R/full_valdata.R

In pik-piam/validation: MAgPIE validation library

#.translate_iso_data <- function(x) {
#  require(moinput)
#  x <- clean_magpie(x)
#  if(all(!(c("SXM","CUW") %in% getRegions(x))) & "ANT"%in%getRegions(x)) {
#    tmp <- x[rep("ANT",2),,]
#    getRegions(tmp) <- c("SXM","CUW")
#    tmp["SXM",,] <- tmp["SXM",,]*1/4
#    tmp["CUW",,] <- tmp["CUW",,]*3/4
#    x <- mbind(tmp,x["ANT",,,invert=TRUE])
#  }
#  #x <- moinput:::.toolISOhistorical(x)
#  x <- moinput:::.toolCountryFill(x, fill = 0)
#}



.write_validation_data <- function(gdx,m,file="validation_data.mif") {
  
  time_horizon <- 1900:2015
  # code to get regional data from gdx
  mreg <- sub(",)",")",sub("(",'(gdx=gdx,level="reg",',m,fixed=TRUE),fixed=TRUE)
  # code to get regional data from validation library
  mhreg <- sub("^([^\\(]*)\\(","getData(\\1,",mreg)
  
  # code to get global data from gdx
  mglo <- sub(",)",")",sub("(",'(gdx=gdx,level="glo",',m,fixed=TRUE),fixed=TRUE)
  # code to get global data from validation library
  mhglo <- sub("^([^\\(]*)\\(","getData(\\1,",mglo)
  
  # write validation data
  
  vdata <- NULL
  
  for(i in 1:length(m)) {
    reg <- eval(parse(text=mhreg[i]))[[1]]$data$historical
    glo <- eval(parse(text=mhglo[i]))[[1]]$data$historical
    sources <- names(glo)
    tmp <- new.magpie(cells_and_regions = c(sort(readGDX(gdx,"i")),"GLO"),years = time_horizon, names = sources,sets = c("region","year","source"))
    tmp <- add_dimension(tmp,dim = 3.1,add = "type",nm = "data")
    tmp <- add_columns(tmp,dim=3.1,addnm = c("lo","up"))
    getSets(tmp,fulldim=FALSE)[3] <- sub(".","_",getSets(tmp,fulldim = FALSE)[3],fixed=TRUE)
    getNames(tmp) <- sub(".","_",getNames(tmp),fixed=TRUE)
    for(j in names(glo)) {
      if(ndata(glo[[j]])==1) {
        getNames(glo[[j]]) <- paste("data",j,sep="_")
      } else {
        glo[[j]] <- collapseNames(glo[[j]])
        getNames(glo[[j]]) <- paste(getNames(glo[[j]]),j,sep="_")
      }
      tmp[getRegions(glo[[j]]),getYears(glo[[j]]),getNames(glo[[j]])] <- glo[[j]]
    }
    for(j in names(reg)) {
      if(ndata(reg[[j]])==1) {
        getNames(reg[[j]]) <- paste("data",j,sep="_")
      } else {
        reg[[j]] <- collapseNames(reg[[j]])
        getNames(reg[[j]]) <- paste(getNames(reg[[j]]),j,sep="_")
      }
      tmp[getRegions(reg[[j]]),getYears(reg[[j]]),getNames(reg[[j]])] <- reg[[j]]
    }
    tmp <- add_dimension(tmp,dim = 3.1,add = "value",nm = names(m)[i])
    vdata <- mbind(vdata,tmp)
  }
  
  #filter non existing entries
  for(i in getNames(vdata)) if(all(is.na(vdata[,,i]))) vdata <- vdata[,,i,invert=TRUE]
  for(i in getYears(vdata, as.integer = TRUE)) if(all(is.na(vdata[,i,]))) vdata <- vdata[,i,,invert=TRUE]
  
  append <- FALSE
  for(i in getNames(vdata,fulldim = TRUE)$type_source) {
    write.report(collapseNames(vdata[,,i],2),file=file,model = i, append=append)
    append <- TRUE
  }
  return(vdata)
}

