search

r-spatialecology / landscapemetrics

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

shape index unexpected behaviour #300

Closed filippo-ferrario closed 1 year ago

filippo-ferrario commented 1 year ago

I think there might be problems with the metric lsm_l_shape_mn using a moving window

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.
# Setup data

Import raster data from Terra package.

landscape <- terra::rast(landscapemetrics::landscape)
terra::plot(landscape)

# polygonized for later plotting
pol_r1<-terra::as.polygons(landscape) 
terra::plot(subset(pol_r1, pol_r1$clumps==1), add=T, border='red', lwd=2)
terra::plot(subset(pol_r1, pol_r1$clumps==2), add=T, border='purple', lwd=2)

# terra::plot(subset(pol_r1, pol_r1$clumps==3), add=T, border='black', lwd=2)

As I need to evaluate the function on just a patch type in a moving window I first mask the other patch clasees.

# select just one class and replace the False with NA (this mimiks the behaviour of raster::rasterize used in  which have a default background to NA, as it has been done in `preprocessing-hyperframe_for_modeling.R`) 
class1<-landscape==1
class1[class1==0]<-NA
terra::plot(class1, col='black')

The function will be tested in both a rectangular and a circle moving window

# Create Rectangular moving window of side length 3
mw<-terra::focalMat(landscape, d=1, type=c('rectangle'), fillNA=T)
mw<-ifelse(mw>0, 1,NA)
mw
#>      [,1] [,2] [,3]
#> [1,]    1    1    1
#> [2,]    1    1    1
#> [3,]    1    1    1

Create Circular moving window of radius 1 (it looks like a cross but it is easier to picture it on the raster)

mwc<-terra::focalMat(landscape, d=1, type=c('circle'), fillNA=T)
mwc<-ifelse(mwc>0, 1,NA)
mwc
#>      [,1] [,2] [,3]
#> [1,]   NA    1   NA
#> [2,]    1    1    1
#> [3,]   NA    1   NA

Rectangular case for Shape index

sh<-window_lsm(class1,what='lsm_l_shape_mn', window=mw)
sh$layer_1$lsm_l_shape_mn
#> class       : SpatRaster 
#> dimensions  : 30, 30, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 30, 0, 30  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source(s)   : memory
#> name        :    clumps 
#> min value   : 0.8333333 
#> max value   : 1.5000000
# plot the patch and values
terra::plot(class1)
terra::plot(subset(pol_r1, pol_r1$clumps==1), add=T, border='red', lwd=2)
terra::text(sh$layer_1$lsm_l_shape_mn, digits=3, cex=0.5)

# highlight region of interest
v1<-terra::vect(rbind(c(18,8),c(21,8),c(21,11),c(18,11)),type='polygons')
v2<-terra::vect(rbind(c(16,19),c(19,19),c(19,22),c(16,22)),type='polygons')
v3<-terra::vect(rbind(c(8,9),c(11,9),c(11,12),c(8,12)),type='polygons')

terra::plot(v1, add=TRUE,border='black', lwd=3)
terra::plot(v2, add=TRUE,border='blue', lwd=3)
terra::plot(v3, add=TRUE, lwd=3, border='green')

Check behaviour

From the help page and the Patton 1975 reference I understad that:

  1. The hypothetical minimum perimeter should be that of a circle in Patton 1975, but I figure that landscapemetric uses a square, can you confirm?
    • looking at the green square the index is 1:
    • consistent with minimum square patch perimeter/minimum square perim = (3 * 4)/(sqrt(9) * 4) =1
    • not consistent with minimum circle: patch perimeter/minimum circonference = (3 * 4)/sqrt(4 * pi*9) =1.1283792
  2. looking at the blue square: is the value 1 correct, or should I read a 0 as per the help page?
    • From the help
      > Equals SHAPE_MN = 0 if all patches are squares
      I don’t know why it should be as described in the help, but anyhow it is not what the function is doing.
  3. looking at the black square I am expecting a value of patch perimeter/minimum square perim = 8/(sqrt(3)*4) =1.1547005 but the function return 1. Why?

Circular case for Shape index

csh<-window_lsm(class1,what='lsm_l_shape_mn', window=mwc)
#> Warning: number of items to replace is not a multiple of replacement length
csh$layer_1$lsm_l_shape_mn
#> class       : SpatRaster 
#> dimensions  : 30, 30, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 30, 0, 30  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source(s)   : memory
#> name        : clumps 
#> min value   :      1 
#> max value   :      2
# plot the patch and values
terra::plot(class1)
terra::plot(subset(pol_r1, pol_r1$clumps==1), add=T, border='red', lwd=2)
terra::text(csh$layer_1$lsm_l_shape_mn, digits=3, cex=0.5)

terra::plot(v1, add=TRUE,border='black', lwd=3)
terra::plot(v2, add=TRUE,border='blue', lwd=3)
terra::plot(v3, add=TRUE, lwd=3, border='green')

Check behaviour

Here I plot the region of interest as squares just as a reference, but the moving window should be a cross inscribed in the colored square.

  1. in the green window
    • I am expecting a value of 12/(sqrt(5)*4)=1.3416408
    • the value returned is 1.2, why?
  2. In the blue square I expect NA, instead there is a 1.
  3. the black seems correct.
  4. the pixel at x=1 y=20 should be 10/(sqrt(4)*4)=1.25 but it is 1.33
    • is it due to some edge correction?

Created on 2023-09-21 with reprex v2.0.2

