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.

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.

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 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.

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

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
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)

Setup

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)

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)
)

# 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
  )

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.

# 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)