Introduction

This document demonstrates how to use the soilDB package to access “soil taxa” mapping (800m grid resolution, CONUS only for now), as part of the ISSR-800 project. Taxa grids were developed from the current SSURGO snapshot using a 800x800 meter grid. Taxa proportions were computed from the geometric intersection of 800m grid and SSURGO polygons (map unit area) and associated component data (component percentage).

Setup

Get required packages.

# run these commands in the R console
install.packages('soilDB', dep = TRUE)
install.packages('mapview', dep = TRUE)
install.packages('sf', dep = TRUE)
install.packages('terra', dep = TRUE)
install.packages('mapview', dep = TRUE)
install.packages('SoilTaxonomy', dep = TRUE)
# load required libraries
library(soilDB)
library(SoilTaxonomy)
library(terra)
library(mapview)
library(sf) 
library(spData)

Prepare state boundaries for CONUS.

# CRS defs used by SEE grids
# EPSG:5070
crs_lower48 <- 5070

# use lower '48 state borders as context
# transform to the CRS of the SEE taxa grids
data("us_states")
us_states <- st_transform(us_states, crs_lower48)

# convert to SpatVect
us_states <- vect(us_states)

# tight bounding box around CONUS
b <- ext(us_states)

Examples

Results are SpatRaster as provided by the terra package.

All vertisols.

taxa <- 'vertisols'
x <- taxaExtent(taxa, level = 'order')
a <- aggregate(x, fact = 5, fun = mean, na.rm = TRUE)

# plot using full extent of CONUS
plot(a, axes = FALSE, col = hcl.colors(25), ext = b, mar = c(1, 1, 3, 4), main = names(a))
lines(us_states, lwd = 2)

Xeralfs.

taxa <- 'xeralfs'
x <- taxaExtent(taxa, level = 'suborder')
a <- aggregate(x, fact = 5, fun = mean, na.rm = TRUE)

# plot using full extent of CONUS
plot(a, axes = FALSE, col = hcl.colors(25), ext = b, mar = c(1, 1, 3, 4), main = names(a))
lines(us_states, lwd = 2)

Durixeralfs.

taxa <- 'durixeralfs'
x <- taxaExtent(taxa, level = 'greatgroup')
a <- aggregate(x, fact = 5, fun = mean, na.rm = TRUE)

# zoomed to taxa extent
plot(a, axes = FALSE, col = hcl.colors(25), mar = c(1, 1, 3, 4), main = names(a))
lines(us_states, lwd = 2)

Abruptic durixeralfs.

taxa <- 'abruptic durixeralfs'
x <- taxaExtent(taxa, level = 'subgroup')
a <- aggregate(x, fact = 5, fun = mean, na.rm = TRUE)

# zoomed to taxa extent
plot(a, axes = FALSE, col = hcl.colors(25), mar = c(1, 1, 3, 4), main = names(a))
lines(us_states, lwd = 2)

Cecil soil series.

s <- seriesExtent('cecil', type = 'raster')
a <- aggregate(s, fact = 5, fun = mean, na.rm = TRUE)

# CONUS
plot(a, axes = FALSE, col = hcl.colors(25), ext = b, mar = c(1, 1, 3, 4), main = names(a))
lines(us_states, lwd = 2)

Formative Elements

x <- taxaExtent('pale', level = 'greatgroup', formativeElement = TRUE)
a <- aggregate(x, fact = 5, fun = mean, na.rm = TRUE)

# zoomed to full extent
plot(a, axes = FALSE, col = hcl.colors(25), ext = b, mar = c(1, 1, 3, 4), main = names(a))
lines(us_states, lwd = 2)

x <- taxaExtent('grossarenic', level = 'subgroup', formativeElement = TRUE)
a <- raster::aggregate(x, fact = 5, fun = mean, na.rm = TRUE)

# zoomed to taxa extent
plot(a, axes = FALSE, col = hcl.colors(25), mar = c(1, 1, 3, 4), main = names(a))
lines(us_states, lwd = 2)

Interactive mapping from R

This is based on the mapview package.

x <- taxaExtent('xeralfs', level = 'suborder')
a <- aggregate(x, fact = 5, fun = mean, na.rm = TRUE)

# requires raster package representation
mapview(raster::raster(a), col.regions = hcl.colors, na.color = NA, use.layer.names = TRUE)
x <- taxaExtent('durixeralfs', level = 'greatgroup')
a <- aggregate(x, fact = 5, fun = mean, na.rm = TRUE)

mapview(raster::raster(a), col.regions = hcl.colors, na.color = NA, use.layer.names = TRUE)
s <- seriesExtent('san joaquin', type = 'raster')
a <- aggregate(s, fact = 5, fun = mean, na.rm = TRUE)

mapview(raster::raster(a), col.regions = hcl.colors, na.color = NA, use.layer.names = TRUE)

Aggregating and mapping multiple soil taxa

This example demonstrates a way of using pattern matching on a soil taxonomic level to filter a list of soil taxa that can be aggregated into one raster and mapped. For example, if we wanted to see all the Andisols and any Andic or Vitrandic subgroups we could use this process. If the pattern match grabs taxa that shouldn’t be include then you can use the same process to exclude taxa that should not be in the aggregated extent raster (ie. ‘Kandic’ elements).

Requires latest terra (1.7-31).

# unique taxa
data("ST_unique_list", package = 'SoilTaxonomy')

# filter to include taxa based on a specified pattern
idx <- grep('and', ST_unique_list$subgroup)
STgroup <- ST_unique_list$subgroup[idx]
length(STgroup)

# filter to remove specified orders and any other based on pattern matching
idx <- grep('ids|ods|ults|ox|els|kandi', STgroup)
STgroup <- STgroup[-idx]
length(STgroup)

# inspect
head(STgroup)


# loop thru list and fetch taxa extent rasters
# not all taxa have been mapped, messages printed
m <- lapply(STgroup, taxaExtent, level = 'subgroup', timeout = 120)

# remove NULL
m <- m[which(!sapply(m, is.null))]

# combine into SpatRasterCollection
# 150 spatRast objects
x <- sprc(m)

# combine into single SpatRaster
x <- mosaic(x, fun = 'sum')

# values > 100 are rounding errors
x[x > 100] <- 100

# give merged raster a name
names(x) <- 'Andisols and Andic Subgroups'

# aggregate to 5x larger grid, sum of cell percent cover
a <- aggregate(x, fact = 5, fun = sum, na.rm = TRUE)

# rescale percent cover to larger grid size
a <- a / 5^2

# CONUS map
plot(a, axes = FALSE, col = hcl.colors(25), ext = b, mar = c(1, 1, 3, 4), main = names(a))
lines(us_states, lwd = 2)


# slippy map
if(requireNamespace("mapview")) {
  mapview::mapview(raster::raster(x), col.regions = hcl.colors, na.color = NA, use.layer.names = TRUE)
}

This document is based on soilDB version 2.8.5.