Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Last active March 1, 2026 11:33
Show Gist options
  • Select an option

  • Save mdsumner/045b19365f605492baed2fe85b15afd5 to your computer and use it in GitHub Desktop.

Select an option

Save mdsumner/045b19365f605492baed2fe85b15afd5 to your computer and use it in GitHub Desktop.
library(wk)
library(geos)
library(PROJ)
## flatten all to POLYGON
map <- wk_flatten(as_wkb(rnaturalearth::ne_countries(scale = 50)))
local_laea <- function(x, merctrans) {
  cent <- geos_centroid(x)
  mercxy <- wk_coords(wk_transform(cent, merctrans))[,c("x", "y")]
  xy <- wk_coords(cent)[, c("x", "y")]
  prj <- sprintf("+proj=laea +lon_0=%f +lat_0=%f +datum=WGS84 +type=crs", xy[1], xy[2])
  trans <- proj_trans_create(wk_crs(x), prj)
  
  out <- wk_transform(x, trans)
  wk_crs(out) <- NA
  wk_transform(out, wk_affine_translate(dx = mercxy[["x"]], dy = mercxy[["y"]]))
}

merctrans <- PROJ::proj_trans_create(wk_crs(map), "EPSG:3857")

x <- do.call(c, lapply(map, local_laea, merctrans))

png("laea-vs-mercator.png", width = 1024, height = 1024)
plot(x, border = "transparent")
plot(wk_transform(map, merctrans), col = "lightgrey", border = "transparent", lty = 1, add = TRUE)

## drop anything "too small"
plot(x[geos_area(x) > 5e8], border = "firebrick", add = TRUE, lwd = 1.5)
dev.off()
@mdsumner
Copy link
Author

mdsumner commented Mar 1, 2026

image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment