Automatically get an adequate projection system based on a bounding box [closed]

Question

The website http://projectionwizard.org/ offers a very easy way to get an appropriate projection system based on a bounding box and choosing whether projection criterion should be equal-area, equi-distant or a comprimise.

It is based on this publication that defines some principles on how to select an adequate projection system:

Jenny, B., Šavrič, B., Arnold, N. D., Marston, B. E. and Preppernau, C. A. (2017). A guide to selecting map projections for world and hemisphere maps. In: M. Lapaine and E. L. Usery (eds), Choosing a Map Projection, Lecture Notes in Geoinformation and Cartography (pp. 213–228). Berlin, Heidelberg, New York: Springer. Doi: 10.1007/978-3-319-51835-0_9

I would like to know if there is any similar tool for R (or Python) that can achieve the same. So e.g. based on a shapefile in WGS84 find an adequate projection system for the shape with equalarea as the main criterion.

Relevance: to calculate areas in R without having to manually find a projection system. Currently I only know of the function area in the raster package. However, this function does not use a projected coordinate system. From it's description:

Compute the approximate surface area of cells in an unprojected (longitude/latitude) Raster object. It is an approximation because area is computed as the height (latitudial span) of a cell (which is constant among all cells) times the width (longitudinal span) in the (latitudinal) middle of a cell. The width is smaller at the poleward side than at the equator-ward side of a cell. This variation is greatest near the poles and the values are thus not very precise for very high latitudes.

Answer 1

I don't know of existing tools, but this is something we can write ourselves. Here are two functions that will give back a projstring based on the bounding box of the (sf) object that is supplied. One is Lambert Azimutal Equal Area, the other one is Two Point Equidistant.

## area projection (lambert azimutal equal area) with projection parameters 
## derived from feature extent
area_proj <- function(x) {
  bb <- sf::st_bbox(x)
  cntr_long <- (bb[3] - bb[1]) * 0.5 + bb[1]
  cntr_lat <- (bb[4] - bb[2]) * 0.5 + bb[2]
  bbm <- bb * 111000
  rng_x <- round(bbm[3] - bbm[1])
  rng_y <- round(bbm[4] - bbm[2])
  paste0("+proj=laea +lat_0=",
         cntr_lat,
         " +lon_0=",
         cntr_long,
         " +x_0=",
         rng_x,
         " +y_0=", 
         rng_y,
         " +a=6371007.181 +b=6371007.181 +units=m +no_defs")
}


## distance projection (tpeqd - two-point eqdistant) with projection parameters 
## derived from feature extent
dist_proj <- function(x) {
  bb <- sf::st_bbox(x)
  paste0("+proj=tpeqd +lat_1=",
         bb[2],
         " +lon_1=", 
         bb[1], 
         " +lat_2=",
         bb[4], 
         " +lon_2=", 
         bb[3],
         " +x_0=0",
         " +y_0=0",
         " +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
}

If you need this for raster data, you will need to adjust the extraction of the bounding box, e.g. bb <- sf::st_bbox(x) to bb <- raster::extent(x) and then use those in the correct order.