From ca355f5d75dd52a8afe063379ec0b2f6eb1ff46e Mon Sep 17 00:00:00 2001 From: Andrew Gene Brown Date: Mon, 1 Apr 2024 07:47:42 -0700 Subject: [PATCH] SDA_spatialQuery: refactor .SDA_geometrySelector for #342 and #343 --- DESCRIPTION | 2 +- R/SDA-spatial.R | 95 +++++++++++------------------------------ man/SDA_spatialQuery.Rd | 5 ++- 3 files changed, 31 insertions(+), 71 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index cd99f5d0..fd709115 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -17,7 +17,7 @@ Suggests: jsonlite, xml2, httr, rvest, odbc, RSQLite, sf, wk, terra, raster, kni Repository: CRAN URL: https://ncss-tech.github.io/soilDB/, https://ncss-tech.github.io/AQP/ BugReports: https://github.com/ncss-tech/soilDB/issues -RoxygenNote: 7.3.0 +RoxygenNote: 7.3.1 Encoding: UTF-8 Language: en-US Roxygen: list(markdown = TRUE) diff --git a/R/SDA-spatial.R b/R/SDA-spatial.R index 1c57c403..2022df12 100644 --- a/R/SDA-spatial.R +++ b/R/SDA-spatial.R @@ -44,79 +44,36 @@ processSDA_WKT <- function(d, g='geom', crs = 4326, p4s = NULL, as_sf = TRUE) { return(sf::as_Spatial(sfobj)) } -## TODO consider adding an 'identity' method for more generic use -## -- what is the intent here? AB 2024/04/01 - # STATSGO features require AND CLIPAREASYMBOL = 'US' to avoid state areasymbol copy # select the right query for SSURGO / STATSGO geometry filters submitted to SDA .SDA_geometrySelector <- function(db, method, geomAcres = TRUE) { + db_table <- switch(db, + SSURGO = "mupolygon", + STATSGO = "gsmmupolygon", + SAPOLYGON = "sapolygon") + db_column <- switch(db, + SSURGO = "mupolygongeo", + STATSGO = "mupolygongeo", + SAPOLYGON = "sapolygongeo") + id_column <- switch(db, + SSURGO = "mukey", + STATSGO = "mukey", + SAPOLYGON = "areasymbol") area_ac_sql <- ", GEOGRAPHY::STGeomFromWKB(geom.STUnion(geom.STStartPoint()).STAsBinary(), 4326).MakeValid().STArea() * 0.000247105 AS area_ac" - # db_table <- switch(db, - # SSURGO = "mupolygon", - # STATSGO = "gsmmupolygon", - # SAPOLYGON = "sapolygon") - # intersection_sql <- sprintf(" - # WITH geom_data (geom, mukey) AS ( - # SELECT mupolygongeo.STIntersection(geometry::STGeomFromText('%%s', 4326)) AS geom, P.mukey - # FROM %s AS P - # WHERE mupolygongeo.STIntersects(geometry::STGeomFromText('%%s', 4326)) = 1 - # ) SELECT geom.STAsText() AS geom, mukey %s - # FROM geom_data", - # db_table, ifelse(geomAcres, area_ac_sql, "")) - # overlap_sql <- sprintf(" - # SELECT mupolygongeo.STAsText() AS geom, P.mukey %s - # FROM %s AS P - # WHERE mupolygongeo.STIntersects(geometry::STGeomFromText('%%s', 4326)) = 1", - # db_table, ifelse(geomAcres, area_ac_sql, "")) - res <- switch(db, - SSURGO = switch(method, - intersection = sprintf(" - WITH geom_data (geom, mukey) AS ( - SELECT mupolygongeo.STIntersection(geometry::STGeomFromText('%%s', 4326)) AS geom, P.mukey - FROM mupolygon AS P - WHERE mupolygongeo.STIntersects(geometry::STGeomFromText('%%s', 4326)) = 1 - ) SELECT geom.STAsText() AS geom, mukey %s - FROM geom_data", - ifelse(geomAcres, area_ac_sql, "")), - overlap = sprintf(" - SELECT mupolygongeo.STAsText() AS geom, P.mukey %s - FROM mupolygon AS P - WHERE mupolygongeo.STIntersects(geometry::STGeomFromText('%%s', 4326)) = 1", - ifelse(geomAcres, area_ac_sql, ""))), - STATSGO = switch(method, - intersection = sprintf(" - WITH geom_data (geom, mukey) AS ( - SELECT mupolygongeo.STIntersection(geometry::STGeomFromText('%%s', 4326)) AS geom, P.mukey - FROM gsmmupolygon AS P - WHERE mupolygongeo.STIntersects(geometry::STGeomFromText('%%s', 4326)) = 1 - AND CLIPAREASYMBOL = 'US' - ) - SELECT geom.STAsText() AS geom, mukey %s - FROM geom_data;", - ifelse(geomAcres, area_ac_sql, "")), - overlap = sprintf(" - SELECT mupolygongeo.STAsText() AS geom, P.mukey %s - FROM gsmmupolygon AS P - WHERE mupolygongeo.STIntersects(geometry::STGeomFromText('%%s', 4326)) = 1 - AND CLIPAREASYMBOL = 'US';", - ifelse(geomAcres, area_ac_sql, ""))), - SAPOLYGON = switch(method, - intersection = sprintf(" - WITH geom_data (geom, areasymbol) AS ( - SELECT - sapolygongeo.STIntersection(geometry::STGeomFromText('%%s', 4326)) AS geom, areasymbol - FROM sapolygon - WHERE sapolygongeo.STIntersects(geometry::STGeomFromText('%%s', 4326)) = 1 - ) - SELECT geom.STAsText() AS geom, areasymbol %s - FROM geom_data;", - ifelse(geomAcres, area_ac_sql, "")), - overlap = sprintf(" - SELECT sapolygongeo.STAsText() AS geom, areasymbol %s - FROM sapolygon - WHERE sapolygongeo.STIntersects(geometry::STGeomFromText('%%s', 4326) ) = 1", - ifelse(geomAcres, area_ac_sql, ""))) - ) + geom_sql <- sprintf(switch(method, + intersection = "%s.STIntersection(geometry::STGeomFromText('%%s', 4326)) AS geom", + overlap = "%s AS geom"), db_column) + res <- sprintf( + "WITH geom_data (geom, %s) AS ( + SELECT %s, %s + FROM %s + WHERE %s.STIntersects(geometry::STGeomFromText('%%s', 4326)) = 1 + AND CLIPAREASYMBOL = 'US' + ) SELECT geom.STAsText() AS geom, %s%s + FROM geom_data", + id_column, geom_sql, id_column, db_table, db_column, id_column, + ifelse(geomAcres, area_ac_sql, "") + ) return(res) } diff --git a/man/SDA_spatialQuery.Rd b/man/SDA_spatialQuery.Rd index b7891d0f..b51a809f 100644 --- a/man/SDA_spatialQuery.Rd +++ b/man/SDA_spatialQuery.Rd @@ -8,6 +8,7 @@ SDA_spatialQuery( geom, what = "mukey", geomIntersection = FALSE, + geomAcres = TRUE, db = c("SSURGO", "STATSGO", "SAPOLYGON"), byFeature = FALSE, idcol = "gid", @@ -20,7 +21,9 @@ SDA_spatialQuery( \item{what}{a character vector specifying what to return. \code{'mukey'}: \code{data.frame} with intersecting map unit keys and names, \code{'mupolygon'} overlapping or intersecting map unit polygons from selected database, \code{'areasymbol'}: \code{data.frame} with intersecting soil survey areas, \code{'sapolygon'}: overlapping or intersecting soil survey area polygons (SSURGO only)} -\item{geomIntersection}{logical; \code{FALSE}: overlapping map unit polygons returned, \code{TRUE}: intersection of \code{geom} + map unit polygons is returned.} +\item{geomIntersection}{logical; \code{FALSE} (default): overlapping map unit polygons returned, \code{TRUE}: intersection of \code{geom} + map unit polygons is returned.} + +\item{geomAcres}{logical; \code{TRUE} (default): calculate acres of result geometry in column \code{"area_ac"} when \code{what} returns a geometry column. \code{FALSE} does not calculate acres.} \item{db}{a character vector identifying the Soil Geographic Databases (\code{'SSURGO'} or \code{'STATSGO'}) to query. Option \var{STATSGO} works with \code{what = "mukey"} and \code{what = "mupolygon"}.}