From bc81e8d07f8a85d06c7f3a4a1119dd0b11c304e0 Mon Sep 17 00:00:00 2001
From: Qianyuli <liqianyu2528@hotmail.com>
Date: Wed, 14 Sep 2022 23:38:44 -0700
Subject: [PATCH 01/14] Add files via upload

Add file that extracts Soilgrids SOC based on user-defined site locations
---
 .../data.land/R/soilgrids_soc_extraction.R    | 149 ++++++++++++++++++
 1 file changed, 149 insertions(+)
 create mode 100644 modules/data.land/R/soilgrids_soc_extraction.R

diff --git a/modules/data.land/R/soilgrids_soc_extraction.R b/modules/data.land/R/soilgrids_soc_extraction.R
new file mode 100644
index 00000000000..ea18a4446e4
--- /dev/null
+++ b/modules/data.land/R/soilgrids_soc_extraction.R
@@ -0,0 +1,149 @@
+###Extraction of total soil organic carbon based on user-defined site location from SoilGrids250m version 2.0 : https://soilgrids.org
+###Input parameter: user-defined site longitude/latitude 
+###Extract variable ocd (organic carbon density)
+###Example:ocs<-ocs_extract(c(-77.904, -119.1944, -121.557), c(40.6658, 37.0678, 44.4524))
+
+rm(list = ls())
+library(terra)
+library(reshape2)
+library(dplyr)
+
+ocs_extract <- function (lon_input, lat_input, outdir=NULL,verbose=TRUE) {
+  #create a variable to store mean and quantile of organic carbon density (ocd) for each soil depth
+  ocdquant <- matrix(NA, nrow = 6, ncol = length(lon_input) * 4) #row represents soil depth, col represents mean, 5%, 50% and 95%-quantile of ocd for all sites 
+  #create a variable for site ID
+  siteid <- matrix(NA, nrow = 6, ncol = length(lon_input) * 4)
+  lonlat <- cbind(lon_input, lat_input)
+  base_data_url <- "/vsicurl?max_retry=30&retry_delay=60&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/ocd/ocd_"
+  depths <- c("0-5cm", "5-15cm", "15-30cm", "30-60cm", "60-100cm", "100-200cm")
+ 
+  # reproject locations to soilgrids projection
+  #Soilgrids data is using Homolosine projection https://www.isric.org/explore/soilgrids/faq-soilgrids 
+  p <- terra::vect(lonlat, crs = "+proj=longlat +datum=WGS84") # Users need to provide lon/lat
+  newcrs <- "+proj=igh +datum=WGS84 +no_defs +towgs84=0,0,0" 
+  p_reproj <- terra::project(p, newcrs) # Transform the point vector to data with Homolosine projection
+  
+  # setup progress bar
+  if (verbose) {
+    j <- 1
+    pb <- txtProgressBar(min = 0, max = length(depths), char="*", width=70, style = 3)
+  }
+  
+  for (dep in seq_along(depths)) {
+    # setup virtual raster URLs for each layer
+    ocd_mean.url <- paste0(base_data_url,depths[dep],"_mean.vrt")
+    ocd_Q0.05.url <- paste0(base_data_url, depths[dep], "_Q0.05.vrt")
+    ocd_Q0.50.url <- paste0(base_data_url, depths[dep], "_Q0.5.vrt")
+    ocd_Q0.95.url <- paste0(base_data_url, depths[dep], "_Q0.95.vrt")
+    
+    # create virtual rasters && extract SOC values - the original unit is hg/m3
+    ocd_mean <- terra::extract(terra::rast(ocd_mean.url), p_reproj)
+    ocd_Q0.05_map <- terra::extract(terra::rast(ocd_Q0.05.url), p_reproj)
+    ocd_Q0.50_map <- terra::extract(terra::rast(ocd_Q0.50.url), p_reproj)
+    ocd_Q0.95_map <- terra::extract(terra::rast(ocd_Q0.95.url), p_reproj)
+    
+    #change the unit to more common kg/m3
+    ocd_mean_real <- ocd_mean[, 2] / 10
+    ocd_Q0.05_real <- ocd_Q0.05_map[, 2] / 10
+    ocd_Q0.50_real <- ocd_Q0.50_map[, 2] / 10
+    ocd_Q0.95_real <- ocd_Q0.95_map[, 2] / 10
+    
+    ocdquant[dep, ] <-c(ocd_mean_real,ocd_Q0.05_real,ocd_Q0.50_real,ocd_Q0.95_real)
+    siteid [dep, ] <- rep(1:length(lon_input), 4)
+
+    ### Display progress to console
+    if (verbose) {
+      setTxtProgressBar(pb, j)
+      j <- j+1
+      flush.console()}
+  
+  # cleanup interim results
+    rm(ocd_mean.url, ocd_Q0.05.url, ocd_Q0.50.url, ocd_Q0.95.url, 
+       ocd_mean, ocd_Q0.05_map, ocd_Q0.50_map, ocd_Q0.95_map,
+       ocd_mean_real, ocd_Q0.05_real, ocd_Q0.50_real, ocd_Q0.95_real)
+  }
+  
+  rownames(ocdquant) <- c("0-5", "5-15", "15-30", "30-60", "60-100", "100-200")
+  colnames(ocdquant) <- c(rep("Mean", length(lon_input)),
+                          rep("0.05", length(lon_input)),
+                          rep("0.5", length(lon_input)),
+                          rep("0.95", length(lon_input)))
+  rownames(siteid) <- c("0-5", "5-15", "15-30", "30-60", "60-100", "100-200")
+  colnames(siteid) <- c(rep("Mean", length(lon_input)),
+                        rep("0.05", length(lon_input)),
+                        rep("0.5", length(lon_input)),
+                        rep("0.95", length(lon_input)))
+  if (verbose) {
+    close(pb)
+  }
+  
+  # parse extracted data and prepare for output
+  ocd_fit <- melt(ocdquant, id.vars = c("Mean"))
+  id_fit <- melt(siteid, id.vars = c("Mean"))
+  colnames(ocd_fit) <- c("Depth", "Quantile", "Value")
+  ocd_fit$Variable <- rep("ocd", length(nrow(ocd_fit)))
+  ocd_fit$siteid <- id_fit$value
+  dat <- split(ocd_fit, list(ocd_fit$siteid, ocd_fit$Depth))
+  
+  #assume the ocd profile follows gamma distribution best
+  cgamma <- function(theta, val, stat) {
+    pred <- rep(NA, 4)
+    names(pred) = stat
+    if ("Mean" %in% stat) {
+      pred["Mean"] <- theta[1] / theta[2]
+    }
+    qstat <- as.numeric(stat)[!is.na(as.numeric(stat))]
+    pred[as.character(qstat)] <- qgamma(qstat, theta[1], theta[2])
+    return(sum((pred - val) ^ 2))
+  }
+  
+  fitQ <- function(x) {
+    val = x$Value
+    stat = as.character(x$Quantile)
+    theta = c(10, 10)
+    fit <-
+      list(Gamma = optim(theta, cgamma, val = val, stat = stat))
+    SS <- sapply(fit, function(f) {
+      f$value
+    })
+    par <- sapply(fit, function(f) {
+      f$par
+    })
+    return(list(par = par, SS = SS))
+  }
+  
+  score <- suppressWarnings(lapply(dat, fitQ))
+  bestPar <- sapply(score, function(f) { f$par })
+  mean <- bestPar[1,] / bestPar[2,]
+  std <- sqrt(bestPar[1,] / bestPar[2,] ^ 2)
+  mean_site <- matrix(mean, length(lon_input), 6)
+  rownames(mean_site) <- paste0("Site_", 1:length(lon_input))
+  colnames(mean_site) <- c("0-5cm",
+                           "5-15cm",
+                           "15-30cm",
+                           "30-60cm",
+                           "60-100cm",
+                           "100-200cm")
+  std_site <- matrix(std, length(lon_input), 6)
+  rownames(std_site) <- paste0("Site_", 1:length(lon_input))
+  colnames(std_site) <- c("0-5cm",
+                          "5-15cm",
+                          "15-30cm",
+                          "30-60cm",
+                          "60-100cm",
+                          "100-200cm")
+  #calculate organic carbon stock (ocs) as the sum of organic carbon density multiplied by layer thickness, the unit of ocs is kg/m2, based on Eq. (6）in paper https://www.sciencedirect.com/science/article/pii/S2215016122000462
+  ocs_sum <- mean_site[,1]*(5-0)*0.01+mean_site[,2]*(15-5)*0.01+mean_site[,3]*(30-15)*0.01+mean_site[,4]*(60-30)*0.01+mean_site[,5]*(100-60)*0.01+mean_site[,6]*(200-100)*0.01 
+  #calculate standard deviation of ocs as the square root of sum of variance of layer-specific ocs, the unit of ocs is kg/m2, based on Eq. (8) in paper https://www.sciencedirect.com/science/article/pii/S2215016122000462, except the correlation term due to the lack of information 
+  ocs_std <- sqrt((std_site[,1]*(5-0)*0.01)^2+(std_site[,2]*(15-5)*0.01)^2+(std_site[,3]*(30-15)*0.01)^2+(std_site[,4]*(60-30)*0.01)^2+(std_site[,5]*(100-60)*0.01)^2+(std_site[,6]*(200-100)*0.01)^2)
+  
+  if (!is.null(outdir)) {
+    PEcAn.logger::logger.info(paste0("Storing results in: ",file.path(outdir,"soilgrids_soc_data.RData")))
+    if (! file.exists(outdir)) dir.create(outdir,recursive=TRUE)
+    soilgrids_soc_data <- list("Total_soc" = ocs_sum, "Sdev_soc" = ocs_std)
+    save(soilgrids_soc_data, file = file.path(outdir,"soilgrids_soc_data.RData"))
+  }
+  # return the results to the terminal as well
+  return(list("Total OCS" = ocs_sum, "Standard deviation of OCS" = ocs_std))
+}
+

