Merge pull request #130 from r-spatialecology/main

v2.2

.github/workflows/R-CMD-check.yaml

@@ -17,7 +17,7 @@ jobs:
     name: ${{ matrix.config.os }} (${{ matrix.config.r }})
 
     strategy:
-      fail-fast: false
+      fail-fast: FALSE
       matrix:
         config:
           - {os: macOS-latest, r: 'release'}

.github/workflows/Test-coverage.yaml

@@ -31,5 +31,24 @@ jobs:
           needs: coverage
 
       - name: Test coverage
-        run: covr::codecov(quiet = FALSE)
+        run: |
+          covr::codecov(
+            quiet = FALSE,
+            clean = FALSE,
+            install_path = file.path(normalizePath(Sys.getenv("RUNNER_TEMP"), winslash = "/"), "package")
+          )
         shell: Rscript {0}
+
+      - name: Show testthat output
+        if: always()
+        run: |
+          ## --------------------------------------------------------------------
+          find ${{ runner.temp }}/package -name 'testthat.Rout*' -exec cat '{}' \; || true
+        shell: bash
+
+      - name: Upload test results
+        if: failure()
+        uses: actions/upload-artifact@v3
+        with:
+          name: coverage-test-failures
+          path: ${{ runner.temp }}/package

DESCRIPTION

@@ -1,11 +1,13 @@
 Type: Package
 Package: shar
 Title: Species-Habitat Associations
-Version: 2.1.1
+Version: 2.2
 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]",
              role = "aut", comment = c(ORCID = "0000-0002-3042-5435")), 
