Merge pull request #1 from MathMarEcol/devel_JDE

Devel jde

Author: jaseeverett

.gitignore

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

DESCRIPTION

@@ -21,7 +21,7 @@ Imports:
     geosphere,
     ggplot2,
     magrittr,
-    parallel,
+    parallelly,
     RColorBrewer,
     rlang,
     sp,
@@ -33,16 +33,18 @@ Suggests:
     knitr,
     mapplots,
     ncdf4,
-    prettydoc,
+    patchwork,
     rasterVis,
     repmis,
     rmarkdown,
-    scales
+    scales,
+    tidyterra,
+    VoCCdata
 URL: https://github.com/JorGarMol/VoCC
 BugReports: https://github.com/JorGarMol/VoCC/issues
 License: AGPL (>= 3)
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 VignetteBuilder: knitr,
     rmarkdown

Meta/vignette.rds

@@ -1,7 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
 export("%>%")
-export(VoCC_get_data)
 export(climPCA)
 export(climPlot)
 export(dVoCC)
@@ -11,10 +10,8 @@ export(shiftTime)
 export(spatGrad)
 export(sumSeries)
 export(tempTrend)
-export(trajClas)
-export(trajLine)
-export(voccTraj)
 importFrom(data.table,":=")
 importFrom(foreach,"%dopar%")
 importFrom(magrittr,"%>%")
 importFrom(rlang,.data)
+importFrom(stats,na.omit)

NAMESPACE

@@ -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/climPCA.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/climPlot.R

@@ -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
@@ -57,8 +58,6 @@
 #' 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
@@ -82,9 +81,9 @@
 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()
   }
@@ -99,17 +98,24 @@ 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)
+  # 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
 
-  doParallel::registerDoParallel(cl)
+  # 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)
 
-  result <- foreach::foreach(x = 1:ncores, .combine = rbind, .packages = c("terra", "gdistance", "geosphere", "data.table"), .multicombine = TRUE) %dopar% {
-    a <- x
-    Dat <- dat[cuts == a, ]
+    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,
@@ -177,7 +183,62 @@ dVoCC <- function(clim, n, tdiff, method = "Single", climTol, geoTol,
     }
     return(resu)
   }
-  parallel::stopCluster(cl)
+  } 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
+      }
+
+      # 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]
+        }
+      }
+    }
+  }
 
   return(result)
 }

R/dVoCC.R

@@ -19,7 +19,7 @@
 #' @author Jorge Garcia Molinos
 #'
 #' @examples
-#'
+#' \dontrun{
 #' HSST <- VoCC_get_data("HSST.tif")
 #' yrSST <- sumSeries(HSST,
 #'   p = "1960-01/2009-12", yr0 = "1955-01-01", l = terra::nlyr(HSST),
@@ -32,6 +32,7 @@
 #'
 #' v <- gVoCC(tr, sg)
 #' terra::plot(v)
