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

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.

• Putting this off-topic is really non-sensical. I do not pretend to program this myself but was rather asking for an existing implementation. This question is on software/plattform advice and you have hundreds of this on this GIS forum. Furthermore it touches a very common problem in the field of cartography. This is NOT a dedicated programming forum and If I had a programming issue I would have gone to SO anyways. The response below from TimSalabim is very nice and you are effectively devaluating his effort. Feb 27, 2019 at 15:01

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.

• should be noted here that these projection system work good for large areas. for smaller ones the code would need to be extended a lot. Mar 7, 2019 at 11:00
• Can you elaborate? Mar 7, 2019 at 11:10