Preserve polygons in orthographic view

ortho_view.R

# Download earth data first
# https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_land.zip

library(sf)
library(lwgeom)
library(dplyr)
library(ggplot2)
library(mapview)

# Read the data
mini_world <- read_sf('data/ne_110m_land/ne_110m_land.shp')

# Define the orthographic projection
# Choose lat_0 with -90 <= lat_0 <= 90 and lon_0 with -180 <= lon_0 <= 180
lat <- 45
lon <- 2
ortho <- paste0('+proj=ortho +lat_0=', lat, ' +lon_0=', lon, ' +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs')

# Define the polygon that will help you finding the "blade"
# to split what lies within and without your projection
circle <- st_point(x = c(0,0)) %>% st_buffer(dist = 6371000) %>% st_sfc(crs = ortho)
  
# Project this polygon in lat-lon
circle_longlat <- circle %>% st_transform(crs = 4326)

# circle_longlat cannot be used as it is
# You must decompose it into a string with ordered longitudes
# Then complete the polygon definition to cover the hemisphere
if(lat != 0) {
    circle_longlat <- st_boundary(circle_longlat)
    
    circle_coords <- st_coordinates(circle_longlat)[, c(1,2)]
    circle_coords <- circle_coords[order(circle_coords[, 1]),]
    circle_coords <- circle_coords[!duplicated(circle_coords),]
    
    # Rebuild line
    circle_longlat <- st_linestring(circle_coords) %>% st_sfc(crs = 4326)
    
    if(lat > 0) {
      rectangle <- list(rbind(circle_coords,
                              c(X = 180, circle_coords[nrow(circle_coords), 'Y']),
                              c(X = 180, Y = 90),
                              c(X = -180, Y = 90),
                              c(X = -180, circle_coords[1, 'Y']),
                              circle_coords[1, c('X','Y')])) %>% 
        st_polygon() %>% st_sfc(crs = 4326)
    } else {
      rectangle <- list(rbind(circle_coords,
                              c(X = 180, circle_coords[nrow(circle_coords), 'Y']),
                              c(X = 180, Y = -90),
                              c(X = -180, Y = -90),
                              c(X = -180, circle_coords[1, 'Y']),
                              circle_coords[1, c('X','Y')])) %>% 
        st_polygon() %>% st_sfc(crs = 4326)
    }
    
    circle_longlat <- st_union(st_make_valid(circle_longlat), st_make_valid(rectangle))
  }
  
# This visualization shows the visible emisphere in red
ggplot() +
    geom_sf(data = mini_world) +
    geom_sf(data = circle_longlat, color = 'red', fill = 'red', alpha = 0.3)
  
# A small negative buffer is necessary to avoid polygons still disappearing in a few pathological cases
# I should not change the shapes too much
  visible <- st_intersection(st_make_valid(mini_world), st_buffer(circle_longlat, -0.09)) %>%
    st_transform(crs = ortho)

# DISCLAIMER: This section is the outcome of trial-and-error and I don't claim it is the best approach 
# Resulting polygons are often broken and they need to be fixed
# Get reason why they're broken
broken_reason <- st_is_valid(visible, reason = TRUE)
  
# First fix NA's by decomposing them
# Remove them from visible for now
na_visible <- visible[is.na(broken_reason),]
visible <- visible[!is.na(broken_reason),]
  
# Open and close polygons
na_visible <- st_cast(na_visible, 'MULTILINESTRING') %>% 
    st_cast('LINESTRING', do_split=TRUE)
na_visible <- na_visible %>% mutate(npts = npts(geometry, by_feature = TRUE))
# Exclude polygons with less than 4 points
na_visible <- na_visible %>%
    filter(npts >=4) %>%
    select(-npts) %>%
    st_cast('POLYGON')

# Fix other broken polygons
broken <- which(!st_is_valid(visible))
for(land in broken) {
    result = tryCatch({
    # visible[land,] <- st_buffer(visible[land,], 0) # Sometimes useful sometimes not
    visible[land,] <- st_make_valid(visible[land,]) %>% 
        st_collection_extract()  
    }, error = function(e) {
      visible[land,] <<- st_buffer(visible[land,], 0)
    })
}

# Bind together the two tables
visible <- rbind(visible, na_visible)

# Final plot
ggplot() +
    geom_sf(data = circle,
        #fill = 'aliceblue') + # if you like the color
        fill = NA) +
    geom_sf(data=st_collection_extract(visible)) +
    coord_sf(crs = ortho)

Thanks a lot for your suggestions @fzenoni. They helped a lot to get what I wanted. I should definitely learn more ggplot!

Did I understand correctly that the graticule can only be behind everything else (in the background) or in front of everything (in the foreground) but not between the layers?

I am glad it worked out for you!

Did I understand correctly that the graticule can only be behind everything else (in the background) or in front of everything (in the foreground) but not between the layers?

After a lot of trial and error, I believe this is correct. But it would be better to ask this question to the ggplot2 maintainers.

This is brilliant! Thanks for sharing!

This is totally brilliant! Thanks a lot! Got this link trough a question a Stack exchange when I tried to solve it on my own.

I haven't really wrapped my head around the effects of using s2 geometry, but I did not get this script to work unless I set sf_use_s2(FALSE).

Thank you @rasmerla, I really appreciate that.

The script was written before sf 1.0, and since then my professional activities shifted a bit, so I haven't had the opportunity yet to adapt it to the new s2 dependency. This vignette gives me the feeling that things got easier though.

Here is an example of ortho view using s2:
https://gist.github.com/rCarto/02dff288244691a53b3c7b55ee44ef86
It's probably not bulletproof though.