Introduction

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 rebuilt when SoilWeb is periodically synced to the official data.

Installation

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('terra', dep=TRUE)
install.packages('sf', dep=TRUE)
install.packages('maps', dep=TRUE)
install.packages('scales', dep=TRUE)
install.packages('mapview', dep=TRUE)

Examples

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

# load required libraries
library(soilDB)
library(sf)
library(terra)
library(maps)
library(rasterVis)

Investigate the structure of soil series extent data.

# the result is an af object
amador <- seriesExtent('amador')

# internal structure / class
class(amador)
## [1] "sf"         "data.frame"
# coordinate reference system in PROJ4 notation
st_crs(amador)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
# check attribute table, details: http://casoilresource.lawr.ucdavis.edu/see/
head(amador)
## Simple feature collection with 1 feature and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -121.193 ymin: 37.183 xmax: -120.078 ymax: 38.588
## Geodetic CRS:  WGS 84
##        series acres gridsize   n                       geometry
## AMADOR AMADOR 41359    0.001 963 MULTIPOLYGON (((-120.489 37...

Compare Amador, Pardee, and San Joaquin Extents

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

# define some nice colors
cols <- c("#E41A1C", "#377EB8", "#4DAF4A")

# 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 sf object
s.extent <- do.call('rbind', s.extent)

# very simple plot method for sf objects
# https://r-spatial.github.io/sf/articles/sf5.html
plot(s.extent['series'], axes = FALSE, key.width = lcm(4.5), pal = cols)

Interactive mapping from R

This is based on the mapview package.

library(mapview)
mapview(s.extent, zcol = 'series', legend = TRUE)

Exporting Series Extent Data

Save soil series extent data in multiple formats.

# save using the coordinate reference system associated with this object (GCS WGS84)
write_sf(s.extent, dsn = 'amador-pardee-san-joaquin-extent.shp')

# save as KML for use in Google Earth
# writeOGR(s.extent, dsn='amador-pardee-san-joaquin-extent.kml', layer='amador-pardee-san-joaquin', driver='KML')

# project to UTM zone 10 NAD83 and save
# s.extent.utm <- spTransform(s.extent, CRS('+proj=utm +zone=10 +datum=NAD83'))
# writeOGR(s.extent.utm, dsn='.', layer='amador-pardee-san-joaquin-extent-utm', driver='ESRI Shapefile')

Graphical Display of Joint Series Extent

Using vector features, works with all SSURGO data.

# need this for alpha()
library(scales)

# define some nice colors
cols <- c("#E41A1C", "#377EB8", "#4DAF4A")

# list of series for investigation of joint occurrence
s <- c('pentz', 'peters', 'pardee', 'amador', 'san joaquin', 'redding', 'auburn')

# get list of extents
e <- lapply(s, seriesExtent)

# flatten to single sf object
e <- do.call('rbind', e)

# map of CA
par(mar = c(1,1,1,1))
map('state', 'ca')

# sf plot method
# polygon extents are "stacked" with transparency
plot(st_geometry(e), border = NA, col = alpha(cols[2], 0.25), add = TRUE)

# finish map
box()
title(main = 'Joint Occurrence', line = 1.25)
mtext(side = 1, text = as.character(Sys.Date()), line = 0)

Join occurrence of Zook and siblings.

# using all siblings of a given soil series
s <- 'zook'
sib <- siblings(s)
sib.names <- c(s, sib$sib$sibling)

# get list of series extent grids
# note these all have different bounding-boxes (raster extents)
g <- lapply(sib.names, function(i) {
  se <- suppressWarnings(seriesExtent(i, type = 'raster'))
  # threshold at 15%
  # se <- se >= 15
  
  return(se)
})


# assemble collection, 44 spatRast objects
g <- sprc(g)

# merge into single SpatRaster
# compute sum whenever there is overlap
m <- mosaic(g, fun = 'sum')

# values > 100% are rounding errors
# truncate at 100%
m[m > 100] <- 100

levelplot(
  m,
  col.regions = hcl.colors,
  main = 'Joint Occurrence',
  margin = FALSE, 
  scales = list(draw = FALSE)
)

# generalize via focal sum
a <- aggregate(m, fact = 5, fun = sum, na.rm = TRUE)

# re-scale to fractions of 5x larger grid size
a <- a / 5^2

levelplot(
  a,
  col.regions = hcl.colors,
  main = 'Joint Occurrence',
  margin = FALSE, 
  scales = list(draw = FALSE)
)

Boomer Soil Series Extent

Lets investigate the spatial extent and MLRA overlap for the Boomer soil series. We will use the SoilWeb SEE data for an estimate of the spatial extent and Soil Data Access for MLRA overlap information.

Get extent and make a simple map.

# define some nice colors
cols <- c("#E41A1C", "#377EB8", "#4DAF4A")

# get extent
boomer <- seriesExtent('boomer')

# basic figure
par(mar=c(1,0,1,0))
map(database='county', regions='california')
plot(st_geometry(boomer), col=cols[2], border=cols[2], add=TRUE)
box()
title(main='Boomer Series Extent', line=1.25)
mtext(side=1, text=as.character(Sys.Date()), line=0)

Investigate MLRA overlap for the Boomer series via SDA (official SSURGO data). These results are based on the assumption that the “muaoverlap” table is correctly populated.

res <- SDA_query("SELECT compname, areasymbol as mlra, areaname as mlra_name, count(compname) as n_components from laoverlap JOIN muaoverlap ON laoverlap.lareaovkey = muaoverlap.lareaovkey JOIN component ON muaoverlap.mukey = component.mukey WHERE compname = 'boomer' AND areatypename = 'MLRA' GROUP BY compname, areasymbol, areaname ORDER BY count(compname) DESC")
knitr::kable(res, row.names = FALSE)
compname mlra mlra_name n_components
Boomer 22A Sierra Nevada and Tehachapi Mountains 71
Boomer 5 Siskiyou-Trinity Area 66
Boomer 15 Central California Coast Range 48
Boomer 18 Sierra Nevada Foothills 44
Boomer 22A Sierra Nevada Mountains 18
Boomer 22B Southern Cascade Mountains 14
Boomer 14 Central California Coastal Valleys 13
Boomer 4B Coastal Redwood Belt 11
Boomer 17 Sacramento and San Joaquin Valleys 10
Boomer 21 Klamath and Shasta Valleys and Basins 6
Boomer 31 Lower Colorado Desert 2
Boomer 20 Southern California Mountains 2

That seems strange, which map units are in MLRA 31?

res <- SDA_query("SELECT component.mukey, muname, compname, comppct_r from laoverlap JOIN muaoverlap ON laoverlap.lareaovkey = muaoverlap.lareaovkey JOIN mapunit ON muaoverlap.mukey = mapunit.mukey JOIN component ON muaoverlap.mukey = component.mukey WHERE compname = 'boomer' AND areatypename = 'MLRA' AND areasymbol = '31'")
knitr::kable(res, row.names = FALSE)
mukey muname compname comppct_r
660462 Hotaw-Crouch-Boomer (s1015) Boomer 1
660462 Hotaw-Crouch-Boomer (s1015) Boomer 15

This document is based on aqp version 2.0, soilDB version 2.7.8 and terra version 1.7-31.