+ Show code- Hide code
library(OpasnetUtils)
openv.setN(0) # use medians instead of whole sampled distributions
objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]] findrest
obstime <- data.frame(Startyear = 2012) # 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 <- TRUE
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:
objects.latest("Op_en7044", code_name = "initiate") # [[Buildings in Basel]]
# 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
head(renovationRate@output)
head(renovationRate@data)
renovationRate <- EvalOutput(renovationRate) * 5 # Rates for 5-year periods
renovationRate@name <- "renovationRate" # This is needed for CheckDecisions
renovationRate@output$renovationRateResult <- renovationRate@output$Result
renovationRate@output$Result <- NULL
#################### 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
colnames(bui)
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" & 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" & 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 = 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
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
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
#DALYs@output <- DALYs@output[DALYs@output$Trait != "Lung cancer" , ] # Has to be removed to avoid double counting.
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 == 2020 & 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)"
)
| |