+#' }
 #'
 gVoCC <- function(tempTrend, spatGrad) {
   VoCC <- tempTrend[[1]] / spatGrad[[1]]

R/gVoCC.R

@@ -21,10 +21,12 @@
 #' @export
 #' @author Jorge Garcia Molinos and Michael T. Burrows
 #' @examples
+#' \dontrun{
 #' HSST <- VoCC_get_data("HSST.tif")
 #' Apr <- shiftTime(HSST, yr1 = 1960, yr2 = 2009, yr0 = 1955, th = 10, m = 4)
 #'
 #' terra::plot(Apr)
+#' }
 shiftTime <- function(r, yr1, yr2, yr0, th, m) {
   # 1. Long term trends in monthly values (e.g. deg/year if temperature)
   m1 <- ((yr1 - yr0) * 12) + m

R/shiftTime.R

@@ -21,11 +21,12 @@
 #' @seealso{\code{\link{tempTrend}}, \code{\link{gVoCC}}}
 #'
 #' @importFrom rlang .data
+#' @importFrom stats na.omit
 #'
 #' @export
 #' @author Jorge Garcia Molinos, David S. Schoeman, and Michael T. Burrows
 #' @examples
-#'
+#' \dontrun{
 #' HSST <- VoCC_get_data("HSST.tif")
 #'
 #' yrSST <- sumSeries(HSST,
@@ -38,13 +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
   sy <- sx <- NSgrad <- WEgrad <- NULL
   clim <- climE <- climN <- climNE <- climNW <- climS <- climSE <- climSW <- climW <- climFocal <- NULL
-  to <- code <- i.to <- LAT <- angle <- Grad <- NULL
+  to <- code <- i.to <- LAT <- angle <- Grad <- .SD <- NULL
 
   if (terra::nlyr(r) > 1) {
     r <- terra::mean(r, na.rm = TRUE)
@@ -55,26 +59,18 @@ spatGrad <- function(r, th = -Inf, projected = FALSE) {
 
   # 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 <- stats::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)
-
-  # TODO JDE added in na.rm = TRUE as I was getting NaN. I can't test if this behaviour has changed from raster....
-  # On second thought I am not sure if NAs are valid here. It gives errors below when calculating weighted means
-  y[, clim := terra::values(r,  na.rm = TRUE)[to]] # Insert values for adjacent cells
+  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[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
 
   # Code cells with positions
-  y[
-    list(
-      code = c("10", "-10", "-11", "-1-1", "11", "1-1", "01", "0-1"),
-      to = c("climE", "climW", "climNW", "climSW", "climNE", "climSE", "climN", "climS")
-    ),
-    on = "code",
-    code := i.to
-  ]
+  y[.(code = c("10", "-10", "-11", "-1-1", "11", "1-1", "01", "0-1"),
+      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
@@ -93,7 +89,8 @@ spatGrad <- function(r, th = -Inf, projected = FALSE) {
   y[, gradWE5 := (climE - climFocal) / (cos(co * CircStats::rad(LAT)) * (d * re[1]))]
   y[, gradWE6 := (climSE - climS) / (cos(co * CircStats::rad(LAT - re[2])) * (d * re[1]))]
 
-  # NS gradients difference in temperatures for each northern and southern pairs divided by the distance between them (111.325 km per degC *re[2] degC)
+  # NS gradients difference in temperatures for each northern and southern pairs divided by
+  # the distance between them (111.325 km per degC *re[2] degC)
   # Positive values indicate an increase in sst from S to N (i.e., in line with the Cartesian y axis)
   y[, gradNS1 := (climNW - climW) / (d * re[2])]
   y[, gradNS2 := (climN - climFocal) / (d * re[2])]
@@ -102,29 +99,16 @@ spatGrad <- function(r, th = -Inf, projected = FALSE) {
   y[, gradNS5 := (climFocal - climS) / (d * re[2])]
   y[, gradNS6 := (climE - climSE) / (d * re[2])]
 
-
-  for (nn in 1:365){
-
-   print(nn)
-
-    print(stats::weighted.mean(y[nn,12:17], w = c(1, 2, 1, 1, 2, 1), na.rm = TRUE))
-
-
-  }
-
-
-  browser()
-  # Calulate NS and WE gradients. 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.
-  # JDE Notes: 1 in apply = operate over rows
-  # Lots of NAs in clim. Can these be removed? Should they be? Chat to Dave S
-  y[, WEgrad := apply(data.table::.SD, 1, function(x) stats::weighted.mean(x, w = c(1, 2, 1, 1, 2, 1), na.rm = TRUE)), .SDcols = gradWE1:gradWE6]
-  y[, NSgrad := apply(data.table::.SD, 1, function(x) stats::weighted.mean(x, c(1, 2, 1, 1, 2, 1), na.rm = T)), .SDcols = 18:23]
+  # Calulate NS and WE gradients.
+  # 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]
   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(data.table::.SD$WEgrad, data.table::.SD$NSgrad), .SDcols = c("WEgrad", "NSgrad")]
+  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))]
@@ -139,5 +123,6 @@ spatGrad <- function(r, th = -Inf, projected = FALSE) {
   rGrad[rGrad[] < th] <- th
   output <- c(rGrad, rAng)
   names(output) <- c("Grad", "Ang")
+
   return(output)
 }