Mapping Dominant Ecological Sites using SSURGO Data

Overview

This vignette demonstrates how to extract and map dominant ecological site information for a given Area of Interest (AOI) using the soilDB package and open-source geospatial tools in R.

For a general introduction to geospatial data in R, see the sf package documentation, terra package documentation, the Geocomputation with R book and the Spatial Data Science with R and terra book.

We will:

Prerequisites

In this vignette we will use sf functions and object types for processing and storing spatial data. This package is used under the hood by soilDB for processing SSURGO data. The terra package can be used with only minor changes in syntax.

library(soilDB)
library(sf)

Load SSURGO Spatial Data

Using Soil Data Access

You can use SDA_spatialQuery() and fetchSDA_spatial() to retrieve spatial data directly from Soil Data Access (SDA), bypassing the need for local files.

Option 1: Use SDA_spatialQuery() to get mukey polygons for an AOI

Define Area of Interest (AOI) as a bounding box sf POLYGON.

aoi <- sf::st_as_sfc(sf::st_bbox(c(
  xmin = -120.85,
  xmax = -120.775,
  ymin = 37.975,
  ymax = 38.0
), crs = 4326))

# Query SDA for map unit polygons cropped to the AOI
soil_polygons <- SDA_spatialQuery(
  aoi, 
  what = "mupolygon"
)

plot(
  soil_polygons["mukey"],
  main = "SSURGO Map Units from SDA_spatialQuery"
)
head(soil_polygons)

You can also set geomIntersection = TRUE so that intersecting geometries are cropped to the AOI. This is convenient if you have a very specific AOI in mind or would like to reduce the amount of data you are going to download.

# Query SDA for map unit polygons cropped to the AOI
soil_polygons <- SDA_spatialQuery(
  aoi, 
  what = "mupolygon", 
  geomIntersection = TRUE
)

plot(
  soil_polygons["mukey"],
  main = "Cropped SSURGO Map Units from SDA_spatialQuery"
)
head(soil_polygons)

Option 2: Use fetchSDA_spatial()

Next, we use fetchSDA_spatial().

This function is different because instead of using an area of interest, it takes a vector of map unit or legend identifiers.

Here, we take the unique map unit keys from the Option 1 result and return the full extent of those map units (not just the intersecting polygons). This is particularly helpful when doing studies that involve the full extent of map unit concepts, as opposed to more site-specific analyses.

mu_ssurgo <- fetchSDA_spatial(
  unique(soil_polygons$mukey), 
  by.col = "mukey",
  add.fields = c("legend.areaname", 
                 "mapunit.muname", 
                 "mapunit.farmlndcl")
)

plot(mu_ssurgo["mukey"], main = "SSURGO Map Units from fetchSDA_spatial")
head(mu_ssurgo)

Note that using add.fields we have included some additional contextual information: area name, map unit name, and map unit farmland classification.

fetchSDA_spatial() geometry sources

fetchSDA_spatial() geom.src argument can be used to return SSURGO map unit polygons and survey area polygons, STATSGO mapunit polygons, and MLRA polygons.

Here is an example using STATSGO. We use SDA_spatialQuery() to fetch the data for our AOI, then fetchSDA_spatial() to get the full extent of those map unit concepts.

statsgo_polygons <- SDA_spatialQuery(
  aoi, 
  what = "mupolygon",
  db = "STATSGO",
  geomIntersection = TRUE
)

plot(statsgo_polygons["mukey"], main = "STATSGO Map Units from SDA_spatialQuery")
head(statsgo_polygons)

mu_statsgo <- fetchSDA_spatial(
  unique(statsgo_polygons$mukey), 
  by.col = "mukey",
  db = "STATSGO",
  add.fields = c("legend.areaname", "mapunit.muname", "mapunit.farmlndcl")
)

plot(mu_statsgo["mukey"], main = "STATSGO Map Units from fetchSDA_spatial")
head(mu_statsgo)
ssas <- SDA_spatialQuery(aoi, what = "areasymbol")
ssas