# write model data that should be validated
.write_model_data <- function(gdx,m,file="magpie_data.mif",scenario="default") {
  # code to get regional data from gdx
  mreg <- sub("(",'(gdx=gdx,level="reg",',m,fixed=TRUE)
  # code to get global data from gdx
  mglo <- sub("(",'(gdx=gdx,level="glo",',m,fixed=TRUE)
  
  t <- readGDX(gdx,"t")
  
  model <- NULL
  
  for(i in 1:length(m)) {
    reg <- eval(parse(text=mreg[i]))
    glo <- eval(parse(text=mglo[i]))
    if(is.null(reg) & is.null(glo)) {
      tmp <- NULL
    } else {
      tmp <- mbind(reg,glo)
      getNames(tmp) <- names(m)[i]
      missing <- t[!(t %in% getYears(tmp))]
      if(length(missing)>0) tmp <- add_columns(tmp,dim = 2.1,addnm = missing)
    }
    model <- mbind(model,tmp)
  }
  write.report(model,file=file,model = "MAgPIE",scenario=scenario)
  return(model)
}


full_valdata <- function(gdx,valfile="validation_data.mif",modelfile="magpie_data.mif") {
  #m <- sapply(mapping,function(x) return(x[[1]]))
  
  m <- c('Land Cover|Cropland (million Ha/yr)' = 'land(types="crop")',
         'Land Cover|Pasture (million Ha/yr)'  = 'land(types="past",siclass="sum")',
         'Emissions|CO2|Land Use (Mt C/yr)' = 'emissions(type="co2_c",cumulative=FALSE)',
         'Emissions|N2O|Land Use (Mt N/yr)' = 'emissions(type="n2o_n",cumulative=FALSE)',
         'Emissions|CH4|Land Use (Mt CH4/yr)' = 'emissions(type="ch4",cumulative=FALSE)',
         'Land Cover|Equipped for irrigation (million Ha/yr)' = 'water_AEI()',
         'Land Cover|Cropland|Irrigated (million Ha/yr)' = 'croparea(crop_aggr=TRUE,water="ir")',
         'Water|Withdrawal|Irrigation (km^3/yr)' = 'water_usage(users="agriculture")',
         'Land Cover|Forest|Unmanaged (million Ha/yr)' = 'land(types="forest",sum=FALSE,siclass="sum")',
         'Land Use|Intensity (Index)' = 'tau(prev_year="y1985")',
         'Price|Agriculture|Non-Energy Crops|Index (Index)' = 'priceIndex(crops=c("kfo", "kli"),chain=FALSE,baseyear="y1995")',
         'Fertilizer Use|Nitrogen (Tg Nr/yr)' = 'nr_fertilizer()',
         'Yield|Cereal (tDM/ha/yr)' = 'yields(crops=c("tece", "maiz", "trce", "rice_pro"),crop_aggr=TRUE,water="sum",unit="DM")')

  crops <- c("Temperate cereals"="tece", "Tropical cereals"="trce", "Maize"="maiz", "Rice"="rice_pro", "Others"="others", "Potatoe"="potato", "Cassava"="cassav_sp", 
             "Pulses"="puls_pro", "Soybean"="soybean", "Rapeseed"="rapeseed", "Groundnut"="groundnut", "Sunflower"="sunflower", 
             "Palm oil"="oilpalm", "Cotton"="cottn_pro", "Sugar beet"="sugr_beet", "Sugar cane"="sugr_cane", "Fodder"="foddr", "Ruminent meat"="livst_rum", 
             "Pig meat"="livst_pig", "Chicken meat"="livst_chick", "Eggs"="livst_egg", "Milk"="livst_milk")

  tmp <- paste0('Production(categories="',crops,'")')
  names(tmp) <- paste0("Production Cost|Agriculture|",names(crops)," (million tDM/yr)")
  m <- c(m,tmp)


  trade_crops <- c("Temperate cereals"="tece", "Tropical cereals"="trce", "Maize"="maiz", "Rice"="rice_pro", "Others"="others", "Potatoe"="potato", "Cassava"="cassav_sp",
                   "Pulses"="puls_pro", "Ruminent meat"="livst_rum", "Pig meat"="livst_pig", "Chicken meat"="livst_chick", "Eggs"="livst_egg", "Milk"="livst_milk")


  tmp <- paste0('trade_flows(commodities="',trade_crops,'")')
  names(tmp) <- paste0("Trade|Exports|Agriculture|Volume|",names(trade_crops)," (million mDT)")
  m <- c(m,tmp)
  
 
  val <- .write_validation_data(gdx,m,file=valfile)
  model <- .write_model_data(gdx,m,file=modelfile)  
  return(list(val=val,model=model))
}