+ Show code- Hide code
# SIirrä Kuopion-datat kässäristä linkin taakse Opasnettiin
# KÄssäriin vain yhteenvetotaulukko joka kaupungista.
# Mieti mitä sanotaan sisäilmasta. Perusmalli toimii ilmankin, ja Matin nostama miljoonan sisäilman hankaluus pitäisi lähinnä keskustella. Käytetäänkä ylileveitä jakaumia?
# Onko järkeä yhdistää kaupungit? Silloin tulisi NA:ta eri päätösten kohdalle, ja tämä pitäisi huomioida kuvissa (muutenkin kannattaisi slaissata data ennen kuvien piirtämistä).
# Tarkista iF jota käytetään: Mikä on iF-summary?
library(OpasnetUtils)
library(ggplot2)
library(rgdal)
library(maptools)
library(RColorBrewer)
library(classInt)
library(OpasnetUtilsExt)
library(RgoogleMaps)
suomenna <- function(ova) {
d <- ova@output
if("Trait" %in% colnames(d)) {
levels(d$Trait)[levels(d$Trait) == "Cancer"] <- "Syöpä"
levels(d$Trait)[levels(d$Trait) == "CHD2"] <- "Sydänkuolema"
levels(d$Trait)[levels(d$Trait) == "Stroke"] <- "Aivohalvaus"
levels(d$Trait)[levels(d$Trait) == "Child's IQ"] <- "Lapsen ÄO"
levels(d$Trait)[levels(d$Trait) == "Tooth defect"] <- "Hammasvaurio"
levels(d$Trait)[levels(d$Trait) == "Dental defect"] <- "Hammasvaurio"
levels(d$Trait)[levels(d$Trait) == "Vitamin D recommendation"] <- "D-vitamiinin saantisuositus"
levels(d$Trait)[levels(d$Trait) == "Dioxin recommendation"] <- "Dioksiinin saantisuositus"
}
if("Exposure_agent" %in% colnames(ova@output)) {
levels(d$Exposure_agent)[levels(d$Exposure_agent) == "Vitamin_D"] <- "D-vitamiini"
levels(d$Exposure_agent)[levels(d$Exposure_agent) == "EPA"] <- "EPA"
levels(d$Exposure_agent)[levels(d$Exposure_agent) == "DHA"] <- "DHA"
levels(d$Exposure_agent)[levels(d$Exposure_agent) == "Omega3"] <- "Omega3"
levels(d$Exposure_agent)[levels(d$Exposure_agent) == "PCDDF"] <- "Dioksiini"
levels(d$Exposure_agent)[levels(d$Exposure_agent) == "PCB"] <- "PCB"
levels(d$Exposure_agent)[levels(d$Exposure_agent) == "TEQ"] <- "TEQ"
levels(d$Exposure_agent)[levels(d$Exposure_agent) == "MeHg"] <- "Metyylielohopea"
}
if("Hedelm" %in% colnames(d)) {
levels(d$Hedelm)[levels(d$Hedelm) == FALSE] <- "Ei"
levels(d$Hedelm)[levels(d$Hedelm) == TRUE] <- "Kyllä"
}
colnames(d)[colnames(d) == "Age"] <- "Ikäryhmä"
colnames(d)[colnames(d) == "Exposure_agent"] <- "Altiste"
colnames(d)[colnames(d) == "Trait"] <- "Vaste"
colnames(d)[colnames(d) == "Ikä"] <- "Ikä"
colnames(d)[colnames(d) == "Silakkamäärä"] <- "Silakkamäärä"
return(d)
}
#language <- "EN"
openv.setN(0) # use medians instead of whole sampled distributions
#openv.setN(1000)
BS <- 18
objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]] findrest
obsyear <- (192:205) * 10 # Observation years must be defined for an assessment.
###################### Decisions
decisions <- opbase.data('Op_en5461', subset = "Decisions") # [[Climate change policies and health in Kuopio]]
DecisionTableParser(decisions)
# Remove previous decisions, if any.
rm(
"buildings",
"stockBuildings",
"buildingTypes",
"changeBuildings",
"stockEfficiencyShares",
"changeEfficiencyShares",
"energyUse",
"heatingShares",
"heatingSharesNew",
"rateBuildings",
"renovationShares",
"fuelTypes",
"year",
envir = openv
)
############################ City-specific data
####!------------------------------------------------
objects.latest("Op_en5417", code_name = "initiate") # [[Population of Kuopio]]
# population: City_area
objects.latest("Op_en5932", code_name = "initiate") # [[Building stock in Kuopio]] Building ovariables:
# buildingStock: Building, Constructed, City_area
# rateBuildings: Age, (RenovationPolicy)
# renovationShares: Renovation
# construction: Building
# constructionAreas: City_area
# buildingTypes: Building, Building2
# heatingShares: Building, Heating, Eventyear
# heatingSharesNew: Building2, Heating
# eventyear: Constructed, Eventyear
# efficiencies: Constructed, Efficiency
# Calculate locations for Kuopio districts
districts <- tidy(opbase.data("Op_en3435.kuopio_city_districts"), widecol = "Location") # [[Exposure to PM2.5 in Finland]]
####i------------------------------------------------
colnames(districts) <- gsub("\\.", "_", colnames(districts))
districts <- Ovariable("districts", data = data.frame(districts, Result = 1))
###################### Energy and emissions
####!------------------------------------------------
objects.latest("Op_en2791", code_name = "initiate") # [[Emission factors for burning processes]]
# emissionFactors: Burner, Fuel, Pollutant
# fuelTypes: Heating, Burner, Fuel
objects.latest("Op_en5488", code_name = "initiate") # [[Energy use of buildings]]
# energyUse: Building, Heating
# changeEfficiencyShares: Efficiency, Constructed
# renovationRatio: Efficiency, Building2, Renovation
####i------------------------------------------------
#### BUILDINGTYPES MUST BE HANDLED OUTSIDE FORMULA!
fuelTypes@data <- merge(fuelTypes@data, data.frame(Year = 1900 + 0:16 * 5))
###### THESE SHOULD HAPPEN WHEN THE OVARIABLES ARE DEFINED
colnames(iF@output)[colnames(iF@output) == "City.area"] <- "Emission_site"
colnames(iF@output)[colnames(iF@output) == "Emissionheight"] <- "Emission_height"
iF@output$Iter <- NULL
colnames(emissionLocations@data) <- gsub(" ", "_", colnames(emissionLocations@data))
# Fill in Heating types and convert from % to fraction.
heatingSharesNew <- findrest(heatingSharesNew, cols = "Heating", total = 100) / 100
changeEfficiencyShares <- findrest(changeEfficiencyShares, cols = "Efficiency", total = 100) / 100
renovationShares <- findrest(renovationShares, cols = "Renovation", total = 100) / 100
################ Transport and fate
####!------------------------------------------------
objects.latest("Op_en3435", code_name = "disperse") # [[Exposure to PM2.5 in Finland]]
# iF: Iter, Emissionheight, City.area ## THESE SHOULD BE UPDATED! (precalculated with N = 1)
# emissionLocations: Heating, Emission site, Emission height
#Summarised Piltti matrix, another copy of the code on a more reasonable page
# Default run: en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=aXDIVDboftr1bTEd
####i------------------------------------------------
fuelTypes <- CheckDecisions(EvalOutput(fuelTypes, verbose = TRUE))
###################### Health assessment
####!------------------------------------------------
objects.latest('Op_en2261', code_name = 'initiate') # [[Health impact assessment]] dose, RR, totcases.
objects.latest('Op_en5917', code_name = 'initiate') # [[Disease risk]] disincidence
objects.latest('Op_en5827', code_name = 'initiate') # [[ERFs of environmental pollutants]] ERF, threshold
#objects.latest('Op_en5453', code_name = 'initiate') # [[Burden of disease in Finland]] BoD
directs <- tidy(opbase.data("Op_en5461", subset = "Direct inputs"), direction = "wide") # [[Climate change policies and health in Kuopio]]
####i------------------------------------------------
frexposed <- 1 # fraction of population that is exposed
bgexposure <- 0 # Background exposure to an agent (a level below which you cannot get in practice)
BW <- 70 # Body weight (is needed for RR calculations although it is irrelevant for PM2.5)
colnames(directs) <- gsub(" ", "_", colnames(directs))
#################### Manage the data before calculating
# This is the population of Kuopio (i.e. population living in the building stock)
population <- EvalOutput(population, verbose = TRUE)
areaweight <- oapply(population, cols = "City_area", FUN = "sum") # Sum across city areas.
areaweight <- population / (1 * areaweight)
buildingStock <- buildingStock * areaweight
construction <- construction * constructionAreas #/ 3 # Statistics are for three years (2010-2012) BUT PER YEAR?!
###################### Actual building model
# The building stock is measured as m^2 floor area.
####!------------------------------------------------
objects.latest("Op_en6289", code_name = "initiate") # [[Building model]] # Generic building model.
# buildings: formula-based
# heatingEnergy: formula-based
####i------------------------------------------------
buildings <- EvalOutput(buildings, verbose = TRUE)
buildings@output$RenovationPolicy <- factor(
buildings@output$RenovationPolicy,
levels = c("BAU", "Active renovation", "Effective renovation"),
ordered = TRUE
)
heatingEnergy <- EvalOutput(heatingEnergy, verbose = TRUE)
heatingEnergy <- unkeep(heatingEnergy, cols = "Building2", sources = TRUE, prevresults = TRUE)
if(exists("server")) # The following code only works from Opasnet server.
{
temp <- buildings * districts
temp@output <- temp@output[temp@output$Year == 2040 , ]
temp <- unkeep(temp, sources = TRUE, prevresults = TRUE)
temp@output <- dropall(temp@output)
temp <- oapply(temp, cols = c("Building", "Heating", "Efficiency", "Renovation"), FUN = "sum", na.rm = TRUE)
MyRmap(
ova2spat(temp, coord = c("E", "N"), proj4string = "+init=epsg:3067"), # National Land Survey uses EPSG:3067 (ETRS-TM35FIN)
plotvar = "Result", legend_title = "Floor area", numbins = 8, pch = 19, cex = 2
)
}
####### Calculate emissions
#emis <- heatingEnergy
#emis@output <- emis@output[emis@output$Year >= 1980 , ]
#emis <- oapply(emis, cols = c("Efficiency", "Renovation", "Building"), FUN = "sum")
#emis@output$Year <- as.numeric(as.character(emis@output$Year))
emissions <- Ovariable("emissions",
dependencies = data.frame(
Name = c(
"heatingEnergy",
"fuelTypes",
"emissionFactors"
)
),
formula = function(...) {
out <- oapply(heatingEnergy, cols = c("Building", "Efficiency"), FUN = sum)
out <- out * fuelTypes * emissionFactors * 3.6 * 1E-9 # convert from kWh /a to MJ /a and mg to ton
out <- unkeep(out * emissionLocations, sources = TRUE, prevresults = TRUE)
out@output$Emission_site <- as.factor(ifelse(
out@output$Emission_site == "At site of consumption",
as.character(out@output$City_area),
as.character(out@output$Emission_site)
))
out <- oapply(
out,
cols = c("Heating", "Burner", "Fuel", "City_area"),
FUN = sum
)
}
)
#"Emission_site", "Emission_height",
emissions <- EvalOutput(emissions)
#############################################
ggplot(emissions@output, aes(x = Year, weight = Result, fill = Emission_site)) + geom_bar() + facet_grid( Pollutant ~ ., scales = "free_y")
#ggplot(health@output) + geom_bar() + theme_gray(base_size = BS) + labs(title = "Health impacts of fuel and renovation policy") +
#aes(x = Year, weight = Result, fill = Heating) + labs(y = "Premature deaths (# /a)") # + facet_grid(FuelPolicy ~ RenovationPolicy)
##############------------------------------
### Use these population and iF values in health impact assessment. Why?
pop <- 5E+5
iF <- 1E-7
exposure <- Ovariable("exposure",
# emis is in ton /a
# iF = conc (g /m3) * pop (#) * BR (m3 /s) / emis (g /s) <=> conc = emis * iF / BR / pop # conc is the exposure concentration
dependencies = data.frame(Name = c("emis2", "iF", "pop")),
formula = function(...) {
BR <- 20 # Nominal breathing rate (m^3 /d)
BR <- BR / 24 / 3600 # m^3 /s
out <- 1E+12 / 365 / 24 / 3600 # Emission scaling from ton /a to ug /s.
out <- (emis2 * out) * iF / BR / pop # the actual equation
out <- unkeep(out, prevresults = TRUE, sources = TRUE)
return(out)
}
)
exposure <- EvalOutput(exposure, verbose = TRUE)
colnames(exposure@output)[colnames(exposure@output) == "Pollutant"] <- "Exposure_agent"
ggplot(exposure@output, aes(x = Year, weight = exposureResult)) + geom_bar() + facet_grid(Exposure_agent ~ FuelPolicy, scales = "free_y")
totcases <- EvalOutput(totcases)
#totcases@output <- totcases@output[!is.na(result(totcases)) , ] # Not needed with new oapply
#cases <- totcases * Ovariable("incre", data = data.frame(Exposcen = c("BAU", "No exposure"), Result = c(1, -1)))
cases <- oapply(totcases, cols = c("Exposcen", "Age", "Sex", "City_area"), FUN = sum)
DW <- Ovariable("DW", data = data.frame(directs["Trait"], Result = directs$DW))
L <- Ovariable("L", data = data.frame(directs["Trait"], Result = directs$L))
#ggplot(cases@output[cases@output$Causes == "PM2.5+" , ], aes(x = Year, weight = Result/10, fill = Disease)) + geom_bar() + facet_grid(FuelPolicy ~ RenovationPolicy, scale = "free_y")
DALYs <- cases * DW * L
#DALYs <- oapply(DALYs, cols = c("Trait"), FUN = sum)
DALYs@output <- DALYs@output[DALYs@output$Trait != "Lung cancer" , ] # Has to be removed to avoid double counting.
#ggplot(DALYs@output[DALYs@output$Causes == "PM2.5+" , ], aes(x = Year, weight = Result/10, fill = Response.metric)) + geom_bar() + facet_grid(FuelPolicy ~ RenovationPolicy, scale = "free_y")
ggplot(DALYs@output, aes(x = Year, weight = Result, fill = Trait)) + geom_bar() + facet_grid(RenovationPolicy ~ FuelPolicy)
############################# Plot graphs
######## Buildings in Kuopio on map
if(exists("server")) # The following code only works from Opasnet server.
{
if(language == "EN") lege2 <- "Age of building" else lege2 <- "Rakennuksen ikä"
shp <- readOGR('PG:host=localhost user=postgres dbname=spatial_db','kuopio_house')
proj4string(shp) <- ("+init=epsg:3067") # The map projection of NLS Finland.
epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
shp2 <- spTransform(shp,epsg4326String) # Convert to longitude-latitude projection.
MyRmap(shp2, plotvar = "ika", legend_title = lege2, numbins = 8, pch = 19, cex = 0.3) # Draw the map.
}
if(exists("server")) BS <- 24 else BS <- 12
### Impact graphs in ENGLISH
ggplot(buildings@output) + geom_bar() + theme_gray(base_size = BS) +
aes(x = Building, weight = buildingsResult/1000000, fill = Heating) + labs(y = "Floor area (M m2)", title = "Building impacts of renovation policy") + coord_flip() # + facet_grid(. ~ RenovationPolicy)
plo <- ggplot(buildings@output) + geom_bar() + facet_grid(. ~ RenovationPolicy) + theme_gray(base_size = BS) +
aes(x = Year, weight = buildingsResult/1000000, fill = Renovation) + labs(y = "Floor area (M m2)",
title = "Building impacts of renovation policy")
print(plo)
ggplot(heatingEnergy@output) + geom_bar() + theme_gray(base_size = BS) + labs(title = "Energy impacts of renovation policy") +
aes(x = Building, weight = heatingEnergyResult/1E+6, fill = Heating) + labs(y = "Heating energy need (GWh /a)") + coord_flip() # + facet_grid(. ~ RenovationPolicy)
plo <- ggplot(heatingEnergy@output) + geom_bar() + facet_grid(. ~ RenovationPolicy) + theme_gray(base_size = BS) +
labs(title = "Energy impacts of renovation policy", y = "Heating energy need (GWh /a)") +
aes(x = Year, weight = heatingEnergyResult/1E+6, fill = Heating)
print(plo)
emis@output <- emis@output[emis@output$Renovation == "BAU" , ]
ggplot(emis@output) + geom_bar() + facet_grid(Pollutant ~ . , scales = "free_y") + theme_gray(base_size = BS) + labs(title = "Emission impacts of biofuel policy") +
aes(x = Heating, weight = Result, fill = Fuel) + labs(y = "Emissions to air (ton /a)") + theme(axis.text.x = element_text(angle = 90, hjust = 1))
plo <- ggplot(emis@output) + geom_bar() + facet_grid(Pollutant ~ FuelPolicy, scales = "free_y") + theme_gray(base_size = BS) +
labs(title = "Emission impacts of biofuel policy", y = "Emissions to air (ton /a)") +
aes(x = Year, weight = Result, fill = Fuel)
print(plo)
# plo <- ggplot(health@output) + geom_bar() + facet_grid(FuelPolicy ~ RenovationPolicy) + theme_gray(base_size = BS) +
# labs(title = "Health impacts of fuel and renovation policy", y = "Premature deaths (# /a)") +
# aes(x = Year, weight = Result, fill = Heating)
#
# print(plo)
### Impact graphs in FINNISH
if(FALSE) {
plo <- ggplot(buildings@output) + geom_bar() + facet_grid(. ~ Remonttipolitiikka) + theme_gray(base_size = BS) +
aes(x = Year, weight = buildingsResult/1000000, fill = Renovation) + labs(y = "Rakennuspinta-ala (M m2)",
title = "Remonttipolitiikan vaikutus rakennuskantaan")
print(plo)
plo <- ggplot(heatingEnergy@output) + geom_bar() + facet_grid(. ~ Remonttipolitiikka) + theme_gray(base_size = BS) +
labs(title = "Remonttipolitiikan vaikutus energiankulutukseen", y = "Lämmitysenergian tarve (GWh /a)") +
aes(x = Year, weight = heatingEnergyResult/1E+6, fill = Heating)
print(plo)
emis@output <- emis@output[emis@output$Remonttipolitiikka == "BAU" , ]
plo <- ggplot(emis@output) + geom_bar() + facet_grid(Pollutant ~ Bioenergiapolitiikka, scales = "free_y") +
theme_gray(base_size = BS) + labs(title = "Bioenergiapolitiikan päästövaikutukset", y = "Päästöt ilmaan (ton /a)") +
aes(x = Year, weight = Result, fill = Fuel)
print(plo)
# plo <- ggplot(health@output) + geom_bar() + facet_grid(Bioenergiapolitiikka ~ Remonttipolitiikka) + theme_gray(base_size = BS) +
# labs(title = "Remontti- ja bioenergiapolitiikan terveysvaikutukset", y = "Ennenaikaisia kuolemia (# /a)") +
# aes(x = Year, weight = Result, fill = Heating)
#
# print(plo)
}
ggplot(dose@output, aes(x = Year, y = doseResult, colour = RenovationPolicy)) + geom_point() + facet_grid(. ~ FuelPolicy)
ggplot(totcases@output, aes(x = Year, y = totcasesResult, colour = Trait)) + geom_point() + facet_grid(Trait ~ FuelPolicy)
ggplot(RR@output, aes(x = RRResult, colour = Trait)) + geom_density()
#ggplot(cases@output, aes(x = Year, weight = Result, fill = Trait)) + geom_bar() + facet_grid(RenovationPolicy ~ #FuelPolicy)
############ Health impact graph corrected. This should be done in a previous stage.
ter <- DALYs
ter2 <- ter
ter2@output <- ter@output[
ter@output$Bioenergiapolitiikka == "BAU" &
ter@output$Remonttipolitiikka == "BAU" ,
]
ter2 <- unkeep(ter2, cols = c("Bioenergiapolitiikka", "Remonttipolitiikka"))
bau <- Ovariable(data = data.frame(
Bioenergiapolitiikka = rep(c("BAU", "Bioenergian lisäys"), each = 2),
Remonttipolitiikka = rep(c("Aktiivinen remontointi", "BAU"), times = 2),
Result = 1
))
bau <- bau * Ovariable(data = data.frame(
Year = c(1980, 1990, 2000, 2010, 2020, 2030, 2040, 2050),
Result = c(1,1,1,1,0,0,0,0)
))
ter <- ter * (1 - bau) + ter2 * bau
levels(ter@output$Heating) <- c("Kaukolämpö", "Sähkö", "Maalämpö", "Öljy", "Puu")
colnames(ter@output)[colnames(ter@output) == "Heating"] <- "Lämmitys"
ggplot(ter@output[ter@output$Tehokkuuspolitiikka == "BAU" , ], aes(x = Year, weight = Result/get("N", envir=openv),
fill = Lämmitys)) + geom_bar() + facet_grid(Remonttipolitiikka ~ Bioenergiapolitiikka) +
labs(y = "Pienhiukkasten terveyshaitta")
ggplot(cases@output, aes(x = Year, weight = totcasesResult, fill = Trait)) + geom_bar() + facet_grid(RenovationPolicy ~ FuelPolicy)
ggplot(DALYs@output, aes(x = Year, weight = Result, fill = Trait)) + geom_bar() + facet_grid(RenovationPolicy ~ FuelPolicy)
############# Building stock in Kuopio
################ EI TARVITTANE?
energyEfficiency@data$Efficiency <- factor(
energyEfficiency@data$Efficiency,
levels = c("Old", "New", "Low-energy", "Passive"),
ordered = TRUE
)
| |