Shortest 3D distance (depth, z) between an underground point and the surface (dtm)

Question

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?

Answer 1

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)
       point_dist_list.append(dist)
       point_dist_list.extend(point)
       dist_list.append(point_dist_list)
   
   dist_min = min(dist_list, key = lambda t:t[0])
   
   with open(r'csvfile.csv','wb') as file:
       file.write('x,y,z,distance')
       file.write('\n')
       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   
           csvfile.triggerRepaint()
       else:
           pass
   

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)