I have a point P1 (see sketch) which is located underneath the surface (DTM). My aim is 1) to find the shortest 3D distance between P1 and the surface (d1 in sketch) and 2) the surface location (P2 in sketch) where the shortest 3D distance leads to. How can I accomplish this using QGIS or Grass or a similiar spatial open source GIS tool?

enter image description here


I have solved this task by using the QGIS Python console (QGIS 2.18.10) and GDAL (2.2.0).

I have a shapefile where P1 is stored called point_shape with an attribute altitude (decimal type).

  1. I've used a DEM (1m) with a size of 1x1km (about 1.000.000 Pixel) and converted it to an CSV-file with the OSGeo4W Shell coming along with the Windows installation:

    gdal2xyz dem.vrt dem.csv
  2. Open your Python console in QGIS and paste following code:

    import numpy as np
    point_array = np.genfromtxt (r'dem.csv', delimiter=" ")
    for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
        if lyr.name() == "point_shape":
            point_shape = lyr
    iter = point_shape.getFeatures()
    for feature in iter:
        geom = feature.geometry()
        attrs = feature.attributes()
        z = attrs[1]
        if geom.type() == QGis.Point:
            xy = geom.asPoint()
            x,y = xy[0], xy[1]
    dist_list = []
    for point in point_array:
        point_dist_list = []
        p1 = np.array(point)
        p2 = np.array([x, y, z])
        dist = np.linalg.norm(p1-p2)
    dist_min = min(dist_list, key = lambda t:t[0])
    with open(r'csvfile.csv','wb') as file:
        file.write(str(dist_min[1]) + ',' + str(dist_min[2]) + ',' + str(dist_min[3]) + ',' + str(dist_min[0]))
    for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
        if lyr.name() == "csvfile":
            csvfile = lyr   

After the first time you have run the code you can add the csv-file as a delimited text layer. After a second time this layer will be refreshed each time you use the code.

In detail:

i) Create a numpy array with the CSV-file

ii) Get the xy-coordinates and the altitude value from the shapefile layer

iii) For every element in the point_array calculate the euclidean distance to the given point p2 with np.linalg.norm(p1-p2)

iv) Append everything (distance and point coordinates) to a list (dist_list)

v) Get the minimum distance value together with the point information of the the dist_list with min(dist_list, key = lambda t:t[0])

Result: [16.202394884691543, 32620649.0, 5710659.0, 257.94]

vi) Write the result to a csv-file with headers x,y,z,distance

Reading the data, calculating the distances and sorting takes about 15 seconds (Windows 7 64 Bit, i7-5500U CPU @ 2,40GHz, 8 GB RAM)

enter image description here

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.