Map Unit Key Web Coverage Service
2024-05-02
Introduction
A web coverage service (WCS) is provided for the gSSURGO and gNATGSO map unit key grids. These grids represent a rasterizaion of the gSSURGO and gNATSGO map unit keys for the conterminous United States at a resolution of 30m, referenced to an Albers Equal Area Conic (NAD83) coordinate reference system (EPSG:5070). Grids are LZW compressed and internally tiled for efficient random access. Cell values are map unit keys, encoded as unsigned 32-bit integers. The grid topology and cell values are identical to the rasters contained within the gSSURGO and gNATSGO file geo-databases (FGDB). The GeoTiff format is used to ensure maximum compatibility. Files are re-created as part of the annual SSURGO refresh cycle.
An experimental WCS is provided for raster soil survey (RSS) products at 10m resolution.
An experimental WCS is provided for STATSGO2 (2016) at 300m resolution.
An experimental WCS is provided for SSURGO in HI (EPSG:6628) and PR (EPSG:32161), at 30m resolution.
A similar WCS is provided for the ISSR-800 soil property grids.
Grid Selection
See the mukey.wcs
and ISSR800.wcs
manual pages for details.
Basic usage.
# 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', ...)
# select various ISSR-800 grids, details below
x <- ISSR800.wcs(aoi = aoi, var = 'paws')gSSURGO
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.
gNATSGO
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.
Raster Soil Survey
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.
STATSGO
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
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.
Caveats / Limitations
- Requests are limited to images sizes of 5000x5000 pixels, this is approximately 1x1 degrees at 30m.
- These instructions are based on the development version of
soilDB, which has recently switched over tosfandterrapackages.
Manually Creating bounding-boxes for WCS requests
The WCS interface functions mukey.wcs and
ISSR800.wcs can automatically generate a bounding-box
(BBOX) from most spatial data formats. However, sometimes it is more
convenient to manually specify a BBOX created from a website or single
point specified in WGS84 coordinates. SoilWeb
provides two convenient keyboard shortcuts, available after clicking on
the map interface:
‘b’ copy bounding-box, returned as: -118.6609 36.4820,-118.6609 36.5972,-118.3979 36.5972,-118.3979 36.4820,-118.6609 36.4820
‘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.
A SoilWeb-style BBOX can be directly converted into an
sf object:
library(sf)
library(soilDB)
library(terra)
# paste your BBOX text here
bb <- '-118.6609 36.4820,-118.6609 36.5972,-118.3979 36.5972,-118.3979 36.4820,-118.6609 36.4820'
# convert text -> WKT -> sf
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')# looks OK
plot(mu, legend = FALSE, axes = FALSE, main = metags(mu, name = 'description'))
# add original BBOX, after transforming to mukey grid CRS
plot(st_transform(x, 5070), add = TRUE)
Metadata
Print attached metadata with theterra::metags() function.
Currently, a short description and “vintage” are included. The gNATSGO,
gSSURGO, and derivatives are typically created 2-3 months after the
annual SSURGO refresh on October 1st of each year (beginning of US
federal fiscal year).
| description | vintage |
|---|---|
| gSSURGO map unit keys | FY2024 |
A buffer applied to a single WGS84 coordinate can create a BBOX:
library(sf)
library(soilDB)
# example point, WGS84 coordinates
p <- st_as_sf(
data.frame(
x = -118.55639,
y = 36.52578
),
coords = c('x', 'y'),
crs = 4326
)
# 1000m buffer applied to WGS84 coordinate
# radius defined in meters
b <- st_buffer(p, dist = units::set_units(1000, m))
b <- st_as_sf(b)
# query WCS
# result is in EPSG:5070
mu <- mukey.wcs(b, db = 'gSSURGO')# looks OK
plot(mu, legend = FALSE, axes = FALSE, main = metags(mu, name = 'description'))
# add buffer, after transforming to mukey grid CRS
plot(st_transform(b, 5070), add = TRUE)
# add original point, after transforming to mukey grid CRS
plot(st_transform(p, 5070), add = TRUE, pch = 16)
Setup
Need the latest version of soilDB and
terra, along with all packages specified below (CRAN
versions fine).
install.packages('soilDB')
remotes::install_github("ncss-tech/soilDB", dependencies = FALSE, upgrade = FALSE, build = FALSE)
# terra
install.packages('terra', repos = 'https://rspatial.r-universe.dev')Load required packages.
# latest from GitHub
# remotes::install_github("ncss-tech/soilDB", dependencies = FALSE)
library(soilDB)
# wrangling polygons and CRS transformations
library(sf)
# raster data / analysis
# need latest: install.packages('terra', repos='https://rspatial.r-universe.dev')
library(terra)
# raster visualization
library(rasterVis)
# figures
library(lattice)
library(tactile)
# soil classification
library(aqp)
# color palettes and manipulation
library(RColorBrewer)
library(colorspace)Map Unit Key Grids
Use the mukey.wcs() function to access chunks of the
CONUS map unit key grids, based on some representation of an area of
interest (AOI). The AOI can be defined manually, as below, or
automatically extracted from sf, sfc,
bbox, SpatRaster, or SpatVector
objects. Backwards-compatibility with Spatial* and
RasterLayer objects is also provided. The resulting grid of
integers isn’t all that useful by itself; join data from Soil Data
Access (SDA) or local files (by map unit key) to create thematic maps.
Note: the gSSURGO and gNATSGO grids are updated annually, and typically
in-sync with the live version of the data hosted by SDA (updated Oct
1st) by early November. Map unit keys can change over time, especially
in soil survey areas that were updated during the last fiscal year.
# 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)## class : SpatRaster
## dimensions : 147, 219, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : -1365495, -1358925, 2869245, 2873655 (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83 / Conus Albers (EPSG:5070)
## source(s) : memory
## varname : file3b28190611e1
## categories : mukey
## name : mukey
## min value : 144983
## max value : 1716001# looks good
plot(
mu,
main = 'gSSURGO map unit keys',
sub = 'Albers Equal Area Projection',
axes = FALSE,
legend = FALSE,
col = hcl.colors(100)
)
SSURGO Polygons from SDA
It is possible to retrieve small areas of vector geometry (SSURGO
polygons) from SDA with the SDA_spatialQuery() function.
These data are stored and delivered in a geographic coordinate system
(WGS84). Overlaying the SSURGO polygons and map unit key grids will
require a simple transformation.
# get intersecting SSURGO polygons from SDA
# result is a SpatVector object
# GCS WGS84
p <- SDA_spatialQuery(mu, what = 'mupolygon', geomIntersection = TRUE)
# transform to AEA coordinate reference system used by gSSURGO / gNATSGO
# this is EPSG:5070
p <- project(p, crs(mu))
par(mar = c(1, 0, 2, 0))
plot(p, main = 'SSURGO Polygons (SDA)', axes = FALSE)
mtext('Albers Equal Area Projection', side = 1, line = -0.5)
Overlay SSURGO polygons and 30m map unit key grid.
par(mar = c(1, 0, 2, 0))
plot(mu, main = 'gSSURGO Grid (WCS)\nSSURGO Polygons (SDA)', col = hcl.colors(50), axes = FALSE, legend = FALSE)
plot(p, add = TRUE, border = 'white')
mtext('CONUS Albers Equal Area Projection (EPSG:5070)', side = 1, line = 1)
Grid Resolution Specification
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 much larger chunk of data by requesting grids at 800m. However, pixels are selected by nearest-neighbor and not generalized to the coarser scale. In short: don’t use the resulting data for anything other than a simple preview.
# part of CA
# make a bounding box and assign a CRS (4326: GCS, WGS84)
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)# OK
plot(
x.800,
main = 'A Preview of gSSURGO Map Unit Keys',
sub = 'Albers Equal Area Projection (800m)\nnearest-neighbor resampling',
axes = FALSE,
legend = FALSE,
col = hcl.colors(100),
)
Raster Soil Survey Data
Note: tabular data for RSS are not yet available via SDA.
# slightly larger than Coweeta RSS BBOX
# specified in EPSG:5070
a <- st_bbox(
c(xmin = 1129000, xmax = 1135000, ymin = 1403000, ymax = 1411000),
crs = st_crs(5070)
)
# convert to sf polygon for later use in figures
a <- st_as_sfc(a)
# gSSURGO grid: 30m resolution
(x <- mukey.wcs(a, db = 'gSSURGO', res = 30))## class : SpatRaster
## dimensions : 267, 200, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 1129005, 1135005, 1402995, 1411005 (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83 / Conus Albers (EPSG:5070)
## source(s) : memory
## varname : file3b2855f94464
## categories : mukey
## name : mukey
## min value : 545800
## max value : 545887# gNATSGO grid: 30m resolution
(y <- mukey.wcs(a, db = 'gNATSGO', res = 30))## class : SpatRaster
## dimensions : 267, 200, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 1129005, 1135005, 1402995, 1411005 (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83 / Conus Albers (EPSG:5070)
## source(s) : memory
## varname : file3b286d6c75a2
## categories : mukey
## name : mukey
## min value : 545800
## max value : 3244759# RSS grid: 10m resolution
(z <- mukey.wcs(a, db = 'RSS', res = 10))## class : SpatRaster
## dimensions : 800, 600, 1 (nrow, ncol, nlyr)
## resolution : 10, 10.0125 (x, y)
## extent : 1129005, 1135005, 1402995, 1411005 (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83 / Conus Albers (EPSG:5070)
## source(s) : memory
## varname : file3b285cf22aa0
## categories : mukey
## name : mukey
## min value : 3244721
## max value : 3244759# graphical comparison
par(mfcol = c(1, 3))
# gSSURGO
plot(x, axes = FALSE, legend = FALSE, main = metags(x, name = 'description'))
# original BBOX
plot(a, add = TRUE)
# gNATSGO
plot(y, axes = FALSE, legend = FALSE, main = metags(y, name = 'description'))
# original BBOX
plot(a, add = TRUE)
# RSS
plot(z, axes = FALSE, legend = FALSE, main = metags(z, name = 'description'), ext = x)
# original BBOX
plot(a, add = TRUE)
STATSGO
Continuing from the example above.
(statsgo <- mukey.wcs(a, db = 'statsgo', res = 300))## class : SpatRaster
## dimensions : 27, 20, 1 (nrow, ncol, nlyr)
## resolution : 300, 300 (x, y)
## extent : 1129005, 1135005, 1402995, 1411095 (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83 / Conus Albers (EPSG:5070)
## source(s) : memory
## varname : file3b283903451
## categories : mukey
## name : mukey
## min value : 659074
## max value : 664845# graphical comparison
par(mfcol = c(1, 2))
# gSSURGO
plot(x, axes = FALSE, legend = FALSE, main = metags(x, name = 'description'))
# STATSGO
plot(statsgo, axes = FALSE, legend = FALSE, main = metags(statsgo, name = 'description'))
Hawaii SSURGO
An experimental, 30m SSURGO mukey WCS based on the EPSG:6628 coordinate reference system. 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 text -> WKT -> sf
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 = 'hi_ssurgo')# make NA (the ocean) blue
plot(mu, legend = FALSE, axes = FALSE, main = metags(mu, name = 'description'), colNA = 'royalblue')
# # check mu names
# .is <- format_SQL_in_statement(cats(mu)[[1]]$mukey)
# .sql <- sprintf("SELECT mukey, muname FROM mapunit WHERE mukey IN %s", .is)
# knitr::kable(SDA_query(.sql))Puerto Rico SSURGO
An experimental, 30m SSURGO mukey WCS based on the EPSG:32161 coordinate reference system. 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 text -> WKT -> sf
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 = 'pr_ssurgo')# make NA (the ocean) blue
plot(mu, legend = FALSE, axes = FALSE, main = metags(mu, name = 'description'), colNA = 'royalblue')
# # check mu names
# .is <- format_SQL_in_statement(cats(mu)[[1]]$mukey)
# .sql <- sprintf("SELECT mukey, muname FROM mapunit WHERE mukey IN %s", .is)
# knitr::kable(SDA_query(.sql))Thematic Mapping
The following example BBOX + resulting gSSURGO mukey grid will be used for the following 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')Map Unit Aggregate Values
The Mapunit Aggregated 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() convenience function, or
write a simple query in SQL and submit via SDA_query().
# copy example grid
mu2 <- mu
# extract RAT for thematic mapping
rat <- cats(mu2)[[1]]
# optionally use convenience function:
# * returns all fields from muagatt table
# * along with map unit name
# tab <- get_SDA_muaggatt(mukeys = as.numeric(rat$mukey), query_string = TRUE)
# the SQL is simple, only need a few columns from a single table
# safely create IN statement from vector of mukey
.sql <- sprintf(
"SELECT mukey, aws050wta, aws0100wta FROM muaggatt WHERE mukey IN %s",
format_SQL_in_statement(as.numeric(rat$mukey))
)
# run query, result is a data.frame
tab <- SDA_query(.sql)
# check
head(tab)## mukey aws050wta aws0100wta
## 1 144983 7.26 14.02
## 2 144984 7.31 14.14
## 3 144985 7.36 14.27
## 4 144986 7.06 13.19
## 5 145005 5.21 9.34
## 6 145009 5.49 9.16# merge aggregate soil data into RAT
rat <- merge(rat, tab, by.x = 'mukey', by.y = 'mukey', sort = FALSE, all.x = TRUE)
# requires that grid cell ID (mukey) be numeric
rat$mukey <- as.integer(rat$mukey)
levels(mu2) <- rat
# convert grid + RAT -> stack of propery grids
# note explicit subset of returned grids, we don't need the mukey grid
aws <- catalyze(mu2)[[c('aws050wta', 'aws0100wta')]]
# simple plot, note scales are different
plot(aws, axes = FALSE, col = hcl.colors(10, 'mako'), main = c('Plant Available Water Storage (cm)\nWeighted Mean over Components, 0-50cm', 'Plant Available Water Storage (cm)\nWeighted Mean over Components, 0-100cm'))
Interpretations for Soil Suitability / Limitation
Use the get_SDA_interpretation() convenience
function.
# 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)## areasymbol musym muname mukey
## 1 MT629 101 McCollum fine sandy loam, 0 to 2 percent slopes 144983
## 2 MT629 102 McCollum fine sandy loam, 2 to 4 percent slopes 144984
## 3 MT629 103 McCollum fine sandy loam, 4 to 8 percent slopes 144985
## 4 MT629 104 McCollum fine sandy loam, gravelly substratum, 0 to 2 percent slopes 144986
## 5 MT629 120 Niarada gravelly loam, 0 to 4 percent slopes 145005
## 6 MT629 123 Niarada gravelly loam, cool, 15 to 30 percent slopes 145009
## rating_ENGConstructionMaterialsRoadfill rating_AWMIrrigationDisposalofWastewater
## 1 1.00 0.00
## 2 0.98 0.05
## 3 0.98 0.68
## 4 1.00 0.15
## 5 0.98 0.19
## 6 0.12 1.00
## class_ENGConstructionMaterialsRoadfill class_AWMIrrigationDisposalofWastewater
## 1 Well suited Not limited
## 2 Moderately well suited Slightly limited
## 3 Moderately well suited Moderately limited
## 4 Well suited Slightly limited
## 5 Moderately well suited Slightly limited
## 6 Poorly suited Very limited
## reason_ENGConstructionMaterialsRoadfill
## 1 <NA>
## 2 Low strength (0.778); Dusty (0.785)
## 3 Low strength (0.778); Dusty (0.785)
## 4 <NA>
## 5 Dusty (0.981); Dusty (0.981); Dusty (0.902)
## 6 Slope (0.08); Dusty (0.975); Slope (0.08); Dusty (0.981); Depth to bedrock (0); Slope (0)
## reason_AWMIrrigationDisposalofWastewater
## 1 <NA>
## 2 Filtering capacity (1); Too steep for surface application (0.08)
## 3 Too steep for surface application (0.68); Too steep for surface application (0.68)
## 4 Filtering capacity (1); Filtering capacity (1); Droughty (0.918)
## 5 Droughty (0.133); Large stones on the surface (0.141); Droughty (0.117); Filtering capacity (1); Droughty (0.918); Slow water movement (0.372)
## 6 Too steep for surface application (1); Too steep for sprinkler application (1); Droughty (0.061); Too steep for surface application (1); Too steep for sprinkler application (1); Large stones on the surface (0.141); Droughty (0.117); Droughty (1); Too steep for surface application (1); Too steep for sprinkler application (1); Depth to bedrock (1); Large stones on the surface (0.684); Filtering capacity (1); Too steep for surface application (1); Droughty (0.918); Too steep for sprinkler application (0.395)# set factor levels
tab$class_ENGConstructionMaterialsRoadfill <- factor(
tab$class_ENGConstructionMaterialsRoadfill,
levels = c('Not suited', 'Poorly suited', 'Moderately suited', 'Moderately well suited', 'Well suited', 'Not Rated'),
ordered = TRUE
)
bwplot(
class_ENGConstructionMaterialsRoadfill ~ rating_ENGConstructionMaterialsRoadfill,
data = tab,
xlab = 'Fuzzy Rating',
main = 'ENG - Construction Materials; Roadfill',
par.settings = tactile.theme()
)
# merge aggregate soil data into RAT
rat <- merge(rat, tab, by.x = 'mukey', by.y = 'mukey', sort = FALSE, all.x = TRUE)
# requires that grid cell ID (mukey) be numeric
rat$mukey <- as.integer(rat$mukey)
levels(mu2) <- rat
# convert grid + RAT -> stack of property grids
# note explicit subset of returned grids, we don't need the mukey grid
vars <- c('rating_ENGConstructionMaterialsRoadfill', 'rating_AWMIrrigationDisposalofWastewater')
rating <- catalyze(mu2)[[vars]]
# simple plot, note scales are different
plot(rating, axes = FALSE, col = hcl.colors(10, 'mako'), main = c('Construction Materials; Roadfill\nWeighted Mean over Components', 'Irrigation Disposal of Wastewater\nWeighted Mean over Components'))
Component level properties
Steel corrosion potential.
# 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)
)
# convert to factor and set levels in logical order
tab$corsteel <- factor(tab$corsteel, levels = c('Low', 'Moderate', 'High'))
# merge aggregate soil data into RAT
rat <- merge(rat, tab, by.x = 'mukey', by.y = 'mukey', sort = FALSE, all.x = TRUE)
# requires that grid cell ID (mukey) be numeric
rat$mukey <- as.integer(rat$mukey)
levels(mu2) <- rat
# set active category to corsteel
activeCat(mu2) <- 'corsteel'
# plot
plot(mu2, axes = FALSE, col = hcl.colors(3), mar = c(0, 0, 0, 6))
Simplified component parent material group.
# 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))
# merge aggregate soil data into RAT
rat <- merge(rat, tab, by.x = 'mukey', by.y = 'mukey', sort = FALSE, all.x = TRUE)
# requires that grid cell ID (mukey) be numeric
rat$mukey <- as.integer(rat$mukey)
levels(mu2) <- rat
# set active category to corsteel
activeCat(mu2) <- 'pmgroupname'
plot(mu2, axes = FALSE, col = hcl.colors(10), mar = c(0, 0, 0, 12))
Hydric Rating.
# 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))
# merge aggregate soil data into RAT
rat <- merge(rat, tab, by.x = 'mukey', by.y = 'mukey', sort = FALSE, all.x = TRUE)
# requires that grid cell ID (mukey) be numeric
rat$mukey <- as.integer(rat$mukey)
levels(mu2) <- rat
# set active category to corsteel
activeCat(mu2) <- 'HYDRIC_RATING'
plot(mu2, axes = FALSE, col = hcl.colors(10), mar = c(0, 0, 0, 6))
Several horizon-level soil properties
The get_SDA_property()
function from soilDB is convenient interface to aggregated
SSURGO/STATSGO tabular data via SDA.
Setup 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')# check: OK
plot(
mu,
legend = FALSE,
axes = FALSE,
col = hcl.colors(100),
main = 'gSSURGO Map Unit Key Grid'
)
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
p <- get_SDA_property(property = vars,
method = "Dominant Component (Numeric)",
mukeys = as.integer(rat$mukey),
top_depth = 0,
bottom_depth = 25)
# check: OK
head(p)## mukey areasymbol musym muname
## 1 623396 OK113 1 Apperson silty clay loam, 1 to 3 percent slopes
## 2 623399 OK113 4 Coyle loam, 1 to 3 percent slopes
## 3 623402 OK113 7 Keokuk very fine sandy loam, 0 to 1 percent slopes, occasionally flooded
## 4 623405 OK113 10 Bethany silt loam, 1 to 3 percent slopes
## 5 623406 OK113 11 Bethany silt loam, 3 to 5 percent slopes
## 6 623407 OK113 12 Bethany-Pawhuska complex, 1 to 5 percent slopes
## dbthirdbar_r awc_r ph1to1h2o_r
## 1 1.45 0.18 6.15
## 2 1.40 0.18 5.70
## 3 1.43 0.17 7.30
## 4 1.35 0.20 6.20
## 5 1.30 0.20 6.30
## 6 1.34 0.20 6.30# convert areasymbol into a factor easy plotting later
p$areasymbol <- factor(p$areasymbol)
# merge aggregate soil data into RAT
rat <- merge(rat, p, by.x = 'mukey', by.y = 'mukey', sort = FALSE, all.x = TRUE)
# requires that grid cell ID (mukey) be numeric
rat$mukey <- as.integer(rat$mukey)
levels(mu) <- rat
# result is a grid of map unit keys + RATPlant available water 0-25cm.
# list variables in the RAT
names(cats(mu)[[1]])## [1] "mukey" "ID" "areasymbol" "musym" "muname" "dbthirdbar_r"
## [7] "awc_r" "ph1to1h2o_r"activeCat(mu) <- 'awc_r'
plot(mu)
Plot aggregated soil properties.
# reset the active category to 1
# this will ensure that all variables in the RAT will be converted by catalyze(mu)
activeCat(mu) <- 1
# convert mukey grid + RAT -> stack of numerical grids
mu.stack <- catalyze(mu)
# keep only properties / remove IDs
mu.stack <- mu.stack[[vars]]
# note implicit simplification via maxpixels
levelplot(
mu.stack[['dbthirdbar_r']],
main = '1/3 Bar Bulk Density (g/cm^3)\nDominant Component\n0-25cm',
margin = FALSE,
scales = list(draw = FALSE),
col.regions = hcl.colors,
maxpixels = 1e5
)
levelplot(
mu.stack[['awc_r']],
main = 'AWC (cm/cm)\nDominant Component\n0-25cm',
margin = FALSE,
scales = list(draw = FALSE),
col.regions = hcl.colors,
maxpixels = 1e5
)
levelplot(
mu.stack[['ph1to1h2o_r']],
main = 'pH 1:1 H2O\nDominant Component\n0-25cm',
margin = FALSE,
scales = list(draw = FALSE),
col.regions = hcl.colors,
maxpixels = 1e5
)
Sand, Silt, and Clay at a Soil Survey Area Boundary
A not so great SSA join (1979 vs 2004).
Setup BBOX and query mukey 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)
# assign GCS WGS84 coordinate reference system
st_crs(x) <- 4326
# get gSSURGO grid here
mu <- mukey.wcs(aoi = x, db = 'gssurgo')# note SSA boundary
plot(mu, legend = FALSE, axes = FALSE, col = hcl.colors(100))
Derive aggregate sand, silt, clay (RV) values from the largest component, taking the weighted mean over 25-50cm depth interval.
# 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)
p <- get_SDA_property(property = vars,
method = "Dominant Component (Numeric)",
mukeys = as.integer(rat$mukey),
top_depth = 25,
bottom_depth = 50)
# check
head(p)## mukey areasymbol musym muname
## 1 691980 MO091 73306 Gressy-Gatewood complex, 3 to 8 percent slopes, rocky
## 2 2502332 MO149 73321 Alred-Gatewood complex, 1 to 8 percent slopes
## 3 2502334 MO149 73322 Alred-Gatewood complex, 8 to 15 percent slopes
## 4 2503322 MO149 76002 Batcave-Farewell complex, 1 to 3 percent slopes, frequently flooded
## 5 2503473 MO149 76046 Secesh silt loam, 1 to 3 percent slopes, rarely flooded
## 6 2503476 MO149 76047 Secesh-Tilk complex, 1 to 3 percent slopes, occasionally flooded
## sandtotal_r silttotal_r claytotal_r
## 1 23.00 59.00 18.00
## 2 28.69 51.11 20.20
## 3 28.69 51.11 20.20
## 4 40.00 40.00 20.00
## 5 25.65 49.00 25.35
## 6 26.91 48.08 25.02# merge aggregate soil data into RAT
rat <- merge(rat, p, by.x = 'mukey', by.y = 'mukey', sort = FALSE, all.x = TRUE)
# requires that grid cell ID (mukey) be numeric
rat$mukey <- as.integer(rat$mukey)
levels(mu) <- rat
# convert mukey grid + RAT -> stack of numerical grids
# retaining only sand, silt, clay via [[vars]]
ssc <- catalyze(mu)[[vars]]
# graphical check
# note implicit simplification via maxpixels
levelplot(
ssc,
main = 'Sand, Silt, Clay (RV)\nDominant Component\n25-50cm',
margin = FALSE,
scales = list(draw = FALSE),
col.regions = hcl.colors,
maxpixels = 1e4,
layout = c(3, 1)
)
Convert to soil texture class of <2mm fraction.
# create a new grid by copy
# we will place soil texture classes into this 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) <- ssc_to_texcl(
sand = values(ssc[['sandtotal_r']]),
clay = values(ssc[['claytotal_r']]),
droplevels = FALSE
)
plot(texture.class, col = hcl.colors(50), axes = FALSE, main = 'Soil Texture Class <2mm Fraction\nDominant Component, 25-50cm (RV)')
TODO SOC Stock Estimates + Expected Variation
Other Relevant WCS
ISSR-800
Available data sets. (TODO: there are more more to add, including NCCPI)
WCS_details(wcs = 'ISSR800')| var | description |
|---|---|
| caco3_kg_sq_m | Total CaCO3 (kg/m^2) |
| cec_025cm | CEC at pH 7 0-25cm depth (cmol[+]/kg) |
| cec_050cm | CEC at pH 7 0-50cm depth (cmol[+]/kg) |
| cec_05cm | CEC at pH 7 0-5cm depth (cmol[+]/kg) |
| clay_025cm | clay percent 0-25cm depth |
| clay_05cm | clay percent 0-5cm depth |
| clay_2550cm | clay percent 25-50cm depth |
| clay_3060cm | clay percent 30-60cm depth |
| drainage_class | Soil Drainage Class |
| ec_025cm | EC 0-25cm depth (dS/m) |
| ec_05cm | EC 0-5cm depth (dS/m) |
| greatgroup | Soil Taxonomy: Greatgroup |
| hydgrp | Hydrologic Soil Group |
| lcc_irrigated | Land Capability Class, irrigated |
| lcc_nonirrigated | Land Capability Class, non-irrigated |
| om_kg_sq_m | Total Soil Organic Matter (kg/m^2) |
| paws | total plant available water storage (cm water) |
| paws_025cm | plant available water storage 0-25cm depth (cm water) |
| paws_050cm | plant available water storage 0-50cm depth (cm water) |
| ph_025cm | pH 1:1 H2O 0-25cm depth |
| ph_05cm | pH 1:1 H2O 0-5cm depth |
| ph_2550cm | pH 1:1 H2O 25-50cm depth |
| ph_3060cm | pH 1:1 H2O 30-60cm depth |
| sand_025cm | sand percent 0-25cm depth |
| sand_05cm | sand percent 0-5cm depth |
| sand_2550cm | sand percent 25-50cm depth |
| sand_3060cm | sand percent 30-60cm depth |
| sar | SAR, entire profile |
| series_name | Soil Series Name |
| silt_025cm | silt percent 0-25cm depth |
| silt_05cm | silt percent 0-5cm depth |
| silt_2550cm | silt percent 25-50cm depth |
| silt_3060cm | silt percent 30-60cm depth |
| soilorder | Soil Taxonomy: Soil Order |
| ssurgo_pct | SSURGO data available, fraction of 800x800m grid cell |
| statsgo_pct | STATSGO data available, fraction of 800x800m grid cell |
| str | Soil Temperature Regime |
| suborder | Soil Taxonomy: Suborder |
| texture_025cm | Soil Texture Class, 0-25cm |
| texture_05cm | Soil Texture Class, 0-5cm |
| texture_2550cm | Soil Texture Class, 25-50cm |
| weg | Wind Erodibility Group |
| wei | Wind Erodibility Index |
Examples for a broad slice of California
# make a bounding box and assign a CRS (4326: GCS, WGS84)
a.CA <- st_bbox(
c(xmin = -121, xmax = -120, ymin = 37, ymax = 38),
crs = st_crs(4326)
)
# convert bbox to sf geometry
a.CA <- st_as_sfc(a.CA)
# floating point grids
paws_total <- ISSR800.wcs(aoi = a.CA, var = 'paws')
paws_025cm <- ISSR800.wcs(aoi = a.CA, var = 'paws_025cm')
paws_050cm <- ISSR800.wcs(aoi = a.CA, var = 'paws_050cm')
pH_05cm <- ISSR800.wcs(aoi = a.CA, var = 'ph_05cm')
pH_3060cm <- ISSR800.wcs(aoi = a.CA, var = 'ph_3060cm')
clay_05cm <- ISSR800.wcs(aoi = a.CA, var = 'clay_05cm')
clay_3060cm <- ISSR800.wcs(aoi = a.CA, var = 'clay_3060cm')
silt_3060cm <- ISSR800.wcs(aoi = a.CA, var = 'silt_3060cm')
sand_3060cm <- ISSR800.wcs(aoi = a.CA, var = 'sand_3060cm')
## NA encoding issues
# 0 is a real value
sar <- ISSR800.wcs(aoi = a.CA, var = 'sar')
# 16bit integer grids
wei <- ISSR800.wcs(aoi = a.CA, var = 'wei')
# 8bit unsigned (BYTE) grids with RAT
drainage_class <- ISSR800.wcs(aoi = a.CA, var = 'drainage_class')
weg <- ISSR800.wcs(aoi = a.CA, var = 'weg')
str <- ISSR800.wcs(aoi = a.CA, var = 'str')
# metadata
metags(sar)## description vintage
## "SAR, entire profile" "FY2024"Check to see if aggregated sand + silt + clay is roughly = 100%. Almost. This isn’t surprising since each component of soil texture is aggregated separately.
z <- sand_3060cm + silt_3060cm + clay_3060cm
plot(z - 100)
Soil texture class 25-50cm.
texture_2550cm <- ISSR800.wcs(aoi = a.CA, var = 'texture_2550cm')
# plot using SoilWeb soil texture class colors
# preset by ISSR800.wcs()
plot(texture_2550cm, axes = FALSE)
Other soil properties of interest.
plot(sar, axes = FALSE, col = hcl.colors(50), main = 'SAR')
plot(wei, axes = FALSE, col = hcl.colors(50), main = 'Wind Erodibility Index')
plot(drainage_class, axes = FALSE, col = hcl.colors(50), main = 'Drainage Class')
plot(weg, axes = FALSE, col = hcl.colors(50), main = 'Wind Erodibility Group')
plot(str, axes = FALSE, col = hcl.colors(50), main = 'Soil Temperature Regime')
Stacking soil property grids for multi-panel figures that share a common color ramp / legend.
# combine into multi-layer SpatRaster
# like raster::stack()
s <- c(pH_05cm, pH_3060cm)
# plot multiple layers via SpatRaster plot method
# common color map / legend
plot(s, col = hcl.colors(50), axes = FALSE, mar = c(2, 2, 2, 4), range = c(4.5, 9))
I like this approach.
# or, plot using rasterVis::levelplot
# note laundering of layer names, required when original layer names contain spaces, punctuation, etc.
names(s) <- make.names(names(s))
# manually specify panel names
# in original ordering specified to c()
# formatting ideas: https://stackoverflow.com/questions/34652580/change-raster-panel-titles-using-levelplot
levelplot(s,
scales = list(draw = FALSE),
col.regions = hcl.colors,
names.attr = c('pH 1:1 H20, 0-5cm', 'pH 1:1 H20, 30-60cm')
)
Another example, with 3 panels.
s <- c(sand_3060cm, silt_3060cm, clay_3060cm)
names(s) <- make.names(names(s))
levelplot(s,
scales = list(draw = FALSE),
col.regions = hcl.colors,
names.attr = c('Percent Sand, 30-60cm', 'Percent Silt, 30-60cm', 'Percent Clay, 30-60cm')
)
Soil / Land Classification
# use terra methods to define a BBOX, from SoilWeb 'b' keypress
# https://casoilresource.lawr.ucdavis.edu/gmap/?loc=47.24008,-99.73114,z11
# generate WKT description of viewport BBOX
bb <- '-100.2976 47.0407,-100.2976 47.4717,-99.3013 47.4717,-99.3013 47.0407,-100.2976 47.0407'
wkt <- sprintf('POLYGON((%s))', bb)
# init spatVector
a <- vect(wkt, crs = 'epsg:4326')
# results are grids + RAT
greatgroup <- ISSR800.wcs(aoi = a, var = 'greatgroup')
# create a large color palette for greatgroup classes
cols <- brewer.pal(12, 'Set3')
cols <- darken(cols, amount = 0.25)
cr <- colorRampPalette(cols, space = 'Lab', interpolate = 'spline')
plot(greatgroup, axes = FALSE, col = cr(50), mar = c(1, 1, 1, 6))
# both LCC data may or may not available in some places
lcc_nonirrigated <- ISSR800.wcs(aoi = a, var = 'lcc_nonirrigated')
lcc_irrigated <- ISSR800.wcs(aoi = a, var = 'lcc_irrigated')
# check NA fraction
global(lcc_nonirrigated, 'isNA') / ncell(lcc_nonirrigated)## isNA
## label 0.000490918global(lcc_irrigated, 'isNA') / ncell(lcc_irrigated)## isNA
## Land Capability Class, irrigated 1# colors for irrigated/non-irrigated LCC
cols <- hcl.colors(50, 'Spectra', rev = TRUE)
cols <- colorspace::darken(cols, amount = 0.1)
cr <- colorRampPalette(cols, space = 'Lab', interpolate = 'spline')
plot(lcc_nonirrigated, axes = FALSE, col = cr(50), mar = c(1, 1, 1, 6), main = 'Non-Irrigated LCC')
Soil Color
# make a bounding box and assign a CRS (4326: GCS, WGS84)
a <- st_bbox(
c(xmin = -121, xmax = -120, ymin = 37, ymax = 38),
crs = st_crs(4326)
)
# convert bbox to sf geometry
a <- st_as_sfc(a)
# moist soil color at 25cm, low-res version
s <- soilColor.wcs(aoi = a, var = 'sc025cm', res = 270)
# metadata
metags(s)## description vintage
## "Moist soil color, 25cm (270m)" "FY2023"# color table is pre-set by soilColor.wcs()
plot(s, legend = FALSE, axes = FALSE, main = metags(s, name = 'description'))
Soil color distribution.
par(mar = c(4.5, 5, 1, 1))
# specific to this AOI
rat <- cats(s)[[1]]
barplot(s, col = rat$col, names.arg = rat$munsell, las = 1, cex.names = 0.75, horiz = TRUE, xlab = 'Cell Count')
# these colors, but frequencies tabulated across all of CONUS
# barplot(rat$freq, col = rat$col, las = 3, names.arg = rat$munsell, cex.names = 0.6)Detailed soil color mapping.
# from SoilWeb
bb <- '-90.1576 41.8579,-90.1576 41.9193,-90.0100 41.9193,-90.0100 41.8579,-90.1576 41.8579'
# use terra methods
# convert text -> WKT -> spatVect
wkt <- sprintf('POLYGON((%s))', bb)
a <- vect(wkt, crs = 'epsg:4326')
# moist soil color at 25cm, low-res version
s <- soilColor.wcs(aoi = a, var = 'sc025cm_hr', res = 30)
# color table is pre-set by soilColor.wcs()
plot(s, legend = FALSE, axes = FALSE, main = metags(s, name = 'description'))
Classification
# part of CA
# make a bounding box and assign a CRS (4326: GCS, WGS84)
a <- st_bbox(
c(xmin = -121, xmax = -120, ymin = 37, ymax = 38),
crs = st_crs(4326)
)
# convert bbox to sf geometry
a <- st_as_sfc(a)
# EC 0-25cm
ec <- ISSR800.wcs(aoi = a, var = 'ec_025cm')
# pH 0-25cm
pH <- ISSR800.wcs(aoi = a, var = 'ph_025cm')
# no ESP grid, use a constant value of 1%
esp <- ec
values(esp) <- 0
names(esp) <- 'ESP'
# copy EC grid, replace values with classification in next step
ss <- ec
names(ss) <- 'FAO Salt Severity'
# resulting raster contains factor levels (integers) and labels
values(ss) <- allocate(EC = values(ec), pH = values(pH), ESP = values(esp), to = 'FAO Salt Severity', droplevels = TRUE)
# note factor details stored in raster attribute table (RAT)
print(ss)## class : SpatRaster
## dimensions : 163, 142, 1 (nrow, ncol, nlyr)
## resolution : 800, 800 (x, y)
## extent : -2178400, -2064800, 1816800, 1947200 (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83 / Conus Albers (EPSG:5070)
## source(s) : memory
## varname : file3b28556b7014
## categories : label
## name : FAO Salt Severity
## min value : extremely saline
## max value : none# plot
plot(ss, axes = FALSE, main = 'FAO Salt Severity', col = hcl.colors(10))
FAO Salt Severity Classification
## TODO: NODATA not correctly handled, possibly need to mask
# combine individual SpatRaster objects into a multi-band SpatRaster
s <- c(pH, ec, esp)
# classification wrapper
.fun <- function(x) {
.pH <- x[1]
.ec <- x[2]
.esp <- x[3]
.res <- allocate(EC = .ec, pH = .pH, ESP = .esp, to = 'FAO Salt Severity', droplevels = FALSE)
return(.res)
}
# apply classification pixel-wise to stack
ss <- app(s, fun = .fun)
## TODO: can this be done within app()?
# manually specified levels
ll <- c("nonsaline", "slightly saline", "moderately saline",
"strongly saline", "very strongly saline", "extremely saline",
"saline-sodic", "slightly sodic", "moderately sodic", "strongly sodic",
"very strongly sodic")
# convert to grid + RAT
levels(ss) <- data.frame(ID = 1:length(ll), class = ll)
plot(ss, axes = FALSE, main = 'FAO Salt Severity', col = hcl.colors(10))PRISM Stack
Others
This document is based on soilDB version 2.8.3 and
terra version 1.7-75.