othrographic plots in r with sf and ggplot

ortho_plot.R

ortho_plot <- function(poly, lat = 55, lon = 20, lwd = 0.25, bgc = "#bfd7ea") {
  require(sf); require(mapview)
  sf_use_s2(FALSE)
  # Define the orthographic projection
  # Choose lat_0 with -90 <= lat_0 <= 90 and lon_0 with -180 <= lon_0 <= 180
  ortho <- sprintf('+proj=ortho +lat_0=%s +lon_0=%s +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs', lat, lon)
  assign(x = "ortho_crs", value = ortho, envir = .GlobalEnv)
  # Define the polygon that will help you finding the "blade"
  # to split what lies within and without your projection
  # buffer is == Semimajor radius of the ellipsoid axis
  circle <- st_sfc(st_buffer(st_point(x = c(0,0)), dist = 6371000),crs = ortho)
  # Project this polygon in lat-lon
  circle_longlat <- st_transform(circle, crs = 4326)
  if(lat != 0) {
    circle_longlat <- st_boundary(circle_longlat)
    circle_coords <- st_coordinates(circle_longlat)[, c(1,2)]
    circle_coords <- circle_coords[order(circle_coords[, 1]),]
    circle_coords <- circle_coords[!duplicated(circle_coords),]
    # Rebuild line
    circle_longlat <- st_sfc(st_linestring(circle_coords), crs = 4326)
    if(lat > 0) {
      rectangle <- list(rbind(circle_coords,
                              c(X = 180, circle_coords[nrow(circle_coords), 'Y']),
                              c(X = 180, Y = 90),
                              c(X = -180, Y = 90),
                              c(X = -180, circle_coords[1, 'Y']),
                              circle_coords[1, c('X','Y')]))
      rectangle <- st_sfc(st_polygon(rectangle), crs = 4326)
    } else {
      rectangle <- list(rbind(circle_coords,
                              c(X = 180, circle_coords[nrow(circle_coords), 'Y']),
                              c(X = 180, Y = -90),
                              c(X = -180, Y = -90),
                              c(X = -180, circle_coords[1, 'Y']),
                              circle_coords[1, c('X','Y')]))
      rectangle <- st_sfc(st_polygon(rectangle), crs = 4326)
    }
    circle_longlat <- st_union(st_make_valid(circle_longlat), st_make_valid(rectangle))
  }
  # A small negative buffer is necessary to avoid polygons still disappearing in a few pathological cases
  visible <- st_transform(st_intersection(st_make_valid(poly), st_buffer(circle_longlat, -0.09)), crs = ortho)
  # DISCLAIMER: This section is the outcome of trial-and-error and I don't claim it is the best approach 
  # Resulting polygons are often broken and they need to be fixed
  # Get reason why they're broken
  broken_reason <- st_is_valid(visible, reason = TRUE)
  # First fix NA's by decomposing them
  # Remove them from visible for now
  na_visible <- visible[is.na(broken_reason),]
  visible <- visible[!is.na(broken_reason),]
  # Open and close polygons
  na_visible <- st_cast(st_cast(na_visible, 'MULTILINESTRING'), "LINESTRING", do_split = TRUE)
  na_visible$npts <- mapview::npts(st_geometry(na_visible), by_feature = TRUE)
  # Exclude polygons with less than 4 points
  na_visible <- na_visible[na_visible$npts >= 4, ]
  na_visible <- subset(na_visible, select = -c(npts))
  na_visible <- st_cast(na_visible, "POLYGON")
  # Fix other broken polygons
  broken <- which(!st_is_valid(visible))
  for(land in broken) {
    result = tryCatch({
      # visible[land,] <- st_buffer(visible[land,], 0) # Sometimes useful sometimes not
      visible[land,] <- st_collection_extract(st_make_valid(visible[land,]))
    }, error = function(e) {
      visible[land,] <<- st_buffer(visible[land,], 0)
    })
  }
  # Bind together the two tables
  visible <- rbind(visible, na_visible)
  # Final plot
  return({ggplot() +
    geom_sf(data = circle, fill = bgc, size = lwd) +
    geom_sf(data = st_collection_extract(visible), size = lwd) +
    coord_sf(crs = ortho)})
}

mask_to_smooth_poly <- function(x, ...) {
  require(smoothr)
  reg <- rast(x)
  poly <- as.polygons(reg)
  poly <- st_as_sf(as.data.frame(poly, geom = TRUE),
                   wkt = "geometry",
                   crs = crs(poly))
  poly <- smoothr::smooth(poly, ...)
  return(poly)
  }