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

[misea]

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 tidyverse/ggplot2#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]

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)

Some further discussion in the context of map plotting and tmap: r-tmap/tmap#457

[kadyb]

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)

[edzer]

#> [1] TRUE # it seems that all geometries are fixed after transformation

Well, "fixed"? - Russia is missing, but should have been shown. This is what sf does too, but it does not what you'd want.

[misea]

From my naive point of view as a user of these packages, what I’d like to see is points on the polygons that are partially in view be remapped to the closest in-view point. That’s probably a lot to ask but I think it would work the way people expect. Certainly better than a difficult to predict/debug exception when one calls plot().