In brownag/NASIStools: Tools for the USDA-NRCS National Soil Information System (NASIS)

knitr::opts_chunk$set(echo = FALSE)
library(DBI)
library(NASIStools)

tabular_path <- "E:/Geodata/FY2021/2021_final_major"
tabular_file <- paste0(tabular_path, ".sqlite")

# create a SSURGO sqlite database from a snapshot
if (!file.exists(tabular_file)) {
  createSSURGO(dsn = tabular_path, output_path = tabular_file)
}

areasymbol <- c("CA649","CA750","CA731")
uprojectid <- "2021-2SON-MLRA-006"
musym_target <- c("7087","BeF","BeFma")
# musym_target <- c("7091","AnE","BeF","BeFma","TbFma","BrF2")
# musym_target <- c("7091","BrF2")
# musym_target <- c("7092","BrF3")
# compare new mapunits to others
# musym_target <- c("7102","BrF2","BrFma")
# musym_target <- c("7102","BrF3")
# musym_target <- c("7103","BrF2")
# musym_target <- c("7103","BrF3")
# musym_target <- c("7103","BrG2")
# musym_target <- c("8162","JdG2","172")
musym_target <- c("8171","JcD2","JbD2")
# musym_target <- c("8172","JeF2","JcF2","JcE2","JbE2")
# musym_target <- c("8173","JcF2","JdG2","BrF2","JeF2","JbF2","BrF3")

lut <- NASIStools:::musym_to_all(musym_target)

# get mus in project of interest, correlated on legend of interest
pmus <- get_projectmapunits_by_uprojectid(uprojectid, areasymbol)
pmus_new <- subset(pmus, seqnum.2 == 500)
pmus_old <- subset(pmus, seqnum.2 == 2021)
pmus_new <- subset(pmus_new, nationalmusym %in% unique(lut$nationalmusym))
pmus_old <- subset(pmus_old, nationalmusym %in% unique(lut$nationalmusym))

if (nrow(pmus_old) == 0) {
  pmus_old <- subset(pmus, seqnum.2 == 2021)
}

# get whole legend of interest
lmus <- get_lmapunit_by_areasymbol(areasymbol)

# connect to sqlite db
con <- dbConnect(RSQLite::SQLite(), tabular_file)

# query cointerps
cointerp <- get_SSURGO_cointerp(con, close = FALSE, ruledepth = 1)

# get key map
keymap <- get_SSURGO_component_keys(con, close = FALSE)

cat("**SSURGO Tabular Source:** <i>", tabular_path, "</i>\n\n")
cat("**SQLite Database:** <i>", tabular_file, "</i>\n\n")

cat("###", areasymbol, "\n")
cat("#### <b>", uprojectid, "</b>-", unique(pmus$projectname), "\n")

# use SWR unique seqnums to identify old vs new
tablecols <- c("areasymbol","nationalmusym","musym","lmapunitiid","muiid","uprojectid","muname")
print(knitr::kable(pmus_new[,tablecols], row.names = FALSE, caption = "NEW MAPUNIT"))
print(knitr::kable(pmus_old[,tablecols], row.names = FALSE, caption = "REPLACED MAPUNIT(S)"))

# or use legend mapunit additional status
# pmus_old[pmus_old$lmapunitiid %in% lmus[lmus$mustatus == 4,]$lmapunitiid,]
# pmus_new[pmus_new$lmapunitiid %in% lmus[lmus$mustatus == 3,]$lmapunitiid,]

# OR use nonrep DMUs in mapunit

# subset cointerps based on old vs new mapunits
cointerp_new <- subset(cointerp, mukey %in% pmus_new$lmapunitiid)
cointerp_old <- subset(cointerp, mukey %in% pmus_old$lmapunitiid)
rules <- unique(cointerp$mrulename)

cat("##### Number of mapunits (before):", length(unique(cointerp_old$mukey)), "\n")
cat("##### Number of components (before):", length(unique(cointerp_old$cokey)), "\n")
cat("##### Number of mapunits (after):", length(unique(cointerp_new$mukey)), "\n")
cat("##### Number of components (after):", length(unique(cointerp_new$cokey)), "\n")

# cat("##### Number of miscellaneous areas (before)", zzz, "\n")
# cat("##### Number of miscellaneous areas (after)", zzz, "\n")

