wcs-ssurgo.R
In soilDB: Soil Database Interface

## ----setup, echo=FALSE, results='hide', warning=FALSE---------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  fig.align = 'center',
  fig.width = 8,
  fig.height = 6, 
  dpi = 36,
  tidy = FALSE,
  verbose = FALSE,
  # run when NASIS is defined or when R_SOILDB_SKIP_LONG_EXAMPLES is FALSE
  eval = isTRUE(try(local_NASIS_defined(), silent = TRUE)) ||
           !as.logical(Sys.getenv("R_SOILDB_SKIP_LONG_EXAMPLES", unset = "TRUE"))
)
options(width = 100, stringsAsFactors = FALSE)

## ----eval = FALSE---------------------------------------------------------------------------------
# install.packages(c('soilDB', 'terra', 'sf'))

## ----eval = FALSE---------------------------------------------------------------------------------
# install.packages(c('soilDB', 'terra', 'sf'),
#   repos = c('https://ncss-tech.r-universe.dev',
#             'https://rspatial.r-universe.dev',
#             'https://r-spatial.r-universe.dev')
# )

## ----eval = FALSE---------------------------------------------------------------------------------
# # select gSSURGO grid, 30m resolution
# x <- mukey.wcs(aoi = aoi, db = 'gssurgo', ...)
# 
# # select gNATSGO grid, 30m resolution
# x <- mukey.wcs(aoi = aoi, db = 'gnatsgo', ...)
# 
# # select RSS grid, 10m resolution
# x <- mukey.wcs(aoi = aoi, db = 'RSS', ...)
# 
# # select STATSGO2 grid, 300m resolution
# x <- mukey.wcs(aoi = aoi, db = 'statsgo', ...)

## ----eval = FALSE---------------------------------------------------------------------------------
# # select various ISSR-800 grids, details below
# x <- ISSR800.wcs(aoi = aoi, var = 'paws')

## ----fig.width = 5, fig.height = 5----------------------------------------------------------------
# library(terra)
# library(soilDB)
# 
# # example point, WGS84 coordinates
# p <- vect(
#   data.frame(
#     lon = -118.55639,
#     lat = 36.52578
#   ),
#   crs = "EPSG:4326"
# )
# 
# # 1000m buffer applied to WGS84 coordinate
# # radius defined in meters
# b <- buffer(p, 1000)
# 
# # query WCS
# # result is in EPSG:5070
# mu <- mukey.wcs(b, db = 'gSSURGO')
# 
# # inspect
# plot(
#   mu,
#   legend = FALSE,
#   axes = FALSE,
#   main = paste0(metags(mu)$value, collapse = " - ")
# )
# 
# # add buffer, after transforming to mukey grid CRS
# plot(project(b, "EPSG:5070"), add = TRUE)
# 
# # add original point, after transforming to mukey grid CRS
# plot(project(p, "EPSG:5070"), add = TRUE, pch = 16)

## ----fig.width = 8, fig.height = 7----------------------------------------------------------------
# library(sf)
# library(soilDB)
# library(terra)
# 
# # paste the five coordinates comprising the BBOX polygon here
# bb <- '-118.6609 36.4820,-118.6609 36.5972,-118.3979 36.5972,-118.3979 36.4820,-118.6609 36.4820'
# 
# # convert WKT string -> sfc geometry
# wkt <- sprintf('POLYGON((%s))', bb)
# x <- st_as_sfc(wkt)
# 
# # set coordinate reference system as GCS/WGS84
# st_crs(x) <- 4326
# 
# # query WCS
# mu <- mukey.wcs(x, db = 'gSSURGO')
# 
# # inspect
# plot(
#   mu,
#   legend = FALSE,
#   axes = FALSE,
#   main = paste0(metags(mu)$value, collapse = " - ")
# )
# 
# # add original BBOX, after transforming to mukey grid CRS
# plot(st_transform(x, 5070), add = TRUE)

