R/ExtractData_Soils.R

In DrylandEcology/rSFSW2: Simulation Framework for SOILWAT2

#------------------------------------------------------------------------------#
#------EXTRACT SOIL CHARACTERISTICS------

#' Preparations for the extraction of external soil datasets
prepare_ExtractData_Soils <- function(SWRunInformation, sim_size, field_sources,
  field_include, how_determine_sources, sw_input_soillayers, sw_input_soils_use,
  sw_input_soils) {

  sites_soils_source <- get_datasource_masterfield(SWRunInformation,
    field_sources, sim_size, how_determine_sources)

  lvars <- c("density", "sand", "clay", "rock", "carbon")
  nvars <- length(lvars)
  coln <- c("i", "depth", paste0(rep(lvars, SFSW2_glovars[["slyrs_maxN"]]),
    "_L", rep(SFSW2_glovars[["slyrs_ids"]], each = nvars)))
  dtemp <- matrix(NA, nrow = sim_size[["runsN_sites"]],
    ncol = 2 + nvars * SFSW2_glovars[["slyrs_maxN"]],
    dimnames = list(NULL, coln))
  vars <- data.frame(input = c("SoilDepth_cm", "Matricd_L", "Sand_L",
                               "Clay_L", "GravelContent_L", "TOC_GperKG_L"),
                     intern = c("depth", lvars),
                     stringsAsFactors = FALSE)

  do_include <- get_datasource_includefield(SWRunInformation, field_include,
    sim_size)

  list(source = sites_soils_source, data = dtemp, idone = vector(),
    use = sw_input_soils_use, input = sw_input_soils,
    input2 = sw_input_soillayers, cn = coln, vars = vars, nvars = nvars,
    do_include = do_include)
}


update_soils_input <- function(MMC, sim_size, digits = 2, i_Done, ldepths_cm,
  lys, fnames_in) {

  #add data to MMC[["input"]] and set the use flags
  #set and save soil layer structure
  temp <- MMC[["data"]][i_Done, MMC[["vars"]][1, "intern"]]
  irow <- sim_size[["runIDs_sites"]][i_Done]
  MMC[["input2"]][irow, MMC[["vars"]][1, "input"]] <-
    round(temp)
  icol <- grep("depth_L", colnames(MMC[["input2"]]))
  temp <- rep(ldepths_cm[lys], sum(i_Done))
  MMC[["input2"]][irow, icol[lys]] <- matrix(temp, nrow = sum(i_Done),
    ncol = length(lys), byrow = TRUE)
  MMC[["input2"]][irow, icol[-lys]] <- NA
  utils::write.csv(MMC[["input2"]], file = fnames_in[["fslayers"]],
    row.names = FALSE)
  unlink(fnames_in[["fpreprocin"]])

  #set and save soil texture
  for (k in 1 + seq_len(MMC[["nvars"]])) {
    icol <- grep(MMC[["vars"]][k, "input"], names(MMC[["use"]]))
    temp <- MMC[["data"]][i_Done, paste0(MMC[["vars"]][k, "intern"], "_L", lys)]
    if (!all(is.na(temp))) {
      MMC[["input"]][irow, icol[lys]] <- round(temp, digits)
      MMC[["use"]][icol[lys]] <- TRUE
      MMC[["use"]][icol[-lys]] <- FALSE
    }
  }

  #write data to disk
  utils::write.csv(reconstitute_inputfile(MMC[["use"]], MMC[["input"]]),
    file = fnames_in[["fsoils"]], row.names = FALSE)
  unlink(fnames_in[["fpreprocin"]])

  MMC
}

adjust_soils_todos <- function(todos, MMC, sim_size) {
  temp <- is.na(MMC[["input2"]][sim_size[["runIDs_sites"]],
    MMC[["vars"]][1, "input"]])

  for (k in 1 + seq_len(MMC[["nvars"]])) {
    temp <- temp | has_nodata(MMC[["input"]][sim_size[["runIDs_sites"]], ],
      MMC[["vars"]][k, "input"])
  }
  todos <- todos & MMC[["do_include"]] & temp
}



