Rasterio – How to Get Elevation from TIFF Using Python

elevationpythonrasterrasterio

I am trying to get the elevation from a TIFF file for each cell of the raster.

Therefore I thought about looping through every cell of the raster and fetching the elevation value of the current cell. My idea was to iterate throught the raster per row an column like a grid:
grid image

I tried this by counting the rows (dem.height) and columns (dem.width) to get the dimension of the raster. To iterate throught every row as well as every column inside a single row I used the while loops. Inside the nested column loop I wanted to fetch the elevation of the current cell via dem_data[row, col]. To have some kind of structure inside the fechted data, I wanted to use an array cell which contains the id of the cell (count), the coordinates (dem.xy(row, col)) and the elevation (dem_data[row, col]). The cell array is appended to the overall array arr to have this kind of structure:

[[id1, (lng1, lat1), elevation1], [id2, (lng2, lat2), elevation2], … ]

This is my complete code to achive this:

import rasterio

# Reading the DEM using rasterio
dem = rasterio.open(r'C:\...\dem.tif')
rows = dem.height
cols = dem.width
dem_data = dem.read(1).astype("float64")
print(rows, cols)
print(dem.crs)
print(dem.bounds)

# Creating empty list to store arrays of values
arr = list()

# Loop through every cell of the raster
count = 0
row = 0
while row < rows:
    col = 0
    while col < cols:
        cell = list()
        cell.append(count)
        # Get coordinates of current cell 
        cell.append(dem.xy(row, col))
        # Get elevation of current cell
        cell.append(dem_data[row, col])
        # Append helper list to overall list
        arr.append(cell)
        col += 1
        count += 1
    row += 1

print(arr)

But strangely I do not receive elevation values, like to ones which are displayed in my GIS. My GIS displays elevation values between 13.97 – 46.52. The array arrdisplays elevation values mostly 0.0 up to around 4

I tried a different approach by using specific points from a SHP-file with this code from this tutorial from YouTube. And using this code I received the correct elevation for these specific points. But I would like to have the elevation not only for specific, manual points but for the whole raster per cell.

Best Answer

You can read the contents of a raster (excluding the nodata values) using GDAL and numpy. See this example which works when pointing to an ArcGIS Pro Python installation.

import numpy as np
from osgeo import gdal
#The input raster
inras = r"Z:\GISpublic\pit.tif"
#Create the gdal object...
ds = gdal.Open(inras)

def get_value(ds):
    '''returns the nodata value, minimum, and maximum values of a single band raster'''
    band = ds.GetRasterBand(1)
    stat = band.GetStatistics(True, True)
    return (band.GetNoDataValue(), stat[0], stat[1])

raster_values = get_value(ds)
no_data = raster_values[0]

#Send the gdal object to a numpy array....
myarray = np.array(ds.GetRasterBand(1).ReadAsArray())
#Iterate the numpy array
for row in myarray: 
    for item in row:
        print(item)
Related Question