## -------------------------------------------------------------------------------------------------
# # make a bounding box and assign a CRS (4326: GCS, WGS84)
# a <- st_bbox(
#   c(xmin = -114.16, xmax = -114.08, ymin = 47.65, ymax = 47.68),
#   crs = st_crs(4326)
# )
# 
# # fetch gSSURGO map unit keys at native resolution (30m)
# mu <- mukey.wcs(aoi = a, db = 'gssurgo')
# 
# # check:
# print(mu)
# 
# plot(
#   mu,
#   main = 'gSSURGO map unit keys',
#   sub = 'Albers Equal Area Projection',
#   axes = FALSE,
#   legend = FALSE
# )

## -------------------------------------------------------------------------------------------------
# # because mu is a SpatRaster, result is a SpatVector object (GCS WGS84)
# p <- SDA_spatialQuery(mu, what = 'mupolygon', geomIntersection = TRUE)

## -------------------------------------------------------------------------------------------------
# p <- project(p, crs(mu))

## ----fig.width = 8, fig.height = 7----------------------------------------------------------------
# plot(mu,
#      main = 'gSSURGO Grid (WCS)\nSSURGO Polygons (SDA)',
#      axes = FALSE,
#      legend = FALSE)
# plot(p, add = TRUE, border = 'white')
# mtext('CONUS Albers Equal Area Projection (EPSG:5070)', side = 1, line = 1)

## -------------------------------------------------------------------------------------------------
# # make a bounding box (in California) and assign a CRS (GCS WGS84 / EPSG:4326)
# a.CA <- st_bbox(c(
#   xmin = -121,
#   xmax = -120,
#   ymin = 37,
#   ymax = 38
# ), crs = st_crs(4326))
# 
# # fetch gSSURGO map unit keys at ~800m
# # nearest-neighbor resampling = this is a "preview"
# # result is a SpatRaster object
# x.800 <- mukey.wcs(aoi = a.CA, db = 'gssurgo', res = 800)
# 
# plot(
#   x.800,
#   main = 'A Preview of gSSURGO Map Unit Keys',
#   sub = 'Albers Equal Area Projection (800m)\nnearest-neighbor resampling',
#   axes = FALSE,
#   legend = FALSE
# )

## ----fig.width=8, fig.height=6--------------------------------------------------------------------
# # Coweeta Hydrologic Laboratory extent; specified in EPSG:5070
# a <- st_bbox(
#   c(xmin = 1129000, xmax = 1135000, ymin = 1403000, ymax = 1411000),
#   crs = st_crs(5070)
# )
# 
# # convert boundary sf polygon
# a <- st_as_sfc(a)
# 
# # gSSURGO grid: 30m resolution
# (x <- mukey.wcs(a, db = 'gSSURGO', res = 30))
# 
# # gNATSGO grid: 30m resolution
# (y <- mukey.wcs(a, db = 'gNATSGO', res = 30))
# 
# # RSS grid: 10m resolution
# (z <- mukey.wcs(a, db = 'RSS', res = 10))
# 
# # graphical comparison
# par(mfcol = c(1, 3))
# 
# 
# # gSSURGO
# plot(
#   x,
#   axes = FALSE,
#   legend = FALSE,
#   main = paste0(metags(x)$value, collapse = " - ")
# )
# plot(a, add = TRUE)
# 
# # gNATSGO
# plot(
#   y,
#   axes = FALSE,
#   legend = FALSE,
#   main = paste0(metags(y)$value, collapse = " - ")
# )
# plot(a, add = TRUE)
# 
# # RSS
# plot(
#   z,
#   axes = FALSE,
#   legend = FALSE,
#   main = paste0(metags(z)$value, collapse = " - "),
#   ext = x
# )
# plot(a, add = TRUE)

## ----fig.width=8, fig.height=6--------------------------------------------------------------------
# (statsgo <- mukey.wcs(a, db = 'statsgo', res = 300))
# 
# # graphical comparison
# par(mfcol = c(1, 2))
# 
# # gSSURGO
# plot(
#   x,
#   axes = FALSE,
#   legend = FALSE,
#   main = paste0(metags(mu)$value, collapse = " - ")
# )
# 
# # STATSGO
# plot(
#   statsgo,
#   axes = FALSE,
#   legend = FALSE,
#   main = paste0(metags(statsgo)$value, collapse = " - ")
# )

