diff --git a/DESCRIPTION b/DESCRIPTION
index d0c5d69..48ec313 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -9,18 +9,17 @@ Authors@R: c(person("Jorge", "Garcia Molinos", email = "jorgegmolinos@arc.hokuda
Description: Functions to calculate the velocity of climate change (VoCC) and related metrics including
both gradient-based (Burrows et al. 2011, Burrows et al. 2014), and distance-based (Hamann et a. 2013,
Garcia Molinos et al. 2017) approaches.
-Depends: R (>= 3.5)
+Depends: R (>= 4.1)
Imports:
assertthat,
CircStats,
cowplot,
data.table,
- doParallel,
dplyr,
- foreach,
+ furrr,
geosphere,
magrittr,
- parallelly,
+ purrr,
RColorBrewer,
rlang,
sf,
@@ -28,15 +27,15 @@ Imports:
terra,
tibble,
tidyr,
- tidyselect,
xts
Suggests:
+ future,
ggplot2,
knitr,
mapplots,
ncdf4,
+ parallelly,
patchwork,
- purrr,
rmarkdown,
scales,
tidyterra,
@@ -46,6 +45,6 @@ BugReports: https://github.com/MathMarEcol/VoCC/issues
License: AGPL (>= 3)
Encoding: UTF-8
LazyData: true
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
VignetteBuilder: knitr,
rmarkdown
diff --git a/NAMESPACE b/NAMESPACE
index 2517e72..09d1214 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -14,7 +14,6 @@ export(trajClas)
export(trajLine)
export(voccTraj)
importFrom(data.table,":=")
-importFrom(foreach,"%dopar%")
importFrom(magrittr,"%>%")
importFrom(rlang,.data)
importFrom(stats,na.omit)
diff --git a/R/angulo.R b/R/angulo.R
index 39744e2..b88010e 100644
--- a/R/angulo.R
+++ b/R/angulo.R
@@ -1,16 +1,16 @@
#' Internal. Angle associated to the spatial gradient
#' @param dx \code{numeric} giving the longitudinal gradient component
#' @param dy \code{numeric} giving the latitudinal gradient component
-#' @author Jorge Garcia Molinos and David S. Schoeman
#' angulo()
angulo <- function(dx, dy) {
- d <- cbind(dx, dy)
- angline <- function(rw) {
- angle <- ifelse(rw[2] < 0, 180 + CircStats::deg(atan(rw[1] / rw[2])),
- ifelse(rw[1] < 0, 360 + CircStats::deg(atan(rw[1] / rw[2])), CircStats::deg(atan(rw[1] / rw[2])))
- )
- return(angle)
- }
- return(apply(d, 1, angline))
+ # OPTIMIZATION: Fully vectorized angle calculation - eliminates apply() loop
+ # Convert atan to degrees once for all values
+ atan_deg <- CircStats::deg(atan(dx / dy))
+
+ # Vectorized conditional logic - much faster than apply()
+ angle <- ifelse(dy < 0, 180 + atan_deg,
+ ifelse(dx < 0, 360 + atan_deg, atan_deg))
+
+ return(angle)
}
diff --git a/R/climPCA.R b/R/climPCA.R
index 7356993..bd45331 100644
--- a/R/climPCA.R
+++ b/R/climPCA.R
@@ -20,10 +20,10 @@
#' @importFrom data.table ":="
#'
#' @export
-#' @author Jorge Garcia Molinos
+#'
#' @examples
-#' \dontrun{
-#' JapTC <- VoCC_get_data("JapTC.tif")
+#'
+#' JapTC <- terra::rast(system.file("extdata", "JapTC.tif", package = "VoCC"))
#'
#' comp <- climPCA(JapTC[[c(1, 3, 5)]], JapTC[[c(2, 4, 6)]],
#' trans = NA, cen = TRUE, sc = TRUE, th = 0.85)
@@ -31,24 +31,38 @@
#' # Create a data frame with the necessary variables in the required order (see climAna? for details)
#' clim <- comp[[2]][, c(2, 4, 3, 5, 1)]
#' clim[, c("x", "y")] <- terra::xyFromCell(JapTC[[1]], clim$cid)
-#' }
#'
climPCA <- function(climp, climf, trans = function(x) log(x), cen = TRUE, sc = TRUE, th = 0.8) {
+
+ .SD <- NULL
+
+ # OPTIMIZATION: Pre-calculate cell count to avoid repeated calls
+ n_cells <- terra::ncell(climp[[1]])
+
# get a data table with the pooled values (current/future) of the clim variables
clim <- data.table::data.table(rbind(terra::values(climp), terra::values(climf)))
- clim <- stats::na.omit(clim[, c("cid", "p") := list(rep(1:terra::ncell(climp[[1]]), times = 2), rep(c("present", "future"), each = terra::ncell(climp[[1]])))])
+ clim <- stats::na.omit(clim[, c("cid", "p") := list(rep(1:n_cells, times = 2), rep(c("present", "future"), each = n_cells))])
+
+ # OPTIMIZATION: Store column indices to avoid repeated column selection
+ data_cols <- !names(clim) %in% c("cid", "p")
# apply transformation if required
if (!is.na(trans)) {
- clim <- trans(clim[, -c("cid", "p")])
+ clim_data <- clim[, .SD, .SDcols = data_cols]
+ clim_data <- trans(clim_data)
+ # Rebuild clim with transformed data
+ clim <- cbind(clim_data, clim[, c("cid", "p")])
}
# apply PCA
- clim.pca <- stats::prcomp(clim[, -c("cid", "p")], center = cen, scale. = sc)
+ clim.pca <- stats::prcomp(clim[, .SD, .SDcols = data_cols], center = cen, scale. = sc)
+
+ # OPTIMIZATION: Vectorized variance calculation
+ sdev_squared <- clim.pca$sdev^2
+ cumvar_prop <- cumsum(sdev_squared) / sum(sdev_squared)
+ a <- which(cumvar_prop >= th)[1]
- # extract numper of components explaining more than th accumulated variance
- a <- which((cumsum((clim.pca$sdev)^2) / sum(clim.pca$sdev^2)) >= th)[1]
- val.pca <- clim.pca$x[, 1:a]
+ val.pca <- clim.pca$x[, 1:a, drop = FALSE]
val <- data.frame(val.pca, cid = clim$cid, p = clim$p)
# put it back in wide form
diff --git a/R/climPlot.R b/R/climPlot.R
index 73f979f..10d1acf 100644
--- a/R/climPlot.R
+++ b/R/climPlot.R
@@ -15,10 +15,9 @@
#' @seealso{\code{\link{dVoCC}}, \code{\link{climPCA}}}
#'
#' @export
-#' @author Jorge Garcia Molinos and Naoki H. Kumagai
#' @examples
-#' \dontrun{
-#' JapTC <- VoCC_get_data("JapTC.tif")
+#'
+#' JapTC <- terra::rast(system.file("extdata", "JapTC.tif", package = "VoCC"))
#'
#' # Plot climate space for the two first variables(annual precipitation and maximum temperature)
#' xy <- stats::na.omit(data.frame(
@@ -37,18 +36,24 @@
#' plot = out, filename = file.path(getwd(), "example_plot.pdf"),
#' width = 17, height = 17, unit = "cm"
#' )
-#' }
+#'
climPlot <- function(xy, x.binSize, y.binSize, x.name = "V1", y.name = "V2") {
xp <- xy[, 1]
yp <- xy[, 3]
xf <- xy[, 2]
yf <- xy[, 4]
+ # OPTIMIZATION: Pre-calculate ranges to avoid repeated min/max calls
+ x_combined <- c(xp, xf)
+ y_combined <- c(yp, yf)
+ x_range <- range(x_combined)
+ y_range <- range(y_combined)
+
# bins per axis
- x.nbins <- floor((abs(range(xp, xf)[2] - range(xp, xf)[1])) / x.binSize)
- y.nbins <- floor((abs(range(yp, yf)[2] - range(yp, yf)[1])) / y.binSize)
- x.bin <- seq(floor(min(cbind(xp, xf))), ceiling(max(cbind(xp, xf))), length = x.nbins)
- y.bin <- seq(floor(min(cbind(yp, yf))), ceiling(max(cbind(yp, yf))), length = y.nbins)
+ x.nbins <- floor(abs(x_range[2] - x_range[1]) / x.binSize)
+ y.nbins <- floor(abs(y_range[2] - y_range[1]) / y.binSize)
+ x.bin <- seq(floor(x_range[1]), ceiling(x_range[2]), length = x.nbins)
+ y.bin <- seq(floor(y_range[1]), ceiling(y_range[2]), length = y.nbins)
# define palette
rf <- grDevices::colorRampPalette(rev(RColorBrewer::brewer.pal(11, "Spectral")))
@@ -75,18 +80,20 @@ climPlot <- function(xy, x.binSize, y.binSize, x.name = "V1", y.name = "V2") {
freq2Dp[freq2Dp > UL] <- UL
freq2Df[freq2Df > UL] <- UL
- # novel (in future but not present, 2), remnant (in both, 1), and dissapearing (in present but not future, 3) climates
- freq2D <- diag(nrow = x.nbins, ncol = y.nbins)
- freq2D[] <- NA
- for (i in 1:x.nbins) {
- for (j in 1:y.nbins) {
- freq2D[i, j] <- ifelse(is.na(freq2Dp[i, j]) & !is.na(freq2Df[i, j]), 1,
- ifelse(!is.na(freq2Dp[i, j]) & is.na(freq2Df[i, j]), 2,
- ifelse(is.na(freq2Dp[i, j]) & is.na(freq2Df[i, j]), NA, 0)
- )
- )
- }
- }
+ # OPTIMIZATION: Vectorized climate classification - eliminates nested loops
+ freq2D <- matrix(NA, nrow = x.nbins, ncol = y.nbins)
+
+ # Vectorized logical operations - much faster than nested loops
+ present_na <- is.na(freq2Dp)
+ future_na <- is.na(freq2Df)
+
+ # Novel climates: in future but not present
+ freq2D[present_na & !future_na] <- 1
+ # Disappearing climates: in present but not future
+ freq2D[!present_na & future_na] <- 2
+ # Remnant climates: in both present and future
+ freq2D[!present_na & !future_na] <- 0
+ # NA remains NA (neither present nor future)
# plot climate space
Freq2Dpf <- rbind(
diff --git a/R/dVoCC.R b/R/dVoCC.R
index d9e7c69..f4916b2 100644
--- a/R/dVoCC.R
+++ b/R/dVoCC.R
@@ -5,7 +5,6 @@
#' other (target) cells in the future by reference to a specified climatic threshold. The function allows for the
#' specification of search distances and incorporates both least-cost path and Great Circle (as-the-crow-flies) distances.
#'
-#' @importFrom foreach %dopar%
#'
#' @param clim \code{data.frame} with the value for the climatic parameters (columns) by cell (rows), arranged as follows (see examples below):
#' The first 2n columns must contain the present and future values for each of the n climatic variables (V1p, V1f, V2p, V2f,...).
@@ -39,10 +38,9 @@
#' @seealso{\code{\link{climPCA}}, \code{\link{climPlot}}}
#'
#' @export
-#' @author Jorge Garcia Molinos
#' @examples
-#' \dontrun{
-#' JapTC <- VoCC_get_data("JapTC.tif")
+#'
+#' JapTC <- terra::rast(system.file("extdata", "JapTC.tif", package = "VoCC"))
#'
#' # Create a data frame with the necessary variables in the required order
#' clim <- stats::na.omit(data.frame(terra::values(JapTC), cid = 1:terra::ncell(JapTC)))
@@ -58,6 +56,7 @@
#' r1[avocc1$focal] <- avocc1$vel
#' terra::plot(r1)
#'
+#' \dontrun{
#' # Cell-specific, distance-unrestricted climate analogue velocity based on least-cost path distances
#' # First, create the conductance matrix (all land cells considered to have conductance of 1)
#' r <- JapTC[[1]]
@@ -81,8 +80,6 @@
dVoCC <- function(clim, n, tdiff, method = "Single", climTol, geoTol,
distfun = "GreatCircle", trans = NA, lonlat = TRUE) {
- geoDis <- climDis <- ang <- vel <- target <- cid <- a <- NULL # Fix devtools check warnings
-
if (distfun == "Euclidean" && lonlat == TRUE) {
print("Error: Euclidean distances specified for unprojected coordinates")
stop()
@@ -99,144 +96,118 @@ dVoCC <- function(clim, n, tdiff, method = "Single", climTol, geoTol,
dat <- stats::na.omit(data.table::data.table(clim))
- # matrix with the future climatic values for all cells
- fut <- dat[, seq(2, (2 * n), by = 2), with = FALSE]
-
- # Determine optimal number of cores, ensuring we don't exceed data rows
- ncores <- parallelly::availableCores(constraints = "connections", omit = 2)
- ncores <- min(ncores, nrow(dat)) # Don't use more cores than data rows
- ncores <- max(ncores, 1) # Ensure at least 1 core
-
- # Only use parallel processing if we have multiple cores and sufficient data
- if (ncores > 1 && nrow(dat) > ncores) {
- cuts <- cut(seq_len(nrow(dat)), ncores, labels = FALSE)
- cl <- parallelly::makeClusterPSOCK(ncores, autoStop = TRUE)
-
- doParallel::registerDoParallel(cl)
-
- result <- foreach::foreach(a = seq_len(ncores),
- .combine = rbind,
- .packages = c("terra", "gdistance", "geosphere", "data.table"),
- .multicombine = TRUE) %dopar% {
-
- Dat <- dat[cuts == a, ]
-
- resu <- data.table::data.table(
- focal = Dat$cid,
- target = as.integer(NA),
- climDis = as.double(NA),
- geoDis = as.double(NA),
- ang = as.double(NA),
- vel = as.double(NA)
- )
-
- i <- 0
- while (i <= nrow(Dat)) {
- i <- i + 1
-
- # for each focal cell subset target cell analogues (within ClimTol)
- pres <- as.numeric(Dat[i, seq(1, (2 * n), by = 2), with = FALSE])
- dif <- data.table::data.table(sweep(fut, 2, pres, "-"))
-
- # Identify future analogue cells
- if (method == "Single") { # Ohlemuller et al 2006 / Hamann et al 2015
- upper <- colnames(dif)
- l <- lapply(upper, function(x) call("<", call("abs", as.name(x)), climTol[grep(x, colnames(dif))]))
- ii <- Reduce(function(c1, c2) substitute(.c1 & .c2, list(.c1 = c1, .c2 = c2)), l)
- anacid <- dat$cid[dif[eval(ii), which = TRUE]] # cids analogue cells
- }
-
- if (method == "Variable") { # Garcia Molinos et al. 2017
- climTol <- as.numeric(Dat[i, ((2 * n) + 1):(3 * n), with = FALSE]) # focal cell tolerance
- upper <- colnames(dif)
- l <- lapply(upper, function(x) call("<", call("abs", as.name(x)), climTol[grep(x, colnames(dif))]))
- ii <- Reduce(function(c1, c2) substitute(.c1 & .c2, list(.c1 = c1, .c2 = c2)), l)
- anacid <- dat$cid[dif[eval(ii), which = TRUE]] # cids analogue cells
- }
-
- # LOCATE CLOSEST ANALOGUE
- if (length(anacid) > 0) {
- # check which of those are within distance and get the analogue at minimum distance
- if (distfun == "Euclidean") {
- d <- stats::dist(cbind(Dat$x[i], Dat$y[i]), cbind(dat$x[dat$cid %in% anacid], dat$y[dat$cid %in% anacid]))
- } # in x/y units
- if (distfun == "GreatCircle") {
- d <- (geosphere::distHaversine(cbind(Dat$x[i], Dat$y[i]), cbind(dat$x[dat$cid %in% anacid], dat$y[dat$cid %in% anacid]))) / 1000
- } # in km
-
- # LeastCost distances not supported - error is thrown at function start
-
- an <- anacid[d < geoTol] # cids analogue cells within search radius
- dis <- d[d < geoTol] # distance to candidate analogues
- if (length(an) > 0) {
- resu[i, target := an[which.min(dis)]] # cid of geographically closest climate analogue
- if (method == "Single") {
- resu[i, climDis := mean(as.numeric(dif[which(anacid == resu[i, target]), ]))]
- } # mean clim difference for the closest analogue
- resu[i, geoDis := min(dis)]
- resu[i, ang := geosphere::bearing(Dat[i, c("x", "y")], dat[cid == resu[i, target], c("x", "y")])]
- resu[i, vel := resu$geoDis[i] / tdiff]
- }
- }
- }
- return(resu)
- }
- } else {
- # Sequential processing for small datasets or limited cores
- result <- data.table::data.table(
- focal = dat$cid,
- target = as.integer(NA),
- climDis = as.double(NA),
- geoDis = as.double(NA),
- ang = as.double(NA),
- vel = as.double(NA)
- )
-
- for (i in seq_len(nrow(dat))) {
- # for each focal cell subset target cell analogues (within ClimTol)
- pres <- as.numeric(dat[i, seq(1, (2 * n), by = 2), with = FALSE])
- dif <- data.table::data.table(sweep(fut, 2, pres, "-"))
-
- # Identify future analogue cells
- if (method == "Single") { # Ohlemuller et al 2006 / Hamann et al 2015
- upper <- colnames(dif)
- l <- lapply(upper, function(x) call("<", call("abs", as.name(x)), climTol[grep(x, colnames(dif))]))
- ii <- Reduce(function(c1, c2) substitute(.c1 & .c2, list(.c1 = c1, .c2 = c2)), l)
- anacid <- dat$cid[dif[eval(ii), which = TRUE]] # cids analogue cells
- }
- if (method == "Variable") { # Garcia Molinos et al. 2017
- climTol <- as.numeric(dat[i, ((2 * n) + 1):(3 * n), with = FALSE]) # focal cell tolerance
- upper <- colnames(dif)
- l <- lapply(upper, function(x) call("<", call("abs", as.name(x)), climTol[grep(x, colnames(dif))]))
- ii <- Reduce(function(c1, c2) substitute(.c1 & .c2, list(.c1 = c1, .c2 = c2)), l)
- anacid <- dat$cid[dif[eval(ii), which = TRUE]] # cids analogue cells
- }
+ # Check if column cid (cell ID exists)
+ if (!"cid" %in% names(dat)) {
+ dat <- dat %>%
+ dplyr::mutate(cid = dplyr::row_number())
+ }
+
+ # OPTIMIZATION: Improved chunking strategy for better load balancing
+ n_cores <- future::availableCores()
+
+ # Calculate optimal chunk size based on data size and cores
+ # Aim for 2-4 chunks per core for good load balancing
+ target_chunks <- n_cores * 3
+ optimal_chunk_size <- max(100, ceiling(nrow(dat) / target_chunks))
+
+ # Ensure we don't create too many tiny chunks
+ actual_chunks <- min(target_chunks, ceiling(nrow(dat) / optimal_chunk_size))
+ final_chunk_size <- ceiling(nrow(dat) / actual_chunks)
+
+ # Create balanced chunks
+ dat_chunks <- split(dat, ceiling(seq_len(nrow(dat)) / final_chunk_size))
+
+ result <- dat_chunks %>%
+ furrr::future_map(pdVoCC,
+ dat_full = dat, # Pass full dataset for analogue search
+ n = n,
+ tdiff = tdiff,
+ method = method,
+ climTol = climTol,
+ geoTol = geoTol,
+ distfun = distfun,
+ trans = trans,
+ lonlat = lonlat,
+ .options = furrr::furrr_options(seed = TRUE),
+ .progress = TRUE) %>%
+ purrr::list_rbind()
+
+ return(result)
+
+}
+
+
- # LOCATE CLOSEST ANALOGUE
- if (length(anacid) > 0) {
- # check which of those are within distance and get the analogue at minimum distance
- if (distfun == "Euclidean") {
- d <- stats::dist(cbind(dat$x[i], dat$y[i]), cbind(dat$x[dat$cid %in% anacid], dat$y[dat$cid %in% anacid]))
- } # in x/y units
- if (distfun == "GreatCircle") {
- d <- (geosphere::distHaversine(cbind(dat$x[i], dat$y[i]), cbind(dat$x[dat$cid %in% anacid], dat$y[dat$cid %in% anacid]))) / 1000
- } # in km
-
- an <- anacid[d < geoTol] # cids analogue cells within search radius
- dis <- d[d < geoTol] # distance to candidate analogues
- if (length(an) > 0) {
- result[i, target := an[which.min(dis)]] # cid of geographically closest climate analogue
- if (method == "Single") {
- result[i, climDis := mean(as.numeric(dif[which(anacid == result[i, target]), ]))]
- } # mean clim difference for the closest analogue
- result[i, geoDis := min(dis)]
- result[i, ang := geosphere::bearing(dat[i, c("x", "y")], dat[cid == result[i, target], c("x", "y")])]
- result[i, vel := result$geoDis[i] / tdiff]
- }
+#' Run parallel dVoCC
+#'
+#' @noRd
+pdVoCC <- function(dat, dat_full, n, tdiff, method = "Single", climTol, geoTol,
+ distfun = "GreatCircle", trans = NA, lonlat = TRUE) {
+
+ geoDis <- climDis <- ang <- vel <- target <- cid <- a <- NULL # Fix devtools check warnings
+
+ # Ensure both datasets are data.tables
+ dat <- data.table::as.data.table(dat) # Focal cells (chunk)
+ dat_full <- data.table::as.data.table(dat_full) # All cells (for analogue search)
+
+ result <- data.table::data.table(
+ focal = dat$cid,
+ target = as.integer(NA),
+ climDis = as.double(NA),
+ geoDis = as.double(NA),
+ ang = as.double(NA),
+ vel = as.double(NA)
+ )
+
+ # FIXED: Use full dataset for analogue search, not just the chunk
+ fut <- dat_full[, seq(2, (2 * n), by = 2), with = FALSE]
+
+ for (i in seq_len(nrow(dat))) {
+
+ # for each focal cell subset target cell analogues (within ClimTol)
+ pres <- as.numeric(dat[i, seq(1, (2 * n), by = 2), with = FALSE])
+ dif <- data.table::data.table(sweep(fut, 2, pres, "-"))
+
+ # Identify future analogue cells
+ if (method == "Single") { # Ohlemuller et al 2006 / Hamann et al 2015
+ upper <- colnames(dif)
+ l <- lapply(upper, function(x) call("<", call("abs", as.name(x)), climTol[grep(x, colnames(dif))]))
+ ii <- Reduce(function(c1, c2) substitute(.c1 & .c2, list(.c1 = c1, .c2 = c2)), l)
+ anacid <- dat_full$cid[dif[eval(ii), which = TRUE]] # FIXED: Use full dataset for analogue search
+ }
+
+ if (method == "Variable") { # Garcia Molinos et al. 2017
+ climTol <- as.numeric(dat[i, ((2 * n) + 1):(3 * n), with = FALSE]) # focal cell tolerance
+ upper <- colnames(dif)
+ l <- lapply(upper, function(x) call("<", call("abs", as.name(x)), climTol[grep(x, colnames(dif))]))
+ ii <- Reduce(function(c1, c2) substitute(.c1 & .c2, list(.c1 = c1, .c2 = c2)), l)
+ anacid <- dat_full$cid[dif[eval(ii), which = TRUE]] # FIXED: Use full dataset for analogue search
+ }
+
+ # LOCATE CLOSEST ANALOGUE
+ if (length(anacid) > 0) {
+ # check which of those are within distance and get the analogue at minimum distance
+ if (distfun == "Euclidean") {
+ d <- stats::dist(cbind(dat$x[i], dat$y[i]), cbind(dat_full$x[dat_full$cid %in% anacid], dat_full$y[dat_full$cid %in% anacid]))
+ } # in x/y units
+ if (distfun == "GreatCircle") {
+ d <- (geosphere::distHaversine(cbind(dat$x[i], dat$y[i]), cbind(dat_full$x[dat_full$cid %in% anacid], dat_full$y[dat_full$cid %in% anacid]))) / 1000
+ } # in km
+
+ an <- anacid[d < geoTol] # cids analogue cells within search radius
+ dis <- d[d < geoTol] # distance to candidate analogues
+ if (length(an) > 0) {
+ result[i, target := an[which.min(dis)]] # cid of geographically closest climate analogue
+ if (method == "Single") {
+ result[i, climDis := mean(as.numeric(dif[which(anacid == result[i, target]), ]))]
+ } # mean clim difference for the closest analogue
+ result[i, geoDis := min(dis)]
+ result[i, ang := geosphere::bearing(dat[i, c("x", "y")], dat_full[cid == result[i, target], c("x", "y")])]
+ result[i, vel := result$geoDis[i] / tdiff]
}
}
- }
+ } # End of the loop
return(result)
}
diff --git a/R/gVoCC.R b/R/gVoCC.R
index 8a9e21d..c9c2737 100644
--- a/R/gVoCC.R
+++ b/R/gVoCC.R
@@ -16,11 +16,11 @@
#' @seealso{\code{\link{tempTrend}}, \code{\link{spatGrad}}}
#'
#' @export
-#' @author Jorge Garcia Molinos
#'
#' @examples
-#' \dontrun{
-#' HSST <- VoCC_get_data("HSST.tif")
+#'
+#' HSST <- terra::rast(system.file("extdata", "HadiSST.tif", package = "VoCCdata"))
+#'
#' yrSST <- sumSeries(HSST,
#' p = "1960-01/2009-12", yr0 = "1955-01-01", l = terra::nlyr(HSST),
#' fun = function(x) colMeans(x, na.rm = TRUE), freqin = "months", freqout = "years"
@@ -37,11 +37,23 @@
gVoCC <- function(tempTrend, spatGrad) {
VoCC <- tempTrend[[1]] / spatGrad[[1]]
- # velocity angles have opposite direction to the spatial climatic gradient if warming and same direction (cold to warm) if cooling
- ind <- which(terra::values(VoCC) > 0)
+ # OPTIMIZATION: Extract values once and use vectorized operations
+ VoCC_values <- terra::values(VoCC)
+ spatGrad_ang_values <- terra::values(spatGrad[[2]])
+
+ # Vectorized angle calculation with proper NA handling
+ warming_cells <- VoCC_values > 0 & !is.na(VoCC_values)
+ VoCCang_values <- spatGrad_ang_values
+ VoCCang_values[warming_cells] <- spatGrad_ang_values[warming_cells] + 180
+
+ # Handle angle wrapping with NA protection
+ needs_wrapping <- !is.na(VoCCang_values) & VoCCang_values >= 360
+ VoCCang_values[needs_wrapping] <- VoCCang_values[needs_wrapping] - 360
+
+ # Create output raster efficiently
VoCCang <- spatGrad[[2]]
- VoCCang[ind] <- spatGrad[[2]][ind] + 180
- VoCCang[] <- ifelse(VoCCang[] >= 360, VoCCang[] - 360, VoCCang[])
+ terra::values(VoCCang) <- VoCCang_values
+
output <- c(VoCC, VoCCang)
names(output) <- c("voccMag", "voccAng")
return(output)
diff --git a/R/resTime.R b/R/resTime.R
index a8e327a..8c6385d 100644
--- a/R/resTime.R
+++ b/R/resTime.R
@@ -15,19 +15,21 @@
#' @seealso{\code{\link{gVoCC}}}
#'
#' @export
-#' @author Jorge Garcia Molinos
#' @examples
#'
#' # Load example Exclusive Economic Zone polygon
-#' \dontrun{
-#' EEZ <- VoCC_get_data("EEZ.gpkg")
-#' HSST <- VoCC_get_data("HSST.tif")
+#'
+#' EEZ <- terra::vect(system.file("extdata", "EEZ.gpkg", package = "VoCC"))
+#'
+#' HSST <- terra::rast(system.file("extdata", "HadiSST.tif", package = "VoCCdata")) %>%
+#' terra::crop(terra::ext(EEZ))
#'
#' yrSST <- sumSeries(HSST,
#' p = "1969-01/2009-12", yr0 = "1955-01-01", l = terra::nlyr(HSST),
#' fun = function(x) colMeans(x, na.rm = TRUE),
#' freqin = "months", freqout = "years"
#' )
+#'
#' tr <- tempTrend(yrSST, th = 10)
#' sg <- spatGrad(yrSST, th = 0.0001, projected = FALSE)
#' v <- gVoCC(tr, sg)
@@ -42,16 +44,16 @@
#' a2
#'
#' # Using a user defined polygon
-#' x_coord <- c(-28, -20, -20.3, -25.5)
-#' y_coord <- c(60, 61, 63, 62)
+#' x_coord <- c(-28, -20, -20.3, -25.5, -28)
+#' y_coord <- c(60, 61, 63, 62, 60)
#' coords <- matrix(c(x_coord, y_coord), ncol = 2)
#' poly_sf <- sf::st_sf(geometry = sf::st_sfc(sf::st_polygon(list(coords))))
#' a3 <- resTime(poly_sf, vel, areapg = NA)
#'
#' terra::plot(vel)
-#' plot(sf::st_geometry(EEZ), add = TRUE)
+#' terra::plot(EEZ, add = TRUE)
#' plot(sf::st_geometry(poly_sf), add = TRUE)
-#' }
+#'
resTime <- function(pg, vel, areapg = NA) {
resTim <- v <- d <- NULL # Fix devtools check warnings
@@ -72,12 +74,12 @@ resTime <- function(pg, vel, areapg = NA) {
# Extract velocity values using terra::vect for efficiency
extracted_values <- terra::extract(vel, pg_vect, fun = mean, na.rm = TRUE)
-
+
# Handle case where extraction returns NULL or empty results
if (is.null(extracted_values) || nrow(extracted_values) == 0) {
stop("No values could be extracted from the raster. Check that polygons overlap with the raster and have the same coordinate system.")
}
-
+
# Extract the mean values, handling potential column name variations
if ("mean" %in% names(extracted_values)) {
RT[, v := extracted_values$mean]
diff --git a/R/shiftTime.R b/R/shiftTime.R
index 5148bc5..e25d98a 100644
--- a/R/shiftTime.R
+++ b/R/shiftTime.R
@@ -19,40 +19,48 @@
#' shifting climate in marine and terrestrial ecosystems. Science, 334, 652-655.
#'
#' @export
-#' @author Jorge Garcia Molinos and Michael T. Burrows
#' @examples
-#' \dontrun{
-#' HSST <- VoCC_get_data("HSST.tif")
+#'
+#' HSST <- terra::rast(system.file("extdata", "HadiSST.tif", package = "VoCCdata"))
#' Apr <- shiftTime(HSST, yr1 = 1960, yr2 = 2009, yr0 = 1955, th = 10, m = 4)
#'
#' terra::plot(Apr)
-#' }
+#'
shiftTime <- function(r, yr1, yr2, yr0, th, m) {
+ # OPTIMIZATION: Pre-calculate constants and month indices
+ yr_offset <- (yr1 - yr0) * 12
+
# 1. Long term trends in monthly values (e.g. deg/year if temperature)
- m1 <- ((yr1 - yr0) * 12) + m
+ m1 <- yr_offset + m
m2 <- ((yr2 - yr0) * 12) + m
r1 <- r[[seq(m1, m2, by = 12)]]
trend <- tempTrend(r1, th)[[1]]
- # 2. seasonal rate of shift in temperature centered on each month (deg/month) = difference in
- # temperature between preceding and following months divided by 2 months (slope) preceding month
- b <- ifelse((m - 1) == 0, 12, (m - 1))
- m1 <- ((yr1 - yr0) * 12) + b
- m2 <- ((yr2 - yr0) * 12) + b
- x2 <- r[[seq(m1, m2, by = 12)]]
+ # OPTIMIZATION: Pre-calculate preceding and following months with bounds checking
+ prev_month <- if (m == 1) 12 else (m - 1)
+ next_month <- if (m == 12) 1 else (m + 1)
+
+ # 2. seasonal rate of shift in temperature centered on each month (deg/month)
+ # preceding month
+ m1_prev <- yr_offset + prev_month
+ m2_prev <- ((yr2 - yr0) * 12) + prev_month
+ x2 <- r[[seq(m1_prev, m2_prev, by = 12)]]
# following month
- b <- ifelse((m + 1) == 13, 1, (m + 1))
- m1 <- ((yr1 - yr0) * 12) + b
- m2 <- ((yr2 - yr0) * 12) + b
- x3 <- r[[seq(m1, m2, by = 12)]]
+ m1_next <- yr_offset + next_month
+ m2_next <- ((yr2 - yr0) * 12) + next_month
+ x3 <- r[[seq(m1_next, m2_next, by = 12)]]
+
+ # OPTIMIZATION: Pre-calculate conversion factor
+ days_per_decade_factor <- 3652.5 / 12 # 304.375
# slope
x4 <- terra::mean((x3 - x2) / 2, na.rm = TRUE)
- # 3. seasonal shifts (month/year) converted to days per decade by multiplying by 10 years, 365.25 days per year and dividing by 12 months
- sShift <- (trend / x4) * (3652.5 / 12)
+ # 3. seasonal shifts (month/year) converted to days per decade
+ sShift <- (trend / x4) * days_per_decade_factor
sShift[sShift == Inf | sShift == -Inf] <- NA
+
r2 <- c(trend, x4, sShift)
names(r2) <- c("mTrend", "seaRate", "seaShift")
return(r2)
diff --git a/R/spatGrad.R b/R/spatGrad.R
index 60b7e87..c88bcf4 100644
--- a/R/spatGrad.R
+++ b/R/spatGrad.R
@@ -16,7 +16,8 @@
#' (Grad in C per km for unprojected rasters and C per spatial unit for projected rasters),
#' and the associated angle (Ang in degrees).
#'
-#' @references \href{http://science.sciencemag.org/content/334/6056/652}{Burrows et al. 2011}. The pace of shifting climate in marine and terrestrial ecosystems. Science, 334, 652-655.
+#' @references \href{http://science.sciencemag.org/content/334/6056/652}{Burrows et al. 2011}.
+#' The pace of shifting climate in marine and terrestrial ecosystems. Science, 334, 652-655.
#'
#' @seealso{\code{\link{tempTrend}}, \code{\link{gVoCC}}}
#'
@@ -24,10 +25,9 @@
#' @importFrom stats na.omit
#'
#' @export
-#' @author Jorge Garcia Molinos, David S. Schoeman, and Michael T. Burrows
#' @examples
-#' \dontrun{
-#' HSST <- VoCC_get_data("HSST.tif")
+#'
+#' HSST <- terra::rast(system.file("extdata", "HadiSST.tif", package = "VoCCdata"))
#'
#' yrSST <- sumSeries(HSST,
#' p = "1969-01/2009-12", yr0 = "1955-01-01", l = terra::nlyr(HSST),
@@ -39,15 +39,16 @@
#' sg <- spatGrad(yrSST, th = 0.0001, projected = FALSE)
#'
#' terra::plot(sg)
-#' }
#'
spatGrad <- function(r, th = -Inf, projected = FALSE) {
# Fix devtools check warnings
"." <- NULL
- gradNS1 <- gradNS2 <- gradNS3 <- gradNS4 <- gradNS5 <- gradNS6 <- gradWE1 <- gradWE2 <- gradWE3 <- gradWE4 <- gradWE5 <- gradWE6 <- NULL
+ gradNS1 <- gradNS2 <- gradNS3 <- gradNS4 <- gradNS5 <- gradNS6 <- NULL
+ gradWE1 <- gradWE2 <- gradWE3 <- gradWE4 <- gradWE5 <- gradWE6 <- NULL
sy <- sx <- NSgrad <- WEgrad <- NULL
- clim <- climE <- climN <- climNE <- climNW <- climS <- climSE <- climSW <- climW <- climFocal <- NULL
+ clim <- climE <- climN <- climNE <- climNW <- NULL
+ climS <- climSE <- climSW <- climW <- climFocal <- NULL
to <- code <- i.to <- LAT <- angle <- Grad <- .SD <- NULL
if (terra::nlyr(r) > 1) {
@@ -57,12 +58,20 @@ spatGrad <- function(r, th = -Inf, projected = FALSE) {
# get resolution of the raster
re <- terra::res(r)
+ # OPTIMIZATION: Pre-extract raster values once to avoid repeated calls
+ spatRaster_values <- terra::values(r)
+ n_cells <- terra::ncell(r)
+
+ # OPTIMIZATION: Pre-calculate row and column indices to avoid repeated terra calls
+ row_indices <- terra::rowFromCell(r, 1:n_cells)
+ col_indices <- terra::colFromCell(r, 1:n_cells)
+
# Create a columns for focal and each of its 8 adjacent cells
y <- data.table::data.table(terra::adjacent(r, 1:terra::ncell(r), directions = 8, pairs = TRUE))
- y <- na.omit(y[, climFocal := terra::values(r)[from]][order(from, to)]) # Get value for focal cell, order the table by raster sequence and omit NAs (land cells)
- y[, clim := terra::values(r)[to]] # Insert values for adjacent cells
- y[, sy := terra::rowFromCell(r, from) - terra::rowFromCell(r, to)] # Column to identify rows in the raster (N = 1, mid = 0, S = -1)
- y[, sx := terra::colFromCell(r, to) - terra::colFromCell(r, from)] # Same for columns (E = 1, mid = 0, W = -1)
+ y <- na.omit(y[, climFocal := spatRaster_values[from]][order(from, to)]) # OPTIMIZED: Use pre-extracted values
+ y[, clim := spatRaster_values[to]] # OPTIMIZED: Use pre-extracted values
+ y[, sy := row_indices[from] - row_indices[to]] # Column to identify rows in the raster (N = 1, mid = 0, S = -1)
+ y[, sx := col_indices[to] - col_indices[from]] # Same for columns (E = 1, mid = 0, W = -1)
y[sx > 1, sx := -1] # Sort out the W-E wrap at the dateline, part I
y[sx < -1, sx := 1] # Sort out the W-E wrap at the dateline, part II
y[, code := paste0(sx, sy)] # Make a unique code for each of the eight neighbouring cells
@@ -72,8 +81,11 @@ spatGrad <- function(r, th = -Inf, projected = FALSE) {
to = c("climE", "climW", "climNW", "climSW", "climNE", "climSE", "climN", "climS")),
on = "code", code := i.to]
y <- data.table::dcast(y[, c("from", "code", "clim")], from ~ code, value.var = "clim")
- y[, climFocal := terra::values(r)[from]] # Put climFocal back in
- y[, LAT := terra::yFromCell(r, from)] # Add focal cell latitude
+ y[, climFocal := spatRaster_values[from]] # OPTIMIZED: Use pre-extracted values
+
+ # OPTIMIZATION: Pre-calculate latitude values to avoid repeated terra calls
+ lat_values <- terra::yFromCell(r, 1:n_cells)
+ y[, LAT := lat_values[from]] # Add focal cell latitude using pre-calculated values
# Calculate individual spatial temperature gradients: grads (degC per km)
# WE gradients difference in temperatures for each western and eastern pairs divided by the distance between the cells in each pair (corrected for latitudinal distortion if unprojected)
@@ -99,19 +111,49 @@ spatGrad <- function(r, th = -Inf, projected = FALSE) {
y[, gradNS5 := (climFocal - climS) / (d * re[2])]
y[, gradNS6 := (climE - climSE) / (d * re[2])]
- # Calulate NS and WE gradients.
+ # OPTIMIZATION: Vectorized weighted mean calculations instead of apply()
+ # Calculate NS and WE gradients using vectorized operations
# NOTE: for angles to work (at least using simple positive and negative values on Cartesian axes),
# S-N & W-E gradients need to be positive)
- y[, WEgrad := apply(.SD, 1, function(x) stats::weighted.mean(x, c(1, 2, 1, 1, 2, 1), na.rm = TRUE)), .SDcols = 12:17]
- y[, NSgrad := apply(.SD, 1, function(x) stats::weighted.mean(x, c(1, 2, 1, 1, 2, 1), na.rm = TRUE)), .SDcols = 18:23]
+
+ # Extract gradient matrices for vectorized operations
+ WE_gradients <- as.matrix(y[, .(gradWE1, gradWE2, gradWE3, gradWE4, gradWE5, gradWE6)])
+ NS_gradients <- as.matrix(y[, .(gradNS1, gradNS2, gradNS3, gradNS4, gradNS5, gradNS6)])
+
+ # Vectorized weighted mean calculation
+ weights <- c(1, 2, 1, 1, 2, 1)
+
+ # Calculate weighted means with proper NA handling (equivalent to na.rm = TRUE)
+ # Replace NA values with 0 for matrix multiplication, but track valid values
+ WE_valid <- !is.na(WE_gradients)
+ NS_valid <- !is.na(NS_gradients)
+
+ # Create matrices with NAs replaced by 0 for calculation
+ WE_clean <- WE_gradients
+ NS_clean <- NS_gradients
+ WE_clean[!WE_valid] <- 0
+ NS_clean[!NS_valid] <- 0
+
+ # Calculate weighted sums only for valid (non-NA) values
+ WE_weighted_sums <- rowSums(WE_clean * rep(weights, each = nrow(WE_clean)))
+ NS_weighted_sums <- rowSums(NS_clean * rep(weights, each = nrow(NS_clean)))
+
+ # Calculate sum of weights for valid values only
+ WE_weight_sums <- rowSums(WE_valid * rep(weights, each = nrow(WE_valid)))
+ NS_weight_sums <- rowSums(NS_valid * rep(weights, each = nrow(NS_valid)))
+
+ # Calculate weighted means (equivalent to stats::weighted.mean with na.rm = TRUE)
+ y[, WEgrad := ifelse(WE_weight_sums > 0, WE_weighted_sums / WE_weight_sums, NA_real_)]
+ y[, NSgrad := ifelse(NS_weight_sums > 0, NS_weighted_sums / NS_weight_sums, NA_real_)]
+
y[is.na(WEgrad) & !is.na(NSgrad), WEgrad := 0L] # Where NSgrad does not exist, but WEgrad does, make NSgrad 0
y[!is.na(WEgrad) & is.na(NSgrad), NSgrad := 0L] # same the other way around
# Calculate angles of gradients (degrees) - adjusted for quadrant (0 deg is North)
y[, angle := angulo(.SD$WEgrad, .SD$NSgrad), .SDcols = c("WEgrad", "NSgrad")]
- # Calculate the vector sum of gradients (C/km)
- y[, Grad := sqrt(apply(cbind((y$WEgrad^2), (y$NSgrad^2)), 1, sum, na.rm = TRUE))]
+ # OPTIMIZATION: Vectorized magnitude calculation instead of apply()
+ y[, Grad := sqrt(WEgrad^2 + NSgrad^2)]
# Merge the reduced file back into the main file to undo the initial na.omit
from <- data.table::data.table(1:terra::ncell(r)) # Make ordered from cells
@@ -121,6 +163,13 @@ spatGrad <- function(r, th = -Inf, projected = FALSE) {
rAng[y$from] <- y$angle
rGrad[y$from] <- y$Grad
rGrad[rGrad[] < th] <- th
+
+ # MEMORY LEAK FIX: Clean up large intermediate objects
+ rm(spatRaster_values, lat_values, row_indices, col_indices,
+ WE_gradients, NS_gradients, WE_clean, NS_clean,
+ WE_valid, NS_valid, WE_weighted_sums, NS_weighted_sums,
+ WE_weight_sums, NS_weight_sums, y, from)
+
output <- c(rGrad, rAng)
names(output) <- c("Grad", "Ang")
diff --git a/R/splitLine.R b/R/splitLine.R
index bd61aea..5671502 100644
--- a/R/splitLine.R
+++ b/R/splitLine.R
@@ -2,8 +2,7 @@
#' @param A \code{numeric} giving coordinates of first point
#' @param B \code{numeric} giving coordinates of second point
#' @param n \code{numeric} number of segments to divide the distance between points with
-#' @author Jorge Garcia Molinos
-
+#'
splitLine <- function(A, B, n) {
# Remove sign change for dateline, if needed
B[, 1] <- with(B, ifelse(abs(B[, 1] - A[, 1]) > 180 & B[, 1] < 0, B[, 1] + 360, B[, 1]))
diff --git a/R/sumSeries.R b/R/sumSeries.R
index 3de447a..2bbfa92 100644
--- a/R/sumSeries.R
+++ b/R/sumSeries.R
@@ -26,12 +26,11 @@
#' of Matrices (and to Vectors). R package version 0.53.1.
#'
#' @export
-#' @author Jorge Garcia Molinos
#' @examples
-#' \dontrun{
+#'
#' # Monthly mean SST (HadISST) data for Europe Jan-1950 to Dec-2010
#'
-#' HSST <- VoCC_get_data("HSST.tif")
+#' HSST <- terra::rast(system.file("extdata", "HadiSST.tif", package = "VoCCdata"))
#'
#' # Calculate mean annual monthly SST
#'
@@ -66,7 +65,6 @@
#' p = "1969-01/2009-12", yr0 = "1950-01-01", l = terra::nlyr(HSST),
#' fun = myf, freqin = "months", freqout = "other"
#' )
-#' }
#'
sumSeries <- function(r, p, yr0, l = terra::nlyr(r), fun = function(x) colMeans(x, na.rm = TRUE), freqin = "months", freqout = "years") {
# construct xts object
@@ -93,18 +91,27 @@ sumSeries <- function(r, p, yr0, l = terra::nlyr(r), fun = function(x) colMeans(
s <- fun(x)
}
- # create raster stack
- for (i in 1:nrow(s)) {
- r2 <- terra::rast(r[[1]])
- r2[] <- as.numeric(s[i, ])
- if (i == 1) {
- r1 <- r2
- } else {
- r1 <- c(r1, r2)
- }
+ # MEMORY LEAK FIX: More efficient raster creation using terra::setValues
+ raster_list <- vector("list", nrow(s))
+ template_raster <- terra::rast(r[[1]])
+
+ # OPTIMIZED: Vectorized raster creation
+ for (i in seq_len(nrow(s))) {
+ r2 <- terra::rast(template_raster) # Use template to avoid repeated rast() calls
+ terra::values(r2) <- as.numeric(s[i, ]) # More efficient than r2[] <-
+ raster_list[[i]] <- r2
}
+
+ # MEMORY LEAK FIX: Create stack once from list instead of incremental growth
+ r1 <- terra::rast(raster_list)
+
+ # Set names before cleaning up variables (x is needed for stats::start)
if (freqout != "other") {
names(r1) <- seq(stats::start(x), length = terra::nlyr(r1), by = freqout)
}
+
+ # Clean up temporary objects after using them
+ rm(raster_list, template_raster, m, ts1, x, s)
+
return(r1)
}
diff --git a/R/tempTrend.R b/R/tempTrend.R
index 638fb53..4613522 100644
--- a/R/tempTrend.R
+++ b/R/tempTrend.R
@@ -17,10 +17,9 @@
#' @seealso{\code{\link{spatGrad}}, \code{\link{gVoCC}}}
#'
#' @export
-#' @author Jorge Garcia Molinos and Christopher J. Brown
#' @examples
-#' \dontrun{
-#' HSST <- VoCC_get_data("HSST.tif")
+#'
+#' HSST <- terra::rast(system.file("extdata", "HadiSST.tif", package = "VoCCdata"))
#'
#' yrSST <- sumSeries(HSST,
#' p = "1969-01/2009-12", yr0 = "1955-01-01", l = terra::nlyr(HSST),
@@ -32,42 +31,60 @@
#' tr <- tempTrend(yrSST, th = 10)
#'
#' terra::plot(tr)
-#' }
+#'
tempTrend <- function(r, th) {
y <- terra::values(r)
ocean <- which(rowSums(is.na(y)) != ncol(y)) # remove land cells
y <- t(y[ocean, ])
- N <- apply(y, 2, function(x) sum(!is.na(x)))
+
+ # OPTIMIZED: Vectorized NA counting using colSums instead of apply
+ N <- colSums(!is.na(y))
ind <- which(N >= th)
- y <- y[, ind] # drop cells with less than th observations
- N <- apply(y, 2, function(x) sum(!is.na(x)))
- x <- matrix(nrow = terra::nlyr(r), ncol = ncol(y))
- x[] <- 1:terra::nlyr(r)
+ y <- y[, ind, drop = FALSE] # drop cells with less than th observations
+ N <- N[ind] # Use pre-calculated N values
+
+ # OPTIMIZED: More efficient x matrix creation using outer product
+ n_layers <- terra::nlyr(r)
+ n_cells <- ncol(y)
+
+ # Create time vector once and replicate efficiently
+ time_vec <- seq_len(n_layers)
+ x <- matrix(rep(time_vec, n_cells), nrow = n_layers, ncol = n_cells)
# put NA values into the x values so they correspond with y
- x1 <- y
- x1[!is.na(x1)] <- 1
- x <- x * x1
+ x[is.na(y)] <- NA
- # calculate the sum terms
- sx <- apply(x, 2, sum, na.rm = T)
- sy <- apply(y, 2, sum, na.rm = T)
- sxx <- apply(x, 2, function(x) sum(x^2, na.rm = T))
- syy <- apply(y, 2, function(x) sum(x^2, na.rm = T))
- xy <- x * y
- sxy <- apply(xy, 2, sum, na.rm = T)
+ # OPTIMIZED: Vectorized sum calculations using colSums instead of apply
+ sx <- colSums(x, na.rm = TRUE)
+ sy <- colSums(y, na.rm = TRUE)
+ sxx <- colSums(x^2, na.rm = TRUE)
+ syy <- colSums(y^2, na.rm = TRUE)
+ sxy <- colSums(x * y, na.rm = TRUE)
- # Estimate slope coefficients and associated standard errors and p-values
- slope <- (N * sxy - (sx * sy)) / (N * sxx - sx^2)
- sres <- (N * syy - sy^2 - slope^2 * (N * sxx - sx^2)) / (N * (N - 2))
- SE <- suppressWarnings(sqrt((N * sres) / (N * sxx - sx^2)))
+ # OPTIMIZED: Vectorized slope calculations with better numerical stability
+ denominator <- N * sxx - sx^2
+ # Add small epsilon to prevent division by zero
+ denominator[denominator == 0] <- .Machine$double.eps
+
+ slope <- (N * sxy - (sx * sy)) / denominator
+ sres <- (N * syy - sy^2 - slope^2 * denominator) / (N * (N - 2))
+ SE <- suppressWarnings(sqrt((N * sres) / denominator))
Test <- slope / SE
- p <- mapply(function(x, y) (2 * stats::pt(abs(x), df = y - 2, lower.tail = FALSE)), x = Test, y = N)
- slpTrends <- sigTrends <- seTrends <- terra::rast(r[[1]])
- slpTrends[ocean[ind]] <- slope
- seTrends[ocean[ind]] <- SE
- sigTrends[ocean[ind]] <- p
+ # OPTIMIZED: Vectorized p-value calculation with bounds checking
+ df_vals <- pmax(N - 2, 1) # Ensure df >= 1
+ p <- 2 * stats::pt(abs(Test), df = df_vals, lower.tail = FALSE)
+
+ # OPTIMIZED: Direct raster creation and assignment
+ template_raster <- terra::rast(r[[1]])
+ slpTrends <- seTrends <- sigTrends <- template_raster
+
+ # Use vectorized assignment
+ ocean_ind <- ocean[ind]
+ slpTrends[ocean_ind] <- slope
+ seTrends[ocean_ind] <- SE
+ sigTrends[ocean_ind] <- p
+
output <- c(slpTrends, seTrends, sigTrends)
names(output) <- c("slpTrends", "seTrends", "sigTrends")
diff --git a/R/trajClas.R b/R/trajClas.R
index f029d29..2a34229 100644
--- a/R/trajClas.R
+++ b/R/trajClas.R
@@ -42,10 +42,13 @@
#' @seealso{\code{\link{voccTraj}}}
#'
#' @export
-#' @author Jorge Garcia Molinos
#' @examples
+#'
#' \dontrun{
-#' HSST <- VoCC_get_data("HSST.tif")
+#' EEZ <- terra::vect(system.file("extdata", "EEZ.gpkg", package = "VoCC"))
+#'
+#' HSST <- terra::rast(system.file("extdata", "HadiSST.tif", package = "VoCCdata")) %>%
+#' terra::crop(terra::ext(EEZ) + 10)
#'
#' # input raster layers
#' yrSST <- sumSeries(HSST,
@@ -61,16 +64,13 @@
#' ang <- v[[2]]
#'
#' # Get the set of starting cells for the trajectories and calculate trajectories
-#' # at 1/4-deg resolution (16 trajectories per 1-deg cell)
-#' mnd <- terra::disagg(mn, 4)
-#' veld <- terra::disagg(vel, 4)
-#' angd <- terra::disagg(ang, 4)
#' lonlat <- stats::na.omit(data.frame(
-#' terra::xyFromCell(veld, 1:terra::ncell(veld)),
-#' terra::values(veld), terra::values(angd), terra::values(mnd)
+#' terra::xyFromCell(vel, 1:terra::ncell(vel)),
+#' terra::values(vel), terra::values(ang), terra::values(mn)
#' ))[, 1:2]
#'
-#' traj <- voccTraj(lonlat, vel, ang, mn, tyr = 50, correct = TRUE)
+#' traj <- voccTraj(lonlat, vel, ang, mn, tstep = 1/4, tyr = 20, seed = 23)
+#'
#'
#' # Generate the trajectory-based classification
#' clas <- trajClas(traj, vel, ang, mn,
@@ -78,35 +78,38 @@
#' Nend = 45, Nst = 15, NFT = 70, DateLine = FALSE
#' )
#'
-#' # Define first the colour palette for the full set of categories
-#' my_col <- c(
-#' "gainsboro", "darkseagreen1", "coral4", "firebrick2", "mediumblue", "darkorange1",
-#' "magenta1", "cadetblue1", "yellow1"
-#' )
-#' # Keep only the categories present in our raster
-#' my_col <- my_col[sort(unique(terra::values(clas[[7]])))]
-#'
#' # Classify raster / build attribute table
#' clasr <- terra::as.factor(clas[[7]])
#' rat_r <- data.frame(ID = sort(unique(terra::values(clas[[7]]))),
#' trajcat = c("N-M", "S-M", "IS", "BS", "Srce",
#' "RS", "Cor", "Div", "Con")[sort(unique(terra::values(clas[[7]])))])
-#' terra::cats(clasr) <- rat_r
-#' # Produce the plot using the rasterVis levelplot function
-#' rasterVis::levelplot(clasr,
-#' col.regions = my_col,
-#' xlab = NULL, ylab = NULL, scales = list(draw = FALSE)
-#' )
+#' levels(clasr) <- rat_r
+#' terra::plot(clasr)
#' }
+#'
trajClas <- function(traj, vel, ang, mn, trajSt, tyr, nmL, smL, Nend, Nst, NFT, DateLine = FALSE) {
ang1 <- ang2 <- ang3 <- ang4 <- d1 <- d2 <- d3 <- d4 <- NULL # Fix devtools check warnings
- isink <- .SD <- .N <- cid <- coastal <- val <- trajIDs <- NULL # Fix devtools check warnings
+ isink <- .SD <- .N <- cid <- coastal <- val <- ID <- NULL # Fix devtools check warnings
+
+ # MEMORY LEAK FIX: Create template raster once and reuse
+ template_raster <- terra::rast(ang)
# Force loading of ang values to avoid lazy loading warnings
invisible(terra::values(ang))
- TrajEnd <- TrajFT <- TrajSt <- IntS <- BounS <- TrajClas <- terra::rast(ang)
+ TrajEnd <- terra::rast(template_raster)
+ terra::values(TrajEnd) <- 0
+ TrajFT <- terra::rast(template_raster)
+ terra::values(TrajFT) <- 0
+ TrajSt <- terra::rast(template_raster)
+ terra::values(TrajSt) <- 0
+ IntS <- terra::rast(template_raster)
+ terra::values(IntS) <- 0
+ BounS <- terra::rast(template_raster)
+ terra::values(BounS) <- 0
+ TrajClas <- terra::rast(template_raster)
+ terra::values(TrajClas) <- 0
# add cell ID to the data frame
traj <- data.table::data.table(traj)
@@ -117,16 +120,14 @@ trajClas <- function(traj, vel, ang, mn, trajSt, tyr, nmL, smL, Nend, Nst, NFT,
terra::values(TrajSt)[!is.na(terra::values(ang))] <- trajSt
# B. Number of traj ending in each cell
- tr <- traj[, data.table::.SD[.N], by = trajIDs] # subset last point of each trajectory
+ tr <- traj[, data.table::.SD[.N], by = ID] # subset last point of each trajectory
enTraj <- tr[, .N, by = cid]
- terra::values(TrajEnd) <- 0
TrajEnd <- terra::mask(TrajEnd, vel)
terra::values(TrajEnd)[enTraj$cid] <- enTraj$N
# C. Number of traj flowing through each cell
- cxtrj <- unique(traj, by = c("trajIDs", "cid"))
+ cxtrj <- unique(traj, by = c("ID", "cid"))
TotTraj <- cxtrj[, .N, by = cid] # total number of trajectories per cell
- terra::values(TrajFT) <- 0
TrajFT <- terra::mask(TrajFT, vel)
terra::values(TrajFT)[TotTraj$cid] <- TotTraj$N
TrajFT <- TrajFT - TrajEnd - TrajSt # subtract traj starting and ending to get those actually transversing the cell
@@ -152,27 +153,39 @@ trajClas <- function(traj, vel, ang, mn, trajSt, tyr, nmL, smL, Nend, Nst, NFT,
res_x/2, res_y/2), # top-right
ncol = 2, byrow = TRUE)
+ # OPTIMIZED: Vectorized corner coordinate calculation
# Get the four cells for each coordinate point
- # Initialize as numeric matrix explicitly
- a <- matrix(as.numeric(NA), nrow = nrow(coords), ncol = 4)
-
- for(i in 1:nrow(coords)) {
- if(!is.na(coords[i,1]) && !is.na(coords[i,2])) {
- # Calculate corner coordinates
- corner_coords <- sweep(offsets, 2, coords[i,], "+")
- # Get cell numbers for each corner - handle one at a time
- for(k in 1:4) {
- cell_num <- terra::cellFromXY(ang, corner_coords[k, , drop = FALSE])
- a[i, k] <- as.numeric(cell_num[1])
- }
+ n_coords <- nrow(coords)
+ a <- matrix(as.numeric(NA), nrow = n_coords, ncol = 4)
+
+ # Only process valid coordinates
+ valid_idx <- !is.na(coords[,1]) & !is.na(coords[,2])
+ valid_coords <- coords[valid_idx, , drop = FALSE]
+
+ if (nrow(valid_coords) > 0) {
+ # Vectorized corner calculation for all valid coordinates at once
+ n_valid <- nrow(valid_coords)
+
+ # Create all corner coordinates in one operation
+ all_corners <- array(NA, dim = c(n_valid, 4, 2))
+ for(k in 1:4) {
+ # Fix dimension mismatch by ensuring proper matrix dimensions
+ offset_matrix <- matrix(rep(offsets[k,], each = n_valid), nrow = n_valid, ncol = 2)
+ all_corners[, k, ] <- offset_matrix + valid_coords
}
- }
- # Sort each row as in original code
- for(i in 1:nrow(a)) {
- a[i,] <- sort(a[i,])
+ # Get cell numbers for all corners at once
+ for(k in 1:4) {
+ corner_matrix <- all_corners[, k, , drop = FALSE]
+ dim(corner_matrix) <- c(n_valid, 2)
+ cell_nums <- terra::cellFromXY(ang, corner_matrix)
+ a[valid_idx, k] <- as.numeric(cell_nums)
+ }
}
+ # OPTIMIZED: Vectorized row sorting using apply
+ a <- t(apply(a, 1, sort, na.last = TRUE))
+
# If date line crossing, correct sequences on date line
if (DateLine == TRUE) {
a[seq(terra::ncol(ang), by = terra::ncol(ang), length = terra::nrow(ang)), ] <- t(apply(a[seq(terra::ncol(ang), by = terra::ncol(ang), length = terra::nrow(ang)), ], 1, function(x) {
@@ -180,21 +193,17 @@ trajClas <- function(traj, vel, ang, mn, trajSt, tyr, nmL, smL, Nend, Nst, NFT,
}))
}
+ # OPTIMIZED: Vectorized angle extraction
# Extract the angles for each group of 4 cells
- # Use direct indexing approach that works better with terra
b <- matrix(NA, nrow = nrow(a), ncol = 4)
# Get all angle values once to avoid repeated calls
ang_values <- terra::values(ang)
+ max_cell <- length(ang_values)
- for(i in 1:nrow(a)) {
- for(j in 1:4) {
- cell_val <- a[i,j] # Should now be numeric
- if(!is.na(cell_val) && cell_val > 0 && cell_val <= length(ang_values)) {
- b[i,j] <- ang_values[cell_val]
- }
- }
- }
+ # Vectorized angle extraction
+ valid_cells <- !is.na(a) & a > 0 & a <= max_cell
+ b[valid_cells] <- ang_values[a[valid_cells]]
# Ensure a and b have the same number of rows before combining
if(nrow(a) != nrow(b)) {
# Pad b to match a if needed
@@ -213,7 +222,6 @@ trajClas <- function(traj, vel, ang, mn, trajSt, tyr, nmL, smL, Nend, Nst, NFT,
# get the cids for the cells contained in the sinks
celdas <- ll[isink == 1, 3:6]
- terra::values(IntS) <- 0
IntS <- terra::mask(IntS, vel)
# Convert data.table columns to vectors and combine
@@ -231,17 +239,22 @@ trajClas <- function(traj, vel, ang, mn, trajSt, tyr, nmL, smL, Nend, Nst, NFT,
# detect coastal cells
coast <- suppressWarnings(terra::boundaries(vel, inner = TRUE)) # terra uses 'inner' parameter instead of 'type'
+ # OPTIMIZED: Pre-extract raster values once to avoid repeated terra::values() calls
+ vel_values <- terra::values(vel)
+ mn_values <- terra::values(mn)
+ coast_values <- terra::values(coast)
+
# make a list of vel values and SST values for each coastal cells and their marine neighbours
- cc <- stats::na.omit(data.table::data.table(cid = 1:terra::ncell(vel), coast = terra::values(coast)))
+ cc <- stats::na.omit(data.table::data.table(cid = 1:terra::ncell(vel), coast = coast_values))
ad <- terra::adjacent(vel, cc$cid, directions = 8, include = TRUE) # matrix with adjacent cells
# Sort the adjacency matrix to ensure consistent ordering (replaces sorted=TRUE from raster package)
ad <- ad[order(ad[, 1], ad[, 2]), ]
ad <- data.table::data.table(
coastal = ad[, 1],
adjacent = ad[, 2],
- cvel = terra::values(vel)[ad[, 1]],
- ctemp = terra::values(mn)[ad[, 1]],
- atemp = terra::values(mn)[ad[, 2]]
+ cvel = vel_values[ad[, 1]],
+ ctemp = mn_values[ad[, 1]],
+ atemp = mn_values[ad[, 2]]
)
# locate the sinks
@@ -249,51 +262,85 @@ trajClas <- function(traj, vel, ang, mn, trajSt, tyr, nmL, smL, Nend, Nst, NFT,
j <- ad[, ifelse(.SD$cvel > 0, all(.SD$ctemp <= .SD$atemp), all(.SD$ctemp >= .SD$atemp)), by = coastal]
data.table::setkey(j)
j <- unique(j)
- terra::values(BounS) <- 0
BounS <- terra::mask(BounS, vel)
boundary_cells <- unique(subset(j$coastal, j$V == TRUE))
if(length(boundary_cells) > 0) {
terra::values(BounS)[as.numeric(boundary_cells)] <- 1
}
+ # OPTIMIZED: Pre-extract ang values and minimize raster operations
+ ang_values <- terra::values(ang)
+
# Total number of trajectories per cell and proportions per cell
TrajTotal <- TrajSt + TrajFT + TrajEnd
- terra::values(TrajTotal)[is.na(terra::values(ang))] <- NA
+ terra::values(TrajTotal)[is.na(ang_values)] <- NA
PropTrajEnd <- (TrajEnd / TrajTotal) * 100
PropTrajFT <- (TrajFT / TrajTotal) * 100
PropTrajSt <- (TrajSt / TrajTotal) * 100
- # reclassify by traj length
- rclM <- matrix(c(0, (nmL / tyr), 1, (nmL / tyr), (smL / tyr), 2, (smL / tyr), Inf, 3), ncol = 3, byrow = TRUE)
- v <- terra::rast(vel)
- terra::values(v) <- abs(terra::values(vel))
- ClassMov <- terra::classify(v, rclM)
+ # OPTIMIZED: Direct classification using pre-extracted values - avoid creating intermediate rasters
+ abs_vel_values <- abs(vel_values)
- # Classify the cells
- terra::values(TrajClas) <- 0
- TrajClas <- terra::mask(TrajClas, vel)
+ # Vectorized classification instead of terra::classify
+ ClassMov_values <- ifelse(abs_vel_values <= (nmL / tyr), 1,
+ ifelse(abs_vel_values <= (smL / tyr), 2, 3))
+ ClassMov_values[is.na(vel_values)] <- NA
+
+ # OPTIMIZED: Pre-extract all raster values for classification
+ IntS_values <- terra::values(IntS)
+ BounS_values <- terra::values(BounS)
+ PropTrajEnd_values <- terra::values(PropTrajEnd)
+ PropTrajSt_values <- terra::values(PropTrajSt)
+ PropTrajFT_values <- terra::values(PropTrajFT)
+
+ # Classify the cells using vectorized operations
+ TrajClas_values <- rep(0, terra::ncell(TrajClas))
+ TrajClas_values[is.na(vel_values)] <- NA
# capture non-moving (1)
- terra::values(TrajClas)[terra::values(ClassMov) == 1] <- 1
+ TrajClas_values[ClassMov_values == 1] <- 1
# capture slow-moving (2)
- terra::values(TrajClas)[terra::values(ClassMov) == 2] <- 2
+ TrajClas_values[ClassMov_values == 2] <- 2
# capture internal (3) and (4) boundary sinks
- terra::values(TrajClas)[terra::values(IntS) == 1] <- 3
- terra::values(TrajClas)[terra::values(BounS) == 1] <- 4
-
- # Classify remaining cells into sources(5), rel sinks(6), corridors(7), divergence(8) and convergence(9)
- d <- stats::na.omit(data.table::data.table(cid = 1:terra::ncell(TrajClas), val = terra::values(TrajClas)))
- d <- d[val == 0, 1]
- d[, Nend := terra::values(PropTrajEnd)[d$cid]]
- d[, Nst := terra::values(PropTrajSt)[d$cid]]
- d[, NFT := terra::values(PropTrajFT)[d$cid]]
- d$clas <- ifelse(d$Nend == 0, 5, ifelse(d$Nend > Nend & d$Nst < Nst, 6, ifelse(d$NFT > NFT, 7, ifelse(d$Nend < d$Nst, 8, 9))))
- terra::values(TrajClas)[d$cid] <- d$clas
+ TrajClas_values[IntS_values == 1] <- 3
+ TrajClas_values[BounS_values == 1] <- 4
+
+ # OPTIMIZED: Vectorized classification for remaining cells
+ remaining_cells <- which(TrajClas_values == 0)
+ if (length(remaining_cells) > 0) {
+ Nend_vals <- PropTrajEnd_values[remaining_cells]
+ Nst_vals <- PropTrajSt_values[remaining_cells]
+ NFT_vals <- PropTrajFT_values[remaining_cells]
+
+ # Vectorized classification logic
+ clas_vals <- ifelse(Nend_vals == 0, 5,
+ ifelse(Nend_vals > Nend & Nst_vals < Nst, 6,
+ ifelse(NFT_vals > NFT, 7,
+ ifelse(Nend_vals < Nst_vals, 8, 9))))
+ TrajClas_values[remaining_cells] <- clas_vals
+ }
+
+ # Set final values
+ terra::values(TrajClas) <- TrajClas_values
+
+ # Create ClassMov raster from the values before cleanup
+ ClassMov <- terra::rast(template_raster)
+ terra::values(ClassMov) <- ClassMov_values
+
+ # MEMORY LEAK FIX: Clean up large intermediate objects before final assembly
+ rm(TrajClas_values, ClassMov_values, IntS_values, BounS_values,
+ PropTrajEnd_values, PropTrajSt_values, PropTrajFT_values,
+ vel_values, mn_values, coast_values, ang_values, template_raster,
+ abs_vel_values)
# return raster
s <- c(PropTrajEnd, PropTrajFT, PropTrajSt, ClassMov, IntS, BounS, TrajClas)
names(s) <- c("PropEnd", "PropFT", "PropSt", "ClassL", "IntS", "BounS", "TrajClas")
+
+ # Force garbage collection
+ gc()
+
return(s)
}
diff --git a/R/trajLine.R b/R/trajLine.R
index 21007b2..a36931f 100644
--- a/R/trajLine.R
+++ b/R/trajLine.R
@@ -7,7 +7,7 @@
#' @param projx \code{CRS} detailing the coordinate reference system of the input data
#' (default geographic CRS).
#'
-#' @return A \code{SpatialLinesDataFrame} with one line per trajectory as specified in x.
+#' @return A \code{sf} object with one line per trajectory as specified in x.
#' To avoid artifacts, trajectories crossing the date line need to be split into two segments.
#' Where the trajectory on one side of the date line is only composed of a single point,
#' the trajectory won't be displayed (no line object created). The function assumes
@@ -16,10 +16,11 @@
#' @seealso{\code{\link{voccTraj}}}
#'
#' @export
-#' @author Jorge Garcia Molinos
#' @examples
-#' \dontrun{
-#' HSST <- VoCC_get_data("HSST.tif")
+#'
+#' EEZ <- terra::vect(system.file("extdata", "EEZ.gpkg", package = "VoCC"))
+#' HSST <- terra::rast(system.file("extdata", "HadiSST.tif", package = "VoCCdata")) %>%
+#' terra::crop(terra::ext(EEZ) + 10)
#'
#' yrSST <- sumSeries(HSST,
#' p = "1969-01/2009-12", yr0 = "1955-01-01", l = terra::nlyr(HSST),
@@ -41,41 +42,54 @@
#' ))[, 1:2]
#'
#' # Calculate trajectories.
-#' traj <- voccTraj(lonlat, vel, ang, mn, tyr = 50, correct = TRUE)
+#' traj <- voccTraj(lonlat, vel, ang, mn, tyr = 50, tstep = 1/12)
#'
#' # create a spatial line data frame from traj
#' lns <- trajLine(x = traj)
#' terra::plot(mn)
-#' terra::plot(lns, add = TRUE)
+#' plot(lns %>% sf::st_geometry(), add = TRUE)
+#'
#'
-#' # Export as ESRI shape file
-#' terra::writeVector(lns, filename = "velTraj", filetype = "ESRI Shapefile")
-#' }
trajLine <- function(x, projx = "EPSG:4326") {
+ x %>%
+ dplyr::group_split(.data$ID) %>%
+ furrr::future_map(get_trajLine, proj_x = projx,
+ .options = furrr::furrr_options(seed = TRUE),
+ .progress = TRUE) %>%
+ purrr::list_rbind() %>%
+ sf::st_sf()
+
+}
+
+
+#' @noRd
+get_trajLine <- function(x, proj_x){
+
# Get distance between first and last points
- d <- dplyr::bind_rows(dplyr::slice(x, 1), dplyr::slice(x, dplyr::n())) %>%
- sf::st_as_sf(coords = c("llon", "llat")) %>%
- sf::st_set_crs(projx) %>%
- sf::st_distance(.) %>%
+ d <- dplyr::bind_rows(dplyr::slice(x, 1), dplyr::slice(x, dplyr::n())) |>
+ sf::st_as_sf(coords = c("lon", "lat")) |>
+ sf::st_set_crs(proj_x) |>
+ sf::st_distance() |>
max()
# Get number of steps; anything <240 means that the tracer died on land
steps <- max(x$Steps)
# Get distance along the tracer
- line_string <- x %>%
- sf::st_as_sf(coords = c("llon", "llat")) %>%
- sf::st_set_crs(projx) %>%
- dplyr::summarise(do_union = FALSE) %>%
- sf::st_cast(to = "LINESTRING") %>%
+ line_string <- x |>
+ sf::st_as_sf(coords = c("lon", "lat")) |>
+ sf::st_set_crs(proj_x) |>
+ dplyr::summarise(do_union = FALSE) |>
+ sf::st_cast(to = "LINESTRING")
+
+# Seperate this out to allow use of sf::st_length(line_string)
+ line_string |>
dplyr::mutate(steps = steps,
line_distance = d,
- line_length = sf::st_length(.),
- cellID = dplyr::first(x$cellIDs),
- lon = dplyr::first(x$llon),
- lat = dplyr::first(x$llat))
-
- return(line_string)
+ line_length = sf::st_length(line_string),
+ ID = dplyr::first(x$ID),
+ lon = dplyr::first(x$lon),
+ lat = dplyr::first(x$lat))
}
diff --git a/R/utils_voccTraj.R b/R/utils_voccTraj.R
new file mode 100644
index 0000000..793b6f2
--- /dev/null
+++ b/R/utils_voccTraj.R
@@ -0,0 +1,369 @@
+# Internal helper functions ----------------------------------------------------
+
+# Function for finding minimum distance
+mfind <- function(rw) {
+ X <- which.min(rw)
+ return(X)
+}
+
+# Function to extract corner coordinates per row
+mplace <- function(rw) {
+ X <- rw[(rw[1] + 1)]
+ Y <- rw[rw[1]]
+ return(c(X, Y))
+}
+
+# Simple circular functions
+deg2rad <- function(deg) {
+ return(deg * pi / 180)
+}
+rad2deg <- function(rad) {
+ return(180 * rad / pi)
+}
+
+
+#' Find distances from departure point to corners of ttcell
+#' OPTIMIZED: Use direct Haversine formula instead of expensive sf operations
+#'
+#' @noRd
+get_dist <- function(y1, x1, y2, x2) {
+ # Convert degrees to radians
+ lat1_rad <- y1 * pi / 180
+ lon1_rad <- x1 * pi / 180
+ lat2_rad <- y2 * pi / 180
+ lon2_rad <- x2 * pi / 180
+
+ # Haversine formula - much faster than sf::st_distance
+ dlat <- lat2_rad - lat1_rad
+ dlon <- lon2_rad - lon1_rad
+
+ a <- sin(dlat/2)^2 + cos(lat1_rad) * cos(lat2_rad) * sin(dlon/2)^2
+ c <- 2 * atan2(sqrt(a), sqrt(1-a))
+
+ # Earth radius in meters (same units as sf::st_distance)
+ R <- 6371000
+ distance <- R * c
+
+ return(distance)
+}
+
+
+#' Function for to find closest cooler/warmer cell within the union of two levels of 8-cell adjacency from the "departed" and "destination" cells
+#' Use when working with the coarse grid!
+#' JDE - Add bfr as an argument to voccTraj?
+#' A from-to cell pair and coordinates of departure point; bfr is a search radius in km for nearest cooler cell, if you end on land
+#'
+#' @noRd
+get_dest_cell_coarse <- function(rw, x_res, y_res, bfr, vel_values, mn_values, vel_raster, mn_raster) {
+
+ # Add error handling for invalid inputs
+ if (any(is.na(rw[1:4]))) {
+ return(NA)
+ }
+
+ # Find clumps of ocean; the rule is that you can't jump from one clump to another,
+ # because this would mean passing over unbroken land
+ pos_depart <- tryCatch({
+ data.frame(x = purrr::pluck(rw, 3),
+ y = purrr::pluck(rw, 4)) %>% # Departure coordinates
+ sf::st_as_sf(coords = c("x", "y"), crs = "EPSG:4326")
+ }, error = function(e) {
+ return(NULL)
+ })
+
+ if (is.null(pos_depart)) {
+ return(NA)
+ }
+
+ clumped <- tryCatch({
+ get_clumps(pos_depart, mn_raster, bfr, x_res, y_res)
+ }, error = function(e) {
+ return(NULL)
+ })
+
+ if (is.null(clumped)) {
+ return(NA)
+ }
+
+ # Which clump did I start in?
+ r1 <- purrr::pluck(rw, 1) # We use this cell a lot, so let's just make it an object
+ from_clump <- terra::extract(clumped, terra::xyFromCell(mn_raster, r1)) %>%
+ unlist()
+
+ # What are the coordinates of cells within the search area that fall in the clump I came from?
+ search_xy <- terra::xyFromCell(clumped, which(clumped[] == from_clump)) %>%
+ as.data.frame()
+
+ # Get the ssts in the cells to search
+ or <- terra::extract(mn_raster, search_xy, cells = TRUE, xy = TRUE) %>%
+ dplyr::rename(sst = 2, x = 4, y = 5)
+
+ # MEMORY LEAK FIX: Use pre-extracted values instead of raster indexing
+ if (vel_values[r1] > 0) {
+ o <- or %>%
+ dplyr::filter(.data$sst < mn_values[r1]) %>%
+ stats::na.omit()
+ if (nrow(o) == 0) {
+ dest_cell <- NA
+ } else {
+ # OPTIMIZED: Use vectorized Haversine instead of sf::st_distance
+ depart_coords <- sf::st_coordinates(pos_depart)
+ distances <- get_dist(depart_coords[1,2], depart_coords[1,1], o$y, o$x)
+ dest_cell <- o$cell[which.min(distances)]
+ }
+ } else { # Otherwise, find nearest warmer cell, if there is one
+ o <- or %>%
+ dplyr::filter(.data$sst > mn_values[r1]) %>%
+ stats::na.omit()
+ if (nrow(o) == 0) {
+ dest_cell <- NA
+ } else {
+ # OPTIMIZED: Use vectorized Haversine instead of sf::st_distance
+ depart_coords <- sf::st_coordinates(pos_depart)
+ distances <- get_dist(depart_coords[1,2], depart_coords[1,1], o$y, o$x)
+ dest_cell <- o$cell[which.min(distances)]
+ }
+ }
+ return(dest_cell)
+}
+
+
+
+
+
+#' Function to find closest cooler/warmer cell within the union of two levels of 8-cell adjacency from the "departed" and "destination" cells
+#'
+#' @noRd
+get_dest_cell_fine <- function(rw, x_res, y_res, bfr, vel_values, mn_values, vel_raster, mn_raster) {
+
+ # Add error handling for invalid inputs
+ if (any(is.na(rw[1:4]))) {
+ return(NA)
+ }
+
+ # Find clumps of ocean; the rule is that you can't jump from one clump to another, because this would mean passing over unbroken land
+ pos_depart <- tryCatch({
+ data.frame(
+ x = purrr::pluck(rw, 3),
+ y = purrr::pluck(rw, 4)
+ )
+ }, error = function(e) {
+ return(NULL)
+ })
+
+ if (is.null(pos_depart)) {
+ return(NA)
+ }
+
+ xy <- tryCatch({
+ terra::xyFromCell(mn_raster, as.vector(as.matrix(rw[1:2]))) %>% # Coordinates of cell centres for start cell and end cell
+ as.data.frame() %>%
+ sf::st_as_sf(coords = c("x", "y"), crs = "EPSG:4326")
+ }, error = function(e) {
+ return(NULL)
+ })
+
+ if (is.null(xy)) {
+ return(NA)
+ }
+
+ clumped <- tryCatch({
+ get_clumps(xy, mn_raster, bfr, x_res, y_res)
+ }, error = function(e) {
+ return(NULL)
+ })
+
+ if (is.null(clumped)) {
+ return(NA)
+ }
+
+ # Which clump did I start in?
+ r1 <- purrr::pluck(rw, 1) # We use this cell a lot, so let's just make it an object
+ from_clump <- purrr::pluck(terra::extract(clumped, terra::xyFromCell(mn_raster, r1)), 1)
+
+ # What are the coordinates of cells within the search area that fall in the clump I came from?
+ search_xy <- terra::xyFromCell(clumped, which(clumped[] == from_clump)) %>%
+ as.data.frame()
+ search_cells <- terra::cellFromXY(mn_raster, search_xy) # Which cells are these
+
+ # Get the ssts in the cells to search
+ or <- terra::extract(mn_raster, search_xy, cells = TRUE, xy = TRUE) %>%
+ dplyr::rename(sst = 2, x = 4, y = 5)
+
+ # MEMORY LEAK FIX: Use pre-extracted values instead of raster indexing
+ if (vel_values[r1] > 0) {
+ # Find all cells in the search area that meet the sst criterion
+ o <- or %>%
+ dplyr::filter(.data$sst < mn_values[r1]) %>%
+ stats::na.omit()
+ if (nrow(o) == 0) {
+ dest_cell <- NA # Set condition with which to ID stuck cells (i.e., if there are no suitable cells in the search area)
+ } else {
+ # OPTIMIZED: Use vectorized Haversine instead of sf::st_distance
+ distances <- get_dist(pos_depart$y, pos_depart$x, o$y, o$x)
+ potential_dest_cells <- o %>%
+ dplyr::mutate(distances = distances) %>%
+ dplyr::arrange(.data$distances)
+ n_cells <- 25 # Number of "close" cells to pick among *** Can change this, or set as an argument (just remove the hard-code here)
+ if (nrow(potential_dest_cells) >= n_cells) {
+ dest_cell <- potential_dest_cells %>%
+ dplyr::slice_head(n = n_cells) %>% # Get the closest cells
+ dplyr::slice_sample(n = 1) %>% # Remove .preserve = TRUE which may cause seed issues
+ dplyr::pull(.data$cell)
+ } else { # If fewer than n_cells, sample 10% of what you have
+ n_sample <- max(1, ceiling(nrow(potential_dest_cells) / 10))
+ dest_cell <- potential_dest_cells %>%
+ dplyr::slice_head(n = n_sample) %>%
+ dplyr::slice_sample(n = 1) %>% # Remove .preserve = TRUE which may cause seed issues
+ dplyr::pull(.data$cell)
+ }
+ }
+ } else { # Otherwise, find nearest cells that is warmer, if there are any
+ # Find all cells in the search area that meet the sst criterion
+ o <- or %>%
+ dplyr::filter(.data$sst > mn_values[r1]) %>%
+ stats::na.omit()
+ if (nrow(o) == 0) {
+ dest_cell <- NA # Set condition with which to ID stuck cells (i.e., if there are no suitable cells in the search area)
+ } else {
+ # OPTIMIZED: Use vectorized Haversine instead of sf::st_distance
+ distances <- get_dist(pos_depart$y, pos_depart$x, o$y, o$x)
+ potential_dest_cells <- o %>%
+ dplyr::mutate(distances = distances) %>%
+ dplyr::arrange(.data$distances)
+
+ n_cells <- 25 # Number of "close" cells to pick among *** Can change this, or set as an argument (just remove the hard-code here)
+ if (nrow(potential_dest_cells) >= n_cells) {
+ dest_cell <- potential_dest_cells %>%
+ dplyr::slice_head(n = n_cells) %>% # Get the closest cells
+ dplyr::slice_sample(n = 1) %>% # Remove .preserve = TRUE which may cause seed issues
+ dplyr::pull(.data$cell)
+ } else { # If fewer than n_cells, sample 10% of what you have
+ n_sample <- max(1, ceiling(nrow(potential_dest_cells) / 10))
+ dest_cell <- potential_dest_cells %>%
+ dplyr::slice_head(n = n_sample) %>%
+ dplyr::slice_sample(n = 1) %>% # Remove .preserve = TRUE which may cause seed issues
+ dplyr::pull(.data$cell)
+ }
+ }
+ }
+ return(dest_cell)
+}
+
+
+# JDE - I think we should rearrange arguments to x, y to be consistent
+
+# Find new destination, given velocity (km/yr), angle (º), time step (yr) and initial coordinates (ºlon, ºlat); max allowed jump is 1 deg
+# vell = vel[fcells] %>% pull(1); angg = ang[fcells] %>% pull(1); timestep = tstep; ll = llold
+#'
+#' @noRd
+destcoords <- function(vell, angg, timestep, ll, y_res, x_res) {
+ # OPTIMIZED: Pre-compute trigonometric values to avoid repeated calculations
+ angg_rad <- deg2rad(angg)
+ cos_angg <- cos(angg_rad)
+ sin_angg <- sin(angg_rad)
+
+ # OPTIMIZED: Vectorized calculations with better numerical stability
+ abs_vell <- abs(vell)
+ latshift <- (abs_vell * timestep * cos_angg) / 111.325 # Calculate shift in lat
+ latnew <- ll[, 2] + latshift # Find new lat...first approximation
+
+ # OPTIMIZED: Pre-compute cos(deg2rad(latnew)) for longitude calculation
+ cos_latnew <- cos(deg2rad(latnew))
+ # Add small epsilon to prevent division by zero
+ cos_latnew[cos_latnew == 0] <- .Machine$double.eps
+
+ lonshift <- (abs_vell * timestep * sin_angg) / (111.325 * cos_latnew) # Shift in lon
+
+ # OPTIMIZED: Vectorized clamping instead of conditional assignment
+ lonshift <- pmax(pmin(lonshift, x_res), -x_res)
+
+ # OPTIMIZED: Vectorized adjustment for extreme longitude shifts
+ needs_adjustment <- (abs(lonshift) == x_res)
+ if (any(needs_adjustment)) {
+ cos_ll_adj <- cos(deg2rad(ll[needs_adjustment, 2]))
+ tan_angg_adj <- tan(angg_rad[needs_adjustment])
+ # Prevent division by zero
+ tan_angg_adj[tan_angg_adj == 0] <- .Machine$double.eps
+ latshift[needs_adjustment] <- (x_res * 111.325 * cos_ll_adj / tan_angg_adj) / 111.325
+ }
+
+ latnew <- ll[, 2] + latshift # Find new lat by adding the adjusted lats
+
+ # OPTIMIZED: Vectorized pole clamping
+ latnew <- pmax(pmin(latnew, 90), -90)
+
+ # Adjust lon
+ lonnew <- ll[, 1] + lonshift # Find new lon...first approximation
+
+ # OPTIMIZED: More efficient dateline adjustment
+ lonnew <- lonnew - (360 * floor((lonnew + 180) / 360))
+
+ # OPTIMIZED: Direct data.frame creation instead of pipe
+ result <- data.frame(dlon = lonnew, dlat = latnew)
+ return(result)
+}
+
+
+#'
+#' @noRd
+get_clumps <- function(xy, mn, bfr, x_res, y_res){
+
+ # MEMORY LEAK FIX: Add error handling and explicit memory management
+ tryCatch({
+ # MEMORY LEAK FIX: Limit buffer size to prevent excessive memory usage
+ max_bfr <- min(bfr, 200) # Cap buffer at 200km to prevent memory issues
+ sp_buffer <- sf::st_buffer(xy, max_bfr * 1000) # Buffer around departure point, remembering that buffer is in metres
+
+ buffer_zone <- terra::extract(mn, sp_buffer, cells = TRUE, xy = TRUE) %>%
+ dplyr::select(-"ID") %>%
+ dplyr::distinct() %>%
+ dplyr::rename(sst = 1) %>% #*** rename "climatology", if needed
+ dplyr::select("x", "y", "sst", "cell") %>%
+ tidyr::drop_na("sst")
+
+ # Check if buffer_zone is empty
+ if (nrow(buffer_zone) == 0) {
+ return(NULL)
+ }
+
+ # MEMORY LEAK FIX: Adaptive sampling based on buffer size and resolution
+ max_cells <- min(10000, ceiling(pi * max_bfr^2 / (x_res * y_res * 111.325^2)))
+ if (nrow(buffer_zone) > max_cells) {
+ # Stratified sampling to maintain spatial distribution
+ buffer_zone <- buffer_zone[sample(nrow(buffer_zone), max_cells), ]
+ }
+
+ # OPTIMIZATION: Pre-calculate raster bounds to avoid repeated min/max calls
+ x_range <- range(buffer_zone$x)
+ y_range <- range(buffer_zone$y)
+
+ clumped_rast <- terra::rast(
+ xmin = x_range[1] - x_res/2,
+ xmax = x_range[2] + x_res/2,
+ ymin = y_range[1] - y_res/2,
+ ymax = y_range[2] + y_res/2,
+ resolution = c(x_res, y_res),
+ crs = "EPSG:4326"
+ )
+
+ # MEMORY LEAK FIX: Create SpatVector more efficiently and clean up immediately
+ buffer_vect <- terra::vect(buffer_zone[, c("x", "y", "sst")], geom = c("x", "y"), crs = "EPSG:4326")
+
+ clumped <- terra::rasterize(
+ x = buffer_vect,
+ y = clumped_rast, # The template to rasterize onto
+ field = "sst") %>% # The data
+ terra::patches(directions = 8, allowGaps = FALSE)
+
+ # MEMORY LEAK FIX: Explicitly clean up temporary objects
+ rm(buffer_vect, clumped_rast, sp_buffer, buffer_zone, x_range, y_range)
+
+ return(clumped)
+
+ }, error = function(e) {
+ warning("Error in get_clumps: ", e$message)
+ return(NULL)
+ })
+}
+
diff --git a/R/voccTraj.R b/R/voccTraj.R
index 7810e9b..1883529 100644
--- a/R/voccTraj.R
+++ b/R/voccTraj.R
@@ -11,24 +11,26 @@
#' @param vel \code{raster} with the magnitude of gradient-based climate velocity.
#' @param ang \code{raster} with velocity angles in degrees.
#' @param mn \code{raster} with the overall mean climatic value over the period of interest.
-#' @param vel_c
-#' @param ang_c
-#' @param mn_c
-#' @param fine_coast
-#' @param lk_up
-#' @param tstep
-#' @param tyr \code{integer} temporal length of the period of interest.
+#' @param x_res Numeric. Resolution of the grid in longitude direction (degrees or km).
+#' @param y_res Numeric. Resolution of the grid in latitude direction (degrees or km).
+#' @param tstep Numeric. Timestep for each trajectory iteration (usually decimal year).
+#' @param tyr Integer. Temporal length of the period of interest (years).
+#' @param bfr Numeric. Distance around a cell to look for a cooler/warmer cell.
+#' @param grid_resolution Character. "coarse" (default) or "fine". Controls how land crossings are handled and allows for higher resolution grids.
+#' @param seed Integer. Random seed for reproducibility. If NULL (default), no seed is set.
#'
-#' @return a \code{data.frame} containing the coordinates ("x", "y") of the constituent
-#' points and identification number ("trajIDs") for each trajectory.
+#'
+#' @return a \code{tibble} containing the coordinates ("lon", "lat") of the constituent
+#' points, time step ("Steps"), identification number ("ID") for each trajectory, and cell IDs for start and end cells.
#'
#' @references \href{https://www.nature.com/articles/nature12976}{Burrows et al. 2014}. Geographical limits to species-range shifts are suggested by climate velocity. Nature, 507, 492-495.
#'
#' @seealso{\code{\link{gVoCC}}, \code{\link{trajClas}}}
#' @export
-#' @author Jorge Garcia Molinos, David S. Schoeman and Michael T. Burrows
#' @examples
#' \dontrun{
+#' HSST <- terra::rast(system.file("extdata", "HadiSST.tif", package = "VoCCdata"))
+#'
#' yrSST <- sumSeries(HSST,
#' p = "1960-01/2009-12", yr0 = "1955-01-01", l = terra::nlyr(HSST),
#' fun = function(x) colMeans(x, na.rm = TRUE),
@@ -56,369 +58,424 @@
#' ))[, 1:2]
#'
#' # Calculate trajectories
-#' # The following throws an error due to the trajectories moving beyond the raster extent
-#' traj <- voccTraj(lonlat, vel, ang, mn, tyr = 50)
-#'
-#' # This accounts for the extent issue
-#' traj <- voccTraj(lonlat, vel, ang, mn, tyr = 50, correct = TRUE)
+#' traj <- voccTraj(lonlat, vel, ang, mn, tstep = 1/4, tyr = 50)
#' }
#'
-voccTraj <- function(lonlat,
- vel, ang, mn,
- # vel_c, ang_c, mn_c,
- #fine_coast, lk_up,
- tstep, tyr = 50) {
+voccTraj <- function(lonlat, # Starting points
+ vel, ang, mn, # Components of velocity: speed, angle and climatology
+ tstep, # Timestep (usually decimal year)
+ x_res = NULL, y_res = NULL, # Resolution of the grid at which velocity was computed
+ tyr = 20, # Number of years to run for
+ bfr = 75,
+ grid_resolution = "coarse", # Set to "fine" if you have disaggregated to original velocity field to a finer resolution
+ seed = NULL) { # Random seed for reproducibility
+
# Setup -------------------------------------------------------------------
- y_res <- x_res <- terra::res(vel)[1] # Set resolution of operations
+ # Validation checks using assertthat
+ assertthat::assert_that(
+ # lonlat must be a data.frame with at least 2 columns
+ is.data.frame(lonlat),
+ msg = "lonlat must be a data.frame"
+ )
+
+ assertthat::assert_that(
+ ncol(lonlat) >= 2,
+ msg = "lonlat must have at least 2 columns (longitude and latitude)"
+ )
+
+ assertthat::assert_that(
+ nrow(lonlat) > 0,
+ msg = "lonlat must contain at least one row"
+ )
+
+ assertthat::assert_that(
+ # vel, ang, and mn must be terra rasters
+ inherits(vel, "SpatRaster"),
+ msg = "vel must be a terra SpatRaster"
+ )
+
+ assertthat::assert_that(
+ inherits(ang, "SpatRaster"),
+ msg = "ang must be a terra SpatRaster"
+ )
+
+ assertthat::assert_that(
+ inherits(mn, "SpatRaster"),
+ msg = "mn must be a terra SpatRaster"
+ )
+
+ assertthat::assert_that(
+ # tstep must be positive numeric
+ is.numeric(tstep) && length(tstep) == 1 && tstep > 0,
+ msg = "tstep must be a positive numeric value"
+ )
+
+ assertthat::assert_that(
+ # x_res and y_res must be positive if provided
+ is.null(x_res) || (is.numeric(x_res) && length(x_res) == 1 && x_res > 0),
+ msg = "x_res must be a positive numeric value if provided"
+ )
+
+ assertthat::assert_that(
+ is.null(y_res) || (is.numeric(y_res) && length(y_res) == 1 && y_res > 0),
+ msg = "y_res must be a positive numeric value if provided"
+ )
+
+ assertthat::assert_that(
+ # When either x_res or y_res is given, the other must also be given
+ (is.null(x_res) && is.null(y_res)) || (!is.null(x_res) && !is.null(y_res)),
+ msg = "When either x_res or y_res is provided, both must be provided"
+ )
+
+ assertthat::assert_that(
+ # tyr must be positive numeric
+ is.numeric(tyr) && length(tyr) == 1 && tyr > 0,
+ msg = "tyr must be a positive numeric value"
+ )
+
+ assertthat::assert_that(
+ # bfr must be positive numeric
+ is.numeric(bfr) && length(bfr) == 1 && bfr > 0,
+ msg = "bfr must be a positive numeric value"
+ )
+
+ assertthat::assert_that(
+ # grid_resolution can only be "coarse" or "fine"
+ is.character(grid_resolution) && length(grid_resolution) == 1 &&
+ grid_resolution %in% c("coarse", "fine"),
+ msg = "grid_resolution must be either 'coarse' or 'fine'"
+ )
+
+ assertthat::assert_that(
+ # When both x_res and y_res are NULL, grid_resolution must be "coarse"
+ !(is.null(x_res) && is.null(y_res)) || grid_resolution == "coarse",
+ msg = "When both x_res and y_res are NULL, grid_resolution must be 'coarse'"
+ )
+
+ assertthat::assert_that(
+ # seed must be numeric if provided
+ is.null(seed) || (is.numeric(seed) && length(seed) == 1),
+ msg = "seed must be a numeric value if provided"
+ )
+
+ if (is.null(x_res) | is.null(y_res)){
+ vel_res <- terra::res(vel)
+
+ x_res <- vel_res[1]
+ y_res <- vel_res[2]
+
+ }
+
+
+ # Set seed for reproducibility if provided
+ if (!is.null(seed)) {
+ set.seed(seed)
+ }
+
+ # A base raster with original resolution and extent
+ r_base <- terra::rast(res = c(x_res, y_res)) %>%
+ terra::crop(vel)
- # Constrain max velocity to avoid stepping over grid squares
- max_vel = 111.325*x_res/tstep
- vel[vel > max_vel] <- max_vel
- vel[vel < -max_vel] <- -max_vel
- # vel_c[vel_c > max_vel] <- max_vel
- # vel_c[vel_c < -max_vel] <- -max_vel
+ # MEMORY LEAK FIX: Don't modify original raster, work with values directly
+ max_vel <- 111.325 * x_res / tstep
+ vel_values <- terra::values(vel)
+ vel_values[vel_values > max_vel] <- max_vel
+ vel_values[vel_values < -max_vel] <- -max_vel
+ # Sort out start points
lonlat <- lonlat %>%
- dplyr::select(tidyselect::all_of(c("x", "y"))) %>% # Collect just lon and lat (in case there's anything else there)
+ dplyr::select("x", "y") %>% # Collect just lon and lat (in case there's anything else there)
as.data.frame()
+ # Get initial descriptors
tcells <- terra::cellFromXY(vel, lonlat) # # Cell IDs of starting cells
n <- nrow(lonlat) # Get number of cells in your sequence
- sflags <- rep(NA, n) # Set a string of cells that change from NA to 1 where the trajectory sticks
-
- # Set up variables to catch results, allocating the right amount of memory
- llon <- numeric((n * tyr / tstep) + n) # Needs to have a set of starting lons, plus one more set for each time step
- llat <- numeric((n * tyr / tstep) + n) # Needs to have a set of starting lats, plus one more set for each time step
- # cellIDs <- rep(tcells, ((tyr / tstep) + 1)) # Needs to have a set of starting cells, plus one more set for each time step—just as a reference
- cellIDs <- rep(1:n, ((tyr / tstep) + 1)) # Needs to have a set of starting cells, plus one more set for each time step—just as a reference
- cellIDend <- numeric((n * tyr / tstep) + n) # Needs to have a set of ending cells, plus one more set for each time step
- # coast <- llon != 0 # Needs to have a set of starting coastal flags, plus one more set for each time step...set up as boolean
- flags <- numeric((n * tyr / tstep) + n) # Needs to have a set of starting flags, plus one more set for each time step
- Steps <- numeric((n * tyr / tstep) + n) # Needs to have a set of starting steps, plus one more set for each time step
-
- # Populate the first n slots with starting points
- llon[1:n] <- lonlat[,1]
- llat[1:n] <- lonlat[,2]
- cellIDend[1:n] <- tcells
- # coast[1:n] <- scoast
- flags[1:n] <- NA
- Steps[1:n] <- 0
-
- # Get coarse cellIDs that exist in fine raster
- # fn_lk_up <- lk_up[] %>%
- # na.omit() %>%
- # unique()
-
- # Helper functions --------------------------------------------------------
-
- # Trajectory helper functions
-
- # Function for to grep min distance
- mfind <- function(rw){
- X <- which.min(rw)
- return(X)
- }
- # Function for to extract corner coordinates per row
- mplace <- function(rw){
- X <- rw[(rw[1]+1)]
- Y <- rw[rw[1]]
- return(c(X, Y))
- }
+ # OPTIMIZATION: Use the already constrained vel_values and extract ang_values
+ ang_values <- terra::values(ang)
+ mn_values <- terra::values(mn) # Pre-extract mn values too
- # Simple circular functions
- deg2rad <- function(deg) return(deg*pi/180)
- rad2deg <- function(rad) return(180*rad/pi)
-
- # A function to find the nearest coastal cell
- get_nearest_coastal_cell <- function(x, y, tccell, lk_up, ...) { # JDE added lk_up to arguments
- t_block <- which(lk_up[] == tccell) # The cellIDs of lk_up
- cst_pts <- terra::xyFromCell(lk_up, t_block) %>%
- as.data.frame() %>%
- sf::st_as_sf(coords = c("x", "y"), crs = terra::crs(terra::rast()))
- # Here, we can't search for corners and add random "fuzz" to we collect 10 random points (if there are 10), and select the nearest of those, instead. If there are 10 or fewer, just pick the nearest
- if(nrow(cst_pts) > 10) {
- cst_pts <- cst_pts[sample(1:nrow(cst_pts), 10, replace = FALSE),]
- } # The unspoken "else" is that we just retain the cst_pts we have
- pt <- sf::st_as_sf(data.frame(x = x, y = y), coords = c("x", "y"), crs = terra::crs(terra::rast())) # The point departed from
- nearest <- sf::st_distance(cst_pts, pt) %>%
- which.min()
- out <- cst_pts[nearest,] %>%
- sf::st_coordinates() %>%
- as.data.frame() %>%
- dplyr::rename(x = 1, y = 2)
- return(out)
- }
+ # MEMORY OPTIMIZATION: Pre-allocate list structure (not content) to avoid dynamic growth
+ # This creates a list of NULL pointers - no contiguous memory required
+ max_steps <- ceiling(tyr / tstep) + 1 # Maximum possible steps + initial
- # Find new destination, given velocity (km/yr), angle (º), time step (yr) and initial coordinates (ºlon, ºlat); max allowed jump is 1 deg
- # vell = vel[fcells] %>% pull(1); angg = ang[fcells] %>% pull(1); timestep = tstep; ll = llold
- destcoords <- function(vell, angg, timestep, ll, y_res, x_res){
-
- latshift <- (abs(vell) * timestep * cos(deg2rad(angg))) / 111.325 # Calculate shift in lat
- latnew <- ll[,2] + latshift # Find new lat...first approximation
- lonshift <- (abs(vell) * timestep * sin(deg2rad(angg))) / (111.325 * cos(deg2rad(latnew))) # Shift in lon
- # Limit large longitudinal jumps at high latitudes
-
- # Because we constrain velocity to be no more than 12 * y_res, all problems will be with lonshift
- # Limit lonshift to at most 1 cell
- lonshift[lonshift > x_res] <- x_res
- lonshift[lonshift < -x_res] <- -x_res
-
- # Now on that basis, adjust latshift
- x_gt <- which(lonshift == x_res) # Indices for adjusted lon shifts
- latshift[x_gt] <- ((x_res*111.325 * cos(deg2rad(ll[x_gt,2])))/tan(deg2rad(angg[x_gt])))/111.325 # Using trig on distances
- x_lt <- which(lonshift == -x_res) # Indices for adjusted lon shifts
- latshift[x_lt] <- ((x_res*111.325 * cos(deg2rad(ll[x_lt,2])))/tan(deg2rad(angg[x_lt])))/111.325 # Using trig on distances
- latnew <- ll[,2] + latshift # Find new lat by adding the adjusted lats
- # Stop new lat from jumping the poles
- latnew[latnew > 90] <- 90
- latnew[latnew < -90] <- -90
- # Adjust lon
- lonnew <- ll[,1] + lonshift # Find new lon...first approximation
- # Adjust for dateline jumps
- lonnew <- lonnew - (360 * floor((lonnew + 180) / 360))
- return(data.frame(lonnew, latnew) %>%
- stats::setNames(c("dlon", "dlat")))
- }
+ # Pre-allocate list structure only - each element will be allocated independently
+ llon <- vector("list", max_steps)
+ llat <- vector("list", max_steps)
+ cellIDend <- vector("list", max_steps)
+ Steps <- vector("list", max_steps)
+ cellIDs <- vector("list", max_steps)
+ trajIDs <- vector("list", max_steps)
- # Function for to find closest cooler/warmer cell within the union of two levels of 8-cell adjacency from the "departed" and "destination" cells
- # TODO JDE - The arguments are not correct. Only 1, but it should include. `mn`
- get_dest_cell <- function(rw) { # A to-from cell pair and the buffer around a cell on land that can be searched for a suitable cell
+ # Populate the first slots with starting points
+ cellIDs[[1]] <- tcells
+ trajIDs[[1]] <- 1:n
+ llon[[1]] <- lonlat[, 1]
+ llat[[1]] <- lonlat[, 2]
+ cellIDend[[1]] <- tcells
+ Steps[[1]] <- rep(0, n)
+ # Set up objects that keep track of things
+ sflags <- rep(NA, n)
+ trj_id <- trajIDs[[1]]
- x_res <- terra::res(mn)[1]
- y_res <- terra::res(mn)[2]
- # Find clumps of ocean; the rule is that you can't jump from one clump to another, because this would mean passing over unbroken land
- xy <- terra::xyFromCell(mn, as.vector(as.matrix(rw))) %>%
- as.data.frame()
+ # Loop --------------------------------------------------------------------
+ # Loop through the trajectories
+ n_steps <- ceiling(tyr / tstep)
- bfr = ((y_res*111.325) + 1) * 1000 # Set buffer to be 1 grid-square width at the equator + 1 km
+ # Add safety check for reasonable number of steps
+ if (n_steps > 1000) {
+ warning("Large number of trajectory steps (", n_steps, "). Consider reducing tyr or increasing tstep.")
+ }
- xy <- xy %>%
- sf::st_as_sf(coords = c("x", "y"), crs = "EPSG:4326")
+ # Track actual steps used for efficient final processing
+ actual_steps_used <- 1 # Start with 1 (initial step)
- sp_buffer <- sf::st_buffer(sf::st_as_sf(xy), bfr) # Remembering that buffer is in metres
+ for (i in seq_len(n_steps)) {
- buffer_zone <- terra::extract(mn, sp_buffer, cells = TRUE, xy = TRUE, touches = TRUE) %>%
- dplyr::select(-"ID") %>%
- dplyr::distinct() %>%
- dplyr::rename(sst = mean) %>% # JDE Changed from climatology to mean
- dplyr::select(tidyselect::all_of(c("x", "y", "sst", "cell"))) %>%
- tidyr::drop_na(.data$sst)
+ # Safety check: if no trajectories remain, break early
+ if (nrow(lonlat) == 0) {
+ message("All trajectories terminated at step ", i-1)
+ break
+ }
- # JDE Changes. When all cells are in a straight line (only 1 unique x or y),
- # rast() will not work. So here we add a more robust approach to creating
- # a rast() regardless of the data
+ llold <- lonlat # Take a copy of lonlat
- # clumped <- buffer_zone %>%
- # dplyr::select(-cell) %>%
- # rast(crs = "EPSG:4326") %>%
- # terra::patches(directions = 8, allowGaps = FALSE)
+ # OPTIMIZATION: Get cell IDs first, then extract values by indexing (much faster)
+ fcells <- terra::cellFromXY(vel, llold) # Get the cells that the trajectories start in
+ # MEMORY LEAK FIX: Add bounds checking to prevent invalid indexing
+ valid_fcells <- !is.na(fcells) & fcells > 0 & fcells <= length(vel_values)
+ if (!any(valid_fcells)) {
+ message("All trajectories moved out of bounds at step ", i)
+ break
+ }
- clumped_rast <- terra::rast(
- xmin = min(buffer_zone$x) - x_res/2,
- xmax = max(buffer_zone$x) + x_res/2,
- ymin = min(buffer_zone$y) - y_res/2,
- ymax = max(buffer_zone$y) + y_res/2,
- resolution = c(x_res, y_res),
- crs = "EPSG:4326"
- )
+ # Extract values by direct indexing (much faster than terra::extract)
+ vell <- vel_values[fcells[valid_fcells]]
+ angg <- ang_values[fcells[valid_fcells]]
- clumped <- terra::rasterize(
- x = buffer_zone %>% # Needs to be SpatVector to add the sst
- dplyr::select(-"cell") %>%
- terra::vect(geom = c("x", "y"), crs = "EPSG:4326"),
- y = clumped_rast, # The template to rasterize onto
- field = "sst") %>% # The data
- terra::patches(directions = 8, allowGaps = FALSE)
+ # Update tracking variables for valid cells only
+ if (!all(valid_fcells)) {
+ llold <- llold[valid_fcells, , drop = FALSE]
+ fcells <- fcells[valid_fcells]
+ trj_id <- trj_id[valid_fcells]
+ sflags <- sflags[valid_fcells]
+ }
- # JDE END CHANGES
+ # Get new locations
+ lonlat <- destcoords(vell, angg, tstep, llold, y_res, x_res) # Extract lon and lat of landing point
+ tcells <- terra::cellFromXY(vel, lonlat) # Get the cells that the trajectories end in
+ sflags[which(is.na(tcells))] <- 1 # Sets the trajectory to "stuck" if it exits the velocity field (i.e., tcells == NA)
+
+ # Remove out-of-bounds cells
+ in_bounds <- which(is.na(sflags))
+ llold <- llold[in_bounds, ]
+ fcells <- fcells[in_bounds]
+ lonlat <- lonlat[in_bounds, ]
+ tcells <- tcells[in_bounds]
+ sflags <- sflags[in_bounds]
+ trj_id <- trj_id[in_bounds]
+
+ # MEMORY LEAK FIX: Use pre-extracted values instead of raster indexing
+ onland <- which(is.na(vel_values[tcells]))
+
+ if (length(onland) > 0) { # Only bother if there is at least one cell that returns a NA velocity = land
+ if (grid_resolution == "fine") { #*** Fine switch
+
+ # Collect the stuff we need here for cells that are onland
+ fpos <- llold[onland, ]
+ tpos <- lonlat[onland, ]
+ SFlags <- sflags[onland]
+ fcell <- fcells[onland]
+ tcell <- tcells[onland]
+
+ ft <- data.frame(fcell = fcell,
+ tcell = tcell,
+ fx = purrr::pluck(fpos, 1),
+ fy = purrr::pluck(fpos, 2),
+ code = paste(fcell, tcell, sep = " "),
+ ref = 1:length(fcell))
+
+ # MEMORY LEAK FIX: Pass values instead of full raster objects
+ ttcell <- apply(ft[, 1:4], 1, get_dest_cell_fine, x_res = x_res, y_res = y_res, bfr = bfr,
+ vel_values = vel_values, mn_values = mn_values, vel_raster = vel, mn_raster = mn)
+
+ # Filter "stuck" flags here
+ stuck <- which(is.na(ttcell)) # This is set in get_dest_cell(), where no cell in the "catchment" has a suitable sst to facilitate movement
+ unstuck <- which(!is.na(ttcell))
+ SFlags[stuck] <- 1 # Adjust flags
+
+ # Make data frame to catch data needed to find new positions
+ ttpos <- data.frame(x = rep(NA, length(onland)), y = rep(NA, length(onland)))
+ ttpos[stuck, ] <- fpos[stuck, ] # If they're stuck, pass on starting points
+ ttcell[stuck] <- fcell[stuck] # If they're stuck, pass on starting cells
+
+ ttpos[unstuck, ] <- terra::xyFromCell(mn, ttcell[unstuck])
+
+ # Collect results
+ lonlat[onland, ] <- ttpos
+ tcells[onland] <- ttcell
+ sflags[onland] <- SFlags
+ } else {
+ # Old onland loop here
+ # Collect the stuff we need here
+ fpos <- llold[onland, ]
+ tpos <- lonlat[onland, ]
+ SFlags <- sflags[onland]
+ fcell <- fcells[onland]
+ tcell <- tcells[onland]
+
+ # ft <- data.frame(fcell = fcell, tcell = tcell, code = paste(fcell, tcell, sep = " "), ref = 1:length(fcell))
+ # ttcell <- apply(ft[,1:2], 1, get_dest_cell)
+ ft <- data.frame(
+ fcell = fcell, tcell = tcell,
+ fx = fpos %>% purrr::pluck(1),
+ fy = fpos %>% purrr::pluck(2),
+ code = paste(fcell, tcell, sep = " "), ref = 1:length(fcell)
+ )
+ # MEMORY LEAK FIX: Pass values instead of full raster objects
+ ttcell <- apply(ft[, 1:4], 1, get_dest_cell_coarse, x_res = x_res, y_res = y_res, bfr = bfr,
+ vel_values = vel_values, mn_values = mn_values, vel_raster = vel, mn_raster = mn)
+
+ # Filter "stuck" flags here
+ stuck <- which(is.na(ttcell)) # This is set in get_dest_cell(), where no cell in the "catchment" has a suitable sst to facilitate movement
+ unstuck <- which(!is.na(ttcell))
+ SFlags[stuck] <- 1 # Adjust flags
+
+ # Make data frame to catch data needed to find new positions #***done up to here
+ ttpos <- data.frame(x = rep(NA, length(onland)), y = rep(NA, length(onland)))
+ ttpos[stuck, ] <- fpos[stuck, ] # If they're stuck, pass on starting points
+ ttcell[stuck] <- fcell[stuck] # If they're stuck, pass on starting cells
+ if (length(unstuck) > 0) {
+ tt_original_cells <- terra::cellFromXY(mn, fpos[unstuck, ]) # Departure cells in the resolution of the original velocity field
+ ttdat <- tibble::tibble(ttcell = ttcell[unstuck]) %>% # Destination cells (nearest cell with appropriate sst)
+ dplyr::mutate(
+ loncent = terra::xFromCell(mn, ttcell), # Coordinates of destination cell
+ latcent = terra::yFromCell(mn, ttcell)
+ ) %>% # Coordinates of destination cell
+ dplyr::mutate(e = NA, w = NA, n = NA, s = NA, dlon = NA, dlat = NA) # To facilitate finding corners of the cells
+ corner_block_size <- 0.25 * x_res # The "corner" is set to 0.25 of the grid square at the original resolution
+
+ # Send trajectory to the nearest corner of the appropriate cell, where corner is a quarter of the grid size. Position is "fuzzed" within this corner at random.
+ # Use reproducible random numbers for package checking
+ n_points <- nrow(ttdat) #TODO Check that this is still 1. Not sure why we need it.
+ ttdat$e <- ttdat$loncent + (0.5 * x_res) - (stats::runif(n_points, -1, 1) * corner_block_size) # The centre of the "corner" +- random fuzz...
+ ttdat$w <- ttdat$loncent - (0.5 * x_res) + (stats::runif(n_points, -1, 1) * corner_block_size) # The centre of the "corner" +- random fuzz...
+ ttdat$n <- ttdat$latcent + (0.5 * y_res) - (stats::runif(n_points, -1, 1) * corner_block_size) # The centre of the "corner" +- random fuzz...
+ ttdat$s <- ttdat$latcent - (0.5 * y_res) + (stats::runif(n_points, -1, 1) * corner_block_size) # The centre of the "corner" +- random fuzz...
+ coords <- with(ttdat, cbind(n, e, n, w, s, w, s, e)) # NE, NW, SW, SE corners' coordinates
+
+ corners <- data.frame(ne = get_dist(fpos[unstuck, 2], fpos[unstuck, 1], coords[, 1], coords[, 2])) %>%
+ dplyr::mutate(
+ nw = get_dist(fpos[unstuck, 2], fpos[unstuck, 1], coords[, 3], coords[, 4]),
+ sw = get_dist(fpos[unstuck, 2], fpos[unstuck, 1], coords[, 5], coords[, 6]),
+ se = get_dist(fpos[unstuck, 2], fpos[unstuck, 1], coords[, 7], coords[, 8])
+ )
+
+ # Select nearest corner
+ cornset <- apply(corners, 1, mfind) * 2 # Identify which corners for each onland cell. Have to mul by 2 to shift along correctly.
+ cornset <- cbind(cornset, coords) # Add in coordinates
+ ttdat[, 8:9] <- data.frame(t(apply(cornset, 1, mplace))) # Add in coordinates of correct corner point
+ ttpos[unstuck, ] <- terra::xyFromCell(mn, ttcell[unstuck])
+ ttcell[unstuck] <- terra::cellFromXY(mn, ttpos[unstuck, ])
+ }
- # Which clump did I start in?
- r1 <- rw[1] %>%
- as.vector() # We use this cell a lot, so let's just make it an object
+ # Collect results
+ lonlat[onland, ] <- ttpos
+ tcells[onland] <- ttcell
+ sflags[onland] <- SFlags
+ }
+ }
- from_clump <- terra::extract(clumped, terra::xyFromCell(mn, r1)) %>%
- unlist() %>%
- as.vector()
+ # Pass on only those cells that are not stuck
+ # if(sum(is.na(lonlat[,1])) > 0) {break}
+ cells_to_keep <- which(is.na(sflags))
+ lonlat <- lonlat[cells_to_keep, ]
+
+ # MEMORY OPTIMIZATION: Direct assignment to pre-allocated list slots
+ # Each element is stored independently - no contiguous memory required
+ step_index <- i + 1
+ llon[[step_index]] <- lonlat[, 1]
+ llat[[step_index]] <- lonlat[, 2]
+ cellIDs[[step_index]] <- fcells[cells_to_keep]
+ cellIDend[[step_index]] <- tcells[cells_to_keep]
+ Steps[[step_index]] <- rep(i, length(tcells[cells_to_keep]))
+ trajIDs[[step_index]] <- trj_id[cells_to_keep]
+
+ # Update tracking variables
+ sflags <- sflags[cells_to_keep]
+ trj_id <- trj_id[cells_to_keep]
+ actual_steps_used <- step_index
+
+ # Progress reporting - only in interactive sessions or when explicitly requested
+ if (interactive() && getOption("VoCC.verbose", FALSE)) {
+ message("Step ", i, "/", tyr/tstep, " (", round(100 * i / (tyr / tstep), 1), "%) - ",
+ length(onland), " cells on land - ", Sys.time())
+ }
- # What are the coordinates of cells within the search area that fall in the clump I came from?
- search_xy <- terra::xyFromCell(clumped, which(clumped[] == from_clump)) %>%
- as.data.frame() #****
- search_cells <- terra::cellFromXY(mn, search_xy) # Which cells are these
+ }
- # Check if any of these are NOT in in EITHER fine coast or mn, and eliminate, if needed
- to_keep <- c(which(!is.na(terra::extract(mn, search_xy, ID = FALSE) %>%
- dplyr::pull(1)))) %>% #, # Coords in the coarse grid
- # which(search_cells %in% fn_lk_up)) %>%
- unique()
+ # MEMORY LEAK FIX: Clean up large objects before final processing
+ rm(vel_values, ang_values, mn_values)
+
+ # MEMORY OPTIMIZATION: Only process used slots and clean up progressively
+ # Trim to actual used length - unused slots remain as NULL (minimal memory)
+ if (actual_steps_used < max_steps) {
+ llon <- llon[1:actual_steps_used]
+ llat <- llat[1:actual_steps_used]
+ Steps <- Steps[1:actual_steps_used]
+ trajIDs <- trajIDs[1:actual_steps_used]
+ cellIDs <- cellIDs[1:actual_steps_used]
+ cellIDend <- cellIDend[1:actual_steps_used]
+ }
- search_xy1 <- search_xy[to_keep,] %>%
- as.data.frame() %>%
- dplyr::distinct()
+ # MEMORY LEAK FIX: Progressive cleanup to minimize peak memory usage
+ # Each unlist operation works on independent memory chunks
+ steps_vec <- unlist(Steps, use.names = FALSE)
+ rm(Steps)
- # Get the ssts in the cells to search
- or <- terra::extract(mn, search_xy1, cells = TRUE, xy = TRUE) %>%
- dplyr::rename(sst = 2, x = 4, y = 5)
+ lon_vec <- unlist(llon, use.names = FALSE)
+ rm(llon)
+ lat_vec <- unlist(llat, use.names = FALSE)
+ rm(llat)
- # If velocity is positive, find nearest cell that is cooler, if there is one
- if(unlist(vel[r1]) > 0) {
- o <- or %>%
- dplyr::filter(.data$sst < unlist(mn[r1])) %>%
- na.omit()
- if(nrow(o) == 0) {
- dest_cell <- NA # Set condition with which to ID stuck cells
- } else {
- pt <- terra::xyFromCell(mn, r1) %>%
- data.frame()
+ id_vec <- unlist(trajIDs, use.names = FALSE)
+ rm(trajIDs)
- dest_cell <- sf::st_distance(sf::st_as_sf(pt, coords = c("x", "y"), crs = terra::crs(terra::rast())),
- sf::st_as_sf(data.frame(o), coords = c("x", "y"), crs = terra::crs(terra::rast()))) %>%
- which.min() %>%
- o$cell[.]
- }
- } else { # Otherwise, find nearest cell that is warmer, if there is one
- o <- or %>%
- dplyr::filter(.data$sst > unlist(mn[r1])) %>%
- na.omit()
- if(nrow(o) == 0) {
- dest_cell <- NA # Set condition with which to ID stuck cells
- } else {
- pt <- terra::xyFromCell(mn, r1) %>%
- data.frame()
+ start_cell_vec <- unlist(cellIDs, use.names = FALSE)
+ rm(cellIDs)
- dest_cell <- sf::st_distance(sf::st_as_sf(pt, coords = c("x", "y"), crs = terra::crs(terra::rast())),
- sf::st_as_sf(o, coords = c("x", "y"), crs = terra::crs(terra::rast()))) %>%
- which.min() %>%
- o$cell[.]
- }
- }
- return(dest_cell)
- }
+ end_cell_vec <- unlist(cellIDend, use.names = FALSE)
+ rm(cellIDend)
- # Loop --------------------------------------------------------------------
+ # Create final result with cleaned vectors
+ result <- tibble::tibble(
+ Steps = steps_vec,
+ lon = lon_vec,
+ lat = lat_vec,
+ ID = id_vec,
+ start_cell = start_cell_vec,
+ end_cell = end_cell_vec
+ )
- # Loop through the trajectories
- for(i in 1:(tyr/tstep)) { # 1:(tyr/tstep)
- llold <- lonlat # Take a copy of lonlat
- fcells <- terra::cellFromXY(vel,llold) # Get the cells that the trajectories start in
+ # Clean up final vectors
+ rm(steps_vec, lon_vec, lat_vec, id_vec, start_cell_vec, end_cell_vec)
- # Pull velocity and angle from the "inland" versions, because we always ensure that we can find destination cells for departures within the coastal "blind spot"
- # TIN: changed this from vel_c to vel and ang_c to ang because we don't have coarse/fine versions
- vc <- vel[fcells] %>% dplyr::pull(1)
- ac <- ang[fcells] %>% dplyr::pull(1)
- llc <- llold
+ # Clean up any remaining large objects from the loop
+ if (exists("lonlat")) rm(lonlat)
+ if (exists("llold")) rm(llold)
+ if (exists("fcells")) rm(fcells)
+ if (exists("tcells")) rm(tcells)
+ if (exists("sflags")) rm(sflags)
+ if (exists("trj_id")) rm(trj_id)
- # Get new locations
- lonlat <- destcoords(vc, ac, tstep, llc, y_res, x_res) # Extract lon and lat of landing point
- tcells <- terra::cellFromXY(vel,lonlat) # Get the cells that the trajectories end in
- sflags[which(is.na(tcells))] <- 1 # Sets the trajectory to "stuck" if it exits the velocity field (i.e., tcells == NA)
-
- # Bounce stuck cells
- stuck_cells <- which(!is.na(sflags))
- lonlat[stuck_cells,] <- llold[stuck_cells,] # Bounce the corresponding coordinates across for "stuck" cells
- tcells[stuck_cells] <- fcells[stuck_cells] # Bounce the corresponding cells across for "stuck" cells
-
- # Deal with cells that end on land
- onland <- which(is.na(vel[tcells])) %>% # Identify which rows of velend are on land, provided that they are not stuck
- dplyr::setdiff(stuck_cells) # Ignoring stuck cells
-
- # if(length(onland) > 0){ # Only bother if there is at least one cell that returns a NA velocity = land
- # Here, you need to check whether it really *IS* on land, or whether it is just in the coastal "blind spot"
- # fn <- terra::extract(fine_coast, lonlat[onland,], ID = FALSE) %>%
- # pull(1)
- # onland <- onland[which(is.na(fn))]
- if(length(onland) > 0) {
- # Collect the stuff we need here for cells that are really onland
- fpos <- llold[onland,]
- tpos <- lonlat[onland,]
- SFlags <- sflags[onland]
- fcell <- fcells[onland]
- tcell <- tcells[onland]
- ft <- data.frame(fcell = fcell, tcell = tcell, code = paste(fcell, tcell, sep = " "), ref = 1:length(fcell))
-
- # Get the destination cells of those that are hitting the coast
- ttcell <- apply(ft[,1:2], 1, get_dest_cell) # TODO using apply in places and map in others. Can we consistently move to purrr?
-
- # Filter "stuck" flags here
- stuck <- which(is.na(ttcell))
- unstuck <- which(!is.na(ttcell))
- SFlags[stuck] <- 1 # Adjust flags
-
- # Make data frame to catch data needed to find new positions
- ttpos <- data.frame(x = rep(NA, length(onland)), y = rep(NA, length(onland)))
- ttpos[stuck,] <- fpos[stuck,] # If they're stuck, pass on starting points
- ttcell[stuck] <- fcell[stuck] # If they're stuck, pass on starting cells
-
- if(length(unstuck) > 0) {
- ttdat <- data.frame(ttcell = ttcell[unstuck],
- loncent = terra::xFromCell(vel, ttcell[unstuck]),
- latcent = terra::yFromCell(vel, ttcell[unstuck])) %>% # Start building coordinates
- dplyr::mutate(e = NA, w = NA, n = NA, s = NA, dlon = NA, dlat = NA) # To facilitate finding corners, if we need them
-
- # Separate cells in the coastal blind spot from those that are not
- # which_coast <- which(ttdat$ttcell %in% fn_lk_up)
- # which_not_coast <- which(!(ttdat$ttcell %in% fn_lk_up))
- # For non-coastal cells
-
- idx <- 1:nrow(ttdat)
-
- if(nrow(ttdat) > 0) {
- corner_block_size <- 0.25*x_res
- nc_ttdat <- ttdat # Copying data
- nc_ttdat$e <- nc_ttdat$loncent + (0.5 * x_res) - (stats::runif(1, 0, 1) * corner_block_size)
- nc_ttdat$w <- nc_ttdat$loncent - (0.5 * x_res) + (stats::runif(1, 0, 1) * corner_block_size)
- nc_ttdat$n <- nc_ttdat$latcent + (0.5 * y_res) - (stats::runif(1, 0, 1) * corner_block_size)
- nc_ttdat$s <- nc_ttdat$latcent - (0.5 * y_res) + (stats::runif(1, 0, 1) * corner_block_size)
- coords <- with(nc_ttdat, cbind(n, e, n, w, s, w, s, e)) # NE, NW, SW, SE corners' coordinates
- # Find distances from departure point to corners of ttcell
- get_dist <- function(y1, x1, y2, x2) {
- pt1 <- sf::st_as_sf(data.frame(x = x1, y = y1), coords = c("x", "y"), crs = terra::crs(terra::rast()))
- pt2 <- sf::st_as_sf(data.frame(x = x2, y = y2), coords = c("x", "y"), crs = terra::crs(terra::rast()))
- out <- sf::st_distance(pt1, pt2, by_element = TRUE) %>%
- as.vector()
- return(out)
- }
-
- corners <- data.frame(ne = get_dist(fpos[idx,2], fpos[idx,1], coords[,1], coords[,2])) %>%
- dplyr::mutate(nw = get_dist(fpos[idx,2], fpos[idx,1], coords[,3], coords[,4]),
- sw = get_dist(fpos[idx,2], fpos[idx,1], coords[,5], coords[,6]),
- se = get_dist(fpos[idx,2], fpos[idx,1], coords[,7], coords[,8]))
- cornset <- apply(corners, 1, mfind)*2 # Identify which corners for each onland cell. Have to mul by 2 to shift along correctly.
- cornset <- cbind(cornset, coords) # Add in coordinates
- ttdat[,8:9] <- data.frame(t(apply(cornset, 1, mplace))) # Extract coordinates of correct corner
- }
- # # For coastal cells
- # if(length(which_coast) > 0) {
- # coastal_dat <- ttdat[which_coast,2:3] %>%
- # mutate(x = fpos[which_coast,1], y = fpos[which_coast,2],
- # tccell = ttdat[which_coast,]$ttcell)
- # # JDE TODO - Argument to get_nearest_coast_cell are not correct - `get_nearest_coastal_cell <- function(x, y, tccell, lk_up, ...)`
- # ttdat[which_coast,8:9] <- pmap_df(coastal_dat, get_nearest_coastal_cell)
- # }
- ttpos[unstuck,] <- ttdat[,8:9]
- ttcell[unstuck] <- terra::cellFromXY(mn, ttpos[unstuck,])
- }
- # Collect results
- lonlat[onland,] <- ttpos
- tcells[onland] <- ttcell
- sflags[onland] <- SFlags
- }
- # }
- llon[((i * n) + 1): ((i * n) + n)] <- lonlat[,1] # Add lon to the list
- llat[((i * n) + 1): ((i * n) + n)] <- lonlat[,2]# Add lat to the list
- cellIDend[((i * n) + 1): ((i * n) + n)] <- tcells # Add cellIDs to the list
- flags[((i * n) + 1): ((i * n) + n)] <- sflags # Add flags to the list
- Steps[((i * n) + 1): ((i * n) + n)] <- i # Capture the time_step
- # print(c(i, length(onland), round(100*i/(tyr/tstep), 3), date())) # Keep a reference of progress
+ # Force garbage collection to free memory immediately
+ gc()
- }
- return(cbind(Steps, llon, llat, cellIDs, cellIDend, flags) %>%
- tibble::as_tibble())
+ return(result)
}
-
diff --git a/data-raw/DATASET.R b/data-raw/DATASET.R
index 2a25e19..e2a0bc0 100644
--- a/data-raw/DATASET.R
+++ b/data-raw/DATASET.R
@@ -4,8 +4,9 @@
# Use the existing ones for now. I will speak to Jorge.
load("data-raw/EEZ.RData")
-EEZ <- terra::vect(EEZ)
-terra::writeVector(EEZ, filename = "inst/extdata/EEZ.gpkg", overwrite = TRUE)
+EEZ <- sf::st_as_sf(EEZ) %>%
+ sf::st_make_valid() %>%
+ sf::st_write("inst/extdata/EEZ.gpkg", append = FALSE)
load("data-raw/HSST.RData")
HSST <- terra::rast(HSST)
diff --git a/docs/articles/VoCC.html b/docs/articles/VoCC.html
index 7792bf8..f65f866 100644
--- a/docs/articles/VoCC.html
+++ b/docs/articles/VoCC.html
@@ -287,17 +287,17 @@
# prepare raster layers
vel <- gv[[1]]
ang <- gv[[2]]
-mn <- app(r, mean, na.rm = T)
+mn <- terra::app(r, mean, na.rm = T)
# generate a velocity layer centered and cropped to study region to extract the initial coordinates for the trajectories from
-x1 <- crop(gv[[1]], ext(-180, 0, -90, 90))
-x2 <- crop(gv[[1]], ext(0, 180, -90, 90))
-ext(x1) <- c(180, 360, -90, 90)
-velc <- merge(x1, x2)
+x1 <- terra::crop(gv[[1]], terra::ext(-180, 0, -90, 90))
+x2 <- terra::crop(gv[[1]], terra::ext(0, 180, -90, 90))
+terra::ext(x1) <- c(180, 360, -90, 90)
+velc <- terra::merge(x1, x2)
# crop to the desired extent
# display restricted to +180 longitude to avoid plotting issues with date line crossing
-velc <- crop(velc, c(90, 180, -32, 33))
+velc <- terra::crop(velc, c(90, 180, -32, 33))
We can now populate the data frame with the cell centroid coordinates
for the trajectories.
@@ -312,14 +312,15 @@
-traj <- voccTraj(lonlat, vel, ang, mn, tyr = 50, tstep = 1/12)
+traj <- voccTraj(lonlat, vel, ang, mn, x_res = terra::res(vel)[1], y_res = terra::res(vel)[2], tyr = 50, tstep = 1/12)
Plot them over the climate velocities and the EEZ polygons from the
EEZs data set (Fig. 3a in Garcia Molinos et al. 2019)
# create the spatial object with the trajectories and plot them together with the EEZ polygons
-lns <- map(traj %>% group_split(cellIDs), trajLine) %>%
- purrr::list_rbind() %>%
- sf::st_sf()
+lns <- trajLine(traj)
+
+# mirai::daemons(0) # Set daemons to 0 to close cluster
+
# Load and simplify polygons to speed plotting up
EEZs <- sf::st_read(system.file("extdata", "EEZs.gpkg", package = "VoCCdata")) %>%
@@ -348,7 +349,6 @@