Skip to content
Snippets Groups Projects
Commit b0e9aab2 authored by R0slyn's avatar R0slyn
Browse files

Changes to yield shock generator file.

parent c8ac4ebf
No related branches found
No related tags found
No related merge requests found
require(data.table) require(data.table)
require(raster) require(raster)
generateShocksTable = function(dt, startYear, endYear, shockMultiplier,makeMapFiles){ generateShocksTable = function(dt, startYear, endYear, shockFrequencyMultiplier,shockMagnitudeMultiplier, makeMapFiles){
shockTable = NULL shockTable = NULL
cotterellRegions =fread(file.path(baseDataDir,"cotterellRegions.csv"), header=TRUE) cotterellRegions =fread(file.path(baseDataDir,"cotterellRegions.csv"), header=TRUE)
noCountriesTab = data.table(table(cotterellRegions[,list(Region)]))[,list(Region=V1,nCountries = N)] noCountriesTab = data.table(table(cotterellRegions[,list(Region)]))[,list(Region=V1,nCountries = N)]
...@@ -11,7 +11,7 @@ generateShocksTable = function(dt, startYear, endYear, shockMultiplier,makeMapFi ...@@ -11,7 +11,7 @@ generateShocksTable = function(dt, startYear, endYear, shockMultiplier,makeMapFi
ratesTab[,shockRate := N/nCountries] ratesTab[,shockRate := N/nCountries]
#probability of getting a shock rate of value x in year y #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] shockProbs = rbind(shockProbs,shockProbs[,list(shockRate=0,prob=1-sum(prob)), by=c('Region','Sector')])[prob != 0]
for(cropType in unique(dt$Sector)){ for(cropType in unique(dt$Sector)){
...@@ -31,7 +31,7 @@ generateShocksTable = function(dt, startYear, endYear, shockMultiplier,makeMapFi ...@@ -31,7 +31,7 @@ generateShocksTable = function(dt, startYear, endYear, shockMultiplier,makeMapFi
for(year in va[,unique(Year)]){ for(year in va[,unique(Year)]){
#randomly choose countries within the region to shock and by what magnitude #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)) 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]) countries = sample(cotterellRegions[Region==region, Country],va[Year==year,numberOfCountriesToShock])
shockTable=rbind(shockTable, data.table(year,mapFilename=countries,crop=cropType,shockFactor=shockMagnitudes, region)) 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" ...@@ -72,19 +72,38 @@ baseOutputDir="/exports/csce/eddie/geos/groups/LURG/models/PLUM/output"
ensembleFile = commandArgs(trailingOnly = TRUE)[1] ensembleFile = commandArgs(trailingOnly = TRUE)[1]
scenarioTable = fread(file.path(plumDirectory,ensembleFile)) scenarioTable = fread(file.path(plumDirectory,ensembleFile))
singleYieldShockFile = commandArgs(trailingOnly = TRUE)[2]
print(singleYieldShockFile)
if("YIELD_SHOCKS_MULTIPLIER" %in% colnames(scenarioTable)){ 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){ if (singleYieldShockFile != 1){
print("Yield shock multiplier different within ensemble, using only the first value") 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, if(length(shockFreqMultiplier) > 1 | length(shockMagMultiplier) > 1){
shockMultiplier=yieldShockMultiplier, makeMapFiles=TRUE) 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]){ for(scenario in scenarioTable[Ensemble==ensemble, Scenario]){
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment