Difference between revisions of "Helsinki energy decision 2015"

#Smaller fuel oil heat plants around Helsinki Smaller gas heat plants around Helsinki
library(OpasnetUtils)
library(ggplot2)

openv.setN(0) # use medians instead of whole sampled distributions
BS <- 24 # b
objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]]

# Download all pre-calculated inputs, e.g. building stock.
objects.latest("Op_en7237", code_name = "intermediates") # [[Helsinki energy decision 2015]]

#plants <- union(as.character(unique(plantParameters@data$Plant)), as.character(unique(energyProcess@data$Plant)))

shutdown <- c(shutdown, "None")

customdecisions <- data.frame(
	Obs = NA, 
	Decision_maker = "Helen", 
	Decision = "PlantPolicy", 
	Option = "Custom", 
	Variable = c("plantParameters", "energyProcess"),
	Cell = c(paste("Plant:", paste(shutdown, collapse = ","), ";Parameter:Max;Time:>=2035"), ""),
	Change = c("Replace", "Identity"),
	Unit = NA,
	Result = 0
)

decisions <- rbind(decisions, customdecisions)

if ("Hanasaari" %in% renovation) {
	customdecisions <- data.frame(
		Obs = NA, 
		Decision_maker = "Helen", 
		Decision = "PlantPolicy", 
		Option = "Custom", 
		Variable = "plantParameters",
		Cell = c("Plant:Hanasaari;Time:>=2018;Time:<=2060;Parameter:Max", "Plant:Hanasaari;Time:>=2018;Parameter:Investment cost"),
		Change = "Replace",
		Unit = NA,
		Result = 100
	)
	decisions <- rbind(decisions, customdecisions)
	customdecisions <- data.frame(
		Obs = NA, 
		Decision_maker = "Helen", 
		Decision = "PlantPolicy", 
		Option = "Custom", 
		Variable = "energyProcess",
		Cell = c("Plant:Hanasaari;Fuel:Biofuel;Time:>=2018", "Plant:Hanasaari;Fuel:Coal;Time:>=2018"),
		Change = "Add",
		Unit = NA,
		Result = c(-0.40, 0.40)
	)
	decisions <- rbind(decisions, customdecisions)
}

DecisionTableParser(decisions)

cat("Plants that will be shut down in the option Plant policy/Custom are:", shutdown, "\n")

fuelUse <- EvalOutput(fuelUse)

fuelUse$Fuel <- factor(
		fuelUse$Fuel, levels = c(
				"Biofuel",
				"Coal",
				"Fuel oil",
				"Gas",
				"Light oil",
				"Wood",
				"Electricity",
				"Heat",
				"Cooling"
		), ordered = TRUE
)

# Consumption in energyBalance is negative, therefore multiply by -1.
#fuelUse <- energyBalance * temperdays * -1 * 3600 * 24 # From MW -> MJ /a. # This is used for continuous energy use.
#fuelUse <- oapply(fuelUse, cols = c("Temperature", "Burner", "Consumable"), FUN = sum)
#emissions <- EvalOutput(emissions)
# Calculate exposure based on average iF.
#exposure <- EvalOutput(exposure) 
#exposure@output <- exposure@output[exposure@output$Area == "Average" , ]
#exposure <- oapply(exposure, cols = c("Plant", "Emission_site", "Emission_height", "Area"), FUN = sum)
#totcases <- EvalOutput(totcases)
#DALYs <- EvalOutput(DALYs)

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

eb <-EnergyNetworkOptim[EnergyNetworkOptim$Process_variable_type == "Activity",]
colnames(eb@output)[colnames(eb@output) == "Process_variable_name"] <- "Plant"
eb$Process_variable_type <- NULL

ebtemp <- eb
ebtemp@output <- subset(eb@output, Time %in% c("2035") & EnergySavingPolicy == "BAU" & PlantPolicy == "BAU")
ebtemp <- truncateIndex(ebtemp, cols = "Plant", bins = 7)