Session info ``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.2.2 (2022-10-31 ucrt) #> os Windows 10 x64 (build 19044) #> system x86_64, mingw32 #> ui RTerm #> language ENG #> collate English_United States.utf8 #> ctype English_United States.utf8 #> tz America/New_York #> date 2023-09-21 #> pandoc 2.18 @ C:/Users/FERRAR~1/AppData/Local/Pandoc/ (via rmarkdown) #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> ! package * version date (UTC) lib source #> cli 3.6.1 2023-03-23 [1] CRAN (R 4.2.3) #> P codetools 0.2-18 2020-11-04 [4] CRAN (R 4.2.2) #> P curl 4.3.3 2022-10-06 [4] CRAN (R 4.2.2) #> P digest 0.6.31 2022-12-11 [4] CRAN (R 4.2.2) #> P evaluate 0.19 2022-12-13 [4] CRAN (R 4.2.2) #> fansi 1.0.4 2023-01-22 [1] CRAN (R 4.2.3) #> P fastmap 1.1.0 2021-01-25 [4] CRAN (R 4.2.2) #> P fs 1.5.2 2021-12-08 [4] CRAN (R 4.2.2) #> P glue 1.6.2 2022-02-24 [4] CRAN (R 4.2.2) #> P highr 0.9 2021-04-16 [4] CRAN (R 4.2.2) #> P htmltools 0.5.4 2022-12-07 [4] CRAN (R 4.2.2) #> P httr 1.4.4 2022-08-17 [4] CRAN (R 4.2.2) #> P knitr 1.41 2022-11-18 [4] CRAN (R 4.2.2) #> landscapemetrics * 2.0.0 2023-09-20 [1] Github (r-spatialecology/landscapemetrics@7d2aa82) #> P lifecycle 1.0.3 2022-10-07 [4] CRAN (R 4.2.2) #> P magrittr 2.0.3 2022-03-30 [4] CRAN (R 4.2.2) #> P mime 0.12 2021-09-28 [4] CRAN (R 4.2.0) #> pillar 1.9.0 2023-03-22 [1] CRAN (R 4.2.3) #> P pkgconfig 2.0.3 2019-09-22 [4] CRAN (R 4.2.2) #> P purrr 0.3.5 2022-10-06 [4] CRAN (R 4.2.2) #> R.cache 0.16.0 2022-07-21 [1] CRAN (R 4.2.3) #> R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.2.2) #> R.oo 1.25.0 2022-06-12 [1] CRAN (R 4.2.2) #> R.utils 2.12.2 2022-11-11 [1] CRAN (R 4.2.3) #> P R6 2.5.1 2021-08-19 [4] CRAN (R 4.2.2) #> P Rcpp 1.0.11 2023-07-06 [4] CRAN (R 4.2.3) #> reprex 2.0.2 2022-08-17 [1] CRAN (R 4.2.3) #> rlang 1.1.1 2023-04-28 [1] CRAN (R 4.2.3) #> P rmarkdown 2.19 2022-12-15 [4] CRAN (R 4.2.2) #> P sessioninfo 1.2.2 2021-12-06 [4] CRAN (R 4.2.2) #> P stringi 1.7.8 2022-07-11 [4] CRAN (R 4.2.1) #> P stringr 1.5.0 2022-12-02 [4] CRAN (R 4.2.2) #> styler 1.10.2 2023-08-29 [1] CRAN (R 4.2.3) #> terra 1.7-46 2023-09-06 [1] CRAN (R 4.2.3) #> tibble 3.2.1 2023-03-20 [1] CRAN (R 4.2.3) #> utf8 1.2.3 2023-01-31 [1] CRAN (R 4.2.3) #> vctrs 0.6.3 2023-06-14 [1] CRAN (R 4.2.3) #> P withr 2.5.0 2022-03-03 [4] CRAN (R 4.2.2) #> xfun 0.40 2023-08-09 [1] CRAN (R 4.2.3) #> P xml2 1.3.3 2021-11-30 [4] CRAN (R 4.2.2) #> P yaml 2.3.6 2022-10-18 [4] CRAN (R 4.2.2) #> #> [1] C:/Users/ferrariof/Documents/R_tests/TEST-landscapemetrics-2.0.0/packrat/lib/x86_64-w64-mingw32/4.2.2 #> [2] C:/Users/ferrariof/Documents/R_tests/TEST-landscapemetrics-2.0.0/packrat/lib-ext/x86_64-w64-mingw32/4.2.2 #> [3] C:/Users/ferrariof/Documents/R_tests/TEST-landscapemetrics-2.0.0/packrat/lib-R/x86_64-w64-mingw32/4.2.2 #> [4] C:/Users/ferrariof/AppData/Local/Programs/R/R-4.2.2/library #> #> P ── Loaded and on-disk path mismatch. #> #> ────────────────────────────────────────────────────────────────────────────── ```
mhesselbarth commented 1 year ago

Hey, here are a few comments that might already help

  1. We just updated the calculation of the shape index, following the definition on the official FRAGSTATS homepage

    SHAPE equals patch perimeter (m) divided by the square root of patch area (m2), adjusted by a constant to adjust for a square standard.

https://fragstats.org/index.php/fragstats-metrics/shape-metrics/p2-shape-index

  1. The documentation of the mean shape index was wrong and should be corrected now. The mean value is shape_mn = 1 if all patches are squares (d2499c7e).

  2. If none of the above points helped, we are currently also re-writting the window function, so maybe give us a few days until thats updated and move on from there?

