The task I'm giving myself is about enabling people with little equipment to produce a map of a forest, or a garden with an important presence of trees.
I participate to openstreetmap and I plan to (enable people to) contribute tree locations for botanic gardens in areas where Bing aerial coverage is not available, or where it is useless because tree canopy is too thick to allow for identification of individual trees.
My hope is that there is a QGIS plug-in (and if it does not exist, to develop one myself) that accepts an amount of reference points distributed around the area of interest, a matrix of distances among points (referenced or not), and populates a vector(points) layer with the calculated 3D coordinates.
I have been looking at trilateration, and I found this 22 years old article.
I suppose the math still holds.
I do not want to work with angles, only distances.
not yet an answer, but I have been working at it, and I have a couple more details about the issue.
- I need my input coordinates in metres (or feet, or millimetres), not in degrees, because I need an isometric coordinate system.
- The math in that article still holds, and numpy helps a lot, specifically, the article suggests solving the overdefined problem by least squares, and numpy contains the function
- I gave up hope to work in 3D and I am going to assume earth is flat, but really flat like the Netherlands, or at least that the digital device I will acquire to measure distances is able to split distances in their horizontal and vertical components.
- I have a working prototype, but I still need to discover how to create a point feature using the coordinates in the coordinate system of the layer.
the code for the coordinates calculation is not at all long, while the difficult part was about choosing which point I can reference given the points I have already calculated:
def find_point_coordinates(points, distances, point_id): import numpy as np from numpy.core.umath_tests import matrix_multiply connected_to = [id for id in sorted(distances[point_id]) if points[id].get('coordinates')] connected_matrix = np.array([points[id]['coordinates'] for id in connected_to]) A = connected_matrix[1:,] - connected_matrix[0,] A = 1.0 * A # make sure we work with floating point values ## squared distances vector, beacon_i to first beacon for which we have distances D_i1_2 = matrix_multiply(A*A, [,]) ## distances of targeted point from used reference points dfb_sel = np.array([distances[point_id][ref_id] for ref_id in connected_to]) r2 = dfb_sel * dfb_sel rhs = ((r2 - r2[1:]).reshape(D_i1_2.shape) + D_i1_2) / 2.0 r1, r2, r3, r4 = np.linalg.lstsq(A, rhs.reshape(rhs.shape[:1])) return connected_matrix[0,] + r1