In bocinsky/gutaker2020: Research Compendium for Gutaker et al. 2020 thermal niche model

This RMarkdown document creates the thermal niche model reported in

Gutaker et al. In review.

# Set the knitting behavior
knitr::opts_chunk$set(echo = TRUE,
                      warning = FALSE,
                      message = FALSE,
                      collapse = TRUE,
                      cache = FALSE,
                      cache.lazy = FALSE,
                      results = 'hold',
                      out.width = "100%",
                      fig.height = 5,  
                      comment = "#>",
                      fig.path = here::here("analysis/figures/")
)

library(gutaker2020)
library(sf)
library(magrittr)
library(fields)

# Set the behavior for printing scientific notation
# options(digits = 3, scipen = -2)
options(scipen = 999,
        knitr.table.format = "html")

# Force Raster to load large rasters into memory
raster::rasterOptions(chunksize = 2e+09,
                      maxmemory = 2e+10)

# Create the plan for (possibly parallel) execution
future::plan(future::multisession, 
             workers = min(parallel::detectCores(), 
                           params$cores),
             gc = TRUE)

raw_data <- here::here("analysis/data/raw_data/")
derived_data <- here::here("analysis/data/derived_data/")
figures <- here::here("analysis/figures/")

# Create output directories
dir.create(raw_data,
           recursive = TRUE,
           showWarnings = FALSE)
dir.create(derived_data,
           recursive = TRUE,
           showWarnings = FALSE)
dir.create(figures,
           recursive = TRUE,
           showWarnings = FALSE)

Set parameters {-}

In this section we define the study area by a bounding box, the calendar years used to define climatology (here, 1961--1990), and a plotting theme for ggplot.

# Set a series of colors for the crops
crop_colors <- RColorBrewer::brewer.pal(6,"Dark2")

# A function to set an objects class
set_class <- function(x, classes){
  class(x) <- classes
  return(x)
}

#Define the study region
ASIA_poly <-
  raster::extent(55,150,-11,60) %>%
  FedData::polygon_from_extent("+proj=longlat") %>%
  sf::st_as_sf()

# Set the calibration period for paleoclimate reconstructions
calibration.years <- 1961:1990

# A ggplot2 theme for Nature Publishing Group
nature_theme <- ggplot2::theme(panel.grid.major = ggplot2::element_line(size = 0.5, 
                                                                        color = "grey"),
                               axis.line = ggplot2::element_line(size = 0.7, 
                                                                 color = "black"),
                               legend.position = c(0.85, 0.7),
                               text = ggplot2::element_text(size = 14))

Load ETOPO5 grids {-}

Indicator kriging takes place on the ETOPO5 digital elevation model. Here, we load the ETOPO5 data from the geomapdata package, then mask out the oceans adn large lakes, as well as islands outside of mainland Asia.

message("Preparing the ETOPO5 grid-aligned dataset")

time_check <-  Sys.time()

ne_dir <- paste0(raw_data,"NaturalEarth/")