Nowosad commented 1 year ago

@mhesselbarth could you add an info about the lsm_p_shape changes to the NEWS file?

Nowosad commented 1 year ago

@filippo-ferrario the circular window now works as you expected. However, please recheck the square one:

# remotes::install_github("r-spatialecology/landscapemetrics@commoncalcs")
devtools::load_all()
#> ℹ Loading landscapemetrics
landscape <- terra::rast(landscapemetrics::landscape)
terra::plot(landscape)

# polygonized for later plotting
pol_r1<-terra::as.polygons(landscape) 
terra::plot(subset(pol_r1, pol_r1$clumps==1), add=T, border='red', lwd=2)
terra::plot(subset(pol_r1, pol_r1$clumps==2), add=T, border='purple', lwd=2)


class1<-landscape==1
class1[class1==0]<-NA
terra::plot(class1, col='black')


mw<-terra::focalMat(landscape, d=1, type=c('rectangle'), fillNA=T)
mw<-ifelse(mw>0, 1,NA)
mw
#>      [,1] [,2] [,3]
#> [1,]    1    1    1
#> [2,]    1    1    1
#> [3,]    1    1    1

mwc<-terra::focalMat(landscape, d=1, type=c('circle'), fillNA=T)
mwc<-ifelse(mwc>0, 1,NA)
mwc
#>      [,1] [,2] [,3]
#> [1,]   NA    1   NA
#> [2,]    1    1    1
#> [3,]   NA    1   NA

sh<-window_lsm(class1,what='lsm_l_shape_mn', window=mw)
sh$layer_1$lsm_l_shape_mn
#> class       : SpatRaster 
#> dimensions  : 30, 30, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 30, 0, 30  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source(s)   : memory
#> name        : clumps 
#> min value   :    1.0 
#> max value   :    1.5

# plot the patch and values
terra::plot(class1)
terra::plot(subset(pol_r1, pol_r1$clumps==1), add=T, border='red', lwd=2)
terra::text(sh$layer_1$lsm_l_shape_mn, digits=3, cex=0.5)

# highlight region of interest
v1<-terra::vect(rbind(c(18,8),c(21,8),c(21,11),c(18,11)),type='polygons')
v2<-terra::vect(rbind(c(16,19),c(19,19),c(19,22),c(16,22)),type='polygons')
v3<-terra::vect(rbind(c(8,9),c(11,9),c(11,12),c(8,12)),type='polygons')

terra::plot(v1, add=TRUE,border='black', lwd=3)
terra::plot(v2, add=TRUE,border='blue', lwd=3)
terra::plot(v3, add=TRUE, lwd=3, border='green')


csh<-window_lsm(class1,what='lsm_l_shape_mn', window=mwc)
csh$layer_1$lsm_l_shape_mn
#> class       : SpatRaster 
#> dimensions  : 30, 30, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 30, 0, 30  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source(s)   : memory
#> name        :   clumps 
#> min value   : 1.000000 
#> max value   : 1.732051

terra::plot(class1)
terra::plot(subset(pol_r1, pol_r1$clumps==1), add=T, border='red', lwd=2)
terra::text(csh$layer_1$lsm_l_shape_mn, digits=3, cex=0.5)

terra::plot(v1, add=TRUE,border='black', lwd=3)
terra::plot(v2, add=TRUE,border='blue', lwd=3)
terra::plot(v3, add=TRUE, lwd=3, border='green')

Created on 2023-10-08 with reprex v2.0.2

mhesselbarth commented 1 year ago

@mhesselbarth could you add an info about the lsm_p_shape changes to the NEWS file?

81501e5d

filippo-ferrario commented 1 year ago

Hi Jakub, Thanks for your work. Can you clarify what I should check to help?

Thanks again Filippo

Il giorno dom 8 ott 2023 alle ore 09:32 Jakub Nowosad < @.***> ha scritto:

@filippo-ferrario https://github.com/filippo-ferrario the circular window now works as you expected. However, please recheck the square one:

@.***")devtools::load_all()#> ℹ Loading landscapemetricslandscape <- terra::rast(landscapemetrics::landscape)terra::plot(landscape)

polygonized for later plottingpol_r1<-terra::as.polygons(landscape) terra::plot(subset(pol_r1, pol_r1$clumps==1), add=T, border='red', lwd=2)terra::plot(subset(pol_r1, pol_r1$clumps==2), add=T, border='purple', lwd=2)

https://camo.githubusercontent.com/dc0859e899d0fd168488e1741aa4d9c213379656f02bbab79e44cb47f3731378/68747470733a2f2f692e696d6775722e636f6d2f4252345457416c2e706e67

class1<-landscape==1class1[class1==0]<-NAterra::plot(class1, col='black')

https://camo.githubusercontent.com/3bc9f2fe7a44910fb0a07cabb01f4bdbce0948850aeec13d1c668d6522423a3a/68747470733a2f2f692e696d6775722e636f6d2f7845536f5a4d352e706e67

