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.
This document is based on
soilDB version 1.8-1.