Run R scripts during scenario creation for yield shocks or shocks.

scripts/CreateYieldShockMaps.R

+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]))
+}
+
+
+
+
+
+

scripts/createScenarios.sh

 #!/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