Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
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]))
}