# area all (proportions of rated) interp classes the same?
x <- lapply(seq_along(rules), function(i) {

  x1 <- subset(cointerp_old, mrulename == rules[i])
  x2 <- subset(cointerp_new, mrulename == rules[i])

  k1 <- subset(keymap, cokey %in% x1$cokey)
  k2 <- subset(keymap, cokey %in% x2$cokey)

  # coarse comparison of proportions in each RV rating class
  before <- table(x1$interphrc)
  before_nr <- table(x1$interphrc, exclude = "Not rated")
  after <-  table(x2$interphrc)
  after_nr <- table(x2$interphrc, exclude = "Not rated")
  res1 <- all.equal(prop.table(before), prop.table(after))
  res2 <- all.equal(prop.table(before_nr),
                    prop.table(after_nr))

  # if no match in coarse interp rating classes
  if (!is.logical(res2)){
    # get the reasons
    reasons <- get_SSURGO_interp_reasons_by_mrulename(con, mrulename = rules[1], n=1,close = FALSE)

    r1 <- subset(reasons, cokey %in% cointerp_new$cokey)
    r2 <- subset(reasons, cokey %in% cointerp_old$cokey)

    newreasons <- table(r1$Reasons)
    oldreasons <- table(r2$Reasons)

    # old reasons that are not in new reasons
    oldnotnew <- oldreasons[!names(oldreasons) %in% names(newreasons)]

    # new reasons that are in old reasons
    newinold <-  newreasons[names(newreasons) %in% names(oldreasons)]

    # new reasons that are not in old reasons
    newnotold <- newreasons[!names(newreasons) %in% names(oldreasons)]

    # get the reasons
    # new <- r1[r1$cokey %in% unique(r1$cokey, names(newnotold)),]
    # old <- r2[r2$cokey %in% unique(r2$cokey, names(newnotold)),]

    cat("#### -----------------------------------------------------------\n")
    cat(rules[i], "\n\n")
    # cat("##### BEFORE\n")
    if (length(before))
      print(knitr::kable(before, 
                       col.names = c("Rating", "# Components"),
                       caption = "Rating Class (Before)"))
    cat("\n")
    # cat("##### AFTER\n")
    if (length(after))
      print(knitr::kable(after, 
                       col.names = c("Rating", "# Components"), 
                       caption = "Rating Class (After)"))
    cat("\n\n###### EXPLANATION\n\n")
    cat("\tNew reasons: ", paste0(names(newnotold), collapse=", "), "\n\n")
    cat("\tExisting reasons: ", paste0(names(newinold), collapse=", "), "\n\n")
    cat("\tDropped reasons: ", ifelse(length(oldnotnew) == 0, "None",
                                    paste0(names(oldnotnew), collapse=", ")), "\n\n")
    .explainerror <- function(x) {
      idx1 <- grep("^Numeric: lengths .* differ$", "Numeric: lengths (2, 1) differ", x)
      x[idx1] <- " - Different number of rating classes in before/after mapunits"
      x
    }
    cat(paste0("\t\t", .explainerror(res1[length(res1)]), "\n"), "\n")
    cat("\tExclude 'Not rated':\n", 
        paste0("\t\t", .explainerror(res2[length(res2)]), "\n"), "\n\n")
    cat("\n")    

    # this should show all new classes, and all not rated
    newrating <- subset(r1[,-6],  !interphrc %in% names(before_nr))
    if (nrow(newrating) > 0){
      print(knitr::kable(newrating, 
                       caption = "Components with New or Not Rated Rating Class",
                       row.names = FALSE))
    } else {
      cat("No new rating classes in 'after' mapunits.\n\n\n")
    }
  }
  res2
})
cat("\n\n--------------\n\n## SUMMARY\n")
cat("<b>", sum(sapply(x, is.logical)), "interpretations are unchanged in bulk rating classes.</b>\n\n")
idx.changed <- which(!sapply(x, is.logical))
cat("<b>",length(idx.changed), "interpretations are changed in bulk rating classes:</b>\n\n")
x <- sapply(rules[idx.changed], function(x) cat(" -", x,"</b>\n\n"))
DBI::dbDisconnect(con)

brownag/NASIStools documentation built on Nov. 21, 2023, 11:34 a.m.