Getting, Plotting, and Saving Soil Series Extent Data from SoilWeb

2013-11-21 Dylan Beaudette

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("soilDB", repos = "http://R-Forge.R-project.org")  # most recent copy from r-forge
install.packages("dismo", dep = TRUE)
install.packages("RColorBrewer", dep = TRUE)
install.packages("maps", dep = TRUE)

Examples

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

# load required libraries
library(soilDB)
library(dismo)

# plot map
seriesExtentAsGmap("amador")

plot of chunk example-1plot of chunk example-1

Illustrate the extent of SSURGO map units associated with several soil series from MLRA 17 and 18.

# load required libraries
library(maps)
library(RColorBrewer)

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

# get the spatial extent of select series
amador <- seriesExtent("amador")
pardee <- seriesExtent("pardee")
san.joaquin <- seriesExtent("san joaquin", timeout = 120)

# 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
plot(amador, col = cols[1], add = TRUE)
plot(pardee, col = cols[2], add = TRUE)
plot(san.joaquin, col = cols[3], add = TRUE)

# add a neatline around the map
box()

# make a simple legend
legend("topright", col = cols, pch = 15, legend = c("Amador", "Pardee", "San Joaquin"))

plot of chunk example-2

Investigate the structure of soil series extent data.

# internal structure / class
str(amador, 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  1 obs. of  2 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 slots

# coordinate reference system in PROJ4 notation
proj4string(amador)
## [1] "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"

Save soil series extent data as a shapefile.

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

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

This document is based on aqp version 1.5-6 and soilDB version 1.2-4.