ggplot(
				ebtemp@output, 
				aes(x = Plant, weight = EnergyNetworkOptimResult, fill = Plant)
		) + geom_bar() +
		facet_wrap( ~ Temperature) + theme_gray(base_size = BS) +
		theme(axis.text.x = element_blank()) + # Turn text and adjust to right
		labs(
				title = "Power plant activity by temperature daily optim \nPlant policy = BAU, Year = 2035",
				#x = "Plant",
				y = "Power output daily average (MW)"
		)

###### oggplot does not work as it should
#oggplot(ebtemp, "Plant", fill = "Plant", turnx = TRUE,
#	labs.title = "Power plant activity by temperature daily optim \nPlant policy = BAU, Year = 2035",
#	labs.y = "Power output daily average (MW)"
#) + facet_wrap( ~ Temperature)

ebtemp <- eb
ebtemp@output <- subset(eb@output, Time %in% c("2035") & EnergySavingPolicy == "BAU" & PlantPolicy == "Custom")
ebtemp <- truncateIndex(ebtemp, cols = "Plant", bins = 7)

ggplot(ebtemp@output, 
				aes(x = Plant, weight = EnergyNetworkOptimResult, fill = Plant))+geom_bar()+
		facet_wrap( ~ Temperature) + theme_gray(base_size = BS) +
		theme(axis.text.x = element_blank()) + # Turn text and adjust to right
		labs(
				title = "Power plant activity by temperature daily optim \nPlant policy = Custom, Year = 2035",
				#x = "Plant",
				y = "Power output daily average (MW)"
		)


ebtemp <- eb
ebtemp@output <- subset(eb@output, Time %in% c("2005") & EnergySavingPolicy == "BAU" & PlantPolicy == "BAU" & Temperature == "(0,3]")
ebtemp <- truncateIndex(ebtemp, cols = "Plant", bins = 7)

ggplot(ebtemp@output, 
				aes(x = Plant, weight = EnergyNetworkOptimResult, fill = Plant))+geom_bar()+
		facet_wrap( ~ Temperature) + theme_gray(base_size = BS) +
		theme(axis.text.x = element_blank()) + # Turn text and adjust to right
		labs(
				title = "Power plant activity by temperature daily optim \nPlant policy = BAU, Year = 2005",
				#x = "Plant",
				y = "Power output daily average (MW)"
		)

#oggplot(subset(ebtemp@output, Temperature == "(0,3]"), x = "Plant", fill = "Plant",
#labs.title = "Daily activity at temperature 0-3 C", y ="Activity of power plant (MW)")

fu <- fuelUse
#result(fu) <- round(result(fu), 12) # Round results smaller than 1e-12 to 0
fu$Burner <- NULL
#fu <- truncateIndex(fuelUse / 3.6E+6, cols = "Plant", bins = 7) # From MJ/a -> GWh/a
futemp <- fu / 3.6E+6 # From MJ/a -> GWh/a
futemp@output <- subset(fu@output, Time %in% c("2015", "2035", "2065") & PlantPolicy == "BAU")
futemp <- truncateIndex(futemp, cols = "Plant", bins = 7)

ggplot(futemp@output, aes(x = Fuel, fill = Plant)) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult > 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack"
		) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult < 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack"
		) +
		geom_hline(aes(yintercept=0)) + 
		facet_grid(Time ~ EnergySavingPolicy) + theme_gray(base_size = BS) +
		theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + # Turn text and adjust to right
		labs(
				title = "Energy commodity flows \n Plant policy = BAU",
				y = "Total annual energy (GWh/a)"
		)

futemp <- fu / 3.6E+6 # From MJ/a -> GWh/a
futemp@output <- subset(fu@output, Time %in% c("2005") & PlantPolicy == "BAU" & EnergySavingPolicy == "BAU")
futemp <- truncateIndex(futemp, cols = "Plant", bins = 7)

