st_wrap_dateline() not working as expected

[pieterprovoost]

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"

[edzer]

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"   

[rsbivand]

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 #1991 help?

[pieterprovoost]

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/")

[pieterprovoost]

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

[jmckenna]

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).

[jonathancallahan]

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" 

[pieterprovoost]

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" 

[edzer]

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

[pieterprovoost]

@edzer Disabling s2 did the trick, thanks!

[kadyb]

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()