search

r-spatial / sf

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

Render fails with specific lon_0 value in st_transform (fails lon_0=27.0, works lon_0=27.1) #2477

Open misea opened 2 days ago

misea commented 2 days ago

Plotting fails when drawing globe at very specific rotation.

To Reproduce

library(rnaturalearth)
library(sf)
#> Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

world <- ne_countries(scale = "medium", returnclass = "sf") |> st_as_sfc()

p1 <- world |> 
  st_transform(crs = st_crs("+proj=ortho +lat_0=0 +lon_0=26.9")) |>
  plot()


p2 <- world |> 
  st_transform(crs = st_crs("+proj=ortho +lat_0=0 +lon_0=27.0")) |>
  plot()

#> Error in polypath(p_bind(L), border = border[i], lty = lty[i], lwd = lwd[i], : Invalid graphics path

Created on 2024-11-11 with reprex v2.1.1

Additional context First saw and reported as occurring via coord_sf https://github.com/tidyverse/ggplot2/issues/6180 , but problem is lower down the call chain

R version 4.4.2 (2024-10-31) Platform: aarch64-apple-darwin20 Running under: macOS Sequoia 15.1 Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0 locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 time zone: America/Los_Angeles tzcode source: internal attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] sf_1.0-18 rnaturalearth_1.0.1 loaded via a namespace (and not attached): [1] jsonlite_1.8.9 highr_0.11 compiler_4.4.2 reprex_2.1.1 Rcpp_1.0.13-1 clipr_0.8.0 callr_3.7.6 yaml_2.3.10 [9] fastmap_1.2.0 R6_2.5.1 classInt_0.4-10 knitr_1.48 tibble_3.2.1 units_0.8-5 DBI_1.2.3 pillar_1.9.0 [17] rlang_1.1.4 utf8_1.2.4 terra_1.7-78 xfun_0.49 fs_1.6.5 cli_3.6.3 withr_3.0.2 magrittr_2.0.3 [25] ps_1.8.1 class_7.3-22 digest_0.6.37 grid_4.4.2 processx_3.8.4 rstudioapi_0.17.1 lifecycle_1.0.4 vctrs_0.6.5 [33] KernSmooth_2.23-24 data.table_1.16.2 proxy_0.4-27 evaluate_1.0.1 glue_1.8.0 codetools_0.2-20 rnaturalearthdata_1.0.0 fansi_1.0.6 [41] e1071_1.7-16 rmarkdown_2.29 httr_1.4.7 tools_4.4.2 pkgconfig_2.0.3 htmltools_0.5.8.1 > 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"
edzer commented 2 days ago

This is well known: you expect too much from the software, at least at the moment. What happens is that st_transform assigns NA values to coordinates of polygons not "visible" on the orthographic projection, and the plotting render mechanism simply drops such polygons or raises an error. I don't see how sf::plot() or st_transform() could be modified to make this work (but am open to suggestions). A single function could do this, and is available e.g. in pkg https://github.com/paleolimbot/s2plot. A script (modified) from here that could form the basis for such a function is:

library(sf)
old <- options(s2_oriented = TRUE) # don't change orientation from here on
countries <- s2::s2_data_countries() |> st_as_sfc()
globe <- st_as_sfc("POLYGON FULL", crs = st_crs(countries))
oceans <- st_difference(globe, st_union(countries))
visible <- st_buffer(st_as_sfc("POINT(-30 -10)", crs = st_crs(countries)), 9800000) # visible half
visible_ocean <- st_intersection(visible, oceans)
visible_countries <- st_intersection(visible, countries)
st_transform(visible_ocean, "+proj=ortho +lat_0=-10 +lon_0=-30") |>
    plot(col = 'lightblue')
st_transform(visible_countries, "+proj=ortho +lat_0=-10 +lon_0=-30") |>
    plot(col = NA, add = TRUE)
options(old)

x

Some further discussion in the context of map plotting and tmap: https://github.com/r-tmap/tmap/issues/457

kadyb commented 1 day ago

I compared this example using terra and it basically works (probably doesn't draw that one missing geometry?).

library(terra)

f = "https://naciscdn.org/naturalearth/50m/cultural/ne_50m_admin_0_countries.zip"
v = vect(paste0("/vsizip/vsicurl/", f), what = "geoms")
v = project(v, "+proj=ortho +lat_0=0 +lon_0=27.0")
#> 1: Point outside of projection domain (GDAL error 1)
#> ...
all(is.valid(v))
#> [1] TRUE # it seems that all geometries are fixed after transformation
plot(v)

Rplot