Selecting coordinates of points below DEM (DTM)

Question

I have a DEM file and a big list of X,Y,Z point (about 1M). Each of them can be below or above the DEM file.

Is there a quick way to select the x,y,z points below the DEM file using Python?

I did not find yet a solution on Google or StackOverflow.

Answer 1

A solution is to get the DEM elevation values of the same points xy coordinates (red lines projections in the figure below) using osgeo.GDAL or rasterio to read the DEM, GeoPandas for the points shapefile and affine to get the elevation value (see Python affine transforms).

With GDAL

from osgeo import gdal
from affine import Affine
import geopandas as gpd

Open the dem

layer = gdal.Open('dem.tif')
gt =layer.GetGeoTransform()
T0 = Affine.from_gdal(*gt)
band = layer.GetRasterBand(1)
elev= band.ReadAsArray()

The shapefile

df = gpd.read_file("points.shp")
df
   ID    z                     geometry
0  18  140  POINT (204439.9038755086 89773.5004016017)
1  37  120  POINT (204404.5155014408 89844.2771497373)
2  25  280  POINT (204333.7387533052 89808.8887756695)
3  43  240  POINT (204262.9620051695 89879.6655238051)
4  60  110  POINT (204510.6806236442 89915.0538978729)
5  56  206  POINT (204227.5736311017 89773.5004016017)

The affine projection to get the DEM elevation values -> col, row = ~T0 * (x, y)

df['x'] = df.geometry.apply(lambda p: p.x)
df['y'] = df.geometry.apply(lambda p: p.y)
xy2rc = lambda x, y: [int(i) for i in [x, y] * ~T0]
def zdem(x,y):
   cc = xy2rc(x,y)
   return elev[cc[1],cc[0]]
df['elev'] = df.apply(lambda p: zdem(p.x,p.y), axis=1)
df
   ID    z                         geometry                 x             y            elev
0  18  140  POINT (204439.9038755086 89773.5004016017)  204439.903876  89773.500402  210.160004
1  37  120  POINT (204404.5155014408 89844.2771497373)  204404.515501  89844.277150  209.169998
2  25  280  POINT (204333.7387533052 89808.8887756695)  204333.738753  89808.888776  208.399994
3  43  240  POINT (204262.9620051695 89879.6655238051)  204262.962005  89879.665524  204.380005
4  60  110  POINT (204510.6806236442 89915.0538978729)  204510.680624  89915.053898  209.410004
5  56  206  POINT (204227.5736311017 89773.5004016017)  204227.573631  89773.500402  205.830002

Points above the DEM

df[df['z'] > df['elev']]
   ID    z                         geometry                 x             y            elev
2  25  280  POINT (204333.7387533052 89808.8887756695)  204333.738753  89808.888776  208.399994
3  43  240  POINT (204262.9620051695 89879.6655238051)  204262.962005  89879.665524  204.380005
5  56  206  POINT (204227.5736311017 89773.5004016017)  204227.573631  89773.500402  205.830002

Points below the DEM

df[df['z'] < df['elev']]
   ID    z                         geometry                 x             y            elev
0  18  140  POINT (204439.9038755086 89773.5004016017)  204439.903876  89773.500402  210.160004
1  37  120  POINT (204404.5155014408 89844.2771497373)  204404.515501  89844.277150  209.169998
4  60  110  POINT (204510.6806236442 89915.0538978729)  204510.680624  89915.053898  209.410004

With rasterio

You don't need the affine library (integrated in rasterio)

import rasterio as rio
from rasterio import Affine
# open the dem
with rio.open('dem.tif') as dataset:
    elev = dataset.read(1)
T0 = dataset.transform
# the shapefile
...... (same)