6

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.

enter image description here

10
  • 2
    Please do not comment on your own post. If you take the Tour, you'll see that our purpose is to compile a list of questions with answers. Thanks and apologies detract from question quality, making it harder to understand what it is that you want. This particular question begs for a graphic to better explain what you want.
    – Vince
    Jun 1, 2017 at 10:29
  • 1
    Possible duplicate of Length of line on specific point (from point layer)
    – mgri
    Jun 1, 2017 at 11:03
  • If I understand your issue, it should be very similar to the question linked above.
    – mgri
    Jun 1, 2017 at 11:03
  • @mgri The lines shown in the figure are not shapefiles but rather a raw representation of the type of distance that I want to calculate. Looking at the question linked above, the poster was looking for a distance along a shapefile polyline to a destination point. I could do that but I have tens of points and I would have to draw several polylines Jun 1, 2017 at 11:25
  • 1
    @AndréLourenço I got it. Maybe, the linked answer doesn't work only when it is possible to draw the distance using a straight line (i.e. without going outside the polygon boundary). In the other cases (for example for the two points on the most left part of the attached image) the problem should be solved as in the linked answer. If you can confirm that, I can slightly edit that code and post it here, otherwise I suggest you adding more information on your original question for better explaining the issue.
    – mgri
    Jun 1, 2017 at 12:43

2 Answers 2

4

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).

0

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: enter image description here using 12500 edges in the graph (500 random points * 25 neighbors)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.