ggplot(futemp@output, aes(x = Fuel, fill = Plant)) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult > 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack"
		) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult < 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack"
		) +
		geom_hline(aes(yintercept=0)) + 
		theme_gray(base_size = BS) +
		theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + # Turn text and adjust to right
		labs(
				title = "Energy commodity flows in 2015 \n Plant policy = BAU, Energy saving policy == BAU",
				y = "Total annual energy (GWh/a)"
		)

#oggplot(subset(ebtemp@output, Temperature == "(0,3]"), x = "Plant", fill = "Plant",
#labs.title = "Daily activity at temperature 0-3 C", y ="Total annual energy 2015 (GWh/a)")

futemp <- fu / 3.6E+6 # From MJ/a -> GWh/a
futemp@output <- subset(fu@output, Time %in% c("2015", "2035", "2065") & PlantPolicy == "Custom")
futemp <- truncateIndex(futemp, cols = "Plant", bins = 7)

ggplot(futemp@output, aes(x = Fuel, fill = Plant)) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult > 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack"
		) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult < 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack"
		) +
		geom_hline(aes(yintercept=0)) + 
		facet_grid(Time ~ EnergySavingPolicy) + theme_gray(base_size = BS) +
		theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + # Turn text and adjust to right
		labs(
				title = "Energy commodity flows \n Plant policy = Custom",
				y = "Total annual energy (GWh/a)"
		)

futemp <- fu * 1E-6 # From MJ/a -> TJ/a
futemp@output <- subset(fu@output, Fuel %in% c("Heat"))
futemp <- truncateIndex(futemp, cols = "Plant", bins = 10)
futemp@output$Time <- as.numeric(as.character(futemp@output$Time))

ggplot(futemp@output, aes(x = Time, fill = Plant)) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult > 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack", 
				binwidth = 5
		) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult < 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack", 
				binwidth = 5
		) +
		geom_hline(aes(yintercept=0)) + #geom_bar(binwidth = 5)+
		facet_grid(EnergySavingPolicy ~ PlantPolicy) + theme_gray(base_size = BS) +
		labs(
				title = "District heat flow",
				#x = "Time",
				y = "Total annual energy (TJ/a)"
		)

futemp <- fu * 1E-6 # From MJ/a -> TJ/a
futemp@output <- subset(fu@output, Fuel %in% c("Electricity"))
futemp <- truncateIndex(futemp, cols = "Plant", bins = 10)
futemp@output$Time <- as.numeric(as.character(futemp@output$Time))

ggplot(futemp@output, aes(x = Time, fill = Plant)) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult > 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack", 
				binwidth = 5
		) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult < 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack", 
				binwidth = 5
		) +
		geom_hline(aes(yintercept=0)) + #geom_bar(binwidth = 5)+
		facet_grid(EnergySavingPolicy ~ PlantPolicy) + theme_gray(base_size = BS) +
		labs(
				title = "Electricity flow",
				#x = "Time",
				y = "Total annual energy (TJ/a)"
		)

futemp <- fu * 1E-6 # From MJ/a -> TJ/a
futemp@output <- subset(fu@output, Fuel %in% c("Cooling"))
futemp <- truncateIndex(futemp, cols = "Plant", bins = 10)
futemp@output$Time <- as.numeric(as.character(futemp@output$Time))

ggplot(futemp@output, aes(x = Time, fill = Plant)) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult > 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack", 
				binwidth = 5
		) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult < 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack", 
				binwidth = 5
		) +
		geom_hline(aes(yintercept=0)) + #geom_bar(binwidth = 5)+
		facet_grid(EnergySavingPolicy ~ PlantPolicy) + theme_gray(base_size = BS) +
		labs(
				title = "Cooling flow",
				#x = "Time",
				y = "Total annual energy (TJ/a)"
		)

EnergyNetworkCost <- EvalOutput(EnergyNetworkCost)

