sf
sf copied to clipboard
r-spatial
Buffering by negative value returns empty polygon
Buffering by negative value returns empty polygon. I think, but not entirely sure, this was not the case before.
To Reproduce
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
bbox <- st_bbox(c(xmin = 16.1, xmax = 16.6, ymax = 48.6, ymin = 47.9), crs = st_crs(4326))
st_buffer(st_as_sfc(bbox), -5)
#> Geometry set for 1 feature (with 1 geometry empty)
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
#> Geodetic CRS: WGS 84
#> POLYGON EMPTY
st_bbox(st_buffer(st_as_sfc(bbox), -5))
#> xmin ymin xmax ymax
#> NA NA NA NA
st_buffer(st_as_sfc(bbox), 5)
#> Geometry set for 1 feature
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 16.092 ymin: 47.89505 xmax: 16.6081 ymax: 48.60544
#> Geodetic CRS: WGS 84
#> POLYGON ((16.12336 47.89515, 16.12608 47.89983,...
st_bbox(st_buffer(st_as_sfc(bbox), 5))
#> xmin ymin xmax ymax
#> 16.09200 47.89505 16.60810 48.60544
Created on 2022-08-17 with reprex v2.0.2
sessionInfo()
R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.4 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8
[6] LC_MESSAGES=C.UTF-8 LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices datasets utils methods base
other attached packages:
[1] reprex_2.0.2 sf_1.0-8 forcats_0.5.1 stringr_1.4.0 dplyr_1.0.9 purrr_0.3.4 readr_2.1.2
[8] tidyr_1.2.0 tibble_3.1.8 ggplot2_3.3.6 tidyverse_1.3.2
loaded via a namespace (and not attached):
[1] httr_1.4.3 jsonlite_1.8.0 modelr_0.1.8 assertthat_0.2.1 sp_1.5-0
[6] highr_0.9 renv_0.15.5 googlesheets4_1.0.0 cellranger_1.1.0 yaml_2.3.5
[11] pillar_1.8.0 backports_1.4.1 lattice_0.20-45 glue_1.6.2 digest_0.6.29
[16] rvest_1.0.2 colorspace_2.0-3 htmltools_0.5.3 pkgconfig_2.0.3 broom_1.0.0
[21] raster_3.5-29 haven_2.5.0 stars_0.5-6 s2_1.1.0 scales_1.2.0
[26] processx_3.7.0 terra_1.6-7 tzdb_0.3.0 proxy_0.4-27 googledrive_2.0.0
[31] generics_0.1.3 ellipsis_0.3.2 withr_2.5.0 cli_3.3.0 magrittr_2.0.3
[36] crayon_1.5.1 readxl_1.4.0 ps_1.7.1 evaluate_0.16 fs_1.5.2
[41] fansi_1.0.3 xml2_1.3.3 lwgeom_0.2-8 class_7.3-20 tools_4.2.1
[46] hms_1.1.1 gargle_1.2.0 lifecycle_1.0.1 munsell_0.5.0 callr_3.7.1
[51] compiler_4.2.1 e1071_1.7-11 RNetCDF_2.6-1 rlang_1.0.4 classInt_0.4-7
[56] units_0.8-0 grid_4.2.1 rstudioapi_0.13 rmarkdown_2.14 wk_0.6.0
[61] gtable_0.3.0 codetools_0.2-18 abind_1.4-5 DBI_1.1.3 R6_2.5.1
[66] ncmeta_0.3.0 lubridate_1.8.0 knitr_1.39 fastmap_1.1.0 utf8_1.2.2
[71] KernSmooth_2.23-20 stringi_1.7.8 parallel_4.2.1 Rcpp_1.0.9 vctrs_0.4.1
[76] dbplyr_2.2.1 tidyselect_1.1.2 xfun_0.31
What was returned before, and before what or when? If anything this would be a change in GEOS, which maybe returned a GEOMETRYCOLLECTION EMPTY before?
I worked on it a few months ago, but I've changed machines since and unfortunately can't tell anymore what version of sf or GEOS I had installed. If I recall correctly, it returned a POLYGON, but I merely used it as a quick way to get a shrunken bbox, which now doesn't work anymore either.
Regardless, I would expect a negative buffer to shrink polygons, but all polygons become empty if I do now. But maybe my assumption is incorrect. See below:
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(units)
#> udunits database from /usr/share/xml/udunits/udunits2.xml
nc_sf <- st_read(system.file("shape/nc.shp", package="sf"), quiet=TRUE)
plot(nc_sf["AREA"])
st_buffer(nc_sf["AREA"], set_units(0.05, degree))
#> Simple feature collection with 100 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: -84.38574 ymin: 33.82673 xmax: -75.38785 ymax: 36.64189
#> Geodetic CRS: NAD27
#> First 10 features:
#> AREA geometry
#> 1 0.114 MULTIPOLYGON (((-81.21884 3...
#> 2 0.061 MULTIPOLYGON (((-80.86545 3...
#> 3 0.143 MULTIPOLYGON (((-80.7824 36...
#> 4 0.070 MULTIPOLYGON (((-75.80691 3...
#> 5 0.153 MULTIPOLYGON (((-77.21177 3...
#> 6 0.097 MULTIPOLYGON (((-77.19596 3...
#> 7 0.062 MULTIPOLYGON (((-76.15519 3...
#> 8 0.091 MULTIPOLYGON (((-76.44185 3...
#> 9 0.118 MULTIPOLYGON (((-77.92781 3...
#> 10 0.124 MULTIPOLYGON (((-80.50494 3...
plot(st_buffer(nc_sf["AREA"], set_units(0.05, degree)))
st_buffer(nc_sf["AREA"], set_units(-0.05, degree))
#> Simple feature collection with 100 features and 1 field (with 100 geometries empty)
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
#> Geodetic CRS: NAD27
#> First 10 features:
#> AREA geometry
#> 1 0.114 POLYGON EMPTY
#> 2 0.061 POLYGON EMPTY
#> 3 0.143 POLYGON EMPTY
#> 4 0.070 POLYGON EMPTY
#> 5 0.153 POLYGON EMPTY
#> 6 0.097 POLYGON EMPTY
#> 7 0.062 POLYGON EMPTY
#> 8 0.091 POLYGON EMPTY
#> 9 0.118 POLYGON EMPTY
#> 10 0.124 POLYGON EMPTY
plot(st_buffer(nc_sf["AREA"], set_units(-0.05, degree)))
#> Error in plot_sf(x, ...): NA value(s) in bounding box. Trying to plot empty geometries?
Created on 2022-08-17 with reprex v2.0.2
Is the outcome the same if st_use_s2() is FALSE? The input geometries are geographic not planar, so s2 may be in play.
should this issue be addressed with a documentation PR on buffer with a warning about sf_us_s2()
Hello everyone,
I am a new user and my intention basically is to create a convex hull and then create a smaller hull using st_buffer and a negative value of buffer distance. However, like in this and others posts, it always return an empty polygon.
I was checking the Simple Features for R vignette and there is this example on the Geometric operations where the st_buffer works fine with negative values and it got me by surprise. The only difference from my convex hull and the example was the projection (WGS 84 vs WGS 84 / Pseudo-Mercator), so i tried to change the projection of the example and voila, the negative value using WGS 84 returns an empty polygon.
Is that a bug?
library(sf)
nc <- st_read(system.file("shape/nc.shp", package="sf"))
nc.web.EPSG <- st_transform(nc, crs= "EPSG:4326")
sel <- c(1,5,14)
geom <- st_geometry(nc.web.EPSG[sel,])
buf <- st_buffer(geom, dist = 30000)
plot(buf, border = 'red')
plot(geom, add = TRUE)
plot(st_buffer(geom, -10000), add = TRUE, border = 'blue')
#Geometry set for 3 features (with 3 geometries empty)
#Geometry type: POLYGON
#Dimension: XY
#Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
#Geodetic CRS: WGS 84
#POLYGON EMPTY
#POLYGON EMPTY
#POLYGON EMPTY
nc.web_mercator <- st_transform(nc, crs= 3857)
sel <- c(1,5,14)
geom <- st_geometry(nc.web_mercator[sel,])
buf <- st_buffer(geom, dist = 30000)
plot(buf, border = 'red')
plot(geom, add = TRUE)
plot(st_buffer(geom, -10000), add = TRUE, border = 'blue')
#Geometry set for 3 features
#Geometry type: POLYGON
#Dimension: XY
#Bounding box: xmin: -9086682 ymin: 4339742 xmax: -8596499 ymax: 4371623
#Projected CRS: WGS 84 / Pseudo-Mercator
#POLYGON ((-9074360 4347826, -9079727 4354112, -...
#POLYGON ((-8632038 4367097, -8607181 4367131, -...
#POLYGON ((-8782440 4342743, -8788570 4342863, -...
sessionInfo()
R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.1 LTS
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.20.so
locale:
[1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8 LC_COLLATE=pt_BR.UTF-8
[5] LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=pt_BR.UTF-8 LC_PAPER=pt_BR.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] sf_1.0-9
loaded via a namespace (and not attached):
[1] styler_1.8.1 gtools_3.9.3 tidyselect_1.2.0 terra_1.6-47 xfun_0.34 purrr_1.0.0
[7] lattice_0.20-45 vctrs_0.5.1 generics_0.1.3 htmltools_0.5.3 s2_1.1.1 yaml_2.3.6
[13] utf8_1.2.2 rlang_1.0.6 R.oo_1.25.0 e1071_1.7-12 pillar_1.8.1 glue_1.6.2
[19] withr_2.5.0 DBI_1.1.3 R.utils_2.12.2 sp_1.5-1 wk_0.7.0 lifecycle_1.0.3
[25] R.cache_0.16.0 gtable_0.3.1 raster_3.6-11 R.methodsS3_1.8.2 codetools_0.2-18 evaluate_0.18
[31] knitr_1.40 callr_3.7.3 fastmap_1.1.0 ps_1.7.2 class_7.3-20 fansi_1.0.3
[37] highr_0.9 Rcpp_1.0.9 KernSmooth_2.23-20 clipr_0.8.0 classInt_0.4-8 fs_1.5.2
[43] gridExtra_2.3 digest_0.6.30 processx_3.8.0 dplyr_1.0.10 dismo_1.3-9 grid_4.2.1
[49] cli_3.4.1 tools_4.2.1 magrittr_2.0.3 proxy_0.4-27 tibble_3.1.8 pkgconfig_2.0.3
[55] reprex_2.0.2 assertthat_0.2.1 rmarkdown_2.18 rstudioapi_0.14 R6_2.5.1 units_0.8-1
[61] compiler_4.2.1
sf::sf_extSoftVersion()
GEOS GDAL proj.4 GDAL_with_GEOS USE_PROJ_H PROJ
"3.10.2" "3.4.1" "8.2.1" "true" "true" "8.2.1"
Please read the comments above and https://r-spatial.org/r/2020/06/17/s2.html. Does sf_use_s2() make a difference for the non-planar case?
Please read the comments above and https://r-spatial.org/r/2020/06/17/s2.html. Does sf_use_s2() make a difference for the non-planar case?
Using sf_use_s2(FALSE) also returns empty polygons as well.
st_buffer(geom, -10000) dist is assumed to be in decimal degrees (arc_degrees). Geometry set for 3 features (with 3 geometries empty) Geometry type: POLYGON Dimension: XY Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA Geodetic CRS: WGS 84 POLYGON EMPTY POLYGON EMPTY POLYGON EMPTY Warning message: In st_buffer.sfc(geom, -10000) : st_buffer does not correctly buffer longitude/latitude data
Which is as expected for a buffer of -10000 degrees, isn't it?
With sf 1.0-10 (devel) and 1.0-9 (CRAN):
sf_use_s2(FALSE)
plot(st_buffer(geom, -0.1), add = TRUE, border = 'blue')