mw<-terra::focalMat(landscape, d=1, type=c('rectangle'), fillNA=T)mw<-ifelse(mw>0, 1,NA)mw#> [,1] [,2] [,3]#> [1,] 1 1 1#> [2,] 1 1 1#> [3,] 1 1 1 mwc<-terra::focalMat(landscape, d=1, type=c('circle'), fillNA=T)mwc<-ifelse(mwc>0, 1,NA)mwc#> [,1] [,2] [,3]#> [1,] NA 1 NA#> [2,] 1 1 1#> [3,] NA 1 NA sh<-window_lsm(class1,what='lsm_l_shape_mn', window=mw)sh$layer_1$lsm_l_shape_mn#> class : SpatRaster #> dimensions : 30, 30, 1 (nrow, ncol, nlyr)#> resolution : 1, 1 (x, y)#> extent : 0, 30, 0, 30 (xmin, xmax, ymin, ymax)#> coord. ref. : #> source(s) : memory#> name : clumps #> min value : 1.0 #> max value : 1.5

plot the patch and valuesterra::plot(class1)terra::plot(subset(pol_r1, pol_r1$clumps==1), add=T, border='red', lwd=2)terra::text(sh$layer_1$lsm_l_shape_mn, digits=3, cex=0.5)

highlight region of interestv1<-terra::vect(rbind(c(18,8),c(21,8),c(21,11),c(18,11)),type='polygons')v2<-terra::vect(rbind(c(16,19),c(19,19),c(19,22),c(16,22)),type='polygons')v3<-terra::vect(rbind(c(8,9),c(11,9),c(11,12),c(8,12)),type='polygons')

terra::plot(v1, add=TRUE,border='black', lwd=3)terra::plot(v2, add=TRUE,border='blue', lwd=3)terra::plot(v3, add=TRUE, lwd=3, border='green')

https://camo.githubusercontent.com/bbdb96de64f03f67a9dfe95ddd41453c111b3b4bac1b15a0beffba65c6efbf1b/68747470733a2f2f692e696d6775722e636f6d2f6c316e345539642e706e67

csh<-window_lsm(class1,what='lsm_l_shape_mn', window=mwc)csh$layer_1$lsm_l_shape_mn#> class : SpatRaster #> dimensions : 30, 30, 1 (nrow, ncol, nlyr)#> resolution : 1, 1 (x, y)#> extent : 0, 30, 0, 30 (xmin, xmax, ymin, ymax)#> coord. ref. : #> source(s) : memory#> name : clumps #> min value : 1.000000 #> max value : 1.732051 terra::plot(class1)terra::plot(subset(pol_r1, pol_r1$clumps==1), add=T, border='red', lwd=2)terra::text(csh$layer_1$lsm_l_shape_mn, digits=3, cex=0.5) terra::plot(v1, add=TRUE,border='black', lwd=3)terra::plot(v2, add=TRUE,border='blue', lwd=3)terra::plot(v3, add=TRUE, lwd=3, border='green')

https://camo.githubusercontent.com/1377e600ddc74df16306569b5f1f00faf1c73f1521e72cb782dce0c834339e76/68747470733a2f2f692e696d6775722e636f6d2f7763706b5a39712e706e67

Created on 2023-10-08 with reprex v2.0.2 https://reprex.tidyverse.org

— Reply to this email directly, view it on GitHub https://github.com/r-spatialecology/landscapemetrics/issues/300#issuecomment-1752029825, or unsubscribe https://github.com/notifications/unsubscribe-auth/AJH74QAA7XGDZXTPMTDTK63X6KTQFAVCNFSM6AAAAAA5B53KASVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTONJSGAZDSOBSGU . You are receiving this because you were mentioned.Message ID: @.***>

Nowosad commented 1 year ago

@filippo-ferrario please go through your example squares one by one, try to recalculate the value by hand, and let us know if there are any differences.

filippo-ferrario commented 1 year ago

@filippo-ferrario please go through your example squares one by one, try to recalculate the value by hand, and let us know if there are any differences.

Hi Jakub, I did manually calculate the values for the square moving window in the black quadrat and I got 1.154701, which I think it is the value you also obtained.

Also, I tried another couple of cases (marked as yellow and purple squares in the attached image)

square_mw_jacub

The focal cell of the yellow square (i.e. the one circled in yellow) I manually calculate a value of 1.25. The focal cell of the purple square (i.e. the one circled in purple) I manually calculate a value of 1.07735.