EnergyNetworkCost@output$Time <- as.numeric(as.character(EnergyNetworkCost@output$Time))
EnergyNetworkCost <- EnergyNetworkCost[EnergyNetworkCost@output$Time >= 2000 , ]

oprint(head(EnergyNetworkCost@output))

ggplot(subset(truncateIndex(EnergyNetworkCost, cols = "Plant", bins = 12)@output, EnergySavingPolicy == "BAU"), 
	aes(x = Time, weight = EnergyNetworkCostResult / 10, fill = Plant))+geom_bar(binwidth = 10, position = "stack") + 
	facet_grid(PlantPolicy ~ . ) + theme_gray(base_size = BS) +
labs(title = "Investment, management and operation costs",
y = "Cost M€/year")

ggplot(subset(truncateIndex(EnergyNetworkCost, cols = "Plant", bins = 12)@output, EnergySavingPolicy == "BAU"), 
	aes(x = PlantPolicy, weight = EnergyNetworkCostResult / 10))+geom_bar(binwidth = 10, position = "stack") + 
	theme_gray(base_size = BS) +
labs(title = "Investment, management and operation costs",
y = "Cost M€/year")

#DALYs <- EvalOutput(DALYs)

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

#ggplot(subset(emis, PlantPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Fuel)) + 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 (Plant policy = BAU)",
#				y = "Emissions (ton /a)"
#		)

## This code is Op_en7237/intermediates on page [[Helsinki energy decision 2015]]
library(OpasnetUtils)
library(ggplot2)
library(rgdal)
library(maptools)
library(RColorBrewer)
library(classInt)
library(RgoogleMaps)

### Technical parameters

openv.setN(0) # use medians instead of whole sampled distributions
BS <- 24 # base_size = font size in graphs
figstofile <- FALSE
saveobjects <- FALSE

########################## Case-specific data and submodels

### Decisions
decisions <- opbase.data("Op_en7237", subset = "Decisions") # [[Helsinki energy decision 2015]]

oprint(decisions)

### Energy production in Helsinki

objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]] findrest
objects.latest("Op_en7311", code_name = "energyProcess") # [[Helsinki energy production]]
objects.latest("Op_en7311", code_name = "plantParameters") # [[Helsinki energy production]]
objects.latest("Op_en4151", code_name = "fuelPrice") # [[Prices of fuels in heat production]]

#objects.latest("Op_en5141", code_name = "energyBalancetest") # [[Energy balance]]
objects.latest("Op_en5141", code_name = "EnergyNetworkOptim") # [[Energy balance]] incl EnergyNetworkCost

### Building data in Helsinki

# Observation years must be defined for an assessment.
obstime <- Ovariable("obstime", data = data.frame(Obsyear = factor(seq(1985, 2065, 10), ordered = TRUE), Result = 1))
heatingShares <- 1 # Shares of different heating types in the building stock (already in the building data)
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)?

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 = "buildingstest") # [[Building model]] # Generic building model.

## Energy use 
#objects.latest("Op_en5488", code_name = "energyUseTemperature")        # [[Energy use of buildings]] temperature version of energyUse
objects.latest("Op_en5488", code_name = "temperene")        # [[Energy use of buildings]] temperature version of energyUse
objects.latest("Op_en5488", code_name = "nontemperene")        # [[Energy use of buildings]] temperature version of energyUse
objects.latest("Op_en7317", code_name = "temperatures") # [[Helsinki energy consumption]]
#objects.latest("Op_en5488", code_name = "energyUseTemperature") # [[Energy use of buildings]]
objects.latest("Op_en5488", code_name = "EnergyConsumerDemand") # [[Energy use of buildings]]

requiredName <- "Heat" # Name of the fuel type that must match for input and output

