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 purged 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
install.packages('soilDB', dep=TRUE) # stable version from CRAN + dependencies
install.packages('dismo', dep=TRUE)
install.packages('RColorBrewer', dep=TRUE)
install.packages('maps', dep=TRUE)

# you will need the latest version of soilDB
devtools::install_github("ncss-tech/soilDB", dependencies=FALSE, upgrade_dependencies=FALSE)

Examples

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

# load required libraries
library(soilDB)
library(rgdal)
library(dismo)
library(maps)
library(RColorBrewer)

Get series extent for the Amador series and plot on a Google Map.

seriesExtentAsGmap('amador')

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
proj4string(amador)

## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

# check attribute table, details: http://casoilresource.lawr.ucdavis.edu/see/
slot(amador, 'data')

##        series acres gridsize   n
## AMADOR AMADOR 23542    0.001 788

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

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

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

# plot the results:
# set figure margins to 0
par(mar=c(0,0,0,0))

# plot select counties from California
map(database='county', regions='california,calaveras|tuolumne|amador|san joaquin|stanislaus|merced|mariposa|sacramento')

# add each soil series extent object
for(i in 1:nrow(s.extent)) {
  plot(s.extent[i, ], col=cols[i], add=TRUE)
}

# add a neatline around the map
box()

# make a simple legend
legend('topright', col=cols, pch=15, legend=s.extent$series)

# save using the coordinate reference system associated with this object (GCS WGS84)
writeOGR(s.extent, dsn='.', layer='amador-pardee-san-joaquin-extent', driver='ESRI Shapefile')

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

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.

boomer <- seriesExtent('boomer')
par(mar=c(1,0,1,0))
map(database='county', regions='california')
plot(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")
kable(res, row.names = FALSE)

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'")
kable(res, row.names = FALSE)

Save extent for Boomer series as shapefile (GCS/WGS84).

writeOGR(boomer, dsn='.', layer='boomer-extent', driver='ESRI Shapefile')

This document is based on aqp version 1.9.15 and soilDB version 1.8-8.