Fixes to yield shock generation.

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)){
 
 
 
+