1

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.

5
  • isn't what you want just dem_data?
    – Ian Turton
    Jul 13 at 13:19
  • Maybe zero is the no-data value for your raster?
    – GBG
    Jul 13 at 14:42
  • @IanTurton I thought, that dem_data is the data of the whole raster and I would like to get the elevation per single cell for every cell of the raster. That´s why I tried to access the cells via dem_data[row, col]
    – jonsken
    Jul 13 at 14:53
  • 2
    but you then loop through the entire raster querying each pixel to store it in a new array - This seems like an XY problem - could you describe what you are trying to do as well as what you have tried so far.
    – Ian Turton
    Jul 13 at 14:55
  • @GBG could be... But when I display the raster in my GIS there seems to be elevation data everywhere inside the raster
    – jonsken
    Jul 13 at 14:57

2 Answers 2

2

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)
2
  • works as described. I am just looking for a way to save the values (item) to an array to cluster everything. But when I define an empty array arr outside the for-loops and try to append the values with arr.append(item) everythings get somehow corrupted. Am I doing something wrong with .append()?
    – jonsken
    Jul 13 at 19:38
  • 1
    myarray is already a numpy array, which is different from a list but probably what you want. You should generally try to avoid using Python loops over numpy arrays, stick to numpy methods if you need to manipulate it further
    – mikewatt
    Jul 13 at 21:55
1

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)]
3
  • I am pretty new to python, so in general: Is numpy a library which I should try to get my hands on and use it in an essential way? Or is it more of a "it makes your life much more easier, so try to use it as much as possible" solution?
    – jonsken
    Jul 14 at 8:21
  • Don't ignore the pure-python basics, but if you plan on working with raster or pointcloud data a lot then it becomes pretty necessary, both from a performance standpoint and because raster and lidar libraries lean on numpy pretty heavily
    – mikewatt
    Jul 14 at 16:28
  • Also note that the format we achieve above format is much less compact than the original raster data, so make sure that's really what you need for whatever processing you're doing. It's essentially treating the raster like point data, which is a bit of a red flag to me
    – mikewatt
    Jul 14 at 16:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.