From f7fd1f5ee3a5da5577812d7f23a19750fb608e75 Mon Sep 17 00:00:00 2001 From: Aritra Dey Date: Mon, 24 Feb 2025 22:46:06 +0530 Subject: [PATCH 1/2] feat: refactor SoilGrids functions to handle multiple properties and improve efficiency --- modules/data.land/R/soilgrids_utils.R | 99 +++++++++++++++++++++++++++ 1 file changed, 99 insertions(+) create mode 100644 modules/data.land/R/soilgrids_utils.R diff --git a/modules/data.land/R/soilgrids_utils.R b/modules/data.land/R/soilgrids_utils.R new file mode 100644 index 00000000000..d2e8004fc0e --- /dev/null +++ b/modules/data.land/R/soilgrids_utils.R @@ -0,0 +1,99 @@ +#' Download SoilGrids raster files +#' +#' @param soil_property The soil property to download (e.g., "ocd", "clay", "sand", etc.) +#' @param depth The depth layer to download (e.g., "0-5cm", "5-15cm", etc.) +#' @param quantile The quantile to download (e.g., "mean", "Q0.05", "Q0.5", "Q0.95") +#' @param extent The extent to crop the raster to (terra::ext object) +#' @param outdir The directory to save the downloaded raster +#' @param overwrite Whether to overwrite existing files (default: FALSE) +#' @return The file path of the downloaded raster +#' @export +download_soilgrids_raster <- function(soil_property, depth, quantile, extent, outdir, overwrite = FALSE) { + soil_properties <- if (is.vector(soil_property)) soil_property else c(soil_property) + + for (property in soil_properties) { + vrt_url <- sprintf( + "/vsicurl/https://files.isric.org/soilgrids/latest/data/%s/%s_%s_%s.vrt", + property, property, depth, quantile + ) + + raster <- terra::rast(vrt_url) |> + terra::crop(extent) + + out_file <- file.path(outdir, sprintf("%s_%s_%s.tif", property, depth, quantile)) + + if (!file.exists(out_file) || overwrite) { + terra::writeRaster(raster, out_file, overwrite = TRUE) + } + } + + return(out_file) +} + +#' Query SoilGrids data +#' +#' @param points A spatial object (e.g., sf or terra::vect) containing the points to query +#' @param soil_property The soil property to query (e.g., "ocd", "clay", "sand", etc.) +#' @param depth The depth layer to query (e.g., "0-5cm", "5-15cm", etc.) +#' @param quantile The quantile to query (e.g., "mean", "Q0.05", "Q0.5", "Q0.95") +#' @param local_raster Optional. The file path of a local raster to query. If NULL, queries the remote VRT. +#' @return A data frame containing the queried data +#' @export +query_soilgrids <- function(points, soil_property, depth, quantile, local_raster = NULL) { + soil_properties <- if (is.vector(soil_property)) soil_property else c(soil_property) + + data_list <- lapply(soil_properties, function(property) { + if (is.null(local_raster)) { + vrt_url <- sprintf( + "/vsicurl/https://files.isric.org/soilgrids/latest/data/%s/%s_%s_%s.vrt", + property, property, depth, quantile + ) + raster <- terra::rast(vrt_url) + } else { + raster <- terra::rast(local_raster) + } + terra::extract(raster, points) + }) + + return(do.call(rbind, data_list)) +} + +#' Estimate SoilGrids uncertainties +#' +#' @param data A data frame containing the SoilGrids data +#' @param quantiles A vector of quantiles to estimate uncertainties (default: c(0.05, 0.5, 0.95)) +#' @return A data frame containing the uncertainties +#' @export +soilgrids_uncertainty <- function(data, quantiles = c(0.05, 0.5, 0.95)) { + uncertainties <- data %>% + dplyr::group_by(site_id) %>% + dplyr::summarize( + mean = mean(value, na.rm = TRUE), + sd = sd(value, na.rm = TRUE), + quantile_05 = quantile(value, probs = quantiles[1], na.rm = TRUE), + quantile_50 = quantile(value, probs = quantiles[2], na.rm = TRUE), + quantile_95 = quantile(value, probs = quantiles[3], na.rm = TRUE) + ) + + return(uncertainties) +} + +#' Calculate SoilGrids carbon stocks +#' +#' @param soc_data A data frame containing the SoilGrids SOC data +#' @param depths A vector of depths corresponding to the SOC data +#' @param thickness A vector of thicknesses for each depth layer +#' @param coarse_fraction Optional. A vector of coarse fractions to include in the calculation +#' @return A data frame containing the calculated carbon stocks +#' @export +soilgrids_carbonstocks <- function(soc_data, depths, thickness, coarse_fraction = NULL) { + carbonstocks <- soc_data %>% + dplyr::group_by(site_id) %>% + dplyr::summarize( + total_soc = sum(value * thickness, na.rm = TRUE) + + if (!is.null(coarse_fraction)) sum(coarse_fraction * thickness, na.rm = TRUE) else 0, + total_soc_0_30 = sum(value[depths <= 30] * thickness[depths <= 30], na.rm = TRUE) + ) + + return(carbonstocks) +} From 7649780830b7197beb6da7530e6d29e1402a0a72 Mon Sep 17 00:00:00 2001 From: Aritra Dey Date: Wed, 26 Feb 2025 01:21:04 +0530 Subject: [PATCH 2/2] feat: Add customizable depth segments for SOC calculations in soilgrids_carbonstocks function --- modules/data.land/R/soilgrids_utils.R | 40 +++++++++++++++++---------- 1 file changed, 26 insertions(+), 14 deletions(-) diff --git a/modules/data.land/R/soilgrids_utils.R b/modules/data.land/R/soilgrids_utils.R index d2e8004fc0e..67c2b8253a0 100644 --- a/modules/data.land/R/soilgrids_utils.R +++ b/modules/data.land/R/soilgrids_utils.R @@ -6,9 +6,10 @@ #' @param extent The extent to crop the raster to (terra::ext object) #' @param outdir The directory to save the downloaded raster #' @param overwrite Whether to overwrite existing files (default: FALSE) +#' @param target_crs Optional. The target coordinate reference system (default: WGS1984) #' @return The file path of the downloaded raster #' @export -download_soilgrids_raster <- function(soil_property, depth, quantile, extent, outdir, overwrite = FALSE) { +download_soilgrids_raster <- function(soil_property, depth, quantile, extent, outdir, overwrite = FALSE, target_crs = "EPSG:4326") { soil_properties <- if (is.vector(soil_property)) soil_property else c(soil_property) for (property in soil_properties) { @@ -37,9 +38,10 @@ download_soilgrids_raster <- function(soil_property, depth, quantile, extent, ou #' @param depth The depth layer to query (e.g., "0-5cm", "5-15cm", etc.) #' @param quantile The quantile to query (e.g., "mean", "Q0.05", "Q0.5", "Q0.95") #' @param local_raster Optional. The file path of a local raster to query. If NULL, queries the remote VRT. +#' @param target_crs Optional. The target coordinate reference system (default: WGS1984) #' @return A data frame containing the queried data #' @export -query_soilgrids <- function(points, soil_property, depth, quantile, local_raster = NULL) { +query_soilgrids <- function(points, soil_property, depth, quantile, local_raster = NULL, target_crs = "EPSG:4326") { soil_properties <- if (is.vector(soil_property)) soil_property else c(soil_property) data_list <- lapply(soil_properties, function(property) { @@ -70,9 +72,9 @@ soilgrids_uncertainty <- function(data, quantiles = c(0.05, 0.5, 0.95)) { dplyr::summarize( mean = mean(value, na.rm = TRUE), sd = sd(value, na.rm = TRUE), - quantile_05 = quantile(value, probs = quantiles[1], na.rm = TRUE), - quantile_50 = quantile(value, probs = quantiles[2], na.rm = TRUE), - quantile_95 = quantile(value, probs = quantiles[3], na.rm = TRUE) + # Simulate uncertainty based on mean and quantiles + simulated_values = replicate(1000, rnorm(n(), mean, sd)), # Example simulation + quantiles = sapply(quantiles, function(q) quantile(simulated_values, probs = q, na.rm = TRUE)) ) return(uncertainties) @@ -83,17 +85,27 @@ soilgrids_uncertainty <- function(data, quantiles = c(0.05, 0.5, 0.95)) { #' @param soc_data A data frame containing the SoilGrids SOC data #' @param depths A vector of depths corresponding to the SOC data #' @param thickness A vector of thicknesses for each depth layer +#' @param segments A list of depth segments (e.g., list(c(0, 30), c(30, 70), c(70, 150))) #' @param coarse_fraction Optional. A vector of coarse fractions to include in the calculation -#' @return A data frame containing the calculated carbon stocks +#' @return A data frame containing the calculated carbon stocks for specified segments #' @export -soilgrids_carbonstocks <- function(soc_data, depths, thickness, coarse_fraction = NULL) { - carbonstocks <- soc_data %>% - dplyr::group_by(site_id) %>% - dplyr::summarize( - total_soc = sum(value * thickness, na.rm = TRUE) + - if (!is.null(coarse_fraction)) sum(coarse_fraction * thickness, na.rm = TRUE) else 0, - total_soc_0_30 = sum(value[depths <= 30] * thickness[depths <= 30], na.rm = TRUE) - ) +soilgrids_carbonstocks <- function(soc_data, depths, thickness, segments, coarse_fraction = NULL) { + carbonstocks <- data.frame(site_id = unique(soc_data$site_id)) + + for (segment in segments) { + lower_bound <- segment[1] + upper_bound <- segment[2] + + total_soc <- soc_data %>% + dplyr::filter(depths >= lower_bound & depths < upper_bound) %>% + dplyr::summarize( + total_soc = sum(value * thickness[depths >= lower_bound & depths < upper_bound], na.rm = TRUE) + + if (!is.null(coarse_fraction)) sum(coarse_fraction * thickness[depths >= lower_bound & depths < upper_bound], na.rm = TRUE) else 0 + ) + + carbonstocks <- dplyr::bind_cols(carbonstocks, total_soc) + colnames(carbonstocks)[ncol(carbonstocks)] <- paste0("total_soc_", lower_bound, "_", upper_bound) + } return(carbonstocks) }