Optimization for Large Files

DESCRIPTION

@@ -9,7 +9,6 @@ License: `use_mit_license()`, `use_gpl3_license()` or friends to pick a
     license
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
 Depends: 
     R (>= 2.10)
 LazyData: true
@@ -23,3 +22,4 @@ Imports:
     magrittr,
     ggplot2
 URL: https://nefsc.github.io/READ_EDAB_Utilities/
+Config/roxygen2/version: 8.0.0

R/convert_2d_longitude_gridded.R

@@ -11,6 +11,9 @@
 #' @export
 convert_2d_longitude_gridded <- function(data.in, write.out = FALSE, output.files = NA) {
 
+  # Boost terra memory limit for in-memory processing fallback
+  terra::terraOptions(memfrac = 0.8)
+
   # 1. Standardized input coercion block
   data.ls = EDABUtilities:::import_data(data.in)
 
@@ -23,18 +26,59 @@ convert_2d_longitude_gridded <- function(data.in, write.out = FALSE, output.file
   out.ls <- lapply(seq_along(data.ls), function(i) {
 
     current_data <- data.ls[[i]]
+    is_file_input <- is.character(current_data)
+
+    # Fast path: File-to-file transformation via GDAL (only for rasters written to disk)
+    if (is_file_input && write.out) {
+      if (!file.exists(current_data)) stop(paste("File does not exist:", current_data))
+
+      # Quickly read extent without loading data into memory
+      temp_rast <- terra::rast(current_data)
+      dat.ext <- terra::ext(temp_rast)
+
+      out_dir <- dirname(output.files[i])
+      if (!dir.exists(out_dir)) dir.create(out_dir, recursive = TRUE)
+
+      if (dat.ext[1] >= -0.001 && dat.ext[2] <= 360.001 && dat.ext[2] > 180.001) {
+        message("Detected longitude range approx 0-360. Fast processing via GDAL warp...")
+
+        # Perform fast GDAL extent shift
+        sf::gdal_utils(
+          util = "warp",
+          source = current_data,
+          destination = output.files[i],
+          options = c(
+            "-t_srs", "EPSG:4326", 
+            "-te", "-180", "-90", "180", "90", 
+            "-wo", "SOURCE_EXTRA=1000", 
+            "--config", "CENTER_LONG", "0"
+          )
+        )
+        return(terra::rast(output.files[i]))
+
+      } else if (!(dat.ext[1] >= -180.001 && dat.ext[2] <= 180.001)) {
+        stop("Longitude out of range. Extent is outside expected boundaries.")
+      } else {
+        message("Already standard format (-180:180). Copying file.")
+        file.copy(current_data, output.files[i], overwrite = TRUE)
+        return(terra::rast(output.files[i]))
+      }
+    }
+
+    # ---------------------------------------------------------
+    # Fallback path: In-memory terra processing
+    # (used if input is already an object, or write.out is FALSE)
+    # ---------------------------------------------------------
 
-    # Safe coercion with file validation
-    if (is.character(current_data)) {
+    if (is_file_input) {
       if (!file.exists(current_data)) stop(paste("File does not exist:", current_data))
       current_data <- current_data |> terra::rast()
     }
 
     dat.ext <- current_data |> terra::ext()
 
-    # Spatial logic checks without redundant variables
     if (dat.ext[1] >= -0.001 && dat.ext[2] <= 360.001 && dat.ext[2] > 180.001) {
-      message("Detected longitude range approximately 0-360. Converting to -180 to +180.")
+      message("Detected longitude range approx 0-360. Converting to -180 to +180 via terra::rotate.")
 
       if (!terra::is.lonlat(current_data)) {
         warning("Object does not have a standard geographic (lat/lon) CRS. Rotation may fail or shift bounds unexpectedly.")

R/crop_nc_2d.R

@@ -15,82 +15,114 @@
 crop_nc_2d <- function(data.in, shp.file, var.name, area.names = NA, write.out = FALSE, output.files = NULL) {
 
   # Data Input Standardization
-  data.ls = EDABUtilities:::import_data(data.in)
+  data.ls = EDABUtilities:::import_data(data.in) 
 
   # Spatial Input Standardization
-  shp.vect = EDABUtilities:::import_shp(shp.file)
-  use.shp = ifelse(class(shp.vect) == 'SpatVector',T,F)
+  shp.vect = EDABUtilities:::import_shp(shp.file) 
+  use.shp = ifelse(class(shp.vect) == 'SpatVector',T,F) 
 
-  if (!use.shp) stop("A valid shp.file must be provided.")
+  if (!use.shp) stop("A valid shp.file must be provided.") 
 
   # Optimized Area Names Filtering
-  if (!is.null(area.names) && !all(is.na(area.names))) {
-    shp.str <- as.data.frame(shp.vect)
+  if (!is.null(area.names) && !all(is.na(area.names))) { 
+    shp.str <- as.data.frame(shp.vect) 
 
-    # Safely find the first column that contains the requested area.names
-    target_col <- NULL
-    for (col in names(shp.str)) {
+    target_col <- NULL 
+    for (col in names(shp.str)) { 
       if (all(area.names %in% shp.str[[col]])) {
-        target_col <- col
-        break
+        target_col <- col 
+        break 
       }
     }
 
-    if (is.null(target_col)) {
-      stop("None of the attributes in shp.file contain all specified area.names.")
+    if (is.null(target_col)) { 
+      stop("None of the attributes in shp.file contain all specified area.names.") 
     }
 
-    # Subset vector directly using terra logic
-    shp.vect <- shp.vect[shp.vect[[target_col]] %in% area.names, ]
+    shp.vect <- shp.vect[shp.vect[[target_col]][,1] %in% area.names , ] 
   }
 
-  out.ls <- list()
+  out.ls <- list() 
 
-  for (i in seq_along(data.ls)) {
+  for (i in seq_along(data.ls)) { 
 
-    if (is.character(data.ls[[i]])) {
-      if (!file.exists(data.ls[[i]])) stop(sprintf("File does not exist: %s", data.ls[[i]]))
-      data.orig <- terra::rast(data.ls[[i]])
+    if (is.character(data.ls[[i]])) { 
+      if (!file.exists(data.ls[[i]])) stop(sprintf("File does not exist: %s", data.ls[[i]])) 
+      data.orig <- terra::rast(data.ls[[i]]) 
     } else {
-      data.orig <- data.ls[[i]]
+      data.orig <- data.ls[[i]] 
     }
+    #Check CRS
+    data.crs = terra::crs(data.orig)
+    shp.crs = terra::crs(shp.vect)
 
-    data.orig <- EDABUtilities::convert_2d_longitude_gridded(data.orig)[[1]]
+    if(!identical(data.crs,shp.crs)){
+      shp.vect = terra::project(shp.vect,data.crs)
+    }
+
+    # ---------------------------------------------------------
+    # 1. Fast Vector Shift (Align Shapefile to Raster)
+    # ---------------------------------------------------------
+    e_rast <- as.vector(terra::ext(data.orig)) 
+    shp_crop <- shp.vect 
 
-    # Edge Case Handling: Numeric Extent intersection check (avoids SpatExtent class mismatch)
-    e1 <- as.vector(terra::ext(data.orig))
-    e2 <- as.vector(terra::ext(shp.vect))
+    # If raster is 0-360 but shapefile has negative longitudes
+    if (e_rast["xmax"] > 180.001 && any(terra::ext(shp_crop)[1:2] < 0)) {
+      # Shift the shapefile geometry 360 degrees East so it overlays on the 0-360 raster
+      shp_crop <- terra::rotate(shp_crop, long = 0,split = T, left = F)
+      #If raster is -180 to 180 but shapefile is 0-360
+    } else if (e_rast["xmax"] <= 180.001 && any(terra::ext(shp_crop)[1:2] > 180)) {
+      # Shift the shapefile geometry 360 degrees West
+      shp_crop <- terra::rotate(shp_crop)
+    }
 
-    # Bounding boxes intersect if they overlap on both axes
-    intersects <- (e1["xmin"] <= e2["xmax"]) && (e1["xmax"] >= e2["xmin"]) && 
-      (e1["ymin"] <= e2["ymax"]) && (e1["ymax"] >= e2["ymin"])
+    # ---------------------------------------------------------
+    # 2. Intersection Check & Cropping
+    # ---------------------------------------------------------
+    # Edge Case Handling: Numeric Extent intersection check (avoids SpatExtent class mismatch)[cite: 2]
+    e1 <- as.vector(terra::ext(data.orig)) 
+    e2 <- as.vector(terra::ext(shp_crop))  
 
-    if (!intersects) {
-      warning(sprintf("Data extent and shapefile extent do not intersect for item %s. Skipping.", i))
-      next
+    # Bounding boxes intersect if they overlap on both axes[cite: 2]
+    intersects <- (e1["xmin"] <= e2["xmax"]) && (e1["xmax"] >= e2["xmin"]) &&  
+      (e1["ymin"] <= e2["ymax"]) && (e1["ymax"] >= e2["ymin"]) 
+
+    if (!intersects) { 
+      warning(sprintf("Data extent and shapefile extent do not intersect for item %s. Skipping.", i)) 
+      next 
     }
 
-    data.crop <- terra::crop(data.orig, shp.vect)
+    # Perform the blazing fast crop on the un-rotated raster
+    data.crop <- terra::crop(data.orig, shp_crop) 
+
+    # ---------------------------------------------------------
+    # 3. Post-Crop Standardize
+    # ---------------------------------------------------------
+    # Convert longitude on the tiny cropped raster instead of the massive global one
+    data.crop <- EDABUtilities::convert_2d_longitude_gridded(data.crop)[[1]]
 
-    if (write.out) {
-      if (is.null(output.files) || length(output.files) != length(data.ls)) {
-        stop("output.files must be provided and match the length of data.in when write.out is TRUE.")
+    # ---------------------------------------------------------
+    # 4. Write Output
+    # ---------------------------------------------------------
+    if (write.out) { 
+      if (is.null(output.files) || length(output.files) != length(data.ls)) { 
+        stop("output.files must be provided and match the length of data.in when write.out is TRUE.") 
       }
 
-      out_dir <- dirname(output.files[i])
-      if (!dir.exists(out_dir)) dir.create(out_dir, recursive = TRUE)
+      out_dir <- dirname(output.files[i]) 
+      if (!dir.exists(out_dir)) dir.create(out_dir, recursive = TRUE) 
 
       terra::writeCDF(data.crop, output.files[i], varname = var.name, overwrite = TRUE)  
     } else {
-      out.ls[[i]] <- data.crop
+      out.ls[[i]] <- data.crop 
     }
   }
 
-  if (write.out == FALSE) {
-    if (is.character(data.in)) {
-      names(out.ls) <- basename(data.in)
+  if (write.out == FALSE) { 
+    if (is.character(data.in)) { 
+      names(out.ls) <- basename(data.in) 
     } else {
-      names(out.ls) <- paste0("layer_", seq_along(out.ls))
+      names(out.ls) <- paste0("layer_", seq_along(out.ls)) 
     }
     return(out.ls)  
   }

R/import_data.R

@@ -3,11 +3,12 @@
 #' This function extracts spatial raster data across specified shapefile regions and aggregates it temporally to produce timeseries summary statistics. It processes inputs by grouping them (e.g., aggregating daily layers into annual time series) and outputs either a list of summarized data frames or writes RDS files directly.
 #'
 #' @param data.in character vector, list, or SpatRaster. Single file path, vector of file paths, single SpatRaster, or list of SpatRasters representing the spatial data.
+#' @param var.name character. Name of variable to be subseted from dataset
 #'
 #' @return list of spatRasters
 #' 
 
-import_data = function(data.in){
+import_data = function(data.in, var.name = NULL){
 
   if (inherits(data.in, "SpatRaster")|inherits(data.in, "SpatRasterDataset")) {
     data.ls <- list(data.in)
@@ -23,5 +24,25 @@ import_data = function(data.in){
     stop("data.in must be a file path, a vector of file paths, a single SpatRaster, or a list of SpatRasters.")
   }
 
+  if(!is.null(var.name)){
+    data.ls = lapply(data.ls,function(x){
+      data.varname = terra::varnames(x)
+      if(length(data.varname > 1)){
+        data = x[[terra::varnames(x) == var.name]]
+      }else{
+
+        data.names = terra::names(data)
+        data.names = data.names[grepl(var.name,data.names)]
+
+        if(length(data.names) ==0){
+          warning('NetCDF file does not contain any fields with var.name=',var.name)
+        }else{
+
+          data = subset(x,data.names)
+        }
+      }
+    })
+  }
+
   return(data.ls)
 }

R/make_2d_anomaly_gridded.R

@@ -49,7 +49,7 @@ make_2d_anomaly_gridded <- function(data.in, climatology, var.name, shp.file = N
     }
 
     # Pre-mask climatology once to avoid doing it N times inside the loop
-    climatology <- terra::mask(climatology[[1]], shp.vect)
+    climatology <- EDABUtilities::crop_nc_2d(climatology[[1]], shp.file = shp.vect,area.names = area.names,var.name = var.name  )[[1]]
   }
 
   out.ls <- list()
@@ -70,7 +70,7 @@ make_2d_anomaly_gridded <- function(data.in, climatology, var.name, shp.file = N
 
     # Align extents and resolutions if mismatched
     if (!(all(terra::res(data) == terra::res(climatology)) && all(terra::ext(data) == terra::ext(climatology)))) {
-      climatology <- terra::crop(climatology,data)
+      climatology <- terra::crop(climatology[[1]],data)
       data <- terra::crop(terra::mask(data, climatology), climatology)
       data <- terra::resample(data, climatology)
     }

R/make_2d_climatology_gridded.R

@@ -63,7 +63,8 @@ make_2d_climatology_gridded <- function(data.in, var.name, agg.time, statistic,
 
     # 3. OPTIMIZATION: Mask significantly fewer aggregated subset layers
     if (use.shp) {
-      data.subset <- terra::mask(data.subset, shp.vect)
+
+      data.subset <- EDABUtilities::crop_nc_2d(data.subset, shp.vect)[[1]]
     }
 
     data.time.agg.ls[[i]] <- data.subset

R/make_2d_deg_day_gridded.R

@@ -97,7 +97,7 @@ make_2d_deg_day_gridded <- function(data.in, var.name, metric, ref.value, type,
 
     # OPTIMIZATION: Apply spatial mask on the final single aggregated layer
     if (use.shp) {
-      data.out <- terra::mask(data.out, shp.vect)
+      data.out <- EDABUtilities::crop_nc_2d(data.out, shp.vect)[[1]]
     }
 
     if (write.out) {

R/make_2d_summary_gridded.R

@@ -72,7 +72,8 @@ make_2d_summary_gridded <- function(data.in, var.name, statistics, agg.time, fil
 
     # OPTIMIZATION: Pre-crop to bounding box once to reduce spatial memory footprint immediately
     if (use.shp) {
-      data <- terra::crop(data, shp.vect)
+      data <- EDABUtilities::crop_nc_2d(data, shp.file = shp.vect,area.names = area.names,var.name = var.name  )[[1]]
+      # data <- terra::crop(data, shp.vect)
     }
 
     data.stat.ls <- list()

R/make_2d_summary_ts.R

@@ -94,6 +94,7 @@ make_2d_summary_ts <- function(data.in, var.name, statistics, agg.time, file.tim
       stop('monthly files not yet implemented')
     }
 
+
     if(terra::crs(data) != terra::crs(shp.vect)){
       data = terra::project(data, terra::crs(shp.vect))
     }
@@ -121,7 +122,7 @@ make_2d_summary_ts <- function(data.in, var.name, statistics, agg.time, file.tim
       }
 
       # OPTIMIZATION: Crop immediately to shapefile bounding box before ANY iterations
-      data <- terra::crop(data, shp.vect)
+      data <- EDABUtilities::crop_nc_2d(data, shp.vect)[[1]]
 
       # OPTIMIZATION: Pull terra::tapp entirely out of the area loop. 
       # Execute once per statistics across the master clipped extent.

R/mask_nc_2d.R

@@ -62,7 +62,8 @@ mask_nc_2d <- function(data.in, var.name, min.value, max.value, write.out = FALS
     data <- if (is.character(item)) terra::rast(item) else item
 
     if (use.shp) {
-      data <- terra::crop(data, shp.vect)
+      data <- EDABUtilities::crop_nc_2d(data, shp.file = shp.vect,area.names = area.names,var.name = var.name  )[[1]]
+      # data <- terra::crop(data, shp.vect)
       data <- terra::mask(data, shp.vect) 
     }