circular window in window_lsm
filippo-ferrario opened this issue · 9 comments
I am new to landscapemetrics and I am experimenting with it.
I want to use window_lsm to calculate metrics within a ‘buffer’ around each pixel, and the first metric I was testing was the mean area of patches of a given class within the buffer.
I would like to use a circular moving window rather than a square one.
However the circle doesn’t seem to work correctly, as opposed to the square.
Here a synthetic reproducible example where I create a raster with only 1 patch of 9 m2, surrounded by NAs.
library(raster)
#> Loading required package: sp
#> Warning: package 'sp' was built under R version 4.1.3
library(landscapemetrics)
# ===================
# create raster
# ===================
r1<-raster(xmn=594124,xmx=594145,ymn=5461881,ymx=5461902, res=1, crs='+proj=utm +zone=19 +datum=WGS84 +units=m +no_defs +type=crs')
r1[10:12,10:12]<-1 # create the patch
r1
#> class : RasterLayer
#> dimensions : 21, 21, 441 (nrow, ncol, ncell)
#> resolution : 1, 1 (x, y)
#> extent : 594124, 594145, 5461881, 5461902 (xmin, xmax, ymin, ymax)
#> crs : +proj=utm +zone=19 +datum=WGS84 +units=m +no_defs
#> source : memory
#> names : layer
#> values : 1, 1 (min, max)
plot(r1)
pol_r1<-rasterToPolygons(r1) # polygonized for later plottingCase 1: using squared matrix (it works)
sq<-matrix(1, ncol=5,nrow=5)
sq
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1 1 1 1 1
#> [2,] 1 1 1 1 1
#> [3,] 1 1 1 1 1
#> [4,] 1 1 1 1 1
#> [5,] 1 1 1 1 1
sq_1class<-window_lsm(r1,what='lsm_l_area_mn', window=sq)
sq_1class$layer_1$lsm_l_area_mn<-sq_1class$layer_1$lsm_l_area_mn*10000 # convert values in meter squared rather than Hecatres.
raster::plot(sq_1class$layer_1$lsm_l_area_mn, colNA='red')
plot(pol_r1, add=T)
text(sq_1class$layer_1$lsm_l_area_mn, digits=3,cex=0.5)
this is looking good.
The patch is actually composed of 9 cells of 1 m2
Problematic cases with circular window.
Case2: Using circular window as created by raster with 0 padding
circ0<-raster::focalWeight(r1, d=2, type=c('circle'), fillNA=F)
circ0
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.00000000 0.00000000 0.07692308 0.00000000 0.00000000
#> [2,] 0.00000000 0.07692308 0.07692308 0.07692308 0.00000000
#> [3,] 0.07692308 0.07692308 0.07692308 0.07692308 0.07692308
#> [4,] 0.00000000 0.07692308 0.07692308 0.07692308 0.00000000
#> [5,] 0.00000000 0.00000000 0.07692308 0.00000000 0.00000000
circ0_1class<-window_lsm(r1,what='lsm_l_area_mn', window=circ0)
#> Warning: no non-missing arguments to max; returning -Inf
circ0_1class$layer_1$lsm_l_area_mn<-circ0_1class$layer_1$lsm_l_area_mn*10000 # convert values in meter squared rather than Hecatres.
raster::plot(circ0_1class$layer_1$lsm_l_area_mn, colNA='red')
plot(pol_r1, add=T)
text(circ0_1class$layer_1$lsm_l_area_mn, digits=3,cex=0.5)
- Why the NA in the central cell?
- The numbers do not make sense
Case 3: using circular window replacing values >0 with ‘1’
circ1<-ifelse(circ0>0, 1,0)
circ1_1class<-window_lsm(r1,what='lsm_l_area_mn', window=circ1)
circ1_1class$layer_1$lsm_l_area_mn<-circ1_1class$layer_1$lsm_l_area_mn*10000 # convert to meters
raster::plot(circ1_1class$layer_1$lsm_l_area_mn, colNA='red')
plot(pol_r1, add=T)
text(circ1_1class$layer_1$lsm_l_area_mn, digits=3, cex=0.5)
- The numbers do not make sense: for example the two cells in the central row reporting a ‘3’ should read ‘7’ in theory. Similarly the corners cells of the patch that read ‘4.5’ should read ‘6’ instead, isn’t in?
However it seems that having values that are only ‘1’ or ‘0’ helps. Do the values in the circular matrix always need to be replaced by “1”?
Case 4: using circular window replacing values >0 with ‘1’ and padding with NA.
circ1NA<-ifelse(circ0>0, 1,NA)
circ1NA_1class<-window_lsm(r1,what='lsm_l_area_mn', window=circ1NA)
#> Warning: data length [13] is not a sub-multiple or multiple of the number of
#> rows [5]
circ1NA_1class$layer_1$lsm_l_area_mn<-circ1NA_1class$layer_1$lsm_l_area_mn*10000 # convert to meters
raster::plot(circ1NA_1class$layer_1$lsm_l_area_mn, colNA='red')
plot(pol_r1, add=T)
text(circ1NA_1class$layer_1$lsm_l_area_mn, digits=3, cex=0.5)
- The numbers do not make sense again…
Created on 2022-09-28 by the reprex package (v2.0.1)
Session info
sessioninfo::session_info()
#> - Session info ---------------------------------------------------------------
#> setting value
#> version R version 4.1.2 (2021-11-01)
#> os Windows 10 x64 (build 19042)
#> system x86_64, mingw32
#> ui RTerm
#> language ENG
#> collate English_United States.1252
#> ctype English_United States.1252
#> tz America/New_York
#> date 2022-09-28
#> pandoc 2.18 @ C:/Users/FERRAR~1/AppData/Local/Pandoc/ (via rmarkdown)
#>
#> - Packages -------------------------------------------------------------------
#> package * version date (UTC) lib source
#> cli 3.3.0 2022-04-25 [1] CRAN (R 4.1.3)
#> codetools 0.2-18 2020-11-04 [2] CRAN (R 4.1.2)
#> curl 4.3.2 2021-06-23 [1] CRAN (R 4.1.2)
#> digest 0.6.29 2021-12-01 [1] CRAN (R 4.1.2)
#> ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.1.2)
#> evaluate 0.14 2019-05-28 [1] CRAN (R 4.1.2)
#> fansi 1.0.2 2022-01-14 [1] CRAN (R 4.1.2)
#> fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.1.2)
#> fs 1.5.2 2021-12-08 [1] CRAN (R 4.1.2)
#> glue 1.6.1 2022-01-22 [1] CRAN (R 4.1.2)
#> highr 0.9 2021-04-16 [1] CRAN (R 4.1.2)
#> htmltools 0.5.2 2021-08-25 [1] CRAN (R 4.1.2)
#> httr 1.4.2 2020-07-20 [1] CRAN (R 4.1.2)
#> knitr 1.37 2021-12-16 [1] CRAN (R 4.1.2)
#> landscapemetrics * 1.5.5 2022-07-26 [1] Github (r-spatialecology/landscapemetrics@acd3fd6)
#> lattice 0.20-45 2021-09-22 [2] CRAN (R 4.1.2)
#> lifecycle 1.0.1 2021-09-24 [1] CRAN (R 4.1.2)
#> magrittr 2.0.2 2022-01-26 [1] CRAN (R 4.1.2)
#> mime 0.12 2021-09-28 [1] CRAN (R 4.1.1)
#> pillar 1.8.0 2022-07-18 [1] CRAN (R 4.1.3)
#> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.1.2)
#> purrr 0.3.4 2020-04-17 [1] CRAN (R 4.1.2)
#> R.cache 0.15.0 2021-04-30 [1] CRAN (R 4.1.3)
#> R.methodsS3 1.8.1 2020-08-26 [1] CRAN (R 4.1.1)
#> R.oo 1.24.0 2020-08-26 [1] CRAN (R 4.1.1)
#> R.utils 2.11.0 2021-09-26 [1] CRAN (R 4.1.3)
#> R6 2.5.1 2021-08-19 [1] CRAN (R 4.1.2)
#> raster * 3.5-15 2022-01-22 [1] CRAN (R 4.1.2)
#> Rcpp 1.0.8.3 2022-03-17 [1] CRAN (R 4.1.3)
#> reprex 2.0.1 2021-08-05 [1] CRAN (R 4.1.2)
#> rgdal 1.5-32 2022-05-09 [1] CRAN (R 4.1.3)
#> rlang 1.0.4 2022-07-12 [1] CRAN (R 4.1.3)
#> rmarkdown 2.11 2021-09-14 [1] CRAN (R 4.1.2)
#> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.1.3)
#> sp * 1.5-0 2022-06-05 [1] CRAN (R 4.1.3)
#> stringi 1.7.6 2021-11-29 [1] CRAN (R 4.1.2)
#> stringr 1.4.0 2019-02-10 [1] CRAN (R 4.1.2)
#> styler 1.7.0 2022-03-13 [1] CRAN (R 4.1.3)
#> terra 1.5-34 2022-06-09 [1] CRAN (R 4.1.3)
#> tibble 3.1.6 2021-11-07 [1] CRAN (R 4.1.2)
#> utf8 1.2.2 2021-07-24 [1] CRAN (R 4.1.2)
#> vctrs 0.3.8 2021-04-29 [1] CRAN (R 4.1.2)
#> withr 2.5.0 2022-03-03 [1] CRAN (R 4.1.3)
#> xfun 0.29 2021-12-14 [1] CRAN (R 4.1.2)
#> xml2 1.3.3 2021-11-30 [1] CRAN (R 4.1.2)
#> yaml 2.2.2 2022-01-25 [1] CRAN (R 4.1.2)
#>
#> [1] C:/Users/ferrariof/Documents/R/win-library/4.1
#> [2] C:/Program Files/R/R-4.1.2/library
#>
#> ------------------------------------------------------------------------------Hi @filippo-ferrario -- I tried to look at your questions, but I am a bit overwhelmed by your code and comments. Would you be able, for example, to prepare a separate (and small) reproducible example for each of the questions (you can post them as separate messages here)? Also -- could you use the {reprex} package (https://reprex.tidyverse.org/) to keep the code and the outputs together (try the wd = . argument)?
@Nowosad I have updated the issue and - hopefully - improved the reproducible example. I could not make separate examples for each question because it is basically only 1 question (i.e. why the circular window does not behave as one would expect?). To make my case I use a square window as a reference case that works, and then present the unexpected result I obtain when using the circular window.
I hope this is now better presented and that you could shed some light into what is happening.
Thank you
Thanks -- it is much cleaner now. I will try to find some time in the next few days to look into it.
@mhesselbarth @bitbacchus -- maybe also you can try to figure it out now?
I will try to have a look at it, however, my next few weeks are pretty busy unfortunately 😭
Hi @filippo-ferrario - I am just writing to (a) confirm the issue and (b) let you know that I found some solution. I will post my code tomorrow.
@filippo-ferrario -- the main issue was that the window_lsm() was implemented for regular (square) windows only. However, with a small change (after some debugging) I was able to add support for circular windows. The only important thing to remember -- the window needs to be set with 1s and NAs only.
@mhesselbarth -- see my related PR at #265
Code example:
library(raster)
library(landscapemetrics)
r1 <- raster(xmn = 594124, xmx = 594145, ymn = 5461881, ymx = 5461902,
res = 1, crs = "EPSG:32619")
r1[10:12, 10:12] <- 1 # create the patch
plot(r1)
pol_r1 <- rasterToPolygons(r1) # polygonized for later plotting
sq <- matrix(1, ncol = 5, nrow = 5)
sq
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1 1 1 1 1
#> [2,] 1 1 1 1 1
#> [3,] 1 1 1 1 1
#> [4,] 1 1 1 1 1
#> [5,] 1 1 1 1 1
sq_1class <- window_lsm(r1, what = "lsm_l_area_mn", window = sq)
sq_1class$layer_1$lsm_l_area_mn <- sq_1class$layer_1$lsm_l_area_mn * 10000 # convert values in meter squared rather than Hecatres.
raster::plot(sq_1class$layer_1$lsm_l_area_mn, colNA = "red")
plot(pol_r1, add = TRUE)
text(sq_1class$layer_1$lsm_l_area_mn, digits = 3, cex = 0.5)
circ0 <- raster::focalWeight(r1, d = 2, type = "circle", fillNA = TRUE)
circ0[!is.na(circ0)] = 1
circ0
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] NA NA 1 NA NA
#> [2,] NA 1 1 1 NA
#> [3,] 1 1 1 1 1
#> [4,] NA 1 1 1 NA
#> [5,] NA NA 1 NA NA
circ0_1class <- window_lsm(r1, what = "lsm_l_area_mn", window = circ0)
circ0_1class$layer_1$lsm_l_area_mn<-circ0_1class$layer_1$lsm_l_area_mn * 10000 # convert values in meter squared rather than Hecatres.
raster::plot(circ0_1class$layer_1$lsm_l_area_mn, colNA = "red")
plot(pol_r1, add = TRUE)
text(circ0_1class$layer_1$lsm_l_area_mn, digits = 3, cex = 0.5)
Created on 2022-10-07 with reprex v2.0.2
@Nowosad I will try to have a look next week since I have a few deadlines until Friday 👍
@Nowosad Thanks for looking into this issue and for finding a solution.
Please, can you tell me how implement this in my scripts? I guess the updated function is not already included in a package update.
What should I do? Copy the function in my R console after loading the package?
thanks
You can install the updated version using remotes::install_github("r-spatialecology/landscapemetrics@263-circular-window-in-window_lsm").