Merge pull request #308 from r-spatialecology/commoncalcs

IIASA 2023 improvements
r-spatialecology · Nov 1, 2023 · b866cbb · b866cbb
2 parents 81501e5 + e266358
commit b866cbb
Show file tree

Hide file tree

Showing 193 changed files with 2,616 additions and 1,761 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: landscapemetrics
 Title: Landscape Metrics for Categorical Map Patterns
-Version: 2.0.0
+Version: 2.1.0
 Authors@R: c(person("Maximilian H.K.", "Hesselbarth",
              role = c("aut", "cre"),
              email = "[email protected]",

diff --git a/NAMESPACE b/NAMESPACE
@@ -15,11 +15,17 @@ export(construct_buffer)
 export(data_info)
 export(extract_lsm)
 export(get_adjacencies)
+export(get_area_patches)
 export(get_boundaries)
 export(get_centroids)
 export(get_circumscribingcircle)
+export(get_class_patches)
+export(get_complexity)
+export(get_enn_patch)
 export(get_nearestneighbour)
 export(get_patches)
+export(get_perimeter_patch)
+export(get_points)
 export(get_unique_values)
 export(landscape_as_list)
 export(list_lsm)
@@ -160,6 +166,7 @@ export(matrix_to_raster)
 export(options_landscapemetrics)
 export(pad_raster)
 export(points_as_mat)
+export(prepare_extras)
 export(proj_info)
 export(raster_to_points)
 export(rcpp_get_nearest_neighbor)

diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,30 @@
+# landscapemetrics 2.1.0
+* Improvements
+    * Many performance improvements. Most visible are in
+    `calculate_lsm` (all metrics: more than 5 times faster with 70% less memory 
+    allocation for `augusta_nlcd`; larger increases were found for smaller data)
+    and `window_lsm` (a single metric: more than 6 times faster for `augusta_nlcd`; 
+    larger increases were found for smaller data)
+    * Some performance improvements are related to the new "extras" mechanism, in which several objects 
+    are precalculated in `calculate_lsm`
+    * Creates an internal `extras_df` object that lists which extras are needed by
+    each metric
+    * Replaces the use of `tibble::tibble()` with `tibble::new_tibble(list())` in most functions.
+    This change is partially responsible for improvements of the `window_lsm` speed
+    * Replaces `raster_to_points` with `get_points` in several places. 
+    The `get_points` function is based on the column and row numbers multiplied by
+    the resolution, not actual coordinates.
+    * Replaces `table` with (faster) `tabulate` in `lsm_p_core`
+* New functions
+    * Adds a few internal helper functions and documents them, including `prepare_extras`,
+    `get_area_patches`, `get_class_patches`, `get_complexity`, `get_enn_patch`, 
+    `get_points`, and `get_perimeter_patch`
+* Bugfixes
+    * Fixes `window_lsm` behaviour for situations with NAs values and non-square windows
+* Various
+    * Fixes several typos and improves documentation in many places
+    * Uses object references in most rcpp functions
+
 # landscapemetrics 2.0.0
 * Improvements
     * `terra` and `sf` instead of `raster` and `sp` as underlying frameworks
@@ -11,7 +38,7 @@
     * The shape index now follows exactly the definition of the FRAGSTATS manual
     * Minor bug in clumpy index fixed
 * Various 
-    * Updated FRAGSTATS reference (thanks to  Oto Kaláb @kalab-oto)
+    * Updated FRAGSTATS reference (thanks to Oto Kaláb @kalab-oto)
     * Update FRAGSTATS tests
 
 # landscapemetrics 1.5.6

diff --git a/R/calculate_lsm.R b/R/calculate_lsm.R
@@ -151,6 +151,7 @@ calculate_lsm_internal <- function(landscape,
                     call. = FALSE)
         }
     }