### Energy and emissions
objects.latest("Op_en7311", code_name = "emissionLocationsPerPlant") # [[Helsinki energy production]] also heatingShares
objects.latest("Op_en2791", code_name = "emissionFactors") # [[Emission factors for burning processes]] 
objects.latest("Op_en2791", code_name = "emissionstest") # [[Emission factors for burning processes]]
objects.latest("Op_en7311", code_name = "fuelShares") # [[Helsinki energy production]]
objects.latest("Op_en7311", code_name = "nondynsupply") # [[Helsinki energy production]]
timelylimit <- 10000 # The largest production (MW) available at each plant (used to describe time-varying limits).

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

#energyUse <- EvalOutput(energyUse)

# energy <- energyUse * 1E-3 / 365.25 / 24 # kWh /a -> MW This is used for annual energy consumption
#energy <- energyUse / 1E+6 # W -> MW # This is used for continuous energy use.

#energyBalance <- EvalOutput(energyBalance)
#EnergyConsumerDemand <- EvalOutput(EnergyConsumerDemand)

fuelUse <- EvalOutput(fuelUse)

fuelUse$Fuel <- factor(
		fuelUse$Fuel, levels = c(
				"Biofuel",
				"Coal",
				"Fuel oil",
				"Gas",
				"Light oil",
				"Wood",
				"Electricity",
				"Heat",
				"Cooling"
		), ordered = TRUE
)

# Consumption in energyBalance is negative, therefore multiply by -1.
#fuelUse <- energyBalance * temperdays * -1 * 3600 * 24 # From MW -> MJ /a. # This is used for continuous energy use.
#fuelUse <- oapply(fuelUse, cols = c("Temperature", "Burner", "Consumable"), FUN = sum)
emissions <- EvalOutput(emissions)
# Calculate exposure based on average iF.
exposure <- EvalOutput(exposure) 
exposure@output <- exposure@output[exposure@output$Area == "Average" , ]
exposure <- oapply(exposure, cols = c("Plant", "Emission_site", "Emission_height", "Area"), FUN = sum)
totcases <- EvalOutput(totcases)
DALYs <- EvalOutput(DALYs)

############################ Output tables and graphs

oprint(head(buildings@output))
#oprint(head(nondynsupply@output))
#oprint(head(energyUse@output))
oprint(head(EnergyConsumerDemand@output))
#oprint(head(energy@output))
#oprint(head(energyBalance@output))
oprint(head(EnergyNetworkOptim@output))
oprint(head(emissions@output))
oprint(head(exposure@output))
oprint(head(totcases@output))
oprint(head(DALYs@output))

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

bui <- oapply(buildings * 1E-6, cols = c("City_area", "buildingsSource"), FUN = sum)
bui <- truncateIndex(bui, cols = "Heating", bins = 4)@output

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(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)"
		)

ggplot(buildings@output, aes(x = Time, weight = buildingsResult, fill = Efficiency))+geom_bar(binwidth = 5)+
		facet_grid(EnergySavingPolicy ~ Renovation)

emissions@output$Time <- as.numeric(as.character(emissions@output$Time))

# Plot energy need and emissions

hea <- EnergyConsumerDemand * temperdays * 24 * 1E-9 # From W -> GWh /a.
hea@output <- hea@output[hea@output$Consumable == "Heating" , ]
hea <- oapply(hea, cols = c("City_area", "buildingsSource"), FUN = sum)
hea <- truncateIndex(hea, cols = "Heating", bins = 4)@output

ggplot(hea, aes(x = Time, weight = Result, fill = Renovation)) + geom_bar(binwidth = 5) + 
		facet_wrap( ~ EnergySavingPolicy) + theme_gray(base_size = BS) +
		labs(
				title = "Energy used in heating in Helsinki",
				x = "Time",
				y = "Heating energy (GWh /a)"
		)

eb <-EnergyNetworkOptim[EnergyNetworkOptim$Process_variable_type == "Activity",]
colnames(eb@output)[colnames(eb@output) == "Process_variable_name"] <- "Plant"
eb$Process_variable_type <- NULL

