gSSURGO: Fix spatial sampling and improve data aggregation accuracy

CHANGELOG.md

@@ -64,6 +64,11 @@ For more information about this file see also [Keep a Changelog](http://keepacha
 - Fixed a bugs and BADM now process both single-site and multi-site settings, detecting the input structure and processing each site independently to generate the correct number of ensemble members per site.
 - Fixed "external pointer is not valid" error and addressed key bugs in `soilgrids_soilC_extract()` function (#3506)
 - Fixed a bug within the `model2netcdf.SIPNET` function where we assumed the constant calculations of `pecan_start_doy` across years (the calculations should vary depending on the last date from the last loop and the start date of the current loop), which will lead to incorrect calculations of the start `sub_dates` and `sub_dates_cf` if we are jumping between years (e.g., from 2012-12-31 to 2013-01-01). The `sipnet2datetime` function is no longer used anywhere and therefore has been removed.
+- `extract_soil_gssurgo()` (#3643)
+  - Replaced point-based WFS queries with raster-based WCS approach using `soilDB::mukey.wcs()`.
+  - Replaced `grid_size`/`grid_spacing` parameters with `radius` (meters) for simpler buffer-based AOI creation.
+  - Switched to single `soilDB::fetchSDA()` call for component-level soil data retrieval, enabling better ensemble uncertainty quantification through within-map-unit variability. Added support for custom AOI polygons. 
+  - This eliminates spatial coverage gaps and reduces network requests while maintaining backward compatibility.
 
 ### Changed
 

modules/data.land/DESCRIPTION

@@ -25,6 +25,7 @@ URL: https://pecanproject.github.io
 BugReports: https://github.com/PecanProject/pecan/issues
 Depends: R (>= 4.1.0)
 Imports:
+    aqp,
     coda,
     curl,
     doSNOW,
@@ -53,6 +54,7 @@ Imports:
     rlang,
     sf,
     sirt,
+    soilDB,
     sp,
     stringr,
     terra,

modules/data.land/NAMESPACE

@@ -35,6 +35,7 @@ export(get.attributes)
 export(get.soil)
 export(get_resource_map)
 export(get_veg_module)
+export(gssurgo_fetch_area)
 export(ic_process)
 export(id_resolveable)
 export(load_veg)

modules/data.land/NEWS.md

@@ -1,5 +1,12 @@
 # PEcAn.data.land 1.9.0.9000
 
+## Changed
+
+* `extract_soil_gssurgo()` now accepts a polygon area of interest or a circular buffer radius,
+     replacing the rectangular area previously specified through the `grid_size`/`grid_spacing` parameters.
+     It also now performs gSSURGO queries through the `soilDB` package, which should be faster and more
+     reliable than the previous hand-rolled XML queries (#3643).
+
 ## Fixed
 
 * `soil2netcdf()` no longer drops depth information for soils with only one layer. (#3785)

modules/data.land/R/gSSURGO_Query.R

@@ -134,3 +134,134 @@ gSSURGO.Query <- function(mukeys,
 }
 
 
+#' Fetch gSSURGO soil data for an area of interest
+#'
+#' Retrieves soil property data from the USDA gSSURGO database using the soilDB package.
+#' This function performs the data retrieval step separately from ensemble generation,
+#' enabling inspection of raw database values and simpler unit testing.
+#'
+#' @param lat Latitude of center point (optional if aoi provided)
+#' @param lon Longitude of center point (optional if aoi provided)
+#' @param aoi Custom area of interest as sf or terra polygon (optional)
+#' @param radius Buffer radius in meters around lat/lon point (default: 500)
+#' @param depths Soil depth breakpoints in meters, must start with 0 
+#'
+#' @return A list containing:
+#'   \item{soilprop}{Data frame with component-level soil properties}
+#'   \item{mukey_counts}{Table of mapunit key pixel counts (for area weighting)}
+#'   \item{depths_cm}{Depth breakpoints in centimeters}
+#'
+#' @export
+#' @author Akash
+#' @examples
+#' \dontrun{
+#'   result <- gssurgo_fetch_area(lat = 40.1, lon = -88.2)
+#'   head(result$soilprop)
+#' }
+gssurgo_fetch_area <- function(lat = NULL, lon = NULL, aoi = NULL, 
+                                radius = 500, depths = c(0, 0.15, 0.30, 0.60)) {
+
+  # Validate inputs
+  if (is.null(aoi) && (is.null(lat) || is.null(lon))) {
+    PEcAn.logger::logger.severe("Must provide either 'aoi' OR both 'lat' and 'lon'")
+  }
+
+  # Create AOI from point + radius if needed
+  if (is.null(aoi)) {
+    aoi <- data.frame(lon = lon, lat = lat) %>%
+      terra::vect(crs = "epsg:4326") %>%
+      terra::buffer(width = radius)
+  }
+
+  # Validate depths parameter
+  if (depths[1] != 0) {
+    PEcAn.logger::logger.severe(
+      "First depth must be 0. Use depths = c(0, 0.15, 0.30, ...) like hist() breaks."
+    )
+  }
+
+  PEcAn.logger::logger.info("Querying gSSURGO Web Coverage Service for map unit keys")
+  mu_raster <- soilDB::mukey.wcs(aoi = aoi, db = 'gSSURGO', res = 30)
+
+  mukey_values <- terra::values(mu_raster)
+  mukey_values <- mukey_values[!is.na(mukey_values)]
+  mukey_counts <- table(mukey_values)
+  mukeys_all <- as.character(names(mukey_counts))
+
+  if (length(mukeys_all) == 0) {
+    PEcAn.logger::logger.severe("No mapunit keys were found for this site.")
+  }
+
+  # Fetch complete soil data via soilDB
+  sda_data <- tryCatch({
+    soilDB::fetchSDA(
+      WHERE = paste0("mukey IN (", paste(mukeys_all, collapse = ","), ")"),
+      duplicates = TRUE,
+      childs = TRUE,
+      nullFragsAreZero = TRUE,
+      rmHzErrors = TRUE
+    )
+  }, error = function(e) {
+    PEcAn.logger::logger.error(paste("Failed to fetch SDA data:", e$message))
+    return(NULL)
+  })
+
+  if (is.null(sda_data)) {
+    PEcAn.logger::logger.error("Could not retrieve soil data from SDA")
+    return(NULL)
+  }
+
+  hz_data <- aqp::horizons(sda_data)
+  site_data <- aqp::site(sda_data)
+
+  # Component-level aggregation by depth
+  depths_cm <- depths * 100
+  all_soil_data <- list()
+
+  for (i in seq_len(length(depths_cm) - 1)) {
+    top_depth <- depths_cm[i]
+    bottom_depth <- depths_cm[i + 1]
+
+    depth_hz <- hz_data %>%
+      dplyr::filter(hzdept_r < bottom_depth & hzdepb_r > top_depth)
+
+    if (nrow(depth_hz) == 0) next
+
+    component_data <- depth_hz %>%
+      dplyr::left_join(site_data[, c("cokey", "comppct_r", "mukey")], by = "cokey") %>%
+      dplyr::mutate(
+        hz_top_adj = pmax(hzdept_r, top_depth),
+        hz_bot_adj = pmin(hzdepb_r, bottom_depth),
+        hz_thickness = hz_bot_adj - hz_top_adj
+      ) %>%
+      dplyr::group_by(mukey, cokey, comppct_r) %>%
+      dplyr::summarise(
+        sandtotal_r = stats::weighted.mean(sandtotal_r, hz_thickness, na.rm = TRUE),
+        silttotal_r = stats::weighted.mean(silttotal_r, hz_thickness, na.rm = TRUE),
+        claytotal_r = stats::weighted.mean(claytotal_r, hz_thickness, na.rm = TRUE),
+        om_r = stats::weighted.mean(om_r, hz_thickness, na.rm = TRUE),
+        dbthirdbar_r = stats::weighted.mean(dbthirdbar_r, hz_thickness, na.rm = TRUE),
+        fragvol_r = stats::weighted.mean(fragvol_r, hz_thickness, na.rm = TRUE),
+        .groups = "drop"
+      ) %>%
+      dplyr::mutate(
+        tex_sum = sandtotal_r + silttotal_r + claytotal_r,
+        sandtotal_r = sandtotal_r / tex_sum * 100,
+        silttotal_r = silttotal_r / tex_sum * 100,
+        claytotal_r = claytotal_r / tex_sum * 100,
+        hzdept_r = top_depth,
+        hzdepb_r = bottom_depth
+      ) %>%
+      dplyr::select(-tex_sum)
+
+    all_soil_data[[i]] <- component_data
+  }
+
+  soilprop <- do.call(rbind, all_soil_data)
+
+  return(list(
+    soilprop = soilprop,
+    mukey_counts = mukey_counts,
+    depths_cm = depths_cm
+  ))
+}

modules/data.land/R/soil_process.R

@@ -58,19 +58,13 @@ soil_process <- function(settings, input, dbfiles, overwrite = FALSE,run.local=T
     names(newfile) <- rep("path", length(newfile))
 
     if(length(newfile)==0){
-      radius <- ifelse(is.null(settings$run$input$soil$radius), 100, 
+      radius <- ifelse(is.null(settings$run$input$soil$radius), 500, 
                        as.numeric(settings$run$input$soil$radius))
-      grid_size <- max(3, ifelse(is.null(settings$run$input$soil$grid_size), 3, 
-                                 as.numeric(settings$run$input$soil$grid_size)))
-
-      grid_extent <- radius * sqrt(pi)
-      grid_spacing <- grid_extent / (grid_size - 1)
       newfile <- extract_soil_gssurgo(
         outfolder, 
         lat = latlon$lat, 
         lon = latlon$lon,
-        grid_size = grid_size,
-        grid_spacing = grid_spacing
+        radius = radius
       )
 
       # register files in DB