aggregate.R created

CHANGELOG.md

@@ -45,6 +45,7 @@ For more information about this file see also [Keep a Changelog](http://keepacha
 - Added new feature of preparing initial conditions for MODIS LAI, AGB, ISCN SOC, and soil moisture across NA anchor sites.
 - Added GEDI AGB preparation workflow.
 - Added new feature of downloading datasets from the NASA DAAC ORNL database.
+- Added 'aggregate' as a new feature for downscaled data
 
 ### Fixed
 

CITATION.cff

@@ -120,7 +120,9 @@ authors:
   - given-names: Eric R. Scott
     affiliation: University of Arizona
     orcid: 'https://orcid.org/0000-0002-7430-7879'
-
+  - given-names: Harunobu Ishii
+    affiliation: Boston University Software & Application Innovation Lab(SAIL)
+
 preferred-citation:
   type: article
   title: Facilitating feedbacks between field measurements and ecosystem models

docker/depends/pecan_package_dependencies.csv

@@ -61,6 +61,7 @@
 "dplyr",">= 1.1.2","base/db","Imports",FALSE
 "ellipse","*","modules/assim.batch","Imports",FALSE
 "emdbook","*","modules/assim.sequential","Suggests",FALSE
+"exactextractr","*","modules/assim.sequential","Suggests",FALSE
 "foreach","*","base/remote","Imports",FALSE
 "foreach","*","modules/data.atmosphere","Suggests",FALSE
 "foreach","*","modules/data.remote","Imports",FALSE

modules/assim.sequential/DESCRIPTION

@@ -32,6 +32,7 @@ Imports:
     stringr
 Suggests:
     corrplot,
+    exactextractr,
     ggrepel,
     emdbook,
     glue,

modules/assim.sequential/NAMESPACE

@@ -21,6 +21,7 @@ export(SDA_control)
 export(SDA_remote_launcher)
 export(SDA_timeseries_plot)
 export(adj.ens)
+export(aggregate)
 export(alltocs)
 export(alr)
 export(assessParams)

modules/assim.sequential/R/aggregate.R

@@ -0,0 +1,59 @@
+#' @title Aggregation Function
+#' @name aggregate
+#' @author Harunobu Ishii
+#'
+#' @param downscale_output  Raster file output from downscale_function.R. Read file in this way if stored locally: \code{downscale_output <- readRDS("xxx.rds")}
+#' @param polygon_data A spatial polygon object (e.g., an `sf` object) that defines the spatial units for aggregation. 
+#'                     This data should be in a coordinate reference system compatible with the raster data (e.g., "EPSG:4326").
+#' @param func A character string specifying the aggregation function to use (e.g., 'mean', 'sum').
+#' @details This function will aggregate previously downscaled carbon flux amount to a spatial unit of choice 
+#'
+#' @return It returns the `polygon_data` with added columns for mean and sum values of the aggregated raster data for each ensemble member.
+#' @import sf
+#' @import exactextractr
+#' @import raster
+#' @export
+#' @examples
+#'        \dontrun{
+#'        # Download a shapefile of U.S. (polygon data)
+#'        url <- "https://www2.census.gov/geo/tiger/GENZ2020/shp/cb_2020_us_state_20m.zip"
+#'        download.file(url, destfile = "polygon/us_states.zip")
+#'        
+#'        # Unzip the downloaded file and save locally
+#'        unzip("polygon/us_states.zip", exdir = "polygon/us_states")
+#'        us_states <- st_read("polygon/us_states/cb_2020_us_state_20m.shp")
+#'        saveRDS(us_states, "polygon/us_states.rds")
+#'        
+#'        # Load the saved polygon data with Massachusetts as an example
+#'        us_states <- readRDS("polygon/us_states.rds")
+#'        state <- "MA"
+#'        polygon_data <- st_transform(us_states[us_states$STUSPS == state, ], crs = "EPSG:4326")
+#'      
+#'        # Load the downscaled raster output
+#'        downscale_output <- readRDS("path/to/downscale_output.rds")
+#'        
+#'        # Slot in as argument to the aggregate function
+#'        result <- aggregate(downscale_output, polygon_data)
+#'        print(result)
+#'        }
+
+aggregate <- function(downscale_output, polygon_data, func = 'mean'){
+  grand_TTL <- 0
+  if (sf::st_crs(downscale_output$maps$ensemble1) != sf::st_crs(polygon_data)) {
+    stop("CRS of downscale_output and polygon_data must match.")
+  }
+
+  # Perform spatial operations on each raster
+  for (name in names(downscale_output$maps)) {
+    raster_data <- downscale_output$maps[[name]]
+    agg_values <- exactextractr::exact_extract(raster_data, polygon_data, fun = func)
+
+    polygon_data[[paste0(name, "_", func)]] <- agg_values
+    grand_TTL = grand_TTL + agg_values
+  }
+  if(func == 'mean'){
+    grand_TTL = grand_TTL/length(downscale_output$maps)
+  }
+  polygon_data[[paste0("TTL_", func)]] <- grand_TTL
+  return (polygon_data)
+}

modules/assim.sequential/tests/testthat/test_aggregation.R

@@ -0,0 +1,60 @@
+library(testthat)
+library(sf)
+library(raster)
+library(exactextractr)
+library(terra)
+source("../../R/aggregate.R")
+test_that("returns aggregated values for RI", {
+    # Load the saved polygon data with Massachusetts as an example
+    us_states <- readRDS("test_aggregation/us_states.rds")
+    state <- "RI"
+    polygon_data <- st_transform(us_states[us_states$STUSPS == state, ], crs = "EPSG:4326")
+
+    # Load the downscaled raster output
+    downscale_output <- list(
+      maps = list(
+        ensemble1 = "test_aggregation/ensemble1.tif",
+        ensemble2 = "test_aggregation/ensemble2.tif",
+        ensemble3 = "test_aggregation/ensemble3.tif"
+      )
+    )
+
+    read_raster <- function(file_path) {
+        rast(file_path)
+    }
+
+    downscale_output$maps <- lapply(downscale_output$maps, read_raster)
+    # Aggregate for RI
+    RI <- aggregate(downscale_output, polygon_data, func = 'mean')
+    comp <- RI$TTL_mean * 10^9
+    comparison_result <- (1.31 < comp & comp < 1.32)
+    expect_true(comparison_result)
+})
+
+test_that("returns error of unmatched CRS", {
+  # Load the saved polygon data with Massachusetts as an example
+  us_states <- readRDS("test_aggregation/us_states.rds")
+  state <- "RI"
+  polygon_data <- st_transform(us_states[us_states$STUSPS == state, ], crs = "EPSG:2222")
+
+  # Load the downscaled raster output
+  downscale_output <- list(
+    maps = list(
+      ensemble1 = "test_aggregation/ensemble1.tif",
+      ensemble2 = "test_aggregation/ensemble2.tif",
+      ensemble3 = "test_aggregation/ensemble3.tif"
+    )
+  )
+
+  read_raster <- function(file_path) {
+    rast(file_path)
+  }
+
+  downscale_output$maps <- lapply(downscale_output$maps, read_raster)
+  expect_error(
+    aggregate(downscale_output, polygon_data, func = 'mean'),
+    "CRS of downscale_output and polygon_data must match."
+  )
+})
+
+