Merge pull request #3537 from dlebauer/clip_and_move_raster

Add function to data.land that will clip a raster file, mask it, and then write it out to a new file

Author: dlebauer

CHANGELOG.md

@@ -3,7 +3,14 @@ All notable changes are kept in this file. All changes made should be added to t
 `Unreleased`. Once a new release is made this file will be updated to create a new `Unreleased`
 section for the next release.
 
-For more information about this file see also [Keep a Changelog](http://keepachangelog.com/) .
+	For more information about this file see also [Keep a Changelog](http://keepachangelog.com/) .
+
+## Unreleased
+
+### Added
+
+* Add function `clip_and_save_raster_file()` for subsetting rasters to match a polygon of interest (#3537).
+
 
 ## [1.9.0] - 2025-05-25
 

docker/depends/pecan_package_dependencies.csv

@@ -681,6 +681,7 @@
 "withr","*","models/sipnet","Suggests",FALSE
 "withr","*","modules/allometry","Suggests",FALSE
 "withr","*","modules/data.atmosphere","Suggests",FALSE
+"withr","*","modules/data.land","Suggests",FALSE
 "xgboost","*","modules/assim.sequential","Suggests",FALSE
 "XML","*","base/workflow","Imports",FALSE
 "XML","*","models/biocro","Imports",FALSE

modules/data.land/DESCRIPTION

@@ -67,7 +67,8 @@ Suggests:
     redland,
     raster,
     reticulate,
-    testthat (>= 1.0.2)
+    testthat (>= 1.0.2),
+    withr
 Remotes:
     github::ropensci/traits
 License: BSD_3_clause + file LICENSE

modules/data.land/NAMESPACE

@@ -10,6 +10,7 @@ export(Read.IC.info.BADM)
 export(Read_Tucson)
 export(Soilgrids_SoilC_prep)
 export(buildJAGSdata_InventoryRings)
+export(clip_and_save_raster_file)
 export(cohort2pool)
 export(dataone_download)
 export(download.SM_CDS)

modules/data.land/NEWS.md

@@ -1,3 +1,8 @@
+# Unreleased
+
+* Add function `clip_and_save_raster_file()` for subsetting rasters to match a polygon of interest (#3537).
+
+
 # PEcAn.data.land 1.8.2
 - Removed unused parameter `machine` from put_veg_module()
 

modules/data.land/R/clip_and_save_raster_file.R

@@ -0,0 +1,58 @@
+#' Clip and Save a Raster File
+#'
+#' Clips a raster to a polygon bounding box, optionally masks to polygon, and saves the 
+#' output in the same format as the input.
+#'
+#' @param input_path Character. Path to the input raster file.
+#' @param polygon An object or file coercible to a `SpatVector` by `terra::vect()`
+#'   (e.g., an `sf` object, a `SpatVector`, or a file path to a vector dataset).
+#'   used for clipping and masking. Must have a valid CRS.
+#' @param out_path Character. Path to save the processed raster.
+#' @param mask Logical: Should pixels outside the polygon but inside its bounding box
+#'   be masked out (TRUE) or included (FALSE)?
+#' @param overwrite Logical: Replace output file if it already exists?
+#' @return Invisibly, the clipped `SpatRaster` object. The raster is also saved to `out_path`.
+#' @export
+#' @author David LeBauer
+clip_and_save_raster_file <- function(input_path, polygon, out_path, mask = TRUE, overwrite = TRUE) {
+
+  # Check that input and output files have same extension
+  # This function is not designed to convert between raster formats
+  if (tools::file_ext(input_path) != tools::file_ext(out_path)) {
+    PEcAn.logger::logger.error("Input and output files must have the same extension.")
+  }
+
+  rast_in <- terra::rast(input_path)
+
+  # Coerce to SpatVector if not already
+  if (inherits(polygon, "SpatVector")) { # NB passing a SpatVector to terra::vect() fails
+    poly_sv <- polygon
+  } else {
+    poly_sv <- terra::vect(polygon)
+  }
+
+  if (terra::crs(poly_sv) == "") {
+    PEcAn.logger::logger.error("Input polygon must have CRS defined.")
+  }
+
+  # Reproject polygon to raster CRS if different
+  if (!terra::same.crs(poly_sv, rast_in)) {
+    poly_sv <- terra::project(poly_sv, terra::crs(rast_in))
+  }
+
+  rast_crop <- terra::crop(rast_in, poly_sv)
+
+  if (mask) {
+    rast_to_write <- terra::mask(rast_crop, poly_sv)
+  } else {
+    rast_to_write <- rast_crop
+  }
+
+  terra::writeRaster(
+    rast_to_write,
+    filename = out_path, 
+    overwrite = overwrite
+  )
+
+  invisible(rast_to_write)
+}

modules/data.land/man/clip_and_save_raster_file.Rd

@@ -0,0 +1,73 @@
+# helper to create a small test raster
+make_raster <- function(outfile, crs = "EPSG:4326") {
+  r <- terra::rast(matrix(1:16, 4, 4),
+    extent = terra::ext(0, 4, 0, 4),
+    crs = crs
+  )
+  terra::writeRaster(r, outfile, filetype = "GTiff", overwrite = TRUE)
+  return(outfile)
+}
+
+test_that("clip & mask works: output clipped to polygon bbox and masked", {
+  in_r <- withr::local_tempfile(fileext = ".tif")
+  out_f <- withr::local_tempfile(fileext = ".tif")
+
+  make_raster(outfile = in_r)
+
+  poly <- terra::as.polygons(
+    terra::ext(1, 3, 1, 3),
+    crs = "EPSG:4326"
+  )
+
+  clip_and_save_raster_file(input_path = in_r, polygon = poly, out_path = out_f, mask = TRUE)
+
+  expect_true(file.exists(out_f))
+
+  r_out <- terra::rast(out_f)
+  expect_equal(terra::ext(r_out), terra::ext(sf::st_bbox(poly)))
+
+  inside_vals <- terra::values(terra::mask(r_out, poly, inverse = FALSE))
+  expect_true(all(!is.na(inside_vals)))
+
+  outside_vals <- terra::values(terra::mask(r_out, poly, inverse = TRUE))
+  expect_true(all(is.na(outside_vals)))
+})
+
+test_that("clip without mask retains all values within bbox", {
+  in_r <- withr::local_tempfile(fileext = ".tif")
+  make_raster(outfile = in_r)
+
+  poly <- sf::st_as_sf(
+    sf::st_as_sfc(
+      sf::st_bbox(c(xmin = 1, ymin = 1, xmax = 3, ymax = 3), crs = sf::st_crs(4326))
+    )
+  )
+  out_f <- withr::local_tempfile(fileext = ".tif")
+
+  clip_and_save_raster_file(in_r, poly, out_f, mask = FALSE)
+  r_out <- terra::rast(out_f)
+  expect_false(any(is.na(terra::values(r_out))))
+})
+
+test_that("preserves CRS and filetype", {
+  in_r_path <- withr::local_tempfile(fileext = ".tif")
+  make_raster(outfile = in_r_path, crs = "EPSG:3857")
+
+  spatvect_raster <- terra::rast(in_r_path)
+
+  poly <- sf::st_as_sf(
+    sf::st_as_sfc(
+      sf::st_bbox(c(xmin = 1, ymin = 1, xmax = 3, ymax = 3), crs = sf::st_crs(3857))
+    )
+  )
+  out_f_path <- withr::local_tempfile(fileext = ".tif")
+
+  clip_and_save_raster_file(input_path = in_r_path, polygon = poly, out_path = out_f_path)
+  r_out <- terra::rast(out_f_path)
+
+  expect_equal(
+    tools::file_ext(terra::sources(r_out)[1]),
+    tools::file_ext(terra::sources(spatvect_raster)[1])
+  )
+  expect_true(terra::same.crs(r_out, spatvect_raster))
+})

modules/data.land/tests/testthat/test-clip_and_save_raster_file.R

@@ -22,20 +22,28 @@ test_that("Species matching works", {
     user = "bety",
     password = "bety",
     host = "localhost",
-    driver = "Postgres")
-  con <- PEcAn.DB::db.open(db_params)
+    driver = "Postgres"
+  )
+  con <- tryCatch(
+    PEcAn.DB::db.open(db_params),
+    error = function(e) NULL
+  )
+
+  skip_if(is.null(con), "No database connection available for species matching tests.")
 
   test_merge(c('ACRU', 'TSCA'), 'usda', con)
   test_merge(c(316L, 261L), 'fia', con)
   test_merge(c('Acer rubrum', 'Tsuga canadensis'), 'latin_name', con)
 
   test_table <- data.frame(
     bety_species_id = c(30L, 1419L),
-    input_code = c('AceRub', 'TsuCan'))
+    input_code = c('AceRub', 'TsuCan')
+  )
 
   test_merge(
     input_codes = test_table$input_code,
     format_name = 'custom',
     bety = con,
-    translation_table = test_table)
+    translation_table = test_table
+  )
 })