if(!params$clean & file.exists(paste0(raw_data,"ASIA_rast_etopo5.tif"))){

  ASIA_rast_etopo5 <- raster::raster(paste0(raw_data,"ASIA_rast_etopo5.tif"))

}else{

  # Get the Natural Earth country lakes data
  dir.create(ne_dir, 
             showWarnings = FALSE, 
             recursive = TRUE)

  FedData::download_data(
    url = "http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_lakes.zip",
    destdir = ne_dir
  )

  unzip(
    paste0(ne_dir,"ne_10m_lakes.zip"), 
    exdir = paste0(ne_dir,"ne_10m_lakes/")
  )

  ne_10m_lakes <-  paste0(ne_dir,"ne_10m_lakes/ne_10m_lakes.shp") %>%
    sf::st_read() %>%
    sf::st_make_valid() %>%
    sf::st_transform(sf::st_crs(ASIA_poly)) %>%
    sf::st_intersection(ASIA_poly)

  # Get the Natural Earth country boundaries data
  dir.create(
    paste0(ne_dir,"ne_10m_admin_0_countries_lakes/"),
    showWarnings = FALSE, 
    recursive = TRUE
  )
  FedData::download_data(
    url = "http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_0_countries_lakes.zip",
    destdir = ne_dir
  )

  unzip(
    paste0(ne_dir,"ne_10m_admin_0_countries_lakes.zip"), 
    exdir = paste0(ne_dir,"ne_10m_admin_0_countries_lakes/")
  )

  ne_10m_admin_0_countries_lakes <- 
    paste0(ne_dir,"ne_10m_admin_0_countries_lakes/ne_10m_admin_0_countries_lakes.shp") %>%
    sf::st_read() %>%
    sf::st_transform(sf::st_crs(ASIA_poly)) %>%
    sf::st_intersection(ASIA_poly) %>%
    dplyr::filter(!(NAME %in% c("Scarborough Reef",
                                "Maldives",
                                "Spratly Is.",
                                "Oman",
                                "United Arab Emirates",
                                "Saudi Arabia")))

  library(geomapdata)
  data("ETOPO5")

  ASIA_rast_etopo5 <- ETOPO5 %>%
    t() %>%
    raster::raster(xmn = 0,
                   xmx = 360,
                   ymn = -90,
                   ymx = 90,
                   crs = sp::CRS("+proj=longlat +ellps=clrk66 +no_defs")) %>%
    raster::crop(ASIA_poly %>%
                   sf::st_transform("+proj=longlat +ellps=clrk66 +no_defs") %>%
                   as("Spatial")) %>%
    raster::mask(ne_10m_admin_0_countries_lakes %>%
                   sf::st_transform("+proj=longlat +ellps=clrk66 +no_defs") %>%
                   as("Spatial")) %>%
    raster::mask(ne_10m_lakes %>%
                   sf::st_transform("+proj=longlat +ellps=clrk66 +no_defs") %>%
                   as("Spatial"),
                 inverse = TRUE) %T>%
    raster::writeRaster(filename = paste0(raw_data,"ASIA_rast_etopo5.tif"),
                        datatype="INT2S",
                        options=c("COMPRESS=DEFLATE", "ZLEVEL=9", "INTERLEAVE=BAND", "PHOTOMETRIC=MINISWHITE"),
                        overwrite=T,
                        setStatistics=FALSE)
  rm(ETOPO5)
}

countries <- 
  paste0(ne_dir,"ne_10m_admin_0_countries_lakes/ne_10m_admin_0_countries_lakes.shp") %>%
  sf::st_read()  %>%
  sf::st_transform(sf::st_crs(ASIA_poly)) %>%
  sf::st_intersection(ASIA_poly) %>%
  dplyr::filter(!(NAME %in% c("Scarborough Reef",
                              "Maldives",
                              "Spratly Is.",
                              "Oman",
                              "United Arab Emirates",
                              "Saudi Arabia"))) %>%
  dplyr::select(NAME, geometry)

message("ETOPO5 grid-aligned dataset preparation complete: ", capture.output(Sys.time() - time_check))

Prepare the GHCN data {-}

Here, we download daily weather records from the Global Historical Climatology Network,

## Downloads and cleans daily climate records from the Global Historical Climate Database.
message("Preparing the daily climate records from the Global Historical Climate Database")
time_check <-  Sys.time()
GHCN.data.final <- prepare_ghcn(region = ASIA_poly %>%
                                  as("Spatial"),
                                label = "ASIA_poly",
                                calibration.years = calibration.years,
                                google_maps_elevation_api_key = params$google_maps_elevation_api_key,
                                raw_dir = paste0(raw_data,"ghcn/"),
                                derived_dir = paste0(derived_data,"ghcn/"),
                                force.redo = params$clean)
message("GHCN preparation complete: ", capture.output(Sys.time() - time_check))
## An example of plotting the GHCN data
# climate_plotter(data = GHCN.data.final, station = "CHM00051334", element = "TMIN")

Prepare the Marcott data {-}

# Run the script that transforms the Marcott et al. 2013 data into standard scores.
message("Preparing the Marcott et al. 2013 data")
time_check <-  Sys.time()
marcott2013 <- prepare_marcott(calibration.years = calibration.years,
                               raw_dir = paste0(raw_data,"marcott2013/"),
                               derived_dir = paste0(derived_data,"marcott2013/"))
