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.
    – joaoal
    Feb 27, 2019 at 15:01

1 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=",
         " +lon_0=",
         " +x_0=",
         " +y_0=", 
         " +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=",
         " +lon_1=", 
         " +lat_2=",
         " +lon_2=", 
         " +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.
    – joaoal
    Mar 7, 2019 at 11:00
  • Can you elaborate?
    – TimSalabim
    Mar 7, 2019 at 11:10

Not the answer you're looking for? Browse other questions tagged or ask your own question.