Skip to content

Optimization for Large Files - #39

Open
jcaracappa1 wants to merge 3 commits into
devfrom
#36/convert_longitude
Open

Optimization for Large Files#39
jcaracappa1 wants to merge 3 commits into
devfrom
#36/convert_longitude

Conversation

@jcaracappa1

Copy link
Copy Markdown
Collaborator

Improved memory management for global files

  • Changed how crop_nc_2d works to crop prior to converting longitude
  • Added crop_nc_2d to summary functions prior to spatial aggregation instead of terra::crop

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remaining comments which cannot be posted as a review comment to avoid GitHub Rate Limit

air

[air] reported by reviewdog 🐶

e1 <- as.vector(terra::ext(data.orig))
e2 <- as.vector(terra::ext(shp_crop))


[air] reported by reviewdog 🐶

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


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶

data.crop <- terra::crop(data.orig, shp_crop)


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶

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.")


[air] reported by reviewdog 🐶

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)


[air] reported by reviewdog 🐶

out.ls[[i]] <- data.crop


[air] reported by reviewdog 🐶

if (write.out == FALSE) {
if (is.character(data.in)) {
names(out.ls) <- basename(data.in)


[air] reported by reviewdog 🐶

names(out.ls) <- paste0("layer_", seq_along(out.ls))


[air] reported by reviewdog 🐶

return(out.ls)


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶

climatology <- EDABUtilities::crop_nc_2d(climatology[[1]], shp.file = shp.vect,area.names = area.names,var.name = var.name )[[1]]


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶

if (!(all(terra::res(data) == terra::res(climatology)) && all(terra::ext(data) == terra::ext(climatology)))) {
climatology <- terra::crop(climatology[[1]],data)


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶

data <- EDABUtilities::crop_nc_2d(data, shp.file = shp.vect,area.names = area.names,var.name = var.name )[[1]]


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶

# OPTIMIZATION: Pull terra::tapp entirely out of the area loop.


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶

data <- EDABUtilities::crop_nc_2d(data, shp.file = shp.vect,area.names = area.names,var.name = var.name )[[1]]


[air] reported by reviewdog 🐶

data <- terra::mask(data, shp.vect)


[air] reported by reviewdog 🐶

Comment on lines 12 to 13
convert_2d_longitude_gridded <- function(data.in, write.out = FALSE, output.files = NA) {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
convert_2d_longitude_gridded <- function(data.in, write.out = FALSE, output.files = NA) {
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)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change

Comment on lines 18 to 19
data.ls = EDABUtilities:::import_data(data.in)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
data.ls = EDABUtilities:::import_data(data.in)
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) {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change


current_data <- data.ls[[i]]
is_file_input <- is.character(current_data)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change

Comment thread R/crop_nc_2d.R
Comment on lines +33 to +34
target_col <- col
break

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
target_col <- col
break
target_col <- col
break

Comment thread R/crop_nc_2d.R
Comment on lines 37 to +39

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.")

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
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.")
if (is.null(target_col)) {
stop(
"None of the attributes in shp.file contain all specified area.names."
)

Comment thread R/crop_nc_2d.R
Comment on lines 41 to +42

# 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 , ]

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
# 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 , ]
shp.vect <- shp.vect[shp.vect[[target_col]][, 1] %in% area.names, ]

Comment thread R/crop_nc_2d.R
Comment on lines 44 to +51

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]])

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
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]])
out.ls <- list()
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]])

