search

r-spatial / sf

Simple Features for R
https://r-spatial.github.io/sf/
Other
1.34k stars 298 forks source link

`st_make_grid` produces a non-overlaping row of hexagons. #2218

Closed albhasan closed 1 year ago

albhasan commented 1 year ago

The documentation states that st_make_grid creates a hexagonal grid covering the bounding box of the given sf object. However, the top row of hexagons produced by the code below doesn't overlap the bounding box of the given sf object.

This issue could be related to #1229 but, as I'm not sure, I'd rather report it.

library(sf)
nc <- st_read(system.file("shapes/", package = "maptools"), "sids")
hex_grid <- st_make_grid(nc, cellsize = 1.5, square = FALSE)
plot(hex_grid, bg = NULL, border = "red")
plot(st_as_sfc(st_bbox(nc)), bg = NULL, lwd = 2, add = TRUE)

In this plot, the black square is nc's bounding box and the grid is in red.

This issue happens only when square = FALSE and at specific cellsizes (i.e. cellsizes 0.5, 1.4, 1.5).

> sessionInfo() R version 4.0.4 (2021-02-15) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Debian GNU/Linux 11 (bullseye) Matrix products: default BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.13.so 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 attached base packages: [1] stats graphics grDevices utils datasets methods base loaded via a namespace (and not attached): [1] Rcpp_1.0.11 fansi_1.0.4 class_7.3-22 utf8_1.2.3 [5] sf_1.0-14 dplyr_1.1.2 grid_4.0.4 R6_2.5.1 [9] lifecycle_1.0.3 DBI_1.1.3 magrittr_2.0.3 e1071_1.7-13 [13] units_0.8-3 pillar_1.9.0 KernSmooth_2.23-22 rlang_1.1.1 [17] cli_3.6.1 generics_0.1.3 vctrs_0.6.3 tools_4.0.4 [21] glue_1.6.2 proxy_0.4-27 compiler_4.0.4 pkgconfig_2.0.3 [25] classInt_0.4-9 tidyselect_1.2.0 tibble_3.2.1 > sf::sf_extSoftVersion() GEOS GDAL proj.4 GDAL_with_GEOS USE_PROJ_H "3.9.0" "3.2.2" "7.2.1" "true" "true" PROJ "7.2.1"
edzer commented 1 year ago

Alber, good to hear from you! I really had some idle hope I would never have to return to that code, but alas. I noted that if you give flat_topped = TRUE then there are no hexagons not intersecting with the bounding box, in your example. Maybe that is a hint for finding the solution, who knows.

albhasan commented 1 year ago

It's good to hear from you too!

You're right, flat_topped fixes it. Thanks!

Cheers!