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 # base_size = font size in graphs
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)?
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_en5488", code_name = "efficiencyShares") # [[Energy use of buildings]] 

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

## Energy use 
objects.latest("Op_en5488", code_name = "energyUseAnnual")        # [[Energy use of buildings]] annual version of energyUse
#objects.latest("Op_en5488", code_name = "efficiencyRatio")  # [[Energy use of buildings]] DEFAULT USED
# objects.latest("Op_en5488", code_name = "renovationRatio")  # [[Energy use of buildings]] DEFAULT USED

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

## Exposure

objects.latest("Op_en5813", code_name = "exposure") # [[Intake fractions of PM]] uses Humbert iF as default.

## Health assessment

objects.latest('Op_en2261', code_name = 'totcases') # [[Health impact assessment]] totcases and dependencies. 
objects.latest('Op_en5461', code_name = 'DALYs') # [[Climate change policies and health in Kuopio]] DALYs, DW, L

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 if the exposure-response function in linear.
# 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$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)

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(bui, aes(x = Time, weight = buildingsResult, fill = Renovation)) + geom_bar(binwidth = 5) + 
		theme_gray(base_size = BS) + facet_grid(. ~ EnergySavingPolicy) +
		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(energyUse * 1E-6, cols = c("City_area", "buildingsSource"), FUN = sum)
hea <- truncateIndex(hea, cols = "Heating", bins = 4)@output

ggplot(hea, aes(x = Time, weight = energyUseResult, 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)"
		)

	emis <- truncateIndex(emissions, cols = "Emission_site", bins = 5)@output

	ggplot(subset(emis, FuelPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Emission_site)) + geom_bar(binwidth = 5) +
			facet_grid(Pollutant ~ EnergySavingPolicy, 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)"
#		)

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

emis <- emissions * koord

emis@output <- subset(emis@output, 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.