Since I did not downloaded any revised function yet, I am basing my comparison on the numbers I can read from the plots in your reply (i.e., 300#issuecomment-1752029825). And since the figures in the plot overlap I can't guarantee that my results are exactly what it is printed in the plot: maybe double check.

I concurrently did some check on the values you obtained with the circular moving window and it seems also ok.

Please let me know when and how I can download a version of the package with these issue resolved (as I can see that my issue was included in a broader revision work).

Thanks again for your work.

mhesselbarth commented 1 year ago

Hey,

just run remotes::install_github("r-spatialecology/landscapemetrics") (after installing the remotes package if you don't already have it).

filippo-ferrario commented 1 year ago

I downloaded the updated package from github as indecated by @mhesselbarth (i.e., remotes::install_github(“r-spatialecology/landscapemetrics”)).
I also updated the packages suggested during the installation.

However, testing the function on my pc, I do not obtain the same results as you do when using the circular moving window.
Look at the blue and black squares in particular, but numbers are different also in other cells.

Follows the code

library(landscapemetrics)
# Setup data

Import raster data from Terra package.

landscape <- terra::rast(landscapemetrics::landscape)
terra::plot(landscape)

# polygonized for later plotting
pol_r1<-terra::as.polygons(landscape) 
terra::plot(subset(pol_r1, pol_r1$clumps==1), add=T, border='red', lwd=2)
terra::plot(subset(pol_r1, pol_r1$clumps==2), add=T, border='purple', lwd=2)

# terra::plot(subset(pol_r1, pol_r1$clumps==3), add=T, border='black', lwd=2)

select just one class and replace the False with NA (this mimiks the behaviour of raster::rasterize used in which have a default background to NA, as it has been done in preprocessing-hyperframe_for_modeling.R)

class1<-landscape==1
class1[class1==0]<-NA
terra::plot(class1, col='black')


# Create Rectangular moving window of side length 3
mw<-terra::focalMat(landscape, d=1, type=c('rectangle'), fillNA=T)
mw<-ifelse(mw>0, 1,NA)
mw
#>      [,1] [,2] [,3]
#> [1,]    1    1    1
#> [2,]    1    1    1
#> [3,]    1    1    1

# Create Circular moving window of radius 1 (it looks like a cross but it is easier to picture it on the raster)
mwc<-terra::focalMat(landscape, d=1, type=c('circle'), fillNA=T)
mwc<-ifelse(mwc>0, 1,NA)
mwc
#>      [,1] [,2] [,3]
#> [1,]   NA    1   NA
#> [2,]    1    1    1
#> [3,]   NA    1   NA

Rectangular case for Shape index

sh<-window_lsm(class1,what='lsm_l_shape_mn', window=mw)
sh$layer_1$lsm_l_shape_mn
#> class       : SpatRaster 
#> dimensions  : 30, 30, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 30, 0, 30  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source(s)   : memory
#> name        : clumps 
#> min value   :    1.0 
#> max value   :    1.5
# plot the patch and values
terra::plot(class1)
terra::plot(subset(pol_r1, pol_r1$clumps==1), add=T, border='red', lwd=2)
terra::text(sh$layer_1$lsm_l_shape_mn, digits=3, cex=0.5)

# highlight region of interest
v1<-terra::vect(rbind(c(18,8),c(21,8),c(21,11),c(18,11)),type='polygons')
v2<-terra::vect(rbind(c(16,19),c(19,19),c(19,22),c(16,22)),type='polygons')
v3<-terra::vect(rbind(c(8,9),c(11,9),c(11,12),c(8,12)),type='polygons')

terra::plot(v1, add=TRUE,border='black', lwd=3)
terra::plot(v2, add=TRUE,border='blue', lwd=3)
terra::plot(v3, add=TRUE, lwd=3, border='green')

Circular case for Shape index

csh<-window_lsm(class1,what='lsm_l_shape_mn', window=mwc)
#> Warning: number of items to replace is not a multiple of replacement length
csh$layer_1$lsm_l_shape_mn
#> class       : SpatRaster 
#> dimensions  : 30, 30, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 30, 0, 30  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source(s)   : memory
#> name        : clumps 
#> min value   :      1 
#> max value   :      2
# plot the patch and values
terra::plot(class1)
terra::plot(subset(pol_r1, pol_r1$clumps==1), add=T, border='red', lwd=2)
terra::text(csh$layer_1$lsm_l_shape_mn, digits=3, cex=0.5)

terra::plot(v1, add=TRUE,border='black', lwd=3)
terra::plot(v2, add=TRUE,border='blue', lwd=3)
terra::plot(v3, add=TRUE, lwd=3, border='green')

Created on 2023-10-13 with reprex v2.0.2

Session info ``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.2.2 (2022-10-31 ucrt) #> os Windows 10 x64 (build 19044) #> system x86_64, mingw32 #> ui RTerm #> language ENG #> collate English_United States.utf8 #> ctype English_United States.utf8 #> tz America/New_York #> date 2023-10-13 #> pandoc 2.18 @ C:/Users/FERRAR~1/AppData/Local/Pandoc/ (via rmarkdown) #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> ! package * version date (UTC) lib source #> cli 3.6.1 2023-03-23 [1] CRAN (R 4.2.3) #> P codetools 0.2-18 2020-11-04 [4] CRAN (R 4.2.2) #> P curl 4.3.3 2022-10-06 [4] CRAN (R 4.2.2) #> P digest 0.6.31 2022-12-11 [4] CRAN (R 4.2.2) #> P evaluate 0.19 2022-12-13 [4] CRAN (R 4.2.2) #> fansi 1.0.5 2023-10-08 [1] CRAN (R 4.2.3) #> P fastmap 1.1.0 2021-01-25 [4] CRAN (R 4.2.2) #> P fs 1.5.2 2021-12-08 [4] CRAN (R 4.2.2) #> P glue 1.6.2 2022-02-24 [4] CRAN (R 4.2.2) #> P highr 0.9 2021-04-16 [4] CRAN (R 4.2.2) #> P htmltools 0.5.4 2022-12-07 [4] CRAN (R 4.2.2) #> P httr 1.4.4 2022-08-17 [4] CRAN (R 4.2.2) #> P knitr 1.41 2022-11-18 [4] CRAN (R 4.2.2) #> landscapemetrics * 2.0.0 2023-10-13 [1] Github (r-spatialecology/landscapemetrics@81501e5) #> P lifecycle 1.0.3 2022-10-07 [4] CRAN (R 4.2.2) #> P magrittr 2.0.3 2022-03-30 [4] CRAN (R 4.2.2) #> P mime 0.12 2021-09-28 [4] CRAN (R 4.2.0) #> pillar 1.9.0 2023-03-22 [1] CRAN (R 4.2.3) #> P pkgconfig 2.0.3 2019-09-22 [4] CRAN (R 4.2.2) #> P purrr 0.3.5 2022-10-06 [4] CRAN (R 4.2.2) #> R.cache 0.16.0 2022-07-21 [1] CRAN (R 4.2.3) #> R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.2.2) #> R.oo 1.25.0 2022-06-12 [1] CRAN (R 4.2.2) #> R.utils 2.12.2 2022-11-11 [1] CRAN (R 4.2.3) #> P R6 2.5.1 2021-08-19 [4] CRAN (R 4.2.2) #> P Rcpp 1.0.11 2023-07-06 [4] CRAN (R 4.2.3) #> reprex 2.0.2 2022-08-17 [1] CRAN (R 4.2.3) #> rlang 1.1.1 2023-04-28 [1] CRAN (R 4.2.3) #> P rmarkdown 2.19 2022-12-15 [4] CRAN (R 4.2.2) #> P sessioninfo 1.2.2 2021-12-06 [4] CRAN (R 4.2.2) #> P stringi 1.7.8 2022-07-11 [4] CRAN (R 4.2.1) #> P stringr 1.5.0 2022-12-02 [4] CRAN (R 4.2.2) #> styler 1.10.2 2023-08-29 [1] CRAN (R 4.2.3) #> terra 1.7-46 2023-09-06 [1] CRAN (R 4.2.3) #> tibble 3.2.1 2023-03-20 [1] CRAN (R 4.2.3) #> utf8 1.2.3 2023-01-31 [1] CRAN (R 4.2.3) #> vctrs 0.6.4 2023-10-12 [1] CRAN (R 4.2.3) #> withr 2.5.1 2023-09-26 [1] CRAN (R 4.2.3) #> xfun 0.40 2023-08-09 [1] CRAN (R 4.2.3) #> P xml2 1.3.3 2021-11-30 [4] CRAN (R 4.2.2) #> P yaml 2.3.6 2022-10-18 [4] CRAN (R 4.2.2) #> #> [1] C:/Users/ferrariof/Documents/R_tests/TEST-landscapemetrics-2.0.0/packrat/lib/x86_64-w64-mingw32/4.2.2 #> [2] C:/Users/ferrariof/Documents/R_tests/TEST-landscapemetrics-2.0.0/packrat/lib-ext/x86_64-w64-mingw32/4.2.2 #> [3] C:/Users/ferrariof/Documents/R_tests/TEST-landscapemetrics-2.0.0/packrat/lib-R/x86_64-w64-mingw32/4.2.2 #> [4] C:/Users/ferrariof/AppData/Local/Programs/R/R-4.2.2/library #> #> P ── Loaded and on-disk path mismatch. #> #> ────────────────────────────────────────────────────────────────────────────── ```
Nowosad commented 1 year ago

Hi @filippo-ferrario -- sorry for the confusion -- the proper installation code is remotes::install_github("r-spatialecology/landscapemetrics@commoncalcs"). These changes are still in their own branch -- I plan to merge it when you confirm that the circular window works correctly. Please recheck it and let me know.

mhesselbarth commented 1 year ago

Oh that's my bad. I thought the changes on main should already address the issue. Didn't realize it's still on another branch,

filippo-ferrario commented 1 year ago

I would say now it works. But see my example for the case when NA is produced.

TEST shape index updates

library(landscapemetrics)
library(terra)
#> Warning: package 'terra' was built under R version 4.2.3
#> terra 1.7.55

landscape <- terra::rast(landscapemetrics::landscape)

terra::plot(landscape)


res_rect<-data.frame(win=c('black','blue', 'green', 'purple','gray','orange'), cell=c(790,288,580,620,272,436), calc=NA, man=NA, check=NA)
res_rect
#>      win cell calc man check
#> 1  black  790   NA  NA    NA
#> 2   blue  288   NA  NA    NA
#> 3  green  580   NA  NA    NA
#> 4 purple  620   NA  NA    NA
#> 5   gray  272   NA  NA    NA
#> 6 orange  436   NA  NA    NA

# Example Rectangular moving windows
mw_list<-list(
        terra::vect(rbind(c(8,2),c(11,2),c(11,5),c(8,5)),type='polygons'),
        terra::vect(rbind(c(16,19),c(19,19),c(19,22),c(16,22)),type='polygons'),
        terra::vect(rbind(c(8,9),c(11,9),c(11,12),c(8,12)),type='polygons'),
        terra::vect(rbind(c(18,8),c(21,8),c(21,11),c(18,11)),type='polygons'),
        terra::vect(rbind(c(0,19),c(3,19),c(3,22),c(0,22)),type='polygons'),
        terra::vect(rbind(c(14,14),c(17,14),c(17,17),c(14,17)),type='polygons')
        )

# Example Circular moving windows

res_circ<-data.frame(win=c('black','blue', 'green', 'purple','gray','orange'), cell=c(790,257,580,652,272,436), calc=NA, man=NA, check=NA)
res_circ
#>      win cell calc man check
#> 1  black  790   NA  NA    NA
#> 2   blue  257   NA  NA    NA
#> 3  green  580   NA  NA    NA
#> 4 purple  652   NA  NA    NA
#> 5   gray  272   NA  NA    NA
#> 6 orange  436   NA  NA    NA

mwc_list<-terra::xyFromCell(landscape,res_circ$cell) |>
apply(MARGIN=1,function(x){
    # browser()
    pl<-terra::vect(rbind( c(x['x']-0.5,x['y']-1.5), c(x['x']+0.5,x['y']-1.5), 
                       c(x['x']+0.5,x['y']-0.5), c(x['x']+1.5,x['y']-0.5),
                       c(x['x']+1.5,x['y']+0.5), c(x['x']+0.5,x['y']+0.5),
                       c(x['x']+0.5,x['y']+1.5), c(x['x']-0.5,x['y']+1.5),
                       c(x['x']-0.5,x['y']+0.5), c(x['x']-1.5,x['y']+0.5),
                       c(x['x']-1.5,x['y']-0.5), c(x['x']-0.5,x['y']-0.5)
                       ) ,type='polygons'  ) 
})

# Create rectangular moving window
mw<-terra::focalMat(landscape, d=1, type=c('rectangle'), fillNA=T)
mw<-ifelse(mw>0, 1,NA)
mw
#>      [,1] [,2] [,3]
#> [1,]    1    1    1
#> [2,]    1    1    1
#> [3,]    1    1    1

# Create circular moving window
mwcir<-terra::focalMat(landscape, d=1, type=c('circle'), fillNA=T)
mwcir<-ifelse(mwcir>0, 1,NA)
mwcir
#>      [,1] [,2] [,3]
#> [1,]   NA    1   NA
#> [2,]    1    1    1
#> [3,]   NA    1   NA

class1<-landscape==1
class1[class1==0]<-NA

class3<-landscape==3
class3[class3==0]<-NA

shape_sqr<-window_lsm(class1,what=c('lsm_l_shape_mn' ), window=mw)
shape_circ<-window_lsm(class3,what=c('lsm_l_shape_mn' ), window=mwcir)

Rectangular case Check

terra::plot(class1,col='yellow')
for (i in 1: length(mw_list) ) {terra::plot(mw_list[[i]], add=T, border=res_rect$win[i], lwd=3)}
terra::text(shape_sqr$layer_1$lsm_l_shape_mn, digits=3, cex=0.5)

Checks

res_out<- res_rect
res_out$man<-round(c('black'= 0,
               'blue' = mean(c(4/(sqrt(1)*4), 4/(sqrt(1)*4))), 
               'green'= 12/(sqrt(9)*4), 
               'purple'= 8/(sqrt(3)*4),
               'gray'=  10/(sqrt(4)*4), 
               'orange'= mean(c(4/(sqrt(1)*4) , 8/(sqrt(3)*4)))  ), 
               3) 
res_out$calc<-round(terra::extract(shape_sqr$layer_1$lsm_l_shape_mn, res_rect$cell), 3)
res_out$check<- with(res_out, man==calc)
res_out
#>      win cell clumps   man clumps
#> 1  black  790     NA 0.000     NA
#> 2   blue  288  1.000 1.000   TRUE
#> 3  green  580  1.000 1.000   TRUE
#> 4 purple  620  1.155 1.155   TRUE
#> 5   gray  272  1.250 1.250   TRUE
#> 6 orange  436  1.077 1.077   TRUE

all ok? NA
All checks (last columns dataframe) ok except the NA

WARNING: when the moving window does not contain any patch the result is NA instead of 0. This is to be corrected before analysing the data.
FYI: the same behaviour noted for lsm_l_np, lsm_l_area_mn, lsm_l_ta, lsm_l_para_mn when calculated on masked classes as in this example (these are the one I know…potentially others).

Circular case Check

terra::plot(class3,col='yellow')
for (i in 1: length(mwc_list) ) {terra::plot(mwc_list[[i]], add=T, border=res_circ$win[i], lwd=3)}
terra::text(shape_circ$layer_1$lsm_l_shape_mn, digits=3, cex=0.5)

Checks

res_out<- res_circ 
res_out$man<-round(c('black'= mean(c( 4/(sqrt(1)*4 ) ) ),
               'blue' = mean(c( 8/(sqrt(2)*4 ) ) ), 
               'green'= 0, 
               'purple'= mean(c( 4/(sqrt(1)*4 ), 4/(sqrt(1)*4 ) ) ),
               'gray'=  4/(sqrt(1)*4), 
               'orange'= mean(c( 8/(sqrt(3)*4 ) ) )  ), 
               3) 
res_out$calc<-round(terra::extract(shape_circ$layer_1$lsm_l_shape_mn, res_circ$cell), 3)
res_out$check<- with(res_out, man==calc)
res_out
#>      win cell clumps   man clumps
#> 1  black  790  1.000 1.000   TRUE
#> 2   blue  257  1.414 1.414   TRUE
#> 3  green  580     NA 0.000     NA
#> 4 purple  652  1.000 1.000   TRUE
#> 5   gray  272  1.000 1.000   TRUE
#> 6 orange  436  1.155 1.155   TRUE

all ok? NA All checks (last columns dataframe) ok except the NA

WARNING: when the moving window does not contain any patch the result is NA instead of 0. This is to be corrected before analysing the data.
FYI: the same behaviour noted for lsm_l_np, lsm_l_area_mn, lsm_l_ta, lsm_l_para_mn when calculated on masked classes as in this example (these are the one I know…potentially others).

Created on 2023-10-26 with reprex v2.0.2

Session info ``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.2.2 (2022-10-31 ucrt) #> os Windows 10 x64 (build 19044) #> system x86_64, mingw32 #> ui RTerm #> language ENG #> collate English_United States.utf8 #> ctype English_United States.utf8 #> tz America/New_York #> date 2023-10-26 #> pandoc 2.18 @ C:/Users/FERRAR~1/AppData/Local/Pandoc/ (via rmarkdown) #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> ! package * version date (UTC) lib source #> cli 3.6.1 2023-03-23 [1] CRAN (R 4.2.3) #> P codetools 0.2-18 2020-11-04 [4] CRAN (R 4.2.2) #> P curl 4.3.3 2022-10-06 [4] CRAN (R 4.2.2) #> P digest 0.6.31 2022-12-11 [4] CRAN (R 4.2.2) #> P evaluate 0.19 2022-12-13 [4] CRAN (R 4.2.2) #> fansi 1.0.5 2023-10-08 [1] CRAN (R 4.2.3) #> P fastmap 1.1.0 2021-01-25 [4] CRAN (R 4.2.2) #> P fs 1.5.2 2021-12-08 [4] CRAN (R 4.2.2) #> P glue 1.6.2 2022-02-24 [4] CRAN (R 4.2.2) #> P highr 0.9 2021-04-16 [4] CRAN (R 4.2.2) #> P htmltools 0.5.4 2022-12-07 [4] CRAN (R 4.2.2) #> P httr 1.4.4 2022-08-17 [4] CRAN (R 4.2.2) #> P knitr 1.41 2022-11-18 [4] CRAN (R 4.2.2) #> landscapemetrics * 2.1.0 2023-10-26 [1] Github (r-spatialecology/landscapemetrics@1e7b0d4) #> P lifecycle 1.0.3 2022-10-07 [4] CRAN (R 4.2.2) #> P magrittr 2.0.3 2022-03-30 [4] CRAN (R 4.2.2) #> P mime 0.12 2021-09-28 [4] CRAN (R 4.2.0) #> pillar 1.9.0 2023-03-22 [1] CRAN (R 4.2.3) #> P pkgconfig 2.0.3 2019-09-22 [4] CRAN (R 4.2.2) #> P purrr 0.3.5 2022-10-06 [4] CRAN (R 4.2.2) #> R.cache 0.16.0 2022-07-21 [1] CRAN (R 4.2.3) #> R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.2.2) #> R.oo 1.25.0 2022-06-12 [1] CRAN (R 4.2.2) #> R.utils 2.12.2 2022-11-11 [1] CRAN (R 4.2.3) #> P R6 2.5.1 2021-08-19 [4] CRAN (R 4.2.2) #> P Rcpp 1.0.11 2023-07-06 [4] CRAN (R 4.2.3) #> reprex 2.0.2 2022-08-17 [1] CRAN (R 4.2.3) #> rlang 1.1.1 2023-04-28 [1] CRAN (R 4.2.3) #> P rmarkdown 2.19 2022-12-15 [4] CRAN (R 4.2.2) #> P sessioninfo 1.2.2 2021-12-06 [4] CRAN (R 4.2.2) #> P stringi 1.7.8 2022-07-11 [4] CRAN (R 4.2.1) #> P stringr 1.5.0 2022-12-02 [4] CRAN (R 4.2.2) #> styler 1.10.2 2023-08-29 [1] CRAN (R 4.2.3) #> terra * 1.7-55 2023-10-13 [1] CRAN (R 4.2.3) #> tibble 3.2.1 2023-03-20 [1] CRAN (R 4.2.3) #> utf8 1.2.4 2023-10-22 [1] CRAN (R 4.2.3) #> vctrs 0.6.4 2023-10-12 [1] CRAN (R 4.2.3) #> withr 2.5.1 2023-09-26 [1] CRAN (R 4.2.3) #> xfun 0.40 2023-08-09 [1] CRAN (R 4.2.3) #> P xml2 1.3.3 2021-11-30 [4] CRAN (R 4.2.2) #> P yaml 2.3.6 2022-10-18 [4] CRAN (R 4.2.2) #> #> [1] C:/Users/ferrariof/Documents/R_tests/TEST-landscapemetrics-2.0.0/packrat/lib/x86_64-w64-mingw32/4.2.2 #> [2] C:/Users/ferrariof/Documents/R_tests/TEST-landscapemetrics-2.0.0/packrat/lib-ext/x86_64-w64-mingw32/4.2.2 #> [3] C:/Users/ferrariof/Documents/R_tests/TEST-landscapemetrics-2.0.0/packrat/lib-R/x86_64-w64-mingw32/4.2.2 #> [4] C:/Users/ferrariof/AppData/Local/Programs/R/R-4.2.2/library #> #> P ── Loaded and on-disk path mismatch. #> #> ────────────────────────────────────────────────────────────────────────────── ```
Nowosad commented 1 year ago

Thanks @filippo-ferrario! @mhesselbarth please take a look at the above examples. Should we postprocess some results or just add an explanation to the help file?

mhesselbarth commented 1 year ago

I would say put an explanation to the help file? I feel like this is such a specific case

Nowosad commented 1 year ago

Thank you @filippo-ferrario for letting us know about the initial issues and providing all of the examples!

@mhesselbarth thanks -- I have done that. This is one of those NA or background situations that are probably better to be decided by the users.

filippo-ferrario commented 1 year ago

I'm glad it helped. I guess in this moment I should still use the remotes::install_github("r-spatialecology/landscapemetrics@commoncalcs") branch, isn't it?

Thank you for your work. Filippo

Nowosad commented 1 year ago

Yes -- you should still use the branch. I think it will be merge in a foreseeable future to the main branch as it brings a lot of performance improvements and other fixes (see https://github.com/r-spatialecology/landscapemetrics/blob/commoncalcs/NEWS.md#landscapemetrics-210).