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)

Hi @rnuske, thank you for testing the script.

Don't worry about the warning, it is a standard print when you manipulate or intersect geometries in lat/lon coordinates. To avoid it, in principle one should:

1. project to some planar CRS
2. do the operation
3. project back to lat/lon.

Unless you are looking after precise results (e.g. compute the centroid of a region with great accuracy), usually it is not worth it.

Concerning the error, I was able to reproduce it.
I have no explanation about why it used to work before, but I have now replaced line 67 with

visible <- st_intersection(st_make_valid(mini_world), st_buffer(circle_longlat, -0.09)) %>%

That should do it.

Thanks for your very fast answer!
With the new line 67 it also works for me.

just a short follow up question: Do you happen to know how to color the oceans?

I am very sorry I missed your message.
To display the color you may uncomment https://gist.github.com/fzenoni/ef23faf6d1ada5e4a91c9ef23b0ba2c1#file-ortho_view-r-L108 and comment https://gist.github.com/fzenoni/ef23faf6d1ada5e4a91c9ef23b0ba2c1#file-ortho_view-r-L109. That will completely hide the graticule, though.

To preserve (and control) the graticule, you need to use the ggplot2's theme() function.

This is inspired by the code I showed in the original Github issue.

ggplot() +
  geom_sf(data = circle,
          fill = 'aliceblue') +
  geom_sf(data=st_collection_extract(visible)) +
  coord_sf(crs = ortho) +
  theme(
    panel.grid.major = element_line(
      color = 'red', linetype = 2, size = 0.5),
    panel.ontop = TRUE,
    panel.background = element_rect(fill = NA)
  )

You might also want to play with the alpha level, so that the graticule goes below the land masses.
It's a quirky solution because one can't perfectly control the color of the sea anymore, nor the color of the graticule.

# Final plot
ggplot() +
  geom_sf(data = circle,
          fill = 'aliceblue', alpha=0.5)+
  geom_sf(data=st_collection_extract(visible), fill='darkgreen', color = 'black') +
  coord_sf(crs = ortho) +
  theme(
    panel.grid.major = element_line(
      color = 'red', linetype = 2, size = 0.5)
  )

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.