Topography of a region - Python

Introduction

https://www.norgeskart.no/?backgroundLayer=toporaster&lon=822368.0472314445&lat=7767892.806770398&rotation=0&zoom=7&printTool=hiking

An example of Alta-joki Alta river.

https://caltopo.com/map.html#ll=68.8546,23.64807&z=8&b=mbt

Mostly needs an API.

• DEM Digital elevation model
• DSM Digital surface model
• DTM Digital terrain model

Python GIS

• folium
• geopy
• GeoPandas

Earth Model

DEM files (likes SRTM data) with

• Rasterio
• GDAL

Norwegian national DEM (Kartverket). The Norwegian Mapping Authority provides very high-resolution DEM data (down to 1 meter). https://hoydedata.no. The files in the .zip package

• .tif: pixel = elevation value
• .tfw (world file). Georeferencing info (pixel size, rotation, coordinates)
• .tif.aux.xml. May include Statistics (min/max elevation), Pyramid layers (for faster display), Projection overrides
• .prj Defines the coordinate reference system (CRS)
• .cpg. character encoding
• .shp. Contains shapes (points, lines, or polygons)
• .shx. Shape index file
• .dbf. Attribute table (like a spreadsheet)

Global DEMs (digital elevation model). No need super high resolution. Common datasets:

• SRTM (~30 m resolution)
• ASTER GDEM (~30 m)

Open portals. You can use

• OpenTopography
• OpenDEM
• NASA EarthData

which gives eg GeoTIFF.

Open Topo Data

https://www.opentopodata.org/

Python and GeoTIFF

import rasterio

with rasterio.open("dem.tif") as dem:
    for lon, lat in track:
        row, col = dem.index(lon, lat)
        elevation = dem.read(1)[row, col]

Downsampling the GeoTIFF

from rasterio.enums import Resampling

data = src.read(
    out_shape=(
        src.count,
        src.height // 10,
        src.width // 10
    ),
    resampling=Resampling.average
)

Path-finding

A depression filling algorithm, such as:

• Planchon–Darboux algorithm
• Wang & Liu algorithm
• Priority-flood (modern, very efficient). Overflow-hydro or RichDEM, WhoteboxTools.

Flow direction using:

• D8 flow direction (simplest, most common)
• D∞ (more accurate, distributes flow)
• MFD (multiple flow direction)

Treat the DEM as a graph. Nodes are raster cells and edges neighbor connections. The Cost function can be

1. Elevation drop (prefer downhill)
2. Slope
3. Accumulated flow (prefer existing rivers)

And use Dijkstra’s algorithm or A* (if you know destination).

Tools

GDAL (basic raster ops) TauDEM (excellent for hydrology) WhiteboxTools (very strong, modern) GRASS GIS (r.watershed, r.fill.dir)

Priority-flood Python

1. Start from the edges (they always drain outward)
2. Use a priority queue (min-heap) ordered by elevation
3. Expand inward
4. If a neighbor is lower → raise it (fill the sink)
5. Otherwise → keep it as-is

import heapq
import numpy as np

def priority_flood(dem):
    """
    Fill depressions in a DEM using Priority-Flood.
    
    Parameters:
        dem (2D numpy array): elevation grid
        
    Returns:
        filled_dem (2D numpy array)
    """
    rows, cols = dem.shape
    filled = dem.copy()
    visited = np.zeros((rows, cols), dtype=bool)
    
    heap = []

    # Step 1: push all edge cells into heap
    for r in range(rows):
        for c in [0, cols - 1]:
            heapq.heappush(heap, (filled[r, c], r, c))
            visited[r, c] = True

    for c in range(cols):
        for r in [0, rows - 1]:
            if not visited[r, c]:
                heapq.heappush(heap, (filled[r, c], r, c))
                visited[r, c] = True

    # 4-neighborhood (can switch to 8 if needed)
    neighbors = [(-1, 0), (1, 0), (0, -1), (0, 1)]

    # Step 2: process cells
    while heap:
        elev, r, c = heapq.heappop(heap)

        for dr, dc in neighbors:
            nr, nc = r + dr, c + dc

            if 0 <= nr < rows and 0 <= nc < cols and not visited[nr, nc]:
                visited[nr, nc] = True

                # If neighbor is lower → fill it
                if filled[nr, nc] < elev:
                    filled[nr, nc] = elev

                heapq.heappush(heap, (filled[nr, nc], nr, nc))

    return filled