library(aqp)
library(soilDB)
library(sf)
library(sharpshootR)
library(SoilTaxonomy)
## possible AOIs defined using a bounding-box via SoilWeb
# WI: many ties
bb <- '-97.0983 39.3808,-97.0983 39.4127,-97.0282 39.4127,-97.0282 39.3808,-97.0983 39.3808'
# KS069
bb <- '-100.5534 37.9177,-100.5534 37.9822,-100.4389 37.9822,-100.4389 37.9177,-100.5534 37.9177'
## assemble AOI polygon into WKT
wkt <- sprintf('POLYGON((%s))', bb)
## init sf polygon
# GCS WGS84
x <- st_as_sfc(wkt, crs = 4326)
## get overlapping map unit keys
# could also use SDA_query() with more elaborate SQL
m <- SDA_spatialQuery(x, what = 'mukey')
## compose SQL to return component details for these map unit keys
# return only:
# * map units overlapping with BBOX
# * major components
# * no misc. areas that might share name with a poorly-named soil series
sql <- sprintf(
"SELECT mukey, cokey, compname, compkind, comppct_r
FROM component
WHERE mukey IN %s
-- AND majcompflag = 'Yes'
AND compkind != 'Miscellaneous area'
", format_SQL_in_statement(as.integer(m$mukey))
)
## send to SDA, result is a data.frame
s <- SDA_query(sql)
## get OSD morphology + extended summaries
osd <- fetchOSD(unique(s$compname), extended = TRUE)
## check out results
str(osd, 1)
## List of 18
## $ SPC :Formal class 'SoilProfileCollection' [package "aqp"] with 9 slots
## $ competing :'data.frame': 61 obs. of 3 variables:
## $ geog_assoc_soils:'data.frame': 75 obs. of 2 variables:
## $ geomcomp :'data.frame': 12 obs. of 9 variables:
## $ hillpos :'data.frame': 12 obs. of 8 variables:
## $ mtnpos : logi FALSE
## $ terrace :'data.frame': 6 obs. of 5 variables:
## $ flats :'data.frame': 7 obs. of 7 variables:
## $ shape_across :'data.frame': 14 obs. of 8 variables:
## $ shape_down :'data.frame': 14 obs. of 8 variables:
## $ pmkind :'data.frame': 31 obs. of 5 variables:
## $ pmorigin :'data.frame': 5 obs. of 5 variables:
## $ mlra :'data.frame': 99 obs. of 4 variables:
## $ ecoclassid :'data.frame': 59 obs. of 5 variables:
## $ climate.annual :'data.frame': 112 obs. of 12 variables:
## $ climate.monthly :'data.frame': 336 obs. of 14 variables:
## $ NCCPI :'data.frame': 14 obs. of 16 variables:
## $ soilweb.metadata:'data.frame': 24 obs. of 2 variables:
# single iteration of hydrologic ordering
h1 <- hydOrder(osd$shape_across, g = 'shape', clust = TRUE)
# return trace log for eval of objective function
# increase max iterations
h2 <- iterateHydOrder(osd$shape_across, 'shape', maxIter = 100, verbose = TRUE, trace = TRUE)
# compare
h1$match.rate
## [1] 0.7857143
h2$match.rate
## [1] 0.7857143
# inspect objective function evolution
tr <- h2$trace
obj <- sapply(tr, '[[', 'obj')
hist(obj, las = 1, xlab = 'Objective Function\nLower is Better')
match.rate <- sapply(tr, '[[', 'match.rate')
hist(match.rate, las = 1, xlab = 'Matching Rate')
Note that profiles deeper than 180cm are marked as truncated via ragged bottoms.
# provide additional arguments to aqp::plotSPC() via options
options(.aqp.plotSPC.args = list(max.depth = 190))
par(mar = c(1, 0, 1, 2))
plotGeomorphCrossSection(osd, type = 'line', maxIter = 100, j.amount = 0.05, verbose = TRUE)
plotGeomorphCrossSection(osd, type = 'bar', maxIter = 100, j.amount = 0.05, verbose = TRUE)
par(mar = c(1, 0, 1, 2))
plotGeomorphCrossSection(osd, type = 'line', clust = FALSE)
plotGeomorphCrossSection(osd, type = 'bar', clust = FALSE)
o <- reconcileOSDGeomorph(osd, 'hillpos')
res <- vizHillslopePosition(o$geom, verbose = TRUE)
print(res$fig)
par(mar = c(1, 0, 1, 1))
idx <- match(hydOrder(o$geom, g = 'hillpos', clust = FALSE), profile_id(o$SPC))
plotSPC(o$SPC, plot.order = idx, width = 0.35, name.style = 'center-center', cex.names = 0.65)