world-slices.R

library(tidyverse)
library(sf)
library(rnaturalearth)

# World map
world <- ne_countries(scale = "medium") |> 
  filter(adm0_a3 != "ATA") |>
  st_transform(crs = "+proj=merc")

# The center of the map is off the coast of Sardegna
center <- c(10, 40)

# Transform center to Mercator coordinates
center_merc <- st_point(center) |> 
  st_sfc(crs = 4326) |> 
  st_transform(crs = "+proj=merc") |> 
  st_coordinates() |> 
  as.numeric()

# Make straight lines in Mercator projection
make_line_mercator <- function(angle, npts = 200, length = 50000000) {  # length in meters
  x_coords <- seq(
    center_merc[1],
    center_merc[1] + length * cos(angle * pi / 180),
    length.out = npts
  )
  y_coords <- seq(
    center_merc[2],
    center_merc[2] + length * sin(angle * pi / 180),
    length.out = npts
  )
  st_linestring(cbind(x_coords, y_coords))
}

# The C angle in the original is ≈34° and B is ≈56°, so we build the angles
# counterclockwise from C
first_slice <- 34

angles <- c(
  0,
  first_slice,
  90,
  180 - first_slice,
  180,
  180 + first_slice,
  270,
  360 - first_slice
)

# Make the actual lines
lines <- angles |>
  map(make_line_mercator) |>
  st_sfc(crs = "+proj=merc") |>
  st_as_sf()

# Crop lines to world boundaries
lines_cropped <- st_crop(lines, st_bbox(world))

# Find the angles in between each of the main angles
mid_angles <- map2_dbl(
  angles,
  lead(angles, default = angles[1] + 360),
  \(x, y) (x + y) / 2
)

# Get the midpoint coordinates in Mercator projection
# Distance is in meters
midpoint_coord_mercator <- function(angle, distance = 5000000) {
  x <- center_merc[1] + distance * cos(angle * pi / 180)
  y <- center_merc[2] + distance * sin(angle * pi / 180)
  tibble(x = x, y = y, angle = angle)
}

# Make a dataframe of all the midpoints
midpoints <- map(mid_angles, midpoint_coord_mercator) |>
  list_rbind() |>
  mutate(label = c("C", "B", "A", "H", "G", "F", "E", "D")) |>
  st_as_sf(coords = c("x", "y"), crs = "+proj=merc")

# Finaly plot these things!
map_thing <- ggplot() +
  geom_sf(data = world, fill = "grey60", color = "white") +
  geom_sf(data = lines_cropped, color = "red", linewidth = 0.5) +
  geom_sf_text(
    data = midpoints,
    aes(label = label),
    size = 6,
    color = "blue",
    fontface = "bold"
  ) +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

map_thing + 
  labs(title = "Original Mercator slices") +
  coord_sf(crs = "+proj=merc")

map_thing + 
  labs(title = "Slices in the Robinson projection") +
  coord_sf(crs = "+proj=robin")

map_thing + 
  labs(title = "Slices in the Equal Earth projection") +
  coord_sf(crs = "+proj=eqearth")