+ Näytä koodi- Piilota koodi
library(OpasnetUtils)
library(ggplot2)
library(xtable)
library(OpasnetUtilsExt)
library(rgdal)
library(maptools)
library(RColorBrewer)
library(classInt)
library(raster)
#Fetch2 <- function(dependencies, evaluate = FALSE, indent = 0, verbose = TRUE, ...) {
# if (nrow(dependencies) > 0) {
# for (i in 1:nrow(dependencies)) {
# if(!exists(as.character(dependencies$Name[i]))) {
# testkey <- if (is.null(dependencies$Key[i])) TRUE else is.na(dependencies$Key[i]) | dependencies$Key[i] == ""
# testid <- if (is.null(dependencies$Ident[i])) TRUE else is.na(dependencies$Ident[i]) | dependencies$Ident[i] == ""
# if (testkey & testid) {
# stop(paste("No key nor ident given for dependent variable: ", dependencies$Name[i], "!", sep = ""))
# }
# if (!testkey) {
# objects.get(dependencies$Key[i]) # Key is the R-tools session identifier (shown at the end of the url)
# }
# if (testkey & !testid) {
# ident <- strsplit(as.character(dependencies$Ident[i]), "/")[[1]] # Ident should be in format <page_id>/<code_name>
# objects.latest(ident[1], ident[2])
# }
# if (evaluate) assign(
# as.character(dependencies$Name[i]),
# EvalOutput(get(as.character(dependencies$Name[i])), ...),
# envir = .GlobalEnv
# )
# if (verbose) cat("\n", rep("-", indent), as.character(dependencies$Name[i]), "fetched successfully!\n")
# }
# }
# }
#} # no need to return anything since variables are written in global memory by objects.get
# Arvioinnin tulos
Tulos <- new("ovariable", name = "Tulos",
dependencies = data.frame(Name = c(
"Terveysvaikutukset",
"Altistuminen",
"Pitoisuus",
"tieliikennepäästöt",
"suorite"
), Ident = c(
"", "", "", "", "Op_fi3192/alustus"
)),
formula = function(dependencies, ...) {
ComputeDependencies(dependencies, ...)
# Leikataan turhat pois
temp <- tieliikennepäästöt@output[tieliikennepäästöt@output$Saaste == "PM", ] # Valitaan vain PM-rivit.
temp <- as.data.frame(as.table(tapply(
temp$tieliikennepäästötResult,
temp[c("LA", "LO", "Vaihtoehto", "Saaste")], # HUOM! Probabilistiset päästöt eivät toimi koska Iter summataan.
sum
)))
temp <- temp[!is.na(temp$Freq), ] # Pudotetaan tyhjät rivit pois.
oprint(temp)
# Käydään päästö läpi rivi kerrallaan ja lasketaan pitoisuus, altistuminen ja terveysvaikutus.
out <- data.frame() # Tähän taulukkoon kerätään data
for(i in 1:nrow(temp)) {
Päästö <- temp$Freq[i]
LA <- as.numeric(as.character(temp$LA[i]))
LO <- as.numeric(as.character(temp$LO[i]))
Altistuminen@output <- data.frame() # Nollataan, jotta malli laskee ne uusiksi uudelle päästölle.
Pitoisuus@output <- data.frame()
# EvalOutput käyttää tilapäisiä, rivikohtaisia tietoja. Kunkin rivin tulos lisätään lopputulokseen.
out <- rbind(
out,
merge(
temp[i, ],
EvalOutput(Terveysvaikutukset, N = N)@output, by = NULL
)
)
}
return(out)
}
)
## Terveysvaikutukset
Terveysvaikutukset <- new(
"ovariable",
name = "Terveysvaikutukset",
dependencies = data.frame(
Name = c("Altistuminen") # , "erf", "bg.mort")
),
formula = function(dependencies, ...) {
ComputeDependencies(dependencies, ...)
bg.mort <- 45182 / 5203826 # same values as used in PILTTI
erf <- 0.0097 # J. T. Tuomisto, A. Wilson, et al. Uncertainty in mortality response to airborne fine
# particulate matter... 2008. unit: m^3 /ug
out <- Altistuminen * erf * bg.mort
return(out)
}
)
## Altistuminen
Altistuminen <- new(
"ovariable",
name = "Altistuminen",
dependencies = data.frame(
Name = c("Pitoisuus", "LO", "LA")
),
formula = function(dependencies, ...) {
ComputeDependencies(dependencies, ...)
out <- GIS.Exposure(Pitoisuus, LO, LA, ...)
return(out)
}
)
## Pitoisuus
Pitoisuus <- new(
"ovariable",
name = "Pitoisuus",
dependencies = data.frame(
Name = c("Päästö", "LO", "LA")
),
formula = function(dependencies, ...) {
ComputeDependencies(dependencies, ...)
temp <- GIS.Concentration.matrix(Päästö, LO, LA, ...)
return(temp)
}
)
Päästö <- 1
LA <- 61.477491
LO <- 21.787756
################# tieliikennepäästöt: funktio tieliikennepäästön laskemiseen
## suorite = ajoneuvojen kulkemat ajokilometrit. Junien osalta ilmoitetaan tonnikilometrit.
tieliikennepäästöt <- new("ovariable",
name = "tieliikennepäästöt",
dependencies = data.frame(Name = c("suorite")),
formula = function(dependencies, ...) {
ComputeDependencies(dependencies, ...)
päästökerroin <- new("ovariable", name = "päästökerroin", ddata = "Op_fi3192") # Haetaan päästökerrointiedot
päästökerroin <- EvalOutput(päästökerroin, N = N, ...)
out <- suorite * päästökerroin # Varsinainen laskentakaava
return(out)
}
)
# suorite
suorite <- new("ovariable",
name = "suorite",
data = {temp <- tidy(opbase.data("Op_fi3357"))
# Luodaan tilapäiset ovariablet, jotta mahdollinen probabilistisuus menee oikein.
lm <- new("ovariable", name = "lm", data = {
temp2 <- temp[colnames(temp) != "Pituus"]
colnames(temp2)[colnames(temp2) == "Liikennemäärä"] <- "Result"
temp2
})
pi <- new("ovariable", name = "pi", data = {
temp3 <- temp[colnames(temp) != "Liikennemäärä"]
colnames(temp3)[colnames(temp3) == "Pituus"] <- "Result"
temp3
})
data <- EvalOutput(lm, N = N) * EvalOutput(pi, N = N) * 365 * 1E-6 # Muutetaan d -> a ja km -> Gm
data@output[!colnames(data@output) %in% c("lmSource", "lmResult", "piSource", "piResult", "Description", "Unit")]
}
)
Tulos <- EvalOutput(Tulos)
if(intermediates) {
cat("Liikennesuorite.\n")
oprint(suorite@data)
cat("Kokonaispienhiukkaspäästöt eri pisteistä.\n")
oprint(tieliikennepäästöt)
cat("Yksityiskohtaiset päästötiedot.\n")
oprint(temp)
cat("Yksityiskohtainen lopputulos.\n")
oprint(Tulos)
}
#if(N > 10) {ggplot(temp, aes(x = Freq, fill = Vaihtoehto)) + geom_density(alpha = 0.2)}
cat("Pienhiukkasten aiheuttamia ylimääräisiä kuolemantapauksia vuodessa.\n")
ggplot(Tulos@output, aes(x = TerveysvaikutuksetResult, fill = Vaihtoehto)) + geom_density()
oprint(Tulos@output[Tulos@output$Iter == 1, ])
#################################
# Draw a concentration map.
par(mfrow=c(6,1), mar=c(3,1,0,1), cex=1.5)
colorstrip <- function(colors, labels) {
count <- length(colors)
m <- matrix(1:count, count, 1)
image(m, col=colors, ylab="", axes=FALSE)
axis(1,approx(c(0, 1), n=length(labels))$y, labels)
}
cat("Esimerkki yhden päästöpisteen aiheuttamasta pitoisuuskentästä.\n")
temp <- EvalOutput(Pitoisuus, N = 1)
temp2 <- gsub("\\(", "", temp@output$LObin)
temp2 <- gsub("\\]", "", temp2)
temp2 <- strsplit(temp2, ",")
temp3 <- c()
for(i in 1:length(temp2)) {
a <- as.numeric(temp2[[i]][1])
b <- as.numeric(temp2[[i]][2])
temp3[i] <- ((a + b) / 2)
}
temp@output$LO <- temp3
temp2 <- gsub("\\(", "", temp@output$LAbin)
temp2 <- gsub("\\]", "", temp2)
temp2 <- strsplit(temp2, ",")
temp3 <- c()
for(i in 1:length(temp2)) {
a <- as.numeric(temp2[[i]][1])
b <- as.numeric(temp2[[i]][2])
temp3[i] <- ((a + b) / 2)
}
temp@output$LA <- temp3
#print(temp)
data <- data.frame(LO=temp@output$LO, LA=temp@output$LA, concentration=temp@output$PitoisuusResult)
data$truncated_concentration <- ifelse(data$concentration > 1, 1, data$concentration)
# Plot the data
coordinates(data)=c("LO","LA")
proj4string(data)<-("+init=epsg:4326")
epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
shp<-spTransform(data,epsg4326String)
#Create blank raster
rast<-raster()
#Set raster extent to that of point data
extent(rast)<-extent(shp)
#Choose number of columns and rows
ncol(rast) <- 41
nrow(rast) <- 41
#cat("<span style='font-size: 1.2em;font-weight:bold;'>PM2.5</span>\n")
#Rasterize point data
rast2<-rasterize(shp, rast, shp$truncated_concentration, fun=mean)
steps <- c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)
colors <- rev(rainbow(length(steps), start=0, end=0.50))
colorstrip(colors, steps)
#Plot data
google.show_raster_on_maps(rast2, col=colors, style="height:500px;")
| |