sf
sf copied to clipboard
r-spatial
st_wrap_dateline() not working as expected
I'm trying to split a polygon at the dateline, but st_wrap_dateline() is not behaving as expected:
st_point(c(170, 0)) %>%
st_sfc(crs = 4326) %>%
st_buffer(2000000) %>%
st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180")) %>%
plot()
Session info:
> sessionInfo()
R version 4.1.2 (2021-11-01)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.4
Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] glue_1.6.2 mapview_2.11.0 sf_1.0-9 dplyr_1.0.9
loaded via a namespace (and not attached):
[1] Rcpp_1.0.9 pillar_1.8.0 compiler_4.1.2 base64enc_0.1-3 class_7.3-20 tools_4.1.2
[7] digest_0.6.29 jsonlite_1.8.0 satellite_1.0.4 lifecycle_1.0.1 tibble_3.1.8 lattice_0.20-45
[13] png_0.1-7 pkgconfig_2.0.3 rlang_1.0.4 DBI_1.1.3 cli_3.3.0 rstudioapi_0.13
[19] crosstalk_1.2.0 yaml_2.3.5 fastmap_1.1.0 terra_1.6-7 e1071_1.7-11 s2_1.1.0
[25] raster_3.5-21 leaflet.providers_1.9.0 generics_0.1.3 vctrs_0.4.1 htmlwidgets_1.5.4 webshot_0.5.3
[31] stats4_4.1.2 classInt_0.4-7 leaflet_2.1.1 grid_4.1.2 tidyselect_1.1.2 R6_2.5.1
[37] fansi_1.0.3 sp_1.5-0 purrr_0.3.4 magrittr_2.0.3 scales_1.2.0 codetools_0.2-18
[43] htmltools_0.5.3 units_0.8-0 assertthat_0.2.1 colorspace_2.0-3 utf8_1.2.2 KernSmooth_2.23-20
[49] proxy_0.4-27 wk_0.6.0 munsell_0.5.0 leafem_0.2.0 crayon_1.5.1
>
> sf::sf_extSoftVersion()
GEOS GDAL proj.4 GDAL_with_GEOS USE_PROJ_H PROJ
"3.11.0" "3.5.2" "9.1.0" "true" "true" "9.1.0"
Interestingly I get
with
> 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=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_GB.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] sf_1.0-8 nvimcom_0.9-92
loaded via a namespace (and not attached):
[1] Rcpp_1.0.9 magrittr_2.0.3 units_0.8-0 tidyselect_1.1.2
[5] R6_2.5.1 rlang_1.0.5 fansi_1.0.3 s2_1.1.0
[9] dplyr_1.0.9 wk_0.6.0 tools_4.2.1 grid_4.2.1
[13] KernSmooth_2.23-20 utf8_1.2.2 cli_3.4.0.9000 e1071_1.7-11
[17] DBI_1.1.3 class_7.3-20 assertthat_0.2.1 tibble_3.1.8
[21] lifecycle_1.0.2 purrr_0.3.4 vctrs_0.4.1 glue_1.6.2
[25] proxy_0.4-27 compiler_4.2.1 pillar_1.8.1 generics_0.1.3
[29] classInt_0.4-7 pkgconfig_2.0.3
> sf_extSoftVersion()
GEOS GDAL proj.4 GDAL_with_GEOS USE_PROJ_H
"3.10.2" "3.4.3" "8.2.1" "true" "true"
PROJ
"8.2.1"
How did you install sf and GDAL? The GDAL version is not that in the CRAN sf binary (3.4.2), but the error is also present in the CRAN sf binary package. It is also present in built from source GDAL 3.5.2 (Fedora). Might https://github.com/r-spatial/sf/pull/1991 help?
GDAL was installed with Homebrew, sf was installed from GitHub:
remotes::install_github("r-spatial/sf", configure.args = "--with-proj-include=/opt/homebrew/include/ --with-proj-lib=/opt/homebrew/lib/", configure.vars = "GDAL_DATA=/opt/homebrew/opt/gdal/share/gdal/")
After experimenting a bit with ogr2ogr it looks like this is a GDAL 3.5.2 issue, with 3.4.3 it wraps fine.
conda install -c conda-forge gdal=3.4.3 # or 3.5.2
ogr2ogr -wrapdateline -nlt multipolygon -datelineoffset 20 out.gml polygon.gml
For the record, with GDAL-master it also wraps fine, with that same file (it seems that this will be fixed in the next GDAL release).
FYI -- I am getting the same horizontal-lines/unsplit-polygon issue with R 4.2.2 and updated packages:
> sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Monterey 12.6
Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] sf_1.0-8
loaded via a namespace (and not attached):
[1] Rcpp_1.0.9 rstudioapi_0.14 magrittr_2.0.3 units_0.8-0
[5] tidyselect_1.2.0 R6_2.5.1 rlang_1.0.6 fansi_1.0.3
[9] dplyr_1.0.10 tools_4.2.2 grid_4.2.2 KernSmooth_2.23-20
[13] utf8_1.2.2 cli_3.4.1 e1071_1.7-12 DBI_1.1.3
[17] class_7.3-20 tibble_3.1.8 lifecycle_1.0.3 vctrs_0.5.0
[21] glue_1.6.2 proxy_0.4-27 compiler_4.2.2 pillar_1.8.1
[25] generics_0.1.3 classInt_0.4-8 pkgconfig_2.0.3
>
>
> sf::sf_extSoftVersion()
GEOS GDAL proj.4 GDAL_with_GEOS USE_PROJ_H PROJ
"3.11.0" "3.5.3" "9.1.0" "true" "true" "9.1.0"
I upgraded GDAL to the newest 3.6.2, but the issue remains.
> sf::sf_extSoftVersion()
GEOS GDAL proj.4 GDAL_with_GEOS USE_PROJ_H PROJ
"3.11.1" "3.6.2" "9.1.1" "true" "true" "9.1.1"
With
library(sf)
# Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE
sf_extSoftVersion()
# GEOS GDAL proj.4 GDAL_with_GEOS USE_PROJ_H
# "3.11.1" "3.6.2" "9.1.1" "true" "true"
# PROJ
# "9.1.1"
sf_use_s2(FALSE) # world as flat, 2-D!
# Spherical geometry (s2) switched off
st_point(c(170, 0)) |>
st_sfc(crs = 4326) |>
st_buffer(units::set_units(15, degrees)) |> # degrees!
st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180")) |>
plot()
# Warning message:
# In st_buffer.sfc(st_sfc(st_point(c(170, 0)), crs = 4326), units::set_units(15, :
# st_buffer does not correctly buffer longitude/latitude data
I get
Note that
- s2 was switched off
- buffer distance was expressed in degrees
This problem has been fixed in GDAL 3.8 (with s2 enabled).
library("sf")
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
st_point(c(170, 0))|>
st_sfc(crs = 4326) |>
st_buffer(2000000) |>
st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180")) |>
plot()
