Which Way Would Rain Flow from Your Roof?

Start from Bhaktapur, then try these.

• Pokhara, Jumla, Darchula. Each drop finds its own river: the Seti and the Gandaki, the Karnali, the Mahakali. The dashed line on the graph is where it crosses into India; Jumla’s and Darchula’s drops then run east across the plains in the Ghaghara and the Ganges.
• Namche. A drop at Namche Bazaar, at 3,150 m, runs down the Dudh Koshi into the Koshi. On the graph, nearly nine-tenths of its fall comes in the first quarter of the way, as for every town: steep mountains, nearly flat plains.
• Pits: As they are. Now the heights are used exactly as measured. Bhaktapur’s drop stops at about 140 m, on the Bagmati’s valley floor near the foot of the hills, in a square with no lower neighbor: a pit. Every town’s drop stops like this, in Nepal, where the measured heights rise a little before they fall again.
• Rain everywhere, with the pits filled, then as they are. With the pits filled, the river lines join into a few great trees, one for each of Nepal’s river systems. As they are, the same lines break off in the plains.
• Tibet, north of the mountains, with As they are. Nearly every drop there stops (462 of 500 random places when this page was made): the plateau is full of flat ground and of real basins with no way out. With Filled, some leave the map to the north, and others cross the Himalaya southwards in the gorges of rivers that rise in Tibet, such as the Arun.

A map made of numbers

The heights come from the Shuttle Radar Topography Mission (SRTM): for 11 days in February 2000, the Space Shuttle Endeavour carried a radar that measured the height of the ground under most of the world. The version used here, NASA’s own, has one height every 30 m or so (one every second of arc). That is far too many numbers for a page, so for the map above each square of about 1.5 km gets one number, and for your program each square of about 3 km.

Which number? Not the average height of the square: a river in a narrow gorge is lower than nearly everything around it, and an average hides it under the hills on both sides. On 3 km squares made with averages, the drop from Bhaktapur escapes the valley the wrong way, west towards the Trishuli. So each square takes the lowest height inside it, the bottom of whatever valley crosses it. With that, the gorge at Chobhar survives at 3 km squares.

The lowest height has one danger: one wrong number wins. Where the radar couldn’t see into deep Himalayan valleys, the copy used here filled a few gaps badly: a patch on Annapurna II reads 606 m where the mountain is about 6,000 m high. So steep-sided hollows over 300 m deep were first filled to where they would spill. Real ground has no such holes.

One rule

For each square, send the water to the lowest of its eight neighbors, if that is lower than the square itself. That is all. A drop is that rule repeated until it reaches the edge of the map or a square with no lower neighbor.

Repeat it for every square at once and you get a flow direction for each square, an arrow to one neighbor. Follow the arrows backwards from any square and you find all the land whose water passes through it, its drainage area or catchment. Where the drainage area is large, there is a river.

Pits, and how to fill them

A square with no lower neighbor is a pit. A few are real: a lake with no outlet, such as Paiku Co in Tibet, north of Langtang. Most are not. Before filling, the map above has 8,112 pits; 6,645 of them are below 300 m, in the plains, and most of the rest are on the Tibetan Plateau, above 3,000 m. Inside Nepal there are 691 pits among 61,010 squares, 465 of them below 300 m.

The fix is to do what water would do: fill each pit until it spills. Imagine the sea rising slowly from every edge of the map, always creeping in at the lowest point it can reach. When it reaches a hollow, the hollow fills to the water’s level, and the sea flows on. This is called Priority-Flood, and it raised 14,851 of the map’s squares, half of them by about a meter or less (about a hundred by more than 100 m), so that every drop finds its way to the edge. The pits you saw with As they are are gone with Filled, and the rivers join up.

How good is it?

To check, the map’s rivers were compared with MERIT Hydro, the world’s flow directions made in 2019 by Dai Yamazaki and his colleagues from the same radar heights, corrected, with rivers seen from space and from OpenStreetMap added. The site’s Rivers of the World is drawn from it.

