Building the shortest path between two points inside polygon using QGIS

Question

With QGIS I need to build the shortest route between two points inside a polygon. I have a vector polygon layer and I need to build a path inside it from the starting point to the ending point. How can I do this?

I found that the vector layer needs to be rasterized, but I don’t understand what to do next. Existing tools are based on vector line layers. I studied so much information and nothing helps. I would write a Python module or use ready modules.

Answer 1

There are probably better approaches available, however I am suggesting the following intuitive workaround for Vector Data.

Let's assume there is a point layer (two features), and a polygon layer (one feature), see the image below.

Step 1. Apply the "Polygon Divider" plugin to break the original polygon feature into pieces of equal areas. Alternatively, one can use the "Create grid" geoalgorithm and overlap it with the initial polygon feature

Step 2. Select all objects from the previous output, without overlap with original points. Afterward, get "Centroids" for the selected pieces only

Step 3. "Extract vertices" of the original polygon

Step 4. Merge original points with outcomes of Steps 2 and 3 with the "Merge vector layers" tool

On this step it is vital to mention, that the next approach should be chosen wisely, depending on precision/tolerance/complexity etc.

---

Approach #1

Step 5. Triangulate points with the "Voronyi/Delaunay triangulation"

Step 6. "Extract by location" triangles, that are within the original polygon

Step 7. Convert these triangles to polylines, with "Polygons to lines"

Step 8. "Explode lines" from the previous step, and also "Delete duplicate geometries"

Step 9. Apply the "Shortest path (point to point)" and get the output like this:

Simplify and/or smooth if needed.

---

Approach #2

Step 5. Use the "Geometry by expression" together with this expression:

collect_geometries(
    array_filter(
        array_foreach(
            overlay_nearest(
                layer:='points_to_test', -- change accordingly
                expression:=@geometry,
                limit:=10 -- adjust properly
                ),
            make_line($geometry, @element)
            ),
        within(@element, aggregate(
                                layer:='polygon', -- change accordingly
                                aggregate:='collect',
                                expression:=$geometry
                                )
            )
        )
    )

Step 6. "Multipart to singleparts", and if required "Delete duplicate geometries"

Step7. Apply the "Shortest path (point to point)" and get the output like this:

Simplify and/or smooth if needed.

---

Approach #3

Step 5. Create all possible line connections between merged points like in this tread Creating all possible line segments between all points using QGIS. And after extract only hub lines that are within the original polygon. Do not forget about the spatial index. And also be careful with polygon edges, do not skip them. Apply the "Multipart to singleparts" if needed

Step 6. Apply the "Shortest path (point to point)" and get the output like this:

---

Once can also compare the results in the end. With some pros&cons one will find a suitable approach.

---

Notes:

• for the Step 1 the area value has to be chosen wisely
• this workflow was only performed on a single polygon and two points, not multiple inputs of each layer
• this workflow can be converted into a QGIS-model or pure PyQGIS code for automized routine work
• mind some issues, that could happen. This one evoked after diving the polygon
• In the first approach, triangulation result may not always correspond to the initial polygon
• In the second approach, some polygon sides can be skipped'

---

References:

• Splitting all lines in layer at vertices using QGIS

Answer 2

Basic solution: shortest path/least cost

A very easy way is a raster least cost analysis. The only cost is distance in this case what makes it very simple:

1. Rasterize the polygon. Use a resolution small enogh that the smallest "bottleneck" of the polygon is still at least 3 to 5 times larger than the pixel size. For A fixed value to burn use 1.

2. Install Least Cost Path plugin and run it (from Menu Processing > Toolbox > Cost Distance analysis) with the raster created before and the start- and end-point.

Variant: If the path should not run along the polygon's boundary, before step 1 apply a (small) negative buffer to the polygon.

Black shape: rasterized shape of the initial polygon; red line: shortest path:

Variant: prefer path in the middle of the polygon

A bit more sophisticated, you can create a path that prefers a path in the middle of the polygon - e.g. to model a ship on a lake/river that always wants to stay in the middle where the water is deeper, not near the shore.

1. Create a polygon grid with the same resolution as the raster. Select by location the cells that are within the polygon, inverse selection, delete selected to keep only grid cells completely inside the polygon.

2. Use Field calculator to create a new attribute cost that measures the inverse distance from each grid cell to the shore with this expression:

    1/distance(
        centroid ($geometry),
        closest_point (
            overlay_within(
                'polygon',
                boundary($geometry)
            )[0],
            centroid ($geometry)
        )
    )
   

3. Rasterize (as in step 1 above), but this time use Field to use for a burn-in value and select the attribute cost.

4. Run Least Cost Path plugin as in step 2 above.

Yellow line shows the result of the refined approach; in the background the raster with values from white (high costs) to blue (low costs), so the algorithm tries to create the shortest distance path, but also tries to avoid the white cell and prefers the darker blue cells; see for comparision the red line created above:

The line simplified and smoothed:

