How to get the ground level altitude from DSM or Point Cloud using Python [closed]

Question

I am trying to get the ground level altitude from either DSM or PC. I am new to GIS and Python and still learning, any pointer or source?

I have a latitude and longitude of the point, for the same I want altitude from DSM or PC.

Answer 1

Usually DSM or DEM rasters are georeferenced Tiff files. You can use GDAL or rasterio library to load them into Numpy Array.

    import numpy as np
    from osgeo import gdal

    raster = gdal.Open("dem.tif")
    elevation_array = np.array(raster.GetRasterBand(1).ReadAsArray())

The georeferenced Tiff files will have information about corner coordinates of the 2D elevation raster, width and height per pixel, row rotation and column rotation as part of its GeoTransform. You can obtain this GeoTransfrom using following method.

transform = raster.GetGeoTransform()
print(transform)

It should print a tuple containing 6 values like

(x-coordinate of the upper-left corner of the upper-left pixel, west-east pixel resolution / pixel width, row rotation (typically zero), y-coordinate of the upper-left corner of the upper-left pixel, column rotation (typically zero),north-south pixel resolution / pixel height (negative value for a north-up image)

Refer documentation of GetGeoTransfrom

Once you have this transform values, you can use it to convert between row,col of the numpy array of the raster to lat,lon and reverse.

So given lon, lat (or x,y)

# x,y to col,row
col = int((x - transform[0]) / transform[1]) 
row = int((y - transform[3]) / transform[5])

now look up the elevation_array with col and row by

elevation = elevation_array[row][col]

this gives you elevation at that location.