Pulling new changes to branch

.Rbuildignore

@@ -10,3 +10,6 @@
 ^data-raw$
 ^README\.Rmd$
 ^LICENSE\.md$
+^_pkgdown\.yml$
+^docs$
+^pkgdown$

.gitignore

@@ -2,3 +2,5 @@
 .Rhistory
 .RData
 .DS_Store
+/doc/
+/Meta/

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: VoCC
 Title: The Velocity of Climate Change and related climatic metrics
-Version: 1.0.1
+Version: 0.0.1
 Authors@R: c(person("Jorge", "Garcia Molinos", email = "[email protected]", role = c("aut", "cre")),
     person("David", "S. Schoeman", email = "", role = "aut"),
     person("Christopher", "J. Brown", email = "", role = "aut"),
@@ -10,39 +10,42 @@ Description: Functions to calculate the velocity of climate change (VoCC) and re
     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)
-Imports: 
+Imports:
     assertthat,
     CircStats,
     cowplot,
     data.table,
     doParallel,
+    dplyr,
     foreach,
-    gdistance,
     geosphere,
-    ggplot2,
     magrittr,
-    parallel,
+    parallelly,
     RColorBrewer,
     rlang,
-    sp,
+    sf,
     stats,
     terra,
+    tibble,
+    tidyr,
+    tidyselect,
     xts
 Suggests:
-    gridExtra,
+    ggplot2,
     knitr,
     mapplots,
     ncdf4,
-    prettydoc,
-    rasterVis,
-    repmis,
+    patchwork,
+    purrr,
     rmarkdown,
-    scales
-URL: https://github.com/JorGarMol/VoCC
-BugReports: https://github.com/JorGarMol/VoCC/issues
+    scales,
+    tidyterra,
+    VoCCdata
+URL: https://mathmarecol.github.io/VoCC/
+BugReports: https://github.com/MathMarEcol/VoCC/issues
 License: AGPL (>= 3)
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 VignetteBuilder: knitr,
     rmarkdown

NAMESPACE

@@ -1,7 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
 export("%>%")
-export(VoCC_get_data)
 export(climPCA)
 export(climPlot)
 export(dVoCC)
@@ -18,3 +17,4 @@ importFrom(data.table,":=")
 importFrom(foreach,"%dopar%")
 importFrom(magrittr,"%>%")
 importFrom(rlang,.data)
+importFrom(stats,na.omit)

R/climPCA.R

@@ -22,7 +22,7 @@
 #' @export
 #' @author Jorge Garcia Molinos
 #' @examples
-#'
+#' \dontrun{
 #' JapTC <- VoCC_get_data("JapTC.tif")
 #'
 #' comp <- climPCA(JapTC[[c(1, 3, 5)]], JapTC[[c(2, 4, 6)]],
@@ -31,6 +31,8 @@
 #' # 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) {
   # 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)))

R/climPlot.R

@@ -17,7 +17,7 @@
 #' @export
 #' @author Jorge Garcia Molinos and Naoki H. Kumagai
 #' @examples
-#'
+#' \dontrun{
 #' JapTC <- VoCC_get_data("JapTC.tif")
 #'
 #' # Plot climate space for the two first variables(annual precipitation and maximum temperature)
@@ -32,7 +32,6 @@
 #'   y.name = "Temperature max (°C)"
 #' )
 #'
-#' \dontrun{
 #' # output plots can be saved as:
 #' ggplot2::ggsave(
 #'   plot = out, filename = file.path(getwd(), "example_plot.pdf"),
