SDA_spatialQuery(
  geom,
  what = "mukey",
  geomIntersection = FALSE,
  db = c("SSURGO", "STATSGO")
)

if(requireNamespace("curl") &
   curl::has_internet() & 
   requireNamespace("sp") &
   requireNamespace("raster") 
   ) {

library(aqp)
library(sp)
library(raster)

## query at a point

# example point
p <- SpatialPoints(cbind(x = -119.72330, y = 36.92204), 
                   proj4string = CRS('+proj=longlat +datum=WGS84'))

# query map unit records at this point
res <- SDA_spatialQuery(p, what = 'mukey')

# convert results into an SQL "IN" statement
# useful when there are multiple intersecting records
mu.is <- format_SQL_in_statement(res$mukey)

# composite SQL WHERE clause
sql <- sprintf("mukey IN %s", mu.is)

# get commonly used map unit / component / chorizon records
# as a SoilProfileCollection object
# confusing but essential: request that results contain `mukey`
# with `duplicates = TRUE`
x <- fetchSDA(sql, duplicates = TRUE)

# safely set texture class factor levels
# by making a copy of this column
# this will save in lieu of textures in the original
# `texture` column
horizons(x)$texture.class <- factor(x$texture, levels = SoilTextureLevels())

# graphical depiction of the result
plotSPC(x, color='texture.class', label='compname', 
        name='hzname', cex.names = 1, width=0.25, 
        plot.depth.axis=FALSE, hz.depths=TRUE, 
        name.style='center-center'
)



## query mukey + geometry that intersect with a bounding box

# define a bounding box: xmin, xmax, ymin, ymax
#
#         +-------------------(ymax, xmax)
#         |                        |
#         |                        |
#     (ymin, xmin) ----------------+
b <- c(-119.747629, -119.67935, 36.912019, 36.944987)

# convert bounding box to WKT
bbox.sp <-as(extent(b), 'SpatialPolygons')
proj4string(bbox.sp) <- '+proj=longlat +datum=WGS84'

# results contain associated map unit keys (mukey)
# return SSURGO polygons, after intersection with provided BBOX
ssurgo.geom <- SDA_spatialQuery(
  bbox.sp, 
  what = 'geom', 
  db = 'SSURGO', 
  geomIntersection = TRUE
)

# return STATSGO polygons, after intersection with provided BBOX
statsgo.geom <- SDA_spatialQuery(
  bbox.sp, 
  what = 'geom', 
  db = 'STATSGO', 
  geomIntersection = TRUE
)

# inspect results
par(mar = c(0,0,3,1))
plot(ssurgo.geom, border = 'royalblue')
plot(statsgo.geom, lwd = 2, border = 'firebrick', add = TRUE)
plot(bbox.sp, lwd = 3, add = TRUE)
legend(
  x = 'topright', 
  legend = c('BBOX', 'STATSGO', 'SSURGO'), 
  lwd = c(3, 2, 1),
  col = c('black', 'firebrick', 'royalblue'),
)


# quick reminder that STATSGO map units often contain many components
# format an SQL IN statement using the first STATSGO mukey
mu.is <- format_SQL_in_statement(statsgo.geom$mukey[1])

# composite SQL WHERE clause
sql <- sprintf("mukey IN %s", mu.is)

# get commonly used map unit / component / chorizon records
# as a SoilProfileCollection object
x <- fetchSDA(sql)

# tighter figure margins
par(mar = c(0,0,3,1))


# organize component sketches by national map unit symbol
# color horizons via awc
# adjust legend title
# add alternate label (vertical text) containing component percent
# move horizon names into the profile sketches
# make profiles wider
groupedProfilePlot(
  x, 
  groups = 'nationalmusym', 
  label = 'compname', 
  color = 'awc_r', 
  col.label = 'Available Water Holding Capacity (cm / cm)',
  alt.label = 'comppct_r',
  name.style = 'center-center',
  width = 0.3
)


mtext(
  'STATSGO (1:250,000) map units contain a lot of components!', 
  side = 1, 
  adj = 0, 
  line = -1.5, 
  at = 0.25, 
  font = 4
)
 }