From 1f824d717daaba30530ac08a58613737cc4d466d Mon Sep 17 00:00:00 2001
From: R0slyn <roslyn.henry.08@aberdeen.ac.uk>
Date: Tue, 6 Aug 2019 15:16:45 +0100
Subject: [PATCH] Fixes to yield shock generation.
---
scripts/CreateYieldShockMaps.R | 138 ++++++++++++++++++---------------
1 file changed, 74 insertions(+), 64 deletions(-)
diff --git a/scripts/CreateYieldShockMaps.R b/scripts/CreateYieldShockMaps.R
index 8df02505..0cc80be2 100644
--- a/scripts/CreateYieldShockMaps.R
+++ b/scripts/CreateYieldShockMaps.R
@@ -3,46 +3,48 @@ require(raster)
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)]
-
- ratesTab = merge(data.table(table(dt[,list(Year,Region,Sector)])),noCountriesTab, by=c('Region'))
- #Shock rate is number of shocks occurring in a region in a given year, shock rate scaled according to number of countries in a region
- 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*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)){
- #magnitudes of shocks take from global distribution per sector at a country level
- sectorHist = hist(dt[Sector==cropType,PropMagnitude],freq=FALSE)
+ if(shockFrequencyMultiplier != 0 & shockMagnitudeMultiplier != 0){
+ cotterellRegions =fread(file.path(baseDataDir,"cotterellRegions.csv"), header=TRUE)
+ noCountriesTab = data.table(table(cotterellRegions[,list(Region)]))[,list(Region=V1,nCountries = N)]
- for(region in unique(dt$Region)){
- #calculate years where shocks could happen using probability of shock occurring
- va = data.table(Year=seq(startYear,endYear,1),
- shockRate=sample(shockProbs[Sector==cropType & Region==region,as.numeric(shockRate)],size=endYear-startYear+1, replace=TRUE,
- prob=shockProbs[Sector==cropType & Region==region,prob]))[shockRate > 0]
- if(nrow(va)>0){
- #calculate for above years how many countries to shock
- va[,numberOfCountriesToShock := sum(rbinom(noCountriesTab[Region==region,nCountries],size=1,prob=shockRate))]
- va=va[numberOfCountriesToShock>0]
-
- 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])*shockMagnitudeMultiplier
- countries = sample(cotterellRegions[Region==region, Country],va[Year==year,numberOfCountriesToShock])
- shockTable=rbind(shockTable, data.table(year,mapFilename=countries,crop=cropType,shockFactor=shockMagnitudes, region))
- }
- }
+ ratesTab = merge(data.table(table(dt[,list(Year,Region,Sector)])),noCountriesTab, by=c('Region'))
+ #Shock rate is number of shocks occurring in a region in a given year, shock rate scaled according to number of countries in a region
+ 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*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)){
+ #magnitudes of shocks take from global distribution per sector at a country level
+ sectorHist = hist(dt[Sector==cropType,PropMagnitude],freq=FALSE)
+
+ for(region in unique(dt$Region)){
+ #calculate years where shocks could happen using probability of shock occurring
+ va = data.table(Year=seq(startYear,endYear,1),
+ shockRate=sample(shockProbs[Sector==cropType & Region==region,as.numeric(shockRate)],size=endYear-startYear+1, replace=TRUE,
+ prob=shockProbs[Sector==cropType & Region==region,prob]))[shockRate > 0]
+ if(nrow(va)>0){
+ #calculate for above years how many countries to shock
+ va[,numberOfCountriesToShock := sum(rbinom(noCountriesTab[Region==region,nCountries],size=1,prob=shockRate))]
+ va=va[numberOfCountriesToShock>0]
+
+ 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])*shockMagnitudeMultiplier
+ countries = sample(cotterellRegions[Region==region, Country],va[Year==year,numberOfCountriesToShock])
+ shockTable=rbind(shockTable, data.table(year,mapFilename=countries,crop=cropType,shockFactor=shockMagnitudes, region))
+ }
+ }
+ }
}
+
+ if (makeMapFiles){
+ makeMaps(shockTable, outputDir)
+ }
+ shockTable[,mapFilename := paste0(mapFilename,'.asc')]
}
-
- if (makeMapFiles){
- makeMaps(shockTable, outputDir)
- }
- shockTable[,mapFilename := paste0(mapFilename,'.asc')]
shockTable
}
@@ -75,51 +77,58 @@ scenarioTable = fread(file.path(plumDirectory,ensembleFile))
singleYieldShockFile = commandArgs(trailingOnly = TRUE)[2]
print(singleYieldShockFile)
-if("YIELD_SHOCKS_MULTIPLIER" %in% colnames(scenarioTable)){
+if("SHOCKS_FREQUENCY_MULTIPLIER" %in% colnames(scenarioTable)){
if (singleYieldShockFile == 1){
- shockFreqMultiplier = unique(scenarioTable[, SHOCKS_FREQUENCY_MULTIPLIER])
- shockMagMultiplier = unique(scenarioTable[, SHOCKS_MAGNITUDE_MULTIPLIER])
+ shockFreqMultiplier = unique(scenarioTable[, SHOCKS_FREQUENCY_MULTIPLIER])
+ shockMagMultiplier = unique(scenarioTable[, SHOCKS_MAGNITUDE_MULTIPLIER])
+
+
+ if(length(shockFreqMultiplier) > 1 | length(shockMagMultiplier) > 1){
+ print("A Yield shock multiplier parameter has more than one value in scenario table, using only the first values")
+ }
+
+ ensembleShockTable=generateShocksTable(dt=fread(file.path(baseDataDir,"historicShocks.csv")), startYear=2010, endYear=2100,
+ shockFrequencyMultiplier=shockFreqMultiplier, shockMagnitudeMultiplier=shockMagMultiplier,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(ensemble in unique(scenarioTable[,Ensemble])){
- if (singleYieldShockFile != 1){
- shockFreqMultiplier = unique(scenarioTable[Ensemble==ensemble, SHOCKS_FREQUENCY_MULTIPLIER])
- shockMagMultiplier = unique(scenarioTable[Ensemble==ensemble, SHOCKS_MAGNITUDE_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)
-
- }
+ if (singleYieldShockFile != 1){
+ shockFreqMultiplier = unique(scenarioTable[Ensemble==ensemble, SHOCKS_FREQUENCY_MULTIPLIER])
+ shockMagMultiplier = unique(scenarioTable[Ensemble==ensemble, SHOCKS_MAGNITUDE_MULTIPLIER])
+
+ if(length(shockFreqMultiplier) > 1 | length(shockMagMultiplier) > 1){
+ print("A Yield shock multiplier parameter is different within a single ensemble, 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]){
ensDir=file.path(ensemble,scenario)
if (dir.exists(file.path(baseOutputDir,ensDir))) {
- print(paste(ensDir, "exists. Creating yield shock file"))
- write.csv(ensembleShockTable,file.path(baseOutputDir,ensDir,"yieldshocks.csv"), row.names=FALSE, quote=FALSE)
+ if(!is.null(ensembleShockTable)){
+ print(paste(ensDir, "exists. Creating yield shock file"))
+ write.csv(ensembleShockTable,file.path(baseOutputDir,ensDir,"yieldshocks.csv"), row.names=FALSE, quote=FALSE)
+ }
+ else{
+ print(paste0("SHOCKS_FREQUENCY_MULTIPLIER or SHOCKS_MAGNITUDE_MULTIPLIER is set to zero, not creating yield shock files for ", ensemble))
+ }
}
else{
print(paste(ensDir, "does not exist. Stopping"))
}
- }
+
+
+ }
}
-}else{
- print(paste("No value for yield shock multiplier, not creating yieldShocks.csv file for ",scenarioTable[i,Ensemble],scenarioTable[i,Scenario]))
+} else{
+ print("No value for shock frequency multiplier in scenario table, not creating yieldShock files")
}
@@ -127,3 +136,4 @@ if("YIELD_SHOCKS_MULTIPLIER" %in% colnames(scenarioTable)){
+
--
GitLab