Skip to content

Instantly share code, notes, and snippets.

@mstrimas
Created May 15, 2020 09:01
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save mstrimas/7634c738d6e34ded3d3cf05b1d9f3b23 to your computer and use it in GitHub Desktop.
Save mstrimas/7634c738d6e34ded3d3cf05b1d9f3b23 to your computer and use it in GitHub Desktop.
Project geometries to CRS centered on longitudes other than 0
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