Writing SDA Queries that Return Geometry

2016-08-17
Dylan Beaudette

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)

Simple Queries

Now that you have the required packages, load them into the current R session.

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.93265560125401 37.6131532349993, -120.93258720477149 37.613285333628, -120.93255434836267 37.613331937167388, -1"| __truncated__ "POLYGON ((-121.06518726313537 37.5862540180816, -121.06513915045151 37.586182939883336, -121.06513479196238 37.586101802698359,"| __truncated__ "POLYGON ((-121.11821940874651 37.593238990536442, -121.11816626749126 37.593204122144805, -121.11814732430042 37.59318232927007"| __truncated__ "POLYGON ((-120.62789456723003 37.63902668716419, -120.62776364240879 37.6390587059678, -120.62769692244254 37.639063902414712, "| __truncated__ ...
##  $ mukey: int  462594 462594 462594 462594 462594 462594 462594 462594 462594 462594 ...
# convert to SPDF
s <- processSDA_WKT(res)

# check
head(s@data)
gid mukey
1 462594
2 462594
3 462594
4 462594
5 462594
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')"

res <- SDA_query(q)
str(res)
## 'data.frame':    262 obs. of  3 variables:
##  $ geom  : chr  "POLYGON ((-121.05686027711923 38.461657600967762, -121.05679087490459 38.4616014427027, -121.05679053944492 38.461583337311417,"| __truncated__ "POLYGON ((-121.03052749990923 38.319438143011652, -121.03057661740183 38.319477873008218, -121.03060545092197 38.31952414109228"| __truncated__ "POLYGON ((-121.10304826277805 38.47214151628399, -121.1029176725373 38.472175211524956, -121.10254082289315 38.472307142353067,"| __truncated__ "POLYGON ((-121.06601716300834 38.372910499366796, -121.06585689999224 38.372808743001038, -121.06551743292594 38.37262501223038"| __truncated__ ...
##  $ mukey : int  461845 461845 461845 461845 461845 461845 461845 461845 461845 461845 ...
##  $ muname: chr  "Amador-Gillender complex, 2 to 15 percent slopes" "Amador-Gillender complex, 2 to 15 percent slopes" "Amador-Gillender complex, 2 to 15 percent slopes" "Amador-Gillender complex, 2 to 15 percent slopes" ...
s <- processSDA_WKT(res)

# check: OK
head(s@data)
gid mukey muname
1 461845 Amador-Gillender complex, 2 to 15 percent slopes
2 461845 Amador-Gillender complex, 2 to 15 percent slopes
3 461845 Amador-Gillender complex, 2 to 15 percent slopes
4 461845 Amador-Gillender complex, 2 to 15 percent slopes
5 461845 Amador-Gillender complex, 2 to 15 percent slopes
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)


This document is based on soilDB version 1.8-1.