+ 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)
}
openv.setN(0) # use medians instead of whole sampled distributions
objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]] findrest
obstime <- data.frame(Startyear = (192:205) * 10) # Observation years must be defined for an assessment.
## Additional index needed in followup of ovariables efficiencyShares and stockBuildings
year <- Ovariable("year", data = data.frame(
Constructed = factor(
c("1799-1899", "1900-1909", "1910-1919", "1920-1929", "1930-1939", "1940-1949",
"1950-1959", "1960-1969", "1970-1979", "1980-1989", "1990-1999",
"2000-2010", "2011-2019", "2020-2029", "2030-2039", "2040-2049"
),
ordered = TRUE
),
Time = c(1880, 1905 + 0:14 * 10),
Result = 1
))
BS <- 24
heating_before <- FALSE
efficiency_before <- TRUE
###################### 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",
"changeBuildings",
"efficiencyShares",
"energyUse",
"heatingShares",
"renovationShares",
"renovationRate",
"fuelShares",
"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
# efficiencyShares: Time, Efficiency
renovationRate <- renovationRate * 5 # Rates for 5-year periods
renovationRate@name <- "renovationRate" # This is needed for CheckDecisions
#################### Energy use (needed for buildings submodel)
####!------------------------------------------------
objects.latest("Op_en5488", code_name = "initiate") # [[Energy use of buildings]]
# energyUse: Building, Heating
# efficiencyShares: Efficiency, Constructed
# renovationRatio: Efficiency, Building2, Renovation
####i------------------------------------------------
###################### 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
)
buildings@output$EfficiencyPolicy <- factor(
buildings@output$EfficiencyPolicy,
levels = c("BAU", "Active efficiency"),
ordered = TRUE
)
bui <- oapply(buildings * 1E-6, cols = c("City_area", "buildingsSource"), FUN = sum)@output
ggplot(subset(bui, EfficiencyPolicy == "BAU"), aes(x = Time, weight = Result, fill = Renovation)) + geom_bar(binwidth = 5) +
facet_grid(. ~ RenovationPolicy) + theme_gray(base_size = BS) +
labs(
title = "Building stock in Kuopio by renovation policy",
x = "Time",
y = "Floor area (M m2)"
)
ggplot(subset(bui, RenovationPolicy == "BAU"), aes(x = Time, weight = Result, fill = Efficiency)) + geom_bar(binwidth = 5) +
facet_grid(. ~ EfficiencyPolicy) + theme_gray(base_size = BS) +
labs(
title = "Building stock in Kuopio by efficiency policy",
x = "Time",
y = "Floor area (M m2)"
)
ggplot(subset(bui, RenovationPolicy == "BAU" & EfficiencyPolicy == "BAU"), aes(x = Time, weight = Result, fill = Heating)) + geom_bar(binwidth = 5) +
theme_gray(base_size = BS) +
labs(
title = "Building stock in Kuopio",
x = "Time",
y = "Floor area (M m2)"
)
ggplot(subset(bui, RenovationPolicy == "BAU" & EfficiencyPolicy == "BAU"), aes(x = Time, weight = Result, fill = Building)) + geom_bar(binwidth = 5) +
theme_gray(base_size = BS) +
labs(
title = "Building stock in Kuopio",
x = "Time",
y = "Floor area (M m2)"
)
###################### Energy and emissions
####!------------------------------------------------
objects.latest("Op_en2791", code_name = "initiate") # [[Emission factors for burning processes]]
# emissionFactors: Burner, Fuel, Pollutant
# fuelShares: Heating, Burner, Fuel
####i------------------------------------------------
#fuelShares <- CheckDecisions(EvalOutput(fuelShares, verbose = TRUE))
heatingEnergy <- EvalOutput(heatingEnergy, verbose = TRUE)
################ Transport and fate
####!------------------------------------------------
iF <- Ovariable("iF", ddata = "Op_en3435", subset = "Intake fractions of PM")
# [[Exposure to PM2.5 in Finland]] Humbert et al 2011 data
emissionLocations <- Ovariable("emissionLocations", ddata = "Op_en3435", subset = "Emission locations")
####i------------------------------------------------
colnames(iF@data) <- gsub("[ \\.]", "_", colnames(iF@data))
iF@data$iFResult <- iF@data$iFResult * 1E-6
colnames(emissionLocations@data) <- gsub("[ \\.]", "_", colnames(emissionLocations@data))
emissionLocations@data$emissionLocationsResult <- 1
# Old data:
# 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
emissions <- EvalOutput(emissions)
class(emissions@output$Time)
emissions@output$Time <- as.numeric(as.character(emissions@output$Time))
# Plot energy need and emissions
ggplot(heatingEnergy@output, aes(x = Time, weight = heatingEnergyResult * 1E-6, fill = Heating)) + geom_bar(binwidth = 5) +
facet_grid(EfficiencyPolicy ~ RenovationPolicy, scale = "free_y") + theme_gray(base_size = BS) +
labs(
title = "Energy used in heating in Kuopio",
x = "Time",
y = "Heating energy (GWh /a)"
)
emis <- truncateIndex(emissions, cols = "Emission_site", bins = 5)@output
ggplot(subset(emis, EfficiencyPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Emission_site)) + geom_bar(binwidth = 5) +
facet_grid(Pollutant ~ RenovationPolicy, scale = "free_y") + theme_gray(base_size = BS) +
labs(
title = "Emissions from heating in Kuopio",
x = "Time",
y = "Emissions (ton /a)"
)
ggplot(subset(emis, EfficiencyPolicy == "BAU" & RenovationPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Fuel)) + geom_bar(binwidth = 5) +
facet_grid(Pollutant ~ FuelPolicy, scale = "free_y") + theme_gray(base_size = BS) +
labs(
title = "Emissions from heating in Kuopio",
x = "Time",
y = "Emissions (ton /a)"
)
ggplot(subset(emis, EfficiencyPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Fuel)) + geom_bar(binwidth = 5) +
facet_grid(Pollutant ~ RenovationPolicy, scale = "free_y") + theme_gray(base_size = BS) +
labs(
title = "Emissions from heating in Kuopio",
x = "Time",
y = "Emissions (ton /a)"
)
## LISÄTTÄVÄ FUELPOLICY, KUNHAN SE SAADAAN TOIMIMAAN.
###################### 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------------------------------------------------
colnames(directs) <- gsub(" ", "_", colnames(directs))
### Use these population and iF values in health impact assessment. Why?
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)
population <- 5E+5
#iF <- 1E-7
exposure <- EvalOutput(exposure, verbose = TRUE)
ggplot(subset(exposure@output, RenovationPolicy == "BAU" & EfficiencyPolicy == "BAU" & FuelPolicy == "BAU"), aes(x = Time, weight = exposureResult, fill = Heating)) + geom_bar(binwidth = 5) + facet_grid(Area ~ Emission_height) + theme_gray(base_size = BS) +
labs(
title = "Exposure to PM2.5 from heating in Kuopio",
x = "Time",
y = "Average PM2.5 (µg/m3)"
)
exposure@output <- exposure@output[exposure@output$Area == "Average" , ] # Kuopio is an average area,
# rather than rural or urban.
ggplot(subset(exposure@output, EfficiencyPolicy == "BAU"), aes(x = Time, weight = exposureResult, fill = Heating)) + geom_bar(binwidth = 5) + facet_grid(FuelPolicy ~ RenovationPolicy) + theme_gray(base_size = BS) +
labs(
title = "Exposure to PM2.5 from heating in Kuopio",
x = "Time",
y = "Average PM2.5 (µg/m3)"
)
totcases <- EvalOutput(totcases)
totcases@output$Time <- as.numeric(as.character(totcases@output$Time))
totcases <- oapply(totcases, cols = c("Age", "Sex"), FUN = sum)
ggplot(totcases@output, aes(x = Time, weight = totcasesResult, fill = Heating))+geom_bar(binwidth = 5) +
facet_grid(Trait ~ RenovationPolicy) +
theme_gray(base_size = BS) +
labs(
title = "Health effects of PM2.5 from heating in Kuopio",
x = "Time",
y = "Health effects (deaths /a)"
)
DW <- Ovariable("DW", data = data.frame(directs["Trait"], Result = directs$DW))
L <- Ovariable("L", data = data.frame(directs["Trait"], Result = directs$L))
DALYs <- totcases * DW * L
DALYs@output <- DALYs@output[DALYs@output$Trait != "Lung cancer" , ] # Has to be removed to avoid double counting.
ggplot(DALYs@output, aes(x = Time, weight = Result, fill = Heating))+geom_bar(binwidth = 5) +
facet_grid(EfficiencyPolicy ~ RenovationPolicy) +
theme_gray(base_size = BS) +
labs(
title = "Health effects in DALYs of PM2.5 from heating in Kuopio",
x = "Time",
y = "Health effects (DALY /a)"
)
ggplot(subset(DALYs@output, Time = 2020), aes(x = FuelPolicy, weight = Result, fill = Heating))+geom_bar() +
facet_grid(EfficiencyPolicy ~ RenovationPolicy) +
theme_gray(base_size = BS) +
labs(
title = "Health effects in DALYs of PM2.5 from heating in Kuopio",
x = "Time",
y = "Health effects (DALY /a)"
)
############################# Plot graphs
if(FALSE) {
BS <- 24
######## Buildings in Kuopio on map
if(exists("server")) # The following code only works from Opasnet server.
{
# 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))
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
)
}
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(binwidth = 5) + 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(binwidth = 5) + 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(binwidth = 5) + 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(binwidth = 5) + 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(binwidth = 5) + 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(binwidth = 5) + 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)
ggplot(emissions@output, aes(x = Year, weight = Result, fill = Emission_site)) + geom_bar(binwidth = 5) + facet_grid( Pollutant ~ ., scales = "free_y")
colnames(exposure@output)[colnames(exposure@output) == "Pollutant"] <- "Exposure_agent"
ggplot(exposure@output, aes(x = Year, weight = exposureResult)) + geom_bar(binwidth = 5) + facet_grid(Exposure_agent ~ FuelPolicy, scales = "free_y")
# plo <- ggplot(health@output) + geom_bar(binwidth = 5) + 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)
#ggplot(health@output) + geom_bar(binwidth = 5) + 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)
#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")
#ggplot(DALYs@output[DALYs@output$Causes == "PM2.5+" , ], aes(x = Year, weight = Result/10, fill = Response.metric)) + geom_bar(binwidth = 5) + facet_grid(FuelPolicy ~ RenovationPolicy, scale = "free_y")
ggplot(DALYs@output, aes(x = Year, weight = Result, fill = Trait)) + geom_bar(binwidth = 5) + facet_grid(RenovationPolicy ~ FuelPolicy)
### Impact graphs in FINNISH
if(FALSE) {
plo <- ggplot(buildings@output) + geom_bar(binwidth = 5) + 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(binwidth = 5) + 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(binwidth = 5) + 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(binwidth = 5) + 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(binwidth = 5) + facet_grid(RenovationPolicy ~ #FuelPolicy)
ggplot(health@output[ter@output$Tehokkuuspolitiikka == "BAU" , ], aes(x = Year, weight = Result/get("N", envir=openv),
fill = Lämmitys)) + geom_bar(binwidth = 5) + facet_grid(Remonttipolitiikka ~ Bioenergiapolitiikka) +
labs(y = "Pienhiukkasten terveyshaitta")
ggplot(cases@output, aes(x = Year, weight = totcasesResult, fill = Trait)) + geom_bar(binwidth = 5) + facet_grid(RenovationPolicy ~ FuelPolicy)
ggplot(DALYs@output, aes(x = Year, weight = Result, fill = Trait)) + geom_bar(binwidth = 5) + facet_grid(RenovationPolicy ~ FuelPolicy)
} # END if(FALSE)
| |