## ----fig.width = 6.5, fig.height=5----------------------------------------------------------------
# # paste your BBOX text here
# bb <- '-159.7426 21.9059,-159.7426 22.0457,-159.4913 22.0457,-159.4913 21.9059,-159.7426 21.9059'
# 
# # convert WKT string -> sfc geometry
# wkt <- sprintf('POLYGON((%s))', bb)
# x <- st_as_sfc(wkt, crs = 4326)
# 
# # query WCS
# mu <- mukey.wcs(x, db = 'hi_ssurgo')
# 
# # make NA (the ocean) blue
# plot(
#   mu,
#   legend = FALSE,
#   axes = FALSE,
#   main = paste0(metags(mu)$value, collapse = " - "),
#   colNA = 'royalblue'
# )

## ----eval=FALSE, include=FALSE--------------------------------------------------------------------
# # # check mu names
# # .is <- format_SQL_in_statement(cats(mu)[[1]]$mukey)
# # .sql <- sprintf("SELECT mukey, muname FROM mapunit WHERE mukey IN %s", .is)
# # knitr::kable(SDA_query(.sql))

## ----fig.width = 6.5, fig.height=5----------------------------------------------------------------
# # paste your BBOX text here
# bb <- '-65.7741 18.1711,-65.7741 18.3143,-65.5228 18.3143,-65.5228 18.1711,-65.7741 18.1711'
# 
# # convert WKT string -> sfc geometry
# wkt <- sprintf('POLYGON((%s))', bb)
# x <- st_as_sfc(wkt, crs = 4326)
# 
# # query WCS
# mu <- mukey.wcs(x, db = 'pr_ssurgo')
# 
# # make missing data (NA; the ocean) blue
# plot(
#   mu,
#   legend = FALSE,
#   axes = FALSE,
#   main = paste0(metags(mu)$value, collapse = " - "),
#   colNA = 'royalblue'
# )

## ----eval=FALSE, include=FALSE--------------------------------------------------------------------
# # # check mu names
# # .is <- format_SQL_in_statement(cats(mu)[[1]]$mukey)
# # .sql <- sprintf("SELECT mukey, muname FROM mapunit WHERE mukey IN %s", .is)
# # knitr::kable(SDA_query(.sql))

## -------------------------------------------------------------------------------------------------
# # make a bounding box and assign a CRS (4326: GCS, WGS84)
# a <- st_bbox(
#   c(xmin = -114.16, xmax = -114.08, ymin = 47.65, ymax = 47.68),
#   crs = st_crs(4326)
# )
# 
# # convert bbox to sf geometry
# a <- st_as_sfc(a)
# 
# # fetch gSSURGO map unit keys at native resolution (~30m)
# mu <- mukey.wcs(aoi = a, db = 'gssurgo')

## ----fig.width=8----------------------------------------------------------------------------------
# # copy example grid
# mu2 <- mu
# 
# # extract raster attribute table for thematic mapping
# (rat <- cats(mu2)[[1]])
# 
# # optionally use convenience function:
# # * returns all fields from muaggatt table
# # * along with map unit name
# # tab <- get_SDA_muaggatt(mukeys = as.numeric(rat$mukey), query_string = TRUE)
# 
# .sql <- paste0(
#   "SELECT mukey, aws050wta, aws0100wta FROM muaggatt WHERE mukey IN ",
#   format_SQL_in_statement(as.numeric(rat$mukey))
# )
# 
# # run query, result is a data.frame
# tab <- SDA_query(.sql)
# 
# # check
# head(tab)
# 
# # set raster categories
# levels(mu2) <- tab
# 
# # convert grid + RAT -> stack of property grids
# aws <- catalyze(mu2)
# 
# # plot, set a common range [0, 20] for both layers
# plot(
#   aws,
#   axes = FALSE,
#   cex.main = 0.7,
#   main = c(
#     'Plant Available Water Storage (cm)\nWeighted Mean over Components, 0-50cm',
#     'Plant Available Water Storage (cm)\nWeighted Mean over Components, 0-100cm'
#   ),
#   range = c(0, 20)
# )