#' \var{\sQuote{CONUS-SOIL}} is a rasterized and controlled
#' \var{\sQuote{STATSGO}} dataset; information for 11 soil are layers available.
#'
#' @param default_TOC_GperKG A numeric value. The default value is
#'   0 g \var{TOC} per kg soil.
#'
#' @references Miller, D. A. and R. A. White. 1998. A conterminous United States
#'  multilayer soil characteristics dataset for regional climate and hydrology
#'  modeling. Earth Interactions 2:1-26.
#' @section Note: It appears that \var{NODATA} values are recorded as 0s.
#' @section Note: We convert what \var{\sQuote{CONUS-SOIL}} calls
#'   \var{\dQuote{bulk density}} to matric density with equation 20 from
#'   Saxton et al. 2006: \eqn{bulkd = matricd * (1 - rockvol) + rockvol * 2.65}
#'   If this variable is indeed \var{\dQuote{bulk density}}, then equation 20
#'   (Saxton et al. 2006) would give negative values
extract_soil_CONUSSOIL <- function(MMC, sim_size, sim_space, dir_ex_soil,
  fnames_in, resume, verbose, default_TOC_GperKG = 0) {

  if (verbose) {
    t1 <- Sys.time()
    temp_call <- shQuote(match.call()[1])
    print(paste0("rSFSW2's ", temp_call, ": started at ", t1))

    on.exit({
      print(paste0("rSFSW2's ", temp_call, ": ended after ",
      round(difftime(Sys.time(), t1, units = "secs"), 2), " s"))
      cat("\n")}, add = TRUE)
  }

  MMC[["idone"]]["CONUSSOIL1"] <- FALSE
  todos <- is.na(MMC[["source"]]) |
    MMC[["source"]] == "CONUSSOILFromSTATSGO_USA"

  if (resume) {
    todos <- adjust_soils_todos(todos, MMC, sim_size)
  }
  names(todos) <- NULL
  n_extract <- sum(todos)

  if (n_extract > 0) {
    if (verbose)
      print(paste("Soil data from 'CONUSSOILFromSTATSGO_USA' will be extracted",
        "for n =", n_extract, "sites"))

    dir.ex.conus <- file.path(dir_ex_soil, "CONUSSoil", "output", "albers")
    stopifnot(file.exists(dir.ex.conus))

    ldepth_CONUS <- c(0, 5, 10, 20, 30, 40, 60, 80, 100, 150, 200, 250)  #in cm
    layer_N <- length(ldepth_CONUS) - 1
    ils <- seq_len(layer_N)

    g <- raster::brick(file.path(dir.ex.conus, "bd.tif"))
    crs_data <- raster::crs(g)

    #locations of simulation runs
    sites_conus <- sim_space[["run_sites"]][todos, ]
    # Align with data crs
    if (!raster::compareCRS(sim_space[["crs_sites"]], crs_data)) {
      sites_conus <- sp::spTransform(sites_conus, CRS = crs_data)
    }

    if (sim_space[["scorp"]] == "point") {
      cell_res_conus <- NULL
      args_extract <- list(y = sites_conus, type = sim_space[["scorp"]])

    } else if (sim_space[["scorp"]] == "cell") {
      cell_res_conus <- align_with_target_res(res_from = sim_space[["sim_res"]],
        crs_from = sim_space[["sim_crs"]],
        sp = sim_space[["run_sites"]][todos, ],
        crs_sp = sim_space[["crs_sites"]], crs_to = crs_data)
      args_extract <- list(y = cell_res_conus, coords = sites_conus,
        method = "block", type = sim_space[["scorp"]])
    }

    #---extract data
    # bulk density -> matric density
    message("NOTE: soil density values extracted from CONUS-soil ",
      "(gridded STATSGO) may be too low!")
    cond30 <- compiler::cmpfun(function(v) ifelse(is.na(v) | v < 30, NA, v))
    ftemp <- file.path(dir.ex.conus, "bd_cond30.tif")
    g <- if (file.exists(ftemp)) {
        raster::brick(ftemp)
      } else {
        # bulk density of less than 0.3 g / cm3 should be treated as no soil
        raster::calc(g, fun = cond30, filename = ftemp)
      }
    temp <- do.call("extract_rSFSW2", args = c(args_extract, x = list(g)))
    MMC[["data"]][todos, grep("density", MMC[["cn"]])[ils]] <- temp / 100

    # soil depth
    # depth in cm >< bedrock from datafile.bedrock, but seems to make more
    # sense?
    cond0 <- compiler::cmpfun(function(v) ifelse(!is.na(v) & v > 0, v, NA))
    ftemp <- file.path(dir.ex.conus, "rockdepm_cond0.tif")
    g <- if (file.exists(ftemp)) {
        raster::raster(ftemp)
      } else {
        # rockdepth of 0 cm should be treated as no soil
        raster::calc(raster::raster(file.path(dir.ex.conus, "rockdepm.tif")),
          fun = cond0, filename = ftemp)
      }
    rockdep_cm <- do.call("extract_rSFSW2", args = c(args_extract, x = list(g)))

    # rock volume: new with v31: rockvol -> gravel vol%
    g <- raster::brick(file.path(dir.ex.conus, "rockvol.tif"))
    temp <- do.call("extract_rSFSW2", args = c(args_extract, x = list(g)))
    temp <- ifelse(is.finite(temp), temp, NA)
    # eq. 7 of Miller et al. 1998
    temp <- pmax(pmin(temp / 100, 1), 0) # volume fraction of bulk = total soil

    # adjust soil depth by layers with 100% rock volume
    solid_rock_nl <- apply(temp >= 1 - SFSW2_glovars[["toln"]], 1, sum,
      na.rm = TRUE)
    solid_rock_nl <- 1 + layer_N - solid_rock_nl
    solid_rock_cm <- ldepth_CONUS[solid_rock_nl]
    MMC[["data"]][todos, grep("rock", MMC[["cn"]])[ils]] <- temp
    MMC[["data"]][todos, "depth"] <- pmin(rockdep_cm, solid_rock_cm)

    lys <- seq_len(max(findInterval(MMC[["data"]][todos, "depth"],
      ldepth_CONUS[-1]), na.rm = TRUE))

    # sand, silt, and clay
    ftemp <- file.path(dir.ex.conus, "sand_cond0.tif")
    g <- if (file.exists(ftemp)) {
        raster::brick(ftemp)
      } else {
        raster::calc(raster::brick(file.path(dir.ex.conus, "sand.tif")),
          fun = cond0, filename = ftemp)
      }
    sand <- do.call("extract_rSFSW2", args = c(args_extract, x = list(g)))

    ftemp <- file.path(dir.ex.conus, "clay_cond0.tif")
    g <- if (file.exists(ftemp)) {
        raster::brick(ftemp)
      } else {
        raster::calc(raster::brick(file.path(dir.ex.conus, "clay.tif")),
          fun = cond0, filename = ftemp)
      }
    clay <- do.call("extract_rSFSW2", args = c(args_extract, x = list(g)))

    ftemp <- file.path(dir.ex.conus, "silt_cond0.tif")
    g <- if (file.exists(ftemp)) {
        raster::brick(ftemp)
      } else {
        raster::calc(raster::brick(file.path(dir.ex.conus, "silt.tif")),
          fun = cond0, filename = ftemp)
      }
    silt <- do.call("extract_rSFSW2", args = c(args_extract, x = list(g)))

    #Normalize to 0-1
    total_matric <- sand + clay + silt # values between 0.99 and 1.01
    total_matric[!is.finite(total_matric)] <- NA
    MMC[["data"]][todos, grep("sand", MMC[["cn"]])[ils]] <- sand / total_matric
    MMC[["data"]][todos, grep("clay", MMC[["cn"]])[ils]] <- clay / total_matric

    # There is no organic carbon data, set all values to a default
    MMC[["data"]][todos, grep("carbon", MMC[["cn"]])[ils]] <- default_TOC_GperKG

    # Determine successful extractions
    MMC[["idone"]]["CONUSSOIL1"] <- TRUE
    i_good <- stats::complete.cases(MMC[["data"]][todos, "depth"])
    MMC[["source"]][which(todos)[!i_good]] <- NA

    if (any(i_good)) {
      i_Done <- rep(FALSE, times = sim_size[["runsN_sites"]])
      # length(i_Done) == length(runIDs_sites) == runsN_sites
      i_Done[which(todos)[i_good]] <- TRUE

      MMC[["source"]][i_Done] <- "CONUSSOILFromSTATSGO_USA"
      MMC <- update_soils_input(MMC, sim_size, digits = 2, i_Done,
        ldepths_cm = ldepth_CONUS[-1], lys, fnames_in)
    }

    if (verbose)
      print(paste("Soil data from 'CONUSSOILFromSTATSGO_USA' was extracted ",
        "for n =", sum(i_good), "out of", n_extract, "sites"))
  }

  MMC
}