• Of the squares where our rule finds a river draining 1,000 km² or more, 93.3% lie on one of MERIT Hydro’s rivers. Above 500 m, in the hills and mountains, all of them do; below 200 m, in the plains, 89.0%.
• Drops started from 2,000 random places in Nepal leave the country by the same river system as MERIT Hydro’s for 94.8% of them. Of those that go wrong, about two in five start below 300 m, on the plains, where only about one in five of all the places are.
• The land each system drains where it leaves Nepal comes out close:

The big rivers are right to a few percent, most likely because their deep valleys leave the water nowhere else to go. The Bagmati probably lost its 8% on the flat Terai, where it runs beside other small rivers and a few meters decide which one a square joins. The Mahakali’s shortfall was tested: much of its land is in India, west of the map’s edge at 80° E; with the map widened to 79° E it came out 4% over. (The Karnali is measured at Chisapani, above where it splits into two channels on the plains.)

Your program works on the Bagmati: from the Kathmandu Valley to the Indian border, 42 rows of 32 squares about 3 km across, each the lowest height in its square. The numbers are in the file bagmati.py, which the programs below import (open it to see them). show_grid is this page’s own way to draw a list of rows as a picture (on your own computer you would use a library such as Pillow).

Step 1: the map

First look at the numbers. Low ground is drawn dark and high ground light; the square root spreads out the low numbers, so the plains don’t all come out the same black.

from bagmati import heights, placesfrom pixels import show_grid ROWS, COLS = len(heights), len(heights[0])print(ROWS, 'rows of', COLS, 'squares')for name, (r, c) in places.items():    print(f'{name}: {heights[r][c]} m')everything = [h for row in heights for h in row]print('lowest', min(everything), 'm, highest', max(everything), 'm') def gray(h):                     # low ground dark, high ground light    return round(40 + 215 * (h / 2100) ** 0.5) show_grid([[gray(h) for h in row] for row in heights])

42 rows of 32 squaresBhaktapur: 1311 mKathmandu: 1276 mChobhar: 1262 mHetauda: 433 mBorder: 62 mlowest 26 m, highest 2012 m

Picture

The Kathmandu Valley is the smooth patch near the top, among patchy hills. Below it, the hills grow lower and darker towards the south, and the plains, all much the same dark gray, fill the bottom half.

Step 2: one drop

downhill looks at the eight neighbors of a square and gives back the lowest one, if it is lower; otherwise None. The 0 <= rr < ROWS test keeps it on the map. drop repeats it and keeps the way in a list.

from bagmati import heights, placesfrom pixels import show_grid ROWS, COLS = len(heights), len(heights[0])NEIGHBORS = [(-1, -1), (-1, 0), (-1, 1),             (0, -1),           (0, 1),             (1, -1),  (1, 0),  (1, 1)] def downhill(grid, r, c):    # The lowest of the eight squares around (r, c), if it is lower; else None.    best, lowest = None, grid[r][c]    for dr, dc in NEIGHBORS:        rr, cc = r + dr, c + dc        if 0 <= rr < ROWS and 0 <= cc < COLS and grid[rr][cc] < lowest:            best, lowest = (rr, cc), grid[rr][cc]    return best def drop(grid, r, c):    way = [(r, c)]    while downhill(grid, r, c) is not None:        r, c = downhill(grid, r, c)        way.append((r, c))    return way way = drop(heights, *places['Bhaktapur'])print(len(way), 'squares')for r, c in way[:6]:    print((r, c), heights[r][c], 'm')r, c = way[-1]print('stops at', (r, c), heights[r][c], 'm') def gray(h):    return round(40 + 215 * (h / 2100) ** 0.5) picture = [[gray(h) for h in row] for row in heights]for r, c in way:    picture[r][c] = 'royalblue'show_grid(picture)

23 squares(6, 14) 1311 m(6, 13) 1300 m(6, 12) 1293 m(6, 11) 1282 m(5, 10) 1276 m(6, 9) 1262 mstops at (23, 16) 113 m

Picture

The first six squares are Bhaktapur, west across the valley floor to Kathmandu at (5, 10) and the gorge at (6, 9): your program found the way out of the valley by itself. Then the drop runs south down the Bagmati, until it stops at row 23, at 113 m, just where the river comes out of the hills onto the plains. The border is still about 45 km to the south.