## -------------------------------------------------------------------------------------------------
# # copy example grid
# mu2 <- mu
# 
# # extract RAT for thematic mapping
# rat <- cats(mu2)[[1]]
# 
# rules <- c('ENG - Construction Materials; Roadfill',
#            'AWM - Irrigation Disposal of Wastewater')
# 
# tab <- get_SDA_interpretation(
#   rulename = rules,
#   method = "Weighted Average",
#   mukeys = as.numeric(rat$mukey)
# )
# 
# # check
# head(tab)
# 
# # set ordered factor levels (for nice label/legend order)
# tab$class_ENGConstructionMaterialsRoadfill <- factor(
#   tab$class_ENGConstructionMaterialsRoadfill,
#   levels = c(
#     'Not suited',
#     'Poorly suited',
#     'Moderately suited',
#     'Moderately well suited',
#     'Well suited',
#     'Not Rated'
#   ),
#   ordered = TRUE
# )
# 
# par(mar = c(4, 12, 3, 3))
# boxplot(
#   rating_ENGConstructionMaterialsRoadfill ~ class_ENGConstructionMaterialsRoadfill,
#   cex.main = 0.7,
#   main = 'ENG - Construction Materials; Roadfill',
#   ylab = "",
#   data = tab,
#   horizontal = TRUE, # fuzzy ratings on X axis
#   las = 1            # rotate axis labels 90 degrees
# )

## ----fig.width=8----------------------------------------------------------------------------------
# vars <- c(
#   'rating_ENGConstructionMaterialsRoadfill',
#   'rating_AWMIrrigationDisposalofWastewater'
# )
# 
# # set raster categories
# levels(mu2) <- tab[, c('mukey', vars)]
# 
# rating <- catalyze(mu2)
# 
# # inspect
# plot(
#   rating,
#   axes = FALSE,
#   cex.main = 0.7,
#   main = c(
#     'Construction Materials; Roadfill\nWeighted Mean over Components',
#     'Irrigation Disposal of Wastewater\nWeighted Mean over Components'
#   )
# )

## ----fig.width = 8, fig.height = 6----------------------------------------------------------------
# # copy example grid
# mu2 <- mu
# 
# # extract RAT for thematic mapping
# rat <- cats(mu2)[[1]]
# 
# tab <- get_SDA_property(property = 'Corrosion of Steel',
#                         method = 'DOMINANT CONDITION',
#                         mukeys = as.integer(rat$mukey),
#                         miscellaneous_areas = TRUE)
# 
# # get soil data viewer standard colors for corsteel
# cols <- get_SDV_legend_elements("attributecolumnname = 'corsteel'")
# 
# # set raster categories
# levels(mu2) <- tab[, c('mukey', 'corsteel')]
# 
# # set active category
# activeCat(mu2) <- 'corsteel'
# 
# # plot
# plot(
#   mu2,
#   col = cols$hex[na.omit(match(unique(tab$corsteel), cols$label))],
#   axes = FALSE,
#   legend = "topleft"
# )

## -------------------------------------------------------------------------------------------------
# # https://casoilresource.lawr.ucdavis.edu/gmap/?loc=36.57666,-96.70175,z14
# # make a bounding box and assign a CRS (4326: GCS, WGS84)
# a <- st_bbox(
#   c(xmin = -96.7696, xmax = -96.6477,
#     ymin = 36.5477, ymax = 36.6139),
#   crs = st_crs(4326)
# )
# 
# # fetch gSSURGO map unit keys at native resolution (~30m)
# mu <- mukey.wcs(aoi = a, db = 'gssurgo')
# 
# plot(
#   mu,
#   legend = FALSE,
#   axes = FALSE,
#   cex.main = 0.7,
#   main = 'gSSURGO Map Unit Key Grid'
# )

## ----fig.width = 8, fig.height = 6----------------------------------------------------------------
# # copy example grid
# mu2 <- mu
# 
# # extract RAT for thematic mapping
# rat <- cats(mu2)[[1]]
# 
# # simplified parent material group name
# tab <- get_SDA_pmgroupname(mukeys = as.integer(rat$mukey),
#                            miscellaneous_areas = TRUE)
# 
# # set raster categories
# levels(mu2) <- tab[, c('mukey', 'pmgroupname')]
# 
# # set active category
# activeCat(mu2) <- 'pmgroupname'
# 
# plot(mu2, legend = "topleft", axes = FALSE)