#' A wrapper for \code{reaggregate_raster} designed to work with the rasters of
#' the \var{\sQuote{ISRIC-WISE}} datasets versions 5-arcmin v1.2 and 30-arcsec
#' v1.0
#'
#' @param i An integer value. The index to select a location from among
#'   \code{sp_sites} and the corresponding resolution \code{res}.
#' @param res A numeric vector of length two or a matrix with two columns. The
#'   x- and y-extent of the rectangle(s) for which to extract values.
#' @param grid A \code{\link[raster:RasterLayer-class]{raster::RasterLayer}}
#'   object with one layer. The raster from which values are extracted.
#' @param sp_sites A \code{\link[sp:SpatialPoints-class]{sp::SpatialPoints}}
#'   object. This object is used to extract the coordinates of the i-th
#'   location.
#' @param att A character string. Which variable in the \var{RAT} table should
#'   be returned. If \code{NULL} then extracted values of \code{grid} are
#'   returned.
#'
#' @seealso \code{\link{reaggregate_raster}}
#'
#' @return A list with four elements \describe{
#'   \item{i}{An integer value. The location index.}
#'   \item{SUIDs_N}{An integer vector. The number of unique values within
#'     the rectangle of \code{x}.}
#'   \item{SUID}{A numeric vector or a vector of character strings. The unique
#'     soil soil identifier values or \code{NA}.}
#'   \item{fraction}{A numeric vector. The relative areas covered by
#'     \code{SUID}.} }
ISRICWISE_extract_SUIDs <- function(i, res = c(0, 0), grid, sp_sites,
  att = NULL) {

  out <- try(reaggregate_raster(x = grid,
        coords = sp::coordinates(sp_sites[i, ]),
        to_res = if (is.null(dim(res))) res else res[i, ],
        with_weights = TRUE,
        method = "block"))

  if (inherits(out, "try-error")) {
    print(out)
    list(i = i, SUIDs_N = -1, SUID = NULL, fraction = NULL)

  } else {
    suids <- if (is.null(att)) {
        temp <- out[[1]][["values"]][[1]]
        ifelse(temp < SFSW2_glovars[["tol"]], NA, temp)
    } else {
        temp <- raster::factorValues(grid, out[[1]][["values"]][[1]], att = att)
        temp <- as.character(unlist(temp))
        ifelse(nchar(temp) == 0L, NA, temp)
    }
    list(i = i, SUIDs_N = out[[1]][["N"]][[1]],
          SUID = suids,
          fraction = out[[1]][["fraction"]][[1]])
  }
}