From 0e1f5cc79c09213957184b726507ab8e760c3c49 Mon Sep 17 00:00:00 2001
From: "Shawn P. Serbin" <sserbin@bnl.gov>
Date: Thu, 15 Sep 2022 14:38:54 -0400
Subject: [PATCH 02/14] A quick update to help guide documentation and changes
 for function

---
 CHANGELOG.md                                  |  1 +
 modules/data.land/DESCRIPTION                 |  2 +
 modules/data.land/NAMESPACE                   |  1 +
 .../data.land/R/soilgrids_soc_extraction.R    | 37 +++++++++++--------
 modules/data.land/man/ocs_extract.Rd          | 27 ++++++++++++++
 5 files changed, 53 insertions(+), 15 deletions(-)
 create mode 100644 modules/data.land/man/ocs_extract.Rd

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 066f3cc801b..d15dfba1bf1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -22,6 +22,7 @@ see if you need to change any of these:
 - TRAEFIK_HTTP_REDIRECT is no longer used, this is the default when you use https.
 
 ### Added
+- Created a new soilgrids function to extract the mean soil organic carbon profile with associated undertainty values at each depth for any lat/long location (#3040). Function was created for the CMS SDA workflow
 
 - `PEcAn.all` gains new function `pecan_version`, which reports for each PEcAn
   package the version that was provided with this release and the version that
diff --git a/modules/data.land/DESCRIPTION b/modules/data.land/DESCRIPTION
index 6da21757643..0e1e0a90dce 100644
--- a/modules/data.land/DESCRIPTION
+++ b/modules/data.land/DESCRIPTION
@@ -45,6 +45,7 @@ Imports:
     PEcAn.utils,
     PEcAn.visualization,
     purrr,
+    reshape2,
     rjags,
     rlang,
     RCurl,
@@ -53,6 +54,7 @@ Imports:
     sirt,
     sp,
     stringr,
+    terra,
     traits,
     XML (>= 3.98-1.4)      
 Suggests:
diff --git a/modules/data.land/NAMESPACE b/modules/data.land/NAMESPACE
index a191fa4da67..d0a7d4538bc 100644
--- a/modules/data.land/NAMESPACE
+++ b/modules/data.land/NAMESPACE
@@ -36,6 +36,7 @@ export(match_pft)
 export(match_species_id)
 export(mpot2smoist)
 export(netcdf.writer.BADM)
+export(ocs_extract)
 export(parse.MatrixNames)
 export(partition_roots)
 export(plot2AGB)
diff --git a/modules/data.land/R/soilgrids_soc_extraction.R b/modules/data.land/R/soilgrids_soc_extraction.R
index ea18a4446e4..fc137610690 100644
--- a/modules/data.land/R/soilgrids_soc_extraction.R
+++ b/modules/data.land/R/soilgrids_soc_extraction.R
@@ -1,14 +1,21 @@
-###Extraction of total soil organic carbon based on user-defined site location from SoilGrids250m version 2.0 : https://soilgrids.org
-###Input parameter: user-defined site longitude/latitude 
-###Extract variable ocd (organic carbon density)
-###Example:ocs<-ocs_extract(c(-77.904, -119.1944, -121.557), c(40.6658, 37.0678, 44.4524))
-
-rm(list = ls())
-library(terra)
-library(reshape2)
-library(dplyr)
-
-ocs_extract <- function (lon_input, lat_input, outdir=NULL,verbose=TRUE) {
+##' ocs_extract function
+##' A function to extract total soil organic carbon for a single or group of 
+##' lat/long locationsbased on user-defined site location from SoilGrids250m 
+##' version 2.0 : https://soilgrids.org
+##' @title ocs_extract
+##' @name ocs_extract
+##' 
+##' @param site_info A dataframe of site info containing the BETYdb site ID, 
+##' site name, latitude, and longitude, e.g. 
+##' (site_id, site_name, lat, lon)
+##' @param outdir Optional. Provide the results as an Rdata file 
+##' (soilgrids_soc_data.RData)
+##' @param verbose Provide progress feedback to the terminal? TRUE/FALSE
+##' 
+##' @export
+##' @author Qianyu Li, Shawn P. Serbin
+##' 
+ocs_extract <- function (site_info, outdir=NULL, verbose=TRUE, ...) {
   #create a variable to store mean and quantile of organic carbon density (ocd) for each soil depth
   ocdquant <- matrix(NA, nrow = 6, ncol = length(lon_input) * 4) #row represents soil depth, col represents mean, 5%, 50% and 95%-quantile of ocd for all sites 
   #create a variable for site ID
@@ -26,7 +33,7 @@ ocs_extract <- function (lon_input, lat_input, outdir=NULL,verbose=TRUE) {
   # setup progress bar
   if (verbose) {
     j <- 1
-    pb <- txtProgressBar(min = 0, max = length(depths), char="*", width=70, style = 3)
+    pb <- utils::txtProgressBar(min = 0, max = length(depths), char="*", width=70, style = 3)
   }
   
   for (dep in seq_along(depths)) {
@@ -53,7 +60,7 @@ ocs_extract <- function (lon_input, lat_input, outdir=NULL,verbose=TRUE) {
 
     ### Display progress to console
     if (verbose) {
-      setTxtProgressBar(pb, j)
+      utils::setTxtProgressBar(pb, j)
       j <- j+1
       flush.console()}
   
@@ -78,8 +85,8 @@ ocs_extract <- function (lon_input, lat_input, outdir=NULL,verbose=TRUE) {
   }
   
   # parse extracted data and prepare for output
-  ocd_fit <- melt(ocdquant, id.vars = c("Mean"))
-  id_fit <- melt(siteid, id.vars = c("Mean"))
+  ocd_fit <- reshape2::melt(ocdquant, id.vars = c("Mean"))
+  id_fit <- reshape2::melt(siteid, id.vars = c("Mean"))
   colnames(ocd_fit) <- c("Depth", "Quantile", "Value")
   ocd_fit$Variable <- rep("ocd", length(nrow(ocd_fit)))
   ocd_fit$siteid <- id_fit$value
diff --git a/modules/data.land/man/ocs_extract.Rd b/modules/data.land/man/ocs_extract.Rd
new file mode 100644
index 00000000000..37b3bb74025
--- /dev/null
+++ b/modules/data.land/man/ocs_extract.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/soilgrids_soc_extraction.R
+\name{ocs_extract}
+\alias{ocs_extract}
+\title{ocs_extract}
+\usage{
+ocs_extract(site_info, outdir = NULL, verbose = TRUE, ...)
+}
+\arguments{
+\item{site_info}{A dataframe of site info containing the BETYdb site ID, 
+site name, latitude, and longitude, e.g. 
+(site_id, site_name, lat, lon)}
+
+\item{outdir}{Optional. Provide the results as an Rdata file 
+(soilgrids_soc_data.RData)}
+
+\item{verbose}{Provide progress feedback to the terminal? TRUE/FALSE}
+}
+\description{
+ocs_extract function
+A function to extract total soil organic carbon for a single or group of 
+lat/long locationsbased on user-defined site location from SoilGrids250m 
+version 2.0 : https://soilgrids.org
+}
+\author{
+Qianyu Li, Shawn P. Serbin
+}

From 1eafb9b686519557103cfb43734419ea11b3edaf Mon Sep 17 00:00:00 2001
From: Qianyuxuan <1404243250@qq.com>
Date: Fri, 23 Sep 2022 22:24:11 -0700
Subject: [PATCH 03/14] Add header and modify return object with site ID, lat
 and lon.

---
 .../data.land/R/soilgrids_soc_extraction.R    | 155 ++++++++++++------
 1 file changed, 104 insertions(+), 51 deletions(-)

diff --git a/modules/data.land/R/soilgrids_soc_extraction.R b/modules/data.land/R/soilgrids_soc_extraction.R
index fc137610690..71517ece3cb 100644
--- a/modules/data.land/R/soilgrids_soc_extraction.R
+++ b/modules/data.land/R/soilgrids_soc_extraction.R
@@ -1,28 +1,78 @@
-##' ocs_extract function
+##' soilgrids_soilC_extract function
 ##' A function to extract total soil organic carbon for a single or group of 
 ##' lat/long locationsbased on user-defined site location from SoilGrids250m 
 ##' version 2.0 : https://soilgrids.org
-##' @title ocs_extract
-##' @name ocs_extract
+##' @title soilgrids_soilC_extract
+##' @name soilgrids_soilC_extract
 ##' 
 ##' @param site_info A dataframe of site info containing the BETYdb site ID, 
 ##' site name, latitude, and longitude, e.g. 
 ##' (site_id, site_name, lat, lon)
-##' @param outdir Optional. Provide the results as an Rdata file 
-##' (soilgrids_soc_data.RData)
+##' @param outdir Optional. Provide the results as a CSV file 
+##' (soilgrids_soilC_data.csv)
 ##' @param verbose Provide progress feedback to the terminal? TRUE/FALSE
 ##' 
+##' ##' @examples
+##' \dontrun{
+##' 
+##' # Example 1 - using the modex.bnl.gov BETYdb and site IDs to extract data
+##' db <- 'betydb'
+##' host_db <- 'modex.bnl.gov'
+##' db_port <- '5432'
+##' db_user <- 'bety'
+##' db_password <- 'bety'
+##' 
+##' bety <- list(user='bety', password='bety', host='modex.bnl.gov',
+##' dbname='betydb', driver=RPostgres::Postgres(),write=FALSE)
+##' 
+##' con <- DBI::dbConnect(drv=bety$driver, dbname=bety$dbname, host=bety$host, 
+##' password=bety$password, user=bety$user)
+##' 
+##' suppressWarnings(site_qry <- glue::glue_sql("SELECT *, ST_X(ST_CENTROID(geometry)) AS lon,
+##' ST_Y(ST_CENTROID(geometry)) AS lat FROM sites WHERE id IN ({ids*})",
+##' ids = c("676","622","678","766","764"), .con = con))
+##' 
+##' suppressWarnings(qry_results.1 <- DBI::dbSendQuery(con,site_qry))
+##' suppressWarnings(qry_results.2 <- DBI::dbFetch(qry_results.1))
+##' DBI::dbClearResult(qry_results.1)
+##' dbDisconnect(con)
+##' 
+##' site_info <- qry_results.2
+##' verbose <- TRUE
+##' system.time(result_soc <- soilgrids_soilC_extract(site_info=site_info, verbose=verbose))
+##' result_soc
+##' 
+##' }
+##' @importFrom reshape2 melt
+##' @importFrom terra vect project rast extract
+##' @importFrom stats qgamma optim
+##' @importFrom utils flush.console setTxtProgressBar txtProgressBar
+##' 
+##' @return a dataframe containing the total soil carbon values  
+##' and the corresponding standard deviation values (uncertainties) for each location 
+##' Output column names are c("Site_ID","Site_Name","Latitude","Longitude",
+##' "Total_soilC","Std_soilC")
+##' 
 ##' @export
 ##' @author Qianyu Li, Shawn P. Serbin
 ##' 
-ocs_extract <- function (site_info, outdir=NULL, verbose=TRUE, ...) {
+soilgrids_soilC_extract <- function (site_info, outdir=NULL, verbose=TRUE, ...) {
+
+    if (is.null(site_info)) {
+    stop("Please provide a BETY DB site list containing at least the site id and PostGIS geometry\
+    as lon and lat")
+  }
+  
+  # prepare site info for extraction
+  internal_site_info <- data.frame(site_info$id, site_info$sitename, site_info$lat,site_info$lon)
   #create a variable to store mean and quantile of organic carbon density (ocd) for each soil depth
-  ocdquant <- matrix(NA, nrow = 6, ncol = length(lon_input) * 4) #row represents soil depth, col represents mean, 5%, 50% and 95%-quantile of ocd for all sites 
+  ocdquant <- matrix(NA, nrow = 6, ncol = length(internal_site_info$site_info.lon) * 4) #row represents soil depth, col represents mean, 5%, 50% and 95%-quantile of ocd for all sites 
   #create a variable for site ID
-  siteid <- matrix(NA, nrow = 6, ncol = length(lon_input) * 4)
-  lonlat <- cbind(lon_input, lat_input)
+  siteid <- matrix(rep(as.numeric(site_info$id),times=4*6), nrow=6,byrow=TRUE)
+  lonlat <- cbind(internal_site_info$site_info.lon, internal_site_info$site_info.lat)
   base_data_url <- "/vsicurl?max_retry=30&retry_delay=60&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/ocd/ocd_"
   depths <- c("0-5cm", "5-15cm", "15-30cm", "30-60cm", "60-100cm", "100-200cm")
+  layer_thick <- c(0.05,0.10,0.15,0.30,0.40,1.00) # in unit m
  
   # reproject locations to soilgrids projection
   #Soilgrids data is using Homolosine projection https://www.isric.org/explore/soilgrids/faq-soilgrids 
@@ -50,14 +100,12 @@ ocs_extract <- function (site_info, outdir=NULL, verbose=TRUE, ...) {
     ocd_Q0.95_map <- terra::extract(terra::rast(ocd_Q0.95.url), p_reproj)
     
     #change the unit to more common kg/m3
-    ocd_mean_real <- ocd_mean[, 2] / 10
-    ocd_Q0.05_real <- ocd_Q0.05_map[, 2] / 10
-    ocd_Q0.50_real <- ocd_Q0.50_map[, 2] / 10
-    ocd_Q0.95_real <- ocd_Q0.95_map[, 2] / 10
+    ocd_mean_real <- ocd_mean[, -1] / 10
+    ocd_Q0.05_real <- ocd_Q0.05_map[, -1] / 10
+    ocd_Q0.50_real <- ocd_Q0.50_map[, -1] / 10
+    ocd_Q0.95_real <- ocd_Q0.95_map[, -1] / 10
     
     ocdquant[dep, ] <-c(ocd_mean_real,ocd_Q0.05_real,ocd_Q0.50_real,ocd_Q0.95_real)
-    siteid [dep, ] <- rep(1:length(lon_input), 4)
-
     ### Display progress to console
     if (verbose) {
       utils::setTxtProgressBar(pb, j)
@@ -70,16 +118,16 @@ ocs_extract <- function (site_info, outdir=NULL, verbose=TRUE, ...) {
        ocd_mean_real, ocd_Q0.05_real, ocd_Q0.50_real, ocd_Q0.95_real)
   }
   
-  rownames(ocdquant) <- c("0-5", "5-15", "15-30", "30-60", "60-100", "100-200")
-  colnames(ocdquant) <- c(rep("Mean", length(lon_input)),
-                          rep("0.05", length(lon_input)),
-                          rep("0.5", length(lon_input)),
-                          rep("0.95", length(lon_input)))
-  rownames(siteid) <- c("0-5", "5-15", "15-30", "30-60", "60-100", "100-200")
-  colnames(siteid) <- c(rep("Mean", length(lon_input)),
-                        rep("0.05", length(lon_input)),
-                        rep("0.5", length(lon_input)),
-                        rep("0.95", length(lon_input)))
+  rownames(ocdquant) <- depths
+  colnames(ocdquant) <- c(rep("Mean", length(internal_site_info$site_info.lon)),
+                          rep("0.05", length(internal_site_info$site_info.lon)),
+                          rep("0.5", length(internal_site_info$site_info.lon)),
+                          rep("0.95", length(internal_site_info$site_info.lon)))
+  rownames(siteid) <- depths
+  colnames(siteid) <- c(rep("Mean", length(internal_site_info$site_info.lon)),
+                        rep("0.05", length(internal_site_info$site_info.lon)),
+                        rep("0.5", length(internal_site_info$site_info.lon)),
+                        rep("0.95", length(internal_site_info$site_info.lon)))
   if (verbose) {
     close(pb)
   }
@@ -90,7 +138,10 @@ ocs_extract <- function (site_info, outdir=NULL, verbose=TRUE, ...) {
   colnames(ocd_fit) <- c("Depth", "Quantile", "Value")
   ocd_fit$Variable <- rep("ocd", length(nrow(ocd_fit)))
   ocd_fit$siteid <- id_fit$value
-  dat <- split(ocd_fit, list(ocd_fit$siteid, ocd_fit$Depth))
+  f1<-factor(ocd_fit$siteid,levels=unique(ocd_fit$siteid))
+  f2<-factor(ocd_fit$Depth,levels=unique(ocd_fit$Depth))
+  #split data by groups of sites and soil depth, while keeping the original order of each group
+  dat <- split(ocd_fit, list(f1, f2))  
   
   #assume the ocd profile follows gamma distribution best
   cgamma <- function(theta, val, stat) {
@@ -123,34 +174,36 @@ ocs_extract <- function (site_info, outdir=NULL, verbose=TRUE, ...) {
   bestPar <- sapply(score, function(f) { f$par })
   mean <- bestPar[1,] / bestPar[2,]
   std <- sqrt(bestPar[1,] / bestPar[2,] ^ 2)
-  mean_site <- matrix(mean, length(lon_input), 6)
-  rownames(mean_site) <- paste0("Site_", 1:length(lon_input))
-  colnames(mean_site) <- c("0-5cm",
-                           "5-15cm",
-                           "15-30cm",
-                           "30-60cm",
-                           "60-100cm",
-                           "100-200cm")
-  std_site <- matrix(std, length(lon_input), 6)
-  rownames(std_site) <- paste0("Site_", 1:length(lon_input))
-  colnames(std_site) <- c("0-5cm",
-                          "5-15cm",
-                          "15-30cm",
-                          "30-60cm",
-                          "60-100cm",
-                          "100-200cm")
+  mean_site <- matrix(mean, length(internal_site_info$site_info.lon), 6)
+  rownames(mean_site) <- as.numeric(internal_site_info$site_info.id)
+  colnames(mean_site) <- depths
+  mean_site.2 <- data.frame(site_id=internal_site_info$site_info.id, 
+                            lat=internal_site_info$site_info.lat, 
+                            lon=internal_site_info$site_info.lon, 
+                            mean_site)
+  colnames(mean_site.2)[4:9] <-  depths 
+
+  std_site <- matrix(std, length(internal_site_info$site_info.lon), 6)
+  rownames(std_site) <- as.numeric(internal_site_info$site_info.id)
+  colnames(std_site) <- depths
+  std_site.2 <- data.frame(site_id=internal_site_info$site_info.id,  
+                            lat=internal_site_info$site_info.lat,
+                            lon=internal_site_info$site_info.lon, 
+                            std_site)
+  colnames(std_site.2)[4:9] <-  depths 
   #calculate organic carbon stock (ocs) as the sum of organic carbon density multiplied by layer thickness, the unit of ocs is kg/m2, based on Eq. (6）in paper https://www.sciencedirect.com/science/article/pii/S2215016122000462
-  ocs_sum <- mean_site[,1]*(5-0)*0.01+mean_site[,2]*(15-5)*0.01+mean_site[,3]*(30-15)*0.01+mean_site[,4]*(60-30)*0.01+mean_site[,5]*(100-60)*0.01+mean_site[,6]*(200-100)*0.01 
+  ocs_sum <- mean_site[,1]*layer_thick[1]+mean_site[,2]*layer_thick[2]+mean_site[,3]*layer_thick[3]+mean_site[,4]*layer_thick[4]+mean_site[,5]*layer_thick[5]+mean_site[,6]*layer_thick[6] 
   #calculate standard deviation of ocs as the square root of sum of variance of layer-specific ocs, the unit of ocs is kg/m2, based on Eq. (8) in paper https://www.sciencedirect.com/science/article/pii/S2215016122000462, except the correlation term due to the lack of information 
-  ocs_std <- sqrt((std_site[,1]*(5-0)*0.01)^2+(std_site[,2]*(15-5)*0.01)^2+(std_site[,3]*(30-15)*0.01)^2+(std_site[,4]*(60-30)*0.01)^2+(std_site[,5]*(100-60)*0.01)^2+(std_site[,6]*(200-100)*0.01)^2)
-  
+  ocs_std <- sqrt((std_site[,1]*layer_thick[1])^2+(std_site[,2]*layer_thick[2])^2+(std_site[,3]*layer_thick[3])^2+(std_site[,4]*layer_thick[4])^2+(std_site[,5]*layer_thick[5])^2+(std_site[,6]*layer_thick[6])^2)
+  soilgrids_soilC_data <- data.frame(internal_site_info$site_info.id,internal_site_info$site_info.sitename,internal_site_info$site_info.lat,internal_site_info$site_info.lon,ocs_sum,ocs_std)
+  colnames(soilgrids_soilC_data)<- c("Site_ID","Site_Name","Latitude","Longitude","Total_soilC","Std_soilC")
+  rownames(soilgrids_soilC_data) <- NULL
+
   if (!is.null(outdir)) {
-    PEcAn.logger::logger.info(paste0("Storing results in: ",file.path(outdir,"soilgrids_soc_data.RData")))
-    if (! file.exists(outdir)) dir.create(outdir,recursive=TRUE)
-    soilgrids_soc_data <- list("Total_soc" = ocs_sum, "Sdev_soc" = ocs_std)
-    save(soilgrids_soc_data, file = file.path(outdir,"soilgrids_soc_data.RData"))
+    PEcAn.logger::logger.info(paste0("Storing results in: ",file.path(outdir,"soilgrids_soilC_data.csv")))
+    utils::write.csv(soilgrids_soilC_data,file=file.path(outdir,"soilgrids_soilC_data.csv"),row.names = FALSE)
   }
   # return the results to the terminal as well
-  return(list("Total OCS" = ocs_sum, "Standard deviation of OCS" = ocs_std))
+  return(soilgrids_soilC_data)
 }
-
+ 

From 98a1223b011e5998c1c9e84c3d241bad1471cefd Mon Sep 17 00:00:00 2001
From: Li <qli1@modex.bnl.gov>
Date: Tue, 27 Sep 2022 12:59:13 -0400
Subject: [PATCH 04/14] Update document files

---
 modules/data.land/NAMESPACE                   | 12 +++-
 modules/data.land/man/ocs_extract.Rd          | 27 --------
 .../data.land/man/soilgrids_soilC_extract.Rd  | 65 +++++++++++++++++++
 3 files changed, 76 insertions(+), 28 deletions(-)
 delete mode 100644 modules/data.land/man/ocs_extract.Rd
 create mode 100644 modules/data.land/man/soilgrids_soilC_extract.Rd

diff --git a/modules/data.land/NAMESPACE b/modules/data.land/NAMESPACE
index d0a7d4538bc..e0da8ccd740 100644
--- a/modules/data.land/NAMESPACE
+++ b/modules/data.land/NAMESPACE
@@ -36,7 +36,6 @@ export(match_pft)
 export(match_species_id)
 export(mpot2smoist)
 export(netcdf.writer.BADM)
-export(ocs_extract)
 export(parse.MatrixNames)
 export(partition_roots)
 export(plot2AGB)
@@ -51,6 +50,7 @@ export(soil.units)
 export(soil2netcdf)
 export(soil_params)
 export(soil_process)
+export(soilgrids_soilC_extract)
 export(subset_layer)
 export(to.Tag)
 export(to.TreeCode)
@@ -58,5 +58,15 @@ export(write_ic)
 export(write_veg)
 importFrom(magrittr,"%>%")
 importFrom(ncdf4,ncvar_get)
+importFrom(reshape2,melt)
 importFrom(rlang,"%||%")
 importFrom(rlang,.data)
+importFrom(stats,optim)
+importFrom(stats,qgamma)
+importFrom(terra,extract)
+importFrom(terra,project)
+importFrom(terra,rast)
+importFrom(terra,vect)
+importFrom(utils,flush.console)
+importFrom(utils,setTxtProgressBar)
+importFrom(utils,txtProgressBar)
diff --git a/modules/data.land/man/ocs_extract.Rd b/modules/data.land/man/ocs_extract.Rd
deleted file mode 100644
index 37b3bb74025..00000000000
--- a/modules/data.land/man/ocs_extract.Rd
+++ /dev/null
@@ -1,27 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/soilgrids_soc_extraction.R
-\name{ocs_extract}
-\alias{ocs_extract}
-\title{ocs_extract}
-\usage{
-ocs_extract(site_info, outdir = NULL, verbose = TRUE, ...)
-}
-\arguments{
-\item{site_info}{A dataframe of site info containing the BETYdb site ID, 
-site name, latitude, and longitude, e.g. 
-(site_id, site_name, lat, lon)}
-
-\item{outdir}{Optional. Provide the results as an Rdata file 
-(soilgrids_soc_data.RData)}
-
-\item{verbose}{Provide progress feedback to the terminal? TRUE/FALSE}
-}
-\description{
-ocs_extract function
-A function to extract total soil organic carbon for a single or group of 
-lat/long locationsbased on user-defined site location from SoilGrids250m 
-version 2.0 : https://soilgrids.org
-}
-\author{
-Qianyu Li, Shawn P. Serbin
-}
diff --git a/modules/data.land/man/soilgrids_soilC_extract.Rd b/modules/data.land/man/soilgrids_soilC_extract.Rd
new file mode 100644
index 00000000000..bb615afb9c9
--- /dev/null
+++ b/modules/data.land/man/soilgrids_soilC_extract.Rd
@@ -0,0 +1,65 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/soilgrids_soc_extraction.R
+\name{soilgrids_soilC_extract}
+\alias{soilgrids_soilC_extract}
+\title{soilgrids_soilC_extract}
+\usage{
+soilgrids_soilC_extract(site_info, outdir = NULL, verbose = TRUE, ...)
+}
+\arguments{
+\item{site_info}{A dataframe of site info containing the BETYdb site ID, 
+site name, latitude, and longitude, e.g. 
+(site_id, site_name, lat, lon)}
+
+\item{outdir}{Optional. Provide the results as a CSV file 
+(soilgrids_soilC_data.csv)}
+
+\item{verbose}{Provide progress feedback to the terminal? TRUE/FALSE
+
+##' @examples
+\dontrun{
+
+# Example 1 - using the modex.bnl.gov BETYdb and site IDs to extract data
+db <- 'betydb'
+host_db <- 'modex.bnl.gov'
+db_port <- '5432'
+db_user <- 'bety'
+db_password <- 'bety'
+
+bety <- list(user='bety', password='bety', host='modex.bnl.gov',
+dbname='betydb', driver=RPostgres::Postgres(),write=FALSE)
+
+con <- DBI::dbConnect(drv=bety$driver, dbname=bety$dbname, host=bety$host, 
+password=bety$password, user=bety$user)
+
+suppressWarnings(site_qry <- glue::glue_sql("SELECT *, ST_X(ST_CENTROID(geometry)) AS lon,
+ST_Y(ST_CENTROID(geometry)) AS lat FROM sites WHERE id IN ({ids*})",
+ids = c("676","622","678","766","764"), .con = con))
+
+suppressWarnings(qry_results.1 <- DBI::dbSendQuery(con,site_qry))
+suppressWarnings(qry_results.2 <- DBI::dbFetch(qry_results.1))
+DBI::dbClearResult(qry_results.1)
+dbDisconnect(con)
+
+site_info <- qry_results.2
+verbose <- TRUE
+system.time(result_soc <- soilgrids_soilC_extract(site_info=site_info, verbose=verbose))
+result_soc
+
+}}
+}
+\value{
+a dataframe containing the total soil carbon values  
+and the corresponding standard deviation values (uncertainties) for each location 
+Output column names are c("Site_ID","Site_Name","Latitude","Longitude",
+"Total_soilC","Std_soilC")
+}
+\description{
+soilgrids_soilC_extract function
+A function to extract total soil organic carbon for a single or group of 
+lat/long locationsbased on user-defined site location from SoilGrids250m 
+version 2.0 : https://soilgrids.org
+}
+\author{
+Qianyu Li, Shawn P. Serbin
+}

From 5489fa5d10592b1e0dba8aa4479768f96914b7f3 Mon Sep 17 00:00:00 2001
From: Li <qli1@modex.bnl.gov>
Date: Thu, 29 Sep 2022 13:49:05 -0400
Subject: [PATCH 05/14] Update pecan depends due to a new library

---
 docker/depends/pecan.depends.R                 | 1 +
 modules/data.land/R/soilgrids_soc_extraction.R | 4 ++--
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/docker/depends/pecan.depends.R b/docker/depends/pecan.depends.R
index ae40f90dc5e..eb725ee81c5 100644
--- a/docker/depends/pecan.depends.R
+++ b/docker/depends/pecan.depends.R
@@ -122,6 +122,7 @@ wanted <- c(
 'stringi',
 'stringr',
 'swfscMisc',
+'terra',
 'testthat',
 'tibble',
 'tictoc',
diff --git a/modules/data.land/R/soilgrids_soc_extraction.R b/modules/data.land/R/soilgrids_soc_extraction.R
index 71517ece3cb..1f4b5223461 100644
--- a/modules/data.land/R/soilgrids_soc_extraction.R
+++ b/modules/data.land/R/soilgrids_soc_extraction.R
@@ -12,7 +12,7 @@
 ##' (soilgrids_soilC_data.csv)
 ##' @param verbose Provide progress feedback to the terminal? TRUE/FALSE
 ##' 
-##' ##' @examples
+##' @examples
 ##' \dontrun{
 ##' 
 ##' # Example 1 - using the modex.bnl.gov BETYdb and site IDs to extract data
@@ -56,7 +56,7 @@
 ##' @export
 ##' @author Qianyu Li, Shawn P. Serbin
 ##' 
-soilgrids_soilC_extract <- function (site_info, outdir=NULL, verbose=TRUE, ...) {
+soilgrids_soilC_extract <- function (site_info, outdir=NULL, verbose=TRUE) {
 
     if (is.null(site_info)) {
     stop("Please provide a BETY DB site list containing at least the site id and PostGIS geometry\

From 34873c2bcb4ccfa672b311010f3517210b79eb70 Mon Sep 17 00:00:00 2001
From: Li <qli1@modex.bnl.gov>
Date: Wed, 5 Oct 2022 00:54:52 -0400
Subject: [PATCH 06/14] Remove dependence of reshaping2 package

---
 .../data.land/R/soilgrids_soc_extraction.R    | 46 ++++++++-----------
 1 file changed, 18 insertions(+), 28 deletions(-)

diff --git a/modules/data.land/R/soilgrids_soc_extraction.R b/modules/data.land/R/soilgrids_soc_extraction.R
index 1f4b5223461..37f77aef508 100644
--- a/modules/data.land/R/soilgrids_soc_extraction.R
+++ b/modules/data.land/R/soilgrids_soc_extraction.R
@@ -43,11 +43,7 @@
 ##' result_soc
 ##' 
 ##' }
-##' @importFrom reshape2 melt
-##' @importFrom terra vect project rast extract
-##' @importFrom stats qgamma optim
-##' @importFrom utils flush.console setTxtProgressBar txtProgressBar
-##' 
+
 ##' @return a dataframe containing the total soil carbon values  
 ##' and the corresponding standard deviation values (uncertainties) for each location 
 ##' Output column names are c("Site_ID","Site_Name","Latitude","Longitude",
@@ -59,7 +55,7 @@
 soilgrids_soilC_extract <- function (site_info, outdir=NULL, verbose=TRUE) {
 
     if (is.null(site_info)) {
-    stop("Please provide a BETY DB site list containing at least the site id and PostGIS geometry\
+      PEcAn.logger::logger.error("No site information found. Please provide a BETY DB site list containing at least the site id and PostGIS geometry\
     as lon and lat")
   }
   
@@ -67,8 +63,6 @@ soilgrids_soilC_extract <- function (site_info, outdir=NULL, verbose=TRUE) {
   internal_site_info <- data.frame(site_info$id, site_info$sitename, site_info$lat,site_info$lon)
   #create a variable to store mean and quantile of organic carbon density (ocd) for each soil depth
   ocdquant <- matrix(NA, nrow = 6, ncol = length(internal_site_info$site_info.lon) * 4) #row represents soil depth, col represents mean, 5%, 50% and 95%-quantile of ocd for all sites 
-  #create a variable for site ID
-  siteid <- matrix(rep(as.numeric(site_info$id),times=4*6), nrow=6,byrow=TRUE)
   lonlat <- cbind(internal_site_info$site_info.lon, internal_site_info$site_info.lat)
   base_data_url <- "/vsicurl?max_retry=30&retry_delay=60&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/ocd/ocd_"
   depths <- c("0-5cm", "5-15cm", "15-30cm", "30-60cm", "60-100cm", "100-200cm")
@@ -110,7 +104,7 @@ soilgrids_soilC_extract <- function (site_info, outdir=NULL, verbose=TRUE) {
     if (verbose) {
       utils::setTxtProgressBar(pb, j)
       j <- j+1
-      flush.console()}
+      utils::flush.console()}
   
   # cleanup interim results
     rm(ocd_mean.url, ocd_Q0.05.url, ocd_Q0.50.url, ocd_Q0.95.url, 
@@ -118,30 +112,23 @@ soilgrids_soilC_extract <- function (site_info, outdir=NULL, verbose=TRUE) {
        ocd_mean_real, ocd_Q0.05_real, ocd_Q0.50_real, ocd_Q0.95_real)
   }
   
-  rownames(ocdquant) <- depths
-  colnames(ocdquant) <- c(rep("Mean", length(internal_site_info$site_info.lon)),
-                          rep("0.05", length(internal_site_info$site_info.lon)),
-                          rep("0.5", length(internal_site_info$site_info.lon)),
-                          rep("0.95", length(internal_site_info$site_info.lon)))
-  rownames(siteid) <- depths
-  colnames(siteid) <- c(rep("Mean", length(internal_site_info$site_info.lon)),
-                        rep("0.05", length(internal_site_info$site_info.lon)),
-                        rep("0.5", length(internal_site_info$site_info.lon)),
-                        rep("0.95", length(internal_site_info$site_info.lon)))
+
   if (verbose) {
     close(pb)
   }
   
-  # parse extracted data and prepare for output
-  ocd_fit <- reshape2::melt(ocdquant, id.vars = c("Mean"))
-  id_fit <- reshape2::melt(siteid, id.vars = c("Mean"))
-  colnames(ocd_fit) <- c("Depth", "Quantile", "Value")
-  ocd_fit$Variable <- rep("ocd", length(nrow(ocd_fit)))
-  ocd_fit$siteid <- id_fit$value
-  f1<-factor(ocd_fit$siteid,levels=unique(ocd_fit$siteid))
-  f2<-factor(ocd_fit$Depth,levels=unique(ocd_fit$Depth))
+ # parse extracted data and prepare for output
+  quantile_name <-c(paste("Mean_",site_info$id,sep=""),paste("0.05_",site_info$id,sep=""),paste("0.5_",site_info$id,sep=""),paste("0.95_",site_info$id,sep=""))
+  colnames(ocdquant) <- quantile_name
+  ocdquant_dep <- cbind(ocdquant,depths)
+  ocdquant_long <- tidyr::pivot_longer(as.data.frame(ocdquant_dep),cols=all_of(quantile_name))
+  ocd_df <- tidyr::separate(ocdquant_long,name, c('Quantile', 'siteid'),"_")
+  colnames(ocd_df) <- c("Depth","Quantile", "Siteid","Value")
+  ocd_df$Value<-as.numeric(ocd_df$Value)
+  f1<-factor(ocd_df$Siteid,levels=unique(ocd_df$Siteid))
+  f2<-factor(ocd_df$Depth,levels=unique(ocd_df$Depth))
   #split data by groups of sites and soil depth, while keeping the original order of each group
-  dat <- split(ocd_fit, list(f1, f2))  
+  dat <- split(ocd_df, list(f1, f2))  
   
   #assume the ocd profile follows gamma distribution best
   cgamma <- function(theta, val, stat) {
@@ -203,6 +190,9 @@ soilgrids_soilC_extract <- function (site_info, outdir=NULL, verbose=TRUE) {
     PEcAn.logger::logger.info(paste0("Storing results in: ",file.path(outdir,"soilgrids_soilC_data.csv")))
     utils::write.csv(soilgrids_soilC_data,file=file.path(outdir,"soilgrids_soilC_data.csv"),row.names = FALSE)
   }
+  else {
+    PEcAn.logger::logger.error("No output directory found.")
+  }
   # return the results to the terminal as well
   return(soilgrids_soilC_data)
 }

From a22b201af5b94368bfdf4ac8b82355ecdb3b01f2 Mon Sep 17 00:00:00 2001
From: Li <qli1@modex.bnl.gov>
Date: Wed, 5 Oct 2022 15:32:06 -0400
Subject: [PATCH 07/14] Update document files

---
 modules/data.land/NAMESPACE                   | 10 ------
 .../data.land/man/soilgrids_soilC_extract.Rd  | 32 +++++++++----------
 2 files changed, 16 insertions(+), 26 deletions(-)

diff --git a/modules/data.land/NAMESPACE b/modules/data.land/NAMESPACE
index b01fd6ef5fc..27380323e31 100644
--- a/modules/data.land/NAMESPACE
+++ b/modules/data.land/NAMESPACE
@@ -57,15 +57,5 @@ export(write_ic)
 export(write_veg)
 importFrom(magrittr,"%>%")
 importFrom(ncdf4,ncvar_get)
-importFrom(reshape2,melt)
 importFrom(rlang,"%||%")
 importFrom(rlang,.data)
-importFrom(stats,optim)
-importFrom(stats,qgamma)
-importFrom(terra,extract)
-importFrom(terra,project)
-importFrom(terra,rast)
-importFrom(terra,vect)
-importFrom(utils,flush.console)
-importFrom(utils,setTxtProgressBar)
-importFrom(utils,txtProgressBar)
diff --git a/modules/data.land/man/soilgrids_soilC_extract.Rd b/modules/data.land/man/soilgrids_soilC_extract.Rd
index bb615afb9c9..a78f02956de 100644
--- a/modules/data.land/man/soilgrids_soilC_extract.Rd
+++ b/modules/data.land/man/soilgrids_soilC_extract.Rd
@@ -4,7 +4,7 @@
 \alias{soilgrids_soilC_extract}
 \title{soilgrids_soilC_extract}
 \usage{
-soilgrids_soilC_extract(site_info, outdir = NULL, verbose = TRUE, ...)
+soilgrids_soilC_extract(site_info, outdir = NULL, verbose = TRUE)
 }
 \arguments{
 \item{site_info}{A dataframe of site info containing the BETYdb site ID, 
@@ -14,9 +14,21 @@ site name, latitude, and longitude, e.g.
 \item{outdir}{Optional. Provide the results as a CSV file 
 (soilgrids_soilC_data.csv)}
 
-\item{verbose}{Provide progress feedback to the terminal? TRUE/FALSE
-
-##' @examples
+\item{verbose}{Provide progress feedback to the terminal? TRUE/FALSE}
+}
+\value{
+a dataframe containing the total soil carbon values  
+and the corresponding standard deviation values (uncertainties) for each location 
+Output column names are c("Site_ID","Site_Name","Latitude","Longitude",
+"Total_soilC","Std_soilC")
+}
+\description{
+soilgrids_soilC_extract function
+A function to extract total soil organic carbon for a single or group of 
+lat/long locationsbased on user-defined site location from SoilGrids250m 
+version 2.0 : https://soilgrids.org
+}
+\examples{
 \dontrun{
 
 # Example 1 - using the modex.bnl.gov BETYdb and site IDs to extract data
@@ -46,19 +58,7 @@ verbose <- TRUE
 system.time(result_soc <- soilgrids_soilC_extract(site_info=site_info, verbose=verbose))
 result_soc
 
-}}
-}
-\value{
-a dataframe containing the total soil carbon values  
-and the corresponding standard deviation values (uncertainties) for each location 
-Output column names are c("Site_ID","Site_Name","Latitude","Longitude",
-"Total_soilC","Std_soilC")
 }
-\description{
-soilgrids_soilC_extract function
-A function to extract total soil organic carbon for a single or group of 
-lat/long locationsbased on user-defined site location from SoilGrids250m 
-version 2.0 : https://soilgrids.org
 }
 \author{
 Qianyu Li, Shawn P. Serbin

From 57d0648c507d0a0cde2a2c7e7d9f006d3ef7362f Mon Sep 17 00:00:00 2001
From: Li <qli1@modex.bnl.gov>
Date: Thu, 6 Oct 2022 00:42:48 -0400
Subject: [PATCH 08/14] Add namespace

---
 modules/data.land/R/soilgrids_soc_extraction.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/modules/data.land/R/soilgrids_soc_extraction.R b/modules/data.land/R/soilgrids_soc_extraction.R
index 37f77aef508..1fdd2f3de8a 100644
--- a/modules/data.land/R/soilgrids_soc_extraction.R
+++ b/modules/data.land/R/soilgrids_soc_extraction.R
@@ -121,7 +121,7 @@ soilgrids_soilC_extract <- function (site_info, outdir=NULL, verbose=TRUE) {
   quantile_name <-c(paste("Mean_",site_info$id,sep=""),paste("0.05_",site_info$id,sep=""),paste("0.5_",site_info$id,sep=""),paste("0.95_",site_info$id,sep=""))
   colnames(ocdquant) <- quantile_name
   ocdquant_dep <- cbind(ocdquant,depths)
-  ocdquant_long <- tidyr::pivot_longer(as.data.frame(ocdquant_dep),cols=all_of(quantile_name))
+  ocdquant_long <- tidyr::pivot_longer(as.data.frame(ocdquant_dep),cols=tidyselect::all_of(quantile_name))
   ocd_df <- tidyr::separate(ocdquant_long,name, c('Quantile', 'siteid'),"_")
   colnames(ocd_df) <- c("Depth","Quantile", "Siteid","Value")
   ocd_df$Value<-as.numeric(ocd_df$Value)

From bc66060edd8247b6c68ad57e204589926d959dbe Mon Sep 17 00:00:00 2001
From: "Shawn P. Serbin" <sserbin@bnl.gov>
Date: Thu, 6 Oct 2022 08:48:43 -0400
Subject: [PATCH 09/14] Removed the reshape2 dependency in the data.land
 DESCRIPTION

---
 modules/data.land/DESCRIPTION | 1 -
 1 file changed, 1 deletion(-)

diff --git a/modules/data.land/DESCRIPTION b/modules/data.land/DESCRIPTION
index 642b9c19a3d..247154b83d3 100644
--- a/modules/data.land/DESCRIPTION
+++ b/modules/data.land/DESCRIPTION
@@ -45,7 +45,6 @@ Imports:
     PEcAn.utils,
     PEcAn.visualization,
     purrr,
-    reshape2,
     rjags,
     rlang,
     sf,

From 635cf6508bc76d96434d602267e656d4965957b7 Mon Sep 17 00:00:00 2001
From: "Shawn P. Serbin" <sserbin@bnl.gov>
Date: Thu, 6 Oct 2022 10:11:56 -0400
Subject: [PATCH 10/14] Added tidyr and tidyselect to imports in
 modules/data.land/DESCRIPTION to fix build check error

---
 modules/data.land/DESCRIPTION | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/modules/data.land/DESCRIPTION b/modules/data.land/DESCRIPTION
index 247154b83d3..64a2ee5e20e 100644
--- a/modules/data.land/DESCRIPTION
+++ b/modules/data.land/DESCRIPTION
@@ -52,6 +52,8 @@ Imports:
     sp,
     stringr,
     terra,
+    tidyr,
+    tidyselect,
     traits,
     XML (>= 3.98-1.4)      
 Suggests:

From 72fd336ab0e0b2617d59150eca6eaa3ca28398c7 Mon Sep 17 00:00:00 2001
From: Li <qli1@modex.bnl.gov>
Date: Thu, 6 Oct 2022 15:11:53 -0400
Subject: [PATCH 11/14] Update parameter in pivot_longer function

---
 modules/data.land/R/soilgrids_soc_extraction.R | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/modules/data.land/R/soilgrids_soc_extraction.R b/modules/data.land/R/soilgrids_soc_extraction.R
index 1fdd2f3de8a..6e87ac5f14a 100644
--- a/modules/data.land/R/soilgrids_soc_extraction.R
+++ b/modules/data.land/R/soilgrids_soc_extraction.R
@@ -121,8 +121,7 @@ soilgrids_soilC_extract <- function (site_info, outdir=NULL, verbose=TRUE) {
   quantile_name <-c(paste("Mean_",site_info$id,sep=""),paste("0.05_",site_info$id,sep=""),paste("0.5_",site_info$id,sep=""),paste("0.95_",site_info$id,sep=""))
   colnames(ocdquant) <- quantile_name
   ocdquant_dep <- cbind(ocdquant,depths)
-  ocdquant_long <- tidyr::pivot_longer(as.data.frame(ocdquant_dep),cols=tidyselect::all_of(quantile_name))
-  ocd_df <- tidyr::separate(ocdquant_long,name, c('Quantile', 'siteid'),"_")
+  ocd_df <- tidyr::pivot_longer(as.data.frame(ocdquant_dep),cols=tidyselect::all_of(quantile_name),names_to=c("Quantile", "Siteid"),names_sep = "_")
   colnames(ocd_df) <- c("Depth","Quantile", "Siteid","Value")
   ocd_df$Value<-as.numeric(ocd_df$Value)
   f1<-factor(ocd_df$Siteid,levels=unique(ocd_df$Siteid))

From 0b46e8baacc64b5bf0f26f400711657491caec4b Mon Sep 17 00:00:00 2001
From: "Shawn P. Serbin" <sserbin@bnl.gov>
Date: Fri, 7 Oct 2022 08:48:13 -0400
Subject: [PATCH 12/14] Added additional dependency. Small update to
 soilgrids_soc_extraction.R example doc

---
 modules/data.land/R/soilgrids_soc_extraction.R   | 10 +++++-----
 modules/data.land/man/soilgrids_soilC_extract.Rd |  6 +++---
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/modules/data.land/R/soilgrids_soc_extraction.R b/modules/data.land/R/soilgrids_soc_extraction.R
index 1fdd2f3de8a..361ffb3b0c1 100644
--- a/modules/data.land/R/soilgrids_soc_extraction.R
+++ b/modules/data.land/R/soilgrids_soc_extraction.R
@@ -22,7 +22,7 @@
 ##' db_user <- 'bety'
 ##' db_password <- 'bety'
 ##' 
-##' bety <- list(user='bety', password='bety', host='modex.bnl.gov',
+##' bety <- list(user='bety', password='bety', host=host_db,
 ##' dbname='betydb', driver=RPostgres::Postgres(),write=FALSE)
 ##' 
 ##' con <- DBI::dbConnect(drv=bety$driver, dbname=bety$dbname, host=bety$host, 
@@ -35,11 +35,11 @@
 ##' suppressWarnings(qry_results.1 <- DBI::dbSendQuery(con,site_qry))
 ##' suppressWarnings(qry_results.2 <- DBI::dbFetch(qry_results.1))
 ##' DBI::dbClearResult(qry_results.1)
-##' dbDisconnect(con)
+##' DBI::dbDisconnect(con)
 ##' 
 ##' site_info <- qry_results.2
 ##' verbose <- TRUE
-##' system.time(result_soc <- soilgrids_soilC_extract(site_info=site_info, verbose=verbose))
+##' system.time(result_soc <- PEcAn.data.land::soilgrids_soilC_extract(site_info=site_info, verbose=verbose))
 ##' result_soc
 ##' 
 ##' }
@@ -138,7 +138,7 @@ soilgrids_soilC_extract <- function (site_info, outdir=NULL, verbose=TRUE) {
       pred["Mean"] <- theta[1] / theta[2]
     }
     qstat <- as.numeric(stat)[!is.na(as.numeric(stat))]
-    pred[as.character(qstat)] <- qgamma(qstat, theta[1], theta[2])
+    pred[as.character(qstat)] <- stats::qgamma(qstat, theta[1], theta[2])
     return(sum((pred - val) ^ 2))
   }
   
@@ -147,7 +147,7 @@ soilgrids_soilC_extract <- function (site_info, outdir=NULL, verbose=TRUE) {
     stat = as.character(x$Quantile)
     theta = c(10, 10)
     fit <-
-      list(Gamma = optim(theta, cgamma, val = val, stat = stat))
+      list(Gamma = stats::optim(theta, cgamma, val = val, stat = stat))
     SS <- sapply(fit, function(f) {
       f$value
     })
diff --git a/modules/data.land/man/soilgrids_soilC_extract.Rd b/modules/data.land/man/soilgrids_soilC_extract.Rd
index a78f02956de..11a98add473 100644
--- a/modules/data.land/man/soilgrids_soilC_extract.Rd
+++ b/modules/data.land/man/soilgrids_soilC_extract.Rd
@@ -38,7 +38,7 @@ db_port <- '5432'
 db_user <- 'bety'
 db_password <- 'bety'
 
-bety <- list(user='bety', password='bety', host='modex.bnl.gov',
+bety <- list(user='bety', password='bety', host=host_db,
 dbname='betydb', driver=RPostgres::Postgres(),write=FALSE)
 
 con <- DBI::dbConnect(drv=bety$driver, dbname=bety$dbname, host=bety$host, 
@@ -51,11 +51,11 @@ ids = c("676","622","678","766","764"), .con = con))
 suppressWarnings(qry_results.1 <- DBI::dbSendQuery(con,site_qry))
 suppressWarnings(qry_results.2 <- DBI::dbFetch(qry_results.1))
 DBI::dbClearResult(qry_results.1)
-dbDisconnect(con)
+DBI::dbDisconnect(con)
 
 site_info <- qry_results.2
 verbose <- TRUE
-system.time(result_soc <- soilgrids_soilC_extract(site_info=site_info, verbose=verbose))
+system.time(result_soc <- PEcAn.data.land::soilgrids_soilC_extract(site_info=site_info, verbose=verbose))
 result_soc
 
 }

From 3924e71ec3a2870fd0500244c439e46fc130a922 Mon Sep 17 00:00:00 2001
From: "Shawn P. Serbin" <sserbin@bnl.gov>
Date: Fri, 7 Oct 2022 09:24:40 -0400
Subject: [PATCH 13/14] Fixed doc line lenght issue

---
 modules/data.land/R/soilgrids_soc_extraction.R   | 3 ++-
 modules/data.land/man/soilgrids_soilC_extract.Rd | 3 ++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/modules/data.land/R/soilgrids_soc_extraction.R b/modules/data.land/R/soilgrids_soc_extraction.R
index 0aed7615823..bb1f6810133 100644
--- a/modules/data.land/R/soilgrids_soc_extraction.R
+++ b/modules/data.land/R/soilgrids_soc_extraction.R
@@ -39,7 +39,8 @@
 ##' 
 ##' site_info <- qry_results.2
 ##' verbose <- TRUE
-##' system.time(result_soc <- PEcAn.data.land::soilgrids_soilC_extract(site_info=site_info, verbose=verbose))
+##' system.time(result_soc <- PEcAn.data.land::soilgrids_soilC_extract(site_info=site_info, 
+##' verbose=verbose))
 ##' result_soc
 ##' 
 ##' }
diff --git a/modules/data.land/man/soilgrids_soilC_extract.Rd b/modules/data.land/man/soilgrids_soilC_extract.Rd
index 11a98add473..175dbc71ee5 100644
--- a/modules/data.land/man/soilgrids_soilC_extract.Rd
+++ b/modules/data.land/man/soilgrids_soilC_extract.Rd
@@ -55,7 +55,8 @@ DBI::dbDisconnect(con)
 
 site_info <- qry_results.2
 verbose <- TRUE
-system.time(result_soc <- PEcAn.data.land::soilgrids_soilC_extract(site_info=site_info, verbose=verbose))
+system.time(result_soc <- PEcAn.data.land::soilgrids_soilC_extract(site_info=site_info, 
+verbose=verbose))
 result_soc
 
 }

From cad65fe3bf1532ee11db33074c454a456113f458 Mon Sep 17 00:00:00 2001
From: "Shawn P. Serbin" <sserbin@bnl.gov>
Date: Fri, 7 Oct 2022 11:13:09 -0400
Subject: [PATCH 14/14] Updated modules/data.land/tests/Rcheck_reference.log to
 remove final build check error

---
 modules/data.land/tests/Rcheck_reference.log | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/modules/data.land/tests/Rcheck_reference.log b/modules/data.land/tests/Rcheck_reference.log
index 426e5c1cc4a..5d741d257c5 100644
--- a/modules/data.land/tests/Rcheck_reference.log
+++ b/modules/data.land/tests/Rcheck_reference.log
@@ -10,12 +10,12 @@
 * checking package namespace information ... OK
 * checking package dependencies ... WARNING
 
-Imports includes 30 non-default packages.
+Imports includes 33 non-default packages.
 Importing from so many packages makes the package vulnerable to any of
 them becoming unavailable.  Move as many as possible to Suggests and
 use conditionally.
 * checking package dependencies ... NOTE
-Imports includes 30 non-default packages.
+Imports includes 33 non-default packages.
 Importing from so many packages makes the package vulnerable to any of
 them becoming unavailable.  Move as many as possible to Suggests and
 use conditionally.