Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

F npi ndc bugfix #164

Merged
merged 14 commits into from
Apr 28, 2020
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
195 changes: 111 additions & 84 deletions scripts/npi_ndc/start_npi_ndc.R
Original file line number Diff line number Diff line change
Expand Up @@ -6,36 +6,34 @@
# | Contact: magpie@pik-potsdam.de
## Functions for the policy targets calculations

get_info <- function(file, grep_expression, sep, pattern="", replacement=""){
if(!file.exists(file)) return("#MISSING#")
file <- readLines(file, warn=FALSE)
tmp <- grep(grep_expression, file, value=TRUE)
tmp <- strsplit(tmp, sep)
tmp <- sapply(tmp, "[[", 2)
tmp <- gsub(pattern, replacement ,tmp)
if(all(!is.na(as.logical(tmp)))) return(as.vector(sapply(tmp, as.logical)))
if (all(!(regexpr("[a-zA-Z]",tmp) > 0))){
tmp <- as.numeric(tmp)
}
return(tmp)
}

## calculates policy for protecting different land pools from land use change
calc_NPI_NDC <- function(policyregions="iso"){
calc_NPI_NDC <- function(policyregions = "iso",
pol_def_file = "policies/policy_definitions.csv",
spatialheader_file = "../../input/spatial_header.rda",
cellmapping = "policies/country2cell.rds",
land_stock_file = "../../modules/10_land/input/avl_land_t_0.5.mz",
spam_file = Sys.glob("../../input/*-to-*_sum.spam"),
outfolder_ad_aolc = c("policies/","../../modules/35_natveg/input/"),
outfolder_aff = c("policies/","../../modules/32_forestry/input/"),
out_ad_file = "npi_ndc_ad_aolc_pol.cs3",
out_aff_file = "npi_ndc_aff_pol.cs3"){

require(magclass)
require(luscale)
require(lucode)
require(madrat)

# load the cell mapping policy
pol_mapping <- readRDS("policies/country2cell.rds")[,policyregions]
pol_mapping <- readRDS(cellmapping)
if(!(policyregions %in% names(pol_mapping))) stop("No cell mapping available for policyregions = ",policyregions)
pol_mapping <- pol_mapping[,policyregions]

##############################################################################
########## Information from the reference observed data ###########
##############################################################################

#read in cellular land cover (stock) from landuse initialization
land_stock <- read.magpie("../../modules/10_land/input/avl_land_t_0.5.mz")
land_stock <- read.magpie(land_stock_file)
getRegions(land_stock) <- paste(pol_mapping,1:59199,sep=".")

forest_stock <- dimSums(land_stock[,,c("primforest","secdforest","forestry")], dim=3)
Expand Down Expand Up @@ -65,7 +63,7 @@ calc_NPI_NDC <- function(policyregions="iso"){
capture.output(print(x, print.gap=3, row.names=FALSE), file=fname, append=TRUE)
}

pol_def <- read.csv("policies/policy_definitions.csv")
pol_def <- read.csv(pol_def_file)

addline("NPI/NDC policies")
addline("MAgPIE")
Expand All @@ -77,28 +75,33 @@ calc_NPI_NDC <- function(policyregions="iso"){
addline("## percentage protection: 0 = no protection, 1 = full protection")
addline("## Ref: BaseYear (1), Baseline (2)")

cat("NPI AD policy\n")
ptm <- proc.time()["elapsed"]
cat("Compute NPI AD policy")
addline("")
addline("##############")
addline("### NPI AD ###")
addline("##############")
npi_ad <- droplevels(subset(pol_def, policy=="npi" & landpool=="forest"))
addtable(npi_ad[,c(-2,-3)])
npi_ad <- calc_policy(npi_ad,forest_stock,pol_type="ad",pol_mapping=pol_mapping)
npi_ad <- calc_policy(npi_ad, forest_stock, pol_type="ad", pol_mapping=pol_mapping,
spatialheader_file=spatialheader_file, spam_file=spam_file)
getNames(npi_ad) <- "npi.forest"
cat(paste0(" (time elapsed: ",format(proc.time()["elapsed"]-ptm,width=6,nsmall=2,digits=2),"s)\n"))

cat("NDC AD policy\n")
cat("Compute NDC AD policy")
addline("")
addline("##############")
addline("### NDC AD ###")
addline("##############")
ndc_ad <- droplevels(subset(pol_def, policy=="ndc" & landpool=="forest"))
addtable(ndc_ad[,c(-2,-3)])
ndc_ad <- calc_policy(ndc_ad,forest_stock,pol_type="ad",pol_mapping=pol_mapping)
ndc_ad <- calc_policy(ndc_ad, forest_stock, pol_type="ad", pol_mapping=pol_mapping,
spatialheader_file=spatialheader_file, spam_file=spam_file)
getNames(ndc_ad) <- "ndc.forest"
#Set all values before 2015 to NPI values; copy the values til 2010 from the NPI data
ndc_ad[,which(getYears(ndc_ad,as.integer=TRUE)<=2020),] <-
npi_ad[,which(getYears(npi_ad,as.integer=TRUE)<=2020),]
cat(paste0(" (time elapsed: ",format(proc.time()["elapsed"]-ptm,width=6,nsmall=2,digits=2),"s)\n"))


addline("")
Expand All @@ -107,24 +110,27 @@ calc_NPI_NDC <- function(policyregions="iso"){
addline("## percentage protection: 0 = no protection, 1 = full protection")
addline("## Ref: BaseYear (1), Baseline (2)")

cat("NPI AOLC policy\n")
cat("Compute NPI AOLC policy")
addline("")
addline("################")
addline("### NPI AOLC ###")
addline("################")
npi_aolc <- droplevels(subset(pol_def, policy=="npi" & landpool=="other"))
addtable(npi_aolc[,c(-2,-3)])
npi_aolc <- calc_policy(npi_aolc,land_stock[,,"other"],pol_type="ad",pol_mapping=pol_mapping)
npi_aolc <- calc_policy(npi_aolc, land_stock[,,"other"], pol_type="ad", pol_mapping=pol_mapping,
spatialheader_file=spatialheader_file, spam_file=spam_file)
getNames(npi_aolc) <- "npi.other"
cat(paste0(" (time elapsed: ",format(proc.time()["elapsed"]-ptm,width=6,nsmall=2,digits=2),"s)\n"))

cat("NDC AOLC policy\n")
cat("Compute NDC AOLC policy")
addline("")
addline("################")
addline("### NDC AOLC ###")
addline("################")
ndc_aolc <- droplevels(subset(pol_def, policy=="ndc" & landpool=="other"))
addtable(ndc_aolc[,c(-2,-3)])
ndc_aolc <- calc_policy(ndc_aolc,land_stock[,,"other"],pol_type="ad",pol_mapping=pol_mapping)
ndc_aolc <- calc_policy(ndc_aolc,land_stock[,,"other"],pol_type="ad",pol_mapping=pol_mapping,
spatialheader_file=spatialheader_file, spam_file=spam_file)
getNames(ndc_aolc) <- "ndc.other"
#Set all values before 2015 to NPI values; copy the values til 2010 from the NPI data
ndc_aolc[,which(getYears(ndc_aolc,as.integer=TRUE)<=2020),] <-
Expand All @@ -135,33 +141,40 @@ calc_NPI_NDC <- function(policyregions="iso"){
none_ad_aolc_pol[] <- 0
getNames(none_ad_aolc_pol) <- c("none.forest","none.other")
ad_aolc_pol <- mbind(none_ad_aolc_pol,npi_ad,npi_aolc,ndc_ad,ndc_aolc)
write.magpie(round(ad_aolc_pol,3), "policies/npi_ndc_ad_aolc_pol.cs3")

adfiles <- paste0(outfolder_ad_aolc, out_ad_file)
write.magpie(ad_aolc_pol, adfiles[1])
if(length(adfiles >1)) for(i in 2:length(adfiles)) file.copy(adfiles[1],adfiles[i], overwrite=TRUE)
cat(paste0(" (time elapsed: ",format(proc.time()["elapsed"]-ptm,width=6,nsmall=2,digits=2),"s)\n"))

addline("")
addline("##----------------------------------------------------------------------------")
addline("## Afforestation - AFF (Mha)")
addline("## Ref: BaseYear (1), Baseline (2)")

cat("NPI AFF policy\n")
cat("Compute NPI AFF policy")
addline("")
addline("###############")
addline("### NPI AFF ###")
addline("###############")
npi_aff <- droplevels(subset(pol_def, policy=="npi" & landpool=="affore"))
addtable(npi_aff[,c(-2,-3)])
npi_aff <- calc_policy(npi_aff,land_stock,pol_type="aff",pol_mapping=pol_mapping,
weight=dimSums(land_stock[,2005,c("crop","past")]))
npi_aff <- calc_policy(npi_aff, land_stock, pol_type="aff", pol_mapping=pol_mapping,
weight=dimSums(land_stock[,2005,c("crop","past")]),
spatialheader_file=spatialheader_file, spam_file=spam_file)
getNames(npi_aff) <- "npi"
cat(paste0(" (time elapsed: ",format(proc.time()["elapsed"]-ptm,width=6,nsmall=2,digits=2),"s)\n"))

cat("NDC AFF policy\n")
cat("Compute NDC AFF policy")
addline("")
addline("###############")
addline("### NDC AFF ###")
addline("###############")
ndc_aff <- droplevels(subset(pol_def, policy=="ndc" & landpool=="affore"))
addtable(ndc_aff[,c(-2,-3)])
ndc_aff <- calc_policy(ndc_aff,land_stock,pol_type="aff",pol_mapping=pol_mapping,
weight=dimSums(land_stock[,2005,c("crop","past")]))
ndc_aff <- calc_policy(ndc_aff, land_stock, pol_type="aff", pol_mapping=pol_mapping,
weight=dimSums(land_stock[,2005,c("crop","past")]),
spatialheader_file=spatialheader_file, spam_file=spam_file)
getNames(ndc_aff) <- "ndc"
#set all values before 2015 to NPI values; copy the values til 2010 from the NPI data
ndc_aff[,which(getYears(ndc_aff,as.integer=TRUE)<=2020),] <-
Expand All @@ -172,105 +185,119 @@ calc_NPI_NDC <- function(policyregions="iso"){
none_aff_pol[] <- 0
getNames(none_aff_pol) <- "none"
aff_pol <- mbind(none_aff_pol,npi_aff,ndc_aff)
write.magpie(round(aff_pol,6), "policies/npi_ndc_aff_pol.cs3")
afffiles <- paste0(outfolder_aff, out_aff_file)
write.magpie(aff_pol, afffiles[1])
if(length(afffiles >1)) for(i in 2:length(afffiles)) file.copy(afffiles[1],afffiles[i], overwrite=TRUE)

#copy files
file.copy("policies/npi_ndc_ad_aolc_pol.cs3",
"../../modules/35_natveg/input/npi_ndc_ad_aolc_pol.cs3",overwrite = TRUE)
file.copy("policies/npi_ndc_aff_pol.cs3",
"../../modules/32_forestry/input/npi_ndc_aff_pol.cs3",overwrite = TRUE)
cat(paste0(" (time elapsed: ",format(proc.time()["elapsed"]-ptm,width=6,nsmall=2,digits=2),"s)\n"))
}

### calc flow function
calc_flows <- function(stock,t_periods) {
flow <- stock
flow[,,] <- 0
for (y in 2:nyears(stock)) {
flow[,y,] <- (setYears(stock[,y-1,],NULL) - stock[,y,])/t_periods[y]
}
return(flow)
calc_tperiods <- function(y) {
t_periods <- c(1,y[3:length(y)]-y[2:(length(y)-1)])
names(t_periods) <- paste0("y",y[2:length(y)])
return(as.magpie(t_periods))
}

### calc annual flow function
calc_flows <- function(stock) {
stock <- stock[,c(1,1:nyears(stock)),]
y <- getYears(stock,as.integer = TRUE)
t_periods <- calc_tperiods(y)
out <- (setYears(stock[,1:(nyears(stock)-1),],getYears(stock)[2:nyears(stock)])
- stock[,2:nyears(stock),])/t_periods
return(out)
}

### calc npi & ndc policy
calc_policy <- function(policy,stock,pol_type="aff",pol_mapping, weight=NULL) {
calc_policy <- function(policy, stock, pol_type="aff", pol_mapping=pol_mapping,
weight=NULL, spatialheader_file = "../../input/spatial_header.rda",
spam_file = Sys.glob("../../input/*-to-*_sum.spam")) {
## pol_type = {"aff","ad"}

require(luscale)
require(lucode)
require(madrat)

#extent stock beyond last observed value with constant values from the last year
ly <- tail(getYears(stock,as.integer=TRUE),1)
ly <- ly - ly%%5
stock <- stock[,seq(1995,ly,5),]
stock_extent <- new.magpie(getCells(stock), seq(ly+5,2150,5),
getNames(stock), stock[,ly,])
stock <- mbind(stock, stock_extent)
rm(stock_extent)
year_extension <- seq(ly+5,2150,5)
stock <- stock[,c(seq(1995,ly,5),rep(ly,length(year_extension))),]
getYears(stock) <- c(seq(1995,ly,5), year_extension)

#the the years and periods
#full years
tp <- getYears(stock, as.integer=TRUE)
t_periods <- c(1, sapply(seq_along(tp[-1]), function(i) tp[i+1] - tp[i]))

#select and filter countries that exist in the chosen policy mapping
policy_countries <- intersect(policy$dummy,unique(pol_mapping))
policy <- policy[policy$dummy %in% policy_countries,]

#create key to distinguish different cases of baseyear, targetyear combinations
policy$key <- paste(policy$baseyear,policy$targetyear)


#set stock to zero or Inf for countries without policies
# (representing no constraint for min and max constraints)
if(pol_type=="ad"){
stock[setdiff(getRegions(stock),policy_countries),,] <- 0
stock[!(sub("\\..*$","",getCells(stock)) %in% policy_countries),,] <- 0
#calculate flows
flow <- calc_flows(stock,t_periods)
#account only for positive flows
flow <- calc_flows(stock)
#account only for positive flows, i.e. deforestation
flow[flow < 0] <- 0
}

#Initialize magpie_policy with 0 (country level)
#This is a return object of this function and contains policy targets at
#cluster level
magpie_policy <- new.magpie(unique(pol_mapping),tp,NULL,0)
for (i in policy_countries) {
cat(i,round(which(policy_countries==i)/length(policy_countries)*100),"%\n")
keys <- unique(policy$key)
for (i in keys) {
#get baseyear and targetyear
baseyear <- policy[which(policy$dummy==i),"baseyear"]
targetyear <- policy[which(policy$dummy==i),"targetyear"]
y_pol <- tp[tp>= baseyear & tp<=targetyear]
y_pol_forever <- tp[tp>= targetyear]
tmp <- as.integer(strsplit(i," ")[[1]])
baseyear <- tmp[1]
targetyear <- tmp[2]
countries <- policy$dummy[policy$key==i]
y_full <- tp[tp >= baseyear]

#set target in targetyear
#percentage: 0 = no reduction, 1 = full reduction of deforestation
#Mha if pol_type=="aff"
magpie_policy[i,targetyear,] <- policy[which(policy$dummy==i),"target"]
magpie_policy[countries,targetyear,] <- policy[policy$key==i,"target"]
#interpolate between baseyear and targetyear
magpie_policy[i,y_pol,] <- time_interpolate(magpie_policy[i,c(baseyear,targetyear),],y_pol)
#set same target for all years after targetyear
magpie_policy[i,y_pol_forever,] <- setYears(magpie_policy[i,targetyear,],NULL)
#get target type
targettype <- policy[which(policy$dummy==i),"targettype"] #1 baseyear target #2 baseline target

if(targetyear==baseyear) {
magpie_policy[countries,y_full,] <- setYears(magpie_policy[countries,baseyear,],NULL)
} else {
magpie_policy[countries,y_full,] <- time_interpolate(magpie_policy[countries,c(baseyear,targetyear),],y_full,
extrapolation_type = "constant")
}
#set reference flow based on target type
if(pol_type=="ad") {
stock[i,tp>targetyear,] <- setYears(stock[i,targetyear,],NULL)
ref_flow <- new.magpie(getCells(stock),getYears(stock),NULL,0)
if (targettype == 1) {
ref_flow[i,,] <- setYears(flow[i,baseyear,],NULL)
} else if (targettype == 2) {
ref_flow[i,,] <- flow[i,,]
} else stop("unknow targettype; needs to be 1 or 2")
c_index <- (sub("\\..*$","",getCells(stock)) %in% countries)
stock[c_index,tp>targetyear,] <- setYears(stock[c_index,targetyear,],NULL)
targettypes <- unique(policy[policy$key==i,"targettype"])
if(!all(targettypes %in% 1:2)) stop("unknow targettype; needs to be 1 or 2")
if (any(targettypes == 1)) {
t1countries <- policy$dummy[policy$key==i & policy$targettype==1]
t1c_index <- (sub("\\..*$","",getCells(flow)) %in% t1countries)
flow[t1c_index,,] <- setYears(flow[t1c_index,baseyear,],NULL)
}
}
}

#calculate the reduction target in absolute numbers
rel <- data.frame(from=pol_mapping,to=paste(pol_mapping,1:length(pol_mapping),sep="."))
rel <- data.frame(from=pol_mapping,to=paste(pol_mapping,1:length(pol_mapping),sep="."), stringsAsFactors = FALSE)
if(pol_type=="aff") {
magpie_policy <- speed_aggregate(x=magpie_policy, rel=rel, weight=weight)
magpie_policy <- madrat::toolAggregate(x=magpie_policy, rel=rel, weight=weight)
} else if(pol_type=="ad") {
magpie_policy <- speed_aggregate(x=magpie_policy, rel=rel)
magpie_policy <- magpie_policy * ref_flow * t_periods + stock
t_periods <- calc_tperiods(c(tp[1],tp))
magpie_policy <- magpie_policy * flow * t_periods + stock
}

load("../../input/spatial_header.rda")
load(spatialheader_file)
getCells(magpie_policy) <- spatial_header

res_out <- get_info("../../input/info.txt","^\\* Output ?resolution:",": ")
res_high <- get_info("../../input/info.txt","^\\* Input ?resolution:",": ")
spam_file <- path("../../input",paste0(res_high,"-to-",res_out,"_sum.spam"))
magpie_policy <- speed_aggregate(magpie_policy,spam_file)

return(magpie_policy)
Expand Down