+ Show code- Hide code
library(OpasnetUtils)
library(OpasnetUtilsExt)
library(ggplot2)
library(rgdal)
library(maptools)
library(RColorBrewer)
library(classInt)
library(RgoogleMaps)
objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]]
# GIS points for emissions.
districts <- tidy(opbase.data("Op_en3435.kuopio_city_districts"), widecol = "Location") # [[Exposure to PM2.5 in Finland]]
districts <- Ovariable("districts", data = data.frame(districts, Result = 1))
cat("PM2.5 intake fractions are being calculated for these locations.\n")
oprint(districts)
dis <- ova2spat(EvalOutput(districts), coord = c("E", "N"), proj4string = "+init=epsg:3067")
# Long-distance iF of PM2.5 for exposures beyond 10 km.
objects.latest("Op_en5813", code_name="initiate") # Long-distance iF for PM2.5 [[Intake fractions of PM]]
iF.PM2.5@data <- iF.PM2.5@data[iF.PM2.5@data$Subcategory == "Large power plants" , ]
iF.PM2.5@data <- iF.PM2.5@data[!colnames(iF.PM2.5@data) %in% c("Obs", "Geographical area", "Year", "PM type", "Source category", "Subcategory")]
# Calculate exposure concentration * population for a unit emission and all emission points.
out <- Ovariable()
for(i in 1:length(dis$City.area))
{
print(paste(i, "\n"))
temp <- GIS.Exposure(GIS.Concentration.matrix(
1,
LA = coordinates(dis)[i, 2],
LO = coordinates(dis)[i, 1],
N = 1
))
out@output <- rbind(out@output, data.frame(City.area = dis$City.area[i], temp@output))
}
out@output <- out@output[out@output$HAVAINTO == "VAESTO" , ]
out@marginal <- !grepl("Result$", colnames(out@output))
out <- unkeep(out, cols = c("KUNTA", "ID_NRO", "XKOORD", "YKOORD", "HAVAINTO", "dx", "dy"), sources = TRUE)
# Large matrix with detailed exposures in grids.
PILTTI.matrix <- out
# This produces an intake fraction if you give PM2.5 emissions as ton /a. GIS.Concentration.matrix takes ton /a and gives ug /m3.
# iF = intake (g /s) per emission (g /s) = concentration (ug /m3) * population (#) * breathing rate (m3 /s) / emission (g /s).
iF <- oapply(out, cols = c("LAbin", "LObin"), FUN = "sum", na.rm = TRUE)
iF <- iF * 20 / (24 * 3600) * 1E-6 # Divide by breathing rate 20 m3 /d and scale from ug to g to get intake fraction.
iF@output <- data.frame(Emissionheight = "Low", iF@output)
iF@output <- orbind(iF, data.frame(Emissionheight = "High", Result = 0))
iF@marginal <- c(TRUE, iF@marginal)
iF@output <- fillna(iF@output, marginals = colnames(iF@output)[iF@marginal])
iF <- iF + iF.PM2.5 * 1E-6 # Scale iF.PM2.5 from ppm to fractions.
colnames(emissionLocations@data) <- gsub("[ \\.]", "_", colnames(emissionLocations@data))
objects.store(PILTTI.matrix, iF)
cat("Objects PILTTI.matrix, emissionLocations and iF saved.\n")
| |