Plot an ellipse between two geographical foci in R

Question

I am trying to find out the geographical area that is equidistant from two points, and to plot this as an ellipse.

I can produce plots for one point easily using st_buffer, and can find numerous R functions that will plot ellipse from a known centroid if I define the axis, but have not been able to find one that will plot an ellipse given two known foci and a defined distance.

The similar question here gets some way towards an answer, but is not readily applicable to geographic situations - Draw an ellipse based on its foci

My code is pretty simple at the moment, and given each coordinate with a 100km radius. However, I would like to find out all the positions that would be reachable by a 200km (or other defined distance) trip between both sites.

library(tidyverse)
library(sf)

#Give Coordinates
citylocations <- tibble::tribble(
  ~city,     ~lon,      ~lat,
  "London",  -0.1276, 51.5072,
  "Birmingham", -1.8904, 52.4862,
)

citydflocations <- as.data.frame(citylocations)

#Convert to SF
citysflocations <- sf::st_as_sf(citydflocations, coords = c("lon","lat" ), crs = 4326)

#Convert location file to National Grid Planar
cityBNGsflocations <- citysflocations %>%
  st_transform(citysflocations, crs = 27700)

#Produce circles with 100km buffer
dat_circles <- st_buffer(cityBNGsflocations, dist = 100000)
join_circles <- st_union(dat_circles) %>%
  st_transform(4326)
plot(join_circles, col = 'lightblue')```

Welcome to Stackoverflow. You might have a look at `nngeo::st_ellipse`. — Chris, Jun 02 '22 at 18:54
Thanks a lot, I've looked at this function, but it requires a centroid location and axes rather than foci and sum of distance from foci, so it doesn't seem to produce what I'm looking for (without some other transformation of the data). — RyMc, Jun 05 '22 at 09:36

score 1 · Accepted Answer · answered Jun 03 '22 at 11:13

The function below should create buffers of varying distances for each of the two points it is given, finds the intersection the two buffers, unions the intersections, and finally returns a convex hull of those intersections. The output should be a near approximation of an ellipse with the two points as foci.

The straight-line(s) distance from one city to any edge of the polygon and then to the other city should equal the distance given in the function (200,000m in the example below).

It works on the data provided, but is fragile as there's no error checking or warning suppression. Make sure the dist argument is greater than the distance between the two points, and that the points have a crs that can use meters as a distance. (lat/lon might not work)

The example below only uses 20 points for the 'ellipse', but changing the function should be relatively straightforward.

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(tidyverse)


#Give Coordinates
citylocations <- tibble::tribble(
  ~city,     ~lon,      ~lat,
  "London",  -0.1276, 51.5072,
  "Birmingham", -1.8904, 52.4862,
)

citydflocations <- as.data.frame(citylocations)

#Convert to SF
citysflocations <- sf::st_as_sf(citydflocations, coords = c("lon","lat" ), crs = 4326)

#Convert location file to National Grid Planar
cityBNGsflocations <- citysflocations %>%
  st_transform(citysflocations, crs = 27700)

#Produce circles with 100km buffer
dat_circles <- st_buffer(cityBNGsflocations, dist = 100000)
join_circles <- st_union(dat_circles) %>%
  st_transform(4326)
#plot(join_circles, col = 'lightblue')

### the ellipse function using 20 buffers ####
ellipse_fn <- function(x_sf, y_sf, distance){
  #set distance argument to meters, get sequence of distances for buffers
  distance = units::set_units(distance, 'm')
  dists_1 <- seq(units::set_units(0, 'm'), distance, length.out = 22)
  
  # create empty sf object to place for loop objects in 
  # purrr would probably be better here
  nrows <- 20
  df <- st_sf(city = rep(NA, nrows), city.1 = rep(NA, nrows), geometry = st_sfc(lapply(1:nrows, function(x) st_geometrycollection())))
  intersections <- for(i in 2:21){
    buff_1 <- st_buffer(cityBNGsflocations[1,], dist = dists_1[i])
    buff_2 <- st_buffer(cityBNGsflocations[2,], dist = distance - dists_1[i])
    intersection <- st_intersection(buff_1, buff_2)
    df[i-1,] <- intersection

  }
  df %>% 
    st_set_crs(st_crs(x_sf)) %>% 
    st_union() %>%
    st_convex_hull()
}
### end ellipse function ###

# Using the ellipse function with 2 points & 200000m distance
ellipse_sf <- ellipse_fn(cityBNGsflocations[1,], cityBNGsflocations[2,], dist = 200000)
# You'll get lots of warnings here about spatial constance...

ggplot() +
  geom_sf(data = ellipse_sf, fill = 'black', alpha = .2) +
  geom_sf(data = cityBNGsflocations, color = 'red')

^{Created on 2022-06-03 by the reprex package (v2.0.1)}

mapview plot of the cities & 'ellipse' on a map:

Thank you very much, this is exactly the type of solution I was looking for, and is certainly beyond my limited skills with R at the moment. It's really appreciated. — RyMc, Jun 05 '22 at 09:38

Plot an ellipse between two geographical foci in R

1 Answers1