Introduction

This is a short tutorial on how to interact with the Soil Data Access (SDA) web-service using R. Queries are written using a dialect of SQL. On first glance SQL appears similar to the language used to write NASIS queries and reports, however, these are two distinct languages. Soil Data Access is a "window" into the spatial and tabular data associated with the current SSURGO snapshot. Queries can contain spatial and tabular filters. If you are new to SDA or SQL, have a look at this page.

If this is your first time using SDA, please see a related tutorial to get started.

Additional tips on advanced spatial queries can be found here.

[details pending]

Follow along with the blocks of code below by copying / pasting into a new R "script" document. Each block of commands can be run by pasting them into the R console, or by "stepping through" lines of code by moving the cursor to the top of a block (in the R script panel) and repeatedly pressing ctrl + enter.

Install Required R Packages

You only need to do this once. If you haven't installed these packages, then copy the code below and paste into the RStudio "console" pane.

# run these commands in the R console
# stable version from CRAN + dependencies
install.packages("httr", dep=TRUE)
install.packages("soilDB", dep=TRUE)
install.packages("rgdal", dep = TRUE)
install.packages("raster", dep = TRUE)
install.packages("rgeos", dep = TRUE)

Critical Note: SSURGO vs. STATSGO

SSURGO (1:24k soil survey) and STATSGO (1:250k soil survey) records are stored together in SDA. Therefore, it is critical that evey query to SDA include some kind of filter for selecting the appropriate records. Filtering strategies include:

  • explicit exclusion of STATSGO records, via legend.areasymbol != 'US' in the WHERE clause
  • implicit exclusion of STATSGO records, via SSURGO areasymbol in the WHERE clause
  • spatial queries using SDA helper functions: e.g. SDA_Get_Mukey_from_intersection_with_WktWgs84()
  • explicit selection of SSURGO / STATSGO records by record ID: e.g. mukey, cokey, etc.

Simple Queries

These queries implicitly filter-out STATSGO via SDA_Get_MupolygonWktWgs84_from_Mukey.

library(soilDB)
library(rgeos)
library(sp)
library(raster)
library(maps)


# get polygons for a single mukey
q <- "SELECT G.MupolygonWktWgs84 as geom, '462594' as mukey from SDA_Get_MupolygonWktWgs84_from_Mukey('462594') as G"
res <- SDA_query(q)

# result is a data.frame, "MupolygonWktWgs84" contains WKT representation of geometry
str(res)
## 'data.frame':    38 obs. of  2 variables:
##  $ geom : chr  "POLYGON ((-120.97991145612417 37.742782233615195, -120.9798283067521 37.742828501994246, -120.97975957521625 37"| __truncated__ "POLYGON ((-120.72864825039305 37.638143059003831, -120.72866384079708 37.638088744767913, -120.72868026919249 3"| __truncated__ "POLYGON ((-120.95867387378898 37.624002099580515, -120.95844789742932 37.624065299195877, -120.95821353920952 3"| __truncated__ "POLYGON ((-120.84277744235267 37.624586493906556, -120.84266579435683 37.6247088698527, -120.84259756464191 37."| __truncated__ ...
##  $ mukey: int  462594 462594 462594 462594 462594 462594 462594 462594 462594 462594 ...
##  - attr(*, "SDA_id")= chr "Table"
# convert to SPDF
s <- processSDA_WKT(res)

# check
head(s@data)
##   gid  mukey
## 1   1 462594
## 2   2 462594
## 3   3 462594
## 4   4 462594
## 5   5 462594
## 6   6 462594
plot(s)

# get polygons associated with map units that contain "amador" as a major component
q <- "select G.MupolygonWktWgs84 as geom, mapunit.mukey, muname
FROM mapunit
CROSS APPLY SDA_Get_MupolygonWktWgs84_from_Mukey(mapunit.mukey) as G
WHERE mukey IN (SELECT DISTINCT mukey FROM component WHERE compname like 'amador%' AND majcompflag = 'Yes')"

# result is a data.frame, "MupolygonWktWgs84" contains WKT representation of geometry
res <- SDA_query(q)
str(res)
## 'data.frame':    396 obs. of  3 variables:
##  $ geom  : chr  "POLYGON ((-121.15252230315829 38.580559769164942, -121.15239322126237 38.5805240627318, -121.15232348375943 38."| __truncated__ "POLYGON ((-121.05686027711923 38.461657600967762, -121.05679087490459 38.4616014427027, -121.05679053944492 38."| __truncated__ "POLYGON ((-121.11379721766438 38.4712833760718, -121.1136671314583 38.471317909610491, -121.11361315117446 38.4"| __truncated__ "POLYGON ((-121.11084293183234 38.4682622025382, -121.11072893724598 38.468321546696139, -121.11059365384207 38."| __truncated__ ...
##  $ mukey : int  461845 461845 461980 461980 461845 461845 461845 461845 461845 461845 ...
##  $ muname: chr  "Amador-Gillender complex, 2 to 15 percent slopes" "Amador-Gillender complex, 2 to 15 percent slopes" "Vleck-Amador-Pits, mine complex, 15 to 50 percent slopes" "Vleck-Amador-Pits, mine complex, 15 to 50 percent slopes" ...
##  - attr(*, "SDA_id")= chr "Table"
# convert to SPDF
s <- processSDA_WKT(res)

# check: OK
head(s@data)
##   gid  mukey                                                   muname
## 1   1 461845         Amador-Gillender complex, 2 to 15 percent slopes
## 2   2 461845         Amador-Gillender complex, 2 to 15 percent slopes
## 3   3 461980 Vleck-Amador-Pits, mine complex, 15 to 50 percent slopes
## 4   4 461980 Vleck-Amador-Pits, mine complex, 15 to 50 percent slopes
## 5   5 461845         Amador-Gillender complex, 2 to 15 percent slopes
## 6   6 461845         Amador-Gillender complex, 2 to 15 percent slopes
# map
par(mar=c(0,0,0,0))
map('county', 'California', xlim=c(-123.25, -118.75), ylim=c(36.5, 39))
plot(s, add=TRUE, col='royalblue', border='royalblue')

This document is based on soilDB version 2.2-6.