Return proj4 CRS object defining a two-point equidistant projection fitting a spatial object

project_tpe.R

require(raster)

project_tpe = function(d, x=0, y=0){
  d = data.frame(t(as.matrix(extent(d))))
  dy = mean(d$y)
  CRS(paste0("+proj=tpeqd +lat_1=", dy, " +lon_1=", d$x[1],
    " +lat_2=", dy, " +lon_2=", d$x[2], " +x_0=", x, " +y_0=", y,
    " +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))
}