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
require(data.table)
cropIntensitySummary = function (luDt, cropType) {
yDt = luDt[cropland > 0,
list( irrig = get(paste0(cropType, '_IQ')),
fert = get(paste0(cropType, '_FQ')),
otherint = get(paste0(cropType, '_OI')),
thisCropFact = get(paste0(cropType, '_A')),
allCropFact = cropland/area,
area = get(paste0(cropType, '_A'))*cropland,
yield = get(paste0(cropType, '_Y'))) ]
# yDt[thisCropFact>0.05 & allCropFact > 0.03, list(cropType = cropType, yield=sum(area * yield)/sum(area),
yDt[thisCropFact>0, list(cropType = cropType, yield=sum(area * yield)/sum(area),
irrig=sum(area * irrig)/sum(area),
fert=sum(area * fert)/sum(area),
otherint=sum(area * otherint)/sum(area),
area = sum(area)) ]
}
summariseLandUseOutput = function(simDir, writeFile=TRUE) {
years = list.files(simDir, pattern="[0-9]+")
cropSumDt = NULL
for (year in years) {
fileName = file.path(simDir, year, 'LandUse.txt')
if (file.exists(fileName)) {
print (paste("Processing", year, "from", fileName))
luDt = data.table(read.table(fileName, header=TRUE))
for (cropType in c('maize', 'rice', 'wheat', 'energycrops', 'oilcrops', 'pulses', 'starchyRoots')) {
cropSummary = cropIntensitySummary(luDt, cropType)
if (nrow(cropSummary) >0)
cropSumDt = rbind(cropSumDt, cbind(cropSummary, year=year))
}
} else {
print(paste("Can't find file, so skipping:", fileName))
}
}
resDt = rbind(cropSumDt, cropSumDt[cropType != 'energycrops', list(cropType='avgExcEC',
yield=sum(yield*area)/sum(area),
irrig=sum(irrig*area)/sum(area),
fert=sum(fert*area)/sum(area),
otherint=sum(otherint*area)/sum(area),
area=sum(area)), by=year])
resDt = rbind(resDt, cropSumDt[, list(cropType='avgIncEC',
yield=sum(yield*area)/sum(area),
irrig=sum(irrig*area)/sum(area),
fert=sum(fert*area)/sum(area),
otherint=sum(otherint*area)/sum(area),
area=sum(area)), by=year])
if (writeFile) write.table(resDt, file.path(simDir, "LandUseSummary.csv"), sep=",", row.names=FALSE, quote=FALSE)
resDt
}
baseOutputDir="/exports/csce/eddie/geos/groups/LURG/models/PLUM/output" # "~/Downloads"
ensemble = commandArgs(trailingOnly = TRUE)[1]
ensDir=file.path(baseOutputDir, ensemble)
if (dir.exists(ensDir)) {
print(paste(ensDir, "exists. Processing"))
summariseLandUseOutput(simDir = ensDir)
} else {
print(paste(ensDir, "does not exist. Stopping"))
}