ebtemp <- eb
ebtemp@output <- subset(eb@output, Time %in% c("2035") & EnergySavingPolicy == "BAU" & PlantPolicy == "BAU")
ebtemp <- truncateIndex(ebtemp, cols = "Plant", bins = 7)

ggplot(
	ebtemp@output, 
	aes(x = Plant, weight = EnergyNetworkOptimResult, fill = Plant)
) + geom_bar() +
	facet_wrap( ~ Temperature) + theme_gray(base_size = BS) +
	theme(axis.text.x = element_blank()) + # Turn text and adjust to right
	labs(
		title = "Power plant activity by temperature daily optim \nPlant policy = BAU, Year = 2035",
		#x = "Plant",
		y = "Power output daily average (MW)"
	)

ebtemp <- eb
ebtemp@output <- subset(eb@output, Time %in% c("2035") & EnergySavingPolicy == "BAU" & PlantPolicy == "Biofuels and process heat")
ebtemp <- truncateIndex(ebtemp, cols = "Plant", bins = 7)

ggplot(ebtemp@output, 
		aes(x = Plant, weight = EnergyNetworkOptimResult, fill = Plant))+geom_bar()+
		facet_wrap( ~ Temperature) + theme_gray(base_size = BS) +
		theme(axis.text.x = element_blank()) + # Turn text and adjust to right
		labs(
				title = "Power plant activity by temperature daily optim \nPlant policy = Biofuels and process heat option, Year = 2035",
				#x = "Plant",
				y = "Power output daily average (MW)"
		)

fu <- fuelUse
#result(fu) <- round(result(fu), 12) # Round results smaller than 1e-12 to 0
fu$Burner <- NULL
#fu <- truncateIndex(fuelUse / 3.6E+6, cols = "Plant", bins = 7) # From MJ/a -> GWh/a
futemp <- fu / 3.6E+6 # From MJ/a -> GWh/a
futemp@output <- subset(fu@output, Time %in% c("2015", "2035", "2065") & PlantPolicy == "BAU")
futemp <- truncateIndex(futemp, cols = "Plant", bins = 7)

ggplot(futemp@output, aes(x = Fuel, fill = Plant)) +
		geom_bar(
			data = subset(futemp@output, fuelUseResult > 0),
			aes(weight = -1 * fuelUseResult), 
			position="stack"
		) +
		geom_bar(
			data = subset(futemp@output, fuelUseResult < 0),
			aes(weight = -1 * fuelUseResult), 
			position="stack"
		) +
		geom_hline(aes(yintercept=0)) + 
		facet_grid(Time ~ EnergySavingPolicy) + theme_gray(base_size = BS) +
		theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + # Turn text and adjust to right
		labs(
				title = "Energy commodity flows \n Plant policy = BAU",
				y = "Total annual energy (GWh/a)"
		)

futemp <- fu / 3.6E+6 # From MJ/a -> GWh/a
futemp@output <- subset(fu@output, Time %in% c("2015", "2035", "2065") & PlantPolicy == "Biofuels and process heat")
futemp <- truncateIndex(futemp, cols = "Plant", bins = 7)

ggplot(futemp@output, aes(x = Fuel, fill = Plant)) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult > 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack"
		) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult < 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack"
		) +
		geom_hline(aes(yintercept=0)) + 
		facet_grid(Time ~ EnergySavingPolicy) + theme_gray(base_size = BS) +
		theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + # Turn text and adjust to right
		labs(
				title = "Energy commodity flows \n Plant policy = Biofuels and process heat",
				y = "Total annual energy (GWh/a)"
		)

futemp <- fu * 1E-6 # From MJ/a -> TJ/a
futemp@output <- subset(fu@output, Fuel %in% c("Heat"))
futemp <- truncateIndex(futemp, cols = "Plant", bins = 10)
futemp@output$Time <- as.numeric(as.character(futemp@output$Time))

