# load the required pacakges
library(hector)
library(dplyr)
library(tidyr)
Because Hector models both SO2 direct and indirect effects we canโt look at the Hector code to get our answer. But we will nedd to look at the ratio of the RF / emissions from Hector output.
# Run hector for rcp45
ini <- system.file('input/hector_rcp45.ini', package = 'hector')
core <- newcore(ini)
run(core)
Hector core: unnamed hector core
Start date: 1745
End date: 2300
Current date: 2300
Input file: /Users/dorh012/Library/R/3.5/library/hector/input/hector_rcp45.ini
rslts <- fetchvars(core, dates = 1990:2100, vars = c(RF_SO2(), RF_OC(), EMISSIONS_OC(), EMISSIONS_SO2()))
rslts %>%
mutate(type = if_else(grepl('emissions', variable), 'emissions', 'RF'),
agent = if_else(grepl('OC', variable), 'OC', 'SO2')) %>%
select(year, value, type, agent) %>%
spread(type, value) %>%
# Divide the RF (W/m2) by the emisisons, if the emissions are in different units you are going to
# need to convert them.
mutate(value = RF/emissions) %>%
group_by(agent) %>%
summarise(value = mean(value))
Check out the units on the emissions to see if you need to convert them.
# What are the units on the OC emissions?
getunits(EMISSIONS_OC())
[1] "Tg"
# What are the units on the SO2 emissions?
getunits(EMISSIONS_SO2())
[1] "Gg S"
However if we look at the units on the emissions we see that the SO2 emissions are in Gg S so that means you are going to have to convert the So2 value from W/m/Gg S to W/m/Tg/.
LS0tCnRpdGxlOiAiT0MgdnMgQ08yIGZvciBKb24iCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCgpgYGB7cn0KIyBsb2FkIHRoZSByZXF1aXJlZCBwYWNha2dlcwpsaWJyYXJ5KGhlY3RvcikKbGlicmFyeShkcGx5cikKbGlicmFyeSh0aWR5cikKYGBgCgoKQmVjYXVzZSBIZWN0b3IgbW9kZWxzIGJvdGggU08yIGRpcmVjdCBhbmQgaW5kaXJlY3QgZWZmZWN0cyB3ZSBjYW4ndCBsb29rIGF0IHRoZSBIZWN0b3IgY29kZSB0byBnZXQgb3VyIGFuc3dlci4gQnV0IHdlIHdpbGwgbmVkZCB0byBsb29rIGF0IHRoZSByYXRpbyBvZiB0aGUgUkYgLyBlbWlzc2lvbnMgZnJvbSBIZWN0b3Igb3V0cHV0LiAKCmBgYHtyfQojIFJ1biBoZWN0b3IgZm9yIHJjcDQ1CmluaSA8LSBzeXN0ZW0uZmlsZSgnaW5wdXQvaGVjdG9yX3JjcDQ1LmluaScsIHBhY2thZ2UgPSAnaGVjdG9yJykKY29yZSA8LSBuZXdjb3JlKGluaSkKcnVuKGNvcmUpCnJzbHRzIDwtIGZldGNodmFycyhjb3JlLCBkYXRlcyA9IDE5OTA6MjEwMCwgdmFycyA9IGMoUkZfU08yKCksIFJGX09DKCksIEVNSVNTSU9OU19PQygpLCBFTUlTU0lPTlNfU08yKCkpKQpgYGAKCmBgYHtyfQpyc2x0cyAlPiUgIAogICAgbXV0YXRlKHR5cGUgPSBpZl9lbHNlKGdyZXBsKCdlbWlzc2lvbnMnLCB2YXJpYWJsZSksICdlbWlzc2lvbnMnLCAnUkYnKSwgCiAgICAgICAgICAgYWdlbnQgPSBpZl9lbHNlKGdyZXBsKCdPQycsIHZhcmlhYmxlKSwgJ09DJywgJ1NPMicpKSAlPiUgCiAgICBzZWxlY3QoeWVhciwgdmFsdWUsIHR5cGUsIGFnZW50KSAlPiUgIAogICAgc3ByZWFkKHR5cGUsIHZhbHVlKSAlPiUgCiAgICAjIERpdmlkZSB0aGUgUkYgKFcvbTIpIGJ5IHRoZSBlbWlzaXNvbnMsIGlmIHRoZSBlbWlzc2lvbnMgYXJlIGluIGRpZmZlcmVudCB1bml0cyB5b3UgYXJlIGdvaW5nIHRvIAogICAgIyBuZWVkIHRvIGNvbnZlcnQgdGhlbS4gCiAgICBtdXRhdGUodmFsdWUgPSBSRi9lbWlzc2lvbnMpICU+JSAKICAgIGdyb3VwX2J5KGFnZW50KSAlPiUgCiAgICBzdW1tYXJpc2UodmFsdWUgPSBtZWFuKHZhbHVlKSkKYGBgCgpDaGVjayBvdXQgdGhlIHVuaXRzIG9uIHRoZSBlbWlzc2lvbnMgdG8gc2VlIGlmIHlvdSBuZWVkIHRvIGNvbnZlcnQgdGhlbS4gCgpgYGB7cn0KIyBXaGF0IGFyZSB0aGUgdW5pdHMgb24gdGhlIE9DIGVtaXNzaW9ucz8KZ2V0dW5pdHMoRU1JU1NJT05TX09DKCkpCgojIFdoYXQgYXJlIHRoZSB1bml0cyBvbiB0aGUgU08yIGVtaXNzaW9ucz8gCmdldHVuaXRzKEVNSVNTSU9OU19TTzIoKSkKYGBgCgpIb3dldmVyIGlmIHdlIGxvb2sgYXQgdGhlIHVuaXRzIG9uIHRoZSBlbWlzc2lvbnMgd2Ugc2VlIHRoYXQgdGhlIFNPMiBlbWlzc2lvbnMgYXJlIGluIEdnIFMgc28gdGhhdCBtZWFucyB5b3UgYXJlIGdvaW5nIHRvIGhhdmUgdG8gY29udmVydCB0aGUgU28yIHZhbHVlIGZyb20gVy9tL0dnIFMgdG8gVy9tL1RnLy4gCgo=