ISRICWISE_get_prids <- function(suid, dat_wise, layer_N, colname_suid) {
  # If is.na(suid) then 'soils' is a 0-row data.frame
  soils <- if (is.na(suid)) {
      dat_wise[0, ]
    } else {
      dat_wise[dat_wise[, colname_suid] == suid, ]
    }
  frac <- unique(soils[, c("PROP", "PRID")])
  depth <- tapply(soils[, "BotDep"], soils[, "PRID"], max)
  idepth <- depth[match(frac$PRID, names(depth))]

  list(PRIDs_N = nrow(soils) / layer_N,
     PRID = frac$PRID,
     fraction = frac$PROP / 100,
     depth = ifelse(idepth > 0, idepth, NA),
     soildat = soils)
}


extract_soillayer_ISRICWISE <- function(dat, soildat_rows, frac) {
  dat_add <- ifelse(is.na(soildat_rows) | soildat_rows < 0, NA, soildat_rows)
  res_frac <- ifelse(is.na(dat_add) | is.na(frac), 0, frac)

  # NAs do not propagate because they are down-weighted by res_frac
  list(
    # weighted mean = sum of values x weights
    soil = sum(dat, dat_add * res_frac, na.rm = TRUE),
    frac = sum(res_frac, na.rm = TRUE))
}


max_depth_byPRIDs <- function(this_soil, var_ids, val_rocks) {
  if (is.null(val_rocks) || length(val_rocks) == 0L) {
    this_soil[["depth"]]

  } else {
    depths <- by(this_soil[["soildat"]][, c("BotDep", var_ids)],
        INDICES = this_soil[["soildat"]][, "PRID"],
        function(x1) {
          temp <- sapply(val_rocks, function(x2) {
            apply(x1[, -1] - x2, 1, function(x3)
              any(abs(x3) < SFSW2_glovars[["tol"]]))})
          itemp <- which(!apply(temp, 1, any))
          if (length(itemp) > 0) x1[itemp[length(itemp)], "BotDep"] else 0
      })
    depths[match(dimnames(depths)[[1]], this_soil[["PRID"]], nomatch = 0)]
  }
}

#' Calculate weighted mean soil variables from one of the
#' \var{\sQuote{ISRIC-WISE}} data bases for one simulation cell
#'
#' Areas with no soil are 'removed' and the remaining data 'scaled' to 1.
#'
#' @param i An integer value. The number of the simulation site/cell location.
#' @param i_sim_cells_SUIDs A named numeric vector. A row of of the data.frame
#'   returned by \code{\link{ISRICWISE_extract_SUIDs}}.
#' @param sim_soils A named numeric vector. First element is 'i' and second is
#'   'depth'. The following elements represent the soil variables (currently,
#'   "density", "sand", "clay", "rock", "carbon"; see
#'   \code{\link{prepare_ExtractData_Soils}}) for each of the maximally possible
#'   soil layers (see \code{SFSW2_glovars[["slyrs_maxN"]]}). Input is an empty
#'   template.
#' @param layer_N An integer value. The number of soil layers in the dataset
#'   (i.e., 5, see \code{\link{extract_soil_ISRICWISE}}).
#' @param layer_Nsim An integer value. The number of soil layers of a
#'   \pkg{rSFSW2} project representing the dataset (i.e., 6, see
#'   \code{\link{extract_soil_ISRICWISE}}).
#' @param ldepth An integer vector. The depth limits of the extracted
#'   \code{rSFSW2} project soil layers including zero where \code{layer_Nsim + 1
#'   == length(ldepth)}.
#' @param dat_wise A data.frame representing the main \var{\sQuote{ISRIC-WISE}}
#'   database.
#' @param nvars An integer value. The number of soil variables extracted from a
#'   \var{\sQuote{ISRIC-WISE}} dataset (currently, 5; see
#'   \code{\link{prepare_ExtractData_Soils}}).
#' @param var_tags A vector of character string. The column names for
#'   \var{\dQuote{suid}}, bulk density (\var{\dQuote{density}}),
#'   \var{\dQuote{sand}}, \var{\dQuote{clay}}, coarse fragments
#'   (\var{\dQuote{rock}}), and organic carbon content (\var{\dQuote{carbon}}).
#' @param val_rocks An integer vector. The (negative) values which correspond to
#'   rocky outcrops and/or rocky subsoils. Soil depth is limited to above layers
#'   that have values in \code{val_rocks}.
#'
calc_cell_ISRICWISE <- function(i, i_sim_cells_SUIDs, sim_soils, layer_N,
  layer_Nsim, ldepth, dat_wise, nvars, var_tags, val_rocks = NULL) {

  sim_soils["i"] <- i

  #Do calculations if any soils in this simulation cell
  if (i_sim_cells_SUIDs["SUIDs_N"] > 0) {
    soil_labels <- names(sim_soils)
    temp <- grep("suid", names(var_tags), invert = TRUE)
    var_labels <- names(var_tags)[temp]
    var_ids <- unlist(var_tags)[temp]

    # Init
    # sim_frac = fraction of how much this simulation cell is covered with
    #   suids and prids that have a depth > 0 cm
    sim_frac <- 0
    # simlyr_frac = fraction of each soil layer x variable covering this
    #   simulation cell
    simlyr_frac <- rep(0, times = length(sim_soils))
    PRIDs_frac <- NULL
    sim_soils["depth"] <- 0

    this_simCell <- c(as.list(i_sim_cells_SUIDs),
      soils = list(t(sapply(i_sim_cells_SUIDs[["SUID"]], ISRICWISE_get_prids,
      dat_wise = dat_wise, layer_N = layer_N,
      colname_suid = var_tags[["suid"]]))))

    # Loop through the suids within this simulation cell; each suid may be
    #   composed of several prids
    for (k in seq_len(this_simCell[["SUIDs_N"]])) {
      this_soil <- this_simCell[["soils"]][k, ]

      if (this_soil[["PRIDs_N"]] > 0) {
        # Vector of the fractions of each prid in relation to the simulation
        #   cell
        prids_frac <- this_soil[["fraction"]] * this_simCell[["fraction"]][k]
        PRIDs_frac <- c(PRIDs_frac, prids_frac)
        temp <- sum(ifelse(is.na(this_soil[["depth"]]), 0, prids_frac))
        sim_frac <- sim_frac + temp
        temp <- sum(
          max_depth_byPRIDs(this_soil, var_ids, val_rocks) * prids_frac,
          na.rm = TRUE)
        sim_soils["depth"] <- sim_soils["depth"] + temp

        if (!all(is.na(this_soil[["depth"]])))
          for (ils in seq_len(layer_Nsim)) {
          # Split wise soil layer 0-20 cm into two layers, 0-10 and 10-20 cm,
          #   to account for lithosols
          lwise <- if (ils == 1) 1 else ils - 1
          # Checks if for each prid, the soils are deeper than this layer.
          #   It also accounts that soil depth for Rock outcrops (RK) is set
          #   to 0 instead of < 0 for such as water and glaciers.
          #   Lithosols (Ix) have depth of 10 cm.
          layer.there <- this_soil[["depth"]] > ldepth[ils]
          pfracl <- prids_frac[layer.there]

          if (sum(layer.there, na.rm = TRUE) > 0) {
            irow <- this_soil[["soildat"]][, "Layer"] == paste0("D", lwise)

            if (sum(irow) > 0) for (iv in seq_along(var_labels)) {
              ids <- which(soil_labels == paste0(var_labels[iv], "_L", ils))
              temp <- extract_soillayer_ISRICWISE(dat = sim_soils[ids],
                soildat_rows = this_soil[["soildat"]][irow, var_ids[iv]],
                frac = pfracl)
              sim_soils[ids] <- temp[["soil"]]
              simlyr_frac[ids] <- simlyr_frac[ids] + temp[["frac"]]
            }
          }
        }
      }
    }

    #Adjust values for area present
    fracs <- c(1, sim_frac, simlyr_frac[- (1:2)])
    sim_soils <- sim_soils / fracs
  }

  sim_soils
}

