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.
See the related tutorial on using the CONUS mukey grids directly.
An experimental WCS is provided for raster soil survey (RSS) products at 10m resolution.
A similar WCS is provided for the ISSR-800 soil property grids.
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 various ISSR-800 grids, details below
x <- ISSR800.wcs(aoi = aoi, var = 'paws')
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.
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.
soilDB
,
which has recently
switched over to sf
and terra
packages.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
m <- mukey.wcs(x, db = 'gSSURGO')
# looks OK
plot(m, legend = FALSE, axes = FALSE, main = attr(m, 'layer name'))
box()
# add original BBOX, after transforming to mukey grid CRS
plot(st_transform(x, 5070), add = TRUE)
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
m <- mukey.wcs(b, db = 'gSSURGO')
# looks OK
plot(m, legend = FALSE, axes = FALSE, main = attr(m, 'layer name'))
box()
# 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)
Need the latest version of soilDB
, along with all
packages specified below (CRAN versions fine).
install.packages('soilDB')
remotes::install_github("ncss-tech/soilDB", dependencies = FALSE, upgrade = FALSE, build = FALSE)
Load required packages.
# latest from GitHub
library(soilDB)
# wrangling polygons and CRS transformations
library(sf)
# raster data / analysis
library(terra)
# raster visualization
library(rasterVis)
library(viridisLite)
# soil classification
library(aqp)
# color palettes and manipulation
library(RColorBrewer)
library(colorspace)
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)
)
# convert bbox to sf geometry
a <- st_as_sfc(a)
# fetch gSSURGO map unit keys at native resolution (30m)
x <- mukey.wcs(aoi = a, db = 'gssurgo')
# check:
print(x)
## class : SpatRaster
## dimensions : 147, 219, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : -1365484, -1358914, 2869259, 2873669 (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83 / Conus Albers (EPSG:5070)
## source : memory
## categories : mukey, ID
## name : mukey
## min value : 144983
## max value : 1716001
# looks good
levelplot(
x,
att = 'ID',
main = 'gSSURGO map unit keys',
sub = 'Albers Equal Area Projection',
margin = FALSE,
colorkey = FALSE,
col.regions = viridis
)
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.
# bounding box as an sf object
# generated from WCS result
# EPSG:5070
x.bbox <- st_as_sfc(st_bbox(x))
# get intersecting SSURGO linework from SDA
# result is an sf object in
# GCS WGS84
p <- SDA_spatialQuery(x.bbox, what = 'mupolygon', geomIntersection = TRUE)
# transform to AEA coordinate reference system used by gSSURGO / gNATSGO
# this is EPSG:5070
p <- st_transform(p, st_crs(x))
par(mar = c(1, 0, 2, 0))
plot(st_geometry(p), main = 'SSURGO Polygons (SDA)')
mtext('Albers Equal Area Projection', side = 1, line = -0.5)