new source data, converting to sf/terra, still some TODOs

MUspatialSummary.R

@@ -1,7 +1,9 @@
 library(soilDB)
+library(sf)
 library(mapview)
 library(latticeExtra)
 library(tactile)
+library(ragg)
 
 source('local-functions.R')
 
@@ -41,40 +43,47 @@ mu.union <- MUspatialSummary(mukey, method = 'union')
 mu.bbox <- MUspatialSummary(mukey, method = 'bbox')
 
 
-mapview(mu.union, fill=NA, color='firebrick', lwd=2, legend=FALSE)
+mapview(mu.union, col.regions = NA, color = 'firebrick', lwd = 2, legend = FALSE)
 
-mapview(mu.bbox, fill='royalblue', alpha.regions=0.25, color='black', lwd=0.5, legend=FALSE)
+mapview(mu.bbox, col.regions = 'royalblue', alpha.regions = 0.25, color = 'black', lwd = 0.5, legend = FALSE)
 
-mapview(mu.bbox, zcol='fd', alpha.regions=0.25, color='black', lwd=0.5)
+mapview(mu.bbox, zcol = 'fd', alpha.regions = 0.25, color = 'black', lwd = 0.5)
 
 
-# get a collection for comparison
+## get a collection for comparison
 mu.summary <- MUspatialSummary(c('1865918', '456137', '162786', '403308', '2534946'), method='bbox')
 
-d <- mu.summary@data
-d$muname_abbv <- gsub(', ', '\n', d$muname)
-d$muname_abbv <- paste(d$muname_abbv, d$mukind, sep = '\n')
+## TODO: these aren't populated as often as we would like
+# map unit / survey order
 
-png(file='mu-spatial-summary-fd.png', width=900, height=400, res = 90)
+table(mu.summary$invesintens)
+table(mu.summary$muname, mu.summary$invesintens)
 
