globe_without_country_distortions.R

#globe map without countries distorted by polygon points being clipped by the "edge of the world"
#This is my hacking around of the fuller script
#https://github.com/r-spatial/sf/blob/master/demo/twitter.R
#linked to from section 8.1 of https://keen-swartz-3146c4.netlify.com/plotting.html

library(sf)
library(maptools)
library(animation)

data(wrld_simpl)
w <- st_as_sf(wrld_simpl)
set.seed(131)
w$f = factor(sample(1:12, nrow(w), replace = TRUE))
circ = function(l = c(-180:180), lon0 = 0, lat0 = 30) {
  deg2rad = pi / 180
  lat = atan(-cos((l - lon0) * deg2rad)/tan(lat0 * deg2rad)) / deg2rad
  xy = if (lat0 == 0) {
    l1 = lon0 - 90
    l2 = lon0 + 90
    rbind(c(l1,-90), c(l2,-90), c(l2,0), c(l2,90), c(l1,90), c(l1,0), c(l1,-90))
  } else if (lat0 > 0) {
    xy = cbind(lon = l, lat = lat)
    rbind(c(-180,90),xy,c(180,90),c(-180,90))
  } else {
    xy = cbind(lon = l, lat = lat)[length(l):1,]
    rbind(c(180,-90), xy, c(-180,-90),c(180,-90))
  }
  st_sfc(st_polygon(list(xy)), crs = st_crs(4326))
}

m = st_make_grid()
m = st_segmentize(m, 4e5)

lat=-16
lon=-160
p4s=paste0("+proj=ortho +lat_0=", lat, " +lon_0=", lon)
plot(st_transform(m, st_crs(p4s), check = TRUE), col = 'lightblue', border = 'lightblue')
crc = circ(lat0 = lat, lon0 = lon)
w0 = suppressWarnings(st_intersection(w, crc))
w0 = st_cast(w0, "MULTIPOLYGON")
fills <- c(rep("darkgreen",2), "#EDC9AF", rep("darkgreen",34), "white", rep("darkgreen",23))
plot(st_transform(w0["f"], st_crs(p4s), check = TRUE), add = TRUE, col=fills, border=NA)