This document demonstrates how to use the soilDB package to access detailed soil series extent maps via SoilWeb. These maps can be directly displayed in R, overlayed on Google Maps, or saved as local files (e.g. shapefiles). Data returned from SoilWeb represent bounding-boxes that enclose SSURGO polygons associated with map units containing the search criteria. Bounding-boxes are snapped to 0.01 degree precision in order to reduce processing time and file size. Note that these files are cached server-side after the first request, and the cache is purged when SoilWeb is periodically synced to the official data.


With a recent version of R, it should be possible to get all of the packages that this tutorial depends on via:

# run these commands in the R console to install required packages
install.packages('soilDB', dep=TRUE)
install.packages('RColorBrewer', dep=TRUE)
install.packages('sf', dep=TRUE)
install.packages('maps', dep=TRUE)
install.packages('scales', dep=TRUE)
install.packages('mapview', dep=TRUE)


Illustrate the extent of SSURGO map units associated with the Amador series.

# load required libraries

Investigate the structure of soil series extent data.

# the result is a SpatialPolygonsDataFrame object
amador <- seriesExtent('amador')

# internal structure / class
str(amador, 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  1 obs. of  4 variables:
##   ..@ polygons   :List of 1
##   ..@ plotOrder  : int 1
##   ..@ bbox       : num [1:2, 1:2] -121.2 37.2 -120.1 38.6
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
# coordinate reference system in PROJ4 notation
## [1] "+proj=longlat +datum=WGS84 +no_defs"
# check attribute table, details:
slot(amador, 'data')
##        series acres gridsize   n
## AMADOR AMADOR 41359    0.001 963

Compare Amador, Pardee, and San Joaquin Extents

# soil series of interest
s <- c('amador', 'pardee', 'san joaquin')

# define some nice colors
cols <- brewer.pal('Set1', n = length(s))

# iterate over vector of soil series names
# and store spatial extents in a list
s.extent <- lapply(s, seriesExtent, timeout=120)

# combine into a single SpatialPolygonsDataFrame
s.extent <-'rbind', s.extent)

# convert into an sf object
s.extent <- st_as_sf(s.extent)

# very simple plot method for sf objects
plot(s.extent['series'], axes = FALSE, key.width = lcm(4.5), pal = cols)