+             person("Chris", "Wudel", email = "[email protected]",
+             role = "aut", comment = c(ORCID = "0000-0003-0446-4665")), 
              person("Zeke", "Marshall", email = "[email protected]",
              role = "ctb", comment = c(ORCID = "0000-0001-9260-7827")), 
              person("Thomas", "Etherington", email = "[email protected]",

NAMESPACE

@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(plot,rd_mar)
+S3method(plot,rd_multi)
 S3method(plot,rd_pat)
 S3method(plot,rd_ras)
 S3method(print,rd_mar)
@@ -15,6 +16,7 @@ export(plot_energy)
 export(randomize_raster)
 export(reconstruct_pattern)
 export(reconstruct_pattern_marks)
+export(reconstruct_pattern_multi)
 export(results_habitat_association)
 export(translate_raster)
 export(unpack_randomized)

NEWS.md

@@ -1,3 +1,7 @@
+# shar 2.2
+* Improvements
+  * Added a new function `reconstruct_pattern_multi()` including several internal functions and methods
+
 # shar 2.1.1
 * Bugfixes (thanks to @baddstats)
 

R/calc_moments.R

@@ -0,0 +1,57 @@
+#' calc_moments
+#'
+#' @description Calculate moments
+#'
+#' @param fn Determination of the weightings of the mark correlation functions.
+#' @param p Defines the initial state of the new ponit pattern.
+#' @param x,y x and y coordinates of the points from the reference point pattern.
+#' @param mark Marks the currently viewed point pattern.
+#' @param kernel Result of the kernel calculation, calculated with the calc_kernels function.
+#' @param rmax_bw Maximum distance at which the summary statistics are
+#' evaluated + Bandwidth with which the kernels are scaled, so that this is the
+#' standard deviation of the smoothing kernel.
+#' @param r Sequence from rmin to rmax in rcount steps.
+#' @param exclude Vector indicating which values not to use.
+#'
+#' @details
+#' Definition of the product-moment function for calculating the contribution
+#' of a point at the coordinates x, y with marking.
+#'
+#' @return matrix
+#'
+#' @aliases calc_moments
+#' @rdname calc_moments
+#'
+#' @keywords internal
+calc_moments <- function(fn,
+                         p,
+                         exclude = NULL,
+                         x,
+                         y,
+                         mark,
+                         kernel,
+                         rmax_bw,
+                         r) {
+
+  d2 <- (p$x-x)^2 + (p$y-y)^2
+  use <-  d2 <= rmax_bw^2
+  use[exclude] <- FALSE
+  z <- crossprod(p$mark[use, , drop = FALSE],
+                 outer(sqrt(d2[use]), r, function(d, r) kernel(r, d)))
+  z[fn$i, , drop = FALSE] * mark[fn$j] + z[fn$j, , drop = FALSE] * mark[fn$i]
+}
+
+calc_moments_full <- function(fn,
+                              p,
+                              kernel,
+                              rmax_bw,
+                              r) {
+
+  f <- 0
+  for (i in seq_len(nrow(p))) {
+    f <- f + calc_moments(fn, p, i:nrow(p), p$x[i], p$y[i], p$mark[i, ],
+                          kernel, rmax_bw, r)
+  }
+  rownames(f) <- paste(names(fn$i), names(fn$j), sep = ":")
+  f
+}

R/calculate_energy.R

@@ -94,18 +94,17 @@ calculate_energy <- function(pattern,
       gest_observed <- spatstat.explore::Gest(X = pattern_observed,
                                               correction = "none", r = r)
 
-      pcf_observed <- spatstat.explore::pcf(X = pattern_observed,
-                                            correction = "none", divisor = "d", r = r)
+      pcf_observed <- spatstat.explore::pcf.ppp(X = pattern_observed, correction = "none", divisor = "d", r = r)
 
       # loop through all reconstructed patterns
       result <- vapply(seq_along(pattern_randomized), function(x) {
 
         gest_reconstruction <- spatstat.explore::Gest(X = pattern_randomized[[x]],
                                                       correction = "none", r = r)
 
-          pcf_reconstruction <- spatstat.explore::pcf(X = pattern_randomized[[x]],
-                                                      correction = "none", divisor = "d",
-                                                      r = r)
+          pcf_reconstruction <- spatstat.explore::pcf.ppp(X = pattern_randomized[[x]],
+                                                          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]]) +

R/compute_statistics.R

@@ -0,0 +1,53 @@
+#' compute_statistics
+#'
+#' @description Compute summary statistics
+#'
+#' @param x,y x and y coordinates of the points from the reference point pattern.
+#' @param k Vector of values k; used only if Dk is included in w_statistics below.
+#' @param xr,yr x and y extension of the observation window (start, end).
+#' @param w_statistics  vector of named weights for optional spatial statistics
+#' from the \code{spatstat} package to be included in the energy calculation. This may
+#' include Dk, K, Hs, pcf.
+#' @param bw,divisor,kernel_arg,r Several parameters related to summary function.
+#'
+#' @details
+#' Compute optional spatial statistics using the spatstat package.
+#'
+#' @return list
+#'
+#' @aliases compute_statistics
+#' @rdname compute_statistics
+#'
+#' @keywords internal
+compute_statistics <- function(x, y, k, xr, yr, w_statistics, bw, divisor, kernel_arg, r) {
+
+  stat <- names(w_statistics)
+  names(stat) <- stat
+  lapply(stat, function(name) switch(name,
+    # Calculation of the Dk(r)-function, if this is to be taken into account for the energy calculation.
+    Dk = {
+      nnd_ <- as.matrix(spatstat.geom::nndist(x, y, k=k))
+      apply(nnd_, 2, function(z) cumsum(graphics::hist(z[z <= max(r)], breaks = c(-Inf, max(r)), plot = FALSE) $ count) / length(z))
+    },
+
+    # Calculation of the K(r)-function, if this is to be taken into account for the energy calculation.
+    K = {
+      kest<-spatstat.explore::Kest(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), rmax=max(r), correction="none")
+      kest$un
+    },
+
+    # Calculation of the pcf(r)-function (spherical contact distribution), if this is to be taken into account for the energy calculation.
+    pcf = {
+      pcfest<-spatstat.explore::pcf.ppp(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), r=c(0,r), kernel=kernel_arg, divisor=divisor, bw=bw, correction="none")
+      pcfest$un
+    },
+    # Calculation of the Hs(r)-function (pair correlation function), if this is to be taken into account for the energy calculation.
+    Hs = {
+      hest<-spatstat.explore::Hest(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), correction="none")
+      hest$raw
+    },
+    # wrong selection
+    stop("unknown statistic")
+    )
+  )
+}

