eval-DSM-data.Rmd

---
title: "Coweeta Raster Soil Survey"
author: "D.E. Beaudette"
date: "`r Sys.Date()`"
output:
  html_document:
    mathjax: null
    jquery: null
    smart: no
    number_sections: yes
    toc: yes
    toc_depth: 3
    toc_float:
      collapsed: true
      smooth_scroll: false
---

```{r setup, echo=FALSE, results='hide', warning=FALSE}
library(knitr, quietly=TRUE)

opts_chunk$set(
  message = FALSE, 
  warning = FALSE, 
  background = '#F7F7F7', 
  fig.align = 'center', 
  dev = 'png', 
  fig.retina = 2,
  tidy = FALSE, 
  verbose = FALSE
)

options(width=100, stringsAsFactors=FALSE)
```

```{r}
library(soilDB)
library(aqp)
library(sharpshootR)
library(terra)
library(rasterVis)
library(sf)
library(stringi)
library(knitr)
library(kableExtra)


# RSS: cell values are map unit keys
# has a simple raster attribute table
r <- rast('grids/rss_utm.tif')

# Coweeta watersheds
b <- read_sf('vect/Coweeta_Hydrologic_Laboratory.shp')

# points
p <- read_sf('vect/Soil_Moisture_Sites_CW.shp')

# reduce to single polygon
x <- st_union(b)

# get intersecting SSURGO linework as SpatialPolygonsDataFrame from SDA
# this uses the bounding-box
mu.poly <- SDA_spatialQuery(x, what = 'mupolygon', geomIntersection = TRUE)
mu.poly <- st_transform(mu.poly, st_crs(crs(r)))

# transform to common CRS
b <- st_transform(b, st_crs(r))
p <- st_transform(p, st_crs(r))
```


Load pre-processed tabular data.
```{r}
load('data/rss-tab-data-raw.rda')
```


Develop thematic data
```{r}

## TODO / BUD: 0-values added to RAT
# grids loaded via raster() have a 0-entry in the RAT


r.mukey <- raster('grids/rss_mukey.tif')
r.coname <- raster('grids/rss_coname.tif')
r.mukind <- raster('grids/rss_mukind.tif')
r.taxpartsize <- raster('grids/rss_taxpartsize.tif')
r.aws050 <- raster('grids/rss_aws050.tif')
```

```{r}


```

```{r}


```






```{r}

# rastervis requires sp objects
mu.sp <- as(mu.poly, 'Spatial')
```

```{r fig.width=8, fig.height=10}
levelplot(
  r.mukind, 
  att = 'category',
  margin = FALSE, 
  col.regions = viridis, 
  main = 'Map Unit Kind', 
  scales = list(draw = FALSE),
  panel = function(...) {
    panel.levelplot(...)
    sp.lines(as(b, 'Spatial'), col = 'white', lwd = 2)
  }
)

levelplot(
  r.coname, 
  att = 'co.names',
  margin = FALSE, 
  col.regions = viridis, 
  main = 'RSS Components', 
  scales = list(draw = FALSE),
  panel = function(...) {
    panel.levelplot(...)
    sp.lines(as(b, 'Spatial'), col = 'white', lwd = 2)
  }
)

levelplot(
  r.taxpartsize, 
  att = 'taxpartsize',
  margin = FALSE, 
  col.regions = viridis, 
  main = 'Particle Size Class\n(Dominant Condition)', 
  scales = list(draw = FALSE),
  panel = function(...) {
    panel.levelplot(...)
    sp.lines(as(b, 'Spatial'), col = 'black', lwd = 2)
  }
)

levelplot(
  r.aws050, 
  margin = FALSE, 
  col.regions = viridis, 
  main = 'PAWS 0-50cm (cm)\n(Dominant Value)', 
  scales = list(draw = FALSE),
  panel = function(...) {
    panel.levelplot(...)
    sp.lines(as(b, 'Spatial'), col = 'white', lwd = 2)
  }
)



levelplot(
  r.coname, 
  att = 'co.names',
  margin = FALSE, 
  col.regions = viridis, 
  main = 'RSS + SSURGO Linework', 
  sub = 'SSURGO (white lines) and Rock outcrop-Cleveland complex (red lines)',
  scales = list(draw = FALSE),
  panel = function(...) {
    panel.levelplot(...)
    # all SSURGO polygons in the BBOX
    sp.lines(mu.sp, col = 'white')
    # Rock outcrop-Cleveland complex, 30 to 95 percent slopes
    sp.lines(mu.sp[which(mu.sp$mukey == '545857'), ], col = 'firebrick')
  }
)
```