Comment thread R/crop_nc_2d.R
data.orig <- terra::rast(data.ls[[i]])
} else {
data.orig <- data.ls[[i]]
data.orig <- data.ls[[i]]

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
data.orig <- data.ls[[i]]
data.orig <- data.ls[[i]]

@jcaracappa1
jcaracappa1 requested a review from atyrell3 August 20, 2026 14:56

@atyrell3 atyrell3 left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Two thoughts...
(1) coordinate reference system: I think we should have a check to make sure that the input nc is in the same CRS as shp.file. We could add a CRS input and project both files to the same CRS. To do it more efficiently: extract the CRS from the nc --> project shp.file to the same as the nc file --> crop nc --> project to the output CRS

(2) rotate vs shift: I think it's safer to use terra::rotate rather than terra::shift ; this is the behavior we want, right?:

Image

changed from shift to rotate

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remaining comments which cannot be posted as a review comment to avoid GitHub Rate Limit

air

[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶

# OPTIMIZATION: Pull terra::tapp entirely out of the area loop.


[air] reported by reviewdog 🐶


[air] reported by reviewdog 🐶

data <- EDABUtilities::crop_nc_2d(data, shp.file = shp.vect,area.names = area.names,var.name = var.name )[[1]]


[air] reported by reviewdog 🐶

data <- terra::mask(data, shp.vect)


[air] reported by reviewdog 🐶

Comment thread R/crop_nc_2d.R
Comment on lines +56 to +60
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)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
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)
data.crs <- terra::crs(data.orig)
shp.crs <- terra::crs(shp.vect)
if (!identical(data.crs, shp.crs)) {
shp.vect <- terra::project(shp.vect, data.crs)

Comment thread R/crop_nc_2d.R
if(!identical(data.crs,shp.crs)){
shp.vect = terra::project(shp.vect,data.crs)
}

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change

Comment thread R/crop_nc_2d.R
Comment on lines +66 to 68
e_rast <- as.vector(terra::ext(data.orig))
shp_crop <- shp.vect

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
e_rast <- as.vector(terra::ext(data.orig))
shp_crop <- shp.vect
e_rast <- as.vector(terra::ext(data.orig))
shp_crop <- shp.vect

Comment thread R/crop_nc_2d.R
# 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)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
shp_crop <- terra::rotate(shp_crop, long = 0,split = T, left = F)
shp_crop <- terra::rotate(shp_crop, long = 0, split = T, left = F)

Comment thread R/crop_nc_2d.R
# 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)) {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
} else if (e_rast["xmax"] <= 180.001 && any(terra::ext(shp_crop)[1:2] > 180)) {
} else if (
e_rast["xmax"] <= 180.001 && any(terra::ext(shp_crop)[1:2] > 180)
) {

@@ -97,7 +97,7 @@ make_2d_deg_day_gridded <- function(data.in, var.name, metric, ref.value, type,

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change

data.out <- terra::mask(data.out, shp.vect)
data.out <- EDABUtilities::crop_nc_2d(data.out, shp.vect)[[1]]
}

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change

@@ -72,7 +72,8 @@ make_2d_summary_gridded <- function(data.in, var.name, statistics, agg.time, fil

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change

# 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]]

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
data <- EDABUtilities::crop_nc_2d(data, shp.file = shp.vect,area.names = area.names,var.name = var.name )[[1]]
data <- EDABUtilities::crop_nc_2d(
data,
shp.file = shp.vect,
area.names = area.names,
var.name = var.name
)[[1]]

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)
}

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change

@jcaracappa1

Copy link
Copy Markdown
Collaborator Author

@atyrell3 Good suggestions. I made those changes. still is relatively faster.

Comment thread R/import_data.R
#' @param var.name character. Name of variable to be subseted from dataset
#'
#' @return list of spatRasters
#'

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
#'
#'

Comment thread R/import_data.R
Comment on lines +11 to 13
import_data = function(data.in, var.name = NULL){

if (inherits(data.in, "SpatRaster")|inherits(data.in, "SpatRasterDataset")) {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
import_data = function(data.in, var.name = NULL){
if (inherits(data.in, "SpatRaster")|inherits(data.in, "SpatRasterDataset")) {
import_data <- function(data.in, var.name = NULL) {
if (
inherits(data.in, "SpatRaster") | inherits(data.in, "SpatRasterDataset")
) {

Comment thread R/import_data.R
@@ -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.")

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
stop("data.in must be a file path, a vector of file paths, a single SpatRaster, or a list of SpatRasters.")
stop(
"data.in must be a file path, a vector of file paths, a single SpatRaster, or a list of SpatRasters."
)

Comment thread R/import_data.R
Comment on lines 26 to +41

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)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
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)
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)

Comment thread R/import_data.R
}
})
}

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change

Comment thread R/make_2d_summary_ts.R
Comment on lines 126 to 127

# OPTIMIZATION: Pull terra::tapp entirely out of the area loop.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
# OPTIMIZATION: Pull terra::tapp entirely out of the area loop.
# OPTIMIZATION: Pull terra::tapp entirely out of the area loop.

Comment thread 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

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change

Comment thread R/mask_nc_2d.R

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]]

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
data <- EDABUtilities::crop_nc_2d(data, shp.file = shp.vect,area.names = area.names,var.name = var.name )[[1]]
data <- EDABUtilities::crop_nc_2d(
data,
shp.file = shp.vect,
area.names = area.names,
var.name = var.name
)[[1]]

Comment thread R/mask_nc_2d.R
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)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change
data <- terra::mask(data, shp.vect)
data <- terra::mask(data, shp.vect)

Comment thread R/mask_nc_2d.R
# data <- terra::crop(data, shp.vect)
data <- terra::mask(data, shp.vect)
}

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[air] reported by reviewdog 🐶

Suggested change

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants