Merge pull request #7 from ncss-tech/crisp-expression1

Improvements to rule and crisp expression parsing

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: InterpretationEngine
 Title: Stand-alone NASIS Interpretation Engine
-Version: 0.1.0
+Version: 0.1.1
 Authors@R: 
     c(person(given = "Dylan",
            family = "Beaudette",

NAMESPACE

@@ -19,6 +19,7 @@ export(extractCrispCurveEval)
 export(extractCrispExpression)
 export(extractEvalCurve)
 export(extractGaussCurveEval)
+export(extractIsNull)
 export(extractLinearCurveEval)
 export(extractPICurveEval)
 export(extractSigmoidCurveEval)

R/EvaluationCurves.R

@@ -132,11 +132,16 @@ extractEvalCurve <- function(evalrec, xlim = NULL, resolution = NULL, sig.scale
     return(res)
   }
 
-  ## ... there are others
-  #   IsNull -- not needed? / not a curve?
+  if (et == "IsNull") {
+    res <- extractIsNull(invert = invert.eval)
+    return(res)
+  }
+
+  warning("extractEvalCurve: curve type (", et, ") not supported", call. = FALSE)
 
-  warning("extractEvalCurve: curve type not yet supported", call. = FALSE)
-  return(function(evalrec) {return(NULL)})
+  return(function(evalrec) {
+    return(NULL)
+  })
 }
 
 #' Function Generators for Interpolating Evaluation Curves
@@ -238,6 +243,26 @@ extractCrispExpression <- function(x, invert = FALSE, asString = FALSE) {
 }
 
 
+#' Extract IsNull Evaluation Logic as R function
+#' 
+#' Default behavior of IsNull evaluation returns value > `0` (1) if `NULL`, inverted behavior returns `0` if `NULL`
+#' 
+#' @param invert invert logic with `!`? Default: `FALSE`
+#'
+#' @return a generated function of an input variable `x` 
+#' 
+#' @details The generated function returns a logical value (converted to numeric) when the relevant property data are supplied.
+#' 
+#' @export
+extractIsNull <- function(invert = FALSE) {
+  if (invert) {
+    function(x) .NULL_HEDGE(x, null.value = 1)
+  } else {
+    function(x) .NULL_HEDGE(x, null.value = 0)
+  }
+}
+
+
 # INTERNAL METHODS ----
 
 # internal method for approxfun() linear interpolation
@@ -260,8 +285,8 @@ extractCrispExpression <- function(x, invert = FALSE, asString = FALSE) {
     crisp.expression <- paste0('domain ', crisp.expression)
 
     # replace logical operators, and add domain vector
-    crisp.expression <- gsub('and', '& domain', crisp.expression)
-    crisp.expression <- gsub('or', '| domain', crisp.expression)
+    crisp.expression <- gsub('and', '& domain', crisp.expression, ignore.case = TRUE)
+    crisp.expression <- gsub('or', '| domain', crisp.expression, ignore.case = TRUE)
   }
 
   if (is.null(xlim)) {
@@ -546,36 +571,84 @@ extractCrispExpression <- function(x, invert = FALSE, asString = FALSE) {
 
 # regex based property crispexpression parser (naive but it works)
 .crispFunctionGenerator <- function(x, invert = FALSE, asString = FALSE) {
+
+  # numeric constants get a short circuit
+  if (grepl("^\\d+$", x)) {
+    res <- sprintf("function(x) { %s }", x)
+    if (asString) return(res)
+    res <- try(eval(parse(text = res)))
+    return(res)
+  }
+
+  # remove empty quotations (some expressions have these trailing with no content)\
+  step0 <- gsub("or +or", "or", gsub("\t", " ", gsub("\"\"", "", gsub("'", "\"", x))))
+
+  if (grepl("[^(]*\\)$", step0)) {
+    step0 <- gsub(")$", "", step0)
+  }
+
+  if (grepl(" TO ", step0, ignore.case = TRUE)) {
+    step0 <- gsub("(.*) TO (.*)", "x >= \\1 & x <= \\2", step0, ignore.case = TRUE)
+  }
+
+  if (grepl("^i?matches", step0, ignore.case = TRUE)) {
+    step0.5 <- gsub("\" *or *i?matches *\"|\" *or *\"|\", \"", "$|^", step0, ignore.case = TRUE)
+  } else {
+    step0.5 <- gsub("\" *or *i?matches *\"|\", \"", "$|^", step0, ignore.case = TRUE)
+  }  
+
   # wildcards matches/imatches
-  step1 <- gsub("i?matches \"([^\"]*)\"", "grepl(\"^\\1$\", x, ignore.case = TRUE)", 
-                gsub("\" or i?matches \"", "$|^", x, ignore.case = TRUE), ignore.case = TRUE)
+  step1 <- gsub(
+    "i?matches +\"([^\"]*)\"",
+    "grepl(\"^\\1$\", x, ignore.case = TRUE)",
+    step0.5,
+    ignore.case = TRUE
+  )
   step2 <- gsub("*", ".*", step1, fixed = TRUE)
 
   # (in)equality  
   step3 <- gsub(" x  grepl", "grepl", gsub("^([><=]*) ?(\")?|(and|or) ([><=]*)? ?(\")?", "\\3 x \\1\\4 \\2\\5", step2))
+  step3 <- gsub("x +x", "x", step3)
+  step3 <- gsub("x  \"", "x == \"", step3)
 
   # convert = to ==
-  step4 <- gsub("x =? ", "x == ", gsub("\" ?(, ?| or ?)\"", "\" | x == \"", step3, ignore.case = TRUE))
+  step4 <- gsub("x [^<>]=? ", "x == ", 
+                gsub("\" ?(, ?| or ?)\"", "\" | x == \"", 
+                     step3, ignore.case = TRUE))
+
+  # convert partial matches to grepl
+  step5 <- gsub("x == +(\"[^\"]*\\.\\*[^\"]*\")", "grepl(\\1, x)", step4, ignore.case = TRUE)
 
   # convert and/or to &/|
-  expr <- trimws(gsub(" or ", " | ", gsub(" and ", " & ", step4)))
+  expr <- trimws(gsub("([^no])or *x?", "\\1 | x ", gsub(" *and *x?", " & x ", step5, ignore.case = TRUE), ignore.case = TRUE))
 
   # various !=
-  expr <- gsub("== != \"|== not \"", "!= \"", expr, ignore.case = TRUE)
+  expr <- gsub("== != *\"|== not *\"", "!= \"", expr, ignore.case = TRUE)
+  expr <- gsub("== !=", "!= ", expr, ignore.case = TRUE)
   expr <- gsub("== \"any class other than ", "!= \"", expr)
 
   # final matches
-  expr <- gsub("== MATCHES ", "== ", expr, ignore.case = TRUE)
+  expr <- gsub("== =", "==", gsub("== MATCHES ", "== ", expr, ignore.case = TRUE))
+
+  # grepl
+  expr <- gsub("x grepl", "grepl", expr, ignore.case = TRUE)
+
+  # not grepl
+  expr <- gsub("x =* *not grepl", "!grepl", expr, ignore.case = TRUE)
 
   # many evals just return the property
   expr[expr == "x =="] <- "x"
 
+  expr <- gsub("x +NOT", "x !=", expr)
+
   # logical expression, possibly inverted, then converted to numeric (0/1)
   # TODO: handle NA via na.rm/na.omit, returning attribute of offending indices
   res <- sprintf("function(x) { as.numeric(%s(%s)) }", 
                  ifelse(invert, "!", ""), expr)
   if (asString) return(res)
-  res <- eval(parse(text = res))
+  res <- try(eval(parse(text = res)))
+  # if (inherits(res, "try-error"))
+  #   browser()
   attr(res, 'CrispExpression') <- x
   res
 }

R/local-functions.R

@@ -119,20 +119,20 @@ getAndCacheData <- function() {
   # get rules, note that "rule" is a reserved word, use [] to protect
   # load ALL rules, even those not ready for use
   rules <- soilDB::dbQueryNASIS(soilDB::NASIS(), 
-                                "SELECT rulename, ruledesign, primaryinterp, notratedphrase, 
+                                "SELECT rulename, ruledesign, primaryinterp, notratedphrase, dataafuse, 
                                         ruledbiidref, ruleiid, [rule]
                                  FROM rule_View_0")
 
   # get all evaluation curves
   evals <- soilDB::dbQueryNASIS(soilDB::NASIS(), 
-                                "SELECT evaliid, evalname, evaldesc, CAST(eval AS text) AS eval,
+                                "SELECT evaliid, evalname, evaldesc, CAST(eval AS text) AS eval, dataafuse,
                                         evaluationtype, invertevaluationresults, propiidref AS propiid
                                  FROM evaluation_View_0")
 
   # get basic property parameters, but not the property definition
   properties <- soilDB::dbQueryNASIS(soilDB::NASIS(), 
                                      "SELECT propiid, propuom, propmin, propmax, propmod, propdefval,
-                                             propname 
+                                             propname, dataafuse
                                       FROM property_View_0")
 
   # property descriptions and CVIR code

misc/dwb-comparison.R

@@ -0,0 +1,54 @@
+library(terra)
+
+projectcode <- "UTGrand250"
+
+paths.dwb <- list.files(path = paste0("inst/extdata/Input/DWB/NASIS Property Rasters/", projectcode), pattern = "tif$", full.names = T)
+
+# all output will be in this crs (5070 = albers AEA conus)
+crs <- crs("EPSG:5070")
+
+sf.studyarea <- vect("demo//Utah_County_Boundaries-shp", "Counties") |>
+  subset(NAME == "GRAND", NSE = TRUE)
+
+paths.dwb <- list.files(path = paste0("inst/extdata/Input/DWB/NASIS Property Rasters/", projectcode), pattern = "tif$", full.names = T)
+
+rast.dwb <- rast(paths.dwb)
+names.dwb <- basename(paths.dwb) |>  
+  gsub(pattern = paste0(projectcode, "_"), replacement = "") |> 
+  gsub(pattern = ".tif", replacement = "")
+names(rast.dwb) <- names.dwb
+
+df.dwb <- as.data.frame(rast.dwb)
+t1 <- Sys.time()
+df.dwb$dwb <- dwb_calc(df.dwb)
+timeDWB <- Sys.time() - t1
+timeDWB
+
+r <- initRuleset("ENG - Dwellings With Basements")
+newnames <- c(
+  "DEPTH.TO.BEDROCK.HARD..BELOW.O.HORIZON",
+  "DEPTH.TO.BEDROCK.SOFT..BELOW.O.HORIZON",
+  "DEPTH.TO.CEMENTED.PAN.THICK..BELOW.O.HORIZON",
+  "DEPTH.TO.CEMENTED.PAN.THIN..BELOW.O.HORIZON",
+  "DRAINAGE.CLASS.IS.NOT.SUBAQUEOUS",
+  "FLOODING.FREQUENCY..Maximum.Frequency.",
+  "FRAGMENTS...75MM.WEIGHTED.AVE..IN.DEPTH.0.100CM",
+  "SUBSIDENCE.DUE.TO.GYPSUM..REV",
+  "COMPONENT.LOCAL.PHASE.IMPACTED",
+  "LEP.WTD_AVG.25.150cm.OR.ABOVE.RESTRICTION..BELOW.O.HORIZON",
+  "USDA.TEXTURE.MODIFIER",
+  "UNIFIED.BOTTOM.LAYER",
+  "DEPTH.TO.PERMAFROST",
+  "USDA.TEXTURE..IN.LIEU.OF.",
+  "PONDING.FREQUENCY",
+  "RESTRICTIVE.FEATURE.HARDNESS",
+  "SLOPE",
+  "SUBSIDENCE.TOTAL",
+  "COMPONENT.LOCAL.PHASE.UNSTABLE.FILL",
+  "HIGH.WATER.TABLE.DEPTH.MINIMUM"
+)
+names(rast.dwb) <- newnames
+rast.dwb$PONDING.DURATION <- rast.dwb$PONDING.FREQUENCY
+df.dwb2 <- as.data.frame(rast.dwb)
+profvis::profvis(res <- interpret(r, rast.dwb, cores = 8))
+plot(res$rating ~ df.dwb$dwb$maxfuzz)

misc/mass-movement-risk.R

@@ -194,8 +194,20 @@ muagg <- comp3[, .SD[which.max(comp_pct)[1], ], by = "mukey"]
 ssurgo_rc <- ssurgo_r
 levels(ssurgo_rc) <- cbind(data.frame(ID = as.numeric(muagg$mukey)), muagg)
 ssurgo_rc <- catalyze(ssurgo_rc)
-plot(ssurgo_rc$Mass_movement_rt, col = hcl.colors(3))
-plot(ssurgo_rc$SLOPE)
+ssurgo_rc_wgs84 <- project(ssurgo_rc, "EPSG:4326")
+
+plot(ssurgo_rc_wgs84$Mass_movement_rt, 
+     col = hcl.colors(20), 
+     range = c(0, 0.8))
+plot(ssurgo_rc_wgs84$rating_new, 
+     col = hcl.colors(20), 
+     range = c(0, 0.8))
+
+plot(ssurgo_rc_wgs84$SLOPE, 
+     col = hcl.colors(20), 
+     range = c(0, 100))
+plot(ssurgo_rc_wgs84$TWO.DIMENSIONAL.SURFACE.MORPHOMETRY, 
+     col = hcl.colors(20))
 cnm <- c("SLOPE", "DEPTH.OF.ISOTROPIC.SOIL", "LIQUID.LIMIT.OF.LAST.LAYER.OR.AT.RESTRICTIVE.LAYER", 
   "CLAY.IN.DIFFERENTIAL.ZONE", "LEP.OF.LAST.LAYER.OR.AT.RESTRICTIVE.LAYER", 
   "LIQUID.LIMIT.IN.DIFFERENTIAL.ZONE", "DEPTH.TO.FIRST.RESTRICTIVE.LAYER..NONE.IS.NULL", 
@@ -220,6 +232,13 @@ rtp <- rast("D:/DebrisFlow/mosquito_rtp.tif")
 test <- c(ssurgo_rc$SLOPE, 
           project(mask(slope, ssurgo_b), ssurgo_rc),
           project(mask(rtp, ssurgo_b), ssurgo_rc))
+test_wgs84 <- project(test, "EPSG:4326")
+
+plot(test_wgs84$slope, 
+     col = hcl.colors(20), 
+     range = c(0, 100))
+plot(test_wgs84$mosquito_rtp, 
+     col = hcl.colors(5))
 # plot(test[[1:2]], range = c(0,100), col = hcl.colors(10))
 # plet(rtp, alpha=0.5, tiles="OpenTopoMap")
 
@@ -231,6 +250,9 @@ rtpk <- k_means(test[[3]], centers = 5)
 ## note that class order depends on seed
 plet(rtpk, alpha=0.5, tiles="OpenTopoMap")
 
+plot(project(rtpk, "EPSG:4326"), 
+     col = hcl.colors(5))
+
 luth <- c("backslope", # 1
           "toeslope",  # 2
           "footslope", # 3 

misc/nccpi.R

@@ -0,0 +1,14 @@
+library(InterpretationEngine)
+
+r <- initRuleset("NCCPI - National Commodity Crop Productivity Index (Ver 3.0)")
+p <- getPropertySet(r)
+p <- merge(p, data.table::rbindlist(propdefByPropname(p$propname)), by = "propiid")
+p2 <- p
+p2$evaluation <- NULL
+View(unique(p2))
+
+cvirrr::cleanCVIR(p$prop[2]) |> 
+  cvirrr::CVIRScript() |> 
+  cvirrr::parseCVIR() |>
+  attr("TSQL") |> 
+  soilDB::dbQueryNASIS(soilDB::NASIS(), q = _)