Circle buffer and projections issue

I’m struggling to figure out the most appropriate CRS to use for computing spheres of influence around arbitrary placed points on the globe. I haven’t found a resource that recommends CRS for this purpose. I’ve made use of epsg.io which has a huge catalogue of CRS.

At the moment my workflow is experimental, I iterate through several CRS to find the best choice. Which is extremely far from efficient. This RMarkdown document shows my current workflow… the version I’ve posted to RPubs (https://rpubs.com/charliejhadley/circle-crs) has the code hidden for readability, the raw code can be found in this gist - https://gist.github.com/charliejhadley/d738166bee8719345d883d379d4f5587

London

Let’s start by storing the location of 221B Baker Street as a tibble with CRS 4236.

## Create an sf object storing the location of 221B Baker Street with EPSG:4326
baker_street <- tibble(lng = -0.158333, lat = 51.523333) %>% 
  st_as_sf(coords = c("lng", "lat"), crs = 4326)

baker_street %>% 
  mapview()

I’ve created circle_buffer(point, crs_for_buffer, distance) which will use st_buffer() to create a circle centered on point after reprojecting the data into crs_for_buffer

circle_buffer <- function(point,
                          crs_for_buffer,
                          distance,
                          crs_type){
  
  crs_for_buffer <- ifelse(str_starts(crs_for_buffer, "[0-9]"), as.numeric(crs_for_buffer), crs_for_buffer)
  
  if(tolower(crs_type) == "geographic"){
    
    point %>% 
    st_transform(crs_for_buffer) %>% 
    st_buffer(dist = distance / 111320) %>% 
    mutate(label = paste("Geographic projection -", crs_for_buffer))
    
  } else {
    
    point %>% 
    st_transform(crs_for_buffer) %>% 
    st_buffer(dist = distance) %>% 
    mutate(label = paste("Geographic projection -", crs_for_buffer))
    
  }
  
}

We’re going to be working with geographic and projected CRS which have different unit measures, degrees and metres respectively. Let’s use 100e5 in projected CRS and 100e3/111320 in geographic CRS - as detailed in this SO answer https://gis.stackexchange.com/a/142327/78822.

Let’s start off with CRS 4326, as expected we get highly eccentric ellipse instead of a circle:

circle_crs_4326 <- baker_street %>% 
  circle_buffer(4326, 100e3, "geographic")
circle_crs_4326 %>% 
  mapview()

I’m going to test the following CRS:

Circle area within each CRS

The table below shows the computed area for each circle within each CRS:

All four of the projected CRS deviate from Pi by the same amount - 1.3%

Circle after reprojection to CRS 4326

Now we reproject each circle into CRS 4326 and recalculate area:

Visualising all of the circes it appears that it’s a draw between: British National Grid, ED50 and the Europe Albers Equal Area Conic CRS.

White House

Let’s do exactly the same but with the White House.

Here are some CRS I chose:

Circle area within each CRS

Once again we compute the area within each CRS:

All four of the projected CRS deviate from Pi by the same amount - 1.35%

Circle after reprojection to CRS 4326

Now we reproject each circle into CRS 4326 and recalculate area:

The least eccentric circle comes from the Pseudo-Mercator projection, all others are clearly eccentric.