knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE, fig.align = 'center', fig.width = 8, fig.height = 6, tidy = FALSE, verbose = FALSE, # run when NASIS is defined or when R_SOILDB_SKIP_LONG_EXAMPLES is FALSE eval = isTRUE(try(local_NASIS_defined(), silent = TRUE)) || !as.logical(Sys.getenv("R_SOILDB_SKIP_LONG_EXAMPLES", unset = "TRUE")) ) options(width = 100, stringsAsFactors = FALSE)
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.
install.packages(c('soilDB', 'terra', 'sf'))
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.
# select various ISSR-800 grids, details below x <- ISSR800.wcs(aoi = aoi, var = 'paws')
This is where we will link to a more detailed ISSR800 vignette in the future.
Excerpt from the gSSURGO documentation.
Excerpt from the gNATSGO documentation.
Excerpt from the RSS documentation.
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.
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 = attr(mu, 'layer name')) box() # 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 = attr(mu, 'layer name')) box() # 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"
).
p <- project(p, crs(mu))
Inspect the result by overlaying SSURGO polygons on the 30m map unit key grid.
plot(mu, main = 'gSSURGO Grid (WCS)\nSSURGO Polygons (SDA)', axes = FALSE, legend = FALSE) plot(p, add = TRUE, border = 'white') mtext('CONUS Albers Equal Area Projection (EPSG:5070)', side = 1, line = 1)
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 = attr(x, 'layer name') ) plot(a, add = TRUE) # gNATSGO plot( y, axes = FALSE, legend = FALSE, main = attr(y, 'layer name') ) plot(a, add = TRUE) # RSS plot( z, axes = FALSE, legend = FALSE, main = attr(z, 'layer name'), 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.
(statsgo <- mukey.wcs(a, db = 'statsgo', res = 300)) # graphical comparison par(mfcol = c(1, 2)) # gSSURGO plot( x, axes = FALSE, legend = FALSE, main = attr(x, 'layer name') ) # STATSGO plot( statsgo, axes = FALSE, legend = FALSE, main = attr(statsgo, 'layer name') )
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 = attr(mu, 'layer name'), 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))
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 = attr(mu, 'layer name'), 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))
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, 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(mu2$awc_r) box()
Plot aggregate soil properties.
plot(mu2[['dbthirdbar_r']], cex.main = 0.7, main = '1/3 Bar Bulk Density (g/cm^3)\nDominant Component\n0-25cm') plot(mu2[['awc_r']], cex.main = 0.7, main = 'AWC (cm/cm)\nDominant Component\n0-25cm') plot(mu2[['ph1to1h2o_r']], cex.main = 0.7, main = 'pH 1:1 H2O\nDominant Component\n0-25cm')
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") )
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.