@@ -44,6 +43,7 @@ climPlot <- function(xy, x.binSize, y.binSize, x.name = "V1", y.name = "V2") {
   yp <- xy[, 3]
   xf <- xy[, 2]
   yf <- xy[, 4]
+
   # 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)
@@ -102,7 +102,7 @@ climPlot <- function(xy, x.binSize, y.binSize, x.name = "V1", y.name = "V2") {
   Freq2D <- data.frame(x = x.bin, y = rep(y.bin, each = length(x.bin)), freq = Freq2D)
   Freq2D <- Freq2D[!is.na(Freq2D$freq), ]
 
-  panelAB <- ggplot2::ggplot(Freq2Dpf, ggplot2::aes(x = x, y = y, fill = freq)) +
+  panelAB <- ggplot2::ggplot(Freq2Dpf, ggplot2::aes(x = .data$x, y = .data$y, fill = freq)) +
     ggplot2::geom_raster() +
     ggplot2::scale_fill_gradientn(
       colors = r,
@@ -121,7 +121,7 @@ climPlot <- function(xy, x.binSize, y.binSize, x.name = "V1", y.name = "V2") {
       strip.text = ggplot2::element_blank()
     )
 
-  panelC <- ggplot2::ggplot(Freq2D, ggplot2::aes(x = x, y = y, fill = freq)) +
+  panelC <- ggplot2::ggplot(Freq2D, ggplot2::aes(x = .data$x, y = .data$y, fill = freq)) +
     ggplot2::geom_raster() +
     ggplot2::scale_fill_manual(values = c("#56B4E9", "#009E73", "#D55E00"), name = "Climate type") +
     ggplot2::labs(x = x.name, y = y.name)

R/dVoCC.R

@@ -14,7 +14,7 @@
 #' These columns are not required if using the "Single" method. The last three columns of the table should contain an identifyier and centroid coordinates of each cell.
 #' @param n \code{integer} defining the number of climatic variables.
 #' @param tdiff \code{integer} defining the number of years (or other temporal unit) between periods.
-#' @param method \code{ch)aracter string} specifying the analogue method to be used. 'Single': a constant, single analogue threshold
+#' @param method \code{character string} specifying the analogue method to be used. 'Single': a constant, single analogue threshold
 #' for each climate variable is applied to all cells (Ohlemuller et al. 2006, Hamann et al. 2015); climate analogy corresponds to target cells
 #' with values below the specified threshold for each climatic variable. 'Variable': a cell-specific climate threshold is used for each climatic variable
 #' to determine the climate analogues associated with each cell by reference to its baseline climatic variability (Garcia Molinos et al. 2017).
@@ -41,6 +41,7 @@
 #' @export
 #' @author Jorge Garcia Molinos
 #' @examples
+#' \dontrun{
 #' JapTC <- VoCC_get_data("JapTC.tif")
 #'
 #' # Create a data frame with the necessary variables in the required order
@@ -53,15 +54,13 @@
 #'   geoTol = 160, distfun = "GreatCircle", trans = NA, lonlat = TRUE
 #' )
 #'
-#' r1 <- JapTC
+#' r1 <- JapTC[[1]]
 #' 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
+#' r <- JapTC[[1]]
 #' r[!is.na(JapTC[[1]])] <- 1
 #' h8 <- gdistance::transition(r, transitionFunction = mean, directions = 8)
 #' h8 <- gdistance::geoCorrection(h8, type = "c")
@@ -73,7 +72,7 @@
 #' )
 #'
 #' # Plot results
-#' r1 <- r2 <- JapTC
+#' r1 <- r2 <- JapTC[[1]]
 #' r1[avocc1$focal] <- avocc1$vel
 #' r2[avocc2$focal] <- avocc2$vel
 #' terra::plot(c(r1, r2))
@@ -82,13 +81,17 @@
 dVoCC <- function(clim, n, tdiff, method = "Single", climTol, geoTol,
                   distfun = "GreatCircle", trans = NA, lonlat = TRUE) {
 
-  geoDis <- climDis <- ang <- vel <- target <- cid <- NULL # Fix devtools check warnings
+  geoDis <- climDis <- ang <- vel <- target <- cid <- a <- NULL # Fix devtools check warnings
 
-  if (distfun == "Euclidean" & lonlat == TRUE) {
+  if (distfun == "Euclidean" && lonlat == TRUE) {
     print("Error: Euclidean distances specified for unprojected coordinates")
     stop()
   }
 
+  if (distfun == "LeastCost") {
+    stop("LeastCost distances are not currently supported. Use 'Euclidean' or 'GreatCircle' instead.")
+  }
+
   assertthat::assert_that(
     all(distfun %in% c("Euclidean", "GreatCircle")),
     is.na(trans)
@@ -99,33 +102,99 @@ dVoCC <- function(clim, n, tdiff, method = "Single", climTol, geoTol,
   # matrix with the future climatic values for all cells
   fut <- dat[, seq(2, (2 * n), by = 2), with = FALSE]
 
-  # set things up for parallel processing
-  cores <- parallel::detectCores()
-  ncores <- cores[1] - 1
-  cuts <- cut(1:nrow(dat), ncores, labels = FALSE)
-  cl <- parallel::makeCluster(ncores)
-
-  doParallel::registerDoParallel(cl)
-
-  result <- foreach::foreach(x = 1:ncores, .combine = rbind, .packages = c("terra", "gdistance", "geosphere", "data.table"), .multicombine = TRUE) %dopar% {
-    a <- x
-    Dat <- dat[cuts == a, ]
-
-    resu <- data.table::data.table(
-      focal = Dat$cid,
+  # 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)
     )
 
-    i <- 0
-    while (i <= nrow(Dat)) {
-      i <- i + 1
-
+    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])
+      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
@@ -137,7 +206,7 @@ dVoCC <- function(clim, n, tdiff, method = "Single", climTol, geoTol,
       }
 
       if (method == "Variable") { # Garcia Molinos et al. 2017
-        climTol <- as.numeric(Dat[i, ((2 * n) + 1):(3 * n), with = FALSE]) # focal cell tolerance
+        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)
@@ -148,36 +217,26 @@ dVoCC <- function(clim, n, tdiff, method = "Single", climTol, geoTol,
       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]))
+          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
+          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
 
-        # TODO transition matrices not possible at the moment as gdistance has not been updated for terra
-        # if(distfun == "LeastCost"){
-        #   SL <- gdistance::shortestPath(trans, cbind(Dat$x[i],Dat$y[i]), cbind(dat$x[dat$cid %in% anacid], dat$y[dat$cid %in% anacid]), output="SpatialLines")
-        #   d <- SpatialLinesLengths(SL, longlat = lonlat)    # in km for longlat TRUE
-        #   # correct for analogues that are within search distance but have no directed path with the focal cell (i.e. conductivity = 0)
-        #   d[which(d == 0 & anacid != Dat$cid[i])] <- Inf
-        # }
-
         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
+          result[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]), ]))]
+            result[i, climDis := mean(as.numeric(dif[which(anacid == result[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]
+          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]
         }
       }
     }
-    return(resu)
   }
-  parallel::stopCluster(cl)
 
   return(result)
 }