Merge pull request #122 from r-spatialecology/main

Update
r-spatialecology · Oct 2, 2023 · b3fedd0 · b3fedd0
2 parents e7c40a0 + 676f5a3
commit b3fedd0
Show file tree

Hide file tree

Showing 63 changed files with 899 additions and 2,737 deletions.
diff --git a/.github/workflows/Render-README.yaml b/.github/workflows/Render-README.yaml
@@ -38,4 +38,4 @@ jobs:
           git config --local user.email "[email protected]"
           git add README.md man/figures/README-*
           git commit -m "Re-build README.md" || echo "No changes to commit"
-          git push origin || echo "No changes to commit"
+          git push -f origin || echo "No changes to commit"
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: shar
 Title: Species-Habitat Associations
-Version: 2.0.4
+Version: 2.1
 Authors@R: c(person("Maximilian H.K.", "Hesselbarth", email = "[email protected]",
              role = c("aut", "cre"), comment = c(ORCID = "0000-0003-1125-9918")),
              person("Marco", "Sciaini", email = "[email protected]",
@@ -30,8 +30,8 @@ Imports:
     grDevices,
     methods,
     spatstat.explore,
-    spatstat.model,
     spatstat.geom,
+    spatstat.model,
     spatstat.random,
     stats,
     terra,

diff --git a/NAMESPACE b/NAMESPACE
@@ -8,7 +8,6 @@ S3method(print,rd_pat)
 S3method(print,rd_ras)
 export(calculate_energy)
 export(classify_habitats)
-export(estimate_pcf_fast)
 export(fit_point_process)
 export(list_to_randomized)
 export(pack_randomized)

diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,12 @@
+# shar 2.1
+* Improvements
+  * Remove `comp_fast` argument
+  * Speed improvements of computation
+  * General updates to code structure 
+* Bugfixes
+  * Removed `n_points` and `window` argument from reconstruction due to methodological issues
+  * Bugfix related to wrap/unwrap raster and printing
+
 # shar 2.0.4
 * Improvements
   * Remove `Sys.sleep` for verbose reconstruction

diff --git a/R/calc_gest.R b/R/calc_gest.R
@@ -0,0 +1,32 @@
+#' calc_gest
+#'
+#' @description Calculate Gest
+#'
+#' @param dist matrix with distance pairs.
+#' @param r vector with distances r.
+#' @param n_points numeric with number of points
+#'
+#' @details
+#' Calculates Gest based on distances created with \code{get_dist_pairs}.
+#'
+#' @seealso
+#' \code{\link{get_dist_pairs}}
+#'
+#' @return data.frame
+#'
+#' @aliases calc_gest
+#' @rdname calc_gest
+#'
+#' @keywords internal
+calc_gest <- function(dist, r, n_points){
+
+  mat <- matrix(nrow = n_points, ncol = n_points, data = Inf)
+  mat[dist[, 1:2]] <- dist[, 3]
+
+  distances_min <- apply(X = mat, MARGIN = 2, FUN = min, na.rm = TRUE)
+
+  hist_min <- graphics::hist(distances_min, breaks = r, plot = FALSE)
+
+  data.frame(r = hist_min$mids, edf = cumsum(hist_min$counts) / n_points)
+
+}
diff --git a/R/calculate_energy.R b/R/calculate_energy.R
@@ -6,18 +6,13 @@
 #' @param weights Vector with weights used to calculate energy.
 #' The first number refers to Gest(r), the second number to pcf(r).
 #' @param return_mean Logical if the mean energy is returned.
-#' @param comp_fast Integer with threshold at which summary functions are estimated
-#' in a computational fast way.
 #' @param verbose Logical if progress report is printed.
 #'
 #' @details
 #' The function calculates the mean energy (or deviation) between the observed
 #' pattern and all reconstructed patterns (for more information see Tscheschel &
 #' Stoyan (2006) or Wiegand & Moloney (2014)). The pair correlation function and the
-#' nearest neighbour distance function are used to describe the patterns. For large
-#' patterns \code{comp_fast = TRUE} decreases the computational demand, because no edge
-#' correction is used and the pair correlation function is estimated based on Ripley's
-#' K-function. For more information see \code{\link{estimate_pcf_fast}}.
+#' nearest neighbour distance function are used to describe the patterns.
 #'
 #' @seealso
 #' \code{\link{plot_energy}} \cr
@@ -54,9 +49,8 @@
 #'
 #' @export
 calculate_energy <- function(pattern,
-                             weights = c(0.5, 0.5),
+                             weights = c(1, 1),
                              return_mean = FALSE,
-                             comp_fast = 1000,
                              verbose = TRUE){
 
   # check if class is correct
@@ -83,7 +77,7 @@ calculate_energy <- function(pattern,
   # calculate r sequence
   r <- seq(from = 0,
            to = spatstat.explore::rmax.rule(W = pattern_observed$window,
-                                         lambda = spatstat.geom::intensity.ppp(pattern_observed)),
+                                            lambda = spatstat.geom::intensity.ppp(pattern_observed)),
            length.out = 250)
 
   if (inherits(x = pattern, what = "rd_pat")) {
@@ -96,71 +90,22 @@ calculate_energy <- function(pattern,
 
     } else {
 
-      # check if weights make sense
-      if (sum(weights) > 1 || sum(weights) == 0) {
-
-        stop("The sum of 'weights' must be 0 < sum(weights) <= 1.", call. = FALSE)
-
-      }
-
-      # check if number of points exceed comp_fast limit
-      if (pattern_observed$n > comp_fast) {
-
-        comp_fast <- TRUE
-
-      } else {
-
-        comp_fast <- FALSE
-
-      }
-
       # calculate summary functions for observed pattern
-      if (comp_fast) {
-
-        gest_observed <- spatstat.explore::Gest(X = pattern_observed, correction = "none",
-                                             r = r)
-
-        pcf_observed <- estimate_pcf_fast(pattern = pattern_observed,
-                                          correction = "none", method = "c",
-                                          spar = 0.5, r = r)
+      gest_observed <- spatstat.explore::Gest(X = pattern_observed,
+                                              correction = "none", r = r)
 
-      } else {
-
-        gest_observed <- spatstat.explore::Gest(X = pattern_observed,
-                                        correction = "han", r = r)
-
-        pcf_observed <- spatstat.explore::pcf(X = pattern_observed,
-                                      correction = "best", divisor = "d", r = r)
-
-      }
+      pcf_observed <- spatstat.explore::pcf(X = pattern_observed,
+                                            correction = "none", divisor = "d", r = r)
 
       # loop through all reconstructed patterns
       result <- vapply(seq_along(pattern_randomized), function(x) {
 
-        # fast computation of summary stats
-        if (comp_fast) {
-
-          gest_reconstruction <- spatstat.explore::Gest(X = pattern_randomized[[x]],
-                                                correction = "none",
-                                                r = r)
-
-          pcf_reconstruction <- estimate_pcf_fast(pattern = pattern_randomized[[x]],
-                                                  correction = "none", method = "c",
-                                                  spar = 0.5, r = r)
-
-        # normal computation of summary stats
-        } else {
-
-          gest_reconstruction <- spatstat.explore::Gest(X = pattern_randomized[[x]],
-                                                correction = "han",
-                                                r = r)
+        gest_reconstruction <- spatstat.explore::Gest(X = pattern_randomized[[x]],
+                                                      correction = "none", r = r)
 
           pcf_reconstruction <- spatstat.explore::pcf(X = pattern_randomized[[x]],
-                                              correction = "best",
-                                              divisor = "d",
-                                              r = r)
-
-        }
+                                                      correction = "none", divisor = "d",
+                                                      r = r)
 
         # difference between observed and reconstructed pattern
         energy <- (mean(abs(gest_observed[[3]] - gest_reconstruction[[3]]), na.rm = TRUE) * weights[[1]]) +

diff --git a/R/data.R b/R/data.R
@@ -34,31 +34,3 @@
 #'
 #' @format A spatstat ppp object.
 "species_b"
-
-#' Gamma test
-#'
-#' Randomized data for species b using the gamma test.
-#'
-#' @format rd_pat object.
-"gamma_test"
-
-#' Reconstruction
-#'
-#' Randomized data for species b using pattern reconstruction.
-#'
-#' @format rd_pat object.
-"reconstruction"
-
-#' Torus trans
-#'
-#' Torus translation of the classified \code{landscape}.
-#'
-#' @format rd_ras object.
-"torus_trans"
-
-#' Random walk
-#'
-#' Randomization of the \code{landscape} using the habitat randomization algorithm.
-#'
-#' @format rd_ras object.
-"random_walk"
diff --git a/R/estimate_pcf_fast.R b/R/estimate_pcf_fast.R
diff --git a/R/get_dist_pairs.R b/R/get_dist_pairs.R
@@ -0,0 +1,26 @@
+#' get_dist_pairs
+#'
+#' @description Distance between points
+#'
+#' @param X ppp object
+#' @param rmax Numeric with maximum distance
+#'
+#' @details
+#' Returns matrix with point pairs and distances between them.
+#'
+#' @seealso
+#' \code{\link{pcf.ppp}}
+#'
+#' @return matrix
+#'
+#' @aliases get_dist_pairs
+#' @rdname get_dist_pairs
+#'
+#' @keywords internal
+get_dist_pairs <- function(X, rmax){
+
+  dist_observed <- spatstat.geom::closepairs(X = X, rmax = rmax, what = "ijd", twice = TRUE)
+
+  cbind(dist_observed$i, dist_observed$j, dist_observed$d)
+
+}
diff --git a/R/pack_randomized.R b/R/pack_randomized.R
@@ -27,8 +27,9 @@
 #' @export
 pack_randomized <- function(raster) {
 
-  # wrap observerd raster
-  raster$observed <- terra::wrap(raster$observed)
+  # check if observed is present
+  # wrap observed raster
+  if (inherits(x = raster$observed, what = "SpatRaster")) raster$observed <- terra::wrap(raster$observed)
 
   # wrap all randomized raster
   raster$randomized <- lapply(X = raster$randomized, FUN = terra::wrap)

diff --git a/R/plot.rd_mar.R b/R/plot.rd_mar.R
@@ -8,8 +8,6 @@
 #' @param n Integer with number or vector of ids of randomized pattern to plot.
 #' See Details section for more information.
 #' @param probs Vector with quantiles of randomized data used for envelope construction.
-#' @param comp_fast Integer with threshold at which summary functions are estimated
-#' in a computational fast way.
 #' @param ask Logical if the user is asked to press <RETURN> before second summary function
 #' is plotted (only used if \code{what = "sf"}).
 #' @param verbose Logical if progress report is printed.
@@ -18,9 +16,6 @@
 #' @details
 #' The function plots the pair correlation function and the nearest neighbour function of
 #' the observed pattern and the reconstructed patterns (as "simulation envelopes").
-#' For large patterns \code{comp_fast = TRUE} decreases the computational demand because no edge
-#' correction is used and the pair correlation function is estimated based on Ripley's
-#' K-function. For more information see \code{\link{estimate_pcf_fast}}.
 #'
 #' It is also possible to plot n randomized patterns and the observed pattern
 #' using \code{what = "pp"}. If \code{n} is a single number, \code{n} randomized
@@ -48,7 +43,7 @@
 #' @rdname plot.rd_mar
 #'
 #' @export
-plot.rd_mar <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975), comp_fast = 1000,
+plot.rd_mar <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975),
                         ask = TRUE, verbose = TRUE, ...) {
 
   # check if class is correct
@@ -70,17 +65,6 @@ plot.rd_mar <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975), comp_
 
   if (what == "sf") {
 
-    # check if number of points exceed comp_fast limit
-    if (x$observed$n > comp_fast) {
-
-      comp_fast <- TRUE
-
-    } else {
-
-      comp_fast <- FALSE
-
-    }
-
     name_unit <- spatstat.geom::unitname(x$observed)[[1]] # unit name for labels
 
     # calculate r
@@ -135,8 +119,7 @@ plot.rd_mar <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975), comp_
     # specify quantums g(r)
     col_kmmr <- ifelse(test = result_observed[, 3] < result_randomized[, 2] |
                          result_observed[, 3] > result_randomized[, 3],
-                       yes =  "#1f78b4",
-                       no = "#b2df8a")
+                       yes =  "#1f78b4", no = "#b2df8a")
 
     # plot results
     graphics::plot(NULL, xlim = range(r), ylim = yrange,
@@ -195,8 +178,11 @@ plot.rd_mar <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975), comp_
       # convert to dataframe
       current_pattern <- as.data.frame(subset_pattern[[i]])
 
+      current_pattern$marks <- ((current_pattern$marks - min(current_pattern$marks)) /
+        (max(current_pattern$marks) - min(current_pattern$marks)) * 1) + 0.25
+
       # plot points
-      graphics::plot(x = current_pattern$x, y = current_pattern$y,
+      graphics::plot(x = current_pattern$x, y = current_pattern$y, cex = current_pattern$marks,
                      type = "p", asp = 1, xlim = x_range, ylim = y_range, axes = FALSE,
                      main  = names_pattern[[i]], xlab = "", ylab = "")