diff --git a/DESCRIPTION b/DESCRIPTION index c2ce3fe..d514981 100644 --- a/DESCRIPTION +++ b/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 diff --git a/R/convert_2d_longitude_gridded.R b/R/convert_2d_longitude_gridded.R index 6866049..26a7fd8 100644 --- a/R/convert_2d_longitude_gridded.R +++ b/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.") diff --git a/R/crop_nc_2d.R b/R/crop_nc_2d.R index 3c18684..bd19062 100644 --- a/R/crop_nc_2d.R +++ b/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) } diff --git a/R/import_data.R b/R/import_data.R index c63cddc..f1eb503 100644 --- a/R/import_data.R +++ b/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) } \ No newline at end of file diff --git a/R/make_2d_anomaly_gridded.R b/R/make_2d_anomaly_gridded.R index abdd4c6..bd40c9e 100644 --- a/R/make_2d_anomaly_gridded.R +++ b/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) } diff --git a/R/make_2d_climatology_gridded.R b/R/make_2d_climatology_gridded.R index 427f4a0..efb814b 100644 --- a/R/make_2d_climatology_gridded.R +++ b/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 diff --git a/R/make_2d_deg_day_gridded.R b/R/make_2d_deg_day_gridded.R index 9d2dd0c..52614de 100644 --- a/R/make_2d_deg_day_gridded.R +++ b/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) { diff --git a/R/make_2d_summary_gridded.R b/R/make_2d_summary_gridded.R index 2f31b17..880ba29 100644 --- a/R/make_2d_summary_gridded.R +++ b/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() diff --git a/R/make_2d_summary_ts.R b/R/make_2d_summary_ts.R index d5f1d1a..5b6ddb5 100644 --- a/R/make_2d_summary_ts.R +++ b/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. diff --git a/R/mask_nc_2d.R b/R/mask_nc_2d.R index 1430ec0..cb585ab 100644 --- a/R/mask_nc_2d.R +++ b/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) } diff --git a/man/import_data.Rd b/man/import_data.Rd index 3183354..17bf57e 100644 --- a/man/import_data.Rd +++ b/man/import_data.Rd @@ -4,10 +4,12 @@ \alias{import_data} \title{Imports data objects of various types and returns a spatRaster list for other functions} \usage{ -import_data(data.in) +import_data(data.in, var.name = NULL) } \arguments{ \item{data.in}{character vector, list, or SpatRaster. Single file path, vector of file paths, single SpatRaster, or list of SpatRasters representing the spatial data.} + +\item{var.name}{character. Name of variable to be subseted from dataset} } \value{ list of spatRasters