|
|
Line 67: |
Line 67: |
| | | |
| :''This model version was used to produce the corrected manuscript in July 2015. | | :''This model version was used to produce the corrected manuscript in July 2015. |
| + | * [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=5YrOTNdv6hO2Gg92 Model run 21.7.2015] runs to the end but emissions are too large exp for wood after 1980. |
| | | |
| <rcode graphics=1 store=0 variables="name:server|type:hidden|default:TRUE"> | | <rcode graphics=1 store=0 variables="name:server|type:hidden|default:TRUE"> |
Line 215: |
Line 216: |
| ####!------------------------------------------------ | | ####!------------------------------------------------ |
| objects.latest("Op_en5488", code_name = "energyUseAnnual") # [[Energy use of buildings]] energyUse | | objects.latest("Op_en5488", code_name = "energyUseAnnual") # [[Energy use of buildings]] energyUse |
− | objects.latest("Op_en2791", code_name = "initiate") # [[Emission factors for burning processes]] | + | objects.latest("Op_en2791", code_name = "emissionstest") # [[Emission factors for burning processes]] |
− | # emissionFactors: Burner, Fuel, Pollutant | + | objects.latest("Op_en2791", code_name = "emissionFactors") # [[Emission factors for burning processes]] |
− | # fuelShares: Heating, Burner, Fuel | + | objects.latest("Op_en7328", code_name = "emissionLocations") # [[Kuopio energy production]] |
| + | objects.latest("Op_en7328", code_name = "fuelShares") # [[Kuopio energy production]] |
| ####i------------------------------------------------ | | ####i------------------------------------------------ |
| | | |
| + | fuelShares <- EvalOutput(fuelShares) |
| energyUse <- EvalOutput(energyUse) | | energyUse <- EvalOutput(energyUse) |
| + | fuelUse <- energyUse * 1E-3 *3600 # kWh -> MJ |
| | | |
− | ################ Transport and fate | + | ## Exposure |
| | | |
− | objects.latest("Op_en5813", code_name = "initiate") # [[Intake fractions of PM]], iF | + | objects.latest("Op_en5813", code_name = "exposure") # [[Intake fractions of PM]] uses Humbert iF as default. |
| + | #objects.latest("Op_en5813", code_name = "initiate") # [[Intake fractions of PM]], iF |
| | | |
| emissions <- EvalOutput(emissions) | | emissions <- EvalOutput(emissions) |
Line 243: |
Line 248: |
| emis <- truncateIndex(emissions, cols = "Emission_site", bins = 5)@output | | emis <- truncateIndex(emissions, cols = "Emission_site", bins = 5)@output |
| | | |
− | ggplot(subset(emis, EfficiencyPolicy == "BAU" & RenovationPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Fuel)) + geom_bar(binwidth = 5) + | + | #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) +
| + | # facet_grid(Pollutant ~ FuelPolicy, scale = "free_y") + theme_gray(base_size = BS) + |
− | labs(
| + | # labs( |
− | title = "Emissions from heating in Kuopio",
| + | # title = "Emissions from heating in Kuopio", |
− | x = "Time",
| + | # x = "Time", |
− | y = "Emissions (ton /a)"
| + | # y = "Emissions (ton /a)" |
− | )
| + | # ) |
| | | |
| if(figstofile) ggsave("Figure5.eps", width = 8, height = 7) | | if(figstofile) ggsave("Figure5.eps", width = 8, height = 7) |
Line 261: |
Line 266: |
| ) | | ) |
| | | |
− | ggplot(subset(emis, EfficiencyPolicy == "BAU" & RenovationPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Fuel)) + geom_bar(binwidth = 5) + | + | #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) +
| + | # facet_grid(Pollutant ~ ., scale = "free_y") + theme_gray(base_size = BS) +#FuelPolicy |
− | labs(
| + | # labs( |
− | title = "Emissions from heating in Kuopio",
| + | # title = "Emissions from heating in Kuopio", |
− | x = "Time",
| + | # x = "Time", |
− | y = "Emissions (ton /a)"
| + | # y = "Emissions (ton /a)" |
− | )
| + | # ) |
| | | |
− | ggplot(subset(emis, EfficiencyPolicy == "BAU" & FuelPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Emission_site)) + geom_bar(binwidth = 5) + | + | ggplot(subset(emis, EfficiencyPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Emission_site)) + geom_bar(binwidth = 5) +# & FuelPolicy == "BAU" |
| facet_grid(Pollutant ~ RenovationPolicy, scale = "free_y") + theme_gray(base_size = BS) + | | facet_grid(Pollutant ~ RenovationPolicy, scale = "free_y") + theme_gray(base_size = BS) + |
| labs( | | labs( |
Line 277: |
Line 282: |
| ) | | ) |
| | | |
− | ggplot(subset(emis, EfficiencyPolicy == "BAU" & FuelPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Fuel)) + geom_bar(binwidth = 5) + | + | ggplot(subset(emis, EfficiencyPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Fuel)) + geom_bar(binwidth = 5) +# & FuelPolicy == "BAU" |
| facet_grid(Pollutant ~ RenovationPolicy, scale = "free_y") + theme_gray(base_size = BS) + | | facet_grid(Pollutant ~ RenovationPolicy, scale = "free_y") + theme_gray(base_size = BS) + |
| labs( | | labs( |
Line 288: |
Line 293: |
| | | |
| ####!------------------------------------------------ | | ####!------------------------------------------------ |
− | objects.latest('Op_en2261', code_name = 'initiate') # [[Health impact assessment]] dose, RR, totcases. | + | #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_en5917', code_name = 'initiate') # [[Disease risk]] disincidence |
− | directs <- tidy(opbase.data("Op_en5461", subset = "Direct inputs"), direction = "wide") # [[Climate change policies and health in Kuopio]] | + | #directs <- tidy(opbase.data("Op_en5461", subset = "Direct inputs"), direction = "wide") # [[Climate change policies and health in Kuopio]] |
| + | objects.latest('Op_en2261', code_name = 'totcases') # [[Health impact assessment]] totcases and dependencies. |
| + | objects.latest('Op_en5461', code_name = 'DALYs') # [[Climate change policies and health in Kuopio]] DALYs, DW, L |
| ####i------------------------------------------------ | | ####i------------------------------------------------ |
| | | |
− | colnames(directs) <- gsub(" ", "_", colnames(directs)) | + | #colnames(directs) <- gsub(" ", "_", colnames(directs)) |
| | | |
| ### Use these population and iF values in health impact assessment. Why? | | ### Use these population and iF values in health impact assessment. Why? |
Line 305: |
Line 312: |
| exposure <- EvalOutput(exposure, verbose = TRUE) | | exposure <- EvalOutput(exposure, verbose = TRUE) |
| | | |
− | ggplot(subset(exposure@output, RenovationPolicy == "BAU" & EfficiencyPolicy == "BAU" & FuelPolicy == "BAU"), aes(x = Time, weight = exposureResult, fill = Heating)) + | + | #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) +
| + | # geom_bar(binwidth = 5) + facet_grid(Area ~ Emission_height) + theme_gray(base_size = BS) + |
− | labs(
| + | # labs( |
− | title = "Exposure to PM2.5 from heating in Kuopio",
| + | # title = "Exposure to PM2.5 from heating in Kuopio", |
− | x = "Time",
| + | # x = "Time", |
− | y = "Average PM2.5 (µg/m3)"
| + | # y = "Average PM2.5 (µg/m3)" |
− | )
| + | # ) |
| | | |
| exposure@output <- exposure@output[exposure@output$Area == "Average" , ] # Kuopio is an average area, | | exposure@output <- exposure@output[exposure@output$Area == "Average" , ] # Kuopio is an average area, |
| # rather than rural or urban. | | # 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) + | + | #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(
| + | # labs( |
− | title = "Exposure to PM2.5 from heating in Kuopio",
| + | # title = "Exposure to PM2.5 from heating in Kuopio", |
− | x = "Time",
| + | # x = "Time", |
− | y = "Average PM2.5 (µg/m3)"
| + | # y = "Average PM2.5 (µg/m3)" |
− | )
| + | # ) |
| | | |
| totcases <- EvalOutput(totcases) | | totcases <- EvalOutput(totcases) |
Line 327: |
Line 334: |
| totcases <- oapply(totcases, cols = c("Age", "Sex"), FUN = sum) | | totcases <- oapply(totcases, cols = c("Age", "Sex"), FUN = sum) |
| | | |
− | ggplot(subset(totcases@output, EfficiencyPolicy == "BAU" & FuelPolicy == "BAU"), aes(x = Time, weight = totcasesResult, fill = Heating))+geom_bar(binwidth = 5) + | + | DALYs <- EvalOutput(DALYs) |
− | facet_grid(Trait ~ RenovationPolicy) +
| + | |
− | theme_gray(base_size = BS) +
| + | #ggplot(subset(totcases@output, EfficiencyPolicy == "BAU" & FuelPolicy == "BAU"), aes(x = Time, weight = totcasesResult, fill = Heating))+geom_bar(binwidth = 5) + |
− | labs(
| + | # facet_grid(Response ~ RenovationPolicy) + |
− | title = "Health effects of PM2.5 from heating in Kuopio",
| + | # theme_gray(base_size = BS) + |
− | x = "Time",
| + | # labs( |
− | y = "Health effects (deaths /a)"
| + | # 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)) | + | #DW <- Ovariable("DW", data = data.frame(directs["Response"], Result = directs$DW)) |
− | L <- Ovariable("L", data = data.frame(directs["Trait"], Result = directs$L)) | + | #L <- Ovariable("L", data = data.frame(directs["Response"], Result = directs$L)) |
| | | |
− | DALYs <- totcases * DW * L | + | #DALYs <- totcases * DW * L |
| | | |
| cat("Total DALYs/a by different combinations of policy options.\n") | | cat("Total DALYs/a by different combinations of policy options.\n") |
Line 346: |
Line 355: |
| temp@output <- subset( | | temp@output <- subset( |
| temp@output, | | temp@output, |
− | as.character(Time) %in% c("2010", "2030") & Trait == "Total mortality" | + | as.character(Time) %in% c("2010", "2030") & Response == "Total mortality" |
| ) | | ) |
| | | |
− | oprint(oapply(temp, INDEX = c("Time", "EfficiencyPolicy", "RenovationPolicy", "FuelPolicy"), FUN = sum)) | + | #oprint(oapply(temp, INDEX = c("Time", "EfficiencyPolicy", "RenovationPolicy", "FuelPolicy"), FUN = sum)) |
| | | |
− | ggplot(subset(DALYs@output, FuelPolicy == "BAU" & Trait == "Total mortality"), aes(x = Time, weight = Result, fill = Heating))+geom_bar(binwidth = 5) + | + | ggplot(subset(DALYs@output, Response == "Total mortality"), aes(x = Time, weight = DALYsResult, fill = Heating))+geom_bar(binwidth = 5) + #FuelPolicy == "BAU" & |
| facet_grid(EfficiencyPolicy ~ RenovationPolicy) + | | facet_grid(EfficiencyPolicy ~ RenovationPolicy) + |
| theme_gray(base_size = BS) + | | theme_gray(base_size = BS) + |
Line 360: |
Line 369: |
| ) | | ) |
| | | |
− | ggplot(subset(DALYs@output, Time == 2030 & Trait == "Total mortality"), aes(x = FuelPolicy, weight = Result, fill = Heating))+geom_bar() + | + | ggplot(subset(DALYs@output, Time == 2030 & Response == "Total mortality"), aes(x = RenovationPolicy, weight = DALYsResult, fill = Heating))+geom_bar() + #x = FuelPolicy |
− | facet_grid(EfficiencyPolicy ~ RenovationPolicy) + | + | facet_grid(EfficiencyPolicy ~ .) + |
| theme_gray(base_size = BS) + | | theme_gray(base_size = BS) + |
| labs( | | labs( |
Line 370: |
Line 379: |
| | | |
| ######## Buildings in Kuopio on map | | ######## Buildings in Kuopio on map |
− | | + | if(FALSE){ |
| # Calculate locations for Kuopio districts | | # Calculate locations for Kuopio districts |
| | | |
Line 401: |
Line 410: |
| cex = 2 | | cex = 2 |
| ) | | ) |
− | | + | } |
| </rcode> | | </rcode> |
| | | |
What are potential climate policies that reach the greenhouse emission targets in the city of Kuopio for years 2010-2030? What are their effects on health and well-being, and what recommendations can be given based on this? The national greenhouse emission target is to reduce greenhouse gas emissions by 20 % between 1990 and 2020; the city of Kuopio has its own, more ambitious target of 40 % for the same time period.
The target of 40 % GHG reduction seems realistic due to reforms in Haapaniemi power plant, assuming that GHG emissions for wood-based fuel is 0. Life-cycle impacts of the wood-based fuel have not yet been estimated.
+ Show code- Hide code
### THIS CODE IS FROM PAGE [[Climate change policies and health in Kuopio]] (Op_en5461, code_name = "")
library(OpasnetUtils)
library(ggplot2)
### Technical parameters
openv.setN(0) # use medians instead of whole sampled distributions
objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]] findrest
BS <- 24 # base_size = font sixe in graphs
figstofile <- FALSE
############################## Case-sepecific data and submodels
obstime <- Ovariable("obstime", data = data.frame(Obsyear = factor(seq(1920, 2050, 10), ordered = TRUE), Result = 1))
## 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, 1910 + 0:14 * 10),
Result = 1
))
###################### Decisions
decisions <- opbase.data('Op_en5461', subset = "Decisions") # [[Climate change policies and health in Kuopio]]
DecisionTableParser(decisions)
# Remove previous decisions, if any.
forgetDecisions <- function() {
for(i in ls(envir = openv)) {
if("dec_check" %in% names(openv[[i]])) openv[[i]]$dec_check <- FALSE
}
return(cat("Decisions were forgotten.\n"))
}
forgetDecisions()
############################ City-specific data
####!------------------------------------------------
objects.latest("Op_en5417", code_name = "initiate") # [[Population of Kuopio]]
# population: City_area
objects.latest("Op_en5932", code_name = "initiatetest") # [[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
#################### 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 = "buildingstest") # [[Building model]] # Generic building model.
# buildings: formula-based
####i------------------------------------------------
renovationRate <- EvalOutput(renovationRate) * 10 # Rates for 10-year periods
renovationRate@marginal[colnames(renovationRate@output) == "Age"] <- TRUE
renovationShares <- EvalOutput(renovationShares)
colnames(renovationShares@output)[colnames(renovationShares@output) == "Startyear"] <- "Obsyear"
stockBuildings <- EvalOutput(stockBuildings)
stockBuildings <- oapply(stockBuildings, cols = c("City_area"), FUN = sum)
changeBuildings <- EvalOutput(changeBuildings)
changeBuildings <- oapply(changeBuildings, cols = c("City_area"), FUN = sum)
buildings <- EvalOutput(buildings)
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, RenovationPolicy == "BAU" & EfficiencyPolicy == "BAU"), aes(x = Time, weight = buildingsResult, 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(buildings@output, aes(x = Time, weight = buildingsResult, fill = Building))+geom_bar()+facet_grid(Efficiency~Heating)
#test <- subset(buildings@output, EfficiencyPolicy == "BAU" & RenovationPolicy == "BAU" & Time == "2000")
if(figstofile) ggsave("Figure3.eps", width = 8, height = 7)
ggplot(subset(bui, EfficiencyPolicy == "BAU"), aes(x = Time, weight = buildingsResult, 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 = buildingsResult, 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 = buildingsResult, 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_en5488", code_name = "energyUseAnnual") # [[Energy use of buildings]] energyUse
objects.latest("Op_en2791", code_name = "emissionstest") # [[Emission factors for burning processes]]
objects.latest("Op_en2791", code_name = "emissionFactors") # [[Emission factors for burning processes]]
objects.latest("Op_en7328", code_name = "emissionLocations") # [[Kuopio energy production]]
objects.latest("Op_en7328", code_name = "fuelShares") # [[Kuopio energy production]]
####i------------------------------------------------
fuelShares <- EvalOutput(fuelShares)
energyUse <- EvalOutput(energyUse)
fuelUse <- energyUse * 1E-3 *3600 # kWh -> MJ
## Exposure
objects.latest("Op_en5813", code_name = "exposure") # [[Intake fractions of PM]] uses Humbert iF as default.
#objects.latest("Op_en5813", code_name = "initiate") # [[Intake fractions of PM]], iF
emissions <- EvalOutput(emissions)
emissions@output$Time <- as.numeric(as.character(emissions@output$Time))
# Plot energy need and emissions
ggplot(subset(energyUse@output, EfficiencyPolicy == "BAU"), aes(x = Time, weight = energyUseResult * 1E-6, fill = Heating)) + geom_bar(binwidth = 5) +
facet_wrap( ~ RenovationPolicy) + theme_gray(base_size = BS) +
labs(
title = "Energy used in heating in Kuopio",
x = "Time",
y = "Heating energy (GWh /a)"
)
if(figstofile) ggsave("Figure4.eps", width = 11, height = 7)
emis <- truncateIndex(emissions, cols = "Emission_site", bins = 5)@output
#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)"
# )
if(figstofile) ggsave("Figure5.eps", width = 8, height = 7)
ggplot(energyUse@output, aes(x = Time, weight = energyUseResult * 1E-6, fill = Heating)) + geom_bar(binwidth = 5) +
facet_grid(EfficiencyPolicy ~ RenovationPolicy) + theme_gray(base_size = BS) +
labs(
title = "Energy used in heating in Kuopio",
x = "Time",
y = "Heating energy (GWh /a)"
)
#ggplot(subset(emis, EfficiencyPolicy == "BAU" & RenovationPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Fuel)) + geom_bar(binwidth = 5) +
# facet_grid(Pollutant ~ ., scale = "free_y") + theme_gray(base_size = BS) +#FuelPolicy
# labs(
# title = "Emissions from heating in Kuopio",
# x = "Time",
# y = "Emissions (ton /a)"
# )
ggplot(subset(emis, EfficiencyPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Emission_site)) + geom_bar(binwidth = 5) +# & FuelPolicy == "BAU"
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"), aes(x = Time, weight = emissionsResult, fill = Fuel)) + geom_bar(binwidth = 5) +# & FuelPolicy == "BAU"
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)"
)
###################### 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
#directs <- tidy(opbase.data("Op_en5461", subset = "Direct inputs"), direction = "wide") # [[Climate change policies and health in Kuopio]]
objects.latest('Op_en2261', code_name = 'totcases') # [[Health impact assessment]] totcases and dependencies.
objects.latest('Op_en5461', code_name = 'DALYs') # [[Climate change policies and health in Kuopio]] DALYs, DW, L
####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 <- 1E+5
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)
DALYs <- EvalOutput(DALYs)
#ggplot(subset(totcases@output, EfficiencyPolicy == "BAU" & FuelPolicy == "BAU"), aes(x = Time, weight = totcasesResult, fill = Heating))+geom_bar(binwidth = 5) +
# facet_grid(Response ~ 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["Response"], Result = directs$DW))
#L <- Ovariable("L", data = data.frame(directs["Response"], Result = directs$L))
#DALYs <- totcases * DW * L
cat("Total DALYs/a by different combinations of policy options.\n")
temp <- DALYs
temp@output <- subset(
temp@output,
as.character(Time) %in% c("2010", "2030") & Response == "Total mortality"
)
#oprint(oapply(temp, INDEX = c("Time", "EfficiencyPolicy", "RenovationPolicy", "FuelPolicy"), FUN = sum))
ggplot(subset(DALYs@output, Response == "Total mortality"), aes(x = Time, weight = DALYsResult, fill = Heating))+geom_bar(binwidth = 5) + #FuelPolicy == "BAU" &
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 == 2030 & Response == "Total mortality"), aes(x = RenovationPolicy, weight = DALYsResult, fill = Heating))+geom_bar() + #x = FuelPolicy
facet_grid(EfficiencyPolicy ~ .) +
theme_gray(base_size = BS) +
labs(
title = "Health effects in DALYs of PM2.5 from heating in Kuopio",
x = "Biofuel policy in district heating",
y = "Health effects (DALY /a)"
)
######## Buildings in Kuopio on map
if(FALSE){
# Calculate locations for Kuopio districts
temp <- buildings
temp@output <- subset(temp@output,
Time == 2030 & EfficiencyPolicy == "BAU" & RenovationPolicy == "BAU"
)
temp <- unkeep(temp, sources = TRUE, prevresults = TRUE)
temp <- oapply(temp, cols = c("Building", "Heating", "Efficiency", "Renovation"), FUN = sum)
####!------------------------------------------------
districts <- tidy(opbase.data("Op_en5932.kuopio_city_districts"), widecol = "Location") # [[Building stock in Kuopio]]
####i------------------------------------------------
colnames(districts) <- gsub("[ \\.]", "_", colnames(districts))
districts <- Ovariable("districts", data = data.frame(districts, Result = 1))
temp <- temp * districts
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
)
}
| |
+ Show code- Hide code
### THIS CODE IS FROM PAGE [[Climate change policies and health in Kuopio]] (Op_en5461, code_name = "")
# 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)
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
figstofile <- FALSE
###################### 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 <- EvalOutput(renovationRate) * 10 # Rates for 10-year periods
#################### 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)
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, RenovationPolicy == "BAU" & EfficiencyPolicy == "BAU"), aes(x = Time, weight = buildingsResult, fill = Heating)) + geom_bar(binwidth = 5) +
theme_gray(base_size = BS) +
labs(
title = "Building stock in Kuopio",
x = "Time",
y = "Floor area (M m2)"
)
if(figstofile) ggsave("Figure3.eps", width = 8, height = 7)
ggplot(subset(bui, EfficiencyPolicy == "BAU"), aes(x = Time, weight = buildingsResult, 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 = buildingsResult, 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 = buildingsResult, 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------------------------------------------------
heatingEnergy <- EvalOutput(heatingEnergy)
################ Transport and fate
objects.latest("Op_en5813", code_name = "initiate") # [[Intake fractions of PM]], iF
emissions <- EvalOutput(emissions)
emissions@output$Time <- as.numeric(as.character(emissions@output$Time))
# Plot energy need and emissions
ggplot(subset(heatingEnergy@output, EfficiencyPolicy == "BAU"), aes(x = Time, weight = heatingEnergyResult * 1E-6, fill = Heating)) + geom_bar(binwidth = 5) +
facet_wrap( ~ RenovationPolicy) + theme_gray(base_size = BS) +
labs(
title = "Energy used in heating in Kuopio",
x = "Time",
y = "Heating energy (GWh /a)"
)
if(figstofile) ggsave("Figure4.eps", width = 11, height = 7)
emis <- truncateIndex(emissions, cols = "Emission_site", bins = 5)@output
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)"
)
if(figstofile) ggsave("Figure5.eps", width = 8, height = 7)
ggplot(heatingEnergy@output, aes(x = Time, weight = heatingEnergyResult * 1E-6, fill = Heating)) + geom_bar(binwidth = 5) +
facet_grid(EfficiencyPolicy ~ RenovationPolicy) + theme_gray(base_size = BS) +
labs(
title = "Energy used in heating in Kuopio",
x = "Time",
y = "Heating energy (GWh /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" & FuelPolicy == "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" & FuelPolicy == "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)"
)
###################### 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
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 <- 1E+5
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(subset(totcases@output, EfficiencyPolicy == "BAU" & FuelPolicy == "BAU"), 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
cat("Total DALYs/a by different combinations of policy options.\n")
temp <- DALYs
temp@output <- subset(
temp@output,
as.character(Time) %in% c("2010", "2030") & Trait == "Total mortality"
)
oprint(oapply(temp, INDEX = c("Time", "EfficiencyPolicy", "RenovationPolicy", "FuelPolicy"), FUN = sum))
ggplot(subset(DALYs@output, FuelPolicy == "BAU" & Trait == "Total mortality"), 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 == 2030 & Trait == "Total mortality"), 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 = "Biofuel policy in district heating",
y = "Health effects (DALY /a)"
)
######## Buildings in Kuopio on map
# Calculate locations for Kuopio districts
temp <- buildings
temp@output <- subset(temp@output,
Time == 2030 & EfficiencyPolicy == "BAU" & RenovationPolicy == "BAU"
)
temp <- unkeep(temp, sources = TRUE, prevresults = TRUE)
temp <- oapply(temp, cols = c("Building", "Heating", "Efficiency", "Renovation"), FUN = sum)
####!------------------------------------------------
districts <- tidy(opbase.data("Op_en5932.kuopio_city_districts"), widecol = "Location") # [[Building stock in Kuopio]]
####i------------------------------------------------
colnames(districts) <- gsub("[ \\.]", "_", colnames(districts))
districts <- Ovariable("districts", data = data.frame(districts, Result = 1))
temp <- temp * districts
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
)
| |