-bwplot(muname_abbv ~ fd, data=d, 
-       scales=list(alternating=3, y=list(cex=0.75), x=list(tick.number=10)), 
-       xlab='Fractal Dimension (polygon complexity)', box.ratio=0.75,
-       par.settings=tactile.theme,
-       panel=function(...) {
+
+# abbreviate map unit names
+mu.summary$muname_abbv <- gsub(', ', '\n', mu.summary$muname)
+mu.summary$muname_abbv <- paste(mu.summary$muname_abbv, mu.summary$mukind, sep = '\n')
+
+agg_png(filename = 'mu-spatial-summary-fd.png', width = 1200, height = 550, scaling = 1.5)
+
+bwplot(muname_abbv ~ fd, data = mu.summary, 
+       scales = list(alternating = 3, y = list(cex = 0.75), x = list(tick.number = 10)), 
+       xlab = 'Fractal Dimension (polygon complexity)', box.ratio = 0.75,
+       par.settings = tactile.theme,
+       panel = function(...) {
          panel.grid(-1, -1)
          panel.bwplot(...)
        })
 
 dev.off()
 
-png(file='mu-spatial-summary-area.png', width=900, height=400, res = 90)
+agg_png(filename = 'mu-spatial-summary-area.png', width = 1200, height = 550, scaling = 1.5)
 
-bwplot(muname_abbv ~ area_ac, data=d, 
-       scales=list(alternating=3, y=list(cex=0.75), x=list(log=10, tick.number=20)), 
-       xscale.components=xscale.components.log10ticks, 
-       xlab='Delineation Area (ac.)', box.ratio=0.75,
-       par.settings=tactile.theme,
+bwplot(muname_abbv ~ area_ac, data = mu.summary, 
+       scales = list(alternating = 3, y = list(cex = 0.75), x = list(log = 10, tick.number = 20)), 
+       xscale.components = xscale.components.log10ticks, 
+       xlab = 'Delineation Area (ac.)', box.ratio = 0.75,
+       par.settings = tactile.theme,
        panel=function(...) {
          panel.grid(-1, -1)
          panel.bwplot(...)
@@ -84,13 +93,35 @@ dev.off()
 
 
 
-xyplot(fd ~ area_ac, groups=muname_abbv, data=d, 
+bwplot(invesintens ~ fd, data = mu.summary, 
+       scales = list(alternating = 3, y = list(cex = 0.75), x = list(tick.number = 10)), 
+       xlab = 'Fractal Dimension (polygon complexity)', box.ratio = 0.75,
+       par.settings = tactile.theme,
+       panel = function(...) {
+         panel.grid(-1, -1)
+         panel.bwplot(...)
+       })
+
+
+bwplot(invesintens ~ area_ac, data = mu.summary, 
+       scales = list(alternating = 3, y = list(cex = 0.75), x = list(log = 10, tick.number = 20)), 
+       xscale.components = xscale.components.log10ticks, 
+       xlab = 'Delineation Area (ac.)', box.ratio = 0.75,
+       par.settings = tactile.theme,
+       panel=function(...) {
+         panel.grid(-1, -1)
+         panel.bwplot(...)
+       })
+
+
+
+xyplot(fd ~ area_ac, groups=muname_abbv, data = mu.summary, 
        type=c('p', 'r'),
        scales=list(y=list(cex=0.75, tick.number=10), x=list(log=10, tick.number=20)), 
        xscale.components=xscale.components.log10ticks, 
        ylab='Fractal Dimension (polygon complexity)', xlab='Delineation Area (ac.)', 
-       par.settings=tactile.theme(),
-       auto.key=list(columns=4, lines=FALSE, points=TRUE, cex=0.75), 
+       par.settings=tactile.theme(superpose.symbol = list(alpha = 0.5)),
+       auto.key=list(columns=3, lines=TRUE, points=FALSE, cex=0.75), 
        panel=function(...) {
          panel.grid(-1, -1)
          panel.xyplot(...)
@@ -101,18 +132,18 @@ xyplot(fd ~ area_ac, groups=muname_abbv, data=d,
 ## does this make any sense?
 
 
-pr <- princomp(cbind(d$fd, log(d$area_ac)), cor = TRUE)
+pr <- princomp(cbind(mu.summary$fd, log(mu.summary$area_ac)), cor = TRUE)
 pc <- predict(pr)
-d$pc1 <- pc[, 1]
-d$pc2 <- pc[, 2]
+mu.summary$pc1 <- pc[, 1]
+mu.summary$pc2 <- pc[, 2]
 
 
-xyplot(pc1 ~ pc2, groups=muname_abbv, data=d, 
+xyplot(pc1 ~ pc2, groups=muname_abbv, data=mu.summary, 
        scales=list(y=list(cex=0.75, tick.number=10), x=list(tick.number=10)), 
        # xscale.components=xscale.components.log10ticks, 
        ylab='PC2', xlab='PC1', 
-       par.settings=tactile.theme(superpose.symbol=list(col=brewer.pal(4, 'Set1'), pch=16, alpha=0.5, cex=0.5)),
-       auto.key=list(columns=4, lines=FALSE, points=TRUE, cex=0.75), 
+       par.settings=tactile.theme(superpose.symbol=list(pch=16, alpha=0.5, cex=0.5)),
+       auto.key=list(columns=3, lines=FALSE, points=TRUE, cex=0.75), 
        panel=function(...) {
          panel.grid(-1, -1)
          panel.xyplot(...)
@@ -121,10 +152,10 @@ xyplot(pc1 ~ pc2, groups=muname_abbv, data=d,
 
 
 
-densityplot(~ pc1, groups=muname_abbv, data=d, pch=NA, bw=0.3,
+densityplot(~ pc1, groups=muname_abbv, data=mu.summary, pch=NA, bw=0.3,
             scales=list(x=list(tick.number=10)), 
-            par.settings=tactile.theme(superpose.line=list(col=brewer.pal(4, 'Set1'), lwd=2)),
-            auto.key=list(columns=4, lines=TRUE, points=FALSE, cex=0.75), 
+            par.settings=tactile.theme(superpose.line=list(lwd=2)),
+            auto.key=list(columns=3, lines=TRUE, points=FALSE, cex=0.75), 
             panel=function(...) {
               panel.grid(-1, -1)
               panel.densityplot(...)

README.md

@@ -1,5 +1,9 @@
 # Intuitive Explanations of Mapping Scale, Survey Order, and MMU
 
+## Key References:
+  * 
+
+
 ## TODO:
    * review D.G. Rossiter's [SSM Lecture Notes](http://www.css.cornell.edu/faculty/dgr2/_static/files/pdf/SSM_LectureNotes2.pdf) on cartographic definitions / thresholds.
    * class accuracy [Three thoughts on soil class maps: (1) Evaluating classification accuracy (2) Taxonomic vs. geographic soil classes (3) Soil geoforms & phenoforms](http://www.css.cornell.edu/faculty/dgr2/_static/files/pdf/PLSCS_seminar_22Mar2018_Rossiter.pdf)

SSA-complexity/MU-entropy.R

@@ -1,15 +1,27 @@
 library(aqp)
 library(latticeExtra)
-library(viridis)
 library(tactile)
 library(hexbin)
-
+library(soilDB)
+library(sf)
 
 ## TODO:
 
 # 1. include misc. areas
 # 2. 
 
+## order 1 mapping?
+# investigate these more
+.sql <- "SELECT mukey, muname, muacres FROM mapunit where invesintens = 'Order 1';"
+SDA_query(.sql)
+
+s <- fetchSDA_spatial('1147457', by.col = 'mukey')
+wk::as_wkt(st_point_on_surface(s))
+
+# https://casoilresource.lawr.ucdavis.edu/soil_web/list_components.php?mukey=1147457
+# https://casoilresource.lawr.ucdavis.edu/gmap/?loc=46.2365,-91.23553
+
+
 
 equal.prob.H <- function(n, b = 2) {
   p <- rep(1, times = n) / n
@@ -30,23 +42,33 @@ epH <- Vectorize(equal.prob.H)
 
 ## how to incorporate area? H/area, by delineation?
 
-
+# created on SoilWeb server
 x <- read.csv('entropy-by-mukey.csv.gz')
 y <- read.csv('entropy-by-mukey-statsgo.csv.gz')
 
-
-x$mukind <- sprintf("SSURGO\n%s", x$mukind)
-y$mukind <- sprintf("STATSGO\n%s", y$mukind)
+# unique mu kinds
+x$u.mukind <- sprintf("SSURGO\n%s", x$mukind)
+y$u.mukind <- sprintf("STATSGO\n%s", y$mukind)
 
 # stack
-vars <- c('areasymbol', 'projectscale', 'mukind', 'mukey', 'entropy', 'n')
+vars <- c('areasymbol', 'projectscale', 'invesintens', 'mukind', 'u.mukind', 'mukey', 'entropy', 'n')
 g <- make.groups(SSURGO = x[, vars], STATSGO = y[, vars])
 
 row.names(g) <- NULL
 
 head(g)
 
-table(g$which, g$mukind, useNA = 'always')
+# set levels
+g$invesintens <- factor(g$invesintens, levels = c('missing', 'Order 1', 'Order 2', 'Order 3', 'Order 4', 'Order 5'))
+g$projectscale <- factor(g$projectscale)
+
+
+
+## TODO: think some more about these
+xtabs( ~ projectscale + mukind, data = g, subset = which == 'SSURGO')
+xtabs( ~ invesintens + mukind, data = g, subset = which == 'SSURGO')
+
+xtabs( ~ which + mukind, data = g)
 
 summary(g$entropy)
 
@@ -63,11 +85,13 @@ z <- g[g$entropy == 0, ]
 head(z)
 
 
-# remove 0-entropy for now
+## remove 0-entropy for now
 g <- subset(g, subset = entropy > 0)
 
 (scale.tab <- sort(round(prop.table(table(g$projectscale)), 3)))
 
+(order.tab <- sort(round(prop.table(table(g$invesintens)), 3)))
+
 
 # equal-prob H
 g$H.max <- epH(g$n)
@@ -85,7 +109,7 @@ g[idx, ]$entropy - g[idx, ]$H.max
 tps <- tactile.theme(plot.symbol = list(cex = 0.5))
 
 
-histogram( ~ entropy | which, data = g, par.settings = tps, breaks = 100, scales = list(alternating = 1, x = list(tick.number = 10), y = list(alternating = 3)), xlab = 'Shannon Entropy (base 2)')
+  histogram( ~ entropy | which, data = g, par.settings = tps, breaks = 100, scales = list(alternating = 1, x = list(tick.number = 10), y = list(alternating = 3)), xlab = 'Shannon Entropy (base 2)')
 
 histogram( ~ entropy / H.max | which, data = g, par.settings = tps, breaks = 100, scales = list(alternating = 1, x = list(tick.number = 10), y = list(alternating = 3)), xlab = 'Shannon Entropy / Equal-Probability Entropy')
 
@@ -95,7 +119,7 @@ bwplot(which ~ entropy, data = g, par.settings = tps, xlab = 'Shannon Entropy (b
 
 
 
-p1 <- bwplot(mukind ~ entropy, 
+p1 <- bwplot(u.mukind ~ entropy, 
              data = g, 
              par.settings = tps, 
              xlab = 'Shannon Entropy (base 2)', 
@@ -109,7 +133,7 @@ p1 <- bwplot(mukind ~ entropy,
 )
 
 
-p2 <- bwplot(mukind ~ entropy / H.max, 
+p2 <- bwplot(u.mukind ~ entropy / H.max, 
              data = g, 
              par.settings = tps, 
              xlab = 'Shannon Entropy / Equal-Probability Entropy', 
@@ -131,6 +155,26 @@ print(p2, split = c(1, 2, 1, 2), more = FALSE)
 dev.off()
 
 
+
+## interesting, mu-level mapping intensity
+
+bwplot(invesintens ~ entropy, 
+       data = g, 
+       par.settings = tps, 
+       xlab = 'Shannon Entropy (base 2)', 
+       main = '', 
+       scales = list(x = list(tick.number = 10)), 
+       varwidth = FALSE, 
+       panel = function(...) {
+         panel.grid(h = 0, v = -1)
+         panel.bwplot(...)
+       }
+)
+
+
+
+
+
 g.sub <- subset(g, subset = areasymbol %in% c('tx027', 'ca630', 'ca113', 'ca792', 'US'))
 g.sub$areasymbol <- factor(g.sub$areasymbol, levels = c('tx027', 'ca630', 'ca113', 'ca792', 'US'))
 
@@ -171,6 +215,19 @@ ggplot(g, aes(x = entropy, y = mukind)) +
   labs(title = 'All Survey Areas FY23', color = 'Interval')
 
 
+ggplot(g, aes(x = entropy, y = invesintens)) +
+  stat_interval(inherit.aes = TRUE, orientation = 'horizontal', size = 5) + 
+  theme_minimal() +
+  theme(legend.position = c(1, 1), legend.justification ='right', legend.direction	
+        = 'horizontal', legend.background = element_rect(fill = 'white', color = NA), axis.text.y = element_text(size = 10, face = 'bold')) + 
+  stat_summary(geom = 'point', fun = median, shape = 21, fill = 'black', col = 'white', cex = 3) +
+  scale_color_brewer(palette = 'Blues') + 
+  scale_x_continuous(n.breaks = 10) +
+  xlab('Shannon Entropy (base 2)\nSingle Component and Misc. Area MU Removed') + ylab('') +
+  labs(title = 'All Survey Areas FY23', color = 'Interval')
+
+
+
 
 ggplot(g, aes(x = n, y = mukind)) +
   stat_interval(inherit.aes = TRUE, orientation = 'horizontal', size = 5) + 
@@ -184,6 +241,19 @@ ggplot(g, aes(x = n, y = mukind)) +
   labs(title = 'All Survey Areas FY23', color = 'Interval')
 
 
+ggplot(g, aes(x = n, y = invesintens)) +
+  stat_interval(inherit.aes = TRUE, orientation = 'horizontal', size = 5) + 
+  theme_minimal() +
+  theme(legend.position = c(1, 1), legend.justification ='right', legend.direction	
+        = 'horizontal', legend.background = element_rect(fill = 'white', color = NA), axis.text.y = element_text(size = 10, face = 'bold')) + 
+  stat_summary(geom = 'point', fun = median, shape = 21, fill = 'black', col = 'white', cex = 3) +
+  scale_color_brewer(palette = 'Blues') + 
+  scale_x_continuous(n.breaks = 10) +
+  xlab('Number of Components\nSingle Component and Misc. Area MU Removed') + ylab('') +
+  labs(title = 'All Survey Areas FY23', color = 'Interval')
+
+
+
 
 
 #

SSA-complexity/combined-analysis.R

@@ -0,0 +1,41 @@
+library(latticeExtra)
+library(hexbin)
+library(tactile)
+library(terra)
+library(soilDB)
+
+
+x <- read.csv('fractal-dimension-test.csv.gz')
+x$invesintens <- factor(x$invesintens, levels = c('missing', 'Order 1', 'Order 2', 'Order 3', 'Order 4', 'Order 5'))
+
+y <- read.csv('entropy-by-mukey.csv.gz')
+
+z <- merge(x, y[, c('mukey', 'mukind', 'entropy')], by = 'mukey', all.x = TRUE, sort = FALSE)
+head(z)
+
+z <- z[which(z$entropy > 0), ]
+
+
+# figure elements and style
+.fy <- '2023'
+.title <- sprintf('FY%s SSURGO\n100k Samples', .fy)
+.cp <- hcl.colors(100, 'zissou1')
+.cpf <- colorRampPalette(.cp)
+
+hexbinplot(
+  fd ~ entropy, 
+  data = z, 
+  xbins = 25, 
+  main = .title, 
+  xlab = 'Shannon Entropy (base 2)', 
+  ylab = 'Fractal Dimension', 
+  trans = log, 
+  inv = exp, 
+  asp = 1, 
+  colramp = .cpf, 
+  type = 'g', 
+  subset = fd < 1.9, 
+  scales = list(tick.number = 6, alterntating = 1), 
+  colorkey = FALSE, 
+  par.settings = tactile.theme()
+)