R/get_SDA_hydric.R

In ncss-tech/soilDB: Soil Database Interface

# Based on ssurgoOnDemand by chad ferguson and jason nemecek
# SDA_hydric.R: translation of SDA_hydric.py into soilDB-style R function by andrew brown
# last update: 2021/04/03

#' Get map unit hydric soils information from Soil Data Access
#'
#' Assess the hydric soils composition of a map unit.
#'
#' The default classes for `method="MAPUNIT"` are as follows:
#' - `'Nonhydric'` - no hydric components
#' - `'Hydric'` - all hydric components
#' - `'Predominantly Hydric'` - hydric component percentage is 50% or more
#' - `'Partially Hydric'` - one or more of the major components is hydric
#' - `'Predominantly Nonhydric'` - hydric component percentage is less than 50%
#'
#' The default result will also include the following summaries of component percentages: `total_comppct`, `hydric_majors` and `hydric_inclusions`.
#'
#' @details Default `method` `"Mapunit"` produces aggregate summaries of all components in the mapunit. Use `"Dominant Component"` and `"Dominant Condition"` to get the dominant component (highest percentage) or dominant hydric condition (similar conditions aggregated across components), respectively. Use `"None"` for no aggregation (one record per component).
#'
#' @param areasymbols vector of soil survey area symbols
#' @param mukeys vector of map unit keys
#' @param WHERE character containing SQL WHERE clause specified in terms of fields in `legend`, `mapunit`, or `component` tables, used in lieu of `mukeys` or `areasymbols`
#' @param method One of: `"Mapunit"`, `"Dominant Component"`, `"Dominant Condition"`, `"None"`
#' @param query_string Default: `FALSE`; if `TRUE` return a character string containing query that would be sent to SDA via `SDA_query`
#' @param dsn Path to local SQLite database or a DBIConnection object. If `NULL` (default) use Soil Data Access API via `SDA_query()`.
#' @author Jason Nemecek, Chad Ferguson, Andrew Brown
#' @return a data.frame
#' @export
get_SDA_hydric <- function(areasymbols = NULL, mukeys = NULL, WHERE = NULL, method = "MAPUNIT", query_string = FALSE, dsn = NULL) {

        method <- match.arg(toupper(method), c("MAPUNIT", "DOMINANT COMPONENT", "DOMINANT CONDITION", "NONE"))

        if (is.null(mukeys) && is.null(areasymbols) && is.null(WHERE)) {
          stop("Please specify one of the following arguments: mukeys, areasymbols, WHERE", call. = FALSE)
        }

        if (!is.null(mukeys)) {
          WHERE <- paste("mapunit.mukey IN", format_SQL_in_statement(as.integer(mukeys)))
        } else if (!is.null(areasymbols)) {
          WHERE <- paste("legend.areasymbol IN", format_SQL_in_statement(areasymbols))
        }
        .h0 <- function(w) .LIMIT_N(paste("SELECT ISNULL(SUM(comppct_r), 0) FROM mapunit AS mu INNER JOIN component AS c ON c.mukey = mu.mukey AND mapunit.mukey = mu.mukey", w), n = 1, sqlite = !is.null(dsn))

        q <- paste0("SELECT mapunit.mukey, areasymbol, musym, mapunit.muname, 
                      (", .h0(""), ") AS total_comppct,
                      (", .h0("AND majcompflag = 'Yes'"), ") AS count_maj_comp,
                      (", .h0("AND hydricrating = 'Yes'"), ") AS all_hydric,
                      (", .h0("AND majcompflag = 'Yes' AND hydricrating = 'Yes'"), ") AS hydric_majors,
                      (", .h0("AND majcompflag = 'Yes' AND hydricrating != 'Yes'"), ") AS maj_not_hydric,
                      (", .h0("AND majcompflag != 'Yes' AND hydricrating = 'Yes'"), ") AS hydric_inclusions,
                      (", .h0("AND hydricrating != 'Yes'"), ") AS all_not_hydric,
                      (", .h0("AND hydricrating IS NULL"), ") AS hydric_null
                    INTO #main_query
                    FROM legend
                    INNER JOIN mapunit ON mapunit.lkey = legend.lkey AND ", WHERE, "
                    SELECT mukey, areasymbol, musym, muname,
                           total_comppct AS total_comppct,
                           hydric_majors AS hydric_majors,
                           hydric_inclusions AS hydric_inclusions,
                    CASE WHEN total_comppct = all_not_hydric + hydric_null THEN 'Nonhydric'
                    WHEN total_comppct = all_hydric THEN 'Hydric'
                    WHEN hydric_majors + hydric_inclusions >= total_comppct / 2 THEN 'Predominantly Hydric'
                    WHEN hydric_majors > 0 THEN 'Partially Hydric'
                    WHEN hydric_majors + hydric_inclusions < total_comppct / 2 THEN 'Predominantly Nonhydric'
                    ELSE 'Error' END AS HYDRIC_RATING
                    FROM #main_query")
                    # TODO: refactor out the temp table and CASE WHEN for HYDRIC_RATING calculated in R

        comp_selection <- ""
        hyd_selection <- ""
        if (method != "MAPUNIT") {
           if (method %in% c("DOMINANT COMPONENT", "DOMINANT CONDITION")) {
                   comp_selection <- sprintf("AND component.cokey = (%s)", .LIMIT_N("SELECT c1.cokey FROM component AS c1
                 INNER JOIN mapunit AS mu1 ON c1.mukey = mu1.mukey AND c1.mukey = mapunit.mukey
                 ORDER BY c1.comppct_r DESC, c1.cokey", n = 1, sqlite = !is.null(dsn)))
           }

           if (method == "DOMINANT CONDITION") {
                   hyd_selection <- sprintf("AND hydricrating = (%s)", .LIMIT_N("SELECT hydricrating FROM mapunit AS mu
                INNER JOIN component AS c ON c.mukey = mu.mukey AND mapunit.mukey = mu.mukey
                GROUP BY hydricrating, comppct_r
                ORDER BY SUM(comppct_r) OVER (PARTITION BY hydricrating) DESC", n = 1, sqlite = !is.null(dsn)))
           }

           q <- sprintf(paste0("SELECT areasymbol, musym, muname, mapunit.mukey, ",
                   ifelse(method == "DOMINANT CONDITION", "", "cokey, compname, compkind, comppct_r, majcompflag, "),
                   "hydricrating
                 FROM legend
                 INNER JOIN mapunit ON mapunit.lkey = legend.lkey
                 INNER JOIN component ON component.mukey = mapunit.mukey %s %s
                 WHERE %s"), comp_selection, hyd_selection, WHERE)
   }

   if (!is.null(dsn)) {
     q <- gsub("ISNULL", "IFNULL", q)
   }

   if (query_string) {
    return(q)
   }

   # execute query
   if (is.null(dsn)) {
     res <- suppressMessages(SDA_query(q))
   } else {
     if (!inherits(dsn, 'DBIConnection')) {
       dsn <- RSQLite::dbConnect(RSQLite::SQLite(), dsn)
       on.exit(DBI::dbDisconnect(dsn), add = TRUE)
     }
     res <- dbGetQuery(dsn, q)
   }

   # stop if bad
   if (inherits(res, 'try-error')) {
    warnings()
    stop(attr(res, 'condition'))
   }

   return(res)
}