## ----fig.width = 8, fig.height = 6----------------------------------------------------------------
# # copy example grid
# mu2 <- mu
# 
# # extract RAT for thematic mapping
# rat <- cats(mu2)[[1]]
# 
# # simplified parent material group name
# tab <- get_SDA_hydric(mukeys = as.integer(rat$mukey))
# 
# levels(mu2) <- tab[, c('mukey', 'HYDRIC_RATING')]
# 
# # set active category
# activeCat(mu2) <- 'HYDRIC_RATING'
# plot(mu2, legend = "topleft", axes = FALSE)

## -------------------------------------------------------------------------------------------------
# # extract RAT for thematic mapping
# rat <- cats(mu)[[1]]
# 
# # variables of interest
# vars <- c("dbthirdbar_r", "awc_r", "ph1to1h2o_r")
# 
# # get / aggregate specific horizon-level properties from SDA
# # be sure to see the manual page for this function
# tab <- get_SDA_property(property = vars,
#                         method = "Dominant Component (Numeric)",
#                         mukeys = as.integer(rat$mukey),
#                         top_depth = 0,
#                         bottom_depth = 25)
# 
# 
# # check
# head(tab)
# 
# # convert areasymbol into a factor easy plotting later
# tab$areasymbol <- factor(tab$areasymbol)
# 
# # set raster categories
# levels(mu) <- tab[, c('mukey', vars)]
# 
# # list variables in the RAT
# names(cats(mu)[[1]])
# 
# # convert categories associated with keys to values
# mu2 <- catalyze(mu)

## ----fig.width = 6, fig.height = 4----------------------------------------------------------------
# plot(mu2$awc_r)

## -------------------------------------------------------------------------------------------------
# plot(mu2[['dbthirdbar_r']], cex.main = 0.7,
#      main = '1/3 Bar Bulk Density (g/cm^3)\nDominant Component\n0-25cm')
# 
# plot(mu2[['awc_r']], cex.main = 0.7,
#      main = 'AWC (cm/cm)\nDominant Component\n0-25cm')
# 
# plot(mu2[['ph1to1h2o_r']], cex.main = 0.7,
#      main = 'pH 1:1 H2O\nDominant Component\n0-25cm')

## -------------------------------------------------------------------------------------------------
# # extract a BBOX like this from SoilWeb by pressing "b"
# bb <- '-91.6853 36.4617,-91.6853 36.5281,-91.5475 36.5281,-91.5475 36.4617,-91.6853 36.4617'
# wkt <- sprintf('POLYGON((%s))', bb)
# 
# # init sf object from WKT
# x <- st_as_sfc(wkt, crs = 4326)
# 
# # get gSSURGO grid here
# mu <- mukey.wcs(aoi = x, db = 'gssurgo')
# 
# # note SSA boundary
# plot(mu, legend = FALSE, axes = FALSE)

## ----fig.width = 8, fig.height = 6----------------------------------------------------------------
# # extract RAT for thematic mapping
# rat <- cats(mu)[[1]]
# 
# # variables of interest
# vars <- c("sandtotal_r", "silttotal_r", "claytotal_r")
# 
# # get thematic data from SDA
# # dominant component
# # depth-weighted average
# # sand, silt, clay (RV)
# tab <-  get_SDA_property(property = vars,
#                          method = "Dominant Component (Numeric)",
#                          mukeys = as.integer(rat$mukey),
#                          top_depth = 25,
#                          bottom_depth = 50)
# 
# # check
# head(tab)
# 
# # set raster categories
# levels(mu) <- tab[, c('mukey', vars)]
# 
# # convert mukey grid + RAT -> stack of numerical grids
# # retaining only sand, silt, clay via [[vars]]
# ssc <- catalyze(mu)
# 
# # create a copy of the grid
# texture.class <- ssc[[1]]
# names(texture.class) <- 'soil.texture'
# 
# # assign soil texture classes for the fine earth fraction
# # using sand and clay percentages
# values(texture.class) <- aqp::ssc_to_texcl(
#   sand = values(ssc[['sandtotal_r']]),
#   clay = values(ssc[['claytotal_r']]),
#   droplevels = FALSE
# )
# r <- c(ssc, texture.class)
# 
# # graphical check
# plot(
#   r,
#   cex.main = 0.7,
#   main = paste0(names(r), " - 25-50cm\nDominant Component")
# )