+ Näytä koodi- Piilota koodi
library(OpasnetUtils)
library(ggplot2)
library(xtable)
library(OpasnetUtilsExt)
library(rgdal)
library(maptools)
library(RColorBrewer)
library(classInt)
library(raster)
# Graph settings?
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)
}
# Päästö
Päästö <- new(
"ovariable",
name = "Päästö",
output = data.frame(PäästöResult = päästö), # unit: Mg /a
marginal = FALSE
)
# Muita tarpeellisia arvoja
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
# Ovariablet
## Pitoisuudet
Pitoisuus <- new(
"ovariable",
name = "Pitoisuus",
formula = function(dependencies, ...) {
ComputeDependencies(dependencies, ...)
temp <- GIS.Concentration.matrix(Päästö, LO, LA, ...)
return(temp)
},
dependencies = data.frame(
Name = c("Päästö", "LO", "LA")
)
)
## Altistuminen
riippuvuudet2 <- data.frame(
Name = c("Pitoisuus", "LO", "LA")
)
funktio2 <- function(dependencies, ...) {
ComputeDependencies(dependencies, ...)
Altistuminen <- GIS.Exposure(Pitoisuus, LO, LA, ...)
return(Altistuminen)
}
Altistuminen <- new(
"ovariable",
name = "Altistuminen",
formula = funktio2,
dependencies = riippuvuudet2
)
## Terveysvaikutukset
riippuvuudet3 <- data.frame(
Name = c("Altistuminen", "erf", "bg.mort")
)
funktio3 <- function(dependencies, ...) {
ComputeDependencies(dependencies, ...)
temp <- Altistuminen * erf * bg.mort
return(temp)
}
Terveysvaikutukset <- new(
"ovariable",
name = "Terveysvaikutukset",
formula = funktio3,
dependencies = riippuvuudet3
)
Terveysvaikutukset <- EvalOutput(Terveysvaikutukset, N = N)
print(xtable(Päästö@output), type = 'html')
print(xtable(Altistuminen@output), type = "html")
print(xtable(Terveysvaikutukset@output), type = "html")
#exposure.Population <- tapply(Population.paasto$Pitoisuus * Population.paasto$Result, Population.paasto[,c("obs")], sum)
#cat("Odotusarvo kuolemille vuodessa:", mean(mort.out), "\n")
#cat("Ohjearvon 40 ugm^-3 mukaisen altistusrajan ylitti", sum(Population.paasto.korjaus$Vaesto[Population.paasto.korjaus$Pitoisuus>40]), "asukasta.\n")
#################################3
# Draw a concentration map.
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)
truncated = c()
for(i in 1:nrow(data)){
if (data$concentration[i] > 1)
{
truncated[i] <- 1
}
else
{
truncated[i] <- data$concentration[i]
}
}
data$truncated_concentration <- truncated
# 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) <- 42
nrow(rast) <- 40
#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;")
| |