Find closest point while staying in a polygon?

Question

I have two cities and one individual. I'd like to find the city that is closest to this individual. To do that, I can use sf::st_nearest_feature(). However, so far it only uses the distance as the crow flies between the individual and each city. I'd like to add the constraint that the path must stay inside a polygon.

In the example below:

individual (red triangle) is closer to city A than to city B if we consider distance as the crow flies;
however if we add the constraint that the individual can only move inside the polygon, then it is closer to city B.

library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(ggplot2)
library(ggrepel)
library(rnaturalearth)

background <- ne_countries(scale = 'small', type = 'map_units', returnclass = 'sf') |>
  subset(name %in% c("England", "Wales")) |> 
  st_union()

cities <- data.frame(
  name = c("A", "B"),
  lon = c(-4.3, -3.3),
  lat = c(51.2, 51.45)
) |> 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

individual <- data.frame(id = 1, lon = -4.3, lat = 51.6) |> 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

ggplot() +
  geom_sf(data = background) +
  geom_sf(data = cities, size = 3) +
  geom_sf(data = individual, color = "red", shape = 17, size = 3) +
  coord_sf(xlim = c(-6, -1), ylim = c(50, 52)) +
  geom_text_repel(
    data = cities, 
    aes(geometry = geometry, label = name),
    stat = "sf_coordinates",
  )
#> Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
#> give correct results for longitude/latitude data

st_nearest_feature() tells me that the nearest city is A:

nearest <- st_nearest_feature(individual, cities)
cities[nearest, "name"]
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -4.3 ymin: 51.2 xmax: -4.3 ymax: 51.2
#> Geodetic CRS:  WGS 84
#>   name          geometry
#> 1    A POINT (-4.3 51.2)

How can I modify my measure of distance so that the nearest city is B? If possible, the solution should scale well to do this for millions of individuals (the number of cities will stay small) in a reasonable time.

There are a couple of solutions to a similar problem here: https://stackoverflow.com/q/68596244/7547327 — mrhellmann, Mar 31 '23 at 10:33

Robert Hacken · Answer 1 · 2023-03-31T13:21:08.330

You can create LINESTRING representations of the paths and use st_contains() to find which are fully inside the polygon. Then select the shortest one after having disqualified those not inside the polygon by setting them to infinity:

paths <- st_cast(st_union(st_geometry(cities), st_geometry(individual)), 
                 'LINESTRING')
path.len <- st_length(paths)
path.len[!st_contains(background, paths, sparse=F)] <- Inf
cities[which.min(path.len), "name"]
# Simple feature collection with 1 feature and 1 field
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: -3.3 ymin: 51.45 xmax: -3.3 ymax: 51.45
# Geodetic CRS:  WGS 84
#   name           geometry
# 2    B POINT (-3.3 51.45)

For more than one individual

Data:

individual <- data.frame(id = 1:3, 
                         lon = -4.3 + c(.7, 0, 1), 
                         lat = 51.6 + c(-.3, 0, 0)) |> 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

ggplot() +
  geom_sf(data = background) +
  geom_sf(data = cities, size = 3) +
  geom_sf(data = individual, color = "red", shape = 17, size = 3) +
  coord_sf(xlim = c(-6, -1), ylim = c(50, 52)) +
  geom_text_repel(
    data = cities, 
    aes(geometry = geometry, label = name),
    stat = "sf_coordinates",
  ) + 
  geom_text_repel(
    data = individual, 
    aes(geometry = geometry, label = id),
    stat = "sf_coordinates",
  )

Basically do the same succesively for each individual:

nearest <- sapply(seq_along(st_geometry(individual)), \(i) {
  paths <- st_cast(st_union(st_geometry(cities), st_geometry(individual)[i]), 
                   'LINESTRING')
  path.len <- st_length(paths)
  path.len[!st_contains(background, paths, sparse=F)] <- Inf
  which.min(path.len)
})
cities[nearest, 'name']
# Simple feature collection with 3 features and 1 field
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: -4.3 ymin: 51.2 xmax: -3.7 ymax: 51.5
# Geodetic CRS:  WGS 84
#     name          geometry
# 1      A POINT (-4.3 51.2)
# 2      B POINT (-3.7 51.5)
# 2.1    B POINT (-3.7 51.5)

It works fine with 1 or 2 individuals, but if I have 3 (or more) then there's an error: `Error in recycle_common(x, as_s2_geography(y)) : Incompatible lengths: 2, 3` — bretauv, Mar 31 '23 at 12:15
@bretauv I see, you didn't indicate there would be more than one individual. See my updated answer. — Robert Hacken, Mar 31 '23 at 13:21
Not explicitly, I understood *points* as cities. It would be best to create example data general enough to make this unambiguous. — Robert Hacken, Mar 31 '23 at 13:32

Find closest point while staying in a polygon?

1 Answers1

For more than one individual