ggplot(futemp@output, aes(x = Time, fill = Plant)) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult > 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack", 
				binwidth = 5
		) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult < 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack", 
				binwidth = 5
		) +
		geom_hline(aes(yintercept=0)) + #geom_bar(binwidth = 5)+
		facet_grid(EnergySavingPolicy ~ PlantPolicy) + theme_gray(base_size = BS) +
		labs(
				title = "District heat flow",
				#x = "Time",
				y = "Total annual energy (TJ/a)"
		)

futemp <- fu * 1E-6 # From MJ/a -> TJ/a
futemp@output <- subset(fu@output, Fuel %in% c("Electricity"))
futemp <- truncateIndex(futemp, cols = "Plant", bins = 10)
futemp@output$Time <- as.numeric(as.character(futemp@output$Time))

ggplot(futemp@output, aes(x = Time, fill = Plant)) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult > 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack", 
				binwidth = 5
		) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult < 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack", 
				binwidth = 5
		) +
		geom_hline(aes(yintercept=0)) + #geom_bar(binwidth = 5)+
		facet_grid(EnergySavingPolicy ~ PlantPolicy) + theme_gray(base_size = BS) +
		labs(
				title = "Electricity flow",
				#x = "Time",
				y = "Total annual energy (TJ/a)"
		)

futemp <- fu * 1E-6 # From MJ/a -> TJ/a
futemp@output <- subset(fu@output, Fuel %in% c("Cooling"))
futemp <- truncateIndex(futemp, cols = "Plant", bins = 10)
futemp@output$Time <- as.numeric(as.character(futemp@output$Time))

ggplot(futemp@output, aes(x = Time, fill = Plant)) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult > 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack", 
				binwidth = 5
		) +
		geom_bar(
				data = subset(futemp@output, fuelUseResult < 0),
				aes(weight = -1 * fuelUseResult), 
				position="stack", 
				binwidth = 5
		) +
		geom_hline(aes(yintercept=0)) + #geom_bar(binwidth = 5)+
		facet_grid(EnergySavingPolicy ~ PlantPolicy) + theme_gray(base_size = BS) +
		labs(
				title = "Cooling flow",
				#x = "Time",
				y = "Total annual energy (TJ/a)"
		)


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

ggplot(subset(emis, PlantPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Fuel)) + 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 (Plant policy = BAU)",
				y = "Emissions (ton /a)"
		)

ggplot(subset(emis, PlantPolicy == "Biofuels and process heat"), 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 (Plant policy = Biofuels and process heat)",
				y = "Emissions (ton /a)"
		)

ggplot(subset(emis, EnergySavingPolicy == "BAU"), aes(x = Time, weight = emissionsResult, fill = Fuel)) + geom_bar(binwidth = 5) +
		facet_grid(Pollutant ~ PlantPolicy, scale = "free_y") + theme_gray(base_size = BS) +
		labs(
				title = "Emissions from heating in Helsinki",
				x = "Time (Energy saving policy = BAU)",
				y = "Emissions (ton /a)"
		)

ggplot(exposure@output, aes(x = Time, weight = exposureResult, fill = Fuel)) + geom_bar(binwidth = 5) + 	
		facet_grid(EnergySavingPolicy ~ Exposure_agent) + 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 ~ .) + 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, PlantPolicy == "BAU"), aes(x = Time, weight = totcasesResult, fill = Fuel))+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 (Plant policy = BAU)",
				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("2015", "2035") & Response == "Total mortality"
)

oprint(
		oapply(temp, INDEX = c("Time", "EnergySavingPolicy", "PlantPolicy"), 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 = Fuel))+geom_bar(binwidth = 5) + 
		facet_grid(EnergySavingPolicy ~ PlantPolicy) +
		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 == 2035), aes(x = EnergySavingPolicy, weight = DALYsResult, fill = Response))+geom_bar() + 
		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)"
		)

