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?
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).
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
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 if geom.type() == QGis.Point: xy = geom.asPoint() x,y = xy, xy 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) with open(r'csvfile.csv','wb') as file: file.write('x,y,z,distance') file.write('\n') file.write(str(dist_min) + ',' + str(dist_min) + ',' + str(dist_min) + ',' + str(dist_min)) 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.
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
iv) Append everything (distance and point coordinates) to a list (
v) Get the minimum distance value together with the point information of the the
min(dist_list, key = lambda t:t)
[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)