Updates to get_SDA_* / SOD-like methods (#185)

* get_SDA_property: Vectorizing `property` argument, add argument to get query string

* get_SDA_interpretation: Vectorizing `rulename`

* get_SDA*: Standardize lower case MUKEY, COKEY, CHKEY, etc.

* get_SDA_interpretations: consistent column names for all methods w/ .cleanRuleColumnName()

* Demo of bottom_depth logic for weighted average, and comparing column names in results

* removing conversion of NULL to 0, why was this in there

* get_SDA_interpretation: add not_rated_value argument, default NA_real_

* Update docs

Co-authored-by: Dylan Beaudette <[email protected]>

R/SDA_hydric.R

@@ -24,47 +24,47 @@ where_clause <- switch(as.character(is.null(areasymbols)),
                        "TRUE" = sprintf("mu.mukey IN %s", mukeys),
                        "FALSE" = sprintf("l.areasymbol IN %s", areasymbols))
 
-q <- sprintf("SELECT AREASYMBOL,
-        MUSYM,
-        MUNAME,
-        mu.mukey/1  AS MUKEY,
+q <- sprintf("SELECT areasymbol,
+        musym,
+        muname,
+        mu.mukey/1 AS mukey,
         (SELECT TOP 1 COUNT_BIG(*)
         FROM mapunit
-        INNER JOIN component ON component.mukey=mapunit.mukey AND mapunit.mukey = mu.mukey) AS comp_count,
+        INNER JOIN component ON component.mukey = mapunit.mukey AND mapunit.mukey = mu.mukey) AS comp_count,
         (SELECT TOP 1 COUNT_BIG(*)
         FROM mapunit
-        INNER JOIN component ON component.mukey=mapunit.mukey AND mapunit.mukey = mu.mukey
+        INNER JOIN component ON component.mukey = mapunit.mukey AND mapunit.mukey = mu.mukey
         AND majcompflag = 'Yes') AS count_maj_comp,
         (SELECT TOP 1 COUNT_BIG(*)
         FROM mapunit
-        INNER JOIN component ON component.mukey=mapunit.mukey AND mapunit.mukey = mu.mukey
+        INNER JOIN component ON component.mukey = mapunit.mukey AND mapunit.mukey = mu.mukey
         AND hydricrating = 'Yes' ) AS all_hydric,
         (SELECT TOP 1 COUNT_BIG(*)
         FROM mapunit
-        INNER JOIN component ON component.mukey=mapunit.mukey AND mapunit.mukey = mu.mukey
+        INNER JOIN component ON component.mukey = mapunit.mukey AND mapunit.mukey = mu.mukey
         AND majcompflag = 'Yes' AND hydricrating = 'Yes') AS maj_hydric,
         (SELECT TOP 1 COUNT_BIG(*)
         FROM mapunit
-        INNER JOIN component ON component.mukey=mapunit.mukey AND mapunit.mukey = mu.mukey
+        INNER JOIN component ON component.mukey = mapunit.mukey AND mapunit.mukey = mu.mukey
         AND majcompflag = 'Yes' AND hydricrating != 'Yes') AS maj_not_hydric,
          (SELECT TOP 1 COUNT_BIG(*)
         FROM mapunit
         INNER JOIN component ON component.mukey=mapunit.mukey AND mapunit.mukey = mu.mukey
         AND majcompflag != 'Yes' AND hydricrating  = 'Yes' ) AS hydric_inclusions,
         (SELECT TOP 1 COUNT_BIG(*)
         FROM mapunit
-        INNER JOIN component ON component.mukey=mapunit.mukey AND mapunit.mukey = mu.mukey
+        INNER JOIN component ON component.mukey = mapunit.mukey AND mapunit.mukey = mu.mukey
         AND hydricrating  != 'Yes') AS all_not_hydric,
          (SELECT TOP 1 COUNT_BIG(*)
         FROM mapunit
-        INNER JOIN component ON component.mukey=mapunit.mukey AND mapunit.mukey = mu.mukey
+        INNER JOIN component ON component.mukey = mapunit.mukey AND mapunit.mukey = mu.mukey
         AND hydricrating  IS NULL ) AS hydric_null
         INTO #main_query
         FROM legend AS l
         INNER JOIN mapunit AS mu ON mu.lkey = l.lkey AND %s
 
 
-        SELECT  AREASYMBOL, MUKEY, MUSYM, MUNAME,
+        SELECT areasymbol, mukey, musym, muname,
         CASE WHEN comp_count = all_not_hydric + hydric_null THEN  'Nonhydric'
         WHEN comp_count = all_hydric  THEN 'Hydric'
         WHEN comp_count != all_hydric AND count_maj_comp = maj_hydric THEN 'Predominantly Hydric'

misc/getSDA_vectorize.R

@@ -0,0 +1,26 @@
+library(soilDB)
+
+# WORKS (dominant condition)
+get_SDA_property(c("taxsuborder", "taxorder", "tfact"),
+                 method = "dominant condition",
+                 areasymbols = "CA630")
+# WORKS (min/max)
+get_SDA_property(c("ksat_l", "ksat_r", "ksat_h"),
+                 method = "min",
+                 areasymbols = "CA630")
+
+# WORKS (weighted average)
+get_SDA_property(c("ksat_l", "ksat_r", "ksat_h"),
+                 method = "weighted average",
+                 areasymbols = 'CA630')
+
+# WORKS (dominant component, numeric -- special case of weighted average)
+q <- get_SDA_property(
+  c("ksat_l", "ksat_r", "ksat_h"),
+  method = "dominant component (numeric)",
+  areasymbols = 'CA630',
+  query_string = TRUE # this just returns the query instead of calling SDA_query
+)
+# cat(q)
+res <- SDA_query(q)
+res

misc/soil-data-aggregation/get_SDA_SOD-tests.R

@@ -0,0 +1,40 @@
+library(soilDB)
+
+# select a single mukey to check that top and bottom depth and different methods produce stable results
+keys <- c(463276)
+
+s1 <- get_SDA_property(property = 'Available Water Capacity - Rep Value', 
+                       method = 'Weighted Average', 
+                       mukeys = keys, 
+                       top_depth = 0, 
+                       bottom_depth = 200)
+s2 <- get_SDA_property(property = 'Available Water Capacity - Rep Value', 
+                       method = 'Weighted Average', 
+                       mukeys = keys, 
+                       top_depth = 0, 
+                       bottom_depth = 1000)
+testthat::expect_equal(s1, s2)
+testthat::expect_equal(s1$awc_r, 0.13)
+
+s3 <- get_SDA_property(property = 'Available Water Capacity - Rep Value', 
+                       method = 'Dominant Component (Numeric)', 
+                       mukeys = keys, 
+                       top_depth = 0, 
+                       bottom_depth = 200)
+s4 <- get_SDA_property(property = 'Available Water Capacity - Rep Value', 
+                       method = 'Dominant Component (Numeric)', 
+                       mukeys = keys, 
+                       top_depth = 0, 
+                       bottom_depth = 1000)
+testthat::expect_equal(s3, s4)
+testthat::expect_equal(s3$awc_r, 0.18)
+
+s5 <- get_SDA_interpretation(rulename = 'FOR - Mechanical Planting Suitability', 
+                       method = 'Weighted Average', 
+                       mukeys = keys)
+testthat::expect_equal(s5$rating_FORMechanicalPlantingSuitability, 0.99)
+
+s6 <- get_SDA_interpretation(rulename = 'FOR - Mechanical Planting Suitability', 
+                       method = 'Dominant Component', 
+                       mukeys = keys)
+testthat::expect_equal(s6$rating_FORMechanicalPlantingSuitability, 0.987)

tests/testthat/test-SDA_interpretations.R

@@ -13,9 +13,10 @@ test_that("SDA interpretations (dominant component) works", {
                                 method = "Dominant Component", areasymbols = target_areas)
   expect_equal(nrow(res), target_area_rows)
 
-  res <- get_SDA_interpretation("FOR - Potential Seedling Mortality",
+  res <- get_SDA_interpretation(c("FOR - Potential Seedling Mortality",
+                                  "FOR - Road Suitability (Natural Surface)"),
                                 method = "Dominant Component", mukeys = target_mukeys)
-  expect_equal(sort(res$MUKEY), sort(target_mukeys))
+  expect_equal(sort(res$mukey), sort(target_mukeys))
 })
 
 test_that("SDA interpretations (dominant condition) works", {
@@ -28,9 +29,10 @@ test_that("SDA interpretations (dominant condition) works", {
   expect_equal(nrow(res), target_area_rows)
 
 
-  res <- get_SDA_interpretation("FOR - Potential Seedling Mortality",
+  res <- get_SDA_interpretation(c("FOR - Potential Seedling Mortality",
+                                  "FOR - Road Suitability (Natural Surface)"),
                                 method = "Dominant Condition", mukeys = target_mukeys)
-  expect_equal(sort(res$MUKEY), sort(target_mukeys))
+  expect_equal(sort(res$mukey), sort(target_mukeys))
 })
 
 test_that("SDA interpretations (weighted average) works", {
@@ -43,20 +45,22 @@ test_that("SDA interpretations (weighted average) works", {
   expect_equal(nrow(res), target_area_rows)
 
 
-  res <- get_SDA_interpretation("FOR - Potential Seedling Mortality",
+  res <- get_SDA_interpretation(c("FOR - Potential Seedling Mortality",
+                                  "FOR - Road Suitability (Natural Surface)"),
                                 method = "Weighted Average", mukeys = target_mukeys)
-  expect_equal(sort(res$MUKEY), sort(target_mukeys))
+  expect_equal(sort(res$mukey), sort(target_mukeys))
 })
 
 test_that("SDA interpretations (no aggregation) works", {
   skip_if_offline()
-  
+
   skip_on_cran()
-
-  res <- get_SDA_interpretation("FOR - Potential Seedling Mortality",
-                                method = "NONE", 
+
+  res <- get_SDA_interpretation(c("FOR - Potential Seedling Mortality",
+                                  "FOR - Road Suitability (Natural Surface)"),
+                                method = "NONE",
                                 areasymbols = target_areas)
   expect_equal(nrow(res), target_area_rows_all)
-  
-  
+
+
 })

tests/testthat/test-SDA_properties.R

@@ -16,7 +16,7 @@ test_that("SDA properties (dominant condition) works", {
                                      method = "Dominant Condition",
                                      areasymbols = target_areas)), target_area_rows)
 
-  expect_equal(get_SDA_property(property = "Taxonomic Suborder",
+  expect_equal(get_SDA_property(property = c("Taxonomic Suborder","Taxonomic Order"),
                                      method = "Dominant Condition",
                                      mukeys = target_mukeys)$mukey, target_mukeys)
 
@@ -30,7 +30,7 @@ test_that("SDA properties (dominant component category) works", {
                                      method = "Dominant Component (Category)",
                                      areasymbols = target_areas)), target_area_rows)
 
-  expect_equal(get_SDA_property(property = "Taxonomic Suborder",
+  expect_equal(get_SDA_property(property = c("Taxonomic Suborder","Taxonomic Order"),
                                      method = "Dominant Component (Category)",
                                      mukeys = target_mukeys)$mukey, target_mukeys)
 })
@@ -49,7 +49,7 @@ test_that("SDA properties (dominant component numeric) works", {
   )), target_area_rows)
 
   expect_equal(get_SDA_property(
-    property = "Very Coarse Sand - Rep Value",
+    property = c("sandvc_l","sandvc_r","sandvc_h"),
     method = "Dominant Component (Numeric)",
     mukeys = target_mukeys,
     top_depth = 25,
@@ -72,7 +72,7 @@ test_that("SDA properties (weighted average) works", {
   )), target_area_rows)
 
   expect_equal(get_SDA_property(
-    property = "Total Clay - Rep Value",
+    property = c("claytotal_l","claytotal_r","claytotal_h"),
     method = "Weighted Average",
     mukeys = target_mukeys,
     top_depth = 25,
@@ -93,7 +93,7 @@ test_that("SDA properties (min/max) works", {
   )), target_area_rows)
 
   expect_equal(get_SDA_property(
-    property = "Saturated Hydraulic Conductivity - Rep Value",
+    property = c("ksat_l","ksat_r","ksat_h"),
     method = "Min/Max",
     mukeys = target_mukeys,
     FUN = "MIN"
@@ -102,17 +102,17 @@ test_that("SDA properties (min/max) works", {
 
 test_that("SDA properties (no aggregation) works", {
   skip_if_offline()
-  
+
   skip_on_cran()
-  
+
   # return results 1:1 with component for component properties
-  expect_equal(nrow(get_SDA_property(property = "Taxonomic Suborder",
+  expect_equal(nrow(get_SDA_property(property = c('rsprod_l','rsprod_r','rsprod_h'),
                                      method = "NONE",
                                      areasymbols = target_areas)), target_area_rows_all)
-  
-  
+
+
   # return results 1:1 with chorizon for horizon properties (includes cokey)
-  expect_equal(nrow(get_SDA_property("Total Sand - Rep Value", 
-                    method = "NONE", 
+  expect_equal(nrow(get_SDA_property(c('sandtotal_l','sandtotal_r','sandtotal_h'),
+                    method = "NONE",
                     areasymbols = target_areas)), target_area_rows_all_chorizon)
 })