Objective: I want to make "globe" maps similar to the attached figure with an automatised work flow (no to little manual tweaking).
Problem: Polygons that are situated outside of the displayed area of the globe get destroyed beyond repair when I transform the vector world map to my desired projection (using the sf-package and the st_transform function). In my attempts to construct the map above, USA became EMPTY until I manually split it into parts and filtered out all colonies leaving only the US mainland and Alaska.
Question: Is there any way to crop polygons that can pose a problem before transforming? I imagine such an operation would need some kind of algorithm that defines the visible area based on the projection, by which polygons could be cropped or removed.
Or are there any other possible solutions that I am missing?
Related previous posts: This question touches on the subject in this post: st-transform-reprojection-issue-with-r-and-sf, but it is not the same problem as all polygons are showing in that case. In addition I do not find the proposed solution completely satisfying as it a) does not explain the origin of the problem, b) why the solution works, and c) demands a change of the centre of the projection that might not be desirable or possible in all cases.
Reproducible example:
library(tmap)
library(sf)
library(spData)
library(dplyr)
globe_crs <- st_crs("+proj=ortho +lat_0=51.470129 +lon_0=-0.452751 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs")
qtm(world) # Hello world
# Transform to the "globe" projection
world_globe <- world %>%
st_transform(crs = globe_crs)
# Make columns indicating the status of different polygons:
world_globe <- world_globe %>% st_make_valid() %>%
mutate(empty = st_is_empty(.),
valid = st_is_valid(.))
# What is the status of the different polygons?
table(world_globe$valid) # All polygons are valid ...
table(world_globe$empty) # ... but quite a few are empty
# Most empty polygons are located in the southern hemisphere
world_globe %>% dplyr::filter(empty == T) %>% pull(name_long)
# Investigate the transformed data:
qtm(world_globe) # Does not work due to empty polygons.
# The plot prints when empty polygons are filtered out:
world_globe %>%
dplyr::filter(empty == F) %>%
qtm()
# (The USA has disappeared)
To investigate if there were any problems related to the (relative) complexity of the world map I also tried to transform a bounding box to my desired projection. However, this result was even worse as it resulted in an invalid polygon.
# Transformation of world bounding box generates invalid polygon:
world %>% st_bbox() %>% st_as_sfc() %>%
st_transform(crs = globe_crs) %>%
st_is_valid()
I believe this result highlights the generality of the problem, as transforming other type of map data could result in even more severe problems.
