Get Spatial Data from Soil Data Access by mukey, nationalmusym or areasymbol

This method facilitates queries to Soil Data Access (SDA) mapunit and survey area geometry. Queries are generated based on map unit key (mukey) and national map unit symbol (nationalmusym) for mupolygon (SSURGO) or gsmmupolygon (STATSGO) geometry OR legend key (lkey) and area symbols (areasymbol) for sapolygon (Soil Survey Area; SSA) geometry).

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.

fetchSDA_spatial(
  x,
  by.col = "mukey",
  method = "feature",
  geom.src = "mupolygon",
  db = "SSURGO",
  add.fields = NULL,
  chunk.size = 10,
  verbose = TRUE,
  as_Spatial = getOption("soilDB.return_Spatial", default = FALSE)
)

Arguments

x: A vector of map unit keys (mukey) or national map unit symbols (nmusym) for mupolygon geometry OR legend keys (lkey) or soil survey area symbols (areasymbol) for sapolygon geometry
by.col: Column name containing map unit identifier "mukey", "nmusym"/"nationalmusym" for geom.src mupolygon OR "areasymbol", "areaname", "mlraoffice", "mouagncyresp" for geom.src sapolygon; default is determined by is.numeric(x) TRUE for mukey or lkey and nationalmusym or areasymbol otherwise.
method: geometry result type: "feature" returns polygons, "bbox" returns the bounding box of each polygon (via STEnvelope()), and "point" returns a single point (via STPointOnSurface()) within each polygon.
geom.src: Either mupolygon (map unit polygons) or sapolygon (soil survey area boundary polygons)
db: Default: "SSURGO". When geom.src is mupolygon, use STATSGO polygon geometry instead of SSURGO by setting db = "STATSGO"
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.
chunk.size: Number of values of x to process per query. Necessary for large results. Default: 10
verbose: Print messages?
as_Spatial: Return sp classes? e.g. Spatial*DataFrame. Default: FALSE.

Value

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.

Querying regions with complex mapping may require smaller chunk.size. Numerically adjacent IDs in the input vector may share common qualities (say, all from same soil survey area or region) which could cause specific chunks to perform "poorly" (slow or error) no matter what the chunk size is. Shuffling the order of the inputs using sample() may help to eliminate problems related to this, depending on how you obtained your set of MUKEY/nationalmusym to query. One could feasibly use muacres as a heuristic to adjust for total acreage within chunks.

Note that STATSGO data are fetched where CLIPAREASYMBOL = 'US' to avoid duplicating state and national subsets of the geometry.

Author

Andrew G. Brown, Dylan E. Beaudette

Examples

# \donttest{
 
    # get spatial data for a single mukey
    single.mukey <- try(fetchSDA_spatial(x = "2924882"))
#> Using 1 chunks...
#> Chunk #1 completed (n = 1; 0.1 secs)
#> Done in 0.1 secs; mean/chunk: 0.1 secs; mean/symbol: 0.13 secs.

    # demonstrate fetching full extent (multi-mukey) of national musym
    full.extent.nmusym <- try(fetchSDA_spatial(x = "2x8l5", by = "nmusym"))
#> Using 1 chunks...
#> Chunk #1 completed (n = 3; 0.8 secs)
#> Done in 1 secs; mean/chunk: 0.8 secs; mean/symbol: 0.33 secs.

    # compare extent of nmusym to single mukey within it
    if (!inherits(single.mukey, 'try-error') && 
        !inherits(full.extent.nmusym, 'try-error')) {
        
        if (requireNamespace("sf")) {
      
         plot(sf::st_geometry(full.extent.nmusym), col = "RED", border = 0)
         plot(sf::st_geometry(single.mukey), add = TRUE, col = "BLUE", border = 0)
       
        }
        
    }


    # demo adding a field (`muname`) to attribute table of result
    head(try(fetchSDA_spatial(x = "2x8l5", by="nmusym", add.fields="muname")))
#> Using 1 chunks...
#> Chunk #1 completed (n = 3; 0.6 secs)
#> Done in 0.8 secs; mean/chunk: 0.6 secs; mean/symbol: 0.26 secs.
#> Simple feature collection with 6 features and 4 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -121.034 ymin: 38.01706 xmax: -120.9596 ymax: 38.24938
#> Geodetic CRS:  WGS 84
#>    mukey areasymbol nationalmusym                                        muname
#> 1 462101      CA077         2x8l5 Pentz-Bellota complex, 2 to 15 percent slopes
#> 2 462101      CA077         2x8l5 Pentz-Bellota complex, 2 to 15 percent slopes
#> 3 462101      CA077         2x8l5 Pentz-Bellota complex, 2 to 15 percent slopes
#> 4 462101      CA077         2x8l5 Pentz-Bellota complex, 2 to 15 percent slopes
#> 5 462101      CA077         2x8l5 Pentz-Bellota complex, 2 to 15 percent slopes
#> 6 462101      CA077         2x8l5 Pentz-Bellota complex, 2 to 15 percent slopes
#>                             geom
#> 1 POLYGON ((-121.0226 38.2276...
#> 2 POLYGON ((-121.0141 38.2325...
#> 3 POLYGON ((-121.0339 38.2325...
#> 4 POLYGON ((-120.9955 38.0674...
#> 5 POLYGON ((-121.0102 38.1891...
#> 6 POLYGON ((-120.9648 38.0208...
# }