if(saveobjects) {
	
	# Clean decisions and previous results not wanted/needed by the half-model
	
	energyProcess@output <- data.frame()
	plantParameters@output <- data.frame()
	EnergyNetworkOptim@output <- data.frame()
	#EnergyConsumerDemand@output <- data.frame()
	
	totcases@output <- data.frame()
	DALYs@output <- data.frame()
	exposure@output <- data.frame()
	
	objects.store(list = ls()) # Save all objects of the global environment.
	cat("All objects of the global environment are stored.\n")

	# Super short version
#	objects.store(
#			BS, buildings, decisions, EnergyConsumerDemand,
#			EnergyConsumerDemandTotal, EnergyFlowHeat, EnergyNetworkDemand, EnergyNetworkOptim, 
#			energyProcess, fuelPrice, fuelUse, plantParameters, temperdays, truncateIndex, EnergyNetworkCost
#	)
	
	#objects.store(
	#		bgexposure, BS, buildings, BW, changeBuildings, DALYs, decisions, disincidence,
	#		dose, dummy, DW, emissionFactors, emissionLocations, emissions, EnergyConsumerDemand,
	#		EnergyConsumerDemandTotal, EnergyFlowHeat, EnergyFlowOther, EnergyNetworkDemand, EnergyNetworkOptim, 
	#		energyProcess, ERF, exposure, findrest, forgetDecisions,
	#		frexposed, fuelPrice, fuelShares, fuelSharesgeneric, fuelUse, heatingShares, iF,
	#		L, nondynsupply, nontemperene, obstime, plantParameters, population, 
	#		renovationRate, renovationRatio, renovationShares, requiredName, RR, stockBuildings,
	#		temperatures, temperdays, temperene, testforrow, threshold, timelylimit,
	#		totcases, truncateIndex
	#)
	#cat("
	#				bgexposure, BS, buildings, BW, changeBuildings, DALYs, decisions, disincidence,
	#		dose, dummy, DW, emissionFactors, emissionLocations, emissions, EnergyConsumerDemand,
	#		EnergyConsumerDemandTotal, EnergyFlowHeat, EnergyFlowOther, EnergyNetworkDemand, EnergyNetworkOptim, 
	#		energyProcess, energyUse, ERF, exposure, findrest, forgetDecisions,
	#		frexposed, fuelPrice, fuelShares, fuelSharesgeneric, fuelUse, heatingShares, iF,
	#		L, nondynsupply, nontemperene, obstime, plantParameters, population, 
	#		renovationRate, renovationRatio, renovationShares, requiredName, RR, stockBuildings,
	#		temperatures, temperdays, temperene, testforrow, threshold, timelylimit,
	#		totcases, truncateIndex.\n")
}

koord <- opbase.data("Op_en7115", subset = "Locations of city areas")
koord <- tidy(koord, widecol = "Location")
colnames(koord)[colnames(koord) == "City_area"] <- "Emission_site"
koord$Result <- 1
koord <- Ovariable("koord", output = koord, marginal = colnames(koord) != "Result")

emis <- emissions * koord

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

EnergyNetworkCost <- EvalOutput(EnergyNetworkCost)

EnergyNetworkCost@output$Time <- as.numeric(as.character(EnergyNetworkCost@output$Time))
EnergyNetworkCost <- EnergyNetworkCost[EnergyNetworkCost@output$Time >= 2000 , ]

ggplot(subset(truncateIndex(EnergyNetworkCost, cols = "Plant", bins = 12)@output, EnergySavingPolicy == "BAU"), 
	aes(x = Time, weight = EnergyNetworkCostResult / 10, fill = Plant))+geom_bar(binwidth = 10, position = "stack") + 
	facet_grid(PlantPolicy ~ . ) +
labs(title = "Investment, management and operation costs",
y = "Cost M€/year")

# Tähän pitäisi laittaa myös kokonaiskustannus diskontattuna. 

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.