Patch labeling inconsistencies between version 1.3 and 1.5.7
Closed this issue · 4 comments
Dear landscapemetrics package authors,
I am experiencing a problem with the way patches are being labeled by the patches() function and the way the calculate_lsm() function labels them in landscapemetrics version 1.5.7. ( I also tried with 1.5.4 and 1.5.3)
The problem arises when I try to use the patch id values produced by lsm_p_area() to reclassify pixels in my landscape data. The patch id values do not seem to correspond to the pixel values in my landscape data, which results in the wrong pixels being processed.
In version 1.3, I could use get_patches(landscape, directions=4) to derive patches, calculate metrics, apply filters, and then reclassify the patches satisfying certain conditions to a new class of interest without issues.
However, in version 1.5.7, the pixel values from the output of lsm_p_area() do not match the ID values, leading to incorrect pixel processing.
Here's a reproducible example using version 1.3 which produces expected results:
#Code using version 1.3
library(terra)
#> terra 1.7.39
suppressPackageStartupMessages(library(raster))
packageVersion("raster")
#> [1] '3.6.23'
library(landscapemetrics)
packageVersion("landscapemetrics")
#> [1] '1.3'
suppressPackageStartupMessages(library(dplyr))
#create test SpatRaster
set.seed(123)
ext <- ext(400304.45, 403304.45, 928602.05, 931602.05)
r <- rast(ext, ncols=100, nrows=100, res=30, vals=sample(1:2, 100*100, replace=TRUE))
(crs(r) <- "EPSG:3118")
#> [1] "EPSG:3118"
#select class 2 and convert into RasterLayer object
cls1 <- raster(segregate(r,classes=2, other=NA))
#extract patches and area metrics
met <- lsm_p_area(cls1, directions=4)
#display metrics
print(met, n=5)
#> # A tibble: 720 × 6
#> layer level class id metric value
#> <int> <chr> <int> <int> <chr> <dbl>
#> 1 1 patch 1 1 area 0.09
#> 2 1 patch 1 2 area 0.9
#> 3 1 patch 1 3 area 3.87
#> 4 1 patch 1 4 area 0.09
#> 5 1 patch 1 5 area 0.18
#> # ℹ 715 more rows
#Extract Patches
pat <- get_patches(cls1,directions=4)[[1]]
#get frequencies
freqs <- as_tibble(freq(pat))
freqs %>% arrange(desc(count)) %>% slice_head(n = 5)
#> # A tibble: 5 × 2
#> value count
#> <dbl> <dbl>
#> 1 NA 5017
#> 2 246 191
#> 3 350 156
#> 4 302 142
#> 5 486 122
# number of patches extracted:
nrow(freqs)
#> [1] 721
#filter by value
met_filt <- met %>% filter(value>10)
met_filt %>% arrange(desc(value))
#> # A tibble: 5 × 6
#> layer level class id metric value
#> <int> <chr> <int> <int> <chr> <dbl>
#> 1 1 patch 1 246 area 17.2
#> 2 1 patch 1 350 area 14.0
#> 3 1 patch 1 302 area 12.8
#> 4 1 patch 1 486 area 11.0
#> 5 1 patch 1 77 area 10.1
#get mask with the patches with area > 10ha
mskd <- subs(pat[[1]], met_filt[,c(4, 1)], by=1, which=2, subsWithNA=TRUE)
freq(mskd)
#> value count
#> [1,] 1 723
#> [2,] NA 9277However, using version 1.5.7, this issue arises:
#Code using version 1.5.7
library(terra)
#> terra 1.7.39
suppressPackageStartupMessages(library(raster))
packageVersion("raster")
#> [1] '3.6.23'
library(landscapemetrics)
#> Starting from v2.0.0, landscapemetrics does not support the 'raster' or 'sp' packages.
#> They are replaced by 'terra' and 'sf', respectively. More information
#> about the 'terra' package can be found here: https://rspatial.org/index.html.
packageVersion("landscapemetrics")
#> [1] '1.5.7'
suppressPackageStartupMessages(library(dplyr))
#create test SpatRaster
set.seed(123)
ext <- ext(400304.45, 403304.45, 928602.05, 931602.05)
r <- rast(ext, ncols=100, nrows=100, res=30, vals=sample(1:2, 100*100, replace=TRUE))
(crs(r) <- "EPSG:3118")
#> [1] "EPSG:3118"
#select class 2 and convert into RasterLayer object
cls1 <- segregate(r,classes=2, other=NA)
#extract area metrics
met <- lsm_p_area(cls1, directions=4)
#display metrics
print(met, n=5)
#> # A tibble: 720 × 6
#> layer level class id metric value
#> <int> <chr> <int> <int> <chr> <dbl>
#> 1 1 patch 1 1 area 0.63
#> 2 1 patch 1 2 area 0.72
#> 3 1 patch 1 3 area 0.36
#> 4 1 patch 1 4 area 0.09
#> 5 1 patch 1 5 area 2.7
#> # ℹ 715 more rows
#Extract patches
pat <- patches(cls1, directions=4, zeroAsNA=TRUE, allowGaps=FALSE)
#get frequencies
freqs <- as_tibble(freq(pat))
freqs %>% arrange(desc(count)) %>% slice_head(n = 5)
#> # A tibble: 5 × 3
#> layer value count
#> <dbl> <dbl> <dbl>
#> 1 1 246 191
#> 2 1 350 156
#> 3 1 302 142
#> 4 1 486 122
#> 5 1 77 112
# number of patches extracted:
nrow(freqs)
#> [1] 720
#filter by value
met_filt <- met %>% filter(value>10)
met_filt %>% arrange(desc(value))
#> # A tibble: 5 × 6
#> layer level class id metric value
#> <int> <chr> <int> <int> <chr> <dbl>
#> 1 1 patch 1 9 area 17.2
#> 2 1 patch 1 331 area 14.0
#> 3 1 patch 1 575 area 12.8
#> 4 1 patch 1 588 area 11.0
#> 5 1 patch 1 114 area 10.1
mskd <- subst(pat, from= met_filt$id, to=1, others=NA)
freq(mskd)
#> layer value count
#> 1 1 1 44Although the metric values seem correct, the ID cannot be used to target the pixels for reclassification, as was possible in the older version. While the terra package has functions like sieve for removing small patches, it doesn't directly allow identifying the pixels or setting intervals, and it doesn't work for other metrics either.
It seems that the issue is present across different environments. For reference, here are the sessionInfo outputs from two different environments where I experienced the problem:
MacOS Environment
sessionInfo()
#> R version 4.2.2 (2022-10-31)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Big Sur ... 10.16
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> loaded via a namespace (and not attached):
#> [1] digest_0.6.33 withr_2.5.0 lifecycle_1.0.3 reprex_2.0.2
#> [5] evaluate_0.21 rlang_1.1.1 cli_3.6.1 rstudioapi_0.15.0
#> [9] fs_1.6.3 rmarkdown_2.23 tools_4.2.2 glue_1.6.2
#> [13] xfun_0.39 yaml_2.3.7 fastmap_1.1.1 compiler_4.2.2
#> [17] htmltools_0.5.5 knitr_1.43- Linux Server Environment
``` r
sessionInfo()
#>R version 4.3.0 (2023-04-21)
#>Platform: x86_64-pc-linux-gnu (64-bit)
#>Running under: Ubuntu 22.04.2 LTS
#>
#>Matrix products: default
#>BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
#>LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#>
#>locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#>[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#>[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#>[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#>[9] LC_ADDRESS=C LC_TELEPHONE=C
#>[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C I appreciate your attention to this matter and look forward to any insights you might have.
Best regards
Jerónimo rodriguez-Escobar
Now, I was able to get the expected results with a slightly modified code:
library(terra)
library(raster)
library(landscapemetrics)
library(dplyr)
#create test SpatRaster
set.seed(123)
ext <- ext(400304.45, 403304.45, 928602.05, 931602.05)
r <- rast(ext, ncols=100, nrows=100, res=30, vals=sample(1:2, 100*100, replace=TRUE))
crs(r) <- "EPSG:3118"
cls1 <- raster(segregate(r,classes=2, other=NA))
met <- lsm_p_area(cls1, directions=4)
met_filt <- met |> filter(value>10)
#Extract Patches
pat <- raster(get_patches(cls1,directions=4)[[1]][[1]])
r <- subs(pat, met_filt[, c(4, 1)], by=1, which=2, subsWithNA=TRUE)
freq(r)
#> value count
#> [1,] 1 723
#> [2,] NA 9277Created on 2023-07-24 with reprex v2.0.2
Here is also terra-only version (raster package is not used):
library(terra)
library(landscapemetrics)
library(dplyr)
#create test SpatRaster
set.seed(123)
ext <- ext(400304.45, 403304.45, 928602.05, 931602.05)
r <- rast(ext, ncols=100, nrows=100, res=30, vals=sample(1:2, 100*100, replace=TRUE))
crs(r) <- "EPSG:3118"
cls1 <- segregate(r,classes=2, other=NA)
met <- lsm_p_area(cls1, directions=4)
met_filt <- met |> filter(value>10)
#Extract Patches
pat <- get_patches(cls1,directions=4)[[1]][[1]]
newr <- subst(pat, met_filt[[4]], met_filt[[1]], others = NA)
freq(newr)
#> layer value count
#> 1 1 1 723Thank you! I was aware that from ver. 1.4 the way patches were labeled had changed, but I was not able to get the id´s to match. So if I am correct the difference would be this line:
pat <-get_patches(cls1,directions=4)[[1]][[1]]and substitute patches() with get_patches() . I will try it right away.
Nevertheless, I see that spatialize_lsm() will suit my needs better.