From 7501bd24c8daea610ce55a7b37e324248c05f8dc Mon Sep 17 00:00:00 2001
From: R0slyn <roslyn.henry.08@aberdeen.ac.uk>
Date: Thu, 23 May 2019 18:12:55 +0100
Subject: [PATCH] Run R scripts during scenario creation for yield shocks or
shocks.
---
scripts/CreateYieldShockMaps.R | 110 +++++++++++++++++++++++++++++++++
scripts/createScenarios.sh | 25 ++++++--
2 files changed, 131 insertions(+), 4 deletions(-)
create mode 100644 scripts/CreateYieldShockMaps.R
diff --git a/scripts/CreateYieldShockMaps.R b/scripts/CreateYieldShockMaps.R
new file mode 100644
index 00000000..a9f12d9f
--- /dev/null
+++ b/scripts/CreateYieldShockMaps.R
@@ -0,0 +1,110 @@
+require(data.table)
+require(raster)
+
+generateShocksTable = function(dt, startYear, endYear, shockMultiplier,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*shockMultiplier/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])
+ 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')]
+ shockTable
+}
+
+##perhaps better to just keep a list of country maps rather than make them each time in this script
+makeMaps=function(countriesTab, outputDir){
+
+ country_codes = fread(file.path(baseDataDir, "country_codes4.csv"), header=TRUE)
+ boundary = raster(file.path(baseDataDir, "halfdeg", "country_boundaries.asc"))
+ crs(boundary) = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
+
+ for(mapFilename in unique(countriesTab[,mapFilename])){
+ b2 = boundary %in% country_codes[Country %in% mapFilename, numCode]
+ values(b2)[values(b2) == 0] = NA
+ fileN = filename=file.path(baseDataDir, "halfdeg","yieldshockmaps",paste0(mapFilename,'.asc'))
+ if (!file.exists(fileN)){
+ writeRaster(b2, filename=fileN, format='ascii', bylayer=TRUE, overwrite=TRUE)
+ }
+
+ }
+}
+
+#plot(raster("C:\\Users\\rhenry2\\eclipse-workspace\\plumv2\\data\\halfdeg\\yieldshockmaps\\United States of America.asc"))
+
+plumDirectory = "/exports/csce/eddie/geos/groups/LURG/models/PLUM"
+baseDataDir = "/exports/csce/eddie/geos/groups/LURG/models/PLUM/roslyns_build_area/plumv2/data"
+baseOutputDir="/exports/csce/eddie/geos/groups/LURG/models/PLUM/output"
+
+ensembleFile = commandArgs(trailingOnly = TRUE)[1]
+scenarioTable = fread(file.path(plumDirectory,ensembleFile))
+
+if("YIELD_SHOCKS_MULTIPLIER" %in% colnames(scenarioTable)){
+
+ for(ensemble in unique(scenarioTable[,Ensemble])){
+
+ yieldShockMultiplier = unique(scenarioTable[Ensemble==ensemble, YIELD_SHOCKS_MULTIPLIER])
+
+ if(length(yieldShockMultiplier) > 1){
+ print("Yield shock multiplier different within ensemble, using only the first value")
+ }
+
+ ensembleShockTable=generateShocksTable(dt=fread(file.path(baseDataDir,"historicShocks.csv")), startYear=2010, endYear=2100,
+ shockMultiplier=yieldShockMultiplier, 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)
+ }
+ 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]))
+}
+
+
+
+
+
+
diff --git a/scripts/createScenarios.sh b/scripts/createScenarios.sh
index 5f859d1d..711b946e 100755
--- a/scripts/createScenarios.sh
+++ b/scripts/createScenarios.sh
@@ -1,8 +1,13 @@
#!/bin/sh
+#$ -o /exports/csce/eddie/geos/groups/LURG/models/PLUM/output/qlogs
+#$ -e /exports/csce/eddie/geos/groups/LURG/models/PLUM/output/qlogs
+#$ -l h_vmem=8G
+#$ -R y
+
-filename="$1"
common_prop_file=/exports/csce/eddie/geos/groups/LURG/models/PLUM/output/common_properties
output_dir=/exports/csce/eddie/geos/groups/LURG/models/PLUM/output
+filename="/exports/csce/eddie/geos/groups/LURG/models/PLUM/$1"
if [[ $* == *-O* ]]; then overwrite=1
else overwrite=0
@@ -51,8 +56,20 @@ do
done < $filename
-echo "Running R shock script"
- module load R
- R < /exports/csce/eddie/geos/groups/LURG/models/PLUM/plumv2/scripts/createShockFiles.R --no-save --args $1
+source /etc/profile
+if [[ $* == *-ys* ]]; then
+ echo "Running R yield shock script"
+ module load R
+ R < /exports/csce/eddie/geos/groups/LURG/models/PLUM/plumv2/scripts/CreateYieldShockMaps.R --no-save --args $1
+else
+ echo "Not running R yield shock script"
+fi
+if [[ $* == *-s* ]]; then
+ echo "Running R shock script"
+ module load R
+ R < /exports/csce/eddie/geos/groups/LURG/models/PLUM/plumv2/scripts/createShockFiles.R --no-save --args $1
+else
+ echo "Not running R shock script"
+fi
\ No newline at end of file
--
GitLab