From b0e9aab2fc02d112d39867848d065b1d5112a11c Mon Sep 17 00:00:00 2001
From: R0slyn <roslyn.henry.08@aberdeen.ac.uk>
Date: Tue, 6 Aug 2019 12:07:19 +0100
Subject: [PATCH] Changes to yield shock generator file.
---
scripts/CreateYieldShockMaps.R | 39 +++++++++++++++++++++++++---------
1 file changed, 29 insertions(+), 10 deletions(-)
diff --git a/scripts/CreateYieldShockMaps.R b/scripts/CreateYieldShockMaps.R
index a9f12d9f..8df02505 100644
--- a/scripts/CreateYieldShockMaps.R
+++ b/scripts/CreateYieldShockMaps.R
@@ -1,7 +1,7 @@
require(data.table)
require(raster)
-generateShocksTable = function(dt, startYear, endYear, shockMultiplier,makeMapFiles){
+generateShocksTable = function(dt, startYear, endYear, shockFrequencyMultiplier,shockMagnitudeMultiplier, makeMapFiles){
shockTable = NULL
cotterellRegions =fread(file.path(baseDataDir,"cotterellRegions.csv"), header=TRUE)
noCountriesTab = data.table(table(cotterellRegions[,list(Region)]))[,list(Region=V1,nCountries = N)]
@@ -11,7 +11,7 @@ generateShocksTable = function(dt, startYear, endYear, shockMultiplier,makeMapFi
ratesTab[,shockRate := N/nCountries]
#probability of getting a shock rate of value x in year y
- shockProbs = data.table(table(ratesTab[N != 0,list(Region,Sector,shockRate)]))[,list(Region, Sector,shockRate, prob=N*shockMultiplier/53)] #53 as that is number of years in historic time period
+ shockProbs = data.table(table(ratesTab[N != 0,list(Region,Sector,shockRate)]))[,list(Region, Sector,shockRate, prob=N*shockFrequencyMultiplier/53)] #53 as that is number of years in historic time period
shockProbs = rbind(shockProbs,shockProbs[,list(shockRate=0,prob=1-sum(prob)), by=c('Region','Sector')])[prob != 0]
for(cropType in unique(dt$Sector)){
@@ -31,7 +31,7 @@ generateShocksTable = function(dt, startYear, endYear, shockMultiplier,makeMapFi
for(year in va[,unique(Year)]){
#randomly choose countries within the region to shock and by what magnitude
bins=with(sectorHist,sample(length(mids),va[Year==year,numberOfCountriesToShock],p=density,replace=TRUE))
- shockMagnitudes=runif(length(bins),sectorHist$breaks[bins],sectorHist$breaks[bins+1])
+ shockMagnitudes=runif(length(bins),sectorHist$breaks[bins],sectorHist$breaks[bins+1])*shockMagnitudeMultiplier
countries = sample(cotterellRegions[Region==region, Country],va[Year==year,numberOfCountriesToShock])
shockTable=rbind(shockTable, data.table(year,mapFilename=countries,crop=cropType,shockFactor=shockMagnitudes, region))
}
@@ -72,19 +72,38 @@ baseOutputDir="/exports/csce/eddie/geos/groups/LURG/models/PLUM/output"
ensembleFile = commandArgs(trailingOnly = TRUE)[1]
scenarioTable = fread(file.path(plumDirectory,ensembleFile))
+singleYieldShockFile = commandArgs(trailingOnly = TRUE)[2]
+print(singleYieldShockFile)
if("YIELD_SHOCKS_MULTIPLIER" %in% colnames(scenarioTable)){
- for(ensemble in unique(scenarioTable[,Ensemble])){
+ if (singleYieldShockFile == 1){
+ shockFreqMultiplier = unique(scenarioTable[, SHOCKS_FREQUENCY_MULTIPLIER])
+ shockMagMultiplier = unique(scenarioTable[, SHOCKS_MAGNITUDE_MULTIPLIER])
- yieldShockMultiplier = unique(scenarioTable[Ensemble==ensemble, YIELD_SHOCKS_MULTIPLIER])
+ if(length(shockFreqMultiplier) > 1 | length(shockMagMultiplier) > 1){
+ print("A Yield shock multiplier is different in ensemble file, using only the first values")
+ }
+
+ ensembleShockTable=generateShocksTable(dt=fread(file.path(baseDataDir,"historicShocks.csv")), startYear=2010, endYear=2100,
+ shockFrequencyMultiplier=shockFreqMultiplier, shockMagnitudeMultiplier=shockMagMultiplier,makeMapFiles=TRUE)
+
+ }
+
+ for(ensemble in unique(scenarioTable[,Ensemble])){
- if(length(yieldShockMultiplier) > 1){
- print("Yield shock multiplier different within ensemble, using only the first value")
- }
+ if (singleYieldShockFile != 1){
+ shockFreqMultiplier = unique(scenarioTable[Ensemble==ensemble, SHOCKS_FREQUENCY_MULTIPLIER])
+ shockMagMultiplier = unique(scenarioTable[Ensemble==ensemble, SHOCKS_MAGNITUDE_MULTIPLIER])
- ensembleShockTable=generateShocksTable(dt=fread(file.path(baseDataDir,"historicShocks.csv")), startYear=2010, endYear=2100,
- shockMultiplier=yieldShockMultiplier, makeMapFiles=TRUE)
+ if(length(shockFreqMultiplier) > 1 | length(shockMagMultiplier) > 1){
+ print("A Yield shock multiplier is different in ensemble file, using only the first values")
+ }
+
+ ensembleShockTable=generateShocksTable(dt=fread(file.path(baseDataDir,"historicShocks.csv")), startYear=2010, endYear=2100,
+ shockFrequencyMultiplier=shockFreqMultiplier, shockMagnitudeMultiplier=shockMagMultiplier,makeMapFiles=TRUE)
+
+ }
for(scenario in scenarioTable[Ensemble==ensemble, Scenario]){
--
GitLab