spin2.R

library(libs2)

# lifted from ggstereo (probably a better way to replicate this in s2)
# this is essentially generating a "cap" but with somewhat irregular
# coordinate spacing
great_circle <- function(x, y, z, length.out = 200) {
  stopifnot(length(x) == 1, length(y) == 1, length(z) == 1)

  # dot-product = 0 for orthogonality
  # VxU = x v1 + y v2 + z v3 = 0

  theta_out <- seq(-pi, pi, length.out = length.out)

  if(z == 0 && y != 0) {
    # solve for v2
    # v2 = (-x * v1 - z * v3) / y
    v1 <- cos(theta_out)
    v3 <- sin(theta_out)
    v2 <- (-x * v1 - z * v3) / y
  } else if(z == 0 && y == 0) {
    # solve for v1
    # v1 <- (-y * v2 - z * v3) / x
    v2 <- cos(theta_out)
    v3 <- sin(theta_out)
    v1 <- (-y * v2 - z * v3) / x
  } else if(z != 0) {
    # solve for v3
    # v3 = (-x v1 - y v2) / z
    v1 <- cos(theta_out)
    v2 <- sin(theta_out)
    v3 <- (-x * v1 - y * v2) / z
  } else {
    stop("Zero-length vector")
  }

  r <- sqrt(v1*v1 + v2*v2 + v3*v3)
  data.frame(x = v1 / r, y = v2 / r, z = v3 / r)
}

# lifted from mdsumner's silicate
# used to detect CCW vs CW ring direction for orthographic crop
tri_area <- function(x, y) {
  colSums(matrix((x[c(3, 1, 2)] + x[1:3]) * (y[c(3, 1, 2)] - y[1:3]), nrow = 3L)) / 2
}

make_hemisphere_latlng <- function(lon, lat, detail, epsilon, precision) {

  ref_latlng <- s2latlng(lat, lon)
  ref_point <- as.data.frame(s2point(ref_latlng))
  great_circle_df <- great_circle(ref_point$x, ref_point$y, ref_point$z, length.out = detail)

  # move these all a tiny bit in the direction of ref_point to prevent
  # invalid geometries at the edges
  great_circle_df <- great_circle_df + (ref_point * epsilon)[rep(1, nrow(great_circle_df)), ]
  lengths <- with(great_circle_df, sqrt(x * x + y * y + z * z))
  great_circle_df <- great_circle_df / data.frame(x = lengths, y = lengths, z = lengths)

  great_circle_point <- s2point(as.matrix(great_circle_df))
  great_circle_latlng <- unique(round(as.data.frame(s2latlng(great_circle_point)), precision))

  # close the ring
  coord_df <- rbind(great_circle_latlng, great_circle_latlng[1, , drop = FALSE])

  # make sure the ring is defined counter clockwise from the perspective of
  # lat, lon
  coord_df_check <- coord_df[c(1, 10, 20), ]
  coord_df_check_proj <- mapproj::mapproject(
    coord_df_check$lng,
    coord_df_check$lat, projection = "orthographic",
    orientation = c(lat, lon, 0)
  )
  if(tri_area(coord_df_check_proj$x, coord_df_check_proj$y) > 0) {
    coord_df[rev(1:nrow(coord_df)), ]
  } else {
    coord_df
  }
}

make_hemisphere_wkt <- function(lon = 0, lat = 0, detail = 10000, epsilon = 0.1, precision = 2) {
  coord_df <- make_hemisphere_latlng(lon, lat, detail = detail, epsilon = epsilon, precision = precision)
  coords <- paste0(coord_df$lng, " ", coord_df$lat, collapse = ", ")
  sprintf("POLYGON ((%s))", coords)
}

plot_hemisphere <- function(geog, lon = 0, lat = 0, rotation = 0) {
  hemisphere <- s2geography(make_hemisphere_wkt(lon, lat))
  geog_hemisphere <- s2_intersection(geog, hemisphere)
  geog_hemisphere_wkb <- wk::new_wk_wkb(s2_asbinary(geog_hemisphere))
  geog_hemisphere_coords <- wk::wkb_coords(geog_hemisphere_wkb, sep_na = TRUE)
  geog_hemisphere_proj <- mapproj::mapproject(
    geog_hemisphere_coords$x,
    geog_hemisphere_coords$y,
    projection = "orthographic",
    orientation = c(lat, lon, rotation)
  )

  withr::with_par(list(mai = c(0, 0, 0, 0), omi = c(0, 0, 0, 0)), {
    plot(
      double(), double(),
      xlim = c(-1, 1),
      ylim = c(-1, 1),
      asp = 1
    )

    polypath(geog_hemisphere_proj, col = "grey90")
  })
}

countries <- s2data_countries()

ragg::agg_png(width = 600, height = 600)
for (lng in seq(0, 360, length.out = 100)[-1]) {
  try(plot_hemisphere(countries, lng, -70))
}
dev.off()

gifski::gifski(list.files(pattern = "Rplot"), width = 600, height  = 600, loop = TRUE, delay = 1/25)
unlink(list.files(pattern = "Rplot"))