"Bottleneck" of the polygon (SSE of the red point in the picture above): the path prefers the darker pixels; pixel size is chosen in a way that even the bottleneck is covered with at least 2 to 3 neighboring pixels:

Answer 3

I don't have a correct answer, but the idea is like this

1. Create a visibility graph
2. Explore all the paths or use A* algo to find the shortest path.
3. Create a line from that in QGIS.

Although the current solution doesn't give correct answer you can get the idea.

import networkx as nx
from qgis.core import *
import processing
from shapely.geometry import Polygon, Point, LineString
from shapely.ops import unary_union
def visibility_graph_2(polygon, start, end):
    # Create graph
    G = nx.Graph()
# Add start and end points
G.add_node((start.x, start.y))
G.add_node((end.x, end.y))

# Add polygon vertices to graph
for vertex in polygon.exterior.coords:
    G.add_node(vertex)

# Connect visible nodes
for node1 in G.nodes:
    for node2 in G.nodes:
        if node1 != node2 and LineString([node1, node2]).intersects(polygon) and not LineString([node1, node2]).crosses(polygon):
            G.add_edge(node1, node2, weight=Point(node1).distance(Point(node2)))

return G


def visibility_graph(polygon, start, end, num_steiner_points=10):
    # Create graph
    G = nx.Graph()
# Convert start and end points to tuples
start = (start.x, start.y)
end = (end.x, end.y)

# Add start and end points
G.add_node(start)
G.add_node(end)

# Add polygon vertices to graph
for vertex in polygon.exterior.coords:
    G.add_node(vertex)

# Add Steiner points to graph
coords = list(polygon.exterior.coords[:-1])
for i in range(len(coords) - 1):
    for j in range(num_steiner_points):
        point = ((1 - j/num_steiner_points) * coords[i][0] + (j/num_steiner_points) * coords[i+1][0],
                 (1 - j/num_steiner_points) * coords[i][1] + (j/num_steiner_points) * coords[i+1][1])
        G.add_node(point)

# Connect visible nodes
for node1 in G.nodes:
    for node2 in G.nodes:
        if node1 != node2:
            line = LineString([node1, node2])
            # Only add the edge if the line is entirely inside the polygon
            if polygon.contains(line):
                G.add_edge(node1, node2, weight=Point(node1).distance(Point(node2)))

return G


def shortest_path(polygon, start, end, path_shortest_path):
    # Create visibility graph
    print('creating  visibility graph')
    G = visibility_graph(polygon, start, end)
# Find shortest path
print('Find shortest path')
path = nx.dijkstra_path(G, source=(start.x, start.y), target=(end.x, end.y), weight='weight')

# Create a new memory layer to store the shortest path
print('Create a new memory layer to store the shortest path')
vl = QgsVectorLayer("LineString", "shortest_path", "memory")
pr = vl.dataProvider()

# Create a new feature
seg = QgsFeature()

# Add the shortest path to the feature as a polyline
seg.setGeometry(QgsGeometry.fromPolyline([QgsPoint(pt[0], pt[1]) for pt in path]))

# Add the feature to the layer
pr.addFeatures([seg])

# Update extent of the layer
vl.updateExtents()

# Write layer to a new shapefile
QgsVectorFileWriter.writeAsVectorFormat(vl, path_shortest_path, "utf-8", vl.crs(), "ESRI Shapefile")

return path


path_polygon = '/Users/manish.sahu/Downloads/polygon.shp'
path_point = '/Users/manish.sahu/Downloads/point.shp'
path_shortest_path = '/Users/manish.sahu/Downloads/path_shortest.shp'
polygon_layer = QgsVectorLayer(path_polygon, 'polygon', 'ogr')
points_layer = QgsVectorLayer(path_point, 'point', 'ogr')
if not layer.isValid():
    print("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(layer)
Assuming 'polygon_layer' and 'points_layer' are your loaded layers
###################################
Get the first (and only) feature from the polygon layer
polygon_feature = polygon_layer.getFeature(0)
Convert the feature to a Shapely Polygon
polygon = Polygon(polygon_feature.geometry().asMultiPolygon()[0])
polygons = [p for p in polygon_feature.geometry().asMultiPolygon()]
Filter out polygons with less than 3 points
selected_polygon = Polygon(polygons[0][0])
Get the features from the points layer
points_features = points_layer.getFeatures()
Convert the features to Shapely Points
points = [Point(feature.geometry().asPoint()) for feature in points_features]
The start and end points are the first and second points
start, end = points[0], points[1]
Calculate the shortest path and write it to a shapefile
path = shortest_path(selected_polygon, start, end, path_shortest_path)
###################################
Create a new layer from the shapefile
output_layer = QgsVectorLayer(path_shortest_path, "Shortest Path", "ogr")
Check if the layer is valid
if not output_layer.isValid():
    print("Layer failed to load!")
else:
    # Add the layer to the project
    QgsProject.instance().addMapLayer(output_layer)