Difference between revisions of "Climate change policies in Helsinki"

### 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 = (192:205) * 10) # Observation years must be defined for an assessment.

obstime <- data.frame(Startyear = c(2010, 2030)) # 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

#year <- Ovariable("year", data = data.frame(
#	Constructed = factor(
#		c("1759-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
#))

###################### 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]]. Somehow only works if done in parts.
#datb <- opbase.data("Op_en7115.stock_details", include = list(Julkisivumateriaali = "")) # 18747 rows
#datc <- opbase.data("Op_en7115.stock_details", exclude = list(Julkisivumateriaali = "")) # 30038 rows
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"
)]
#dat <- rbind(datb, datc)[ , c(
#	"Rakennus ID",
#	"Sijainti",
#	"Valmistumisaika",
#	"Julkisivumateriaali",
#	"Käyttötarkoitus",
#	"Lämmitystapa",
#	"Polttoaine",
#	"Rakennusaine",
#	"Varusteena koneellinen ilmanvaihto",
#	"Perusparannus",
#	"Kunta rakennuttajana",
#	"Energiatehokkuusluokka",
#	"Varusteena aurinkopaneeli",
#	"Tilavuus",
#	"Result"
#)]

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
#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.