From fd3c2053c295891b8a774554e7c8fa8dd08a435c Mon Sep 17 00:00:00 2001 From: "Andrew G. Brown" Date: Thu, 16 Nov 2023 11:18:42 -0800 Subject: [PATCH] misc cleanup --- misc/NASIS-RV-copedon.R | 140 ------------- misc/NASIS-lite.R | 263 ------------------------ misc/ca750-newpedons.R | 88 -------- misc/dm-basics.R | 96 --------- misc/get_NSM_current.R | 51 ----- misc/get_OSD_JSON-demo.R | 279 -------------------------- misc/nasis-lab-morphology-to-SQLite.R | 40 ---- misc/revdepcheck.R | 2 + misc/rhubcheck.R | 2 + 9 files changed, 4 insertions(+), 957 deletions(-) delete mode 100644 misc/NASIS-RV-copedon.R delete mode 100644 misc/NASIS-lite.R delete mode 100644 misc/ca750-newpedons.R delete mode 100644 misc/dm-basics.R delete mode 100644 misc/get_NSM_current.R delete mode 100644 misc/get_OSD_JSON-demo.R delete mode 100644 misc/nasis-lab-morphology-to-SQLite.R create mode 100644 misc/revdepcheck.R create mode 100644 misc/rhubcheck.R diff --git a/misc/NASIS-RV-copedon.R b/misc/NASIS-RV-copedon.R deleted file mode 100644 index 9ada11852..000000000 --- a/misc/NASIS-RV-copedon.R +++ /dev/null @@ -1,140 +0,0 @@ -# make RV copedon snapshot work with fetchNASIS / SQLite backend - -library(soilDB) - -SS <- FALSE -nullFragsAreZero <- TRUE -stringsAsFactors <- FALSE - -dsn <- "C:/Geodata/soils/NASIS-data_patched.sqlite" -copedon_rv_data <- "C:/Geodata/soils/copedon-data.rds" -copedon <- readRDS(copedon_rv_data) - -### -### BEGIN -- PATCHES TO THE NASIS-data.sqlite "Pedon Snapshot" for fetchNASIS -### -# res1 <- dbQueryNASIS(soilDB:::.openNASISchannel(dsn), -# "SELECT geomfname, geomfiid, geomftiidref FROM geomorfeat") -# res2 <- dbQueryNASIS(soilDB:::.openNASISchannel(dsn), -# "SELECT geomftname, geomftiid FROM geomorfeattype") -# -# # # get the geomorphic features tables in the SQLite snapshot -# library(DBI) -# con <- DBI::dbConnect(RSQLite::SQLite(), dsn) -# DBI::dbWriteTable(con, "geomorfeat", res1, overwrite = TRUE) -# DBI::dbWriteTable(con, "geomorfeattype", res2, overwrite = TRUE) -# DBI::dbDisconnect(con) - -# # fix the standard lat/long/gps names truncated to 16 characters -# con <- DBI::dbConnect(RSQLite::SQLite(), dsn) -# sitetab <- DBI::dbReadTable(con, "site") -# sitetab$longstddecimaldegrees <- sitetab$longstddecimalde -# sitetab$latstddecimaldegrees <- sitetab$latstddecimaldeg -# sitetab$gpspositionalerror <- sitetab$gpspositionalerr -# DBI::dbWriteTable(con, "site", sitetab, overwrite = TRUE) -# DBI::dbDisconnect(con) - -# # fix the petaxhistory record ID name truncated to 16 characters -# con <- DBI::dbConnect(RSQLite::SQLite(), dsn) -# petaxmtab <- DBI::dbReadTable(con, "petaxhistmoistcl") -# petaxmtab$pedtaxhistoryiidref <- petaxmtab$pedtaxhistoryiid -# DBI::dbWriteTable(con, "petaxhistmoistcl", petaxmtab, overwrite = TRUE) -# petaxotab <- DBI::dbReadTable(con, "petxhistfmother") -# petaxotab$pedtaxhistoryiidref <- petaxotab$pedtaxhistoryiid -# DBI::dbWriteTable(con, "petxhistfmother", petaxotab, overwrite = TRUE) -# DBI::dbDisconnect(con) - -# # fix the different names in MetadataDomainDetail -# con <- DBI::dbConnect(RSQLite::SQLite(), dsn) -# metatab <- DBI::dbReadTable(con, "MetadataDomainDetail") -# metatab$DomainID <- metatab$domain_id -# metatab$ChoiceValue <- metatab$choice_id -# metatab$ChoiceName <- metatab$choice_label -# DBI::dbWriteTable(con, "MetadataDomainDetail", metatab, overwrite = TRUE) -# DBI::dbDisconnect(con) - -### -### END PATCHES TO PEDON SNAPSHOT for fetchNASIS -### - -system.time(f <- fetchNASIS( - SS = SS, - dsn = dsn, - rmHzErrors = FALSE - )) -# TODO: handle uncoding options - -library(aqp) -aqp_df_class(f) <- "data.table" -f <- rebuildSPC(f) -save(f, file = "C:/Geodata/soils/fetchNASIS-data-2.rda") -# load("C:/Geodata/soils/fetchNASIS-data-2.rda") - -system.time(good.ids <- checkHzDepthLogic(f, fast = TRUE)) -save(good.ids, file = "C:/Geodata/soils/fetchNASIS-data-goodids-2.rda") - -f.sub <- subset(f, good.ids$valid) -all(good.ids$valid) -site(f.sub) <- good.ids -all(f.sub$valid) - -## from full set, you can do subset operations on any site level var -mollisols <- subset(f, f$taxorder == "mollisols") - -## here, we match peiid against a lookup table of RV component pedon peiids -f.cp <- subset(f.sub, profile_id(f.sub) %in% unique(copedon$peiid)) - -library(dplyr) - -site(f.cp) %>% - count(taxgrtgroup) %>% - filter(n > 30) - -# calculate the expected MTR for ~27k RV copedons -res <- profileApply(f.cp, mollic.thickness.requirement, clay.attr = 'clay') -f.cp$mollic_thickness_requirement <- res - -## 27790 elements -save(f.cp, file = "C:/Geodata/soils/fetchNASIS-rv-copedon-2.rda") - -copedon_rv_data <- "C:/Geodata/soils/copedon-data.rds" -copedon <- readRDS(copedon_rv_data) -copedon$peiid <- as.character(copedon$peiid) -copedon$phiid <- as.character(copedon$phiid) - -# add mixed moist color information -horizons(f.cp)$phiid <- as.character(horizons(f.cp)$phiid) -horizons(f.cp) <- copedon[,c("phiid","peiid", - "mxhue_moist","mxvalue_moist","mxchroma_moist", - "mxhue_dry","mxvalue_dry","mxchroma_dry")] - -f.cp$m_hue <- f.cp$mxhue_moist -f.cp$m_value <- f.cp$mxvalue_moist -f.cp$m_chroma <- f.cp$mxchroma_moist -f.cp$d_hue <- f.cp$mxhue_dry -f.cp$d_value <- f.cp$mxvalue_dry -f.cp$d_chroma <- f.cp$mxchroma_dry -f.cp$soil_color <- with(horizons(f.cp), munsell2rgb(m_hue, m_value, m_chroma)) -save(f.cp, file = "C:/Geodata/soils/fetchNASIS-rv-copedon-2.rda") - -f.cp$is_mollic_color <- hasDarkColors(f.cp, d_value = NA) - -darkColorInterval <- function(p) { - mss <- getMineralSoilSurfaceDepth(p) - p.sub <- glom(p, mss, estimateSoilDepth(p)) - if (!inherits(p.sub, 'SoilProfileCollection')) { - return(data.frame(peiid = profile_id(p), - mineral_surface = mss, - darkness_depth = NA)) - } - return(data.frame(peiid = profile_id(p), - mineral_surface = mss, - darkness_depth = getSurfaceHorizonDepth(p.sub, pattern = "TRUE", - hzdesgn = "is_mollic_color" - ))) -} -dcdf <- profileApply(f.cp, darkColorInterval, frameify = TRUE) -f.cp$mineral_surface <- NULL -f.cp$darkness_depth <- NULL -site(f.cp) <- dcdf -save(f.cp, file = "C:/Geodata/soils/fetchNASIS-rv-copedon-2.rda") diff --git a/misc/NASIS-lite.R b/misc/NASIS-lite.R deleted file mode 100644 index 6f066e31f..000000000 --- a/misc/NASIS-lite.R +++ /dev/null @@ -1,263 +0,0 @@ -# the following script has been deprecated in favor of -# NASIS-RV-copedon.R -# MollicEpipedonThickness.Rmd (SoilTaxonomy/misc/presentations) - -# devtools::install() - -library(soilDB) - -# SS <- FALSE -# nullFragsAreZero <- TRUE -# stringsAsFactors <- FALSE -# -# nasislite_path <- "C:/Geodata/soils/NASIS-data.sqlite" -# copedon_rv_data <- "C:/Geodata/soils/copedon-data.rds" -# copedon <- readRDS(copedon_rv_data) - -# res1 <- dbQueryNASIS(soilDB:::.openNASISchannel(nasislite_path), -# "SELECT geomfname, geomfiid, geomftiidref FROM geomorfeat") -# res2 <- dbQueryNASIS(soilDB:::.openNASISchannel(nasislite_path), -# "SELECT geomftname, geomftiid FROM geomorfeattype") - -# # quick hack to get the geomorphic features tables in there -# library(DBI) -# con <- DBI::dbConnect(RSQLite::SQLite(), "C:/Geodata/soils/NASIS-data.sqlite") -# DBI::dbWriteTable(con, "geomorfeat", res1, overwrite = TRUE) -# DBI::dbWriteTable(con, "geomorfeattype", res2, overwrite = TRUE) -# DBI::dbDisconnect(con) - -# extended_data <- get_extended_data_from_NASIS_db(SS = SS, -# dsn = nasislite_path, -# nullFragsAreZero = nullFragsAreZero, -# stringsAsFactors = stringsAsFactors) -# f <- fetchNASIS(SS = SS, dsn = nasislite_path) -# -# library(aqp) -# aqp_df_class(f) <- "data.table" -# f <- rebuildSPC(f) -# save(f, file = "C:/Geodata/soils/fetchNASIS-data.rda") -# load("C:/Geodata/soils/fetchNASIS-data.rda") - -# good.ids <- checkHzDepthLogic(f) -# save(good.ids, file = "C:/Geodata/soils/fetchNASIS-data-goodids.rda") - -# f.sub <- subset(f, good.ids$valid) -# all(good.ids$valid) - -## from full set, you can do subset operations on any site level var -# mollisols <- subset(f, f$taxorder == "Mollisols") - -## here, we match peiid against a lookup table of RV component pedon peiids -# f.cp <- subset(f, profile_id(f) %in% unique(copedon$peiid)) - -# library(dplyr) -# -# site(f.cp) %>% -# count(taxgrtgroup) %>% -# filter(n > 30) - -# calculate the expected MTR for ~27k RV copedons -# res <- profileApply(f.cp, mollic.thickness.requirement, clay.attr = 'clay') -# f.cp$mollic_thickness_requirement <- res - -## 27790 elements -# save(f.cp, file = "C:/Geodata/soils/fetchNASIS-rv-copedon.rda") -# -# -# color analysis -# library(aqp) -load("C:/Geodata/soils/fetchNASIS-rv-copedon.rda") -# -# copedon_rv_data <- "C:/Geodata/soils/copedon-data.rds" -# copedon <- readRDS(copedon_rv_data) -# copedon$peiid <- as.character(copedon$peiid) -# copedon$phiid <- as.character(copedon$phiid) -# -# # add mixed moist color information -# horizons(f.cp)$phiid <- as.character(horizons(f.cp)$phiid) -# horizons(f.cp) <- copedon[,c("phiid","peiid", -# "mxhue_moist","mxvalue_moist","mxchroma_moist", -# "mxhue_dry","mxvalue_dry","mxchroma_dry")] - -# f.cp$m_hue <- f.cp$mxhue_moist -# f.cp$m_value <- f.cp$mxvalue_moist -# f.cp$m_chroma <- f.cp$mxchroma_moist -# f.cp$d_hue <- f.cp$mxhue_dry -# f.cp$d_value <- f.cp$mxvalue_dry -# f.cp$d_chroma <- f.cp$mxchroma_dry -# f.cp$soil_color <- with(horizons(f.cp), munsell2rgb(m_hue, m_value, m_chroma)) -# save(f.cp, file = "C:/Geodata/soils/fetchNASIS-rv-copedon.rda") - -# f.cp$is_mollic_color <- hasDarkColors(f.cp, d_value = NA) -# -# darkColorInterval <- function(p) { -# mss <- getMineralSoilSurfaceDepth(p) -# p.sub <- glom(p, mss, estimateSoilDepth(p)) -# if (!inherits(p.sub, 'SoilProfileCollection')) { -# return(data.frame(peiid = profile_id(p), -# mineral_surface = mss, -# darkness_depth = NA)) -# } -# return(data.frame(peiid = profile_id(p), -# mineral_surface = mss, -# darkness_depth = getSurfaceHorizonDepth(p.sub, pattern = "TRUE", hzdesgn = "is_mollic_color" -# ))) -# } -# dcdf <- profileApply(f.cp, darkColorInterval, frameify = TRUE) -# f.cp$mineral_surface <- NULL -# f.cp$darkness_depth <- NULL -# site(f.cp) <- dcdf -# save(f.cp, file = "C:/Geodata/soils/fetchNASIS-rv-copedon.rda") - -load("C:/Geodata/soils/fetchNASIS-rv-copedon.rda") - -library(aqp) -library(sf) -library(ggplot2) -library(rnaturalearth) - -dat <- site(f.cp) - -vplotdat <- subset(dat, dat$mollic_thickness_requirement >= 18) -vioplot::vioplot(vplotdat$mollic_thickness_requirement ~ vplotdat$taxmoistcl) - -do.call('rbind', lapply(split(dat, dat$taxmoistcl), function(x) { - quantile(x$mollic_thickness_requirement, na.rm = TRUE, probs = c(0.4,0.5,0.6,0.7,0.8,0.9)) - })) - -library(dplyr, warn.conflicts = FALSE) -dat <- dat %>% mutate(mtr_group = case_when(mollic_thickness_requirement < 18 ~ "10 cm", - mollic_thickness_requirement == 18 ~ "18 cm", - mollic_thickness_requirement > 18 & - mollic_thickness_requirement < 25 ~ "18 to 25 cm", - mollic_thickness_requirement == 25 ~ "25 cm")) -dsplit <- split(dat, dat$mtr_group) - -datfilt <- dat[!is.na(dat$mtr_group), ] -datfilt$met_mollic_moist <- (datfilt$darkness_depth >= datfilt$mollic_thickness_requirement) -datfilt$met_20cm_moist <- (datfilt$darkness_depth >= 20) -datfilt$met_25cm_moist <- (datfilt$darkness_depth >= 25) -datfilt$met_30cm_moist <- (datfilt$darkness_depth >= 30) - -dat.sp <- datfilt[,c("peiid", "x_std", "y_std", "mtr_group", - "met_mollic_moist", - "met_20cm_moist", "met_25cm_moist", "met_30cm_moist")] - -dat.sp <- st_as_sf(dat.sp[complete.cases(dat.sp),], coords = c("x_std", "y_std"), crs = 4326) - -plot(subset(dat.sp[dat.sp$peiid %in% datfilt$peiid,'met_mollic_moist'], met_mollic_moist == TRUE)) - -world <- ne_countries(scale = "medium", returnclass = "sf") -world$.id <- 1 -world <- merge(world, data.frame(id = 1, mtr_group = unique(datfilt$mtr_group))) - -usa <- subset(world, admin == "United States of America") -dat.sp <- st_crop(dat.sp, usa) - -ggplot(data = datfilt, aes(group = mtr_group)) + - facet_wrap(~ mtr_group, nrow = 2) + - scale_color_viridis_d(direction = -1) + - ggtitle(sprintf('Mollic Epipedon Thickness Requirements in NASIS Component Representative Pedons (n = %s)', - nrow(dat.sp))) + - geom_sf(data = usa, fill = "#93dfb8") + - geom_sf(data = dat.sp[dat.sp$peiid %in% dsplit[[1]]$peiid,], - size = 0.005, aes(color = met_mollic_moist)) + - geom_sf(data = dat.sp[dat.sp$peiid %in% dsplit[[2]]$peiid,], - size = 0.005, aes(color = met_mollic_moist)) + - geom_sf(data = dat.sp[dat.sp$peiid %in% dsplit[[3]]$peiid,], - size = 0.005, aes(color = met_mollic_moist)) + - geom_sf(data = dat.sp[dat.sp$peiid %in% dsplit[[4]]$peiid,], - size = 0.005, aes(color = met_mollic_moist)) + - guides(colour = guide_legend(override.aes = list(size = 4))) + - labs(color = "Meets Moist Color\n") + - coord_sf(crs = st_crs(2163), xlim = c(-2500000, 2500000), ylim = c(-2300000, 730000)) - - -datfilt_wasmollic <- filter(datfilt, met_mollic_moist == TRUE) -dat.sp_wasmollic <- filter(dat.sp, dat.sp$peiid %in% datfilt_wasmollic$peiid) -ggplot(data = datfilt_wasmollic, aes(group = mtr_group)) + - facet_wrap(~ mtr_group, nrow = 2) + - scale_color_viridis_d(direction = -1) + - ggtitle(sprintf('Mollic Epipedon Thickness Requirements in NASIS Component Representative Pedons (n = %s)', - nrow(dat.sp_wasmollic))) + - geom_sf(data = usa, fill = "#93dfb8") + - geom_sf(data = dat.sp_wasmollic[dat.sp_wasmollic$peiid %in% dsplit[[1]]$peiid,], - size = 0.005, aes(color = met_25cm_moist)) + - geom_sf(data = dat.sp_wasmollic[dat.sp_wasmollic$peiid %in% dsplit[[2]]$peiid,], - size = 0.005, aes(color = met_25cm_moist)) + - geom_sf(data = dat.sp_wasmollic[dat.sp_wasmollic$peiid %in% dsplit[[3]]$peiid,], - size = 0.005, aes(color = met_25cm_moist)) + - geom_sf(data = dat.sp_wasmollic[dat.sp_wasmollic$peiid %in% dsplit[[4]]$peiid,], - size = 0.005, aes(color = met_25cm_moist)) + - guides(colour = guide_legend(override.aes = list(size = 4))) + - labs(color = "Meets Moist Color (25cm)\n") + - coord_sf(crs = st_crs(2163), xlim = c(-2500000, 2500000), ylim = c(-2300000, 730000)) - -res <- subset(datfilt_wasmollic, !met_25cm_moist) -sort(table(res$taxgrtgroup)) -to.update <- subset(f.cp, peiid %in% res$peiid) - -taxonnames <- sort(table(to.update$taxonname)) - - -plot(subset(to.update, taxonname %in% names(taxonnames[taxonnames >= 2]))) - -plot_series <- function(series_name) { - print(series_name) - if (series_name != "SND") { - ful <- subset(f.cp, taxonname == series_name) - upd <- subset(to.update, taxonname == series_name) - osd <- try(fetchOSD(series_name)) - if (!is.null(osd)) { - ser <- aqp::combine(osd, subset(ful, !(peiid %in% upd$peiid)), upd) - ser$threshold[ser$darkness_depth >= 25] <- "Meets Requirement" - ser$threshold[ser$darkness_depth < 25] <- "Does Not Meet Requirement" - ser$threshold[is.na(ser$threshold)] <- "OSD" - groupedProfilePlot(ser, groups = "threshold", cex.names = 0.6, max.depth = 200, print.id=FALSE) - abline(h = 25, lty=3, lwd=1, col="red") - title(sprintf("%s (%s)", series_name, osd$family)) - } else { - stop(osd) - } - } -} - -lapply(tail(names(taxonnames[taxonnames >= 2]), 20), plot_series) - -update_mukeys <- unique(copedon[copedon$peiid %in% to.update$peiid, 'nationalmusym']) - -mupoly <- read_sf("E:/Geodata/soils/MLRA_2_SON_FY2021.gdb", "mupolygon") - -mupoly_sub <- subset(mupoly, mupoly$NATMUSYM %in% update_mukeys) -area <- st_area(st_union(mupoly_sub$Shape)) -units(area) <- "acres" -area - -plot(mupoly_sub$Shape) - -mupoly_sub <- merge(mupoly_sub, copedon, all.x=TRUE) -res <- subset(f.cp, peiid %in% mupoly_sub$peiid) -site(res) <- mupoly_sub[,c("compname","comppct_r","peiid")] - -res$compname -res_moll <- subset(res, taxorder == "Mollisols") - -res_moll$mollic_thickness_requirement - res_moll$darkness_depth -(25 - res_moll$darkness_depth) > 0 - -res_moll20 <- subset(res_moll, (20 - res_moll$darkness_depth) > 0) - -plot(res_moll20, label = "compname") -abline(h=18) - - -res_inc <- subset(res, taxorder == "Inceptisols") - -res_inc$mollic_thickness_requirement - res_inc$darkness_depth -(25 - res_inc$darkness_depth) > 0 - -res_inc20 <- subset(res_inc, (20 - res_inc$darkness_depth) > 0) - -plot(res_inc20, label = "compname") -abline(h=18) -res_inc20$taxgrtgroup diff --git a/misc/ca750-newpedons.R b/misc/ca750-newpedons.R deleted file mode 100644 index ddab8fd70..000000000 --- a/misc/ca750-newpedons.R +++ /dev/null @@ -1,88 +0,0 @@ -# new CA750 data -library(soilDB) -library(sharpshootR) -library(sf) -library(raster) - -# use DBI connection to create a NASIS table list object and select some columns -s.sub <- soilDB::createStaticNASIS("site_View_1")[["site_View_1"]][, c( - "usiteid", - "plssmeridian", - "plsstownship", - "plssrange", - "plsssection", - "plsssdetails", - "elev", - "slope", - "aspect" -)] - -s.sub <- subset(s.sub, grepl("^05", usiteid)) - -# rename some columns to match PLSS2LL() data.frame format -colnames(s.sub) <- c("id","m","t","r","s","plsssdetails", - "elev", "slope", "aspect") - -# add meridian section style etc -s.sub$plsssdetails[grepl("meters", s.sub$plsssdetails)] <- "" -s.sub$qq <- ""#sapply(strsplit(s.sub$plsssdetails, ","), function(x) trimws(x[1])) -s.sub$q <- ""#sapply(strsplit(s.sub$plsssdetails, ","), function(x) trimws(x[2])) - -s.sub$m <- "CA21" -s.sub$type <- "SN" - -is_corner <- grepl("corner$|section$|corner of$", s.sub$plsssdetails) -corners <- gsub("(.*) corner$|(.*) section$|(.*) corner of$", "\\1\\2\\3", s.sub$plsssdetails) - -s.complete <- s.sub[complete.cases(s.sub) & - !s.sub$t == "NA" & - !s.sub$r == "NA" & !(s.sub$r == "E") & - !s.sub$s == "NA" & - !s.sub$q == "NA",] - -s.complete$t <- paste0("T", trimws(s.complete$t)) -s.complete$r <- paste0("R", trimws(s.complete$r)) - -spad <- paste0(ifelse(nchar(s.complete$s) == 0, "", "00"), s.complete$s) - -s.complete$s <- substr(spad, nchar(spad)-1, nchar(spad)) -s.complete$plssid <- sharpshootR::formatPLSS(p = s.complete) -s.complete$plssid <- gsub("(.*A)\\d+ meters.*", "\\1", s.complete$plssid) -s.complete <- s.complete[!grepl("NA", s.complete$plssid),] - -# fix(s.complete) -res <- PLSS2LL(s.complete) - -# convert the longlat points from PLSS to sf -pts <- st_as_sf(res[complete.cases(res),], - coords = c("lon", "lat"), - crs = st_crs(4326)) - -# get SSA polygon -bdy <- fetchSDA_spatial("CA750", geom.src = "sapolygon") - -# get raster of mukeys -mukeyras <- mukey.wcs(bdy) - -# get data from SDA -ssurgo <- SDA_query(sprintf( - "SELECT * FROM legend - INNER JOIN mapunit ON legend.lkey = mapunit.lkey - WHERE mukey IN %s", format_SQL_in_statement(unique(values(mukeyras))))) - -# get just target survey area -ssurgo.sub <- subset(ssurgo, areasymbol == "CA750") - -# get the numbers only out -ssurgo.sub$musym <- as.numeric(gsub("(\\d+)|.*","\\1", ssurgo.sub$musym)) - -# new raster attribute table based on numeric musyms used in (most of) CA750 -rat <- merge(levels(mukeyras)[[1]], ssurgo.sub[,c('mukey',"musym")], - by.x = 'ID', by.y = 'mukey', - sort = FALSE, all.x = TRUE, incomparables = NA) -levels(mukeyras) <- rat - -# make a map -plot(mukeyras, "musym") -plot(st_transform(st_as_sf(bdy), crs(mukeyras))$geometry, add = TRUE) -plot(st_transform(pts, crs(mukeyras))$geometry, add = TRUE) diff --git a/misc/dm-basics.R b/misc/dm-basics.R deleted file mode 100644 index 686492364..000000000 --- a/misc/dm-basics.R +++ /dev/null @@ -1,96 +0,0 @@ -# working with a "static" (non-transactional, not authenticated) -# instance of the NASIS database schema/tables with {soilDB}+{dm} -# -# a demonstration of building the project data model -# -# last update: 2020/03/26 -# andrew g. brown -# -# -library(soilDB) -library(dm) -library(daff) -library(dplyr) - -# get projects from selected set using createStaticNASIS -res <- createStaticNASIS(get_NASIS_table_name_by_purpose("project", SS = TRUE)) - -# list project related tables manually -con <- soilDB::NASIS() -x <- DBI::dbListTables(con) -DBI::dbDisconnect(con) - -# use full join manually -projects <- left_join(res[["project_View_1"]], - by = c("projectiid" = "projectiidref"), - res[["projectmapunit_View_1"]]) - -# get lookup table for specified tables -tkey <- get_NASIS_table_key_by_name(names(res)) -fkeys <- tkey$fkey -names(fkeys) <- tkey$table - -# specify important tables that need primary keys -tnprj <- c("project_View_1", "projectlandcatbreakdown_View_1", "projectmilestone_View_1") -tnprj_pkey <- .get_NASIS_pkey_by_name(tnprj) - -# determine child tables of each of the above (will have a foreign key) -project_tables <- subset(tkey, pkeyref == tnprj_pkey[1] & pkeyref != fkey) -lcb_tables <- subset(tkey, pkeyref == tnprj_pkey[2] & pkeyref != fkey) -pmilestone_tables <- subset(tkey, pkeyref == tnprj_pkey[3] & pkeyref != fkey) - -# use dm to set primary keys -resdm <- as_dm(res) %>% - dm_add_pk(!!tnprj[1], !!tnprj_pkey[1]) %>% - dm_add_pk(!!tnprj[2], !!tnprj_pkey[2]) %>% - dm_add_pk(!!tnprj[3], !!tnprj_pkey[3]) - -# use dm to set foreign keys for project table -# e.g projectstaff, projectmapunit, projectmilestone are children of project -fkey.sub <- fkeys[which(fkeys == unique(project_tables$fkey))] -for(f in seq_along(fkey.sub)) { - resdm <- resdm %>% - dm_add_fk(!!names(fkey.sub)[f], !!fkey.sub[f], ref_table = !!tnprj[1]) -} - -# set projectmappingprogress child of projectlandcatbreakdown -fkey.sub <- fkeys[which(fkeys == unique(lcb_tables$fkey))] -for(f in seq_along(fkey.sub)) { - resdm <- resdm %>% - dm_add_fk(!!names(fkey.sub)[f], !!fkey.sub[f], ref_table = !!tnprj[2]) -} - -# set projectmilestoneprogress child of projectmilestone -fkey.sub <- fkeys[which(fkeys == unique(pmilestone_tables$fkey))] -for(f in seq_along(fkey.sub)) { - resdm <- resdm %>% - dm_add_fk(!!names(fkey.sub)[f], !!fkey.sub[f], ref_table = !!tnprj[3]) -} - -# inspect a wimpy schema (but a fine demo) -resdm %>% - dm_draw() - -# try flattening using {dm} and defined PKEY/FKEY relationship -projects_dm <- resdm %>% - dm_flatten_to_tbl(start = "projectmapunit_View_1") - -# compare column names -colnames(projects) -colnames(projects_dm) - -# count number of project mapunits by project (basic join) -test1a <- projects %>% - group_by(uprojectid) %>% - count() %>% - arrange(desc(n)) -test1a - -# count number of project mapunits by project (using {dm}) -test1b <- projects_dm %>% - group_by(uprojectid) %>% - count() %>% - arrange(desc(n)) -test1b - -diff_data(test1a, test1b) diff --git a/misc/get_NSM_current.R b/misc/get_NSM_current.R deleted file mode 100644 index 49da2b5b1..000000000 --- a/misc/get_NSM_current.R +++ /dev/null @@ -1,51 +0,0 @@ -#' Download nationalsoilmoisture.com CONUS Soil Moisture Percentile ASCII Grids -#' -#' @return data.frame containing file name root, grid label, source URL and download status -#' @export -#' -#' @examples -#' # library(terra) -#' -#' # options(timeout = 300) -#' # nsm_products <- get_NSM_current() -#' # x <- rast(nsm_products$filename) -#' -#' ## as-downloaded (2022/09/18) missing rows at end of file, fixed manually -#' # plot(x$nldas_05_current) -#' # plot(x$nldas_50_current) -#' # writeRaster(x, "~/Geodata/NationalSoilMoisture.tif", overwrite = TRUE) -#' # file.remove(list.files(pattern = "\\.asc")) -#' -#' # x <- rast("~/Geodata/NationalSoilMoisture.tif") -#' # plot(x) -get_NSM_current <- function() { - x <- read.table( - header = TRUE, - text = 'root label - rk_05 "RK 5cm" - rk_20 "RK 20cm" - rk_50 "RK 50cm" - last7days "Precipitation last seven days" - nldas_05 "NLDAS Soil Moisture Layer 1" - nldas_20 "NLDAS Soil Moisture Layer 2" - nldas_50 "NLDAS Soil Moisture Layer 3" - smap_05 "SMAP L3" - smap3_05 "SMAP L3 (3-Day Mean)" - blend_05 "All Blend 5cm" - alluncertaintymax_05 "All Blend 5cm Maximum Uncertainty" - alluncertaintymean_05 "All Blend 5cm Mean Uncertainty" - nkblend_05 "NLDAS/RK 5cm Blend" - nkuncertainty_05 "NLDAS/RK 5cm Uncertainty" - nkblend_20 "NLDAS/RK 20cm Blend" - nkuncertainty_20 "NLDAS/RK 20cm Uncertainty" - nkblend_50 "NLDAS/RK 50cm Blend" - nkuncertainty_50 "NLDAS/RK 50cm Uncertainty"' - ) - x$filename <- paste0(x$root, "_current.asc") - x$url <- - paste0("http://nationalsoilmoisture.com/", x$filename) - x$downloaded <- try(download.file(x$url, x$filename)) - x -} - -# see also: http://nationalsoilmoisture.com/in_situ_05_current.txt etc. diff --git a/misc/get_OSD_JSON-demo.R b/misc/get_OSD_JSON-demo.R deleted file mode 100644 index 0224351fe..000000000 --- a/misc/get_OSD_JSON-demo.R +++ /dev/null @@ -1,279 +0,0 @@ -# There are many analyses that benefit from having direct access to the contents of OSDs parsed by section. -# -# The following code snippet builds a ~30MB .rda file containing OSD data -- a _data.frame_ with one column per "standard OSD section" and one row per series (n=24379). This is too big to deliver as part of the package but it offers significant capabilities over the SoilWeb Postgres fulltext search (soilDB::OSD_query) when one is interested in the details of the OSD narrative rather than just a vector of series names matching some search. -# -# The parsed OSD table from SoilKnowledgeBase is derived from the set of series specific JSON files in https://github.com/ncss-tech/SoilKnowledgeBase. This dataset gives complementary info to that from the SC database and those derived from series/component names/SSURGO: https://github.com/ncss-tech/SoilTaxonomy/blob/master/inst/extdata/SC-database.csv.gz and https://github.com/ncss-tech/SoilTaxonomy/blob/master/inst/extdata/series-stats.csv.gz. -# -# I want to consider some avenues for making the FULL OSD data readily available and _queryable_. This is currently possible if one has an instance of SoilKnowledgeBase / JSON files installed locally. -# -# An interesting new option is https://phiresky.github.io/blog/2021/hosting-sqlite-databases-on-github-pages/ which allows a simple SQLite database to be hosted (read only) via static GitHub Pages. -# -# library(soilDB) -# -# path <- "../SoilKnowledgeBase/inst/extdata/OSD" -# -# # list all -# series <- gsub("(.*).json|(.*log$)", "\\1", basename(list.files(path, recursive = TRUE))) -# series <- series[nchar(series) > 0] -# -# # this uses a _local_ cloned instance of SKB repo (pretty fast) stored in same parent directory as soilDB -# # - base_url = NULL or missing will use GitHub -# res <- soilDB::get_OSD_JSON(series, base_url = path) -# -# ST_series <- res -# save(ST_series, file = "misc/ST_series.rda") -# -# colnames(ST_series) -# #> [1] "SERIES" "STATUS" -# #> [3] "BYREV" "REVDATE" -# #> [5] "STATES" "OVERVIEW" -# #> [7] "TAXONOMIC.CLASS" "TYPICAL.PEDON" -# #> [9] "TYPE.LOCATION" "RANGE.IN.CHARACTERISTICS" -# #> [11] "COMPETING.SERIES" "GEOGRAPHIC.SETTING" -# #> [13] "GEOGRAPHICALLY.ASSOCIATED.SOILS" "DRAINAGE.AND.PERMEABILITY" -# #> [15] "USE.AND.VEGETATION" "DISTRIBUTION.AND.EXTENT" -# #> [17] "REGIONAL.OFFICE" "ORIGIN" -# #> [19] "REMARKS" -# -# nrow(ST_series) -# #> [1] 24379 - - -library(soilDB) -library(tibble) -library(SoilTaxonomy) - -path <- "../SoilKnowledgeBase/inst/extdata/OSD" - -# list all -series <- gsub("(.*).json|(.*log$)", "\\1", basename(list.files(path, - recursive = TRUE))) -series <- series[nchar(series) > 0] - -# this uses a _local_ cloned instance of SKB repo (pretty fast) -# - base_url = NULL or missing will use GitHub -res <- get_OSD_JSON(series, base_url = path) - -ST_series <- res -save(ST_series, file = "misc/ST_series.rda") - -# 2023 series are "incomplete" in one or more sections -tibble(series = series[match(res$SERIES, toupper(series))], - complete = complete.cases(res)) %>% - subset(!complete) - - -# tabulate "missing" sections -apply(res, 2, function(x) sum(is.na(x))) - - -(tibble(res) %>% - subset(is.na(TAXONOMIC.CLASS)))$SERIES - -(tibble(res) %>% - subset(is.na(TYPICAL.PEDON)))$SERIES - - -# most common states -allstates <- table(unlist(strsplit(res$STATES, ","))) - -# count by state -allstates[order(allstates, decreasing = TRUE)] - -# proportion by state -round(prop.table(allstates[order(allstates, decreasing = TRUE)]), 2) - -# how many albolls? -idx <- grep("albolls$", res$TAXONOMIC.CLASS) -res_albolls <- res[idx,] -zzx <- lapply(res_albolls$GEOGRAPHIC.SETTING, function(x) {cat("\n\n"); cat(x)}) - -# examine the geographic setting of albolls -sum(sapply(res_albolls$GEOGRAPHIC.SETTING, function(x) length(grep("lake|lacustrine", x)) > 0)) -sum(sapply(res_albolls$GEOGRAPHIC.SETTING, function(x) length(grep("plain", x)) > 0)) -sum(sapply(res_albolls$GEOGRAPHIC.SETTING, function(x) length(grep("loess", x)) > 0)) -sum(sapply(res_albolls$GEOGRAPHIC.SETTING, function(x) length(grep("till", x)) > 0)) -sum(sapply(res_albolls$GEOGRAPHIC.SETTING, function(x) length(grep("alluvium|stream|terrace", x)) > 0)) - -# mollic epipedon thickness in range in characteristics - -x <- data.table::rbindlist(lapply(seq_along(res$SERIES), function(i) - data.frame( - series = res$SERIES[i], - taxonomy = res$TAXONOMIC.CLASS[i], - ric.content = as.character(strsplit(res$RANGE.IN.CHARACTERISTICS[i], "\\. |\\n")[[1]]) - ))) - -x.mollic <- subset(x, grepl('[^n][^o][^t] ?thick', ric.content, ignore.case = TRUE) & - grepl('mollic epipedon', ric.content, ignore.case = TRUE) | - grepl('Thickness of mollic epipedon[ :\\-]*', ric.content, ignore.case = TRUE)) -x.mollic$units <- NA_character_ -x.mollic$numbers <- gsub(".*[ \\-]+(\\d+ to \\d+) (in|inches|cm|centimeters).*", "\\1;\\2", x.mollic$ric.content) - -x.mollic.split <- strsplit(x.mollic$numbers, " to |;") -names(x.mollic.split) <- x.mollic$series -c("HUNTINGTON","TANEY","TILMA") %in% x$series -taxonfamily <- gsub("TAXONOMIC CLASS: ?", "", x.mollic$taxonomy) -names(taxonfamily) <- x.mollic$series - -x.mollic.result <- data.table::rbindlist(lapply(seq_along(x.mollic.split), function(i) { - if (length(x.mollic.split[[i]]) == 3) { - aSeries <- names(x.mollic.split)[[i]] - data.frame( - series = aSeries, - taxonomy = taxonfamily[[aSeries]], - top = x.mollic.split[[i]][1], - bottom = x.mollic.split[[i]][2], - units = x.mollic.split[[i]][3] - ) - } - })) - -families <- SoilTaxonomy::parse_family(x.mollic.result$taxonomy) -families$series <- x.mollic.result$series - -x.mollic.result <- merge(x.mollic.result, families, all.x=TRUE, sort=FALSE, incomparables=NA) - -inidx <- x.mollic.result[, .I[units %in% c("in","inches")]] - -x.mollic.result$classes_split <- NULL -x.mollic.result$taxonomy <- NULL - -x.mollic.result$topt <- as.numeric(x.mollic.result$top) -x.mollic.result$topt[inidx] <- x.mollic.result$topt[inidx]*2.54 - -x.mollic.result$bottomt <- as.numeric(x.mollic.result$bottom) -x.mollic.result$bottomt[inidx] <- x.mollic.result$bottomt[inidx]*2.54 - -plot(density(round(x.mollic.result$topt), na.rm=T)) -plot(density(round(x.mollic.result$bottomt), na.rm=T)) - -sort(table(round(x.mollic.result$topt))) -sort(table(round(x.mollic.result$bottomt))) - -parenttaxa <- data.table::rbindlist(lapply(SoilTaxonomy::getParentTaxa(code = x.mollic.result$subgroup_code), - function(x) { - if (length(x) == 3) { - data.frame(order = x[1], - suborder = x[2], - greatgroup = x[3]) - } else { - data.frame(order = NA, - suborder = NA, - greatgroup = NA) - } - }), - fill = TRUE) - -x.mollic.result <- cbind(parenttaxa, x.mollic.result) - -min_lt18 <- subset(x.mollic.result, topt < 17.5) -write.csv(min_lt18, file = "misc/osd_mollic_min_thickness_lt18.csv") - -length(unique(min_lt18$greatgroup)) -# [1] 21 - -sort(table(min_lt18$greatgroup), decreasing = TRUE)[1:5] -# top 5 affected: -# Haplustolls Argiustolls Argixerolls Haploxerolls Calcixerolls -# 14 7 4 4 2 - -max_lt25 <- subset(x.mollic.result, bottomt < 24.5) -write.csv(max_lt25, file = "misc/osd_mollic_max_thickness_lt25.csv") -# great groups affected -length(unique(max_lt25$greatgroup)) -# [1] 13 - -sort(table(max_lt25$greatgroup), decreasing = TRUE)[1:5] -# top 5 affected: -# Haplustolls Calciustolls Argixerolls Calciargids Calcixerolls -# 4 2 1 1 1 - -min_is18 <- subset(x.mollic.result, round(topt) == 18) -write.csv(min_is18, file = "misc/osd_mollic_min_thickness_18.csv") -sort(table(min_is18$greatgroup), decreasing = TRUE)[1:5] -# top 5 affected: -# Argiustolls Haploxerolls Haplustolls Calciaquolls Hapludolls -# 73 69 64 43 43 -# - -# looking at other horizon designations and diagnostics -# chunk.idx <- makeChunks(x.mollic.result$series, 50) -# parsedosd <- aqp::combine(lapply(unique(chunk.idx), function(i) fetchOSD(x.mollic.result$series[chunk.idx == i]))) -# save(parsedosd, file = "misc/parsedosds.rda") -load("misc/parsedosds.rda") - -library(aqp) -# 4153/4377 do not have sandy textures through the upper 25 cm -parsedosd025 <- trunc(parsedosd, 0, 25) -parsedosd025$isSandy <- grepl("sand$", as.character(parsedosd025$texture_class)) & - !grepl('very fine', as.character(parsedosd025$texture_class)) -parsedosd025 <- mutate_profile(parsedosd025, any_isSandy = any(isSandy)) -table(parsedosd025$any_isSandy) -sandy_sub <- subset(site(parsedosd025), !any_isSandy) -nrow(sandy_sub) - -# 2753/4377 - have no secondary carbonates -carbonates_sub <- depthOf(parsedosd, pattern = "k", hzdesgn = "hzname") -length(unique(carbonates_sub$id[!complete.cases(carbonates_sub)])) - -# 4316/4377 - no sec. gypsum -gypsum_sub <- depthOf(parsedosd, pattern = "y", hzdesgn = "hzname") -length(unique(gypsum_sub$id[!complete.cases(gypsum_sub)])) - -# 4375/4377 - no fragipan -fragipan_sub <- depthOf(parsedosd, pattern = "x", hzdesgn = "hzname") -length(unique(fragipan_sub$id[!complete.cases(fragipan_sub)])) - -# 3950/4377 - no oxic -oxic_sub <- depthOf(parsedosd, pattern = "o", hzdesgn = "hzname") -length(unique(oxic_sub$id[!complete.cases(oxic_sub)])) - -# 4376/4377 - no spodic -spodic_sub <- depthOf(parsedosd, pattern = "h|s", hzdesgn = "hzname") -length(unique(spodic_sub$id[!complete.cases(spodic_sub)])) - -# 2223/4377 no natric/argillic/kandic -argi_sub <- depthOf(parsedosd, pattern = "t|tn", hzdesgn = "hzname") -length(unique(argi_sub$id[!complete.cases(argi_sub)])) - -# 2894/4377 - no cambic -cambi_sub <- depthOf(parsedosd, pattern = "w|Bg|B2[^t]*", hzdesgn = "hzname") -length(unique(cambi_sub$id[!complete.cases(cambi_sub)])) - -# 2795/4377 - do not have bedrock, duripan, or densic, or petrocalcic -restriction_sub <- depthOf(parsedosd, "Cr|R|Cd|m") -length(unique(restriction_sub$id[!complete.cases(restriction_sub)])) - -# 3976/4377 are not fluv- or cumulic- -fluventic_ids <- site(parsedosd)$id[!grepl("fluv|cumulic", parsedosd$subgroup)] -length(fluventic_ids) - -countdiags <- table(c(sandy_sub$id[!complete.cases(sandy_sub)], - carbonates_sub$id[!complete.cases(carbonates_sub)], - gypsum_sub$id[!complete.cases(gypsum_sub)], - oxic_sub$id[!complete.cases(oxic_sub)], - spodic_sub$id[!complete.cases(spodic_sub)], - argi_sub$id[!complete.cases(argi_sub)], - cambi_sub$id[!complete.cases(cambi_sub)], - restriction_sub$id[!complete.cases(restriction_sub)], - fluventic_ids)) - -nodiags <- subset(x.mollic.result, x.mollic.result$series %in% names(countdiags[countdiags == 9])) -write.csv(nodiags, "mollic_no_diags.csv") - -# series with an OXIC referenced in RIC -res$SERIES[grep("\\b[Oo]xic\\b", res$RANGE.IN.CHARACTERISTICS)] - -idx <- grepl("\\b[Ff]ragipan\\b", res$RANGE.IN.CHARACTERISTICS) & - grepl("olls$", res$TAXONOMIC.CLASS) - -# taxa with FRAGIPAN referenced in RIC that are mollisols -fragitaxa <- SoilTaxonomy::parse_family(res$TAXONOMIC.CLASS[idx]) - -res$RANGE.IN.CHARACTERISTICS[idx] -fragitaxa$SERIES <- res$SERIES[idx] - -fragitaxa[which(getTaxonAtLevel(fragitaxa$subgroup, "order") == "mollisols"),]$SERIES - - diff --git a/misc/nasis-lab-morphology-to-SQLite.R b/misc/nasis-lab-morphology-to-SQLite.R deleted file mode 100644 index ea2452dd6..000000000 --- a/misc/nasis-lab-morphology-to-SQLite.R +++ /dev/null @@ -1,40 +0,0 @@ -# create Lab Morphology SQLite from Geodatabase -# - -library(soilDB) - -exdir <- "D:/Geodata/soils/NASIS_Morphological" -soilDB:::.dumpSSURGOGDB(dsn = "D:/Geodata/soils/NASIS_Morphological_20200116.gdb/NASIS_Morphological_20200116.gdb", - exdir = exdir) - -# throw out the metadata included with the GDB, add lookup tables and area table -for (xn in get_NASIS_table_name_by_purpose(c("metadata","lookup","area"))) { - d <- dbQueryNASIS(NASIS(), paste("SELECT * FROM", xn)) - write.table( - d, - file = file.path(exdir, "tabular", paste0(xn, ".txt")), - sep = "|", - qmethod = "double", - col.names = TRUE, - row.names = FALSE - ) -} - -# Pipe-delimited TXT (with headers) -> SQLite -createSSURGO("D:/Geodata/soils/NASIS_Morphological/nasis_morph.sqlite", - "D:/Geodata/soils/NASIS_Morphological/", header = TRUE) - -# convert names and labels to numeric codes -> standard NASIS static db (storing coded values) -soilDB:::.recode_db("D:/Geodata/soils/NASIS_Morphological/nasis_morph.sqlite") -file.copy("D:/Geodata/soils/NASIS_Morphological/nasis_morph.sqlite", - "D:/Geodata/soils/NASIS_Morphological/nasis_morph_decoded.sqlite") - -# test that it works -f <- fetchNASIS(dsn = "D:/Geodata/soils/NASIS_Morphological/nasis_morph.sqlite", SS = FALSE, mixColors = FALSE) - -# convert numeric codes to names -> a "decoded" NASIS static db (storing user readable values) -soilDB:::.decode_db("D:/Geodata/soils/NASIS_Morphological/nasis_morph_decoded.sqlite") - -# test that it works -options(soilDB.NASIS.skip_uncode = TRUE) -f <- fetchNASIS(dsn = "D:/Geodata/soils/NASIS_Morphological/nasis_morph_decoded.sqlite", SS = FALSE, mixColors = FALSE) diff --git a/misc/revdepcheck.R b/misc/revdepcheck.R new file mode 100644 index 000000000..3ad189ca7 --- /dev/null +++ b/misc/revdepcheck.R @@ -0,0 +1,2 @@ +revdepcheck::revdep_reset() +revdepcheck::revdep_check(num_workers=4) diff --git a/misc/rhubcheck.R b/misc/rhubcheck.R new file mode 100644 index 000000000..46c3e3115 --- /dev/null +++ b/misc/rhubcheck.R @@ -0,0 +1,2 @@ +rhub::check_for_cran() +rhub::check_on_macos()