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:
Load SSURGO spatial data
Query Soil Data Access (SDA) for dominant ecological site assignments
Join tabular and spatial data
Export and visualize the results
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.
You can use SDA_spatialQuery()
and
fetchSDA_spatial()
to retrieve spatial data directly from
Soil Data Access (SDA), bypassing the need for local files.
SDA_spatialQuery()
to get mukey polygons
for an AOIDefine 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.
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 sourcesfetchSDA_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)
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.
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
).
As we did above for fetchSDA_spatial()
we are going to
get the unique set of map units in our AOI.
get_SDA_coecoclass()
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.
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
).
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.
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
).
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
.
Finally, we can write our result out to a spatial file, such as a
GeoPackage, using sf
or terra
.