#' Wrapping \code{\link{calc_cell_ISRICWISE}} in \code{\link[base]{try}} to
#' handle errors without stopping a \pkg{rSFSW2} run.
try_cell_ISRICWISE <- function(i, sim_cells_SUIDs, template_simulationSoils,
  layer_N, layer_Nsim, ldepth, dat_wise = dat_wise, nvars, var_tags,
  val_rocks) {

  if (i %% 1000 == 0)
    print(paste(Sys.time(), "'try_cell_ISRICWISE'",
      "done:", i))

  temp <- try(calc_cell_ISRICWISE(i,
    i_sim_cells_SUIDs = sim_cells_SUIDs[i, ],
    sim_soils = template_simulationSoils, layer_N = layer_N,
    layer_Nsim = layer_Nsim, ldepth = ldepth, dat_wise = dat_wise,
    nvars = nvars, var_tags = var_tags, val_rocks = val_rocks))

  if (inherits(temp, "try-error")) template_simulationSoils else temp
}


#' Extract soil data from one of the \var{\sQuote{ISRIC-WISE}} datasets
#'
#' @param dataset A character string. Identifies the \var{\dQuote{ISRIC-WISE}}
#'   dataset from which soil data are extracted. See details.
#'
#' @section Details: \code{dataset} is current implemented for values:
#'   \describe{
#'   \item{"ISRICWISEv12"}{Dataset \var{\sQuote{ISRIC-WISE 5-arcmin v1.2}}
#'     (Batjes 2012)}
#'   \item{"ISRICWISE30secV1a"}{Dataset \var{\sQuote{ISRIC-WISE 30-arcsec v1.0}}
#'     (Batjes 2015, 2016)}
#'  }
#'
#' @references Batjes, N. H. 2012. ISRIC-WISE derived soil properties on a
#'  5 by 5 arc-minutes global grid (version 1.2). Report 2012/01 (with data set,
#'  available at \url{www.isric.org}). ISRIC-World Soil Information, Wageningen,
#'  The Netherlands. \url{http://library.wur.nl/WebQuery/wurpubs/443845}
#' @references Batjes N.H. 2016. Harmonised soil property values for broad-scale
#'   modelling (WISE30sec) with estimates of global soil carbon stocks.
#'   Geoderma 269, 61-68
#'   (\url{http://dx.doi.org/10.1016/j.geoderma.2016.01.034}).
#' @references Batjes N.H. 2015. World soil property estimates for broad-scale
#'   modelling (WISE30sec, version 1.0). Report 2015/01, ISRIC-World Soil
#'   Information, Wageningen [available at ISRIC Soil Data Hub],
#'   with addendum and corrigendum.
#'
#' @section Note: Cells with negative soil values indicate non-soils and are
#'   eliminated before aggregations: \enumerate{
#'   \item ISRIC-WISE 5-arcmin v1.2 (Batjes 2012) \itemize{
#'       \item -1: Oceans and inland waters (\var{SUIDS} %in% c(0, 1972, 6997))
#'       \item -2: Glaciers and snow caps (\var{SUIDS} %in% 6998)
#'       \item -3: No/insufficient data
#'       \item -7: Rock outcrops (or shallow subsoils) (\var{SUIDS} %in% 6694)
#'     }
#'   \item ISRIC-WISE 30-arcsec v1.0 (Batjes 2015) \itemize{
#'       \item -1: Oceans and inland waters
#'       \item -2: Glaciers and snow caps
#'       \item -3: Rock outcrops
#'       \item -4: Dunes/Shifting sands
#'       \item -5: Salt flats
#'       \item -7: Rocky subsoils as for Leptosols
#'       \item -9: Remaining miscellaneous units
#'     }
#'   }
extract_soil_ISRICWISE <- function(MMC, sim_size, sim_space,
  dir_ex_soil, fnames_in, dataset = c("ISRICWISEv12", "ISRICWISE30secV1a"),
  resume, verbose) {

  dataset <- match.arg(dataset)

  if (verbose) {
    t1 <- Sys.time()
    temp_call <- shQuote(match.call()[1])
    print(paste0("rSFSW2's ", temp_call, ": started at ", t1, " for dataset ",
      shQuote(dataset)))

    on.exit(enable_debug_dump(file_tag = match.call()[[1]]), add = FALSE)

    on.exit({
      print(paste0("rSFSW2's ", temp_call, ": ended after ",
      round(difftime(Sys.time(), t1, units = "secs"), 2), " s"))
      cat("\n")}, add = TRUE)
  }

  MMC[["idone"]][dataset] <- FALSE
  todos <- is.na(MMC[["source"]]) |
    MMC[["source"]] == paste0(dataset, "_Global")

  if (resume) {
    todos <- adjust_soils_todos(todos, MMC, sim_size)
  }
  names(todos) <- NULL
  n_extract <- sum(todos)

  if (n_extract > 0) {
    if (verbose)
      print(paste("Soil data from", shQuote(dataset),
        "will be extracted for n =", n_extract, "sites"))

    ldepth_WISE <- if (dataset == "ISRICWISEv12") {
        c(0, 10, 20, 40, 60, 80, 100)  #in cm
      } else if (dataset == "ISRICWISE30secV1a") {
        c(0, 10, 20, 40, 60, 80, 100, 150, 200) # in cm
      }

    # WISE contains five soil layers for each prid; I added one layer to account
    #   for lithosols (Ix), which have a depth of 10 cm; for all other soil
    #   types, my layers 0-10 cm and 10-20 cm contain the same wise information
    layer_Nsim <- length(ldepth_WISE) - 1L
    layer_N <- layer_Nsim - 1L

    #run_sites_wise of simulation runs
    run_sites_wise <- sim_space[["run_sites"]][todos, ]
    is_ToDo <- seq_along(run_sites_wise)

    #---extract data
    if (dataset == "ISRICWISEv12") {
      rat_att <- NULL
      var_tags <- list(suid = "SUID", density = "BULK", sand = "SDTO",
        clay = "CLPC", rock = "CFRAG", carbon = "TOTC")
      val_rocks <- -7

      dir.ex.dat <- file.path(dir_ex_soil, "WISE", "wise5by5min_v1b")
      fwise_grid <- file.path(dir.ex.dat, "Grid", "smw5by5min")
      fwise_table <- file.path(dir.ex.dat, "WISEsummaryFile.csv")

    } else if (dataset == "ISRICWISE30secV1a") {
      rat_att <- "NEWSUID"
      var_tags <- list(suid = "NEWSUID", density = "BULK", sand = "SDTO",
        clay = "CLPC", rock = "CFRAG", carbon = "ORGC")
      val_rocks <- c(-3, -7)

      dir.ex.dat <- file.path(dir_ex_soil, "WISE", "WISE30sec_v1a")
      fwise_grid <- file.path(dir.ex.dat, "GISfiles", "wise30sec_fin")
      fwise_table <- file.path(dir.ex.dat, "Interchangeable_format",
        "HW30s_FULL.txt")
    }

    stopifnot(file.exists(fwise_grid), file.exists(fwise_table))
    grid_wise <- raster::raster(fwise_grid)
    dat_wise <- utils::read.csv(file = fwise_table, stringsAsFactors = FALSE)

    #- List all the wise cells that are covered by the grid cell or point
    #  location
    if (verbose) {
      print(paste0("rSFSW2's ", temp_call, " for dataset ", shQuote(dataset),
        ": extract WISE dataset identifiers for simulation cell/point ",
        "locations"))
    }

    if (sim_space[["scorp"]] == "point") {
      cell_res_wise <- NULL
      suids <- raster::extract(grid_wise, run_sites_wise)
      if (!is.null(rat_att)) {
        suids <- as.character(raster::factorValues(grid_wise, suids,
          att = rat_att)[[1]])
      }
      sim_cells_SUIDs <- data.frame(i = is_ToDo, SUIDs_N = 1, SUID = suids,
        fraction = 1, stringsAsFactors = FALSE)

    } else if (sim_space[["scorp"]] == "cell") {
      cell_res_wise <- align_with_target_res(res_from = sim_space[["sim_res"]],
        crs_from = sim_space[["sim_crs"]], sp = run_sites_wise,
        crs_sp = sim_space[["crs_sites"]], crs_to = raster::crs(grid_wise))

      if (SFSW2_glovars[["p_has"]]) {

        #call the simulations depending on parallel backend
        if (identical(SFSW2_glovars[["p_type"]], "mpi")) {

          sim_cells_SUIDs <- Rmpi::mpi.applyLB(is_ToDo, ISRICWISE_extract_SUIDs,
            res = cell_res_wise, grid = grid_wise, sp_sites = run_sites_wise,
            att = rat_att)

        } else if (identical(SFSW2_glovars[["p_type"]], "socket")) {
          sim_cells_SUIDs <- parallel::clusterApplyLB(SFSW2_glovars[["p_cl"]],
            x = is_ToDo, fun = ISRICWISE_extract_SUIDs, res = cell_res_wise,
            grid = grid_wise, sp_sites = run_sites_wise, att = rat_att)

        } else {
          sim_cells_SUIDs <- data.frame(i = is_ToDo, SUIDs_N = 0, SUID = NA,
            fraction = 1)
        }

        clean_SFSW2_cluster()

      } else {
        sim_cells_SUIDs <- lapply(is_ToDo, FUN = ISRICWISE_extract_SUIDs,
          res = cell_res_wise, grid = grid_wise, sp_sites = run_sites_wise,
          att = rat_att)
      }

      sim_cells_SUIDs <- do.call(rbind, sim_cells_SUIDs)
    }

    if ("i" %in% dimnames(sim_cells_SUIDs)[[1]]) {
      sim_cells_SUIDs <- t(sim_cells_SUIDs)
    }
    stopifnot("i" %in% dimnames(sim_cells_SUIDs)[[2]])
    sim_cells_SUIDs <- sim_cells_SUIDs[order(unlist(sim_cells_SUIDs[, "i"])), ]

    #- Calculate simulation cell wide weighted values based on each PRID
    #   weighted by SUID.fraction x PRIP.PROP
    if (verbose) {
      print(paste0("rSFSW2's ", temp_call, " for dataset ", shQuote(dataset),
        ": calculate simulation cell/point location weighted values of soil ",
        "characteristics"))
    }

    template_simulationSoils <- unlist(MMC[["data"]][1, ])
    template_simulationSoils[] <- NA
    template_simulationSoils["depth"] <- 0

    if (SFSW2_glovars[["p_has"]]) {
      #call the simulations depending on parallel backend
      if (identical(SFSW2_glovars[["p_type"]], "mpi")) {

        ws <- Rmpi::mpi.applyLB(is_ToDo,
          try_cell_ISRICWISE,
          sim_cells_SUIDs = sim_cells_SUIDs,
          template_simulationSoils = template_simulationSoils,
          layer_N = layer_N, layer_Nsim = layer_Nsim, ldepth = ldepth_WISE,
          dat_wise = dat_wise, nvars = MMC[["nvars"]], var_tags = var_tags,
          val_rocks = val_rocks)

      } else if (identical(SFSW2_glovars[["p_type"]], "socket")) {

        ws <- parallel::clusterApplyLB(SFSW2_glovars[["p_cl"]], x = is_ToDo,
          fun = try_cell_ISRICWISE,
          sim_cells_SUIDs = sim_cells_SUIDs,
          template_simulationSoils = template_simulationSoils,
          layer_N = layer_N, layer_Nsim = layer_Nsim, ldepth = ldepth_WISE,
          dat_wise = dat_wise, nvars = MMC[["nvars"]], var_tags = var_tags,
          val_rocks = val_rocks)
      }

      clean_SFSW2_cluster()

    } else {
      ws <- lapply(is_ToDo, FUN = try_cell_ISRICWISE,
        sim_cells_SUIDs = sim_cells_SUIDs,
        template_simulationSoils = template_simulationSoils,
        layer_N = layer_N, layer_Nsim = layer_Nsim, ldepth = ldepth_WISE,
        dat_wise = dat_wise, nvars = MMC[["nvars"]], var_tags = var_tags,
        val_rocks = val_rocks)
    }

    ws <- do.call(rbind, ws)
    ws <- ws[order(ws[, "i"]), , drop = FALSE]
    # convert percent to fraction
    icol <- grepl("rock", colnames(ws)) |
      grepl("sand", colnames(ws)) | grepl("clay", colnames(ws))
    ws[, icol] <- ws[, icol] / 100

    # Determine successful extractions
    MMC[["idone"]][dataset] <- TRUE

    i_good <- rep(FALSE, n_extract)
    ids <- seq_len(2 + layer_Nsim * MMC[["nvars"]])
    # i is index for todos
    i_good[ws[stats::complete.cases(ws[, ids, drop = FALSE]), "i"]] <- TRUE
    MMC[["source"]][which(todos)[!i_good]] <- NA

    if (verbose) {
      print(paste0("rSFSW2's ", temp_call, " for dataset ", shQuote(dataset),
        ": soil data was extracted for n = ", sum(i_good), " out of ",
        n_extract, " sites"))
    }

    if (any(i_good)) {
      i_Done <- rep(FALSE, times = sim_size[["runsN_sites"]])
      # length(i_Done) == length(runIDs_sites) == runsN_sites
      i_Done[which(todos)[i_good]] <- TRUE
      MMC[["data"]][todos, seq_len(dim(ws)[2])] <- ws

      MMC[["source"]][i_Done] <- paste0(dataset, "_Global")

      if (verbose) {
        print(paste0("rSFSW2's ", temp_call, " for dataset ", shQuote(dataset),
          ": update soil input files"))
      }
      MMC <- update_soils_input(MMC, sim_size, digits = 2, i_Done,
        ldepths_cm = ldepth_WISE[-1], lys = seq_len(layer_Nsim), fnames_in)
    }
  }

  if (verbose) {
    # Remove debug dumping but not other 'on.exit' expressions before returning
    # without error
    oe <- sys.on.exit()
    oe <- remove_from_onexit_expression(oe, "enable_debug_dump")
    on.exit(eval(oe), add = FALSE)
  }

  MMC
}



