New functions for interacting with AWDB API: Snow Courses

R/AWDB-requests-and-util.R

@@ -0,0 +1,53 @@
+##
+##
+##
+
+## 
+.AWDB_elementSets <- list(
+  'SCAN' = c(
+    "STO", "SMS", "RDC", "SAL", "PREC", "TOBS", "SRADV", 
+    "WDIRV", "WSPDV", "WSPDX", "DPTP", "PRCP", 
+    "PVPV", "RHUM", "RHUMN", "RHUMX", "SVPV", "SNWD", 
+    "NTRDV", "WDIR", "WTEQ"
+  ), 
+  'SNTL' = c(
+    "STO", "SMS", "TOBS", "RDC", "SAL", 
+    "PREC", "WTEQ", "SNWD", "PTEMP", 
+    "WDIRV", "WSPDX", "WSPDV", "SRADV", "RHUMV", "SRADX", "RHUMX", 
+    "RHUM", "SRAD", "SRADN", "RHUMN", "PRES", "SWINV", "SWOTV", "WDIRZ", 
+    "LWINV", "LWOTV", "PVPV"
+  ), 
+  'SNOW' = c("SNDN", "SNWD", "WTEQ"), 
+  'SNTLT' = c(
+    "TOBS", "SNWD", "STO", "SMS", "WSPDV", 
+    "WDIRV", "WSPDX", "PREC", "RDC", "SAL", "WDIRZ", 
+    "SMV", "SRADV", "STV", "AWDC", "RHUM", 
+    "RHUMN", "RHUMV", "RHUMX", "SRAD", "SRADN", "SRADX", "WDIR"
+  )
+)
+
+
+
+#' @title Get AWDB Element (Sensor) Metadata
+#' @description This function retrieves all element (sensor) metadata from the AWDB webservice.
+#' 
+#' @param ... Additional arguments to `.soilDB_curl_get_JSON()`, such as `timeout = 120`
+#' @return A `data.frame` unless network or server error, then `NULL`
+#' 
+get_AWDB_elements <- function(...) {
+
+  .u <- 'https://wcc.sc.egov.usda.gov/awdbRestApi/services/v1/reference-data?referenceLists=elements&visibility=all'
+  .res <- .soilDB_curl_get_JSON(.u, gzip = FALSE, quiet = TRUE, ...)
+
+  # errors result in NULL
+  if(inherits(.res, 'data.frame')) {
+    .res <- .res$elements
+  }
+
+  return(.res)
+}
+
+
+
+
+

R/SDA-spatial.R

@@ -27,7 +27,7 @@ processSDA_WKT <- function(d,
   if (!is.null(p4s)) {
     .Deprecated(msg = "Passing PROJ4 strings via `p4s` is deprecated. SDA interfaces in soilDB use the WGS84 Geographic Coordinate System (EPSG:4326) by default. Use the `crs` argument to customize.")
   }
-
+  
   if (inherits(d,  'try-error')) {
     message("Invalid SDA WKT result, returning try-error")
     return(d)
@@ -70,15 +70,15 @@ processSDA_WKT <- function(d,
                              intersection = "%s.STIntersection(geometry::STGeomFromText('%%s', 4326)) AS geom",
                              overlap = "%s AS geom"), db_column)
   res <- sprintf(
-      "WITH geom_data (geom, %s) AS (
+    "WITH geom_data (geom, %s) AS (
           SELECT %s, %s
           FROM %s 
           WHERE %s.STIntersects(geometry::STGeomFromText('%%s', 4326)) = 1 %s
         ) SELECT geom_data.%s, geom.STAsText() AS geom%s
         FROM geom_data",
-      id_column, geom_sql, id_column, db_table, db_column, clip_sql, id_column,
-      ifelse(geomAcres, area_ac_sql, "")
-    )
+    id_column, geom_sql, id_column, db_table, db_column, clip_sql, id_column,
+    ifelse(geomAcres, area_ac_sql, "")
+  )
 
   # handle non-polygon results
   if (db == "SSURGO" && what %in% c("mupoint", "muline", "featpoint", "featline")) {
@@ -323,7 +323,7 @@ SDA_spatialQuery <- function(geom,
     }))
     return(res)
   }
