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 = 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

##########################!------------------------------------------------
## Decisions
decisions <- opbase.data("Op_en7063", subset = "Decisions") # [[Climate change policies in Helsinki]]

## Building stock
objects.latest("Op_en7115", code_name = "stockBuildings") # [[Building stock in Helsinki]]
objects.latest("Op_en7115", code_name = "changeBuildings") # [[Building stock in Helsinki]]
objects.latest("Op_en7115", code_name = "renovationRate") # [[Building stock in Helsinki]]
objects.latest("Op_en7115", code_name = "renovationShares") # [[Building stock in Helsinki]]

objects.latest("Op_en6289", code_name = "initiate") # [[Building model]] # Generic building model.

## Energy use (needed for buildings submodel)
objects.latest("Op_en5488", code_name = "initiate") # [[Energy use of buildings]] 

## Energy and emissions
objects.latest("Op_en7311", code_name = "emissionLocations") # [[Helsinki energy production]] also heatingShares
objects.latest("Op_en2791", code_name = "initiate") # [[Emission factors for burning processes]] 
objects.latest("Op_en7311", code_name = "fuelShares") # [[Helsinki energy production]]

## Exposure

## 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_en5461', code_name = 'DALYs') # [[Climate change policies and health in Kuopio]] DALYs, DW, L

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.

#########################i--------------------------------------------------


################# Actual model and some intermediate processing.

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

buildings <- EvalOutput(buildings)

#buildings@output$RenovationPolicy <- factor(
#		buildings@output$RenovationPolicy,
#		levels = c("BAU", "Active renovation", "Total renovation"),
#		ordered = TRUE
#)

buildings@output$EnergySavingPolicy <- factor(
		buildings@output$EnergySavingPolicy,
		levels = c("BAU", "Energy saving moderate", "Energy saving total"),
		ordered = TRUE
)

exposure <- EvalOutput(exposure)

exposure@output <- exposure@output[exposure@output$Area == "Average" , ] # Helsinki is an average area,
# rather than rural or urban.

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"

DALYs <- EvalOutput(DALYs)

####################### Post-processing and graphs

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, EnergySavingPolicy == "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(renovationRate@output, aes(x = Time, weight = renovationRateResult)) + geom_bar(binwidth = 5) + 
		facet_grid(. ~ EnergySavingPolicy) + theme_gray(base_size = BS) +
		labs(
				title = "Renovation rate by energy saving policy",
				x = "Time",
				y = "Floor area (M m2)"
		)

ggplot(bui, aes(x = Time, weight = buildingsResult, fill = Efficiency)) + geom_bar(binwidth = 5) + 
		facet_grid(. ~ EnergySavingPolicy) + theme_gray(base_size = BS) +
		labs(
				title = "Building stock in Helsinki by efficiency policy",
				x = "Time",
				y = "Floor area (M m2)"
		)

ggplot(subset(bui, EnergySavingPolicy == "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, EnergySavingPolicy == "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)"
		)

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

ggplot(hea, aes(x = Time, weight = heatingEnergyResult, fill = Heating)) + geom_bar(binwidth = 5) + # Tuplamuunnos *1e-6 pois
		facet_wrap( ~ EnergySavingPolicy) + theme_gray(base_size = BS) +
		labs(
				title = "Energy used in heating in Helsinki",
				x = "Time",
				y = "Heating energy (GWh /a)"
		)

if("Emission_site" %in% colnames(emissions@output)) {
	emis <- truncateIndex(emissions, cols = "Emission_site", bins = 5)@output

	ggplot(subset(emis, EnergySavingPolicy == "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(emissions@output, EnergySavingPolicy == "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(emissions@output, EnergySavingPolicy == "BAU" & FuelPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Fuel)) + geom_bar(binwidth = 5) +
		facet_grid(Pollutant ~ . , scale = "free_y") + theme_gray(base_size = BS) +
		labs(
				title = "Emissions from heating in Helsinki",
				x = "Time",
				y = "Emissions (ton /a)"
		)

#exposure <- EvalOutput(exposure)

ggplot(subset(exposure@output, EnergySavingPolicy == "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)"
		)

ggplot(subset(exposure@output, EnergySavingPolicy == "BAU"), aes(x = Time, weight = exposureResult, fill = Heating)) + geom_bar(binwidth = 5) + 	facet_grid(FuelPolicy ~ .) + theme_gray(base_size = BS) +
		labs(
				title = "Exposure to PM2.5 from heating in Helsinki",
				x = "Time",
				y = "Average PM2.5 (µg/m3)"
		)

ggplot(subset(totcases@output, EnergySavingPolicy == "BAU" & FuelPolicy == "BAU"), aes(x = Time, weight = totcasesResult, fill = Heating))+geom_bar(binwidth = 5) + 
		facet_grid(Response ~ .) +
		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)

#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", "EnergySavingPolicy", "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, Response == "Total mortality"), aes(x = Time, weight = DALYsResult, fill = Heating))+geom_bar(binwidth = 5) + 
		facet_grid(EnergySavingPolicy ~ FuelPolicy) +
		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), aes(x = FuelPolicy, weight = DALYsResult, fill = Heating))+geom_bar() + 
		facet_grid(EnergySavingPolicy ~ Response) +
		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 <- opbase.data("Op_en7115", subset = "Locations of city areas")
koord <- tidy(koord, widecol = "Location")
colnames(koord)[colnames(koord) == "City_area"] <- "Emission_site"
#colnames(koord)[colnames(koord) == "N"] <- "X"
#colnames(koord)[colnames(koord) == "E"] <- "Y"
koord$Result <- 1
koord <- Ovariable("koord", output = koord, marginal = colnames(koord) != "Result")

if("Emission_site" %in% colnames(emissions@output)) {

	emis <- emissions * koord

	emis@output <- subset(emis@output, RenovationPolicy == "BAU" & EnergySavingPolicy == "BAU" & 
		FuelPolicy == "BAU" & Pollutant == "PM2.5")
	emis <- oapply(emis, INDEX = c("Emission_site", "N", "E"), FUN = sum)

	MyRmap(
		ova2spat(
			emis, 
			coord = c("E", "N"), 
			proj4string = "+init=epsg:3067"
		), # National Land Survey Finland uses ETRS89 / ETRS-TM35FIN
		#http://spatialreference.org/ref/epsg/etrs89-etrs-tm35fin/
		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.