message("Marcott et al. 2013 data preparation complete: ", capture.output(Sys.time() - time_check))

Modulating climatology by SD {-}

#### Calculating growing degree days ####
message("Modulating local climatology by Marcott et al. 2013 data")
time_check <-  Sys.time()
# How often to sample GDD, in z-space, for model tuning
# Here, we sample from -20 to 20 SD, at 1 SD interval
sample.points <- -20:20

# Read in data on different crop GDD needs
crop_GDD <- 
  readr::read_csv(here::here("inst/crops.csv"),
                  col_types = readr::cols(
                    cultivar_long = readr::col_character(),
                    cultivar = readr::col_character(),
                    crop_long = readr::col_character(),
                    crop = readr::col_character(),
                    t_base = readr::col_double(),
                    min_gdd = readr::col_integer()
                  )) %T>%
  readr::write_csv(paste0(derived_data,"/crops.csv"))

# Transform GHCN data to GDDs of each base, and modulate to Marcott
GDDs <- sort(unique(crop_GDD$t_base))

GHCN.GDD.incremented.sd <-
  furrr::future_map(GDDs,
                    # purrr::map(GDDs,
                    .options = furrr::future_options(
                      packages = c("magrittr",
                                   "gutaker2020"),
                      globals = c("sample.points",
                                  "GHCN.data.final")),
                    .f = function(base){
                      purrr::map_dfr(sample.points,
                                     function(change){

                                       GHCN.GDDs <-
                                         purrr::map_dbl(GHCN.data.final$climatology,
                                                        function(station){
                                                          return(sdModulator(data.df = station,
                                                                             temp.change.sd = change,
                                                                             t.base = base,
                                                                             t.cap = 30))
                                                        })


                                       names(GHCN.GDDs) <- names(GHCN.data.final$climatology)

                                       return(tibble::tibble(SD_change = change,
                                                             ID = names(GHCN.GDDs),
                                                             GDD = GHCN.GDDs))

                                     }) %>%
                        dplyr::left_join(GHCN.data.final$spatial %>%
                                           tibble::as_tibble() %>%
                                           magrittr::set_names(c("ID",
                                                                 "NAME",
                                                                 "elevation",
                                                                 "x",
                                                                 "y")),
                                         by = "ID")

                    })

names(GHCN.GDD.incremented.sd) <- GDDs

# stop the cluster (will free memory)
# stopCluster(cl)
message("Modulation of local climatology by Marcott et al. 2013 data complete: ", capture.output(Sys.time() - time_check))

Generate the crop niche models {-}

message("Calculating indicator Krige models")
time_check <-  Sys.time()
# Create a spatialPointsDataFrame of the etopo5 data, and convert to WGS84 ellipsoid
ASIA_rast_etopo5.sp <-
  ASIA_rast_etopo5 %>%
  magrittr::set_names("elevation") %>%
  raster::rasterToPoints(spatial = T) %>%
  sp::spTransform(sp::CRS(raster::projection(GHCN.data.final$spatial)))

# A function that generates the kriging model, then predicts
krige_and_predict <- function(dt){
  model <- fields::mKrig(x = dt[,c("x","y")],
                         y = dt$GDD_thresh,
                         Z = dt$elevation,
                         Covariance = "Exponential",
                         Distance = "rdist.earth")

  prediction <- ASIA_rast_etopo5.sp %>%
    tibble::as_tibble() %>%
    dplyr::mutate(chunk = rep(1:ceiling(nrow(.)/chunk_size),length.out = nrow(.)) %>%
                    sort()
    ) %>%
    dplyr::group_by(chunk) %>%
    dplyr::do(prediction = fields::predict.mKrig(model,
                                                 xnew = .[,c("x","y")],
                                                 Z = .$elevation) %>%
                as.vector()) %$%
    prediction %>%
    unlist()

  return(list(model = model, prediction = prediction))

}

# Calculate gdd kriging models for each crop
gdd_model_files <- paste0(derived_data,"models/",crop_GDD$cultivar,"_models.rds")