#' Extract soil characteristics
#' @export
ExtractData_Soils <- function(exinfo, SFSW2_prj_meta, SFSW2_prj_inputs,
  opt_parallel, resume, verbose = FALSE) {

  field_sources <- "SoilTexture_source"
  field_include <- "Include_YN_SoilSources"

  #--- SET UP PARALLELIZATION
  setup_SFSW2_cluster(opt_parallel,
    dir_out = SFSW2_prj_meta[["project_paths"]][["dir_prj"]],
    verbose = opt_verbosity[["verbose"]],
    print.debug = opt_verbosity[["print.debug"]])
  on.exit(exit_SFSW2_cluster(verbose = opt_verbosity[["verbose"]]),
    add = TRUE)
  on.exit(set_full_RNG(SFSW2_prj_meta[["rng_specs"]][["seed_prev"]],
    kind = SFSW2_prj_meta[["rng_specs"]][["RNGkind_prev"]][1],
    normal.kind = SFSW2_prj_meta[["rng_specs"]][["RNGkind_prev"]][2]),
    add = TRUE)


  MMC <- prepare_ExtractData_Soils(SFSW2_prj_inputs[["SWRunInformation"]],
    sim_size = SFSW2_prj_meta[["sim_size"]], field_sources = field_sources,
    field_include = field_include,
    how_determine_sources =
      SFSW2_prj_meta[["opt_input"]][["how_determine_sources"]],
    sw_input_soillayers = SFSW2_prj_inputs[["sw_input_soillayers"]],
    sw_input_soils_use = SFSW2_prj_inputs[["sw_input_soils_use"]],
    sw_input_soils = SFSW2_prj_inputs[["sw_input_soils"]])

  if (exinfo$ExtractSoilDataFromCONUSSOILFromSTATSGO_USA) {
    MMC <- extract_soil_CONUSSOIL(MMC,
      sim_size = SFSW2_prj_meta[["sim_size"]],
      sim_space = SFSW2_prj_meta[["sim_space"]],
      dir_ex_soil = SFSW2_prj_meta[["project_paths"]][["dir_ex_soil"]],
      fnames_in = SFSW2_prj_meta[["fnames_in"]], resume, verbose)
  }

  if (exinfo$ExtractSoilDataFromISRICWISE30secV1a_Global) {
    MMC <- extract_soil_ISRICWISE(MMC,
      sim_size = SFSW2_prj_meta[["sim_size"]],
      sim_space = SFSW2_prj_meta[["sim_space"]],
      dir_ex_soil = SFSW2_prj_meta[["project_paths"]][["dir_ex_soil"]],
      fnames_in = SFSW2_prj_meta[["fnames_in"]], dataset = "ISRICWISE30secV1a",
      resume, verbose)
  }

  if (exinfo$ExtractSoilDataFromISRICWISEv12_Global) {
    MMC <- extract_soil_ISRICWISE(MMC,
      sim_size = SFSW2_prj_meta[["sim_size"]],
      sim_space = SFSW2_prj_meta[["sim_space"]],
      dir_ex_soil = SFSW2_prj_meta[["project_paths"]][["dir_ex_soil"]],
      fnames_in = SFSW2_prj_meta[["fnames_in"]], dataset = "ISRICWISEv12",
      resume, verbose)
  }

  SFSW2_prj_inputs[["SWRunInformation"]] <- update_datasource_masterfield(MMC,
    sim_size = SFSW2_prj_meta[["sim_size"]],
    SFSW2_prj_inputs[["SWRunInformation"]],
    SFSW2_prj_meta[["fnames_in"]], field_sources, field_include)

  SFSW2_prj_inputs[["sw_input_soillayers"]] <- MMC[["input2"]]
  SFSW2_prj_inputs[["sw_input_soils_use"]] <- MMC[["use"]]
  SFSW2_prj_inputs[["sw_input_soils"]] <- MMC[["input"]]

  oe <- sys.on.exit()
  oe <- remove_from_onexit_expression(oe, "exit_SFSW2_cluster")
  on.exit(eval(oe), add = FALSE)

  SFSW2_prj_inputs
}


#------END OF SOIL CHARACTERISTICS------
#-----------------------------------------------------------------------------#