R/data.R

@@ -4,7 +4,7 @@
 #' generated with the \code{NLMR::nlm_fbm()} algorithm.
 #'
 #' @format A SpatRaster object.
-#' @source Simulated neutral landscape model with R. https://github.com/ropensci/NLMR/
+#' @source Simulated neutral landscape model with R. <https://github.com/ropensci/NLMR/>
 "landscape"
 
 #' Species a

R/dummy_transf.R

@@ -0,0 +1,29 @@
+#' dummy_transf
+#'
+#' @description Tranfsorm to dummy variables
+#'
+#' @param f Result of the calc_moments_full function which represents
+#' product-moment contribution of a point at coordinates x, y with marks,
+#' for the whole point pattern.
+#'
+#' @details
+#' Function for the transformation of variables to dummy variables and back
+#'
+#' @return matrix
+#'
+#' @aliases dummy_transf
+#' @rdname dummy_transf
+#'
+#' @keywords internal
+to_dummy <- function(f) {
+  x <- matrix(0, length(f), nlevels(f), dimnames=list(names(f), levels(f)))
+  x[cbind(seq_along(f), as.integer(f))] <- 1
+  x
+}
+
+from_dummy <- function(x, levels=colnames(x)) {
+  f <- as.integer(x %*% seq_along(levels))
+  levels(f) <- levels
+  class(f) <- "factor"
+  f
+}

R/energy_fun.R

@@ -0,0 +1,49 @@
+#' energy_fun
+#'
+#' @description Energy function
+#'
+#' @param f Result of the calc_moments_full function which represents
+#' product-moment contribution of a point at coordinates x, y with marks,
+#' for the whole new ponit pattern.
+#' @param f0  Column sums of the weights of the brand correlation functions of
+#' the new point pattern.
+#' @param statistics Results of the compute_statistics function for the
+#' new point pattern (calculation of optional spatial statistics).
+#' @param fn Determination of the weightings of the mark correlation functions.
+#' @param p Defines the initial state of the new ponit pattern.
+#' @param p_ Reference point pattern.
+#' @param Lp Distance measure for the calculation of the energy function
+#' (Lp distance, 1 <= p <Inf).
+#' @param w_statistics Vector of named weights for optional spatial statistics
+#' from the \code{spatstat} package to be included in the energy calculation.This may
+#' include Dk, K, Hs, pcf.
+#'
+#' @details
+#' Defining the Energy_fun function to calculate the "energy" of the pattern
+#' (where a lower energy indicates a better match).
+#'
+#' @return vector
+#'
+#' @aliases energy_fun
+#' @rdname energy_fun
+#'
+#' @keywords internal
+#'
+Energy_fun <- function(f, f0, statistics, f_, f0_, statistics_, fn, p, p_, Lp, w_statistics) {
+    result <- c(
+      f = sum(fn$w * rowMeans(abs(
+        f / nrow(p) -
+        f_ / nrow(p_)
+      )^Lp)),
+      f0 = sum(fn$w0 * abs(
+        f0 / nrow(p) -
+        f0_ / nrow(p_)
+      )^Lp),
+      if (length(w_statistics))
+        sapply(seq_along(w_statistics), function(i) w_statistics[i] *
+          mean(abs(statistics[[i]] - statistics_[[i]])^Lp, na.rm = TRUE),
+          USE.NAMES=FALSE
+        )
+    )
+    c(energy = sum(result), result)
+}