Saving memory for TRI computations #69
Description
Hello Jeffrey and spatialEco contributors, first off thanks for this great package. I recently had to compute a few topographical metrics for a big collection of large rasters for a workflow and ran into some speed/memory usage problems. Instead of slicing everything into smaller chunks to complete the workflow, I dove into the code and found that if we use modified weights in a terra::focalMat to set the centre cell to 0, it allows downstream computations to strictly rely terra's C++ enabled functions. The end result result is identical to both exact = TRUE and FALSE in spatialEco::tri :
# Import libraries
library(bench)
#> Warning: package 'bench' was built under R version 4.5.2
library(terra)
#> Warning: package 'terra' was built under R version 4.5.1
#> terra 1.8.60
library(spatialEco)
#> Warning: package 'spatialEco' was built under R version 4.5.2
#>
#> Attaching package: 'spatialEco'
#> The following objects are masked from 'package:terra':
#>
#> is.empty, shift, sieve
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following object is masked from 'package:spatialEco':
#>
#> combine
#> The following objects are masked from 'package:terra':
#>
#> intersect, union
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(tidyr)
#>
#> Attaching package: 'tidyr'
#> The following object is masked from 'package:terra':
#>
#> extract
library(ggplot2)
#> Warning: package 'ggplot2' was built under R version 4.5.1
# Define the area of interest
aoi <- as.polygons(
x = ext(-1641555.5, -1633868.0, 340608.7, 348892.2),
crs = "EPSG:3979"
)
# Import NRCAN elevation dataset
dem <- rast(paste0("/vsicurl/", "https://canelevation-dem.s3.ca-central-1.amazonaws.com/mrdem-30/mrdem-30-dsm.tif")) %>%
crop(aoi, mask = TRUE)
# Fast TRI function
fast_tri <- function(x, dist, window = c("circle", "rectangle"), ...) {
# Check window type
if (!window %in% c("circle", "rectangle")) {
stop("`type` must be either 'circle' or 'rectangle'.")
}
# Define window, set weights to 1s and 0s
window <- terra::focalMat(x = x, d = dist, type = window)
window <- ifelse(window > 0, 1, 0)
# Set window centre cell to NA
window[((length(window) + 1)/2)] <- 0
# Sum of neighbours
sum_neigh <- terra::focal(x = x, w = window, fun = "sum", ...)
# Sum of squared neighbours
sum_sq_neigh <- terra::focal(x = x^2, w = window, fun = "sum", ...)
# Number of neighbours
n_neigh <- sum(window) # Even faster
# n_neigh <- terra::focal(
# x = terra::ifel(x > 0, 1, NA),
# w = window,
# fun = "sum",
# ...
# )
# TRI formula (expanded):
tri <- sqrt(sum_sq_neigh - 2 * x * sum_neigh + n_neigh * x^2)
# Return
return(tri)
}
# Compute fast TRI
tri_fast <- fast_tri(dem, dist = res(dem)[1], window = "rectangle")
tri_fast <- round(tri_fast, 10)
names(tri_fast) <- "tri"
# Compute spatialEco exact = TRUE
tri_exact <- tri(dem, s = 3, exact = TRUE)
tri_exact <- round(tri_exact, 10)
names(tri_exact) <- "tri"
# Compute spatialEco exact = FALSE
tri_algebra <- tri(dem, s = 3, exact = FALSE)
tri_algebra <- round(tri_algebra, 10)
names(tri_algebra) <- "tri"
# Check if identical up to 10 decimals (overkill)
identical(x = tri_exact, y = tri_fast)
#> [1] TRUE
identical(x = tri_exact, y = tri_algebra)
#> [1] TRUE
# Plot results side by side
par(mfrow = c(1, 3))
plot(tri_fast, main = "Fast TRI")
plot(tri_exact, main = "spatialEco (exact = TRUE)")
plot(tri_algebra, main = "spatialEco (exact = FALSE)")
# Test speed
benchmark <- bench::mark(
fast_tri(dem, dist = 30, window = "rectangle"),
tri(dem, s = 3, exact = FALSE),
tri(dem, s = 3, exact = TRUE),
check = FALSE,
iterations = 100
)
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
# Plot speed and memory
autoplot(benchmark, type = "boxplot")
ggplot(benchmark, aes(x = expression, y = mem_alloc)) +
geom_col() +
scale_y_bench_bytes() +
coord_flip() +
labs(x = "Expression", y = "Total Memory Allocated")
Created on 2025-12-15 with reprex v2.1.1
Session info
sessionInfo()
#> R version 4.5.0 (2025-04-11 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 22631)
#>
#> Matrix products: default
#> LAPACK version 3.12.1
#>
#> locale:
#> [1] LC_COLLATE=English_Canada.utf8 LC_CTYPE=English_Canada.utf8
#> [3] LC_MONETARY=English_Canada.utf8 LC_NUMERIC=C
#> [5] LC_TIME=English_Canada.utf8
#>
#> time zone: America/Toronto
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] ggplot2_3.5.2 tidyr_1.3.1 dplyr_1.1.4 spatialEco_2.0-3
#> [5] terra_1.8-60 bench_1.1.4
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 compiler_4.5.0 reprex_2.1.1 tidyselect_1.2.1
#> [5] Rcpp_1.0.14 xml2_1.3.8 scales_1.4.0 yaml_2.3.10
#> [9] fastmap_1.2.0 R6_2.6.1 generics_0.1.4 curl_6.3.0
#> [13] knitr_1.50 tibble_3.3.0 pillar_1.10.2 RColorBrewer_1.1-3
#> [17] rlang_1.1.6 xfun_0.52 fs_1.6.6 cli_3.6.5
#> [21] withr_3.0.2 magrittr_2.0.3 digest_0.6.37 grid_4.5.0
#> [25] rstudioapi_0.17.1 lifecycle_1.0.4 vctrs_0.6.5 evaluate_1.0.3
#> [29] glue_1.8.0 farver_2.1.2 codetools_0.2-20 profmem_0.7.0
#> [33] rmarkdown_2.29 purrr_1.0.4 tools_4.5.0 pkgconfig_2.0.3
#> [37] htmltools_0.5.8.1Perhaps this is something you've already thought a great deal about and I'm late to the party, but figured I'd contribute my experiment here in case it is useful in any way. An added bonus to the decrease in memory use is that the use of terra::focalMat internally in the function allows the user to easily specify circular windows, which are more intuitive to me when producing spatial indices. One thing I remain uncertain about is the use of the custom r.sqrt function to return NA in spatialEco:tri, perhaps there is no way around that for certain use cases, but I haven't yet encountered the need for it in my work so far.