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 geodatabases (FGDB). The GeoTiff format is used to ensure maximum compatibility. Files are re-created as part of the annual SSURGO refresh cycle.

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
x <- mukey.wcs(aoi = aoi, db = 'gssurgo', ...)

# select gNATSGO grid
x <- mukey.wcs(aoi = aoi, db = 'gnatsgo', ...)

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

Thematic Mapping

Thematic mapping or analysis of soil information requires connecting the grids to our tabular data sources, either using local files (e.g. gSSURGO / gNATSGO) or Soil Data Access (SDA). Note that map unit keys associated with raster soil surveys (included with gNATSGO) are not currently available in SDA.

Caveats / Limitations

Requests are limited to images sizes of 5000x5000 pixels, this is approximately 1x1 degrees at 30m.

Remaianing TODOs

  • Ensure that data are not resampled unless a resolution other than 30m is requested.
  • Currently the data come back close to 30m, usually 29.998-30.111. I’m not sure if this is something that has to be solved in the creation of the BBOX or if it is a Mapserver configuration issue. Still tinkering.
  • Consider converting source mukey grids to COG-GeoTiff, this will require a newer version of GDAL.
  • Simplify / unify soil data aggregation within SDA: this is much more generic than this WCS idea, but is a perfect use case.

A Simple Demonstration

Setup

Need the latest version of soilDB, along with all packages specified below (CRAN versions fine).

remotes::install_github("ncss-tech/soilDB", dependencies = FALSE, upgrade = FALSE, build = FALSE)
# latest from GitHub
library(soilDB)

# wrangling polygons and CRS transformations
library(sp)
library(rgdal)

# geometric calculations
library(rgeos)

# this is the future, will eventually replace much of 
# sp, rgdal, geos methods used here
library(sf)

# raster data visualization
library(raster)
library(rasterVis)
library(viridis)

# soil classification
library(aqp)
# 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')

# OK
levelplot(
  x, 
  att = 'ID',
  main = 'gSSURGO map unit keys',
  sub = 'Albers Equal Area Projection',
  margin = FALSE, 
  colorkey = FALSE, 
  col.regions = viridis
  )

Fetching Polygons from SDA

## TODO: update to sf

# convert raster extent into vector 
g <- as(extent(x), 'SpatialPolygons')
proj4string(g) <- projection(x)

# get intersecting SSURGO linework as SpatialPolygonsDataFrame from SDA
p <- SDA_spatialQuery(g, what = 'geom', geomIntersection = TRUE)

# transform to AEA CRS
p <- spTransform(p, CRS(projection(x)))

# compute area and convert square meters to acres
rgeos::gArea(p) * 0.000247105
## [1] 7140.472
par(mar = c(1, 0, 2, 0))
plot(p, main = 'SSURGO Polygons (SDA)')
mtext('Albers Equal Area Projection', side = 1, line = -0.5)

Overlay on gSSURGO grid.

levelplot(
  x, att = 'ID', 
  main = 'gSSURGO Grid (WCS)\nSSURGO Polygons (SDA)',
  margin = FALSE, 
  colorkey = FALSE,
  col.regions = viridis,
  panel = function(...) {
    panel.levelplot(...)
    sp.lines(p, col = 'white')
  }
)