Middle-earth map projecton

gistfile1.txt

library(sf)
library(ggplot2)
library(rnaturalearth)
library(units)
sf_use_s2(TRUE)
lon0 <- 174
lat0 <- -41
tilt <-   0
proj <- "+proj=eqearth +lon_0=0 +datum=WGS84 +units=m +no_defs"
rot_xy <- function(m, lon0, lat0, tilt = 0, inverse = FALSE) {
  th <- lat0 * pi / 180
  ta <- tilt * pi / 180
  if (!inverse) {
    lam <- (m[, 1] - lon0) * pi / 180
    phi <- m[, 2] * pi / 180
    x <- cos(phi) * cos(lam); y <- cos(phi) * sin(lam); z <- sin(phi)
    x2 <-  x * cos(th) + z * sin(th); y2 <- y; z2 <- -x * sin(th) + z * cos(th)
    x3 <- x2; y3 <- y2 * cos(ta) - z2 * sin(ta); z3 <- y2 * sin(ta) + z2 * cos(ta)
    m[, 1] <- atan2(y3, x3) * 180 / pi
    m[, 2] <- asin(pmax(-1, pmin(1, z3))) * 180 / pi
  } else {
    lam <- m[, 1] * pi / 180
    phi <- m[, 2] * pi / 180
    x3 <- cos(phi) * cos(lam); y3 <- cos(phi) * sin(lam); z3 <- sin(phi)
    x2 <- x3; y2 <- y3 * cos(ta) + z3 * sin(ta); z2 <- -y3 * sin(ta) + z3 * cos(ta)
    x <- x2 * cos(th) - z2 * sin(th); y <- y2; z <- x2 * sin(th) + z2 * cos(th)
    lo <- atan2(y, x) * 180 / pi + lon0
    m[, 1] <- ((lo + 180) %% 360) - 180
    m[, 2] <- asin(pmax(-1, pmin(1, z))) * 180 / pi
  }
  m
}
rot_rec <- function(g, ...) {
  if (is.matrix(g)) return(rot_xy(g, ...))
  if (is.numeric(g) && is.null(dim(g)))
    return(drop(rot_xy(matrix(g, nrow = 1), ...)))
  if (is.list(g)) return(lapply(g, rot_rec, ...))
  g
}
rotate_sf <- function(x, lon0, lat0, tilt = 0) {
  g2 <- st_sfc(
    lapply(st_geometry(x), function(gi)
      structure(rot_rec(unclass(gi), lon0, lat0, tilt), class = class(gi))),
    crs = 4326
  )
  st_set_geometry(x, g2)
}
seam_line <- st_sfc(
  st_linestring(rot_xy(cbind(180, seq(-89.9, 89.9, 0.25)),
                       lon0, lat0, tilt, inverse = TRUE)),
  crs = 4326
)
seam_buf <- st_buffer(seam_line, 4000, max_cells = 20000)   # ~4 km geodesic
reframe <- function(x) {
  x |>
    st_make_valid() |>
    st_difference(seam_buf) |> 
    st_segmentize(set_units(100, km)) |> 
    rotate_sf(lon0, lat0, tilt) |>
    st_transform(proj)
}
world_rot <- ne_countries(scale = "medium", returnclass = "sf") |> reframe()
lonseq <- seq(-180, 180, 0.5)
latseq <- seq(-90, 90, 0.5)
outline <- st_sf(geometry = st_sfc(st_polygon(list(rbind(
  cbind(-180,          latseq),            # up the western edge
  cbind(lonseq,          90),              # along the top
  cbind( 180,      rev(latseq)),           # down the eastern edge
  cbind(rev(lonseq),    -90),              # along the bottom
  cbind(-180,           -90)               # close
))), crs = 4326)) |>
  st_transform(proj)
grat <- st_sf(geometry = st_sfc(c(
  lapply(seq(-150, 180, 30), function(l)
    st_linestring(cbind(l, seq(-90, 90, 0.5)))),
  lapply(seq(-60, 60, 30), function(b)
    st_linestring(cbind(seq(-180, 180, 0.5), b)))
), crs = 4326)) |>
  st_transform(proj)
ggplot() +
  geom_sf(data = outline,   fill = "aliceblue", colour = "grey40", linewidth = 0.3) +
  geom_sf(data = grat,      colour = "grey75",  linewidth = 0.15) +
  geom_sf(data = world_rot, fill = "grey88",    colour = "grey35", linewidth = 0.15) +
  coord_sf(expand = FALSE) +
  theme_void()