Last active
October 29, 2021 05:01
-
-
Save scbrown86/1eca941b6193765889adbbf6dd4e8bd1 to your computer and use it in GitHub Desktop.
othrographic plots in r with sf and ggplot
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment