Find minimum number of polygons connecting two points

Question

Given GIS data of polygons and two points, is there a way to find the minimum number of polygons you need to traverse to walk between the two points? Even an approximate approach would be fine for my purposes if the problem cannot always be solved exactly.

The actual application I am using this for uses tax parcels and pipelines, with the goal of finding minimum number of parcels I need to traverse to lay a pipeline between two points -- so thousands of polygons and thousands of pairs of points. But to provide an example here is some code in R using the sf package. I load North Carolina counties and put two points on it. The straight line path between the two points goes through 5 counties. But if you look at the map you can see there's a path that goes through only 4 counties. I want an approach to calculate that 4 in this case, and more generally.

# install.packages(c('sf','dplyr')) # if needed
library(sf)
library(dplyr)

# use NC counties polygons
polys <- read_sf(system.file("shape/nc.shp", package = "sf"))

# make two points, goal is to find the minimum # of counties to get between them
points <- data.frame(name=c('start','end'),LAT=c(34.623858,35.257017),LONG=c(-78.851954,-78.237195))
points <- st_as_sf(points, coords = c("LONG", "LAT"), crs = 4326)

plot(polys$geometry)
plot(points$geometry, add=TRUE,col='blue')
plot(st_cast(summarize(points,name=first(name)),'LINESTRING'),add=TRUE,col='blue')

Answer 1

Helpful comments to my question led me to study minimum spanning trees, which aren't exactly the solution -- spanning trees connect all nodes in the minimum cost way, while I want to connect two given nodes in the minimum cost way. But the insight to consider the data as a graph carries over and leads to a simple solution. Once the GIS data is converted to an adjacency graph there are ready algorithms to find the shortest path between two polygons along the graph. (For example Dijkstra's algorithm, although for my particular problem which has no cost weights it appears that R's igraph package uses a faster algorithm.)

Steps for my solution are:

1. take the polygon dataset and convert it to a graph using adjacency as the definition of edge, and polygon as the definition of vertex.

2. Given the two points (start and end) find the vertices that contain them.

3. Then apply a shortest path algorithm for the two containing vertices.

R code and packages implementing this for the North Carolina example below.

# install.packages(c('sf','dplyr','spdep','igraph')) # if needed
library(sf)
library(dplyr)
library(spdep)
library(igraph)

# use NC counties polygons for an example
polys <- read_sf(system.file("shape/nc.shp", package = "sf"))

# pick two points, goal is to find the minimum # of counties to get between them
# for these two points, the shortest path between the points (straight line) actually crosses 5 counties, but there are paths of length 4 counties
points <- data.frame(name=c('start','end'),LAT=c(34.623858,35.257017),LONG=c(-78.851954,-78.237195))
points <- st_as_sf(points, coords = c("LONG", "LAT"), crs = 4267)

# show the problem:
plot(polys$geometry)
plot(points$geometry, add=TRUE,col='blue')
plot(st_cast(summarize(points,name=first(name)),'LINESTRING'),add=TRUE,col='blue')

# find the indices of the polygons of the two points
indices <- as.integer(st_intersects(points,polys))
indices # points are in county index 94 for start, and 62 for end

# convert to graph of polygons (polygons = vertices; edges = adjacency)
nc.nb <- spdep::poly2nb(polys)
g <- igraph::graph_from_adj_list(nc.nb)

plot(polys$geometry)
plot(polys[94,]$geometry,add=TRUE,col='red') # starting county
plot(polys[62,]$geometry,add=TRUE,col='red') # ending county

path <- igraph::shortest_paths(g, from = 94, to = 62) # find shortest path
pathint <- as.integer(path$vpath[[1]])
plot(polys[pathint,]$geometry,add=TRUE,col='red') # here's a shortest path plotted in red

length(pathint) # gives 4 -- the right answer