local_proj.md

Apply a local projection to every feature in a longlat data set, and plot.

https://gis.stackexchange.com/questions/335570/row-wise-reprojection-of-spdf-in-r

NOTE: there's a grander problem here in that the centroid used is in longlat, so it's not really the centre in the final map - but probably good enough. There are trig methods to get an angular centre that I don't bother with here.

  local_proj <- function(lonlat) {
  sprintf("+proj=laea +lon_0=%f +lat_0=%f +datum=WGS84", lonlat[1], lonlat[2])
}

local_reproj <- function(x) {
  cc <- sf::st_coordinates(st_centroid(x))
  sf::st_transform(x, local_proj(cc))
}
library(sf)
#> Linking to GEOS 3.7.1, GDAL 2.4.2, PROJ 5.2.0
w <- rnaturalearth::ne_countries(returnclass = "sf")


## easiest way to row-wise is split
listof <- lapply(split(w, 1:nrow(w)), local_reproj)

par(mfrow = n2mfrow(length(listof)), mar = rep(0, 4))

## use walk to ignore all the side effect cruft
purrr::walk(listof, ~plot(st_geometry(.x)))

^{Created on 2019-09-17 by the reprex package (v0.3.0)}

do this equal circles thing, something to update tissot with (and compare how that's doing it)

local_proj <- function(lonlat) {
  sprintf("+proj=laea +lon_0=%f +lat_0=%f +datum=WGS84", lonlat[1], lonlat[2])
}

local_rebuff <- function(x) {
  cc <- sf::st_coordinates(st_centroid(x))
  xx <- sf::st_buffer(sf::st_transform(x, local_proj(cc)), 2e5)
  sf::st_transform(xx, sf::st_crs(x))
}
library(sf)
#> Linking to GEOS 3.7.1, GDAL 2.4.2, PROJ 5.2.0
#w <- rnaturalearth::ne_countries(returnclass = "sf")

w <- sfheaders::sf_cast(sf::st_sf(g = sf::st_sfc(sf::st_multipoint(raster::coordinates(raster::raster(raster::extent(-170, 170, -80, 80), res = 15))))), "POINT")
w <-   sf::st_set_crs(w, 4326)                      

## easiest way to row-wise is split
listof <- lapply(split(w, 1:nrow(w)), local_rebuff)
plot(sf::st_transform(do.call(rbind, listof), "+proj=laea"))

better version

local_proj <- function(lonlat) {
  sprintf("+proj=laea +lon_0=%f +lat_0=%f +datum=WGS84", lonlat[1], lonlat[2])
}

local_rebuff <- function(x) {
  cc <- sf::st_coordinates(st_centroid(x))
  xx <- sf::st_buffer(sf::st_transform(x, local_proj(cc)), 3e5)
  sf::st_transform(xx, sf::st_crs(x))
}
library(sf)
sf_use_s2(FALSE)
w <- sfheaders::sf_cast(sf::st_sf(g = sf::st_sfc(sf::st_multipoint(
  as.matrix(expand.grid(seq(-170, 170, by = 20), seq(-80, 80, by = 20)))
)), crs = "OGC:CRS84"), "POINT")


x <- do.call(rbind, lapply(split(w, 1:nrow(w)), local_rebuff))
p1 <- sf::st_transform(x, "+proj=laea")
p2 <- sf::st_transform(x, "+proj=eck4")
op <- par(mfrow = c(2, 1), mar = rep(0, 4))

plot(st_geometry(p1))
plot(st_geometry(p2))
par(op)

much better: https://github.com/hypertidy/earthcircle