if(params$clean){
  unlink(paste0(derived_data,"models/"), recursive = TRUE, force = TRUE)
  dir.create(paste0(derived_data,"models/"), recursive = TRUE, showWarnings = F)
  crop_GDD_run <- crop_GDD
}else{
  dir.create(paste0(derived_data,"models/"), recursive = TRUE, showWarnings = F)
  crop_GDD_run <- crop_GDD[!file.exists(gdd_model_files),]
}

if(nrow(crop_GDD_run) == 0){
  message("All indicator Krige models have already been calculated. Continuing.")
}else message("Calculating indicator Krige models for ",
              nrow(crop_GDD_run),
              " cultivars:\n",
              paste0(capture.output(crop_GDD_run), collapse = "\n"))

# A function to reduce the size of a loess model by half
skinny.loess <- function(x){
  x[c("fitted",
      "residuals",
      "enp",
      "one.delta",
      "two.delta",
      "trace.hat",
      "call",
      "terms",
      "xnames")] <- NULL
  return(x)
}

# A function of correct the indication predictions and estimate a smooth
# monotonic function
# This first uses isotonic regression, then loess smoothing with a degree of 1
smooth.preds <- function(y){
  y[y<0] <- 0
  y[y>1] <- 1
  y <- loess(isoreg(y~sample.points)$yf~sample.points, span=0.1, degree=1) %>%
    skinny.loess()
  return(y)
}

if(nrow(crop_GDD_run) > 0){
  options(dplyr.show_progress = FALSE)
  chunk_size <- 10000

  gdd.models <-
    furrr::future_map(1:nrow(crop_GDD_run),
                      .options = furrr::future_options(
                        packages = c("fields",
                                     "dplyr",
                                     "magrittr",
                                     "readr",
                                     "gutaker2020")),
                      .f = function(crop){

                        # Threshold for indicator kriging
                        GHCN.GDD.incremented.sd[[as.character(crop_GDD_run[crop,"t_base"])]] %>%
                          dplyr::mutate(GDD_thresh = {GDD >= as.numeric(crop_GDD_run[crop,"min_gdd"])}) %>%
                          dplyr::group_by(SD_change) %>%
                          dplyr::do(out_preds = krige_and_predict(.)) %$%
                          out_preds %>%
                          sapply("[[","prediction") %>%
                          apply(1,smooth.preds) %>%
                          tibble::tibble(model = .) %>%
                          maptools::spCbind(ASIA_rast_etopo5.sp,.) %>%
                          readr::write_rds(paste0(derived_data,"models/",
                                                  crop_GDD_run[crop,"cultivar"],
                                                  "_models.rds"), 
                                           compress = "xz")

                        return(paste0(derived_data,"models/",
                                      crop_GDD_run[crop,"cultivar"],
                                      "_models.rds"))

                      })

}

rm(ASIA_rast_etopo5,
   ASIA_rast_etopo5.sp)

message("Calculation indicator Krige models complete: ", capture.output(Sys.time() - time_check))

Predict crop niches through time {-}

## Predicting crop niche from smoothed Krige models
# Calculate niches for each crop using the Marcott et al. 2013.
# create the cluster for parallel computation
message("Generating niche reconstructions")
time_check <-  Sys.time()
if(params$clean){
  unlink(paste0(derived_data,"recons/"), recursive = TRUE, force = TRUE)
}
dir.create(paste0(derived_data,"recons/"), recursive = TRUE, showWarnings = F)

# Create the plan for (possibly parallel) execution
future::plan(future::multisession, 
             workers = max(
               floor(
                 min(
                   future::availableCores(), 
                   params$cores)/2), 
               1),
             gc = TRUE)