-
+  
   res <- .SDA_spatialQuery(
     geom = geom,
     what = what,
@@ -352,23 +352,23 @@ SDA_spatialQuery <- function(geom,
   # check for required packages
   if (!requireNamespace('sf', quietly = TRUE))
     stop('please install the `sf` package', call. = FALSE)
-
+  
   if (!requireNamespace('wk', quietly = TRUE))
     stop('please install the `wk` package', call. = FALSE)
-
+  
   what <- tolower(what)
   db <- toupper(db)
-
+  
   return_sf <- FALSE
   return_terra <- FALSE
-
+  
   # raster support
   if (inherits(geom, c('RasterLayer', 'RasterBrick', 'RasterStack'))) {
     if (!requireNamespace('terra'))
       stop("packages terra is required", call. = FALSE)
     geom <- terra::rast(geom)
   }
-
+  
   if (inherits(geom, 'SpatRaster') || inherits(geom, 'SpatVector')) {
     # terra support
     return_terra <- TRUE
@@ -386,70 +386,70 @@ SDA_spatialQuery <- function(geom,
   } else {
     stop('`geom` must be an sf, terra, or Spatial* object', call. = FALSE)
   }
-
+  
   # backwards compatibility with old value of what argument
   if (what == 'geom') {
     message("converting what='geom' to what='mupolygon'")
     what <- "mupolygon"
   }
-
+  
   # determine if requested data type is allowed
   if (!what %in% c('mukey', 'mupolygon', 'mupoint', 'muline', 'featpoint', 'featline', 'areasymbol', 'sapolygon')) {
     stop("query type (argument `what`) must be either 'mukey' / 'areasymbol' (tabular result) OR 'mupolygon', 'mupoint', 'muline', 'featpoint', 'featline', 'sapolygon' (geometry result)", call. = FALSE)
   }
-
+  
   # areasymbol is allowed with db = "SSURGO" (default) and db = "SAPOLYGON"
   if (what %in% c('areasymbol', 'sapolygon')) {
     db <- 'SAPOLYGON' # geometry selector uses db argument to specify sapolygon queries
   }
-
+  
   db <- match.arg(db)
-
+  
   if (what == 'areasymbol' && db == 'STATSGO') {
     stop("query type 'areasymbol' for 'STATSGO' is not supported", call. = FALSE)
   }
-
+  
   # geom must have a valid CRS
   if (is.na(sf::st_crs(geom)$wkt)) {
     stop('`geom` must have a valid CRS', call. = FALSE)
   }
-
+  
   # CRS conversion if needed
   target.prj <- sf::st_crs(4326)
   if (suppressWarnings(sf::st_crs(geom)) != target.prj) {
     geom <- sf::st_transform(geom, target.prj)
   }
-
+  
   # WKT encoding
   # use a geometry collection
   wkt <- wk::wk_collection(wk::as_wkt(geom))
-
+  
   # returning geom + mukey or geom + areasymbol
   if (what %in% c("mupolygon", "sapolygon", "mupoint", "muline", "featpoint", "featline")) {
-
+    
     if (what %in% c("mupoint", "muline", "featpoint", "featline")) {
       geomAcres <- FALSE
     }
 
     # return intersection + area
     if (geomIntersection) {
-
+      
       # select the appropriate query
       .template <- .SDA_geometrySelector(db = db, what = what, method = 'intersection', addFields = addFields, geomAcres = geomAcres)
       q <- sprintf(.template, as.character(wkt), as.character(wkt))
-
+      
     } else {
       # return overlapping
-
+      
       # select the appropriate query
       .template <- .SDA_geometrySelector(db = db, what = what, method = 'overlap', addFields = addFields, geomAcres = geomAcres)
       q <- sprintf(.template, as.character(wkt))
     }
-
+    
     if (query_string) {
       return(q)
     }
-
+    
     # single query for all of the features
     # note that row-order / number of rows in results may not match geom
     res <- suppressMessages(SDA_query(q))
@@ -465,7 +465,7 @@ SDA_spatialQuery <- function(geom,
     }
 
   } else {
-      
+
     if (what == 'mukey') {
       if (db == "SSURGO") {
         q <- sprintf("WITH geom_data (mukey, muname) AS (SELECT mukey, muname
@@ -484,7 +484,7 @@ SDA_spatialQuery <- function(geom,
       } else {
         stop("query type 'mukey' for 'SAPOLYGON' is not supported", call. = FALSE)
       }
-
+      
     } else if (what == 'areasymbol') {
       # SSURGO only
       q <- sprintf("WITH geom_data (areasymbol) AS (SELECT areasymbol

R/fetchSDA_spatial.R

@@ -5,15 +5,25 @@
 #' A Soil Data Access query returns geometry and key identifying information about the map unit or area of interest. Additional columns from the map unit or legend table can be included; see `add.fields` argument.
 #'
 #' @param x A vector of map unit keys (`mukey`) or national map unit symbols (`nationalmusym`) for `mupolygon`, `muline` or `mupoint`; feature keys (`featkey`) for `featpoint` and `featline`; legend keys (`lkey`) or soil survey area symbols (`areasymbol`) for `sapolygon` geometry. If `geom.src="mlrapolygon"` then `x` refers to `MLRARSYM` (major land resource area symbols).
+#' 
 #' @param by.col Column name containing map unit identifier `"mukey"`, `"nationalmusym"`, or `"ecoclassid"` for `geom.src` `mupolygon` OR `"areasymbol"`, `"areaname"`, `"mlraoffice"`, `"mouagncyresp"` for `geom.src` `sapolygon`; default is determined by `isTRUE(is.numeric(x))` for `mukey`, `featkey` or `lkey`, using `nationalmusym` or `areasymbol` otherwise.
+#' 
 #' @param method geometry result type: `"feature"` returns polygons, `"bbox"` returns the bounding box of each polygon (via `STEnvelope()`), `"point"` returns a single point (via `STPointOnSurface()`) within each polygon, `"extent"` returns an aggregate bounding box (the extent of all polygons, `geometry::EnvelopeAggregate()`) ), `"convexhull"` (`geometry::ConvexHullAggregate()`) returns the aggregate convex hull around all polygons, `"union"` (`geometry::UnionAggregate()`) and `"collection"` (`geometry::CollectionAggregate()`) return a `MULTIPOLYGON` or a `GEOMETRYCOLLECTION`, respectively, for each `mukey`, `nationalmusym`, or `areasymbol `. In the case of the latter four aggregation methods, the groups for aggregation depend on `by.col` (default by `"mukey"`).
+#' 
 #' @param geom.src Either `mupolygon` (map unit polygons), `muline` (map unit lines), `mupoint` (map unit points), `featpoint` (feature points), `featline` (feature lines), `sapolygon` (soil survey area boundary polygons), or `mlrapolygon` (major land resource area boundary polygons)
+#' 
 #' @param db Default: `"SSURGO"`. When `geom.src` is `mupolygon`, use STATSGO polygon geometry instead of SSURGO by setting `db = "STATSGO"`
+#' 
 #' @param add.fields Column names from `mapunit` or `legend` table to add to result. Must specify parent table name as the prefix before column name e.g. `mapunit.muname`.
+#' 
 #' @param chunk.size Number of values of `x` to process per query. Necessary for large results. Default: `10`
+#' 
 #' @param verbose Print messages?
+#' 
 #' @param as_Spatial Return sp classes? e.g. `Spatial*DataFrame`. Default: `FALSE`.
+#' 
 #' @return an `sf` data.frame corresponding to SDA spatial data for all symbols requested. If `as_Spatial=TRUE` returns a `Spatial*DataFrame` from the sp package via `sf::as_Spatial()` for backward compatibility. Default result contains geometry with attribute table containing unique feature ID, symbol and area symbol plus additional fields in result specified with `add.fields`.
+#' 
 #' @details
 #'
 #' This function automatically "chunks" the input vector (using `makeChunks()`) of map unit identifiers to minimize the likelihood of exceeding the SDA data request size. The number of chunks varies with the `chunk.size` setting and the length of your input vector. If you are working with many map units and/or large extents, you may need to decrease this number in order to have more chunks.

data-raw/scan-snotel-current/dump-all-sites.R

@@ -0,0 +1,12 @@
+## Optionally download all SCAN, SNOTEL, and associated cooperator site metadata
+## J.M. Skovlin and D.E. Beaudette
+## 2025-03-06
+
+
+
+# all stations and cooperators
+u <- 'https://wcc.sc.egov.usda.gov/awdbRestApi/services/v1/stations'
+
+
+
+