summariseUKlanduse.R

# Script to summarise land use for the UK only, PLUM post processing 
# 
# Author: rhenry2
##########################################

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')),
                   area = get(paste0(cropType, '_A'))*cropland, 
                   yield = get(paste0(cropType, '_Y'))) ]
  
  ##irrigation rate in in l/m2 and fertiliser in kg/ha 
  yDt[thisCropFact>0, list(cropType = cropType, 
                           irrig=sum(area * irrig)/sum(area), 
                           fert=sum(area * fert)/sum(area),
                           otherint=sum(area * otherint)/sum(area),
                           area = sum(area))]
}

summariseLandUseOutput = function(simDir) {
  
  countryCells = data.table(read.table("/exports/csce/eddie/geos/groups/LURG/models/PLUM/plumv2/data/halfdeg/ukCells.txt", header=TRUE))
  years = list.files(simDir, pattern="[0-9]+")
  resDt = NULL
  
  for (year in years) {
    cropSumDt = NULL
    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))
      
      luDt = merge(luDt, countryCells, by=c('Lon','Lat'))
      
      for (cropType in c('maize', 'rice', 'wheat', 'energycrops', 'oilcrops', 'pulses', 'starchyRoots', "fruitveg", "sugar")) {
        cropSummary = cropIntensitySummary(luDt, cropType)
        if (nrow(cropSummary) >0)
          cropSumDt = rbind(cropSumDt, cropSummary)
      }
      #weighted intensity rates by crop type
      intensitySum = cropSumDt[,list(IrrigCrops =sum(irrig*area)/(sum(area)), FertCrops=sum(fert*area)/(sum(area)))]

      landCoverSum = luDt[,list(Cropland=sum(cropland), Pasture=sum(pasture), ManForest = sum(managedForest), 
                                UnmanForest=sum(unmanagedForest),Natural=sum(otherNatural))]
      
      landCoverSum[, AbPasture := luDt[pasture_A > 0 & pasture_FQ == 0 & pasture_IQ == 0 & pasture_OI == 0, sum(pasture)]]
      
      resDt = rbind(cbind(Year=year,intensitySum,landCoverSum), resDt)
      
      
    } else {
      print(paste("Can't find file, so skipping:", fileName))
    }
     
  }
 
  write.table(resDt, file.path(simDir, "ukLC.csv"), sep=",", row.names=FALSE, quote=FALSE)
  resDt
}

baseOutputDir="/exports/csce/eddie/geos/groups/LURG/models/PLUM/output" 

for (directory in list.dirs(full.names = TRUE, recursive = FALSE)){
  for (scenario in basename(list.dirs(path = directory, full.names = TRUE, recursive = FALSE))){
    ensDir=file.path(directory,scenario)
    if (dir.exists(ensDir)) {
      print(paste(ensDir, "exists.  Processing"))
      summariseLandUseOutput(simDir = ensDir)
    } else {
      print(paste(ensDir, "does not exist.  Stopping"))
    }
  }
}