Merge pull request #8 from MathMarEcol/devel_JDE

Devel jde

Author: jaseeverett

DESCRIPTION

@@ -9,34 +9,33 @@ Authors@R: c(person("Jorge", "Garcia Molinos", email = "[email protected]
 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,
     stats,
     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

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)

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

R/climPCA.R

@@ -20,35 +20,49 @@
 #' @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)
 #' summary(comp[[1]]) # first two components explain >90% of variance
 #' # 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

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(