In PrincetonUniversity/agroEcoTradeoff: Landuse Tradeoff analysis model

Overview

The following code runs a demonstration of how one might try to balance different land use priorities. We test different weights that reflect the level of priority we might attach to each of the following four goals:

1. Maximizing the potential yield we might get from different crops
2. Minimizing the loss of carbon resulting from land use conversion
3. Minimizing the loss of biodiversity and protected areas
4. Minimizing transport time, a proxy for cost

The model balances these priorities according to different permutations of weights as it sets out to satisfy crop production targets that we set. These targets can be thought of as scenarios for future agricultural demand.

In the following code, we are going to try evaluate which set of weights is optimal, given an across the board increase that needs to be met for 9 different crops.

In this first chunk of code, you are getting to set just two options:

1. The multiple by which existing production must be increased to meet future demand across all crops;
2. Which of the four priorities we are going to care about. You can choose just two, three, or all four. For each pairwise permutation there are 5 combinations, so the model runs a bit longer depending on how many priorities you choose.

# set the target multiplier
targ <- 4  

# choose which of the priorities you want to evaluate. "Ag" = agriculture, or 
# the desire to maximixe potential agricultural productivity, "C" = carbon, 
# or the desire to minimize carbon loss, "bd" = biodiversity, or the desire to 
# minimize protected area loss or higher biodiversity unprotected areas, 
# "cost" = the desire to minimize transport costs. 
# you can simply run all four to get all possible combinations, but could also
# just run a subset (e.g. the commented out one below)
# cnames <- c("Ag", "C")
cnames <- c("Ag", "C", "cost")
# cnames <- c("Ag", "C", "bd", "cost")

Back to top

This chunk will execute the model and plot the outputs, which takes a little while to run

library(agroEcoTradeoff)
yblist <- list(yb1 <- c(1, 1))
step <- 0.25
ot <- pareto(cnames = cnames, step = step, yblist = yblist, targ = targ)
bcode <- ot$bcode
save(bcode, file = paste0("external/output/batch/dt/bcode.rda"))
save(ot, file = paste0("external/output/batch/dt/", bcode, "/ot.rda"))

Weights evaluated

# Load the output tables
library(agroEcoTradeoff)
setwd(set_base_path())
load(paste0("external/output/batch/dt/bcode.rda"))
dnm <- full_path("external/output/batch/dt", bcode)
load(paste0(dnm, "/out_tables.rda"))
load(paste0(dnm, "/parms.rda"))
load(paste0(dnm, "/ot.rda"))

parms[, cnames]

Back to top

Tradeoff plots

Sorted by converted area and carbon lost

The first column are the results sorted by total land converted, and the second by carbon loss. Both are ordered from least to most. The weighting scenario numbers are listed at the bottom of the bars for reference.

ott <- data.frame(ot$table)
ott$colors <- rainbow(length(ott$ind))
ott$ind <- seq(1, length(ott$ind))

odfs <- lapply(c("land", "carbon", "biodiversity", "cost"), function(x) {
  if(x %in% colnames(ott)) {
    odf <- ott[order(ott[, x]), ]
  } else{ 
    odf <- NULL
  }
  odf
})

pvec <- c("land", "carbon", "biodiversity", "cost")
tvec <- c("by area", "by carbon", "by PA loss", "by cost")
yvec <- c("Area (ha X 1000)", "C (ton X 10000)", "PAs (ha X 1000)", "Time")
div <- c(1 / 1000, 1 / 10000, 1 / 1000, 60 * 24 * 7)
dfs <- which(sapply(odfs, is.data.frame))
if(length(dfs) == 2) {
 par(mfcol = c(2, 2), oma = c(2, 5, 1, 0), mar = c(1, 1, 1, 1) + 0.1, 
    xaxs = "i", mgp = c(1, 0.1, 0), tcl = -0.2)
} else if(length(dfs) == 3) {
  par(mfcol = c(3, 2), oma = c(2, 5, 1, 0), mar = c(1, 1, 1, 1) + 0.1, 
    xaxs = "i", mgp = c(1, 0.1, 0), tcl = -0.2)
} else if(length(dfs) == 4) {
  par(mfcol = c(4, 2), oma = c(2, 5, 1, 0), mar = c(1, 1, 1, 1) + 0.1, 
    xaxs = "i", mgp = c(1, 0.1, 0), tcl = -0.2)
}
for(i in 1:2) {  #i <- 1
  d <- odfs[[i]]
  if(!is.null(d)) {
   for(j in 1:length(pvec)) {  # j <- 1
     if(pvec[[j]] %in% colnames(d)) {
      barplot(round(d[, pvec[j]] * div[j]), names.arg = d$ind, col = d$colors,
              yaxt = "n", las = 2, cex.axis = 0.3)
      if(j == 1) mtext(tvec[i], side = 3)
      if(i == 1) mtext(yvec[j], side = 2, cex = 0.8, line = 3)
      if(i == 1) {
       axis(2, las = 2, cex = 0.6)
      } else {
       axis(2, labels = FALSE)
      }
     }
   } 
  }
}

Back to top

Sorted by protected area loss and total cost

The first column are the results sorted by area of protected areas converted (a good proxy for biodiversity impacts, and the second by transportation costs. The weighting scenario numbers are listed at the bottom of the bars for reference.

if(length(dfs) == 3) {
  par(mfcol = c(3, 2), oma = c(2, 5, 1, 0), mar = c(1, 1, 1, 1) + 0.1, 
    xaxs = "i", mgp = c(1, 0.1, 0), tcl = -0.2)
} else if(length(dfs) == 4) {
  par(mfcol = c(4, 2), oma = c(2, 5, 1, 0), mar = c(1, 1, 1, 1) + 0.1, 
    xaxs = "i", mgp = c(1, 0.1, 0), tcl = -0.2)
}
for(i in 3:length(pvec)) {  #i <- 1
  d <- odfs[[i]]
  if(!is.null(d)) {
   for(j in 1:length(pvec)) {  # j <- 1
     if(pvec[[j]] %in% colnames(d)) {
      barplot(round(d[, pvec[j]] * div[j]), names.arg = d$ind, col = d$colors,
              yaxt = "n", 
              las = 2, cex.axis = 0.3)
      if(j == 1) mtext(tvec[i], side = 3)
      if(i == 3) mtext(yvec[j], side = 2, cex = 0.8, line = 3)
      if(i == 3) {
       axis(2, las = 2, cex = 0.6)
      } else {
       axis(2, labels = FALSE)
      }
     }
   } 
  }
}

Back to top

PrincetonUniversity/agroEcoTradeoff documentation built on May 8, 2019, 3:12 p.m.