+ Show code- Hide code
### THIS CODE IS FROM PAGE [[Climate change policies in Helsinki]] (Op_en7063, code_name = "")
library(OpasnetUtils)
library(ggplot2)
library(rgdal)
library(maptools)
library(RColorBrewer)
library(classInt)
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 = c((198:206) * 10)) # Observation years must be defined for an assessment.
BS <- 24
heating_before <- FALSE # Should heatingShares be calculated before renovate and timepoints (or after)?
efficiency_before <- TRUE # Should efficiencyShares be calculated before renovate and timepoints (or after)?
dummy <- 1
figstofile <- FALSE
## Additional index needed in followup of ovariables efficiencyShares and stockBuildings
############################ City-specific data
####!------------------------------------------------
objects.latest("Op_en7115", code_name = "initiate") # [[Building stock in Helsinki]]
###################### Decisions
#
decisions <- opbase.data("Op_en7063", subset = "Decisions") # [[Climate change policies in Helsinki]]
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()
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]]
####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.
####i------------------------------------------------
buildings <- EvalOutput(buildings)
buildings@output$RenovationPolicy <- factor(
buildings@output$RenovationPolicy,
levels = c("BAU", "Active renovation", "Total 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)
bui <- truncateIndex(bui, cols = "Heating", bins = 4)@output
#levels(bui$Heating)[levels(bui$Heating) == "Long-distance heating"] <- "District heating"
ggplot(subset(bui, EfficiencyPolicy == "BAU" & RenovationPolicy == "BAU"), aes(x = Time, weight = buildingsResult, fill = Heating)) + geom_bar(binwidth = 5) +
theme_gray(base_size = BS) +
labs(
title = "Building stock in Helsinki by heating",
x = "Time",
y = "Floor area (M m2)"
)
if(figstofile) ggsave("Figure6.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 Helsinki 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 Helsinki by efficiency policy",
x = "Time",
y = "Floor area (M m2)"
)
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 Helsinki",
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 Helsinki",
x = "Time",
y = "Floor area (M m2)"
)
###################### Energy and emissions
#############################
# Toinen korjaus
#############################
energyUse <- EvalOutput(energyUse)
levels(energyUse@output$Building)[levels(energyUse@output$Building) == "Detached houses"] <- "Residential"
temp <- energyUse@output[energyUse@output$Heating == "Oil",]
temp$Heating <- "Other"
energyUse@output <- rbind(energyUse@output, temp)
#############################
####!------------------------------------------------
objects.latest("Op_en2791", code_name = "initiate") # [[Emission factors for burning processes]]
####i------------------------------------------------
heatingEnergy <- EvalOutput(heatingEnergy)
################ 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
####i------------------------------------------------
colnames(iF@data) <- gsub("[ \\.]", "_", colnames(iF@data))
iF@data$iFResult <- iF@data$iFResult * 1E-6
#################################
# Kolmas korjaus
# fuelShares:sta puuttuu Heating-Other lokaatio
# emissionLocations:sta puuttuu Heating-Oil,Other lokaatiot (myös "District heating")
# heatingEnergy:n Heating-District vastaa fuelShares:ssa jotain turvepohjaista tuotantoa
# emissions@formula:ssa yksi rivi bugaa
#################################
temp <- as.character(heatingEnergy@output$Heating)
temp[temp == "District"] <- "District heating"
heatingEnergy@output$Heating <- temp
emissionLocations <- EvalOutput(emissionLocations)
#levels(emissionLocations@output$Building)[levels(emissionLocations@output$Building) == "Detached houses"] <- "Residential"
temp <- emissionLocations@output[emissionLocations@output$Heating == "District",]
temp$Heating <- "District heating"
emissionLocations@output <- rbind(emissionLocations@output, temp)
temp <- emissionLocations@output[emissionLocations@output$Heating == "Heating oil",]
temp$Heating <- "Oil"
emissionLocations@output <- rbind(emissionLocations@output, temp)
temp <- emissionLocations@output[emissionLocations@output$Heating == "Other sources",]
temp$Heating <- "Other"
emissionLocations@output <- rbind(emissionLocations@output, temp)
temp <- fuelShares@data[1,]
temp$Heating <- "Other"
temp$Burner <- "*"
temp$Fuel <- "*"
temp$Time <- "*"
temp$fuelSharesResult <- 0
fuelShares@data <- rbind(fuelShares@data, temp)
temp <- function(...){
out <- oapply(heatingEnergy, cols = c("Building", "Efficiency"), FUN = sum)
out <- out * fuelShares * emissionFactors * 3.6 * 1e-09
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("Burner", "City_area"), FUN = sum)
return(out)
}
#################################
emissions <- EvalOutput(emissions)
emissions@output$Time <- as.numeric(as.character(emissions@output$Time))
# Plot energy need and emissions
hea <- oapply(heatingEnergy * 1E-6, cols = c("City_area", "buildingsSource"), FUN = sum)
hea <- truncateIndex(hea, cols = "Heating", bins = 4)@output
#levels(hea$Heating)[levels(hea$Heating) == "Long-distance heating"] <- "District heating"
ggplot(hea, aes(x = Time, weight = heatingEnergyResult, fill = Heating)) + geom_bar(binwidth = 5) + # Tuplamuunnos *1e-6 pois
facet_wrap( ~ RenovationPolicy) + theme_gray(base_size = BS) +
labs(
title = "Energy used in heating in Helsinki",
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 Helsinki",
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 Helsinki",
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 Helsinki",
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 <- 623732 # Contains only the Helsinki city, i.e. assumes no exposure outside city. (Wikipedia)
# Note: the population size does NOT affect the health impact as it cancels out. However, it DOES affect
# exposure estimates.
exposure <- EvalOutput(exposure)
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 Helsinki",
x = "Time",
y = "Average PM2.5 (µg/m3)"
)
exposure@output <- exposure@output[exposure@output$Area == "Urban" , ] # Helsinki is an urban 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 Helsinki",
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)
totcases <- truncateIndex(totcases, cols = "Heating", bins = 5)
levels(totcases@output$Heating)[levels(totcases@output$Heating) == "Long-distance heating"] <- "District heating"
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 Helsinki",
x = "Time",
y = "Health effects (deaths /a)"
)
if(figstofile) ggsave("Figure8.eps", width = 11, height = 7)
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),
caption = "Table 1: Total DALYs/a by different combinations of policy options.",
caption.placement = "top",
include.rownames = FALSE
)
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 Helsinki",
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 Helsinki 2030",
x = "Biofuel policy in district heating",
y = "Health effects (DALY /a)"
)
koord <- tidy(opbase.data("Op_en7115", subset = "Locations of postal codes"))
colnames(koord) <- c("Emission_site", "X", "Y")
koord$Result <- 1
koord <- Ovariable("koord", output = koord, marginal = c(TRUE, TRUE, TRUE, FALSE))
emis <- emissions * koord
emis@output <- subset(emis@output, RenovationPolicy == "BAU" & EfficiencyPolicy == "BAU" & FuelPolicy == "BAU" & Pollutant == "PM2.5")
emis <- oapply(emis, INDEX = c("Emission_site", "X", "Y"), FUN = sum)
MyRmap(
ova2spat(
emis,
coord = c("X", "Y"),
proj4string = "+init=epsg:21781"
), # Swiss Land Survey uses CH1903
# http://spatialreference.org/ref/epsg/21782/
# http://en.wikipedia.org/wiki/Swiss_coordinate_system
plotvar = "Result",
legend_title = "PM2.5 emissions (ton/a)",
numbins = 4,
pch = 19,
cex = sqrt(result(emis)) * 3
)
# Map saved manually to .eps with width = 1280, height = 960 px.
| |
Climate policies Helsinki data
Helsingin_rakennuskanta_Facta_ajo_huhtikuu_2015
Bioher model