```{r}


## approximate acreage
## pixel-based estimates of area (largest, non-misc. area)
mu.area <- na.omit(as.data.frame(freq(r.mukey)))
mu.area$ac <- mu.area$count * (mean(res(r.mukey))^2) * 0.000247105

a <- merge(mu.area, rss.co[rss.co$compkind == 'Series', ], by.x = 'value', by.y = 'mukey', sort = FALSE, all.x = TRUE)
a$comp.ac <- a$ac * (a$comppct_r / 100)

co.area <- aggregate(comp.ac ~ compname, data = a, FUN = sum)
co.area <- co.area[order(co.area$comp.ac, decreasing = TRUE), ]

# proportion
co.area$proportion <- co.area$comp.ac / sum(co.area$comp.ac)

# investigate named components
osds <- fetchOSD(co.area$compname)

# encode horizon boundary distinctness
osds$hzd <- hzDistinctnessCodeToOffset(osds$distinctness)

# # hmm...
# kable_styling(kable(co.area, row.names = FALSE, digits = 2), full_width = FALSE, font_size = 10)
```

```{r}
# join to SPC
co.area$id <- toupper(co.area$compname)
site(osds) <- co.area[, c('id', 'comp.ac', 'proportion')]

s <- site(osds)
kable_styling(kable(s[, c('id', 'family', 'drainagecl', 'proportion')], digits = 2), full_width = FALSE, font_size = 10)
```


```{r fig.width=10, fig.height=6}
SoilTaxonomyDendrogram(osds, width = 0.35, name.style = 'center-center', y.offset = 0.4, hz.distinctness.offset = 'hzd', shrink = TRUE)


# get plotting details from aqp environment
lsp <- get('last_spc_plot', envir = aqp.env)

# sensor depths
.depths <- (c(10, 20, 50, 80, 100) * lsp$scaling.factor) + lsp$y.offset[1]

abline(h = .depths, lty = 2, col = 'green')
```

```{r fig.width=10, fig.height=6}
o.sub <- osds

idx <- order(o.sub$comp.ac, decreasing = TRUE)

ragg::agg_png(filename = 'RSS-OSDs.png', width = 1800, height = 900, scaling = 1.5)

par(mar = c(0, 0, 0, 0.5))


plotSPC(o.sub, width = 0.35, name.style = 'center-center', hz.distinctness.offset = 'hzd', shrink = TRUE, hz.depths = TRUE, fixLabelCollisions = TRUE, hz.depths.offset = 0.05, id.style = 'top', cex.names = 0.75, plot.order = idx, cex.id = 0.66, max.depth = 225, axis.line.offset = -4, cex.depth.axis = 1, n.depth.ticks = 10)

# planned depths for soil sensors
.depths <- c(10, 20, 50, 80, 100)
segments(x0 = 0.5, y0 = .depths, x1 = length(o.sub) + 1, y1 = .depths, lty = 2, col = 'blue', lwd = 1.5)

axis(side = 1, at = 1:length(o.sub), labels = round(o.sub$proportion[idx], 2), line = -5, cex.axis = 1)

mtext('Approximate Area Proportion within Coweeta', side = 1, line = -2.5, at = 0.5, adj = 0, font = 2)

dev.off()

```