furrr::future_walk(crop_GDD$cultivar,    
                   .options = furrr::future_options(
                     packages = c("fields",
                                  "dplyr",
                                  "magrittr",
                                  "sp",
                                  "readr",
                                  "gutaker2020")),
                   function(crop){

                     if(file.exists(paste0(derived_data,"recons/",
                                           crop,
                                           "_recons.rds")))
                       return()

                     if(!file.exists(paste0(derived_data,"models/",
                                            crop,
                                            "_models.rds"))) 
                       stop("Models for ",
                            crop,
                            " are missing! Aborting.")

                     Zs <- c("Z_Lower","Z","Z_Upper")

                     crop.models <- readr::read_rds(paste0(derived_data,"models/",
                                                           crop,
                                                           "_models.rds"))

                     purrr::map(Zs, 
                                 function(z){
                                   suppressWarnings(
                                     crop.models@data %$%
                                       model %>%
                                       purrr::map(
                                         function(x){
                                           x %>%
                                             predict(newdata = marcott2013[[z]]) %>%
                                             magrittr::multiply_by(100) %>%
                                             round() %>%
                                             as.integer()
                                         }) %>%
                                       do.call(rbind, .) %>%
                                       tibble::as_tibble() %>%
                                       magrittr::set_colnames(marcott2013$YearBP) %>%
                                       new("SpatialPointsDataFrame",
                                           data = .,
                                           coords.nrs = [e-mail],
                                           coords = crop.models@coords,
                                           bbox = crop.models@bbox,
                                           proj4string = crop.models@proj4string) %>%
                                       as("SpatialPixelsDataFrame") %>%
                                       raster::brick() %>%
                                       raster::setZ(marcott2013$YearBP,
                                                    name="Years BP")
                                   )
                                 }) %>%
                       magrittr::set_names(Zs) %>%
                       purrr::map(raster::readAll) %>%
                       readr::write_rds(paste0(derived_data,"recons/",
                                               crop,
                                               "_recons.rds"),
                                        compress = "xz",
                                        compression = 9)

                     return()
                   })

message("Generation of niche reconstructions complete: ", capture.output(Sys.time() - time_check))

Combine similar crop niches {-}

## Combining crop niches from similar crops by taking an arithmatic mean
message("Combining like crop niches")
time_check <-  Sys.time()

# Force Raster to load large rasters into memory
raster::rasterOptions(chunksize = 2e+10,
                      maxmemory = 2e+11)

crop_GDD %>%
  split(as.factor(crop_GDD$crop)) %>%
  purrr::walk(function(x){

    if(file.exists(paste0(derived_data,"recons/all_",x$crop[[1]],".rds"))) return()
    n_crops <- length(x$crop)

    if(n_crops == 1){
      file.copy(paste0(derived_data,"recons/",
                       x$cultivar[[1]],
                       "_recons.rds"),
                paste0(derived_data,"recons/all_",
                       x$crop[[1]],
                       ".rds"))
      return()
    }

    x %$%
      cultivar %>%
      purrr::map(function(cultivar){
        paste0(derived_data,"recons/",cultivar,"_recons.rds") %>%
          readr::read_rds()
      }) %>%
      purrr::transpose() %>%
      purrr::map(function(x){
        x %>%
          Reduce(f = "+", x = .) %>%
          magrittr::divide_by(n_crops) %>%
          round() %>%
          raster::readAll() %>%
          magrittr::set_names(names(x[[1]]))
      }) %>%
      readr::write_rds(paste0(derived_data,"recons/all_",
                              x$crop[[1]],
                              ".rds"),
                       compress = "xz",
                       compression = 9)
  })

message("Combining like crop niches complete: ", capture.output(Sys.time() - time_check))

## Plotting cultivar niche
# create the cluster for parallel computation
message("Plotting cultivar niche reconstructions")
time_check <-  Sys.time()

dir.create(paste0(figures,"cultivar_niches/"),
           recursive = TRUE,
           showWarnings = FALSE)

breaks <- seq(0, 100, 10)

pal <- viridisLite::viridis(length(breaks) - 1)

