Created
May 15, 2020 09:01
-
-
Save mstrimas/7634c738d6e34ded3d3cf05b1d9f3b23 to your computer and use it in GitHub Desktop.
Project geometries to CRS centered on longitudes other than 0
This file contains hidden or 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
library(sf) | |
library(rnaturalearth) | |
# transform an sf object, moving the dateline seam to edge of projection | |
# | |
# x: sf containing spatial features | |
# crs: character; proj4 string for new projection | |
# | |
# return: | |
# sf containing spatial features in new projection | |
recenter_sf <- function(x, crs) { | |
stopifnot(inherits(x, c("sf", "sfc"))) | |
# find central meridian | |
center <- as.numeric(stringr::str_extract(crs, "(?<=lon_0=)[-.0-9]+")) | |
if (is.na(center) || length(center) != 1 || !is.numeric(center)) { | |
stop("CRS has no central meridian term (+lon_0).") | |
} | |
# edge is offset from center by 180 degrees | |
edge <- ifelse(center < 0, center + 180, center - 180) | |
# create an very narrow sliver to clip out | |
delta <- 1e-5 | |
clipper <- sf::st_bbox(c(xmin = edge - delta, xmax = edge + delta, | |
ymin = -90, ymax = 90), | |
crs = 4326) | |
clipper <- sf::st_as_sfc(clipper) | |
clipper <- suppressWarnings(smoothr::densify(clipper, max_distance = 1e-3)) | |
clipper <- sf::st_transform(clipper, crs = sf::st_crs(x)) | |
# cut then project | |
x_proj <- sf::st_difference(x, clipper) | |
sf::st_transform(x_proj, crs = crs) | |
} | |
# world map | |
countries <- ne_countries(returnclass = "sf") %>% | |
st_geometry() | |
proj <- "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs" | |
# bad :( | |
countries_proj <- st_transform(countries, crs = proj) | |
ggplot(countries_proj) + | |
geom_sf(color = "transparent", fill = "#e6e6e6") + | |
theme_void() | |
# good :) | |
countries_recenter <- recenter_sf(countries, crs = proj) | |
ggplot(countries_recenter) + | |
geom_sf(color = "transparent", fill = "#e6e6e6") + | |
theme_void() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment