Topography of a region - Python

Introduction

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

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

<syntaxhighlight lang="python"> 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

<syntaxhighlight>