CreateYieldShockMaps.R

require(data.table)
require(raster)

generateShocksTable = function(dt, startYear, endYear, shockFrequencyMultiplier,shockMagnitudeMultiplier, makeMapFiles){
  shockTable = NULL
  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)]
    
    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')]
  }
  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))
singleYieldShockFile = commandArgs(trailingOnly = TRUE)[2]
print(singleYieldShockFile)

if("SHOCKS_FREQUENCY_MULTIPLIER" %in% colnames(scenarioTable)){
  
  if (singleYieldShockFile == 1){
    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)
    
  }
  
  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 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))) {
        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("No value for shock frequency multiplier in scenario table, not creating yieldShock files")
}