ssa <- fetchSDA_spatial(
  ssas$areasymbol, 
  by.col = "areasymbol",
  geom.src = "sapolygon",
  add.fields = c("legend.areaname")
)

plot(ssa["areasymbol"], main = "SSURGO Soil Survey Area from SDA")
head(ssa)

SDA_spatialQuery() vs. fetchSDA_spatial()

  • SDA_spatialQuery() is ideal for spatial queries where you have a specific, possibly complex, area of interest.

  • fetchSDA_spatial() returns the full extent of the specified map unit concepts, optionally including more legend and map unit attributes via add.fields argument.

Local Data Sources

If working with large extents, it is generally better to be use a local data source. You can download SSURGO data from Web Soil Survey using downloadSSURGO(). Then, you can create GeoPackage or other database types using createSSURGO() or the SSURGOPortal tools. This process is described in detail in the Local SSURGO Databases vignette. See also the SSURGOPortal R package.

Assuming you have a GeoPackage from SSURGOPortal ("soilmu_a.gpkg"), then you can read it with sf::st_read() or terra::vect()

ssurgo_path <- "data/soilmu_a.gpkg"  # Replace with your actual path
soil_polygons <- sf::st_read(ssurgo_path)
head(soil_polygons)

This soil_polygons object should have the standard set of columns we would expect for a SSURGO map unit data source, including: map unit key (mukey), area symbol (areasymbol), map unit symbol (musym).

Query SDA for Dominant Ecological Site Info

Extract Unique Map Unit Keys

As we did above for fetchSDA_spatial() we are going to get the unique set of map units in our AOI.

mukeys <- unique(soil_polygons$mukey)

Use get_SDA_coecoclass()

eco_data <- get_SDA_coecoclass(mukeys = mukeys)
head(eco_data)

This function returns a data frame including mukey, cokey, ecoclassid, and ecoclassname.

There are several methods for aggregating from component to map unit level available in get_SDA_coecoclass(). The default aggregation method is "none" which will return 1 record per map unit component, so many map units will have more than one record.

Using method="dominant component"

Most often, users want “typical” conditions that can apply to a whole map unit.

Sometimes it is helpful to be able to point to a specific component, so that you do not have to reason over mathematical aggregation of distinct components within map unit concepts. The most common method to select one component per map unit is "dominant component"

eco_data_domcond <- get_SDA_coecoclass(
  mukeys = mukeys,
  method = "dominant component"
)
head(eco_data_domcond)

Note that this output includes information for the dominant component (comppct_r and compname).

Using method="dominant condition"

The method "dominant condition" is convenient because it accounts for the possibility that multiple components have the same class (ecological site) and can have their component percentages summed.

For a property like ecological site, it is fairly common for multiple components in a map unit to have the same site assigned, so the dominant condition makes up a higher percentage of the area of the map unit than the dominant component alone.

eco_data_domcond <- get_SDA_coecoclass(
  mukeys = mukeys,
  method = "dominant condition"
)
head(eco_data_domcond)

Note that this result includes the aggregate ecological site composition ecoclasspct_r.

The cokey and comppct_r are from the dominant component can be used to link to a specific, common component for other more detailed information.

Join Tabular and Spatial Data

We want to have a 1:1 relationship between our map unit polygons and the thematic variable we are mapping, so we will use the map unit dominant condition ecological sites (eco_data_domcond).

soil_polygons <- merge(
  soil_polygons, 
  eco_data_domcond, 
  by = "mukey"
)

Visualize the Result

plot(
  soil_polygons["ecoclassid"],
  main = "Dominant Ecological Site by Map Unit"
)

This same process for merging in aggregate map unit level information generalizes to other SSURGO tables. You can easily create thematic maps from any data you can aggregate to the point that it is 1:1 with the mukey.

Export Spatial File

Finally, we can write our result out to a spatial file, such as a GeoPackage, using sf or terra.

sf::st_write(
  soil_polygons, 
  "ecosite_dominant.gpkg", 
  delete_dsn = TRUE
)

References

SSURGO Data and Tools

soilDB functions

Spatial Analysis in R