gdd.recons <- 
  purrr::map_chr(1:nrow(crop_GDD),
                 function(n){

                   cultivar <- crop_GDD[n,]$cultivar

                   title <- stringr::str_c(crop_GDD[n,]$cultivar_long,
                                           " \u2014 Required GDD: ",
                                           crop_GDD[n,]$min_gdd,
                                           " at a ",
                                           crop_GDD[n,]$t_base,
                                           "°C base")

                   rasts <- readr::read_rds(paste0(derived_data,"recons/",cultivar,"_recons.rds")) %>%
                     purrr::map(function(x){
                       x %>%
                         magrittr::extract2(which(.@z$`Years BP` > 1000)) %>%
                         # raster:::readAll() %>%
                         magrittr::extract2(raster::nlayers(.):1) 
                     })


                   years <- rasts[[1]] %>%
                     names() %>%
                     gsub(pattern = "X",
                          replacement = "",
                          x = .) %>%
                     as.numeric()

                   if(!file.exists(paste0(figures,"cultivar_niches/",cultivar,".pdf")))
                     gutaker2020:::space_time_plot(
                       the_brick = rasts$Z,
                       the_brick_lower = rasts$Z_Lower,
                       the_brick_upper = rasts$Z_Upper,
                       out_file = paste0(figures,"cultivar_niches/",cultivar,".pdf"),
                       title = title,
                       time = years,
                       timelim = c(max(years),min(years)),
                       timeaxis =  seq(from = max(years)-500,
                                       to = min(years),
                                       by = -500),
                       timelab = "Years BP",
                       zbreaks = breaks,
                       zlab = "Probability of being in niche",
                       zaxis = seq(0,100,10),
                       zcolors = pal
                     )

                   if(!file.exists(paste0(figures,"cultivar_niches/",cultivar,".mp4")))
                     gutaker2020:::space_time_video(                       
                       the_brick = rasts$Z,
                       the_brick_lower = rasts$Z_Lower,
                       the_brick_upper = rasts$Z_Upper,
                       out_file = paste0(figures,"cultivar_niches/",cultivar,".mp4"),
                       title = title,
                       time = years,
                       timelim = c(max(years),min(years)),
                       timeaxis =  seq(from = max(years)-500,
                                       to = min(years),
                                       by = -500),
                       timelab = "Years BP",
                       zbreaks = breaks,
                       zlab = "Probability of being in niche",
                       zaxis = seq(0,100,10),
                       zcolors = pal
                     )

                   return(paste0(figures,"cultivar_niches/",cultivar,".pdf"))

                 })

message("Plotting of cultivar niche reconstructions complete: ", capture.output(Sys.time() - time_check))

Generate niche videos {-}

# Create a static plot for paper publication
# A function to extract data from a crop raster brick
tidyCrop <- function(x, years){
  x %>%
    paste0(derived_data,"recons/all_",.,".rds") %>%
    readr::read_rds() %$%
    Z %>%
    magrittr::extract2(stringr::str_c("X",years)) %>%
    as("SpatialPixelsDataFrame") %>%
    tibble::as_tibble() %>%
    tidyr::gather(Year, Niche, -x:-y) %>%
    dplyr::rename(Longitude = x,
                  Latitude = y) %>%
    dplyr::mutate(Year = stringr::str_remove(Year, "X") %>%
                    factor(levels = years) %>%
                    forcats::fct_relabel(function(x){stringr::str_c(x," cal. BP")}),
                  Crop = x)
}

breaks <- seq(0, 100, 10)

pal <- (RColorBrewer::brewer.pal(10, "Spectral") %>%
          rev() %>%
          colorRampPalette(.))(length(breaks) - 1)

# A function to create a plot of years and crops
facet_niche <- function(crops, years){

  these_rasters <- crops %>%
    purrr::map(tidyCrop, years = years) %>%
    dplyr::bind_rows() %>%
    dplyr::mutate(Crop = dplyr::recode_factor(Crop,
                                              rice = "Rice"),
                  Niche = Niche  %>%
                    cut(breaks = breaks,
                        include.lowest = TRUE)) 

  p <- these_rasters %>%
    ggplot2::ggplot() +
    ggplot2::geom_raster(mapping = ggplot2::aes(x = Longitude,
                                                y = Latitude,
                                                fill = Niche),
                         show.legend = TRUE) +
    ggplot2::coord_quickmap() +
    ggplot2::facet_grid(Crop ~ Year,
                        switch = "y")

  return(p)
}

message("Plotting of crop niche reconstructions complete: ", capture.output(Sys.time() - time_check))

Are GDD differences between Japan and Thailand stable across the reconstruction? {-}