Step 3: the pits

A pit is a square with no lower neighbor, that isn’t on the edge of the map (water on the edge leaves the map).

from bagmati import heightsfrom pixels import show_grid ROWS, COLS = len(heights), len(heights[0])NEIGHBORS = [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)] def downhill(grid, r, c):    best, lowest = None, grid[r][c]    for dr, dc in NEIGHBORS:        rr, cc = r + dr, c + dc        if 0 <= rr < ROWS and 0 <= cc < COLS and grid[rr][cc] < lowest:            best, lowest = (rr, cc), grid[rr][cc]    return best def on_edge(r, c):    return r == 0 or c == 0 or r == ROWS - 1 or c == COLS - 1 pits = []for r in range(ROWS):    for c in range(COLS):        if not on_edge(r, c) and downhill(heights, r, c) is None:            pits.append((r, c))print(len(pits), 'pits, at heights:')print(sorted(heights[r][c] for r, c in pits)) def gray(h):    return round(40 + 215 * (h / 2100) ** 0.5) picture = [[gray(h) for h in row] for row in heights]for r, c in pits:    picture[r][c] = 'crimson'show_grid(picture)

24 pits, at heights:[43, 45, 49, 52, 56, 57, 58, 59, 59, 62, 63, 64, 64, 64, 64, 65, 65, 66, 67, 67, 69, 71, 74, 113]

Picture

Of the 24 pits, 23 are below 75 m, on the plains. The last, at 113 m, is the square where your drop stopped, where the hills end: its lowest neighbor is 116 m, just 3 m higher. Not one pit is up in the hills, where your drop ran without trouble.

Step 4: let the sea rise

Now fill the pits, as the sea rising from the edges would. The sea’s shore is kept in a heap (Python’s heapq, the priority queue of Grade 11): a list from which heappop always takes out the lowest square. The sea starts on every edge square; each time, it takes the lowest square on its shore and reaches its neighbors. A neighbor lower than the water is a hollow: it is raised to 1 cm (0.01 m) above the water, so its water still flows on, out the way the sea came in.

import heapqfrom bagmati import heights, placesfrom pixels import show_grid ROWS, COLS = len(heights), len(heights[0])NEIGHBORS = [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)] def downhill(grid, r, c):    best, lowest = None, grid[r][c]    for dr, dc in NEIGHBORS:        rr, cc = r + dr, c + dc        if 0 <= rr < ROWS and 0 <= cc < COLS and grid[rr][cc] < lowest:            best, lowest = (rr, cc), grid[rr][cc]    return best def drop(grid, r, c):    way = [(r, c)]    while downhill(grid, r, c) is not None:        r, c = downhill(grid, r, c)        way.append((r, c))    return way def on_edge(r, c):    return r == 0 or c == 0 or r == ROWS - 1 or c == COLS - 1 def fill(grid):    filled = [row[:] for row in grid]    shore = []                      # a heap of (height, r, c)    reached = set()    for r in range(ROWS):        for c in range(COLS):            if on_edge(r, c):                heapq.heappush(shore, (filled[r][c], r, c))                reached.add((r, c))    while shore:        level, r, c = heapq.heappop(shore)      # the lowest square the sea has reached        for dr, dc in NEIGHBORS:            rr, cc = r + dr, c + dc            if 0 <= rr < ROWS and 0 <= cc < COLS and (rr, cc) not in reached:                reached.add((rr, cc))                if filled[rr][cc] <= level:     # a hollow: raise it to just above the water                    filled[rr][cc] = level + 0.01                heapq.heappush(shore, (filled[rr][cc], rr, cc))    return filled filled = fill(heights)raised = [(r, c) for r in range(ROWS) for c in range(COLS) if filled[r][c] > heights[r][c]]print(len(raised), 'squares raised')way = drop(filled, *places['Bhaktapur'])r, c = way[-1]print(len(way), 'squares; ends at', (r, c), 'on the edge:', on_edge(r, c))print('passes the border:', places['Border'] in way) def gray(h):    return round(40 + 215 * (h / 2100) ** 0.5) picture = [[gray(h) for h in row] for row in heights]for r, c in raised:    picture[r][c] = 'lightskyblue'for r, c in way:    picture[r][c] = 'royalblue'show_grid(picture)