+    landscape <- terra::as.int(landscape)
 
     # get name of metrics
     metrics <- list_lsm(level = level, metric = metric, name = name,
@@ -162,18 +163,14 @@ calculate_lsm_internal <- function(landscape,
     # how many metrics need to be calculated?
     number_metrics <- length(metrics_calc)
 
-    # get coordinates of cells
-    points <- raster_to_points(landscape)[, 2:4]
-
-    # resolution of original raster
+    # prepare extras
     resolution <- terra::res(landscape)
-
-    # convert to matrix
     landscape <- terra::as.matrix(landscape, wide = TRUE)
+    extras <- prepare_extras(metrics, landscape, directions, neighbourhood,
+                                        ordered, base, resolution)
 
     result <- do.call(rbind, lapply(seq_along(metrics_calc), FUN = function(current_metric) {
-
-        # print progess using the non-internal name
+        # print progress using the non-internal name
         if (progress) {
             cat("\r> Progress metrics: ", current_metric, "/", number_metrics)
         }
@@ -185,25 +182,31 @@ calculate_lsm_internal <- function(landscape,
         arguments <- names(formals(foo))
 
         # run function
-        tryCatch(do.call(what = foo,
+        #start_time = Sys.time()
+        resultint <- tryCatch(do.call(what = foo,
                          args = mget(arguments, envir = parent.env(environment()))),
                  error = function(e){
                      message("")
                      stop(e)})
+
+        #end_time = Sys.time()
+        #resultint$time <- as.numeric(difftime(end_time, start_time, units = "secs"))
+        resultint
         })
     )
 
     if (full_name == TRUE) {
 
         col_ordering <- c("level", "class", "id", "metric", "value",
-                          "name", "type", "function_name")
+                          "name", "type", "function_name"#,"time"
+                          )
 
         result <- merge(x = result,
                         y = lsm_abbreviations_names,
                         by = c("level", "metric"),
                         all.x = TRUE, sort = FALSE, suffixes = c("", ""))
 
-        result <- tibble::as_tibble(result[,col_ordering])
+        result <- tibble::as_tibble(result[, col_ordering])
     }
 
     if (progress) {

diff --git a/R/construct_buffer.R b/R/construct_buffer.R
@@ -2,15 +2,15 @@
 #'
 #' @description Internal function to construct plot area around coordinates
 #'
-#' @param coords SpatialPoints or 2-column matrix with coordinates of sample points
+#' @param coords SpatVector, sf object or 2-column matrix with coordinates of sample points
 #' @param shape String specifying plot shape. Either "circle" or "square"
 #' @param size Size of sample plot. Equals the radius for circles or the
-#' side-length for squares in mapunits
-#' @param return_vec If true, vector objects are returned.
+#' side-length for squares in map units
+#' @param return_vec If TRUE, vector objects are returned.
 #' @param verbose Print warning messages.
 #'
 #' @return
-#' matrix or sf objecct
+#' matrix or SpatVector object
 #'
 #' @examples
 #' coords <- matrix(c(10, 5, 25, 15, 5, 25), ncol = 2, byrow = TRUE)

diff --git a/R/data.R b/R/data.R
@@ -26,7 +26,6 @@
 #' @source http://maps.elie.ucl.ac.be/CCI/viewer/
 "podlasie_ccilc"
 
-
 #' Tibble of abbreviations coming from FRAGSTATS
 #'
 #' A single tibble for every abbreviation of every metric that is

diff --git a/R/data_info.R b/R/data_info.R
@@ -39,5 +39,5 @@ data_info <- function(landscape){
                         yes = "integer",
                         no = "non-integer"))
 
-    tibble::tibble(class = class, n_classes = length(landscape_values))
+    tibble::new_tibble(list(class = class, n_classes = length(landscape_values)))
 }
diff --git a/R/get_boundaries.R b/R/get_boundaries.R
@@ -17,7 +17,7 @@
 #' cell or a cell with a different value than itself. Non-boundary cells only
 #' neighbour cells with the same value than themself.
 #'
-#' @return List with RasterLayer or matrix
+#' @return List with SpatRaster or matrix
 #'
 #' @examples
 #' landscape <- terra::rast(landscapemetrics::landscape)

diff --git a/R/get_centroids.R b/R/get_centroids.R
@@ -51,7 +51,7 @@ get_centroids <- function(landscape, directions = 8, cell_center = FALSE,
 
     if (return_vec) {
 
-        result <- terra::vect(result, geom=c("x", "y"), crs = crs)
+        result <- terra::vect(result, geom = c("x", "y"), crs = crs)
     }
 
     return(result)
@@ -73,14 +73,14 @@ get_centroids_calc <- function(landscape, directions, cell_center, verbose) {
     # all values NA
     if (all(is.na(landscape))) {
 
-        return(tibble::tibble(level = "patch",
+        return(tibble::new_tibble(list(level = "patch"),
                               class = as.integer(NA),
                               id = as.integer(NA),
                               y = as.double(NA),
                               y = as.double(NA)))
     }
 
-    # get uniuqe class id
+    # get unique class id
     classes <- get_unique_values_int(landscape, verbose = verbose)
 
     centroid <- do.call(rbind,
@@ -100,11 +100,11 @@ get_centroids_calc <- function(landscape, directions, cell_center, verbose) {
                             # set ID from class ID to unique patch ID
                             points[, 3] <- landscape_labeled[!is.na(landscape_labeled)]
 
-                            # # conver to tibble
+                            # # convert to tibble
                             points <- stats::setNames(object = data.frame(points),
                                                       nm = c("x", "y", "id"))
 
-                            # calcuale the centroid of each patch (mean of all coords)
+                            # calculate the centroid of each patch (mean of all coords)
                             centroid_temp <- stats::aggregate(points[, c(1, 2)],
                                                               by = list(id = points[, 3]),
                                                               FUN = mean)
@@ -159,9 +159,9 @@ get_centroids_calc <- function(landscape, directions, cell_center, verbose) {
         }
     }
 
-    tibble::tibble(level = "patch",
+    tibble::new_tibble(list(level = rep("patch", nrow(centroid)),
                    class = as.integer(centroid$class),
-                   id = as.integer(id),
+                   id = as.integer(centroid$id),
                    x = as.double(centroid$x),
-                   y = as.double(centroid$y))
+                   y = as.double(centroid$y)))
 }
diff --git a/R/get_circumscribingcircle.R b/R/get_circumscribingcircle.R
@@ -100,12 +100,12 @@ get_circumscribingcircle_calc <- function(landscape, level, directions) {
                           )
 
         # resulting tibble
-        circle <- tibble::tibble(level = "patch",
+        circle <- tibble::new_tibble(list(level = rep("patch", nrow(circle)),
                                  class = as.integer(circle$class),
                                  id = as.integer(seq_len(nrow(circle))),
                                  value = circle$circle_diameter,
                                  center_x = circle$circle_center_x,
-                                 center_y = circle$circle_center_y)
+                                 center_y = circle$circle_center_y))
     }
 
     # class level (no labeling)
@@ -115,12 +115,12 @@ get_circumscribingcircle_calc <- function(landscape, level, directions) {
         circle_class <- rcpp_get_circle(landscape, resolution_xy = resolution[1])
 
         # resulting tibble
-        circle <- tibble::tibble(level = "class",
+        circle <- tibble::new_tibble(list(level = rep("class", nrow(circle_class)), 
                                  class = as.integer(circle_class$patch_id),
-                                 id = as.integer(NA),
+                                 id = rep(as.integer(NA), nrow(circle_class)),
                                  value = circle_class$circle_diameter,
                                  center_x = circle_class$circle_center_x,
-                                 center_y = circle_class$circle_center_y)
+                                 center_y = circle_class$circle_center_y))
     }
 
     # shift the coordinates to the original coordinate system

diff --git a/R/get_nearestneighbour.R b/R/get_nearestneighbour.R
@@ -46,14 +46,16 @@ get_nearestneighbour <- function(landscape, return_id = FALSE) {
 
 }
 
-get_nearestneighbour_calc <- function(landscape, return_id,
+get_nearestneighbour_calc <- function(landscape, return_id, resolution,
                                       points = NULL) {
 
+    if (missing(resolution)) resolution <- terra::res(landscape)
+
     # convert to matrix
     if (!inherits(x = landscape, what = "matrix")) {
 
         # get coordinates and values of all cells
-        points <- raster_to_points(landscape)[, 2:4]
+        points <- get_points(landscape, resolution = resolution)
 
         # convert to matrix
         landscape <- terra::as.matrix(landscape, wide = TRUE)
@@ -79,12 +81,12 @@ get_nearestneighbour_calc <- function(landscape, return_id,
     num <- seq_along(ord)
     rank <- match(num, ord)
 
-    res <- rcpp_get_nearest_neighbor(terra::as.matrix(points, wide= TRUE)[ord, ])
+    res <- rcpp_get_nearest_neighbor(as.matrix(points)[ord, ])
 
-    min_dist <- tibble::tibble(cell = num,
+    min_dist <- tibble::new_tibble(list(cell = num,
                                dist = res[rank, 1],
                                id_focal = points[, 3],
-                               id_neighbour = res[rank, 2])
+                               id_neighbour = res[rank, 2]))
 
     min_dist_aggr <- stats::setNames(stats::aggregate(x = min_dist$dist,
                                                       by = list(min_dist$id_focal),

diff --git a/R/get_patches.R b/R/get_patches.R
@@ -31,7 +31,7 @@
 #' algorithm based on immersion simulations. IEEE Transactions on Pattern
 #' Analysis and Machine Intelligence. 13 (6), 583-598
 #'
-#' @return List
+#' @return List of SpatRaster
 #'
 #' @examples
 #' landscape <- terra::rast(landscapemetrics::landscape)

diff --git a/R/get_unique_values.R b/R/get_unique_values.R
@@ -56,7 +56,7 @@ get_unique_values <- function(x, simplify = FALSE, verbose = TRUE) {
     return(result)
 }
 
-get_unique_values_int <- function(landscape, verbose) {
+get_unique_values_int <- function(landscape, verbose = FALSE) {
 
     if (inherits(x = landscape, what = "SpatRaster")) {
 

diff --git a/R/landscape_as_list.R b/R/landscape_as_list.R
@@ -21,7 +21,6 @@ landscape_as_list <- function(landscape) UseMethod("landscape_as_list")
 #' @name landscape_as_list
 #' @export
 landscape_as_list.SpatRaster <- function(landscape) {
-
     landscape <- terra::as.list(landscape)
 
     return(landscape)

diff --git a/R/landscapemetrics-package.R b/R/landscapemetrics-package.R
@@ -15,4 +15,4 @@
 #' @keywords internal
 "_PACKAGE"
 
-globalVariables(c("label", "lsm_abbreviations_names", "metric_1", "metric_2", "value", "values", "x", "y"))
+globalVariables(c(".data", "label", "lsm_abbreviations_names", "metric_1", "metric_2", "value", "values", "x", "y"))
diff --git a/R/lsm_c_ai.R b/R/lsm_c_ai.R
@@ -54,28 +54,29 @@ lsm_c_ai <- function(landscape) {
     tibble::add_column(result, layer, .before = TRUE)
 }
 
-lsm_c_ai_calc <- function(landscape) {
+lsm_c_ai_calc <- function(landscape, extras = NULL) {
 
-    # convert to raster to matrix
-    if (!inherits(x = landscape, what = "matrix")) {
+    if (is.null(extras)){
+        metrics <- "lsm_c_ai"
         landscape <- terra::as.matrix(landscape, wide = TRUE)
+        extras <- prepare_extras(metrics, landscape_mat = landscape)
     }
 
     # all values NA
     if (all(is.na(landscape))) {
-        return(tibble::tibble(level = "class",
+        return(tibble::new_tibble(list(level = "class",
                               class = as.integer(NA),
                               id = as.integer(NA),
                               metric = "ai",
-                              value = as.double(NA)))
+                              value = as.double(NA))))
     }
 
     # get coocurrence matrix of like_adjacencies
     like_adjacencies <- rcpp_get_coocurrence_matrix_diag(landscape,
                                                          directions = as.matrix(4)) / 2
 
     # get number of cells each class
-    cells_class <- rcpp_get_composition_vector(landscape)
+    cells_class <- extras$composition_vector
 
     # calculate maximum adjacencies
     n <- trunc(sqrt(cells_class))
@@ -96,9 +97,9 @@ lsm_c_ai_calc <- function(landscape) {
     # max_adj can be zero if only one cell is present; set to NA
     ai[is.nan(ai)] <- NA
 
-    return(tibble::tibble(level = "class",
+    return(tibble::new_tibble(list(level = rep("class", length(ai)),
                           class = as.integer(names(like_adjacencies)),
-                          id = as.integer(NA),
-                          metric = "ai",
-                          value = as.double(ai)))
+                          id = rep(as.integer(NA), length(ai)),
+                          metric = rep("ai", length(ai)),
+                          value = as.double(ai))))
 }