Skip to content

Commit

Permalink
SDA_spatialQuery: allow toggle of area_ac for #343
Browse files Browse the repository at this point in the history
  • Loading branch information
@brownag
brownag committed Apr 1, 2024
1 parent 0f5510b commit deaef38
Showing 1 changed file with 83 additions and 93 deletions.
176 changes: 83 additions & 93 deletions R/SDA-spatial.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,96 +44,80 @@ 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

## TODO consider adding an 'identity' method for more generic use
# STATSGO features require AND CLIPAREASYMBOL = 'US' to avoid state areasymbol copy
# select the right query for SSURGO / STATSGO geometry filters submitted to SDA
# this is important because STATSGO features require the additional
# AND CLIPAREASYMBOL = 'US'
.SDA_geometrySelector <- function(db, method) {

.SDA_geometrySelector <- function(db, method, geomAcres = TRUE) {
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 = {
"
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,
GEOGRAPHY::STGeomFromWKB(
geom.STUnion(geom.STStartPoint()).STAsBinary(), 4326).MakeValid().STArea() * 0.000247105 AS area_ac
FROM geom_data;
"
},

overlap = {
"
SELECT
mupolygongeo.STAsText() AS geom, P.mukey
FROM mupolygon AS P
WHERE mupolygongeo.STIntersects( geometry::STGeomFromText('%s', 4326) ) = 1;
"
}
)
},

STATSGO = {
switch(method,
intersection = {
"
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,
GEOGRAPHY::STGeomFromWKB(
geom.STUnion(geom.STStartPoint()).STAsBinary(), 4326).MakeValid().STArea() * 0.000247105 AS area_ac
FROM geom_data;
"
},

overlap = {
"
SELECT
mupolygongeo.STAsText() AS geom, P.mukey
FROM gsmmupolygon AS P
WHERE mupolygongeo.STIntersects( geometry::STGeomFromText('%s', 4326) ) = 1
AND CLIPAREASYMBOL = 'US' ;
"

}
)
},
SAPOLYGON = { switch(method,
intersection = "
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,
GEOGRAPHY::STGeomFromWKB(geom.STUnion(geom.STStartPoint()).STAsBinary(), 4326).MakeValid().STArea() * 0.000247105 AS area_ac
FROM geom_data;
",

overlap = "SELECT
sapolygongeo.STAsText() AS geom, areasymbol
FROM sapolygon
WHERE sapolygongeo.STIntersects(geometry::STGeomFromText('%s', 4326) ) = 1")}
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, "")))
)

return(res)

}

#' Query Soil Data Access by spatial intersection with supplied geometry
Expand All @@ -145,14 +129,17 @@ FROM geom_data;
#' SSURGO (detailed soil survey, typically 1:24,000 scale) and STATSGO (generalized soil survey, 1:250,000 scale) data are stored together within SDA. This means that queries that don't specify an area symbol may result in a mixture of SSURGO and STATSGO records. See the examples below and the [SDA Tutorial](http://ncss-tech.github.io/AQP/soilDB/SDA-tutorial.html) for details.
#'
#' @aliases SDA_spatialQuery
#'
#' @param geom an `sf` or `Spatial*` object, with valid CRS. May contain multiple features.
#' @param what a character vector specifying what to return. `'mukey'`: `data.frame` with intersecting map unit keys and names, `'mupolygon'` overlapping or intersecting map unit polygons from selected database, `'areasymbol'`: `data.frame` with intersecting soil survey areas, `'sapolygon'`: overlapping or intersecting soil survey area polygons (SSURGO only)
#' @param geomIntersection logical; `FALSE`: overlapping map unit polygons returned, `TRUE`: intersection of `geom` + map unit polygons is returned.
#' @param geomIntersection logical; `FALSE` (default): overlapping map unit polygons returned, `TRUE`: intersection of `geom` + map unit polygons is returned.
#' @param geomAcres logical; `TRUE` (default): calculate acres of result geometry in column `"area_ac"` when `what` returns a geometry column. `FALSE` does not calculate acres.
#' @param db a character vector identifying the Soil Geographic Databases (`'SSURGO'` or `'STATSGO'`) to query. Option \var{STATSGO} works with `what = "mukey"` and `what = "mupolygon"`.
#' @param byFeature Iterate over features, returning a combined data.frame where each feature is uniquely identified by value in `idcol`. Default `FALSE`.
#' @param idcol Unique IDs used for individual features when `byFeature = TRUE`; Default `"gid"`
#' @param query_string Default: `FALSE`; if `TRUE` return a character string containing query that would be sent to SDA via `SDA_query`
#' @param as_Spatial Return sp classes? e.g. `Spatial*DataFrame`. Default: `FALSE`.
#'
#' @return A `data.frame` if `what = 'mukey'`, otherwise an `sf` object. A `try-error` in the event the request cannot be made or if there is an error in the query.
#' @note Row-order is not preserved across features in \code{geom} and returned object. Use `byFeature` argument to iterate over features and return results that are 1:1 with the inputs. Polygon area in acres is computed server-side when `what = 'mupolygon'` and `geomIntersection = TRUE`.
#' @author D.E. Beaudette, A.G. Brown, D.R. Schlaepfer
Expand Down Expand Up @@ -294,6 +281,7 @@ FROM geom_data;
SDA_spatialQuery <- function(geom,
what = 'mukey',
geomIntersection = FALSE,
geomAcres = TRUE,
db = c("SSURGO", "STATSGO", "SAPOLYGON"),
byFeature = FALSE,
idcol = "gid",
Expand All @@ -318,6 +306,7 @@ SDA_spatialQuery <- function(geom,
geom = geom,
what = what,
geomIntersection = geomIntersection,
geomAcres = geomAcres,
db = db,
query_string = query_string
)
Expand All @@ -330,11 +319,12 @@ SDA_spatialQuery <- function(geom,
}

.SDA_spatialQuery <- function(geom,
what = 'mukey',
geomIntersection = FALSE,
db = c("SSURGO", "STATSGO", "SAPOLYGON"),
query_string = FALSE) {

what = 'mukey',
geomIntersection = FALSE,
geomAcres = TRUE,
db = c("SSURGO", "STATSGO", "SAPOLYGON"),
query_string = FALSE) {

# check for required packages
if (!requireNamespace('sf', quietly = TRUE))
stop('please install the `sf` package', call. = FALSE)
Expand Down Expand Up @@ -419,14 +409,14 @@ SDA_spatialQuery <- function(geom,
if (geomIntersection) {

# select the appropriate query
.template <- .SDA_geometrySelector(db = db, method = 'intersection')
.template <- .SDA_geometrySelector(db = db, method = 'intersection', geomAcres = geomAcres)
q <- sprintf(.template, as.character(wkt), as.character(wkt))

} else {
# return overlapping

# select the appropriate query
.template <- .SDA_geometrySelector(db = db, method = 'overlap')
.template <- .SDA_geometrySelector(db = db, method = 'overlap', geomAcres = geomAcres)
q <- sprintf(.template, as.character(wkt))
}

Expand Down

0 comments on commit deaef38

Please sign in to comment.