Get elevation from TIFF using rasterio (Python)

Question

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:

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.

Answer 1

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)

Answer 2

Here's how you could achieve your desired format without Python loops:

import numpy as np
import rasterio

path = 'example.tif'

with rasterio.open(path) as f:
    arr = f.read(1)
    mask = (arr != f.nodata)
    elev = arr[mask]
    col, row = np.where(mask)
    x, y = f.xy(col, row)
    uid = np.arange(f.height * f.width).reshape((f.height, f.width))[mask]

result = np.rec.fromarrays([uid, x, y, elev],
                           names=['id', 'x', 'y', 'elev'])

print(result.dtype)
print(result[:10])

1. Read all data into a 2d array
2. Create a binary mask indicating where data is valid
3. Mask off 2d array down to 1d array of valid values
4. Grab column and row indices of valid values
5. Convert column and row indices to map coordinates
6. Create a sequential ID for each pixel, mask it down to valid values
7. Merge all of that into a single array

The result is a numpy recarray that looks like:

(numpy.record, [('id', '<i8'), ('x', '<f8'), ('y', '<f8'), ('elev', 'u1')])

[(0, -8047.8756175, 10524.2996175, 40)
 (1, -8046.1748525, 10524.2996175, 40)
 (2, -8044.4740875, 10524.2996175, 44)
 (3, -8042.7733225, 10524.2996175, 43)
 (4, -8041.0725575, 10524.2996175, 35)
 (5, -8039.3717925, 10524.2996175, 35)
 (6, -8037.6710275, 10524.2996175, 36)
 (7, -8035.9702625, 10524.2996175, 37)
 (8, -8034.2694975, 10524.2996175, 38)
 (9, -8032.5687325, 10524.2996175, 41)]