+ 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
###################### Decisions
#
#decisions <- opbase.data('Op_en5461', subset = "Decisions") # [[Climate change policies and health in Kuopio]]
#DecisionTableParser(decisions)
############################ City-specific data
#### First download Basel data and use that as a default. Then, replace Helsinki-specific parts.
####!------------------------------------------------
objects.latest("Op_en7115", code_name = "initiate") # [[Building stock in Helsinki]]
# [[Building stock in Helsinki]]
dat <- opbase.data("Op_en7115.stock_details")[ , c(
# "Rakennus ID",
"Sijainti",
"Valmistumisaika",
# "Julkisivumateriaali",
"Käyttötarkoitus",
# "Lämmitystapa",
"Polttoaine",
# "Rakennusaine",
# "Varusteena koneellinen ilmanvaihto",
# "Perusparannus",
# "Kunta rakennuttajana",
# "Energiatehokkuusluokka",
# "Varusteena aurinkopaneeli",
# "Tilavuus",
"Result" # Rakennusala m2
)]
colnames(dat) <- c("City_area", "Time", "Building types in Facta", "Heating types in Facta", "stockBuildingsResult")
dat$Time <- as.numeric(substring(dat$Time, nchar(as.character(dat$Time)) - 3))
dat$Time <- as.numeric(as.character((cut(dat$Time, breaks = c(0, 1885 + 0:26*5), labels = as.character(1885 + 0:26*5)))))
dat$stockBuildingsResult <- as.numeric(as.character(dat$stockBuildingsResult))
build <- tidy(opbase.data("Op_en7115.building_types"))
colnames(build)[colnames(build) == "Result"] <- "Building"
heat <- tidy(opbase.data("Op_en7115.heating_types"))
colnames(heat)[colnames(heat) == "Result"] <- "Heating"
######################
# Korjaus
########################
temp <- as.character(heat$Heating)
temp[temp == "District heating"] <- "District"
temp[temp == "Light oil"] <- "Oil"
temp[temp == "Fuel oil"] <- "Oil"
heat$Heating <- temp
########################################
dat <- merge(merge(dat, build), heat)[c("City_area", "Time", "Building", "Heating", "stockBuildingsResult")]
###################### Decisions
decisions <- opbase.data('Op_en5480') # [[Climate change policies in Basel]]
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()
temp <- aggregate(dat["stockBuildingsResult"], by = dat[c("Time", "Building", "Heating")], FUN =sum)
temp <- temp[!is.na(temp$stockBuildingsResult) , ]
stockBuildings <- Ovariable("stockBuildings", data = temp) # Replace Basel building data with Helsinki data
# Construction rate is assumed to be 2 % /a from the year 2010 building stock.
# changeBuildings is defined as in Basel but only created now to match Helsinki data.
changeBuildings <- stockBuildings
changeBuildings@name <- "changeBuildings"
colnames(changeBuildings@data)[colnames(changeBuildings@data) == "stockBuildingsResult"] <- "changeBuildingsResult"
changeBuildings@data$changeBuildingsResult <- changeBuildings@data$changeBuildingsResult * 0.02
changeBuildings@data$Time <- NULL
changeBuildings@data <- merge(changeBuildings@data, data.frame(Time = 2015 + 0:7 * 5))
renovationRate <- EvalOutput(renovationRate) * 20 # Rates for 20-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 * 1E-6, fill = Heating)) + geom_bar(binwidth = 5) +
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
#objects.latest('Op_en5827', code_name = 'initiate') # [[ERFs of environmental pollutants]] ERF, threshold # ERF is picked automatically
#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 <- 623732 # Contains only the Helsinki city, i.e. assumes no exposure outside city. (Wikipedia)
# population <- 700000 # Contains the Helsinki metropolitan area, as that is exposed to PM2.5.
# 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))
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