41 squares raised41 squares; ends at (41, 8) on the edge: Truepasses the border: True

Picture

Only 41 squares had to be raised, all on the plains or where the hills end, none by more than about 20 m. The drop now runs to the bottom edge, in India, crossing the border in the same square as MERIT Hydro’s Bagmati.

Why a heap? The sea must always grow from its lowest point, or it would fill a hollow to the height of some higher rim it reached first. With a plain list you would have to search the whole shore for its lowest square every time. A heap finds it at once, so the fill takes one pass over the map.

Step 5: how much land drains through each square

Now every square except those on the edge has a square downhill. Turn the arrows round: uphill[(r, c)] lists the squares whose water comes straight to (r, c). Then the land draining through a square is the square itself plus all the land draining through each of its uphill squares, a recursive definition, like a folder’s size in Grade 11. The dictionary land remembers each answer, so no square is worked out twice.

import heapqfrom bagmati import heights, placesfrom pixels import show_grid ROWS, COLS = len(heights), len(heights[0])NEIGHBORS = [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)]SQUARE_KM2 = 9.9                    # one square is about 3.0 km by 3.3 km def downhill(grid, r, c):    best, lowest = None, grid[r][c]    for dr, dc in NEIGHBORS:        rr, cc = r + dr, c + dc        if 0 <= rr < ROWS and 0 <= cc < COLS and grid[rr][cc] < lowest:            best, lowest = (rr, cc), grid[rr][cc]    return best def on_edge(r, c):    return r == 0 or c == 0 or r == ROWS - 1 or c == COLS - 1 def fill(grid):    filled = [row[:] for row in grid]    shore, reached = [], set()    for r in range(ROWS):        for c in range(COLS):            if on_edge(r, c):                heapq.heappush(shore, (filled[r][c], r, c))                reached.add((r, c))    while shore:        level, r, c = heapq.heappop(shore)        for dr, dc in NEIGHBORS:            rr, cc = r + dr, c + dc            if 0 <= rr < ROWS and 0 <= cc < COLS and (rr, cc) not in reached:                reached.add((rr, cc))                if filled[rr][cc] <= level:                    filled[rr][cc] = level + 0.01                heapq.heappush(shore, (filled[rr][cc], rr, cc))    return filled filled = fill(heights)uphill = {}for r in range(ROWS):    for c in range(COLS):        if not on_edge(r, c):            uphill.setdefault(downhill(filled, r, c), []).append((r, c)) land = {}def drains(r, c):    # how many squares' water passes through (r, c), its own included    if (r, c) not in land:        land[(r, c)] = 1 + sum(drains(rr, cc) for rr, cc in uphill.get((r, c), []))    return land[(r, c)] for name in ['Chobhar', 'Border']:    n = drains(*places[name])    print(f'{name}: {n} squares, about {n * SQUARE_KM2:,.0f} km²') def gray(h):    return round(40 + 215 * (h / 2100) ** 0.5) picture = []for r in range(ROWS):    row = []    for c in range(COLS):        n = drains(r, c)        if n >= 50:            row.append((20, 70, 180))       # a river: 500 km² or more        elif n >= 10:            row.append((90, 150, 230))      # a stream: 100 km² or more        else:            row.append(gray(heights[r][c]))    picture.append(row)show_grid(picture)

Chobhar: 34 squares, about 337 km²Border: 380 squares, about 3,762 km²

Picture

The rivers in the picture were never drawn: they are only the squares that a lot of land drains through. The dark line down the middle is the Bagmati; the one crossing the top right corner is the Sun Koshi, on its way east to the Koshi. By the border, your program says, the Bagmati carries the water of about 3,760 km². MERIT Hydro, with 90 m squares and the rivers seen from space, says 3,803 km² at the same place. Your 3 km squares, one rule and a fill of 41 squares got within about 1%. Not everything is as close: at Chobhar your program finds about 337 km², and MERIT Hydro about 590 km². With squares 3 km across, the hills round the valley are only a square or two thick, too rough to say exactly where its rim is.

