Enhance SDA_spatialQuery(what="geom") to query from either SSURGO or (new) STATSGO?

[dschlaep]

The function SDA_spatialQuery (soilDB v2.5.7) currently returns geometries and mukey values from SSURGO and not from STATSGO. As long as SSURGO remains incomplete and SDA doesn't connect to merged version such as gNATSGO, it would be great to have the functionality to query from either database for a complete workflow coverage.

A call with the argument what = "mukey" queries the tabular table mapunit and relies on a SQL statement based on SDA_Get_Mukey_from_intersection_with_WktWgs84() to handle the spatial conversion; this, however, restricts the output to SSURGO. Until there is a comparable function for STATSGO not much can be done for this case apparently.

However, a call with the argument what = "geom" relies on T-SQL intrinsic spatial function STIntersects that queries the spatial table mupolygon of SSURGO. According to a comment from Dec 2017 made in response to an earlier, different issue, a corresponding table exists for STATSGO #36 (comment), i.e., gsmmupolygon.

I created a local copy of SDA_spatialQuery for my own use with a new argument db = c("SSURGO", "STATSGO") and few minor code changes, which appear to work well (see below).

As a simple user, I don't know whether querying gsmmupolygon is a supported activity. If it is, however, then it would be great if you would consider whether the suggested enhancement could be a good fit for soilDB! Let me know if you would like me to create a pull-request with the suggestion incorporated.

Thanks!

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

  db <- match.arg(db)

  if (what == "mukey" && db == "STATSGO") {
    stop("query type 'mukey' for 'STATSGO' is not supported", call. = FALSE)
  }
  
  ...

  if(what == 'geom') {
    # as suggested in https://github.com/ncss-tech/soilDB/issues/36#issuecomment-349054067
    db_table <- switch(db, SSURGO = "mupolygon", STATSGO = "gsmmupolygon")

    # return intersection
    if(geomIntersection) {
      q <- sprintf("
               SELECT
                 mupolygongeo.STIntersection( geometry::STGeomFromText('%s', 4326) ).STAsText() AS geom, P.mukey
                 FROM %s AS P
                 WHERE mupolygongeo.STIntersects( geometry::STGeomFromText('%s', 4326) ) = 1;", wkt, db_table, wkt)
    } else {
      # return overlapping
      q <- sprintf("
               SELECT
                 mupolygongeo.STAsText() AS geom, P.mukey
                 FROM %s AS P
                 WHERE mupolygongeo.STIntersects( geometry::STGeomFromText('%s', 4326) ) = 1;", db_table, wkt)
    }
  ...
  ...

[dylanbeaudette]

Excellent idea! I would be happy to review a PR and merge ASAP.

Dylan

[dschlaep]

Wonderful, thanks! Working on it.

Do you have a preference for how to name the new argument? Not sure whether there is a standard short terminology for "Soil Geographic Databases". I used db in my code snippet above

[dylanbeaudette]

Thanks for the input, suggestions, and PR!

I've made some updates since the PR in the documentation and "SDA spatial" code (167200a). I'm checking on the efficiency question related to SDA_Get_Mukey_from_intersection_with_WktWgs84, it could be that we can skip the wrapper and hit the base tables so that SDA_spatialQuery(what = 'mukey', db = 'STATSGO) works as expected.

Also, I've added a new section to the examples demonstrating how to get both SSURGO / STATSGO geometry.

I've also deprecated some older interfaces to SDA, so it is probably time to re-visit the related tutorials. Honestly, all of the tutorials should be properly organized via bookdown.