japan_thai <- countries %>%
  dplyr::filter(NAME %in% c("Japan","Thailand"))

# # Get the pixel counts
# crop <- "tropical"
# y.vx <- velox::velox(readr::read_rds(paste0("./data/derived_data/recons/all_",crop,".rds"))[[1]])
# y.vx$extract(japan_thai,
#              fun = function(i){sum(!is.na(i))})[,1]

japan_thai_data <-
  c("Tropical Japonica" = "tropical_japonica",
    "Temperate Japonica" = "temperate_japonica",
    "Indica" = "indica") %>%
  purrr::map(
    function(crop){
      readr::read_rds(paste0(derived_data,"recons/all_",crop,".rds")) %>%
        purrr::map(
          function(y){
            y.vx <- velox::velox(y)
            y.vx$extract(japan_thai,
                         fun = function(i){mean(i, na.rm = TRUE)}) %>%
              t() %>%
              tibble::as_tibble() %>%
              magrittr::set_names(japan_thai$NAME) %>%
              dplyr::mutate(`Years BP` = 
                              names(y) %>% 
                              stringr::str_remove("X") %>% 
                              as.integer()) %>%
              dplyr::select(`Years BP`,
                            dplyr::everything())
          }) 
    }) %>%
  purrr::map(dplyr::bind_rows,
             .id = "Series") %>%
  dplyr::bind_rows(.id = "Crop") %>%
  tidyr::gather(Country, `Probability of being in niche`,-Crop:-`Years BP`) %>%
  tidyr::spread(Series, `Probability of being in niche`)

japan_thai_data %>%
  ggplot2::ggplot(ggplot2::aes(x = `Years BP`,
                               color = Country)) +
  ggplot2::geom_ribbon(ggplot2::aes(ymin = Z_Lower,
                                    ymax = Z_Upper,
                                    fill = Country),
                       color = NA,
                       alpha = 0.3) +
  ggplot2::geom_line(ggplot2::aes(y = Z)) +
  ggplot2::scale_x_reverse(limits = c(5510,1010),
                           breaks = seq(5510,1010,-500),
                           expand = c(0,0)) +
  ggplot2::scale_y_continuous(name = "Probability of being in niche",
                              limits = c(0,100),
                              expand = c(0,0)) +

  ggplot2::facet_wrap(facets = "Crop",
                      ncol = 1) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position="bottom",
                 plot.margin = ggplot2::margin(t = 0, r = 10, b = 0, l = 0, unit = "pt"))

ggplot2::ggsave(paste0(figures,"supplementary_figure_20b.pdf"),
                width = 6.5,
                height = 6,
                device = cairo_pdf)


smooth.preds <- function(y){
  y[y<0] <- 0
  y[y>1] <- 1
  y <- loess(isoreg(y~sample.points)$yf~sample.points, span=0.1, degree=1) %>%
    skinny.loess()
  return(y)
}

GHCN.data.final$spatial %>%
  sf::st_as_sf() %>%
  sf::st_transform(sf::st_crs(japan_thai)) %>%
  sf::st_intersection(japan_thai) %>%
  dplyr::select(ID) %>%
  tibble::as_tibble() %>%
  sf::st_as_sf() %>%
  dplyr::left_join(
    GHCN.GDD.incremented.sd[[1]] %>%
      dplyr::group_by(ID) %>%
      dplyr::summarise(`GDD Model` = list(loess(isoreg(GDD ~ SD_change)$yf ~ sample.points, span=0.1, degree=1)))
  ) %>%
  dplyr::rowwise() %>%
  dplyr::mutate(`GDD` = list(
    marcott2013 %>%
      dplyr::select(`Years BP` = YearBP,
                    Z_Lower,
                    Z,
                    Z_Upper) %>%
      dplyr::mutate_at(.vars = vars(Z_Lower,
                                    Z,
                                    Z_Upper),
                       .funs = function(x){
                         predict(`GDD Model`,newdata = x)
                       })
  )) %>%
  sf::st_as_sf() %>%
  dplyr::select(GDD) %>%
  dplyr::ungroup() %>%
  sf::st_intersection(japan_thai) %>%
  dplyr::rename(Country = NAME) %>%
  sf::st_drop_geometry() %>%
  tidyr::unnest(cols = c(GDD)) %>%
  dplyr::group_by(Country,`Years BP`) %>%
  dplyr::summarise_all(mean, na.rm = TRUE) %>%
  ggplot2::ggplot(ggplot2::aes(x = `Years BP`,
                               color = Country)) +
  ggplot2::geom_ribbon(ggplot2::aes(ymin = Z_Lower,
                                    ymax = Z_Upper,
                                    fill = Country),
                       color = NA,
                       alpha = 0.3) +
  ggplot2::geom_line(ggplot2::aes(y = Z)) +
  ggplot2::scale_x_reverse(limits = c(5510, 1010),
                           breaks = seq(5510, 1010, -500),
                           expand = c(0,0)) +
  ggplot2::scale_y_continuous(name = "Growing Degree Days",
                              limits = c(0,8000),
                              expand = c(0,0)) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position="bottom",
                 plot.margin = ggplot2::margin(t = 10, r = 10, b = 0, l = 0, unit = "pt"))

