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)