Commit ff628d02 authored by Kristine Karstens's avatar Kristine Karstens
Browse files

remove old preprocessing + set cachetype to def

parent 09fd34d0
# (C) 2008-2016 Potsdam Institute for Climate Impact Research (PIK),
# authors, and contributors see AUTHORS file
# This file is part of MAgPIE and licensed under GNU AGPL Version 3
# or later. See LICENSE file or go to http://www.gnu.org/licenses/
# Contact: magpie@pik-potsdam.de
# *********************************************************************
# *** This script combines magpie inputs from a high ***
# *** resolution and a low resolution dataset ***
# *********************************************************************
# Version 1.00
library(landuse)
library(magpie)
############################# BASIC CONFIGURATION #######################################
if(!exists("source_include")) {
input_spam_file<- NULL
output_spam_file <-NULL
folder <- NULL
hcells <- as.vector(as.matrix(read.csv("MAGPIE_Zellnummern_RSF+Rand.csv",head=TRUE,sep=","))) # All the cells, that shall be represented in high resolution
readArgs('input_spam_file','output_spam_file','folder','hcells')
}
#########################################################################################
# read in the spam file:
spam<-as.matrix(read.spam(path(folder,input_spam_file)))
# Affected clusters and relation to cells
cluster_high_res_cells<-c()
for(i in 1:length(hcells)){
cluster_high_res_cells<-c(cluster_high_res_cells,which(spam[,hcells[i]]>0))
}
cell_info<-cbind(cell=hcells,cluster=cluster_high_res_cells)
cluster.numbers<-sort(unique(cluster_high_res_cells))
# Detection of clusters that are completely disaggregated
fully_disaggregated<-c(rep(NA,length(cluster.numbers)))
for(i in 1:length(cluster.numbers)){
fully_disaggregated[i]<-identical(length(which(spam[cluster.numbers[i],]>0)),length(which(cluster_high_res_cells==cluster.numbers[i])))
}
cluster_info<-cbind(number=cluster.numbers,fully_disaggregated=fully_disaggregated)
# order from first to last affected cluster
cluster_info[order(cluster_info[,1]),]
############################## Create new spam-file #####################################
# Cells of the considered region and the remaining cells
# from disaggregated clusters, located also in the particular region,
# are listed.
# Get the cluster cell relationship
tmp<-new.magpie(paste("GLO",1:dim(spam)[1],sep="."),"y1995","cluster")
tmp[,1,1]<-1:dim(spam)[1]
cluster_cells<-as.vector(speed_aggregate(tmp, t(spam))[,1,1])
suppressWarnings(try(rm(final.spam)))
# Fill all cluster up to the first affected
if(cluster_info[1,1]>1){
final.spam<-spam[1:(cluster_info[1,1]-1),]
delete_first<-FALSE
} else{
final.spam<-spam[1,]
delete_first=TRUE
}
for(cluster in 1:dim(cluster_info)[1]){
# Fill the remainder of the cluster if it exists
if(cluster_info[cluster,2]==FALSE){
tmp<-spam[cluster_info[cluster,1],]
tmp[cell_info[which(cell_info[,2]==cluster_info[cluster,1]),1]]<-0
final.spam<-rbind(final.spam,tmp)
}
# Now add the cells that are extracted from cluster
cells<-cell_info[which(cell_info[,2]==cluster_info[cluster,1]),1]
tmp<-array(0,dim=c(length(cells),dim(spam)[2]))
for(i in 1:length(cells)){
tmp[i,cells[i]]<-1
}
final.spam<-rbind(final.spam,tmp)
# Now add the clusters, that lie between cluster and the next cluster
if(cluster==dim(cluster_info)[1]){
next_cluster<-dim(spam)[1]+1
} else {
next_cluster<-cluster_info[cluster+1,1]
}
if(next_cluster>cluster_info[cluster,1]+1){
tmp<-spam[(cluster_info[cluster,1]+1):(next_cluster-1),]
final.spam<-rbind(final.spam,tmp)
}
}
if(delete_first){
final.spam<-final.spam[2:dim(final.spam)[1],]
}
final.spam<-as.spam(final.spam)
write.spam(final.spam,file_name=path(folder,output_spam_file))
#########################################################################################
\ No newline at end of file
# (C) 2008-2016 Potsdam Institute for Climate Impact Research (PIK),
# authors, and contributors see AUTHORS file
# This file is part of MAgPIE and licensed under GNU AGPL Version 3
# or later. See LICENSE file or go to http://www.gnu.org/licenses/
# Contact: magpie@pik-potsdam.de
# *********************************************************************
# *** This script calculates low resolution input data under use of ***
# *** a high resolution base data set ***
# *********************************************************************
aggregation <- function(input_file = "path/input.tgz", # path to the data that should be used for aggregation
regionmapping = "../config/regionmappingH12.csv", # path to regionmapping
output_file = "path/output.tgz", # file path to which the data should be written
rev = 23, # MAgPIE revision
res_high = 0.5, # Input high resolution
res_low = "n100", # input low resolution
hcells = NULL,
weight = NULL,
nrepeat = 5,
nredistribute = 0,
sum_spam_file = NULL,
debug = FALSE,
seed = NULL) {
require(magclass)
require(luscale)
require(lucode)
require(luplot)
require(digest)
require(tools)
require(moinput)
cat("Aggregate data:\n")
umask <- Sys.umask("2")
on.exit(Sys.umask(umask))
# read spatial header information (spatial_header and corresponding regionscode)
spatial_header <- spatial_header(regionmapping)
regionscode <- regionscode(regionmapping)
wkey <- ifelse(is.null(weight), "", gsub(".","",paste0("_",names(weight),weight,collapse=""),fixed=TRUE))
res_out <- paste0(res_low,ifelse(is.null(hcells),"",digest(hcells)),wkey)
#########################################################################################
if(rev<24) stop('Converter is not written for revisions < 24. Please use an older revision of the script instead!')
if(is.numeric(res_out)) stop('Resolution "',res_out, '" is currently not supported as only cluster based aggregation is available (e.g. "n100").')
res_high <- format(as.numeric(res_high), nsmall = 1)
finput <- paste0("tmp/input_",format(Sys.time(),"%Y%m%d%H%M%S"))
dir.create(finput,recursive=TRUE)
cat("Extract inputs to ",finput,".\n")
file.copy(input_file,paste0(finput,"/data.tgz"))
cwd <- getwd()
setwd(finput)
trash <- system("tar -xvf data.tgz", intern = TRUE)
unlink("data.tgz")
setwd(cwd)
foutput <- paste0("tmp/aggregation_",format(Sys.time(),"%Y%m%d%H%M%S"))
dir.create(foutput,recursive=TRUE)
cat("Use temp dir ",foutput,".\n")
# create info file
writeInfo <- function(f_in, f_out, res_high, res_out, regionscode, nrepeat, nredistribute, input_file, output_file) {
functioncall <- paste(deparse(sys.call(-1)), collapse = "")
info <- readLines(f_in)
info <- c(info,'','aggregation settings:',
paste('* Input resolution:',res_high),
paste('* Output resolution:',res_out),
paste('* Input file:',input_file),
paste('* Output file:',output_file),
paste('* Regionscode:',regionscode),
paste('* (clustering) n-repeat:',nrepeat),
paste('* (clustering) n-redistribute:',nredistribute),
paste('* Call:', functioncall),
'')
cat(info,file=f_out,sep='\n')
}
writeInfo(path(finput,'info.txt'),path(foutput,'info.txt'), res_high, res_out, regionscode, nrepeat, nredistribute, input_file, output_file)
# Create new grid
if(is.null(sum_spam_file)){
spam <- clusterspam(lr=res_low,hr=res_high,ifolder=finput,ofolder=foutput,cfiles=c("lpj_yields_0.5", "lpj_airrig", rep("transport_distance",16)),spatial_header=spatial_header, weight=weight, seed=seed)
} else {
if(!file.exists(sum_spam_file)) stop("Spam file", sum_spam_file," not found")
file.copy(sum_spam_file,sub("res_in",res_high,sub("res_out",res_low,path(foutput,"res_in-to-res_out_sum.spam"))))
}
# Create new grid combining low_res and high_res if necessary
if(exists("hcells")) if(!is.null(hcells)) {
source_include <- TRUE
input_spam_file<-sub("res_in",res_high,sub("res_out",res_low,"res_in-to-res_out_sum.spam"))
output_spam_file<-sub("res_in",res_high,sub("res_out",res_out,"res_in-to-res_out_sum.spam"))
folder<-foutput
sys.source("combine.resolutions.R",envir=new.env())
write.magpie(as.magpie(hcells),path(foutput,"highres_cells.mz"))
file.remove(path(foutput,input_spam_file))
}
try(plotspam(sub("res_in",res_high,sub("res_out",res_out,path(foutput,"res_in-to-res_out_sum.spam"))),name=paste("spamplot",paste0(res_low,wkey),regionscode,sep="_"),folder=foutput,color="random"))
################################# Create SPAM files #####################################
f <- list()
#file name templates (just replace fname and fspam to get the right file name)
tinput <- sub("res_in",res_high,path(finput,"fname_res_in.mz"))
tspam <- sub("res_in",res_high,sub("res_out",res_out,path(foutput,"res_in-to-res_out_fspam.spam")))
### area_weighted_mean ###
if (rev >= 31){
x <- rowSums(read.magpie(sub("fname","avl_land_t",tinput))[,1,])
} else {
x <- rowSums(read.magpie(sub("fname","avl_land",tinput)))
}
rel <- sub("fspam","sum",tspam)
fname <- sub("fspam","area_weighted_mean",tspam)
create_spam(x,rel,fname=fname)
cat("SPAM area_weighted_mean created!\n")
### crop_area_weighted_mean ###
if(rev >= 31){
x <- rowSums(read.magpie(sub("fname","avl_land_t",tinput))[,1,"crop"])
} else {
x <- rowSums(read.magpie(sub("fname","avl_land",tinput))[,1,"crop"])
}
rel <- sub("fspam","sum",tspam)
fname <- sub("fspam","crop_area_weighted_mean",tspam)
create_spam(x,rel,fname=fname)
cat("SPAM crop_area_weighted_mean created!\n")
### soilc_weighted_mean ###
x <- read.magpie(sub("fname","lpj_carbon_stocks",tinput))[,1,"other.soilc"]
rel <- sub("fspam","sum",tspam)
fname <- sub("fspam","soilc_weighted_mean",tspam)
create_spam(x,rel,fname=fname)
cat("SPAM soilc_weighted_mean created!\n")
### bph_mask ###
if (rev >= 44){
x <- read.magpie(sub("fname","f32_bph_mask",tinput))
rel <- sub("fspam","sum",tspam)
fname <- sub("fspam","bph_mask",tspam)
create_spam(x,rel,fname=fname)
cat("SPAM bph_mask created!\n")
}
################################### Aggregate data ######################################
f <- NULL
# USAGE:
# f[<name of data set>] <- <aggregation to use>
f["aff_unrestricted"] <- "area_weighted_mean"
f["aff_noboreal"] <- "area_weighted_mean"
f["aff_onlytropical"] <- "area_weighted_mean"
f["koeppen_geiger"] <- "area_weighted_mean"
f["avl_land_si"] <- "sum"
f["avl_irrig"] <- "sum"
f["protect_area"] <- "sum"
f["lpj_airrig"] <- "crop_area_weighted_mean"
f["lpj_yields"] <- "crop_area_weighted_mean"
f["transport_distance"] <- "area_weighted_mean"
f["lpj_carbon_stocks"] <- "area_weighted_mean"
f["lpj_watavail_grper"] <- "sum"
f["lpj_envflow_grper"] <- "sum"
f["watdem_nonagr_grper"] <- "sum"
f["f59_som_initialisation_pools"] <- "sum"
if (rev >= 25 & rev < 34 | rev >= 38) {
f["rr_layer"] <- "area_weighted_mean"
f["luh2_side_layers"] <- "area_weighted_mean"
}
if (rev >= 26 & rev < 34) {
f["f38_croparea_initialisation"] <- "sum"
} else if (rev >= 34){
f["f30_croparea_initialisation"] <- "sum"
}
if (rev >= 29 & rev < 34) {
f["forestageclasses"] <- "sum"
}
if (rev >= 31) {
f["avl_land_t"] <- "sum"
} else {
f["avl_land"] <- "sum"
}
if (rev >= 32){
f["npi_ndc_ad_aolc_pol"] <- "sum"
f["npi_ndc_aff_pol"] <- "sum"
} else if (rev >= 28) {
f["npi_ndc_ad_pol"] <- "sum"
f["npi_ndc_aff_pol"] <- "sum"
f["npi_ndc_emis_pol"] <- "sum"
} else {
f["indc_ad_pol"] <- "sum"
f["indc_aff_pol"] <- "sum"
f["indc_emis_pol"] <- "sum"
}
if (rev >= 33) {
f["lpj_carbon_topsoil"] <- "area_weighted_mean"
}
if (rev >= 36) {
f["cshare_released"] <- "area_weighted_mean"
}
if (rev >= 37) {
f["runoff"] <- "sum"
}
if (rev >= 39) {
f["f30_croparea_w_initialisation"] <- "sum"
}
if (rev >= 40) {
f["avl_irrig_luh_t"] <- "sum"
}
if (rev >= 42) {
f["f50_NitrogenFixationRateNatural"] <- "area_weighted_mean"
f["f50_AtmosphericDepositionRates"] <- "area_weighted_mean"
}
if (rev >= 44) {
f["f58_peatland_degrad"] <- "sum"
f["f58_peatland_intact"] <- "sum"
f["f32_bph_effect"] <- "bph_mask"
f["f32_bph_mask"] <- "area_weighted_mean"
}
for(n in names(f)) {
input_file <- path(finput,paste(n,"_",res_high,".mz",sep=""))
spam_file <- path(foutput,paste(res_high,
"-to-", res_out,
"_",f[n],".spam", sep = ""))
outfile <- path(foutput,paste(n,"_",res_out,".mz",sep=""))
# replace spatial names with spatial header (to make sure that the right region names are used)
input <- read.magpie(input_file)
if(!is.null(spatial_header)) getCells(input) <- spatial_header
if (n == "aff_unrestricted" | n == "aff_noboreal" | n == "aff_onlytropical") {
write.magpie(round(speed_aggregate(input,spam_file,fname=NULL)),file_name=outfile)
} else speed_aggregate(input,spam_file,fname=outfile)
}
cat(paste("aggregated", res_low, "data\n"))
################################### Copy data ######################################
f<-NULL
if (rev >= 31) {
f <- c(f,"avl_land_t_0.5.mz")
} else {
f <- c(f,"avl_land_0.5.mz")
}
if (rev >= 37) {
f <- c(f,"runoff_0.5.mz")
}
if (rev >= 42) {
f <- c(f,"f50_NitrogenFixationRateNatural_0.5.mz")
f <- c(f,"f50_AtmosphericDepositionRates_0.5.mz")
}
if (rev >= 44) {
f <- c(f,"f32_bph_mask_0.5.mz")
f <- c(f,"f32_bph_effect_0.5.mz")
}
file.copy(paste(finput,f,sep="/"),paste(foutput,f,sep="/"))
reduce_time <- function(f,finput,foutput) {
x <- read.magpie(paste(finput,f,sep="/"))
x <- setYears(x[,1,],NULL)
getCells(x) <- paste("GLO",1:59199,sep=".")
write.magpie(x,paste(foutput,f,sep="/"))
}
reduce_time("lpj_yields_0.5.mz",finput,foutput)
reduce_time("lpj_carbon_stocks_0.5.mz",finput,foutput)
if (rev >= 38) reduce_time("lpj_carbon_topsoil_0.5.mz",finput,foutput)
cat(paste("copied required high res data\n"))
cwd <- getwd()
setwd(foutput)
trash <- system("tar -czf data.tgz *", intern = TRUE)
setwd(cwd)
# create output folder
if(!dir.exists(dirname(output_file))) dir.create(dirname(output_file),recursive=TRUE)
file.copy(paste0(foutput,"/data.tgz"),output_file)
unlink(paste0(foutput,"/data.tgz"))
if(!debug) {
unlink(finput, recursive = TRUE)
unlink(foutput, recursive = TRUE)
}
}
# (C) 2008-2016 Potsdam Institute for Climate Impact Research (PIK),
# authors, and contributors see AUTHORS file
# This file is part of MAgPIE and licensed under GNU AGPL Version 3
# or later. See LICENSE file or go to http://www.gnu.org/licenses/
# Contact: magpie@pik-potsdam.de
calibrate_lai <- function(input_file = "/p/projects/landuse/data/input/lpj_input/GLUES2/sresa2/constant_co2/miub_echo_g/pft_harvest.pft",
area_file = "/p/projects/landuse/data/input/raw_data/GLUES2-sresa2-constant_co2-miub_echo_g_rev21.1001/0.5_set/calibrated_area_0.5.mz",
cntr_LAI_file = "cntr.LAI.RData",
start_year = 1901, nbands = 32, avg_range = 8, lai_range = 1:7) {
#########################################################
#### Determines for each country and each crop the best fitting LAI max value ####
#########################################################
####################################################################################
#load LPJ inputs for all seven LAI, convert tem to MAgPIe objects, convert to MAgPIE structure (crops and cellbelongings), convert to tdm/ha and store them in a list
require(ludata)
require(lpjclass)
require(magclass)
require(lucode)
keepObjects<-c(ls(),"keepObjects")
LPJ.yields<-list()
Warnings<-list()
names.lai<-vector()
for(lai in lai_range){
names.lai<-c(names.lai,paste("LAI_",lai,sep=""))
laiyield<-readLPJ(paste(input_file,"_lai",lai,".bin",sep=""),wyears=1995,syear=start_year,averaging_range=avg_range,bands=nbands,soilcells=TRUE, ncells = 67420)
#transform it to a MAgPIe object and change cellordering
#get correct MAGPie crops including surrogates
#gC/m? -> tDM/ha
laiyield<-cft.transform(as.magpie(laiyield),rev=rev,na.value=1)*0.01/0.45
# Check for yields below 0 and set them to 0
if(any(laiyield<0)){
Warnings[[lai]]<-which(laiyield<0,arr.ind=T)
laiyield[laiyield<0]<-0
warning(paste("Some negative yields set to 0."))
}
LPJ.yields[[lai]]<-list(rainfed=laiyield[,,"rainfed"],irrigated=laiyield[,,"irrigated"])
getNames(LPJ.yields[[lai]]$rainfed) <- sub("\\.rainfed$","",getNames(LPJ.yields[[lai]]$rainfed))
getNames(LPJ.yields[[lai]]$irrigated) <- sub("\\.irrigated$","",getNames(LPJ.yields[[lai]]$irrigated))
}
names(LPJ.yields)<-names.lai
rmAllbut(LPJ.yields,names.lai,input_file,list=keepObjects)
####################################################################################
####################################################################################
####################################################################################
# load the area dataset for aggregation
cell.area <-read.magpie(area_file)
cell.area_rf<-as.array(cell.area[,,"rainfed"])[,1,]
cell.area_ir<-as.array(cell.area[,,"irrigated"])[,1,]
dimnames(cell.area_rf)[[2]] <- sub(".rainfed","",dimnames(cell.area_rf)[[2]],fixed=TRUE)
dimnames(cell.area_ir)[[2]] <- sub(".irrigated","",dimnames(cell.area_ir)[[2]],fixed=TRUE)
rmAllbut("LPJ.yields","cell.area_rf","cell.area_ir","cell.area",input_file,names.lai,list=keepObjects)
####################################################################################
####################################################################################
####################################################################################
#calculate the country mean of the yield for all crops using the cell.area as a weight.
#get all MAgPIE crops except for bioenergy crops
require(ludata)
magpie.crops<-as.vector(ludata$cftrel$magpie)
magpie.crops<-magpie.crops[-which(magpie.crops=="pasture")]
magpie.crops<-magpie.crops[-which(magpie.crops=="begr")]
magpie.crops<-magpie.crops[-which(magpie.crops=="betr")]
#rainfed, irrigated and total harvested area in each country for each magpie cft according to MIRCA harmonized with FAO
cntr.area_rf<-array(NA,dim=c(length(getCountries()),length(magpie.crops)),dimnames=list(getCountries(),magpie.crops))
cntr.area_ir<-array(NA,dim=c(length(getCountries()),length(magpie.crops)),dimnames=list(getCountries(),magpie.crops))
cntr.area<-array(NA,dim=c(length(getCountries()),length(magpie.crops)),dimnames=list(getCountries(),magpie.crops))
#mean LPJ yield for each country for all LAI's
cntr.LPJyield<-array(NA,dim=c(length(getCountries()),length(magpie.crops),length(names.lai)),dimnames=list(getCountries(),magpie.crops,names.lai))
for(cntr in getCountries()){
#get the code of that country
cntr.code<-ludata$countryregions$country.code[which(ludata$countryregions$country.name==cntr)]
#get all cells belonging to that country
cntr.cells<-which(ludata$cellbelongings$country.code==cntr.code)
if(length(cntr.cells)==0){ # no area in this country
cntr.area_rf[cntr,crop]<-NA
cntr.area_ir[cntr,crop]<-NA
cntr.area[cntr,]<-NA #no harvested area possible
cntr.LPJyield[cntr,,]<-NA #no yield possible
}
else{
for(crop in magpie.crops){
#calculate total harvested area of all mirca crops belonging to the cft.
cntr.area_rf[cntr,crop]<-sum(cell.area_rf[cntr.cells,crop])
cntr.area_ir[cntr,crop]<-sum(cell.area_ir[cntr.cells,crop])
cntr.area[cntr,crop]<-(cntr.area_rf+cntr.area_ir)[cntr,crop]
#calculate the yield on country level
if(cntr.area[cntr,crop]>0){
for(lai in names.lai){
cntr.LPJyield[cntr,crop,lai]<-sum(cell.area_rf[cntr.cells,crop]*LPJ.yields[[lai]][["rainfed"]][cntr.cells,"y1995",crop]+cell.area_ir[cntr.cells,crop]*LPJ.yields[[lai]][["irrigated"]][cntr.cells,"y1995",crop])/cntr.area[cntr,crop]
}
}
else{
cntr.LPJyield[cntr,crop,]<--1
}
}
}
}
rmAllbut("cntr.LPJyield","cntr.area",input_file,magpie.crops,list=keepObjects)
##########################################################################################
##########################################################################################
###########################################################################################
# Determine the best LAI max on country level .
# If there is no data available, take the (rounded) regional average calculated with cell.area area as weight.
# If no data for the region is available, take the world average.
cntr.LAI<-array(NA,dim=c(length(getCountries()),length(magpie.crops)),dimnames=list(getCountries(),magpie.crops))
#load fao yields on country level
cntr.FAOyield<-as.matrix(getYields(level="country",physical_area=TRUE))
for(cntr in getCountries()){
#get the code of that country
cntr.code<-ludata$countryregions$country.code[which(ludata$countryregions$country.name==cntr)]
#get all cells belonging to that country
cntr.cells<-which(ludata$cellbelongings$country.code==cntr.code)
if(length(cntr.cells)==0){ # no area in this country
cntr.LAI[cntr,]=NA #no LAImax to be determined
}
else{
for(crop in magpie.crops){
if(!(is.na(cntr.FAOyield[cntr,crop]))&&!(-1 %in% cntr.LPJyield[cntr,crop,])){ #fao as well as MIRCA data are present
#determine the lai from the difference between fao and lpj
lai<-as.numeric(which.min(abs(cntr.LPJyield[cntr,crop,]-cntr.FAOyield[cntr,crop])))
cntr.LAI[cntr,crop]<-lai
}
else{
cntr.LAI[cntr,crop]<-0
}
}
}
}
rmAllbut("cntr.area","cntr.LAI",input_file,magpie.crops,list=keepObjects)
################################################################################################
# determine the region averages and world average of LAImax
magpie.regions<-as.vector(sort(unique(ludata$countryregions$MAgPIE.Region)))
reg.LAI<-array(NA,dim=c(length(magpie.regions)+1,length(magpie.crops)),dimnames=list(c(magpie.regions,"WORLD"), magpie.crops))
for(crop in magpie.crops){
for(region in magpie.regions){
#determine all countries in that region which have a calibrated LAImax
countries<-as.vector(ludata$countryregions$country.name[which(ludata$countryregions$MAgPIE.Region==region)])
countries<-countries[which(cntr.LAI[countries,crop] %in% 1:7)]
if(length(countries)>0){
reg.LAI[region,crop]<-round(sum(cntr.LAI[countries,crop]*cntr.area[countries,crop])/sum(cntr.area[countries,crop]),digits=0)
}
else{
# where no data is available in the region, fill with NA's, -1's respectively
reg.LAI[region,crop]<-0
}
}
# now determine the world average
#determine all countries in the world which have a calibrated LAImax
countries<-getCountries()[which(cntr.LAI[getCountries(),crop] %in% 1:7)]
reg.LAI["WORLD",crop]<-round(sum(cntr.LAI[countries,crop]*cntr.area[countries,crop])/sum(cntr.area[countries,crop]),digits=0)
}
rmAllbut(cntr.LAI,reg.LAI,magpie.crops,input_file,list=keepObjects)
################################################################################################
# replace the missing values in cntr.LAI by region (if available) or world average.
for(cntr in getCountries()){
region<-as.vector(ludata$countryregions$MAgPIE.Region[which(ludata$countryregions$country.name==cntr)])