ggplot2::ggsave(paste0(figures,"supplementary_figure_20a.pdf"),
                width = 6.5,
                height = 6.5/1.618,
                device = cairo_pdf)

Export figure source data for Nature Plants

# Force Raster to load large rasters into memory
raster::rasterOptions(chunksize = 2e+09,
                      maxmemory = 2e+10)

dir.create(paste0(figures,"source_data/"))

crop_GDD %$%
  cultivar %>%
  purrr::map(function(cultivar){
    rasts <- 
      readr::read_rds(paste0(derived_data,"recons/",cultivar,"_recons.rds")) %>%
      purrr::map_dfr(function(x){
        # the_brick <- 
        # x %>%
        # magrittr::extract2(which(.@z$`Years BP` > 1000))

        mean.temporal <- 
          raster::cellStats(x, mean, na.rm = T) %>%
          round(digits = 2)

        ci.temporal <- 
          raster::cellStats(x,
                            stats::quantile,
                            probs = c(0.25, 0.5, 0.75),
                            na.rm = T)

        tibble::tibble(`Years BP` = x@z$`Years BP`,
                       `Spatial Mean` = mean.temporal,
                       `Spatial Lower IQR` = ci.temporal[1,],
                       `Spatial Median` = ci.temporal[2,],
                       `Spatial Upper IQR` = ci.temporal[3,])

      },
      .id = "Marcott et al. 2013 Reconstruction") %>%
      dplyr::mutate(`Marcott et al. 2013 Reconstruction`  = `Marcott et al. 2013 Reconstruction` %>%
                      factor(levels = c("Z","Z_Lower","Z_Upper"),
                             ordered = TRUE) %>%
                      forcats::fct_recode("Mean" = "Z",
                                          "1σ Lower" = "Z_Lower",
                                          "1σ Upper" = "Z_Upper")) %>%
      dplyr::select(`Years BP`, dplyr::everything())

  }) %>%
  magrittr::set_names(crop_GDD$cultivar_long) %>%
  magrittr::extract(.,sort(names(.))) %>%
  writexl::write_xlsx(paste0(figures,"source_data/figure_3e_cultivar_graph_data.xlsx"))

## Export raster data for Figure 3
readr::read_rds(paste0(derived_data,"recons/tropical_japonica_120_recons.rds")) %$%
  Z %>%
  magrittr::extract2(paste0("X", c(4410, 3550))) %>%
  raster::unstack() %>%
  magrittr::set_names(.,c(4410, 3550)) %>%
  purrr::iwalk(~raster::writeRaster(.x, paste0(figures,"source_data/figure_3_raster_data_",.y,".asc")))

Colophon {-}

This report was generated on r Sys.time() using the following computational environment and dependencies:

# which R packages and versions?
devtools::session_info()

The current Git commit details are:

# what commit is this file at? You may need to change the path value
# if your Rmd is not in analysis/paper/
git2r::repository("..")

bocinsky/gutaker2020 documentation built on Sept. 25, 2020, 3:10 a.m.