The rule worked in the hills and went wrong where the ground flattens, at every scale: the map above, at 1.5 km, has thousands of pits on the plains. There, heights alone can’t say which way water goes.

Notice, too, the short streams on the plains running side by side and never joining. On an even slope each line of squares passes its own water straight on; real streams wander and join. The parallel lines are a sign of flat ground, which is why the picture marks a stream only once 100 km² drains through it.

How real tools do better

• Flats. Your 1 cm steps tilt a filled hollow towards its outlet; real tools tilt flats by the land around them, so water crosses them as it would.
• Rivers the radar can’t see. A river known from maps or from space can be burned into the heights, lowered a few meters so the water must follow it. MERIT Hydro used water seen from space and OpenStreetMap’s rivers.
• One direction is crude. On a broad slope water spreads out. Other rules share a square’s water among several lower neighbors, by how steep each is.

Try it yourself

• The steepest slope. O’Callaghan and Mark’s rule divides each fall by the distance: a diagonal neighbor is about 1.4 times as far away as one beside it. Change downhill to choose the steepest slope instead of the lowest square. Does the Bagmati move?
• Bring back the lake. Raise the square at Chobhar (places['Chobhar']) and its neighbors by 200 m before filling, and count how many squares the fill raises now, and how high. That is a lake in the Kathmandu Valley. How much land does it hold under water?
• An open question. The rule has no rain in it at all: the same heights in a desert and in Pokhara give the same rivers. What would you need to add to say how much water each river carries, and when it floods?

Three challenges, each checked on the page: write your function in the program above the checker (press Edit), then press Check my program. The examples show what your function should give; there are more tests you can’t see. In each, grid is a list of rows of heights, like heights.

★ Count the pits

Write count_pits(grid): how many squares, not on the edge, have no neighbor lower than themselves.

NEIGHBORS = [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)] def count_pits(grid):    # Count the squares, not on the edge, with no lower neighbor.    return 0 print(count_pits([[5, 5, 5], [5, 1, 5], [5, 5, 5]]))

★★ Whose rain passes here?

Write drains_through(grid, r, c): how many squares’ drops pass through the square (r, c), counting its own. Follow the drop from every square of the grid, by the lowest-neighbor rule, and count those that pass through (r, c). Don’t fill the pits: a drop stops in a pit.

NEIGHBORS = [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)] def drains_through(grid, r, c):    # How many squares' drops pass through (r, c), its own included?    return 0 print(drains_through([[9, 8, 7], [8, 5, 4], [7, 4, 1]], 2, 2))

★★★ How much would be under water?

Write lake_squares(grid): if every pit filled with water up to the level where it spills over, how many squares would be under water? Water leaves from the edge squares. A square is under water when the water stands higher than the square; a lake’s surface is flat (no 1 cm step).

import heapq NEIGHBORS = [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)] def lake_squares(grid):    # How many squares are under water once every pit fills to its brim?    return 0 print(lake_squares([[5, 5, 5, 5], [5, 1, 2, 5], [5, 5, 3, 5], [5, 5, 4, 5]]))

• Rivers of the World: the real rivers on a globe, from MERIT Hydro. Zoom in on Nepal and tap a river to see the land it drains.
• Grade 6 सामाजिक, नेपालका नदी प्रवाह क्षेत्र: the Koshi, Gandaki and Karnali systems, as the school books teach them.
• Grade 9, A photo is a table of numbers: a grid of numbers as a picture, the start of this topic’s map.
• Grade 10, Painting a picture, ring by ring: a search that spreads over a grid, as the sea does in Step 4.
• Grade 11, Flood fill: the paint bucket and Heaps: the most urgent first: the recursion of Step 5 and the heap of Step 4.
• What Does Water Do to Rock?: what the rivers your program found do to the mountains.
• How Much Electricity Can a Small Stream Make?: the height a river falls, turned into power.