Merge pull request #37 from r-spatialecology/master

Update v0.3
r-spatialecology · Mar 22, 2019 · 9f87277 · 9f87277
2 parents f14ebc5 + 828ee61
commit 9f87277
Show file tree

Hide file tree

Showing 97 changed files with 2,079 additions and 474 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: shar
 Title: Species-Habitat Associations
-Version: 0.2
+Version: 0.3
 Authors@R: c(person("Maximillian H.K.", "Hesselbarth", email = "[email protected]",
              role = c("aut", "cre"), comment = c(ORCID = "0000-0003-1125-9918")),
              person("Marco", "Sciaini", email = "[email protected]",
@@ -18,12 +18,13 @@ URL: https://r-spatialecology.github.io/shar
 BugReports: https://github.com/r-spatialecology/shar/issues
 Depends: R (>= 3.1)
 Imports:
-  classInt,
-  graphics,
-  raster,
-  spatstat, 
-  stats,
-  utils
+    classInt,
+    graphics,
+    raster,
+    spatstat,
+    stats,
+    utils,
+    Rcpp
 RoxygenNote: 6.1.1
 Suggests:
     covr,
@@ -34,3 +35,5 @@ Suggests:
 VignetteBuilder: knitr
 Encoding: UTF-8
 LazyData: true
+LinkingTo: 
+    Rcpp
diff --git a/NAMESPACE b/NAMESPACE
@@ -7,8 +7,12 @@ export(estimate_pcf_fast)
 export(extract_points)
 export(fit_point_process)
 export(plot_randomized_pattern)
+export(plot_randomized_raster)
 export(randomize_raster)
+export(rcpp_sample)
 export(reconstruct_marks)
 export(reconstruct_pattern)
 export(results_habitat_association)
 export(translate_raster)
+importFrom(Rcpp,sourceCpp)
+useDynLib(shar, .registration = TRUE)
diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,24 @@
+# shar 0.3
+* Improvements
+  * `plot_randomized_pattern()` now uses envelopes to plot randomized summary functions
+  * `plot_randomized_pattern()` includes a quantum bar
+  * `plot_randomized_pattern()` now can return plots after each other (Press <Enter>)
+  * `calculate_engery()` can also calculate the energy for marked reconstructions
+  * Improved warning messages
+* Bugfixes
+  * Explicitly C++11 compiler
+* New functionality
+  * `plot_randomized_pattern()` to plot randomized rasters
+
+# shar 0.2.1
+* Improvements
+  * minor speed improvment for `reconstruct_pattern()`, `reconstruct_marks()` and `calculate_energy()`
+    * The starting pattern is now identical for all n_random and only the relocation process differs between randomizations
+    * All summary functions are only calculated for 250 steps from 0 to rmax
+* Bugfixes
+* New functionality
+  * `rcpp_sample()` as a faster Rcpp implementation of `sample()`
+
 # shar 0.2
 * Improvements
   * Replaced `cat()` with `message()` for all printing to console

diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -0,0 +1,25 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#' rcpp_sample
+#'
+#' @description Rcpp sample function
+#'
+#' @param x Vector of elements to sample from.
+#' @param n Size of the sample.
+#' @param replace Sample with replacement.
+#'
+#' @details
+#' \code{Rcpp} implementation of the \code{sample} function.
+#'
+#' @seealso
+#' \code{\link{sample}}
+#'
+#' @return vector
+#'
+#' @name rcpp_sample
+#' @export
+rcpp_sample <- function(x, n, replace = FALSE) {
+    .Call(`_shar_rcpp_sample`, x, n, replace)
+}
+
diff --git a/R/SHAR.R b/R/SHAR.R
@@ -1,3 +1,18 @@
+#' @title shar
+#'
+#' @description
+#' Analyse species-habitat associations in R. Therefore, information about the
+#' location of the species is needed and about the environmental conditions. To test
+#' for significance habitat associations, one of the two components is randomized.
+#' Methods are mainly based on Plotkin et al. (2000) <doi:10.1006/jtbi.2000.2158> and
+#' Harms et al. (2001) <doi:10.1111/j.1365-2745.2001.00615.x>.
+#'
+#' @name shar
+#' @docType package
+#' @useDynLib shar, .registration = TRUE
+#' @importFrom Rcpp sourceCpp
+"_PACKAGE"
+
 # Global variables
 globalVariables(c("count",
                   "habitat",
@@ -9,3 +24,6 @@ globalVariables(c("count",
                   "theo",
                   "type",
                   "x_r"))
+
+
+
diff --git a/R/calculate_energy.R b/R/calculate_energy.R
@@ -3,7 +3,8 @@
 #' @description Calculate mean energy
 #'
 #' @param pattern List with reconstructed patterns.
-#' @param return_mean Return the mean energy
+#' @param return_mean Return the mean energy.
+#' @param method String to specifiy if spatial pattern or marks were reconstructed.
 #' @param comp_fast If pattern contains more points than threshold, summary functions are estimated in a computational fast way.
 #' @param verbose Print progress report.
 #'
@@ -23,10 +24,16 @@
 #' @return numeric
 #'
 #' @examples
-#' pattern_random <- fit_point_process(species_a, n_random = 39)
+#' pattern_random <- fit_point_process(species_a, n_random = 19)
 #' calculate_energy(pattern_random)
 #' calculate_energy(pattern_random, return_mean = TRUE)
 #'
+#' \dontrun{
+#' marks_sub <- spatstat::subset.ppp(species_a, select = dbh)
+#' marks_recon <- reconstruct_marks(pattern_random[[1]], marks_sub, n_random = 19, max_runs = 1000)
+#' calculate_energy(marks_recon, return_mean = FALSE, method = "marks")
+#' }
+#'
 #' @aliases calculate_energy
 #' @rdname calculate_energy
 #'
@@ -38,81 +45,161 @@
 #' in ecology. Boca Raton: Chapman and Hall/CRC Press.
 
 #' @export
-calculate_energy <- function(pattern, return_mean = FALSE, comp_fast = 1000, verbose = TRUE){
+calculate_energy <- function(pattern,
+                             return_mean = FALSE,
+                             method = "spatial",
+                             comp_fast = 1000,
+                             verbose = TRUE){
 
   # check if randomized and observed is present
   if(!all(c(paste0("randomized_", seq_len(length(pattern) - 1)), "observed") == names(pattern)) || is.null(names(pattern))) {
     stop("Input must named 'randomized_1' to 'randomized_n' and includ 'observed' pattern.",
          call. = FALSE)
   }
 
-  # check if number of points exceed comp_fast limit
-  if(pattern$observed$n > comp_fast) {
-    comp_fast <- TRUE
-  }
-
-  else {
-    comp_fast <- FALSE
-  }
-
-  pattern_observed <- pattern[names(pattern) == "observed"] # extract observed pattern
-
-  pattern_reconstructed <- pattern[names(pattern) != "observed"] # extract randomized patterns
+  # extract observed pattern
+  pattern_observed <- pattern$observed
 
-  # calculate summary functions for observed pattern
-  if(comp_fast) {
+  # extract randomized patterns
+  pattern_reconstructed <- pattern[names(pattern) != "observed"]
 
-    gest_observed <- spatstat::Gest(X = pattern_observed[[1]], correction = "none")
+  # calculate r sequence
+  r <- seq(from = 0,
+           to = spatstat::rmax.rule(W = pattern_observed$window,
+                                    lambda = spatstat::intensity.ppp(pattern_observed)),
+           length.out = 250)
 
-    pcf_observed <- shar::estimate_pcf_fast(pattern = pattern_observed[[1]],
-                                            correction = "none",
-                                            method = "c",
-                                            spar = 0.5)
-  }
+  if (method == "spatial") {
 
-  else{
+    if(verbose) {
+      # check if pattern is marked
+      if(spatstat::is.marked(pattern_observed) || all(vapply(pattern_reconstructed,
+                                                             spatstat::is.marked,
+                                                             FUN.VALUE = logical(1)))) {
 
-    gest_observed <- spatstat::Gest(X = pattern_observed[[1]], correction = "han")
+        warning("Only energy of spatial summary functions are considered.", call. = FALSE)
+      }
+    }
 
-    pcf_observed <- spatstat::pcf(X = pattern_observed[[1]],
-                                  correction = "best", divisor = "d")
-  }
+    # check if number of points exceed comp_fast limit
+    if(pattern_observed$n > comp_fast) {
+      comp_fast <- TRUE
+    }
 
-  # loop through all reconstructed patterns
-  result <- vapply(seq_along(pattern_reconstructed), function(x) {
+    else {
+      comp_fast <- FALSE
+    }
 
-    # fast computation of summary stats
+    # calculate summary functions for observed pattern
     if(comp_fast) {
 
-      gest_reconstruction <- spatstat::Gest(X = pattern_reconstructed[[x]], correction = "none")
+      gest_observed <- spatstat::Gest(X = pattern_observed,
+                                      correction = "none",
+                                      r = r)
 
-      pcf_reconstruction <- shar::estimate_pcf_fast(pattern = pattern_reconstructed[[x]],
-                                                    correction = "none",
-                                                    method = "c",
-                                                    spar = 0.5)
+      pcf_observed <- shar::estimate_pcf_fast(pattern = pattern_observed,
+                                              correction = "none",
+                                              method = "c",
+                                              spar = 0.5,
+                                              r = r)
     }
 
-    # normal computation of summary stats
     else{
 
-      gest_reconstruction <- spatstat::Gest(X = pattern_reconstructed[[x]], correction = "han")
+      gest_observed <- spatstat::Gest(X = pattern_observed,
+                                      correction = "han",
+                                      r = r)
 
-      pcf_reconstruction <- spatstat::pcf(X = pattern_reconstructed[[x]],
-                                          correction = "best", divisor = "d")
+      pcf_observed <- spatstat::pcf(X = pattern_observed,
+                                    correction = "best",
+                                    divisor = "d",
+                                    r = r)
     }
 
-    # difference between observed and reconstructed pattern
-    energy <- mean(abs(gest_observed[[3]] - gest_reconstruction[[3]]), na.rm = TRUE) +
-      mean(abs(pcf_observed[[3]] - pcf_reconstruction[[3]]), na.rm = TRUE)
+    # loop through all reconstructed patterns
+    result <- vapply(seq_along(pattern_reconstructed), function(x) {
 
-    # print progress
-    if(verbose) {
-      message("\r> Progress: ", x, "/", length(pattern_reconstructed), appendLF = FALSE)
+      # fast computation of summary stats
+      if(comp_fast) {
+
+        gest_reconstruction <- spatstat::Gest(X = pattern_reconstructed[[x]],
+                                              correction = "none",
+                                              r = r)
+
+        pcf_reconstruction <- shar::estimate_pcf_fast(pattern = pattern_reconstructed[[x]],
+                                                      correction = "none",
+                                                      method = "c",
+                                                      spar = 0.5,
+                                                      r = r)
+      }
+
+      # normal computation of summary stats
+      else{
+
+        gest_reconstruction <- spatstat::Gest(X = pattern_reconstructed[[x]],
+                                              correction = "han",
+                                              r = r)
+
+        pcf_reconstruction <- spatstat::pcf(X = pattern_reconstructed[[x]],
+                                            correction = "best",
+                                            divisor = "d",
+                                            r = r)
+      }
+
+      # difference between observed and reconstructed pattern
+      energy <- mean(abs(gest_observed[[3]] - gest_reconstruction[[3]]), na.rm = TRUE) +
+        mean(abs(pcf_observed[[3]] - pcf_reconstruction[[3]]), na.rm = TRUE)
+
+      # print progress
+      if(verbose) {
+        message("\r> Progress: ", x, "/", length(pattern_reconstructed), appendLF = FALSE)
+      }
+
+      return(energy)
+
+    }, FUN.VALUE = numeric(1))
+  }
+
+  else if( method == "marks") {
+
+    # check if pattern is marked
+    if(!spatstat::is.marked(pattern_observed) || !all(vapply(pattern_reconstructed,
+                                                             spatstat::is.marked,
+                                                             FUN.VALUE = logical(1)))) {
+
+      stop("Please provide pattern with reconstruced marks.", call. = FALSE)
     }
 
-    return(energy)
+    # calculate summary functions
+    kmmr_observed <- spatstat::markcorr(pattern_observed,
+                                        correction = "Ripley",
+                                        r = r)
+
+    result <- vapply(seq_along(pattern_reconstructed), function(x) {
 
-  }, FUN.VALUE = numeric(1))
+      # calculate summary functions
+      kmmr_reconstruction <- spatstat::markcorr(pattern_reconstructed[[x]],
+                                                correction = "Ripley",
+                                                r = r)
+
+      # difference between observed and reconstructed pattern
+      energy <- mean(abs(kmmr_observed[[3]] - kmmr_reconstruction[[3]]), na.rm = TRUE) +
+        mean(abs(kmmr_observed[[3]] - kmmr_reconstruction[[3]]), na.rm = TRUE)
+
+      # print progress
+      if(verbose) {
+        message("\r> Progress: ", x, "/", length(pattern_reconstructed), appendLF = FALSE)
+      }
+
+      return(energy)
+
+    }, FUN.VALUE = numeric(1))
+  }
+
+  else {
+    stop("Please select either 'method = spatial' or 'method = marks'.",
+         call. = FALSE)
+  }
 
   # return mean for all reconstructed patterns
   if(return_mean) {

diff --git a/R/estimate_pcf_fast.R b/R/estimate_pcf_fast.R
@@ -35,7 +35,7 @@
 #' @export
 estimate_pcf_fast <- function(pattern, ...){
 
-  k_fun <- spatstat::Kest(X = pattern, ...) # estimate K-fct
+  k_fun <- suppressMessages(spatstat::Kest(X = pattern, ...)) # estimate K-fct
 
   result <- spatstat::pcf.fv(X = k_fun, ...) # estimate pcf from K-fct