Nothing
# new tools for querying lab data from SDA Lab Data Mart
#
#' Query data from Kellogg Soil Survey Laboratory Data Mart via Soil Data Access or local SQLite snapshot
#'
#' LDM model diagram: \url{https://jneme910.github.io/Lab_Data_Mart_Documentation/Documents/SDA_KSSL_Data_model.html}
#'
#' @param x A vector of values to find in column specified by `what`, default `NULL` uses no constraints on `what`
#' @param what A single column name from tables: `lab_combine_nasis_ncss`, `lab_webmap`, `lab_site`, `lab_pedon` or `lab_area`
#' @param bycol A single column name from `lab_layer` used for processing chunks; default: `"pedon_key"`
#' @param tables A vector of table names; Default is `"lab_physical_properties"`, `"lab_chemical_properties"`, `"lab_calculations_including_estimates_and_default_values"`, and `"lab_rosetta_Key"`. May also include one or more of: `"lab_mir"`, `"lab_mineralogy_glass_count"`, `"lab_major_and_trace_elements_and_oxides"`, `"lab_xray_and_thermal"` but it will be necessary to select appropriate `prep_code` and `analyzed_size_frac` for your analysis (see _Details_).
#' @param WHERE character. A custom SQL WHERE clause, which overrides `x`, `what`, and `bycol`, such as `CASE WHEN corr_name IS NOT NULL THEN LOWER(corr_name) ELSE LOWER(samp_name) END = 'musick'`
#' @param chunk.size Number of pedons per chunk (for queries that may exceed `maxJsonLength`)
#' @param ntries Number of tries (times to halve `chunk.size`) before returning `NULL`; default `3`
#' @param layer_type Default: `"horizon"`, `"layer"`, and `"reporting layer"`
#' @param prep_code Default: `"S"` and `""`. May also include one or more of: `"F"`, `"HM"`, `"HM_SK"` `"GP"`, `"M"`, `"N"`, or `"S"`
#' @param analyzed_size_frac Default: `"<2 mm"` and `""`. May also include one or more of: `"<0.002 mm"`, `"0.02-0.05 mm"`, `"0.05-0.1 mm"`, `"0.1-0.25 mm"`, `"0.25-0.5 mm"`, `"0.5-1 mm"`, `"1-2 mm"`, `"0.02-2 mm"`, `"0.05-2 mm"`
#' @param dsn Data source name; either a path to a SQLite database, an open DBIConnection or (default) `NULL` (to use `soilDB::SDA_query`)
#'
#' @details If the `chunk.size` parameter is set too large and the Soil Data Access request fails, the algorithm will re-try the query with a smaller (halved) `chunk.size` argument. This will be attempted up to 3 times before returning `NULL`
#'
#' Currently the `lab_area` tables are joined only for the "Soil Survey Area" related records.
#'
#' When requesting data from `"lab_major_and_trace_elements_and_oxides"`, `"lab_mineralogy_glass_count"`, or `"lab_xray_and_thermal"` multiple preparation codes (`prep_code`) or size fractions (`analyzed_size_frac`) are possible. The default behavior of `fetchLDM()` is to attempt to return a topologically valid (minimal overlaps) _SoilProfileCollection_. This is achieved by setting `prep_code="S"` ("sieved") and `analyzed_size_frac="<2 mm"`. You may specify alternate or additional preparation codes or fractions as needed, but note that this may cause "duplication" of some layers where measurements were made with different preparation or on fractionated samples
#'
#' @return a `SoilProfileCollection` for a successful query, a `try-error` if no site/pedon locations can be found or `NULL` for an empty `lab_layer` (within sites/pedons) result
#' @export
#' @examplesIf curl::has_internet()
#' @examples
#' \dontrun{
#' # fetch by ssa_key
#' res <- fetchLDM(8297, what = "ssa_key")
#'
#' # physical properties correlated as taxonomic subgroup "Typic Argialbolls"
#' res <- fetchLDM(x = "Typic Argialbolls",
#' what = "corr_taxsubgrp",
#' tables = "lab_physical_properties")
#'
#' # fetch by area_code (SSA only)
#' res <- fetchLDM("CA630", what = "area_code")
#' }
#' @importFrom aqp `depths<-` `site<-`
#' @importFrom data.table rbindlist
fetchLDM <- function(x = NULL,
what = "pedlabsampnum",
bycol = "pedon_key",
tables = c(
"lab_physical_properties",
"lab_chemical_properties",
"lab_calculations_including_estimates_and_default_values",
# # remove alternate prep/fractionated data from default set of tables
# "lab_major_and_trace_elements_and_oxides",
# "lab_mineralogy_glass_count",
# "lab_xray_and_thermal",
"lab_rosetta_Key"
# "lab_mir"
),
WHERE = NULL,
chunk.size = 1000,
ntries = 3,
layer_type = c("horizon","layer","reporting layer"),
prep_code = c("S", ""), # , `"F"`, `"HM"`, `"HM_SK"` `"GP"`, `"M"`, `"N"`, or `"S"`
analyzed_size_frac = c("<2 mm", ""),# optional: "<0.002 mm", "0.02-0.05 mm", "0.05-0.1 mm", "0.1-0.25 mm", "0.25-0.5 mm", "0.5-1 mm", "1-2 mm", "0.02-2 mm", "0.05-2 mm"
dsn = NULL) {
# set up data source connection if needed
if (inherits(dsn, 'DBIConnection')) {
# allow any existing DBI connection to be passed via dsn argument
con <- dsn
} else if (!is.null(dsn) && is.character(dsn)) {
# if it is a path to data source, try to connect with RSQLite
con <- RSQLite::dbConnect(RSQLite::SQLite(), dsn)
} else {
# otherwise we are using SDA_query
con <- NULL
}
lab_combine_nasis_ncss <- c("pedon_key", "site_key", "pedlabsampnum", "pedoniid", "upedonid",
"labdatadescflag", "priority", "priority2", "samp_name", "samp_class_type",
"samp_classdate", "samp_classification_name", "samp_taxorder",
"samp_taxsuborder", "samp_taxgrtgroup", "samp_taxsubgrp", "samp_taxpartsize",
"samp_taxpartsizemod", "samp_taxceactcl", "samp_taxreaction",
"samp_taxtempcl", "samp_taxmoistscl", "samp_taxtempregime", "samp_taxminalogy",
"samp_taxother", "samp_osdtypelocflag", "corr_name", "corr_class_type",
"corr_classdate", "corr_classification_name", "corr_taxorder",
"corr_taxsuborder", "corr_taxgrtgroup", "corr_taxsubgrp", "corr_taxpartsize",
"corr_taxpartsizemod", "corr_taxceactcl", "corr_taxreaction",
"corr_taxtempcl", "corr_taxmoistscl", "corr_taxtempregime", "corr_taxminalogy",
"corr_taxother", "corr_osdtypelocflag", "SSL_name", "SSL_class_type",
"SSL_classdate", "SSL_classification_name", "SSL_taxorder", "SSL_taxsuborder",
"SSL_taxgrtgroup", "SSL_taxsubgrp", "SSL_taxpartsize", "SSL_taxpartsizemod",
"SSL_taxceactcl", "SSL_taxreaction", "SSL_taxtempcl", "SSL_taxmoistscl",
"SSL_taxtempregime", "SSL_taxminalogy", "SSL_taxother", "SSL_osdtypelocflag",
"siteiid", "usiteid", "site_obsdate", "latitude_decimal_degrees",
"longitude_decimal_degrees", "country_key", "state_key", "county_key",
"mlra_key", "ssa_key", "npark_key", "nforest_key", "note", "samp_taxfamhahatmatcl",
"corr_taxfamhahatmatcl", "SSL_taxfamhahatmatcl", "pedobjupdate",
"siteobjupdate")
lab_area <- c("area_key", "area_type", "area_sub_type", "parent_area_key",
"parent_org_key", "area_code", "area_name", "area_abbrev", "area_desc")
lab_webmap <- c("wmiid", "Series", "User_pedon_ID", "pedon_Key", "peiid", "Soil_Classification")
lab_site <- c("site_key", "user_site_id")
lab_pedon <- c("pedon_key", "pedlabsampnum", "observation_date")
if (!is.null(x)) {
if (what %in% lab_combine_nasis_ncss)
what <- paste0("lab_combine_nasis_ncss.", what)
if (what %in% lab_area)
what <- paste0("lab_area.", what)
if (what %in% lab_webmap)
what <- paste0("lab_webmap.", what)
if (what %in% lab_site)
what <- paste0("lab_site.", what)
if (what %in% lab_pedon)
what <- paste0("lab_pedon.", what)
what <- match.arg(what, choices = c(paste0("lab_combine_nasis_ncss.", lab_combine_nasis_ncss),
paste0("lab_area.", lab_area),
paste0("lab_webmap.", lab_webmap),
paste0("lab_site.", lab_site),
paste0("lab_pedon.", lab_pedon)))
}
# TODO: set up arbitrary area queries by putting area table into groups:
# country, state, county, mlra, ssa, npark, nforest
if (!is.null(x) && (missing(WHERE) || is.null(WHERE))) {
WHERE <- sprintf("LOWER(%s) IN %s", what, format_SQL_in_statement(tolower(x)))
}
# get site/pedon/area information
site_query <- paste0(
"SELECT * FROM lab_combine_nasis_ncss
LEFT JOIN lab_webmap ON
lab_combine_nasis_ncss.pedon_key = lab_webmap.pedon_key
LEFT JOIN lab_site ON
lab_combine_nasis_ncss.site_key = lab_site.site_key
LEFT JOIN lab_pedon ON
lab_combine_nasis_ncss.site_key = lab_pedon.site_key ",
ifelse(is.null(x) && is.null(WHERE), "", paste("WHERE", WHERE)))
if (inherits(con, 'DBIConnection')) {
# query con using (modified) site_query
# if (FALSE) {
# # fixes for old version of SQLite db
# site_query <- gsub("\\blab_|\\blab_combine_", "", gsub("lab_(site|pedon)", "nasis_\\1", site_query))
# }
sites <- try(DBI::dbGetQuery(con, site_query))
} else {
# the lab_area table allows for overlap with many different area types
# for now we only offer the "ssa" (soil survey area) area_type
site_query_ssaarea <- gsub("WHERE",
"LEFT JOIN lab_area ON
lab_combine_nasis_ncss.ssa_key = lab_area.area_key
WHERE", site_query)
sites <- suppressMessages(SDA_query(site_query_ssaarea))
}
if (!inherits(sites, 'try-error') && !is.null(sites)) {
# TODO: this shouldn't be needed
sites <- sites[,unique(colnames(sites))]
if (is.null(chunk.size) || nrow(sites) < chunk.size) {
# get data for lab layers within pedon_key returned
hz <- .get_lab_layer_by_pedon_key(x = sites[[bycol]],
con = con,
bycol = bycol,
tables = tables,
layer_type = layer_type,
prep_code = prep_code,
analyzed_size_frac = analyzed_size_frac)
if (inherits(hz, 'try-error') && grepl("Layer data query failed", hz[1])) {
return(hz)
}
} else {
hz <- try(stop("Failed to access Web Service or query returned no results", call. = FALSE), silent = TRUE)
}
.do_chunk <- function(con, size) {
chunk.idx <- makeChunks(sites[[bycol]], size)
as.data.frame(data.table::rbindlist(lapply(unique(chunk.idx),
function(i) {
keys <- sites[[bycol]][chunk.idx == i]
res <- .get_lab_layer_by_pedon_key(x = keys,
con = con,
bycol = bycol,
tables = tables,
prep_code = prep_code,
analyzed_size_frac = analyzed_size_frac)
if (inherits(res, 'try-error'))
return(NULL)
res
})))
}
ntry <- 0
while ((inherits(hz, 'try-error') || is.null(hz)) && ntry < ntries) {
if (is.null(chunk.size)) {
stop("query failed and chunk.size argument is NULL", call. = FALSE)
}
hz <- .do_chunk(con, chunk.size)
# repeat as long as there is a try error/NULL, halving chunk.size with each iteration
chunk.size <- pmax(floor(chunk.size / 2), 1)
ntry <- ntry + 1
}
# close connection (if we opened it)
if (is.character(dsn) && inherits(con, "DBIConnection")) {
DBI::dbDisconnect(con)
}
if (!is.null(hz) && nrow(hz) > 0) {
# local SQLite: sometimes prep_code differs for bulk density with no difference in value
# e.g. pedon_key 10010
# prep_code is also repeated across multiple tables
# hz$prep_code <- NULL
# site_key used in multiple tables
hz$site_key <- NULL
# hacks to deal with problems in the various databases
# hz <- unique(hz[,unique(colnames(hz))])
# build SoilProfileCollection
depths(hz) <- pedon_key ~ hzn_top + hzn_bot
hzn <- aqp::horizonNames(hz)
hzn <- hzn[hzn != idname(hz)]
site(hz) <- sites[, !colnames(sites) %in% hzn]
hzdesgnname(hz) <- 'hzn_desgn'
return(hz)
} else {
message("no horizon data in result!")
return(NULL)
}
} else {
# return try-error
return(sites)
}
NULL
}
.get_lab_layer_by_pedon_key <- function(x,
con = NULL,
bycol = "pedon_key",
tables = c(
"lab_physical_properties",
"lab_mineralogy_glass_count",
"lab_chemical_properties",
"lab_major_and_trace_elements_and_oxides",
"lab_xray_and_thermal",
"lab_calculations_including_estimates_and_default_values",
"lab_rosetta_Key",
"lab_mir"
),
layer_type = c("horizon", "layer", "reporting layer"),
prep_code = c("GP", "M", "N", "S", "HM"),
analyzed_size_frac = c("0.02-0.05 mm", "0.05-0.1 mm", "1-2 mm",
"0.5-1 mm", "0.25-0.5 mm", "0.05-2 mm",
"<2 mm", "0.02-2 mm", "0.1-0.25 mm",
"<0.002 mm")) {
flattables <- c(
"lab_physical_properties",
"lab_chemical_properties",
"lab_calculations_including_estimates_and_default_values",
"lab_major_and_trace_elements_and_oxides",
"lab_rosetta_Key",
"lab_mir"
)
fractables <- c("lab_mineralogy_glass_count",
"lab_xray_and_thermal")
tables <- match.arg(tables, c(flattables, fractables), several.ok = TRUE)
fractablejoincriteria <- list(
"lab_mineralogy_glass_count" = "LEFT JOIN lab_mineralogy_glass_count ON
lab_layer.labsampnum = lab_mineralogy_glass_count.labsampnum",
"lab_xray_and_thermal" = "LEFT JOIN lab_xray_and_thermal ON
lab_layer.labsampnum = lab_xray_and_thermal.labsampnum"
)
tablejoincriteria <- list(
"lab_physical_properties" = "LEFT JOIN lab_physical_properties ON
lab_layer.labsampnum = lab_physical_properties.labsampnum",
"lab_chemical_properties" = "LEFT JOIN lab_chemical_properties ON
lab_layer.labsampnum = lab_chemical_properties.labsampnum",
"lab_major_and_trace_elements_and_oxides" = "LEFT JOIN lab_major_and_trace_elements_and_oxides ON
lab_layer.labsampnum = lab_major_and_trace_elements_and_oxides.labsampnum",
"lab_xray_and_thermal" = "LEFT JOIN lab_xray_and_thermal ON
lab_layer.labsampnum = lab_xray_and_thermal.labsampnum",
"lab_calculations_including_estimates_and_default_values" = "LEFT JOIN lab_calculations_including_estimates_and_default_values ON
lab_layer.labsampnum = lab_calculations_including_estimates_and_default_values.labsampnum",
"lab_rosetta_Key" = "LEFT JOIN lab_rosetta_Key ON
lab_layer.layer_key = lab_rosetta_Key.layer_key",
"lab_mir" = "LEFT JOIN lab_mir ON lab_layer.layer_key = lab_mir.layer_key
LEFT JOIN lab_mir_wavelength ON lab_mir.d_wavelength_array_id = lab_mir_wavelength.d_wavelength_array_id"
)
layer_type <- match.arg(layer_type, c("horizon", "layer", "reporting layer"), several.ok = TRUE)
if (any(tables %in% flattables)) {
layer_query <- sprintf(
"SELECT * FROM lab_layer %s WHERE lab_layer.layer_type IN %s %s AND %s",
paste0(sapply(flattables[flattables %in% tables], function(a) tablejoincriteria[[a]]), collapse = "\n"),
format_SQL_in_statement(layer_type),
ifelse(is.null(x), "", paste0(" AND ", bycol, " IN ", format_SQL_in_statement(x))),
paste0(paste0(sapply(flattables[flattables %in% tables[!tables %in% c("lab_rosetta_Key", "lab_mir")]],
function(b) paste0("IsNull(",b,".prep_code, '')")),
" IN ", format_SQL_in_statement(prep_code)), collapse = " AND "))
} else {
layer_query <- sprintf(
"SELECT * FROM lab_layer WHERE lab_layer.layer_type IN %s %s",
format_SQL_in_statement(layer_type),
ifelse(is.null(x), "", paste0(" AND ", bycol, " IN ", format_SQL_in_statement(x))))
}
if (any(tables %in% fractables)) {
layer_fraction_query <- sprintf(
"SELECT * FROM lab_layer %s WHERE lab_layer.layer_type IN %s AND %s IN %s AND %s AND %s",
paste0(sapply(fractables[fractables %in% tables], function(a) fractablejoincriteria[[a]]), collapse = "\n"),
format_SQL_in_statement(layer_type),
bycol,
format_SQL_in_statement(x),
paste0(paste0(sapply(fractables[fractables %in% tables],
function(b) paste0("IsNull(",b,".prep_code, '')")),
" IN ", format_SQL_in_statement(prep_code)), collapse = " AND "),
paste0(paste0(sapply(fractables[fractables %in% tables],
function(b) paste0("IsNull(",b,".analyzed_size_frac, '')")),
" IN ", format_SQL_in_statement(analyzed_size_frac)), collapse = " AND "))
} else {
layer_fraction_query <- NULL
}
if (inherits(con, 'DBIConnection')) {
layerdata <- try(DBI::dbGetQuery(con, gsub("IsNull", "IFNULL", layer_query)))
} else {
layerdata <- suppressWarnings(SDA_query(layer_query))
}
if (inherits(layerdata, 'try-error')) {
return(try(stop("Layer data query failed", call. = FALSE)))
}
# TODO: this shouldn't be needed
layerdata <- layerdata[,unique(colnames(layerdata))]
if (!is.null(layer_fraction_query)) {
layerfracdata <- suppressWarnings(SDA_query(layer_fraction_query))
layerfracdata <- layerfracdata[,unique(colnames(layerfracdata))]
if (!inherits(layerdata, 'try-error') && !inherits(layerfracdata, 'try-error')) {
if (nrow(layerfracdata) == 0)
message(sprintf('no fractionated samples found for selected prep_code (%s) and analyzed_size_frac (%s)',
paste0(prep_code, collapse = ", "), paste0(analyzed_size_frac, collapse = ", ")))
layerdata <- merge(layerdata, layerfracdata[,c("labsampnum", colnames(layerfracdata)[!colnames(layerfracdata) %in% colnames(layerdata)])], by = "labsampnum", all.x = TRUE, incomparables = NA)
}
}
layerdata$prep_code[is.na(layerdata$prep_code)] <- ""
layerdata
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.