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 np.linalg.lstsq.
  • 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, [[1],[1]])
    ## 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[0] - 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
  • 1
    you are essentially solving the same problem as a GPS receiver. You know the distances to multiple points, and you know the exact locations of those points. I suggest you read up on how a GPS receiver derives its location
    – tomfumb
    Jan 3, 2017 at 17:42
  • Surveying Calculation plug-in can solve arc section (using hundred year old formulas) in 2D. For a 3D solution zenith angles are necessary.
    – Zoltan
    Jan 3, 2017 at 17:51
  • I am reading about it, and writing code that I can test, and understanding details of the problem. As of now, I think that I will look for ways to work in 2D. The question would be if a "Digital Handheld Laser Distance Meter" will give me decent horizontal estimations given a measurement taken at an angle.
    – mariotomo
    Jan 3, 2017 at 21:59
  • in my case, I needed transf = QgsCoordinateTransform(QgsCoordinateReferenceSystem(3117),QgsCoordinateReferenceSystem(4326)) and then transf.transform(QgsPoint(point[0], point[1]))
    – mariotomo
    Jan 4, 2017 at 20:16
  • @Zoltan, I have installed the surveying calculation plugin, it looks a lot more powerful than what I need, I will reserve some time to study it. thank you for the hint.
    – mariotomo
    Jan 5, 2017 at 0:33

1 Answer 1


I'm just sketching an answer, I will add to it bit by bit, as I work at DistanceMatrixToCoords, a small project on github, it defines a QGIS plugin that attempts to answer this question.

@tomfumb is right, the problem I want to solve is very similar to what a GPS receiver does, with the difference that my reference points, being distributed on a small portion of earth surface, are almost exactly coplanar, which makes altitude determination close to impossible.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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