Skip to content
Permalink
Browse files

updates to Rangeland Health soil description summary report

  • Loading branch information
brownag committed Nov 7, 2019
1 parent 1f972d3 commit ad8fdd48925a2b656da9cdacf8e923ec9978f561
@@ -0,0 +1,314 @@
---
title: "Rangeland Health EVAL - Soil Reference Description Generator"
author: "Andrew Brown"
output: html_document
params:
ecositeid:
label: "Ecological Site ID:"
value: "R018XI163CA"
input: text
majorcompswitch:
label: "Major components only?"
value: TRUE
ignorecompname:
label: "Ignore component names:"
value: "Peters,Cometa,Red Bluff,Redding,Yokohl,Rocklin,Sesame,Whitney,Mokelumne,Miltonhills"
input: text
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(knitr)
library(tidyverse)
library(aqp)
library(soilDB)
library(FedData)
library(raster)
library(viridis)
# logic for handling major component switch
flag_majcmp <- 'Yes'
if(!params$majorcompswitch)
flag_majcmp <- 'No'
# logic for ignoring specific component names (e.g. errors in assignment in SSURGO)
to_ignore <- strsplit(params$ignorecompname, ",")
if(!length(to_ignore))
to_ignore <- ''
# parse ecosite (for case of multiple ecosites of interest)
list_ecositeid <- unlist(strsplit(params$ecositeid, ","))
```

```{r fetchdata, message=FALSE}
# construct WHERE clause for fetchSDA_component()
q.where.comp <- paste0("ecoclassid IN ", format_SQL_in_statement(list_ecositeid), " AND
areasymbol != 'US' AND
majcompflag = '", flag_majcmp, "' AND
NOT compname IN ", format_SQL_in_statement(to_ignore))
# get component sets that are 1:1 with nmusym (duplicates=FALSE)
comp_nd <- suppressMessages(fetchNASISWebReport(projectname = "ES - MLRA 18 - Thermic Low Rolling Hills 14 - 20 PZ R018XI163CA", stringsAsFactors = FALSE))
# remove components on ignore list
comp_nd$spc <- comp_nd$spc[which(!comp_nd$spc$compname %in% to_ignore[[1]]),]
# remove additionals, misc areas and minor components
comp_nd$spc <- comp_nd$spc[comp_nd$spc$dmuiid %in% comp_nd$mapunit[comp_nd$mapunit$mustatus == 'correlated',]$dmuiid,]
comp_nd$spc <- comp_nd$spc[comp_nd$spc$compkind != "miscellaneous area",]
comp_nd$spc <- comp_nd$spc[which(comp_nd$spc$majcompflag == '1'),]
# query additional horizon attributes (low and high of key attributes)
# create a horizon-level attribute with component name info (to allow filtering of horizons by compname)
comp_nd$spc$parentname <- aqp::denormalize(comp_nd$spc, 'compname')
# get full list of unique nmusyms (~unique DMUs) and MUKEYs (needed for spatial relates within individual SSAs)
nmusym.list <- unique(comp_nd$mapunit$nationalmusym[comp_nd$mapunit$dmuiid %in% comp_nd$spc$dmuiid])
# do a second query with duplicates=TRUE to get all MUKEYs (to be able to reference full spatial extent)
#mukey.list <- unique(suppressMessages(fetchSDA_component(WHERE = q.where.comp, duplicates = TRUE)$mukey))
mukey.list <- unique(comp_nd$mapunit$lmapunitiid[comp_nd$mapunit$dmuiid %in% comp_nd$spc$dmuiid])
# debug
# print(paste0("Found ", length(mukey.list), " MUKEYs representing ",
# length(nmusym.list), " National Mapunit Symbols with ",
# if(params$majorcompswitch) "major ",
# "components correlated to ", paste0(list_ecositeid, collapse=", ")))
# get legend-level information on components returned above (i.e. meeting ecosite constraint)
muall <- suppressMessages(get_mapunit_from_SDA(WHERE=paste0('nationalmusym IN ', format_SQL_in_statement(nmusym.list))))
legall <- suppressMessages(get_legend_from_SDA(WHERE=paste0('nationalmusym IN ', format_SQL_in_statement(nmusym.list))))
leg <- merge(muall, legall, by=c("areasymbol", "lkey"))
comp <- merge(site(comp_nd$spc), comp_nd$mapunit[,c('dmuiid','lmapunitiid')], "dmuiid")
leg.ex <- merge(comp, leg, by.x = "lmapunitiid", by.y="mukey", all.x = TRUE)
leg.ex <- leg.ex[!is.na(leg.ex$areasymbol),]
# # debug
# print(paste0("Cumulative area correlated to ", paste0(list_ecositeid, collapse=", "),": ",
# round(sum(leg$muacres)), " acres in ",
# length(unique(leg$areasymbol)), " soil survey area(s)."))
```

```{r custom}
# create uniform label for bedrock, and set constant bottomdepth
comp_nd$spc$hzname[grepl(comp_nd$spc$hzname, pattern="R|Cr", ignore.case = TRUE)] <- "R"
comp_nd$spc$hzname[grepl(comp_nd$spc$texture, pattern="br|uwb|wb", ignore.case = TRUE)] <- "R"
# custom massaging for R163 - TODO: implment in include file
comp_nd$spc$hzname[grepl(comp_nd$spc$hzname, pattern="^A|H1", ignore.case = TRUE)] <- "A"
comp_nd$spc$hzname[grepl(comp_nd$spc$hzname, pattern="H[2-9]|B", ignore.case = TRUE)] <- "B"
comp_nd$spc$hzname[grepl(comp_nd$spc$hzname, pattern="^C", ignore.case = TRUE)] <-
"C"
comp_nd$spc$hzname[is.na(comp_nd$spc$texture) & !comp_nd$spc$hzdept_r == 0 & is.na(comp_nd$spc$sandtotal_r)] <- "R"
```

### Ecological Site Acres by Soil Survey Area

```{r acres-by-ssa}
# compute some summaries by soil survey area
res1 <- leg.ex %>%
group_by(areasymbol) %>%
summarise(site_acres = round(sum(muacres * (comppct_r / 100), na.rm=TRUE)),
area_acres = areaacres[1],
site_pct_of_area = round(sum(muacres * (comppct_r / 100) / areaacres * 100, na.rm=TRUE)))
# calculate totals across all SSAs
res1$group <- "All"
res2 <- res1 %>%
group_by(group) %>%
summarise(site_acres = round(sum(site_acres, na.rm=TRUE)),
area_acres = round(sum(area_acres, na.rm=TRUE)),
site_pct_of_area = round((sum(site_acres, na.rm=TRUE) / sum(area_acres, na.rm=TRUE) * 100))) %>%
arrange(desc(site_acres))
# format for table
names(res2)[1] <- "areasymbol"
res3 <- rbind(res1[,-5], res2)
area.lut <- unique(data.frame(leg.ex$areasymbol, leg.ex$areaname))
res3$areaname <- c(as.character(area.lut[match(res3$areasymbol[1:(nrow(res3)-1)], area.lut[,1]),2]), "All")
res3 <- res3[,c(1,5,2,3,4)]
names(res3) <- c("Area Symbol","Area Name","Site Acres","Area Acres", "Site % of Area")
# output table
knitr::kable(res3,
caption = paste0("Acres in Ecological Site [", paste0(list_ecositeid, collapse=", "),
"] tabulated by Soil Survey Area"))
```

### Ecological Site Acres by Component Name

```{r acres-by-component}
# Acres by component name
cmp1 <- leg.ex %>%
group_by(compname) %>%
summarise(comp_acres = round(sum(muacres, na.rm=TRUE)),
n_comp = length(muacres)) %>%
arrange(desc(comp_acres))
# output table
names(cmp1) <- c("Component","Acres", "# of Components")
knitr::kable(cmp1, caption = paste0("Acres in Ecological Site [",
paste0(list_ecositeid, collapse=", "), "] tabulated by Component Name"))
```

```{r fetchspatial}
# spatial plot, colored by national mapunit symbol
chunk_SDA_spatial <- function(mukey.list, attr = c("nationalmusym"), nchunk = 10) {
mukey.chunk <- (1:length(mukey.list) %% nchunk) + 1
s <- NULL
for(i in 1:max(mukey.chunk)) {
idx <- which(mukey.chunk == i)
q <- paste0("SELECT G.MupolygonWktWgs84 as geom, mapunit.mukey", if(nzchar(attr)) ", ",
paste0(attr, collapse=", ")," FROM mapunit
CROSS APPLY SDA_Get_MupolygonWktWgs84_from_Mukey(mapunit.mukey) as G
WHERE mukey IN ", format_SQL_in_statement(mukey.list[idx]))
#NB: FedData has a different (much simpler, but not equivalent) definition of SDA_query
# it also uses the post/rest interface
sp.res.sub <- suppressMessages(soilDB::SDA_query(q))
if(!length(sp.res.sub))
s.sub <- s[0,]
else
s.sub <- soilDB::processSDA_WKT(sp.res.sub)
if(is.null(s)) {
s <- s.sub
} else {
s <- rbind(s, s.sub)
}
}
return(s)
}
s <- chunk_SDA_spatial(mukey.list)
s$ecosite <- paste(list_ecositeid, collapse = "_")
rgdal::writeOGR(s, dsn = ".", layer = paste(list_ecositeid, collapse = "_"), overwrite=TRUE,
driver="ESRI Shapefile")
```

### Ecological Site Spatial Extent Map

```{r spatialplot}
# create color lookup table by nmusym
color.lut <- viridis::viridis(length(unique(leg$nationalmusym)))
names(color.lut) <- unique(leg$nationalmusym)
# create plot with county boundaries
par(mar=c(0.1,0.1,0.1,0.1))
plot(FedData::polygon_from_extent(s), border=0)
caplabels <- unlist(lapply(as.list(maps::map('county', 'ca', names=T, plot=F)),
function(x) tools::toTitleCase(as.character(strsplit(x, ",")[[1]][2]))))
maps::map('county', 'ca', add=T)
plot(s, add=T, col=color.lut[s$nationalmusym], border=color.lut[s$nationalmusym],
main=paste0("Unique National Mapunit Symbols with ",
if(params$majorcompswitch) "major ",
"components\ncorrelated to ", paste0(list_ecositeid, collapse=", ")))
maps::map.text('county', 'ca', add=T, labels=caplabels, cex=1.03, font=3)
text(x=mean(raster::extent(s)[2:1]), y=raster::extent(s)[3]+0.05, cex=1.2, font=2,
paste0("Unique mapunits (National MUSYMs) with ", if(params$majorcompswitch) "major ",
"components\ncorrelated to ", paste0(list_ecositeid, collapse=", ")))
```

### Component Data Summary

```{r component-data}
# depth to restrictions
depth_to_br <- profileApply(comp_nd$spc, top='hzdept_r', bottom='hzdepb_r', estimateSoilDepth, p="Cr|R|Cd|qm")
depth_to_br[!is.finite(depth_to_br) & depth_to_br < 0] <- 0
res <- as.data.frame(quantile(depth_to_br,
probs = c(0,0.05,0.25,0.5,0.75,0.95,1)))
colnames(res) <- c("Depth to Bedrock, cm")
par(mar=c(4.1,2.5,2,2.5))
plot(density(depth_to_br, kernel="epanechnikov", na.rm=TRUE),
main="Density Plot of Depth to Bedrock",
xlab="Depth, cm",
ylim=c(0, 0.035))
points(x=jitter(depth_to_br, 1.5), y=rep(0, length(depth_to_br)), pch=4)
abline(v=res[c(1,7),1], lwd=2, lty=2, col="RED")
abline(v=res[c(2,6),1], lwd=2, lty=2, col="BLUE")
abline(v=res[4,1], lwd=2, col="GREEN")
legend('topright', legend = c("Median","5th/95th Percentile", "Minimum/Maximum", "Component Bedrock Depth"),
lwd=c(2,2,2,NA), lty=c(1,2,2, NA), pch=c(NA,NA,NA,4), col=c("GREEN","BLUE","RED","BLACK"))
kable(t(res), colnames = FALSE)
```

### Component Horizon Data Summary

```{r horizon-data}
# variables to summarize
vars <- c('hzdept_r','hzdepb_r', 'texture',
'claytotal_l','claytotal_r','claytotal_h',
'sandtotal_l','sandtotal_r','sandtotal_h',
'ph1to1h2o_l','ph1to1h2o_r','ph1to1h2o_h',
'ksat_l','ksat_r','ksat_h',
'fragvoltot_l','fragvoltot_r','fragvoltot_h')# needfragvol_l & fragvol_h
# better names, used in final tables / figures
var.names <- c('clay', 'sand', 'ph1to1h2o', 'fragvol')
h.i <- horizons(comp_nd$spc)[,c("chiid","parentname","hzname", vars)]
getmode <- function(v, na.rm=TRUE) {
if(na.rm==TRUE)
v <- v[!is.na(v)]
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}
suppressWarnings(h.i.ranges <- h.i %>%
group_by(parentname, hzname) %>%
summarise(depth_q5 = round(quantile(hzdept_r, probs=c(0.05), na.rm = T)),
middepth = median(((hzdepb_r - hzdept_r) / 2) + hzdept_r),
depth_q95 = round(quantile(hzdepb_r, probs=c(0.95), na.rm = T)),
claylow = round(quantile(claytotal_l, probs=c(0.05), na.rm = T)),
clayrv = round(getmode(claytotal_r), 1),
clayhigh = round(quantile(claytotal_h, probs=c(0.95), na.rm = T)),
pHlow = round(min(ph1to1h2o_l, na.rm = T), 1),
pHrv = round(getmode(ph1to1h2o_r), 1),
pHhigh = round(max(ph1to1h2o_h, na.rm = T), 1),
rflow = round(min(fragvoltot_l, na.rm = T)),
rfrv = round(getmode(fragvoltot_r), 1),
rfhigh = round(max(fragvoltot_h, na.rm = T)),
n = length(claytotal_r)) %>%
filter(hzname != "R"))
suppressWarnings(h.i.classes <- h.i %>%
group_by(parentname, hzname) %>%
summarise(textures = paste0(unique(sort(toupper(texture))), collapse=", "),
n = length(claytotal_r)) %>%
filter(hzname != "R") %>%
na.omit())
h.i.ranges[is.na(h.i.ranges)] <- "-"
h.i.classes[is.na(h.i.classes)] <- "-"
knitr::kable(h.i.ranges,
caption=paste0("Representative Values for Horizon Properties (all components correlated to ", paste0(list_ecositeid, collapse=", "),")"))
knitr::kable(h.i.classes,
caption=paste0("Representative Classes for Horizons (all components correlated to ", paste0(list_ecositeid, collapse=", "),")"))
# h.i.long <- h.i %>%
# gather_(key="variable", value = "value", vars) %>%
# filter(!is.na(hzname) & !is.na(value))
```

### Modern Series Concepts & Taxonomy
```{r soiltax}
osd <- fetchOSD(soils = unique(comp_nd$spc$compname))
plot(osd, name="hzname")
# family-level differentia
kable(cbind(site(osd),getSoilDepthClass(osd, top = "top", bottom = "bottom")[,c('depth','depth.class')])[order(osd$tax_tempcl, osd$subgroup, osd$tax_partsize, osd$drainagecl), c('id','depth.class','drainagecl','tax_partsize','tax_minclass','tax_ceactcl','tax_tempcl','subgroup')], row.names = FALSE)
```
@@ -12,7 +12,7 @@ params:
value: TRUE
ignorecompname:
label: "Ignore component names:"
value: ""
value: "Peters,Cometa"
input: text
---

@@ -182,7 +182,7 @@ chunk_SDA_spatial <- function(mukey.list, attr = c("nationalmusym"), nchunk = 10
s <- chunk_SDA_spatial(mukey.list)
rgdal::writeOGR(s, dsn = ".", layer = paste(list_ecositeid,collapse = "_"), driver="ESRI Shapefile")
rgdal::writeOGR(s, dsn = ".", layer = paste(list_ecositeid,collapse = "_"), driver="ESRI Shapefile", overwrite=TRUE)
```

### Ecological Site Spatial Extent Map

0 comments on commit ad8fdd4

Please sign in to comment.
You can’t perform that action at this time.