PecanProject · mdietze · Sep 27, 2022 · Jul 1, 2022 · Jul 3, 2022 · Jul 17, 2022
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -51,6 +51,7 @@ convert data for a single PFT fixed (#1329, #2974, #2981)
   (#2989; @nanu1605)
 - Fixed a bug with ED2 where ED2IN tags supplied in `settings` that were not in the ED2IN template file were not getting added to ED2IN config files (#3034)
 - Fixed a bug where warnings were printed for file paths on remote servers even when they did exist (#3020)
+- `PEcAn.data.land::gSSURGO.Query` has been updated to work again after changes to the gSSURGO API.
 
 ### Changed
 

diff --git a/Makefile.depends b/Makefile.depends
@@ -14,7 +14,7 @@ $(call depends,modules/assim.sequential): | .install/base/db .install/base/logge
 $(call depends,modules/benchmark): | .install/base/db .install/base/logger .install/base/remote .install/base/settings .install/base/utils .install/modules/data.land
 $(call depends,modules/data.atmosphere): | .install/base/db .install/base/logger .install/base/remote .install/base/utils
 $(call depends,modules/data.hydrology): | .install/base/logger .install/base/utils
-$(call depends,modules/data.land): | .install/modules/benchmark .install/modules/data.atmosphere .install/base/db .install/base/logger .install/base/remote .install/base/settings .install/base/utils .install/base/visualization
+$(call depends,modules/data.land): | .install/modules/benchmark .install/modules/data.atmosphere .install/base/db .install/base/logger .install/base/remote .install/base/utils .install/base/visualization .install/base/settings
 $(call depends,modules/data.remote): | .install/base/db .install/base/utils .install/base/logger .install/base/remote
 $(call depends,modules/emulator): | .install/base/logger
 $(call depends,modules/meta.analysis): | .install/base/utils .install/base/db .install/base/logger .install/base/settings

diff --git a/docker/depends/pecan.depends.R b/docker/depends/pecan.depends.R
@@ -101,7 +101,6 @@ wanted <- c(
 'randtoolbox',
 'raster',
 'rcrossref',
-'RCurl',
 'REddyProc',
 'redland',
 'reshape',

diff --git a/modules/data.land/DESCRIPTION b/modules/data.land/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: PEcAn.data.land
 Type: Package
 Title: PEcAn Functions Used for Ecological Forecasts and Reanalysis
-Version: 1.7.2
+Version: 1.7.2.9000
 Date: 2021-10-04
 Authors@R: c(person("Mike", "Dietze", role = c("aut", "cre"),
                      email = "dietze@bu.edu"),
@@ -29,9 +29,9 @@ Imports:
     dplyr,
     dplR,
     fs,
+    httr,
     lubridate,
     magrittr,
-    maptools,
     mvtnorm,
     ncdf4 (>= 1.15),
     neonUtilities,
@@ -42,14 +42,11 @@ Imports:
     PEcAn.DB,
     PEcAn.logger,
     PEcAn.remote,
-    PEcAn.settings,
     PEcAn.utils,
     PEcAn.visualization,
     purrr,
-    RCurl,
     rjags,
     rlang,
-    RPostgreSQL,
     sf,
     sirt,
     sp,
@@ -58,9 +55,12 @@ Imports:
     XML (>= 3.98-1.4)      
 Suggests:
     dataone,
+    maptools,
+    PEcAn.settings,
     redland,
     raster,
     rgdal,
+    RPostgreSQL,
     testthat (>= 1.0.2)
 License: BSD_3_clause + file LICENSE
 Copyright: Authors

diff --git a/modules/data.land/NEWS.md b/modules/data.land/NEWS.md
@@ -0,0 +1,11 @@
+# PEcAn.data.land 1.7.2.9000
+
+## Fixed
+
+* `gSSURGO.Query()` now always returns all the columns requested, even ones that are all NA. It also now always requires `mukeys` to be specified.
+* Updated `gSSURGO.Query()` and `extract_soil_gssurgo()` to work again after formatting changes in the underlying gSSURGO API
+
+# PEcAn.data.land 1.7.1
+
+* All changes in 1.7.1 and earlier were recorded in a single file for all of the PEcAn packages; please see 
+https://github.com/PecanProject/pecan/blob/v1.7.1/CHANGELOG.md for details.
diff --git a/modules/data.land/R/extract_soil_nc.R b/modules/data.land/R/extract_soil_nc.R
@@ -8,7 +8,8 @@
 #' @param depths  Standard set of soil depths in m to create the ensemble of soil profiles with.
 #'
 #' @return It returns the address for the generated soil netcdf file
-#' @export
+#'
+#' @importFrom rlang .data
 #'
 #' @examples
 #' \dontrun{
@@ -18,6 +19,7 @@
 #'    PEcAn.data.land::extract_soil_gssurgo(outdir, lat, lon)
 #' }
 #' @author Hamze Dokoohaki
+#' @export
 #' 
 extract_soil_gssurgo<-function(outdir, lat, lon, size=1, radius=500, depths=c(0.15,0.30,0.60)){
   # I keep all the ensembles here 
@@ -27,93 +29,59 @@ extract_soil_gssurgo<-function(outdir, lat, lon, size=1, radius=500, depths=c(0.
   # Basically I think of this as me going around and taking soil samples within 500m of my site.
   #https://sdmdataaccess.nrcs.usda.gov/SpatialFilterHelp.htm
   mu.Path <- paste0(
-    "https://sdmdataaccess.nrcs.usda.gov/Spatial/SDMWGS84Geographic.wfs?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetFeature&TYPENAME=MapunitPoly&FILTER=<Filter><DWithin><PropertyName>Geometry</PropertyName><gml:Point>%20<gml:coordinates>",
-    lon ,
-    ",",
-    lat,
-    "</gml:coordinates></gml:Point><Distance%20units=%27m%27>",radius,"</Distance></DWithin></Filter>"
+    "https://sdmdataaccess.nrcs.usda.gov/Spatial/SDMWGS84Geographic.wfs?",
+    "SERVICE=WFS",
+    "&VERSION=1.1.0",
+    "&REQUEST=GetFeature&TYPENAME=MapunitPoly",
+    "&FILTER=",
+      "<Filter>",
+        "<DWithin>",
+          "<PropertyName>Geometry</PropertyName>",
+          "<gml:Point>",
+            "<gml:coordinates>", lon, ",", lat, "</gml:coordinates>",
+          "</gml:Point>",
+          "<Distance%20units=%27m%27>", radius, "</Distance>",
+        "</DWithin>",
+      "</Filter>",
+    "&OUTPUTFORMAT=XMLMukeyList"
   )
-  # We are trying to find the mapunit key here. Either using rgdal if the driver is defined or parsing it's gml
 
-  if ("GML" %in% rgdal::ogrDrivers()$name) {
-
-    suppressMessages({
+  xmll <- curl::curl_download(
+    mu.Path,
+    ssl.verifyhost = FALSE,
+    ssl.verifypeer = FALSE)
 
-    #disambiguateFIDs if TRUE, and FID values are not unique, they will be set to unique values 1:N for N features; problem observed in GML files
-    #idk why but gssurgo api seems to fail for no reason, that'w why I try 3 times.
-      for(i in 1:3){
-       # try to read the polygon to get the mukey
-       sarea <-try(rgdal::readOGR(mu.Path, disambiguateFIDs=T), silent = TRUE)
-       if( class(sarea) != "try-error" ) break;
-       PEcAn.logger::logger.warn(paste0(i, " attemp was unsuccessful"))
-       Sys.sleep(1)
-     }
+  mukey_str <- XML::xpathApply(
+    doc = XML::xmlParse(xmll),
+    path = "//MapUnitKeyList",
+    fun = XML::xmlValue)
+  mukeys <- strsplit(mukey_str, ",")[[1]]
 
-      # flipping the coordinates 
-      # gdal reads the gSSUGO layers with fliped coordinateds
-      for (i in seq_along(sarea@polygons)){
-        for (j in seq_along(sarea@polygons[[i]]@Polygons)){
-          # flip the coordinates
-          sarea@polygons[[i]]@Polygons[[j]]@coords <- sarea@polygons[[i]]@Polygons[[j]]@coords[, c(2,1)]
-        }
-      }
-
-      areasf <-sf::st_as_sf(sarea)
-      # getting the site point ready
-      site = sf::st_as_sf(data.frame(long=lon, lat=lat), coords=c("long","lat"))
-
-      #buffer the radius around site / and clip the study area based on buffer
-      site_buffer = sf::st_buffer(site, (radius/111000)) # converting radius m to dgree - each degree is about 111 Km
-      site_area = sf::st_intersection(site_buffer, areasf)
-      # calculating areas again for the clipped regions
-     mukey_area <- data.frame(Area=raster::area(x= as(site_area, 'Spatial')),
-                              mukey=site_area$mukey)
-
-    })
-
-    mukeys <- mukey_area$mukey %>% as.character()
-  }else{
-    #reading the mapunit based on latitude and longitude of the site
-    #the output is a gml file which need to be downloaded and read as a spatial file but I don't do that.
-    #I just read the file as a text and parse it out and try to find the mukey==mapunitkey
-    xmll <-
-      curl::curl_download(mu.Path,
-                    ssl.verifyhost = FALSE,
-                    ssl.verifypeer = FALSE
-      )
-
-    startp <- regexpr('<ms:mukey>', xmll)
-    stopp <- regexpr('</ms:mukey>', xmll)
-
-    #if you found the mapunit key
-    if (startp == -1 |
-        stopp == -1)
-      PEcAn.logger::logger.error("There was no mapunit keys found for this site.")
-    mukeys <-substr(xmll, startp %>% as.numeric + 10, stopp %>% as.numeric - 1)  
+  if (length(mukeys) == 0) {
+    PEcAn.logger::logger.error("No mapunit keys were found for this site.")
   }
-  
+
   # calling the query function sending the mapunit keys
-  soilprop <- gSSURGO.Query(mukeys,c("chorizon.sandtotal_r",
-                                    "chorizon.silttotal_r",
-                                    "chorizon.claytotal_r",
-                                    "chorizon.hzdept_r"))
+  soilprop <- gSSURGO.Query(
+    mukeys,
+    c("chorizon.sandtotal_r",
+      "chorizon.silttotal_r",
+      "chorizon.claytotal_r",
+      "chorizon.hzdept_r"))
 
   soilprop.new <- soilprop %>%
-    dplyr::arrange(hzdept_r) %>%
-    dplyr::select(-comppct_r) %>%
-    `colnames<-`(
-      c(
-        "fraction_of_sand_in_soil",
-        "fraction_of_silt_in_soil",
-        "fraction_of_clay_in_soil",
-        "soil_depth",
-        "mukey"
-      )
-    )
-  #unit conversion
-  soilprop.new [, c("fraction_of_sand_in_soil", "fraction_of_silt_in_soil" , "fraction_of_clay_in_soil" ,
-                    "soil_depth")] <- soilprop.new [, c("fraction_of_sand_in_soil", "fraction_of_silt_in_soil" ,
-                                                        "fraction_of_clay_in_soil" , "soil_depth")]/100
+    dplyr::arrange(.data$hzdept_r) %>%
+    dplyr::select(
+      fraction_of_sand_in_soil = .data$sandtotal_r,
+      fraction_of_silt_in_soil = .data$silttotal_r,
+      fraction_of_clay_in_soil = .data$claytotal_r,
+      soil_depth = .data$hzdept_r,
+      mukey = .data$mukey) %>%
+    dplyr::mutate(dplyr::across(
+        c(dplyr::starts_with("fraction_of"),
+          "soil_depth"),
+        function(x) x / 100))
+
   soilprop.new <- soilprop.new[ complete.cases(soilprop.new) , ]
   #converting it to list
   soil.data.gssurgo <- names(soilprop.new)[1:4] %>%

diff --git a/modules/data.land/R/find.land.R b/modules/data.land/R/find.land.R
@@ -11,15 +11,22 @@
 ##' @export
 ##' @author David LeBauer
 find.land <- function(lat, lon, plot = FALSE) {
-  data("wrld_simpl",package="maptools",envir = environment())
+  if (!requireNamespace("maptools", quietly = TRUE)) {
+    PEcAn.logger::logger.error(
+      "PEcAn.data.land::find.land needs package `maptools`.",
+      "Please install it and try again.")
+  }
+  # If maptools used lazy data, this could just be `maptools::wrld_simpl`...
+  wrld_simpl <- NULL
+  utils::data("wrld_simpl",package="maptools",envir = environment())
 
   ## Create a SpatialPoints object
   points <- expand.grid(lon, lat)
   colnames(points) <- c("lat", "lon")
   pts <- sp::SpatialPoints(points, proj4string = sp::CRS(sp::proj4string(wrld_simpl)))
 
   ## Find which points fall over land
-  landmask <- cbind(points, data.frame(land = !is.na(over(pts, wrld_simpl)$FIPS)))
+  landmask <- cbind(points, data.frame(land = !is.na(sp::over(pts, wrld_simpl)$FIPS)))
   if (plot) {
     plot(wrld_simpl)
     landmask[, points(lon, lat, col = 1 + land, pch = 16)]

diff --git a/modules/data.land/R/gSSURGO_Query.R b/modules/data.land/R/gSSURGO_Query.R
@@ -4,46 +4,48 @@
 #' @param mukeys map unit key from gssurgo
 #' @param fields a character vector of the fields to be extracted. See details and the default argument to find out how to define fields.
 #'
-#' @return a dataframe with soil properties. units can be looked up from database documentation
-#' @export
+#' @return a dataframe with soil properties. Units can be looked up from database documentation
 #'
 #' @details 
 #' Full documention of available tables and their relationships can be found here \url{www.sdmdataaccess.nrcs.usda.gov/QueryHelp.aspx}
 #' There have been occasions where NRCS made some minor changes to the structure of the API which this code is where those changes need
 #' to be implemneted here.
 #' Fields need to be defined with their associate tables. For example, sandtotal is a field in chorizon table which needs to be defined as chorizon.sandotal_(r/l/h), where 
-#' r stands for the representative value, l stand for low and h stands for high. At the momeent fields from mapunit, component, muaggatt, and chorizon tables can be extracted. 
+#' r stands for the representative value, l stands for low and h stands for high. At the moment fields from mapunit, component, muaggatt, and chorizon tables can be extracted.
 #'
 #' @examples
 #' \dontrun{
 #'  PEcAn.data.land::gSSURGO.Query(
+#'    mukeys = 2747727,
 #'    fields = c(
 #'      "chorizon.cec7_r", "chorizon.sandtotal_r",
 #'      "chorizon.silttotal_r","chorizon.claytotal_r",
 #'      "chorizon.om_r","chorizon.hzdept_r","chorizon.frag3to10_r",
 #'      "chorizon.dbovendry_r","chorizon.ph1to1h2o_r",
 #'      "chorizon.cokey","chorizon.chkey"))
 #' }
-gSSURGO.Query <- function(mukeys=2747727,
-                        fields=c("chorizon.sandtotal_r",
-                                 "chorizon.silttotal_r",
-                                 "chorizon.claytotal_r")){
-  #browser()
-  #               ,
-  ######### Reteiv soil
-  headerFields <-
-    c(Accept = "text/xml",
-      Accept = "multipart/*",
-      'Content-Type' = "text/xml; charset=utf-8",
-      SOAPAction = "http://SDMDataAccess.nrcs.usda.gov/Tabular/SDMTabularService.asmx/RunQuery")
+#' @export
+#'
+gSSURGO.Query <- function(mukeys,
+                          fields = c("chorizon.sandtotal_r",
+                                     "chorizon.silttotal_r",
+                                     "chorizon.claytotal_r")) {
+
+  ######### Retrieve soil
+
+  # Avoids duplicating fields that are always included in the query
+  fixed_fields <- c("mapunit.mukey", "component.cokey", "component.comppct_r")
+  qry_fields <- unique(fields[!(fields %in% fixed_fields)])
 
   body <- paste('<?xml version="1.0" encoding="utf-8"?>
                <soap:Envelope xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:xsd="http://www.w3.org/2001/XMLSchema" xmlns:soap="http://schemas.xmlsoap.org/soap/envelope/">
                <soap:Body>
                <RunQuery xmlns="http://SDMDataAccess.nrcs.usda.gov/Tabular/SDMTabularService.asmx">
                <Query>
-               SELECT mapunit.mukey, component.cokey, component.mukey, component.comppct_r, ',paste(fields, collapse = ", "),',
-               muaggatt.aws050wta from mapunit
+               SELECT ',
+                 paste(fixed_fields, collapse = ", "),
+                 paste(qry_fields, collapse = ", "),
+                 ' from mapunit
                join muaggatt on mapunit.mukey=muaggatt.mukey
                join component on mapunit.mukey=component.mukey
                join chorizon on component.cokey=chorizon.cokey
@@ -52,14 +54,23 @@ gSSURGO.Query <- function(mukeys=2747727,
                </RunQuery>
                </soap:Body>
                </soap:Envelope>')
-  reader <- RCurl::basicTextGatherer()
-  out <- RCurl::curlPerform(url = "https://SDMDataAccess.nrcs.usda.gov/Tabular/SDMTabularService.asmx",
-                          httpheader = headerFields,  postfields = body,
-                          writefunction = reader$update
-  )
+
+  out <- httr::POST(
+    url = "https://SDMDataAccess.nrcs.usda.gov/Tabular/SDMTabularService.asmx",
+    config = list(
+      httr::accept("text/xml"),
+      httr::accept("multipart/*"),
+      httr::add_headers(
+        SOAPAction = "http://SDMDataAccess.nrcs.usda.gov/Tabular/SDMTabularService.asmx/RunQuery")),
+    httr::content_type("text/xml; charset=utf-8"), # I expected this to belong inside `config`, but doesn't seem to work there...
+    encode="multipart",
+    body = body)
+  httr::stop_for_status(out)
+  result <- httr::content(out, "text")
+
   suppressWarnings(
     suppressMessages({
-      xml_doc <- XML::xmlTreeParse(reader$value())
+      xml_doc <- XML::xmlTreeParse(result)
       xmltop  <- XML::xmlRoot(xml_doc)
       tablesxml <- (xmltop[[1]]["RunQueryResponse"][[1]]["RunQueryResult"][[1]]["diffgram"][[1]]["NewDataSet"][[1]])
     })
@@ -72,25 +83,16 @@ gSSURGO.Query <- function(mukeys=2747727,
         tables <- XML::getNodeSet(tablesxml,"//Table")
 
         ##### All datatables below newdataset
-        # This method leaves out the variables are all NAs  - so we can't have a fixed naming scheme for this df
-        dfs <- tables %>%
-          purrr::map_dfr(function(child){
-            #converting the xml obj to list
-            allfields <- XML::xmlToList(child)
-            remov <- names(allfields) %in% c(".attrs")
-            #browser()
-            names(allfields)[!remov] %>%
-              purrr::map_dfc(function(nfield){
-                #browser()
-                outv <- allfields[[nfield]] %>% unlist() %>% as.numeric
-                ifelse(length(outv) > 0, outv, NA)
-              })%>%
-              as.data.frame() %>%
-              `colnames<-`(names(allfields)[!remov])
-          })%>%
-          dplyr::select(comppct_r:mukey) %>%
-          dplyr::select(-aws050wta)
-
+        dfs <- purrr::map_dfr(
+            tables,
+            function(tbl){
+              lst <- purrr::map(
+                XML::xmlToList(tbl),
+                function(v)ifelse(is.null(v), NA, v)) #avoid dropping empty columns
+
+              lst[names(lst) != ".attrs"]}
+          )
+          dfs <- dplyr::mutate(dfs, dplyr::across(dplyr::everything(), as.numeric))
       })
     )