r-spatialecology / landscapemetrics

Landscape Metrics for Categorical Map Patterns 🗺️ in R
https://r-spatialecology.github.io/landscapemetrics
GNU General Public License v3.0
227 stars 43 forks source link

Potential spatialize_lsm (direction = 4) bug #317

Closed MStuhlmacher closed 4 months ago

MStuhlmacher commented 4 months ago

First off, thank you for a great package! It's so nice to use compared to the FRAGSTATS GUI and I really appreciate it. I might have found a bug while running the "spatialize_lsm" command on some proprietary data. The spatialize_lsm output raster doesn't match the tabular output from "calculate_lsm" but only when I use the "directions = 4", the "directions = 8" did not have this problem. Below is a re-producible example of this with a toy raster:

library(raster)
library(landscapemetrics)

#Create a 10 by 10 raster in which all cells are 1
r = raster(ncol=10, nrow=10)
values(r) = 1

#Randomly assign some of the values to NA
set.seed(123)  # For reproducibility
na_matrix = matrix(runif(100) < 0.8, nrow=10, ncol=10)

# Apply NA values to the raster based on the NA matrix
values(r)[na_matrix] = NA

#Assign crs
crs(r) = "EPSG:32633"
plot(r)

#Run the tabular
check_landscape(r)
tabular = calculate_lsm(r, what = c('lsm_p_area'),verbose=T,directions = 4)

#Run the spatialize
spatialize = spatialize_lsm(r,what = "lsm_p_area", directions = 4)
spatializeR = spatialize$layer_1$lsm_p_area #Extract the raster
plot(spatializeR) 

Here's the resulting plot:

image

The circled green and grey pixel values look like they are wrong. Based on the tabular results they should be 0.06 and 0.1944 respectively (and the ones that are yellow should be 0.1296). I exported this raster to check that it's not a visualization issue as well. This roughly mirrors what I was running into with the proprietary data as well. Let me know if it would be useful to provide more info on that.

Here is contents of utils::sessionInfo() for my machine: R version 4.2.1 (2022-06-23 ucrt) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale: [1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8 LC_MONETARY=English_United States.utf8 [4] LC_NUMERIC=C LC_TIME=English_United States.utf8

attached base packages: [1] stats graphics grDevices utils datasets methods base

other attached packages: [1] landscapemetrics_2.1.1 raster_3.6-26 sp_2.1-2

loaded via a namespace (and not attached): [1] Rcpp_1.0.12 magrittr_2.0.3 units_0.8-5 tidyselect_1.2.0 lattice_0.20-45 R6_2.5.1 rlang_1.1.3 fansi_1.0.6
[9] dplyr_1.1.4 tools_4.2.1 grid_4.2.1 KernSmooth_2.23-20 utf8_1.2.4 terra_1.7-65 cli_3.6.2 e1071_1.7-14
[17] DBI_1.2.1 class_7.3-20 tibble_3.2.1 lifecycle_1.0.4 sf_1.0-15 vctrs_0.6.5 codetools_0.2-18 glue_1.7.0
[25] proxy_0.4-27 compiler_4.2.1 pillar_1.9.0 generics_0.1.3 classInt_0.4-10 pkgconfig_2.0.3

I appreciate you looking into this! Hopefully I'm not missing something super obvious.

mhesselbarth commented 4 months ago

Hey,

Thanks for finding this. I think the argument does not get passed on properly. I will try to fix this

mhesselbarth commented 4 months ago

Okay, should be fixed in c5336d1e

library(terra)
library(landscapemetrics)

#Create a 10 by 10 raster in which all cells are 1
r = rast(ncol=10, nrow=10)
values(r) = 1

#Randomly assign some of the values to NA
set.seed(123)  # For reproducibility
na_matrix = matrix(runif(100) < 0.8, nrow=10, ncol=10)

# Apply NA values to the raster based on the NA matrix
values(r)[na_matrix] = NA

#Assign crs
crs(r) = "EPSG:32633"

# show landscape patch ids
show_patches(landscape = r, class = 1, directions = 4, labels = TRUE)
#> $layer_1


# calculate the area
area_df <- lsm_p_area(landscape = r, directions = 4)

# get raster layer with values
area_ras <- spatialize_lsm(r, what = "lsm_p_area", directions = 4)

# those should be identical and four unique values
unique(area_df$value)
#> [1] 0.2592 0.1296 0.0648 0.1944
unique(area_ras$layer_1$lsm_p_area[])
#>       value
#> [1,]     NA
#> [2,] 0.1296
#> [3,] 0.0648
#> [4,] 0.2592
#> [5,] 0.1944

# this seems to be correct now
plot(area_ras$layer_1$lsm_p_area, col = viridis::viridis(n = 4, option = "E"))


# compare to function by the package
show_lsm(landscape = r, what = "lsm_p_area", class = 1, directions = 4, labels = TRUE)
#> $layer_1

Created on 2024-02-29 with reprex v2.1.0

Could you please try to install the Github version and see if results are correct now.

MStuhlmacher commented 4 months ago

Great, thank you! It's on my list to try the updated version and report back

MStuhlmacher commented 4 months ago

I installed the Github version and re-ran my code from above and it now works as expected. Thanks for solving that so quickly!

mhesselbarth commented 4 months ago

Cool 👍 Will close this, please re-open if needed