R/soilgrids.R

In bmcnellis/RSFIA: FIA Tree Mortality Analysis

#' @include generics.R
#' @include nc_helper.R
#'
#' @title Class, methods, and functions for pulling/analyzing soilgrids data
#'
#' @author Brandon McNellis
#'
#' @name soilgrids
#' @rdname soilgrids
NULL
#'
#' An S4 class for soilgrids data
#'
#' @slot base_url Length 4 character, requires 2 NAs.
#' @slot var_set 1-length character identifying which variable set is present.
#'
#' @rdname soilgrids
soilgrids <- setClass(
  'soilgrids',
  slots = list(
    base_url = 'character',
    var_set = 'character'
  ),
  contains = 'nc_helper'
)
#' @export
setValidity('soilgrids', function(object) {
  errors <- character()

  # base_url
  bt <- object@base_url
  if (table(is.na(bt))['TRUE'] != 2) {
    msg <- paste0('Missing two NAs from base_url that refer to lon/lat')
    errors <- c(errors, msg)
  }
  if (!is.character(bt)) {
    msg <- paste0('base_url must be a character vector')
    errors <- c(errors, msg)
  }

  # var_set
  if (object@var_set != 'major') {
    msg <- paste0('Bad variable set input')
    errors <- c(errors, msg)
  }

  # returns
  if (length(errors) == 0) {
    TRUE
  } else {
    errors
  }
})
#' @rdname soilgrids
#' @export
setMethod('initialize',
          signature(.Object = 'soilgrids'),
          function (.Object, ...) {
            params <- list(...)

            # base_url
            if ('base_url' %in% names(params)) {
              .Object@base_url <- params$base_url
            } else {
              .Object@base_url <- c('https://rest.soilgrids.org/query?lon=', NA, '&lat=', NA)
            }

            # var_set
            .Object@var_set <- 'major'

            # variables
            meta <- SoilGridMeta()
            if (.Object@var_set == 'major') {
              vs <- meta$vars[meta$major_var]
              .Object@variables <- vs
            }

            # time_units
            .Object@time_units <- 'depth'

            # time
            .Object@time <- c(0, 5, 15, 30, 60, 100, 200)

            # returns
            .Object <- callNextMethod()
            mt <- validObject(.Object)
            if (isTRUE(mt)) {
              return(.Object)
            } else {
              return(mt)
            }
          }
)
#' @rdname soilgrids
#' @export
setMethod('AggregateByDimension',
          signature(object = 'soilgrids'),
          function(object, dim, FUN, ...) {
            stopifnot(dim == 'time', missing('FUN'))
            stop('broken')
            object <- callNextMethod()
          }
)
#' @rdname soilgrids
#' @export
DownloadSoilGrids <- function(object) {

  stopifnot(inherits(object, 'soilgrids'), validObject(object))
  vs <- object@variables

  # Magic numbers for Soil Grids data format
  n_layers <- 7
  nval_norm <- 6
  nval_OCSTHA <- 7
  nval_BDTICM <- 1

  burl <- object@base_url
  fname <- paste0(object@nc_dir, '/', object@file_name)
  c0 <- CoordVecsToList(object@coords)
  nn0 <- object@n_fill
  nn <- nn0
  td <- tempdir()

  if (!(file.exists(fname))) {
    SetupDataFile(object)
  }
  soils_nc <- ncdf4::nc_open(fname, write = T, readunlim = T)

  cat('\n')
  # work loop:
  for (i in c(nn0:nrow(c0))) {
    i0 <- c0[i, 1]
    i1 <- c0[i, 2]
    bb <- burl
    bb[which(is.na(bb))[1]] <- i0
    bb[which(is.na(bb))[1]] <- i1
    bb <- paste0(bb, collapse = '')
    f0 <- paste0()
    fname <- paste0(td, '/url_', i, '_tmp.json')

    rb <- 1L
    repeat {
      rb <- rb + 1L
      try(download.file(url = bb, destfile = fname, quiet = T), silent = T)
      if (file.exists(fname)) {
        break
      } else {
        Sys.sleep(60)
      }
      if (rb > 5) {
        Sys.sleep(120)
      }
      if (rb > 10) {
        message('downloading failed, returning updated object')
        object@n_fill <- nn
        ncdf4::nc_close(soils_nc)
        return(object)
      }
    }

    ijs <- jsonlite::read_json(fname)
    idf <- data.frame(matrix(ncol = length(vs), nrow = n_layers), stringsAsFactors = F)
    colnames(idf) <- vs

    for (j in seq_along(vs)) {
      jj <- vs[j]
      vals <- as.numeric(unlist(ijs$properties[[jj]])[-1])
      if (!(length(vals) %in% c(nval_norm, nval_OCSTHA, nval_BDTICM))) {
        next
      }
      if (jj == 'OCSTHA') {
        idf[c(2:n_layers), jj] <- vals
      } else if (jj == 'BDTICM') {
        idf[1, jj] <- vals
      } else {
        idf[, jj] <- vals
      }
    } # end j

    time <- object@time
    if (length(time) != nrow(idf)) stop('bad row count?')
    sample <- rep(object@sample[i], nrow(idf))
    idf <- data.frame(time, sample, idf, stringsAsFactors = F)

    cat('\r', round(i / nrow(c0) * 100, 2), '%  sample:', i)
    FillArray(object, df = idf, nc = soils_nc)

    nn <- nn + 1L
    object@n_fill <- nn
  }

  ncdf4::nc_close(soils_nc)
  return(object)
}
#' @rdname soilgrids
#' @export
SoilGridMeta <- function() {
  vars <- c(
    'ALUM3S', 'AWCh1', 'AWCh2', 'AWCh3', 'AWCtS',
    'BDRICM', 'BDRLOG', 'BDTICM', 'BLDFIE',
    'CECSOL', 'CLYPPT', 'CRFVOL',
    'DRAINFAO',
    'EALKCL', 'ECAX', 'EMGX', 'ENAX', 'EXKX',

    'GLC100m',

    'NTO',
    'OCSTHA', 'ORCDRC',
    'PHIHOX', 'PHIKCL', 'PREMRG',

    'SLTPPT', 'SNDPPT',
    'TAXNWRB', 'TAXOUSDA', 'TEXMHT', 'TMDMOD_2001', 'TMDMOD_2011', 'TMNMOD_2011', 'TMNMOD_2011',

    'WWP')
  units <- c(
    '', '', '', '', '',
    '', '', '', '',
    '', '', '',
    '',
    '', '', '', '', '',

    '',

    '',
    'tonnes ha-1', '',
    '', '', '',

    '', '',
    '', '', '', '', '', '', '',

    ''
  )
  # mv is the 'major' set
  mv <- c('AWCtS', 'BDTICM', 'BLDFIE', 'CECSOL', 'CLYPPT',
          'CRFVOL', 'OCSTHA', 'PHIHOX', 'SLTPPT', 'SNDPPT', 'TEXMHT')
  major_var <- rep(FALSE, length(vars))
  major_var <- vars %in% mv
  out_df <- data.frame(vars, major_var, stringsAsFactors = F)
  return(out_df)
}
#' @rdname soilgrids
#' @export
SumSoilCarbonStocks <- function(x, d1, d2) {
  dmin <- c(0, 5, 15, 30, 60, 100)
  dmax <- c(5, 15, 30, 60, 100, 200)
  if (!(all(d1 %in% dmin, d2 %in% dmax))) {
    stop('Error in summing carbon stocks, bad depth argument')
  }
  b1 <- which(dmin == d1)
  b2 <- which(dmax == d2)
  z <- sum(x[b1:b2])
  return(z)
}
#' @rdname soilgrids
#' @export
TrapezoidalIntegrateSoil <- function(ref, d1, d2) {
  # d1/d2 are the min/max depths of desired range, respectively, while
  # ref is the reference values at 0, 5, 15, 30, 60, 100, 200 cm respectively
  d <- c(0, 5, 15, 30, 60, 100, 200)
  nl <- which(d %in% c(d1:d2))
  z <- 0
  for (k in 1:(length(nl) - 1)) {
    dd <- d[nl[k + 1]] - d[nl[k]]
    val <- x[nl[k]] + x[nl[k + 1]]
    z <- z + (dd * val)
  }
  l <- 0.5 * (1 / (d2 - d1))
  return(round(z * l, 2))
}