Skip to content

Instantly share code, notes, and snippets.

@scbrown86
Last active September 21, 2022 02:41
Show Gist options
  • Save scbrown86/5003d8253d5873c293ca68655ba919fe to your computer and use it in GitHub Desktop.
Save scbrown86/5003d8253d5873c293ca68655ba919fe to your computer and use it in GitHub Desktop.
rotated projections
# Small function for correctly plotting a global map whith a rotated projection (i.e. non standard centre longitude)
## from here - https://stackoverflow.com/a/68313482
rotate_prj <- function(x, crs) {
stopifnot(inherits(x, what = "sf"))
stopifnot(inherits(st_crs(crs), "crs"))
# make the data valid before doing anything else!
x <- st_make_valid(x)
# determine the rotated/centre longitude from crs
lon <- sapply(strsplit(as.character(st_crs(crs)[2]), "\n"), trimws)
lon <- lon[which(grepl(pattern = "Longitude of natural origin", x = lon))]
lon <- as.numeric(sapply(strsplit(lon, ","),"[", 2))
if (length(lon) == 0) {
lon <- 0
}
# calculate an offset
offset <- 180 - lon
# make a polygon that covers the centre longitude
polygon <- st_sfc(st_polygon(x = list(rbind(
c(-0.0001 - offset, 90),
c(0 - offset, 90),
c(0 - offset, -90),
c(-0.0001 - offset, -90),
c(-0.0001 - offset, 90)))), crs = 4326)
# trim anything the intersect the centre longitude
x2 <- st_difference(x, polygon)
# project to requested crs
x3 <- st_transform(x2, crs)
return(x3)
}
TEST <- FALSE
if (TEST) {
library(sf)
library(rnaturalearth)
world <- ne_countries(110, returnclass = "sf")
# centre lon
lon <- 160
# rotated equal earth projection
## could also use Robinson etc.
prjCRS <- sprintf("+proj=eqearth +lon_0=%s +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs", lon)
# polygons stretch
plot(st_transform(world, crs = prjCRS)[1], main = "polygons stretch across hemisphere")
# polygons don't stretch
plot(rotate_prj(world, prjCRS)[1], main = "polygons don't stretch across hemisphere")
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment