Calculating shortest distance between points within polygon using QGIS?

Question

I performed a transect in which I sampled several animals. The location of these animals are represented by a point shapefile. I am interested in calculating a pairwise matrix of distance (in meters) between all sampled animals. These distances will be used subsequently in spatial autocorrelation analyses to assess animal's dispersal capability (using genetic data).

However, calculating the euclidean distance is not adequate for my case because there is a river at north of these locations and it is considered to be a barrier to movement for the sampled animals. Thus, I rendered a shapefile polygon containing all these points and I intend to calculate the pairwise distance of points within the borders of this polygon. I uploaded a figure to help you understand exactly what I want. This figure contains 4 sampled locations (blue dots) and a yellow polygon from which I want to calculate the minimum distances. The red line represents the euclidean distance between the pair of points displayed. However, this distance to is not realistic because the river is in-between. The green line within the polygon is roughly what I want to estimate (basically a minimum path distance).

Any idea of how to do this with shapefiles or is it only possible to do with rasters?

There is other topic "Length of line on specific point (from point layer)" with a similar problem but that does not do exactly what I want. I have many more locations (besides the four shown) and calculate distances along the border of the polygon would not be realist.

Answer 1

I have no idea on how to do it with shapefile with existing tools, but an algorithm that iteratively adds vertices of the polygon to the line that joins the two points (until the new line is fully inside the polygon) should work. This would be the opposite of Douglas-Peucker.

With existing tools, I suggest that you convert your polygon to raster and assign a value equal to your pixel size for all the pixels inside the polygon (no data or VERY large value outside of the polygon). Then you can compute the minimum cost path using either GRASS or SAGA tools in the processing toolbox (see Cost Distance Analyst In QGIS? ). This will yield an exact but unprecise value, and you can increase the precision by reducing the cell size (which will also increase the processing time).

Answer 2

I have a REALLY rough solution to this question by doing the following steps:

1. run "Random points in polygons" processing algorithm on your with a small min distance between created points (the more points the better the final path, but more points means slower)

2. run "Join attributes by nearest" on your randomly generated points joining the auto-generated id field ("rand_point_id") to the x nearest neighbors (the more neighbors the better the final path, but more neighbors means slower)

3. run the following script in the console:

import networkx as nx

random_points = QgsProject.instance().mapLayersByName("Random points in polygons")[0]
points = QgsProject.instance().mapLayersByName("points")[0]
neighbored_points = QgsProject.instance().mapLayersByName("Joined layer")[0]
poly = QgsProject.instance().mapLayersByName("poly")[0]

f_start_point, f_end_point = list(points.getFeatures())

f_poly = next(poly.getFeatures(), None)
geom_poly = f_poly.geometry()
poly_eng = QgsGeometry.createGeometryEngine(geom_poly.constGet())
poly_eng.prepareGeometry()
all_neighbored_points = list(neighbored_points.getFeatures())

list_edge = []
print(len(all_neighbored_points))
for point in all_neighbored_points:
    pt1 = next(random_points.getFeatures("rand_point_id = {}".format(point["rand_point_id"])), None).geometry().asPoint()
    pt2 = next(random_points.getFeatures("rand_point_id = {}".format(point["rand_point_id_2"])), None).geometry().asPoint()
    line_geom = QgsGeometry.fromPolylineXY([pt1,pt2])
    if not poly_eng.contains(line_geom.constGet()):
        continue
    list_edge.append(((pt1.x(), pt1.y()), (pt2.x(), pt2.y()), {"weight": line_geom.length()}))

qsi = QgsSpatialIndex(random_points.getFeatures())

g = nx.Graph()
g.add_edges_from(list_edge)

p1 = f_start_point.geometry().asGeometryCollection()[0].asPoint()
n1_f = random_points.getFeature(qsi.nearestNeighbor(p1, 1)[0])
n1 = n1_f.geometry().asGeometryCollection()[0].asPoint()
n1 = (n1.x(), n1.y())

p2 = f_end_point.geometry().asGeometryCollection()[0].asPoint()
n2_f = random_points.getFeature(qsi.nearestNeighbor(p2, 1)[0])
n2 = n2_f.geometry().asGeometryCollection()[0].asPoint()
n2 = (n2.x(), n2.y())

print(n1)
print(n2)

list_nodes = list(nx.bidirectional_dijkstra(g, n1, n2, weight='weight'))[-1]
geom = QgsGeometry.fromPolylineXY([p1] + [QgsPointXY(*node) for node in list_nodes] + [p2])
print(geom.asWkt())

This can be really slow if you generate a lot of points because the number of edges in the final graph is based on the first two steps: number_of_edges = nb_random_points * nb_neighbors

example of result: using 12500 edges in the graph (500 random points * 25 neighbors)