This vignette will explore how to use the soilDB package to create thematic maps of the Soil Survey Geographic Database (SSURGO) via Soil Data Access (SDA) and SoilWeb coverage services.
A web coverage service (WCS) is provided for the gSSURGO and gNATGSO
map unit
key grids by the UCDavis California Soil
Resource Lab SoilWeb server. These grids are a raster representation
of the gSSURGO and gNATSGO map unit keys (mukey) for the
conterminous United States at a resolution of 30m. The GeoTIFF format is
used to ensure maximum compatibility. Cell values are map unit keys,
encoded as unsigned 32-bit integers. The standard spatial reference for
the grids is Albers Equal Area Conic (NAD83) coordinate reference system
("EPSG:5070"). The grids are ‘LZW’ compressed and
internally tiled for efficient random access. The grid topology and cell
values are identical to the rasters contained within the gSSURGO and
gNATSGO File Geodatabases (FGDB).
The files backing each WCS are re-created annually after the SSURGO annual refresh on October 1st. They are typically in-sync with the official version of the data hosted by Soil Data Access by early November. Map unit keys can change over time, especially in soil survey areas that were updated during the last fiscal year.
In addition to the standard CONUS gSSURGO and gNATSGO grids, these web coverage services are also provided:
Raster Soil Survey (RSS) products at 10m resolution (EPSG:5070).
STATSGO2 (2016) at 300m resolution (EPSG:5070).
SSURGO in HI (EPSG:6628) and PR (EPSG:32161) at 30m resolution.
Get the latest CRAN version of soilDB,
terra and sf for the following examples.
terra is required for handling of raster data.
sf or terra may be used for handling vector
data for as inputs to web coverage service functions, and examples
demonstrate usage of both packages for this purpose. In general for
soilDB functions that take spatial inputs, if
terra objects are used as input, terra objects
are returned.
Consider also installing the latest development versions from GitHub or r-universe:
install.packages(c('soilDB', 'terra', 'sf'),
repos = c('https://ncss-tech.r-universe.dev',
'https://rspatial.r-universe.dev',
'https://r-spatial.r-universe.dev')
)Here are some basic usage examples for the coverage services that can
be accessed with mukey.wcs().
See the mukey.wcs
manual pages for details.
# select gSSURGO grid, 30m resolution
x <- mukey.wcs(aoi = aoi, db = 'gssurgo', ...)
# select gNATSGO grid, 30m resolution
x <- mukey.wcs(aoi = aoi, db = 'gnatsgo', ...)
# select RSS grid, 10m resolution
x <- mukey.wcs(aoi = aoi, db = 'RSS', ...)
# select STATSGO2 grid, 300m resolution
x <- mukey.wcs(aoi = aoi, db = 'statsgo', ...)ISSR800.wcs
is a similar function that makes available a variety of pre-aggregated
800m resolution properties derived from gNATSGO.
This is where we will link to a more detailed ISSR800 vignette in the future.
Excerpt from the gSSURGO documentation.
The gSSURGO Database is derived from the official Soil Survey Geographic (SSURGO) Database. SSURGO generally has the most detailed level of soil geographic data developed by the National Cooperative Soil Survey (NCSS) in accordance with NCSS mapping standards. The tabular data represent the soil attributes and are derived from properties and characteristics stored in the National Soil Information System (NASIS). The gSSURGO data were prepared by merging the traditional vector-based SSURGO digital map data and tabular data into statewide extents, adding a statewide gridded map layer derived from the vector layer, and adding a new value-added look up table (Valu1) containing “ready to map” attributes. The gridded map layer is a file geodatabase raster in an ArcGIS file geodatabase. The raster and vector map data have a statewide extent. The raster map data have a 10-meter cell size that approximates the vector polygons in an Albers Equal Area projection. Each cell (and polygon) is linked to a map unit identifier called the map unit key. A unique map unit key is used to link the raster cells and polygons to attribute tables. Due to file size, the raster layer for the conterminous United States is only available in a 30-meter resolution.
Excerpt from the gNATSGO documentation.
The gNATSGO databases contain a raster of the soil map units and 70 related tables of soil properties and interpretations. They are designed to work with the SPSD gSSURGO ArcTools. Users can create full coverage thematic maps and grids of soil properties and interpretations for large geographic areas, such as the extent of a State or the conterminous United States. Please note that the State-wide geodatabases contain a 10 meter raster and the CONUS database contains a 30 meter raster.
The gNATSGO database is composed primarily of SSURGO data, but STATSGO2 data was used to fill in the gaps. The RSSs are newer product with relatively limited spatial extent. These RSSs were merged into the gNATSGO after combining the SSURGO and STATSGO2 data. The extent of RSS is expected to increase in the coming years.
Excerpt from the RSS documentation.
Raster Soil Survey is a reference to the products of soil survey work completed using digital soil mapping methodologies. Digital soil mapping is the production of georeferenced soil databases based on the quantitative relationships between soil measurements made in the field or laboratory and environmental data and may be represented as either discrete classes or continuous soil properties. Both digital and traditional soil mapping use a conceptual soil-landscape model as a means for organizing environmental information into discrete divisions. The primary difference between these two approaches is that digital methods exploit quantitative relationships of the environmental information, while traditional methods utilize a more subjective approach and the approximate relationships of the environmental information to spatially represent where the divisions are represented.
An experimental, 300m gridded representation of STATSGO 2 is provided by the SoilWeb web coverage service. This is not an official USDA-NRCS product.
Excerpt from STATSGO2 Documentation.
The Digital General Soil Map of the United States or STATSGO2 is a broad-based inventory of soils and non-soil areas that occur in a repeatable pattern on the landscape and that can be cartographically shown at the scale mapped of 1:250,000 in the continental U.S., Hawaii, Puerto Rico, and the Virgin Islands and 1:1,000,000 in Alaska. The level of mapping is designed for broad planning and management uses covering state, regional, and multi-state areas. The U.S. General Soil Map is comprised of general soil association units and is maintained and distributed as a spatial and tabular dataset.
Thematic mapping or analysis of soil information requires connecting
the grids to our tabular data sources, either using local files or Soil
Data Access (SDA) web-service. The soilDB package provides
many convenient interfaces to SDA. Note that SDA does not yet contain
tabular data for the raster soil surveys.
A buffer applied to a single WGS84 coordinate can be used to create a bounding box (BBOX):
We can create a terra SpatVector containing the single point, then
create a 1000m radius circular polygon around the point with
buffer():
library(terra)
library(soilDB)
# example point, WGS84 coordinates
p <- vect(
data.frame(
lon = -118.55639,
lat = 36.52578
),
crs = "EPSG:4326"
)
# 1000m buffer applied to WGS84 coordinate
# radius defined in meters
b <- buffer(p, 1000)
# query WCS
# result is in EPSG:5070
mu <- mukey.wcs(b, db = 'gSSURGO')
# inspect
plot(mu, legend = FALSE, axes = FALSE, main = metags(mu)['description'])
# add buffer, after transforming to mukey grid CRS
plot(project(b, "EPSG:5070"), add = TRUE)
# add original point, after transforming to mukey grid CRS
plot(project(p, "EPSG:5070"), add = TRUE, pch = 16)Sometimes it is convenient to specify a BBOX created from a website or single point specified in WGS84 coordinates. Arbitrary spatial objects can be used as input. SoilWeb provides two keyboard shortcuts in the map interface:
Press 'b' to copy the current bounding box
coordinates, returned as: -118.6609
36.4820,-118.6609 36.5972,-118.3979 36.5972,-118.3979 36.4820,-118.6609
36.4820. These five coordinates pairs form a rectangular
polygon, the geometry of which can be represented as Well-Known
Text.
Press 'p' copy link to center coordinate, returned
as:
https://casoilresource.lawr.ucdavis.edu/gmap/?loc=36.53964,-118.52943,z13
Right-clicking anywhere in the map interface will also generate a link to those coordinates and zoom level.
library(sf)
library(soilDB)
library(terra)
# paste the five coordinates comprising the BBOX polygon here
bb <- '-118.6609 36.4820,-118.6609 36.5972,-118.3979 36.5972,-118.3979 36.4820,-118.6609 36.4820'
# convert WKT string -> sfc geometry
wkt <- sprintf('POLYGON((%s))', bb)
x <- st_as_sfc(wkt)
# set coordinate reference system as GCS/WGS84
st_crs(x) <- 4326
# query WCS
mu <- mukey.wcs(x, db = 'gSSURGO')
# inspect
plot(mu, legend = FALSE, axes = FALSE, main = metags(mu)['description'])
# add original BBOX, after transforming to mukey grid CRS
plot(st_transform(x, 5070), add = TRUE)Use the mukey.wcs() function to access chunks of the
CONUS “map unit key” grids. The area of interest (AOI) can be defined
manually, as below, or automatically extracted from sf,
sfc, bbox, SpatRaster,
SpatVector, Spatial* or
RasterLayer objects.
The resulting grid of integers (mukey; or map unit key)
represents unique map units within specific soil survey areas; the grid
isn’t all that useful by itself.
To make the grids more useful, join data from Soil Data Access (SDA)
or local files to create thematic maps based on the map unit key. SSURGO
data for specific soil components, depths, and properties within a map
unit can be aggregated so the aggregate values are 1:1 with
mukey, then the resulting values can be used to symbolize
the map.
# make a bounding box and assign a CRS (4326: GCS, WGS84)
a <- st_bbox(
c(xmin = -114.16, xmax = -114.08, ymin = 47.65, ymax = 47.68),
crs = st_crs(4326)
)
# fetch gSSURGO map unit keys at native resolution (30m)
mu <- mukey.wcs(aoi = a, db = 'gssurgo')
# check:
print(mu)
plot(
mu,
main = 'gSSURGO map unit keys',
sub = 'Albers Equal Area Projection',
axes = FALSE,
legend = FALSE
)It is possible to retrieve vector geometries from SDA with the
SDA_spatialQuery() function. Standard SSURGO map unit
polygons can be obtained in addition to soil survey area polygons and
STATSGO map unit polygons.
The vector data are stored and delivered in a geographic coordinate
system (WGS84), whereas the WCS grids generally use a locally relevant
projected coordinate system ("EPSG:5070" in CONUS).
Overlaying SSURGO polygons and map unit key grids will therefore require
a simple transformation.
First get intersecting SSURGO polygons from SDA with
SDA_spatialQuery()
# because mu is a SpatRaster, result is a SpatVector object (GCS WGS84)
p <- SDA_spatialQuery(mu, what = 'mupolygon', geomIntersection = TRUE)Then transform to AEA coordinate reference system used by CONUS
gSSURGO / gNATSGO ("EPSG:5070").
Inspect the result by overlaying SSURGO polygons on the 30m map unit key grid.
Requesting map unit key grids at a resolution other than 30m is possible, but only suitable for a quick “preview” of the data.
For example, it is possible to get a larger areal extent of data by requesting grids at coarser resolution (e.g. 800m).
However, the pixels represent categories (unique map units) and are selected by nearest-neighbor; there is no other generalization used to convert the source 30m grid to the coarser scale. A coarser representation of data can be used for inspection of general patterns. Detailed analysis should be based on derived property data sets aggregated up from 30m results.
# make a bounding box (in California) and assign a CRS (GCS WGS84 / EPSG:4326)
a.CA <- st_bbox(c(
xmin = -121,
xmax = -120,
ymin = 37,
ymax = 38
), crs = st_crs(4326))
# fetch gSSURGO map unit keys at ~800m
# nearest-neighbor resampling = this is a "preview"
# result is a SpatRaster object
x.800 <- mukey.wcs(aoi = a.CA, db = 'gssurgo', res = 800)
plot(
x.800,
main = 'A Preview of gSSURGO Map Unit Keys',
sub = 'Albers Equal Area Projection (800m)\nnearest-neighbor resampling',
axes = FALSE,
legend = FALSE
)The specific RSS (state-level) data sets can be downloaded (map unit
key grids, tabular data) on Box:
https://nrcs.app.box.com/v/soils. Please note that tabular
data for Raster Soil Surveys are not yet available via Soil Data
Access.
# Coweeta Hydrologic Laboratory extent; specified in EPSG:5070
a <- st_bbox(
c(xmin = 1129000, xmax = 1135000, ymin = 1403000, ymax = 1411000),
crs = st_crs(5070)
)
# convert boundary sf polygon
a <- st_as_sfc(a)
# gSSURGO grid: 30m resolution
(x <- mukey.wcs(a, db = 'gSSURGO', res = 30))
# gNATSGO grid: 30m resolution
(y <- mukey.wcs(a, db = 'gNATSGO', res = 30))
# RSS grid: 10m resolution
(z <- mukey.wcs(a, db = 'RSS', res = 10))
# graphical comparison
par(mfcol = c(1, 3))
# gSSURGO
plot(
x,
axes = FALSE,
legend = FALSE,
main = metags(x)['description']
)
plot(a, add = TRUE)
# gNATSGO
plot(
y,
axes = FALSE,
legend = FALSE,
main = metags(y)['description']
)
plot(a, add = TRUE)
# RSS
plot(
z,
axes = FALSE,
legend = FALSE,
main = metags(z)['description'],
ext = x
)
plot(a, add = TRUE)Continuing from the example above, we can use
db='statsgo' to compare gSSURGO product with the Digital
General Soil Map of the United States (STATSGO2). STATSGO data are
provided at 10x the nominal resolution of gSSURGO (300m v.s. 30m) to
reflect the relative generality of this product.
A new 30m SSURGO map unit key WCS based on the
"EPSG:6628" coordinate reference system has been added for
Hawaii.
The example bounding box is centered on the southern coast of Kauai.
# paste your BBOX text here
bb <- '-159.7426 21.9059,-159.7426 22.0457,-159.4913 22.0457,-159.4913 21.9059,-159.7426 21.9059'
# convert WKT string -> sfc geometry
wkt <- sprintf('POLYGON((%s))', bb)
x <- st_as_sfc(wkt, crs = 4326)
# query WCS
mu <- mukey.wcs(x, db = 'hi_ssurgo')
# make NA (the ocean) blue
plot(
mu,
legend = FALSE,
axes = FALSE,
main = metags(mu)['description'],
colNA = 'royalblue'
)A new 30m SSURGO map unit key WCS based on the
"EPSG:32161" coordinate reference system has been added for
Puerto Rico.
The example bounding box is centered on the eastern coast of Puerto Rico.
# paste your BBOX text here
bb <- '-65.7741 18.1711,-65.7741 18.3143,-65.5228 18.3143,-65.5228 18.1711,-65.7741 18.1711'
# convert WKT string -> sfc geometry
wkt <- sprintf('POLYGON((%s))', bb)
x <- st_as_sfc(wkt, crs = 4326)
# query WCS
mu <- mukey.wcs(x, db = 'pr_ssurgo')
# make missing data (NA; the ocean) blue
plot(
mu,
legend = FALSE,
axes = FALSE,
main = metags(mu)['description'],
colNA = 'royalblue'
)The following example BBOX + resulting gSSURGO map unit key grid will be used for thematic mapping examples:
# make a bounding box and assign a CRS (4326: GCS, WGS84)
a <- st_bbox(
c(xmin = -114.16, xmax = -114.08, ymin = 47.65, ymax = 47.68),
crs = st_crs(4326)
)
# convert bbox to sf geometry
a <- st_as_sfc(a)
# fetch gSSURGO map unit keys at native resolution (~30m)
mu <- mukey.wcs(aoi = a, db = 'gssurgo')The “Mapunit Aggregate Attribute” table records a variety of soil attributes and interpretations that have been aggregated from the component level to a single value at the map unit level. They have been aggregated by one or more appropriate means in order to express a consolidated value or interpretation for the map unit as a whole.
Use the get_SDA_muaggatt() function, or write a query in
SQL and submit via SDA_query().
# copy example grid
mu2 <- mu
# extract raster attribute table for thematic mapping
(rat <- cats(mu2)[[1]])
# optionally use convenience function:
# * returns all fields from muaggatt table
# * along with map unit name
# tab <- get_SDA_muaggatt(mukeys = as.numeric(rat$mukey), query_string = TRUE)
.sql <- paste0(
"SELECT mukey, aws050wta, aws0100wta FROM muaggatt WHERE mukey IN ",
format_SQL_in_statement(as.numeric(rat$mukey))
)
# run query, result is a data.frame
tab <- SDA_query(.sql)
# check
head(tab)
# set raster categories
levels(mu2) <- tab
# convert grid + RAT -> stack of property grids
aws <- catalyze(mu2)
# plot, set a common range [0, 20] for both layers
plot(
aws,
axes = FALSE,
cex.main = 0.7,
main = c(
'Plant Available Water Storage (cm)\nWeighted Mean over Components, 0-50cm',
'Plant Available Water Storage (cm)\nWeighted Mean over Components, 0-100cm'
),
range = c(0, 20)
)You can use the get_SDA_interpretation() function to
return interpretation ratings for specific components, or aggregated up
to the map unit level. Here, we produce maps using two interpretation
rules: 'ENG - Construction Materials; Roadfill', and
'AWM - Irrigation Disposal of Wastewater'.
# copy example grid
mu2 <- mu
# extract RAT for thematic mapping
rat <- cats(mu2)[[1]]
rules <- c('ENG - Construction Materials; Roadfill',
'AWM - Irrigation Disposal of Wastewater')
tab <- get_SDA_interpretation(
rulename = rules,
method = "Weighted Average",
mukeys = as.numeric(rat$mukey)
)
# check
head(tab)
# set ordered factor levels (for nice label/legend order)
tab$class_ENGConstructionMaterialsRoadfill <- factor(
tab$class_ENGConstructionMaterialsRoadfill,
levels = c(
'Not suited',
'Poorly suited',
'Moderately suited',
'Moderately well suited',
'Well suited',
'Not Rated'
),
ordered = TRUE
)
par(mar = c(4, 12, 3, 3))
boxplot(
rating_ENGConstructionMaterialsRoadfill ~ class_ENGConstructionMaterialsRoadfill,
cex.main = 0.7,
main = 'ENG - Construction Materials; Roadfill',
ylab = "",
data = tab,
horizontal = TRUE, # fuzzy ratings on X axis
las = 1 # rotate axis labels 90 degrees
)From above graph we can see that the different suitability rating
classes class_ENGConstructionMaterialsRoadfill each
correspond to a range of fuzzy values
(rating_ENGConstructionMaterialsRoadfill).
Next, we can view the ratings as a thematic map:
vars <- c(
'rating_ENGConstructionMaterialsRoadfill',
'rating_AWMIrrigationDisposalofWastewater'
)
# set raster categories
levels(mu2) <- tab[, c('mukey', vars)]
rating <- catalyze(mu2)
# inspect
plot(
rating,
axes = FALSE,
cex.main = 0.7,
main = c(
'Construction Materials; Roadfill\nWeighted Mean over Components',
'Irrigation Disposal of Wastewater\nWeighted Mean over Components'
)
)Soil “components” are the members (land types) of a map unit. A map unit may contain several distinct soil and non-soil areas.
A widely-used property which is calculated as a standard part of
SSURGO soil map unit components is the “steel corrosion potential”. Here
we also use get_SDV_legend_elements() to get the standard
Soil Data Viewer colors for the selected property/interpretation.
# copy example grid
mu2 <- mu
# extract RAT for thematic mapping
rat <- cats(mu2)[[1]]
tab <- get_SDA_property(property = 'Corrosion of Steel',
method = 'DOMINANT CONDITION',
mukeys = as.integer(rat$mukey),
miscellaneous_areas = TRUE)
# get soil data viewer standard colors for corsteel
cols <- get_SDV_legend_elements("attributecolumnname = 'corsteel'")
# set raster categories
levels(mu2) <- tab[, c('mukey', 'corsteel')]
# set active category
activeCat(mu2) <- 'corsteel'
# plot
plot(
mu2,
col = cols$hex[na.omit(match(unique(tab$corsteel), cols$label))],
axes = FALSE,
legend = "topleft"
)Another example is thematic mapping of the “simplified component parent material group”. First, set up a new AOI for the following examples:
# https://casoilresource.lawr.ucdavis.edu/gmap/?loc=36.57666,-96.70175,z14
# make a bounding box and assign a CRS (4326: GCS, WGS84)
a <- st_bbox(
c(xmin = -96.7696, xmax = -96.6477,
ymin = 36.5477, ymax = 36.6139),
crs = st_crs(4326)
)
# fetch gSSURGO map unit keys at native resolution (~30m)
mu <- mukey.wcs(aoi = a, db = 'gssurgo')
plot(
mu,
legend = FALSE,
axes = FALSE,
cex.main = 0.7,
main = 'gSSURGO Map Unit Key Grid'
)We use get_SDA_pmgroupname() to obtain the tabular
parent material information to relate to map unit keys:
# copy example grid
mu2 <- mu
# extract RAT for thematic mapping
rat <- cats(mu2)[[1]]
# simplified parent material group name
tab <- get_SDA_pmgroupname(mukeys = as.integer(rat$mukey),
miscellaneous_areas = TRUE)
# set raster categories
levels(mu2) <- tab[, c('mukey', 'pmgroupname')]
# set active category
activeCat(mu2) <- 'pmgroupname'
plot(mu2, legend = "topleft", axes = FALSE)We can also inspect a mapunit-level hydric rating derived from the
default aggregation method in get_SDA_hydric().
# copy example grid
mu2 <- mu
# extract RAT for thematic mapping
rat <- cats(mu2)[[1]]
# simplified parent material group name
tab <- get_SDA_hydric(mukeys = as.integer(rat$mukey))
levels(mu2) <- tab[, c('mukey', 'HYDRIC_RATING')]
# set active category
activeCat(mu2) <- 'HYDRIC_RATING'
plot(mu2, legend = "topleft", axes = FALSE)The get_SDA_property()
function from soilDB is a general interface to aggregated
SSURGO/STATSGO tabular data via SDA. It was used to obtain a
component-level property previously (steel corrosion), and now we will
use it to aggregate several horizon-level property values for a specific
depth interval.
Derive aggregate soil properties, merge with raster attribute table (RAT).
# extract RAT for thematic mapping
rat <- cats(mu)[[1]]
# variables of interest
vars <- c("dbthirdbar_r", "awc_r", "ph1to1h2o_r")
# get / aggregate specific horizon-level properties from SDA
# be sure to see the manual page for this function
tab <- get_SDA_property(property = vars,
method = "Dominant Component (Numeric)",
mukeys = as.integer(rat$mukey),
top_depth = 0,
bottom_depth = 25)
# check
head(tab)
# convert areasymbol into a factor easy plotting later
tab$areasymbol <- factor(tab$areasymbol)
# set raster categories
levels(mu) <- tab[, c('mukey', vars)]
# list variables in the RAT
names(cats(mu)[[1]])
# convert categories associated with keys to values
mu2 <- catalyze(mu)Inspect just the plant available water 0-25cm.
Plot aggregate soil properties.
Here is an example of not so great exact join between soil survey areas. In this case the one soil survey was published in 1979 and the other in 2004.
First, we setup BBOX and query map unit key WCS.
# extract a BBOX like this from SoilWeb by pressing "b"
bb <- '-91.6853 36.4617,-91.6853 36.5281,-91.5475 36.5281,-91.5475 36.4617,-91.6853 36.4617'
wkt <- sprintf('POLYGON((%s))', bb)
# init sf object from WKT
x <- st_as_sfc(wkt, crs = 4326)
# get gSSURGO grid here
mu <- mukey.wcs(aoi = x, db = 'gssurgo')
# note SSA boundary
plot(mu, legend = FALSE, axes = FALSE)Then we derive aggregate sand, silt, clay (RV) values from the largest component, taking the weighted mean over 25-50cm depth interval. We also will take the sand and clay values to calculate the surface texture class for comparison.
# extract RAT for thematic mapping
rat <- cats(mu)[[1]]
# variables of interest
vars <- c("sandtotal_r", "silttotal_r", "claytotal_r")
# get thematic data from SDA
# dominant component
# depth-weighted average
# sand, silt, clay (RV)
tab <- get_SDA_property(property = vars,
method = "Dominant Component (Numeric)",
mukeys = as.integer(rat$mukey),
top_depth = 25,
bottom_depth = 50)
# check
head(tab)
# set raster categories
levels(mu) <- tab[, c('mukey', vars)]
# convert mukey grid + RAT -> stack of numerical grids
# retaining only sand, silt, clay via [[vars]]
ssc <- catalyze(mu)
# create a copy of the grid
texture.class <- ssc[[1]]
names(texture.class) <- 'soil.texture'
# assign soil texture classes for the fine earth fraction
# using sand and clay percentages
values(texture.class) <- aqp::ssc_to_texcl(
sand = values(ssc[['sandtotal_r']]),
clay = values(ssc[['claytotal_r']]),
droplevels = FALSE
)
r <- c(ssc, texture.class)
# graphical check
plot(
r,
cex.main = 0.7,
main = paste0(names(r), " - 25-50cm\nDominant Component")
)