Is the reservoir my town depends on rising or falling — and who relies on it?

A reservoir is a big artificial lake held back by a dam. It’s where a town stores its drinking and irrigation water. Two numbers describe its health: how high the water sits (its surface elevation, in metres) and how much area it covers from above. A satellite called SWOT now measures both directly from orbit, for lakes all over the world. So you can watch the water your town drinks from rise or fall even where there’s no public gauge on the shore.

Which data to use, and how

Use this dataset: SWOT lakes (short name SWOT_L2_HR_LakeSP_D). SWOT (a NASA/CNES radar mission launched December 2022) is the one instrument that gives you a reservoir’s water-surface elevation (wse, in metres) and its surface area (area_total, in km²) directly. Before it existed, you had to stitch together cloudy optical photos or hope a gauge station was nearby. The catch: SWOT files come as zipped shapefiles (a map-data format you open with the geopandas library), not simple grids, so this one needs more handling than most.

Then do five steps (the runnable code further down does all of this — this is just the idea):

1. Pick your reservoir (a small latitude/longitude box around it) and a handful of dates spread across the years you care about.
2. Grab each pass. For each date, download the SWOT lake file that covers your box. A “pass” is one flyover; SWOT revisits roughly every 10–21 days.
3. Find your lake and clean it up. Keep only real, fully-seen water bodies inside your box and throw away junk readings (more on that in the gotchas below).
4. Read the two numbers. Record the elevation (wse) and area for the biggest lake in your box on each date.
5. See the trend. Line the dates up and look: is the water level going up, down, or holding steady? Compare against the long-term history and count who lives nearby.

That’s the whole method. Below, “How a scientist answers this” names the exact tools, and the Run it cell does it live on synthetic data so you can see each step work.

What you can find out

• Whether the reservoir’s level is rising or falling. Track SWOT elevation (wse) across passes to see the trend, including in places with no public gauge.
• How its current size compares to the long-term normal. Overlay the free JRC Global Surface Water baseline (1984–2021) to see how today’s extent sits against decades of history.
• How much surface area it has now. SWOT reports the observed lake area on each pass.
• Roughly how many people depend on it. Overlay free WorldPop population data on the towns or metro it serves, and name the area with geoBoundaries.
• Whether a drought or a wet year moved it. The SWOT record picks up the seasonal draw-down and multi-year swings. For Lake Mead, this workflow shows the reservoir held steady and partly recovered across 2023–2025, matching the Colorado River’s better snow years. (Don’t assume the answer. Let the data tell you.)

What it can’t tell you

• The exact stored volume (acre-feet or m³). SWOT gives surface elevation and area, not volume. Turning those into “how much water is in there” needs the reservoir’s bathymetry, the depth-to-capacity curve that only the dam operator has.
• A clean reading from every single pass. SWOT only sees the strip of ground it flies over, so on some passes it catches just part of a big lake. Those partial passes report a too-small area and a noisy elevation, and you have to filter them out before trusting any one point.
• The absolute level against a local gauge. SWOT’s elevation is measured from a global reference surface (a geoid datum, roughly “average sea level, modelled worldwide”), which can sit tens of metres off a local town gauge’s zero point. Use the change over time, not the raw number.
• Who is legally served, or who holds the water rights. WorldPop just counts residents in an area. Actual served population follows utility service boundaries and downstream water allocation, not the lake’s footprint.
• Water quality. Level and area say nothing about salinity, sediment, or contamination.

Gotchas to watch for

• Each pass may see only part of the lake (partial coverage). SWOT crosses a slice of a large lake on most flyovers. If you read the level off a half-seen lake you’ll think it shrank when it didn’t. The fix: keep only passes where the observed area is large (well-observed), and drop the rest. The code filters on area > 300 km² for Lake Mead, a threshold you’d tune per reservoir.
• It’s a datum, not a gauge. The elevation (wse) is referenced to that global geoid surface, so it can differ from a local gauge by tens of metres. Always track the change between dates, never the absolute number.
• The record is young. SWOT science data starts mid-2023, so you get a few years of trend, not decades. For the long view, lean on the JRC baseline below.
• Most rows in a file are empty placeholders. A SWOT lake file lists every lake it might see from a prior catalogue (the Prior Lake Database), and about 70% of those rows are unobserved on any given pass. They carry null geometry and a fill value of -999999999999. Drop the empty ones before doing anything; the code does this for you.
• Level isn’t volume, and volume isn’t supply. Even a confident rise in level doesn’t equal more usable water. That depends on the operator’s capacity curve, allocation rules, and downstream demand.
• Small reservoirs slip through. SWOT resolves lakes down to roughly a few hundred metres across; farm ponds and narrow canyon reservoirs may be missed entirely.

The real-data code

The Run it cell above runs this method on synthetic data with no login. Below is the whole recipe against the real archive: search SWOT passes, open the zipped lake shapefiles, pull wse/area for the biggest well-observed lake per window, read off the level change, and plot it — then overlay the free, no-login population and boundary layers. It needs a free Earthdata Login.

Verified locally. SWOT lake products are zipped shapefiles; read them with geopandas (zip://). Each lake has a stable lake_id from the Prior Lake Database; area_total and wse are the surface area (km²) and elevation (m). Fill value is -999999999999, and ~70% of rows in a file are unobserved prior lakes with null geometry. For Lake Mead the id is 7720003253.

import earthaccess, os, requests, numpy as np
import geopandas as gpd, pandas as pd, rasterio
import matplotlib.pyplot as plt
from shapely.geometry import box, Point
from rasterio.windows import from_bounds

earthaccess.login(strategy="netrc")
sess = earthaccess.get_requests_https_session()   # authed, to fetch single assets

aoi = (-114.95, 35.95, -113.9, 36.65)             # Lake Mead, USA, as (W, S, E, N)
os.makedirs("swot", exist_ok=True)

# 1-3. Search SWOT passes per window, open each zipped Prior-lake shapefile, drop the
#       empty placeholder rows + partial passes, keep the biggest well-observed lake.
def reservoir_series(time_windows):
    rows = []
    for win in time_windows:
        results = earthaccess.search_data(short_name="SWOT_L2_HR_LakeSP_D",
                                          bounding_box=aoi, temporal=win, count=60)
        best = None
        for g in results:
            prior = [l for l in g.data_links() if "_Prior_" in l and l.endswith(".zip")]
            if not prior:
                continue
            fp = f"swot/{os.path.basename(prior[0])}"
            if not os.path.exists(fp):
                open(fp, "wb").write(sess.get(prior[0], timeout=120).content)
            gdf = gpd.read_file(f"zip://{fp}")
            gdf = gdf[~gdf.geometry.isna()]                 # drop empty placeholder rows
            gdf["area"] = pd.to_numeric(gdf["area_total"], errors="coerce")
            gdf["wse_m"] = pd.to_numeric(gdf["wse"], errors="coerce")
            # QC: real water-bodies in the AOI, observed area > 300 km² (drop partial passes)
            obs = gdf[gdf.intersects(box(*aoi)) & (gdf["area"] > 300) & (gdf["wse_m"] > -1000)]
            if len(obs):
                b = obs.loc[obs["area"].idxmax()]
                if best is None or b["area"] > best["area"]:
                    best = {"date": pd.to_datetime(win[0]), "lake_id": b["lake_id"],
                            "area_km2": float(b["area"]), "wse_m": float(b["wse_m"])}
        if best:
            rows.append(best)
    return pd.DataFrame(rows).sort_values("date").reset_index(drop=True)

windows = [("2023-08-01", "2023-08-20"), ("2024-01-01", "2024-01-25"),
           ("2024-06-01", "2024-06-25"), ("2024-12-01", "2024-12-25"),
           ("2025-04-01", "2025-04-25")]
ts = reservoir_series(windows)            # verified: area ~380-406 km², wse ~359-363 m

# 4. The answer: change in surface elevation, first well-observed pass to last.
change = ts["wse_m"].iloc[-1] - ts["wse_m"].iloc[0]
slope = np.polyfit((ts["date"] - ts["date"].iloc[0]).dt.days, ts["wse_m"], 1)[0] * 365
verdict = "RISING" if change > 0.2 else "FALLING" if change < -0.2 else "HOLDING STEADY"
print(f"Lake Mead is {verdict}: {change:+.2f} m over the record "
      f"({slope:+.2f} m/yr), area now {ts['area_km2'].iloc[-1]:.0f} km².")

# 5. Plot both numbers across the SWOT passes.
fig, ax = plt.subplots()
ax.plot(ts["date"], ts["wse_m"], "o-", color="tab:blue", label="elevation (m)")
ax.set_ylabel("SWOT water-surface elevation (m, geoid)")
ax2 = ax.twinx()
ax2.plot(ts["date"], ts["area_km2"], "s--", color="tab:green", label="area (km²)")
ax2.set_ylabel("surface area (km²)")
ax.set_title("Lake Mead — SWOT level & area"); fig.autofmt_xdate()
fig.tight_layout(); plt.show()

# --- Who depends on it: free WorldPop + geoBoundaries (no NASA login) ---
served = (-115.4, 35.9, -114.9, 36.4)      # Las Vegas valley (one served metro)
meta = requests.get("https://www.worldpop.org/rest/data/pop/wpic1km?iso3=USA").json()
pop_url = next(f for f in meta["data"][-1]["files"] if f.endswith(".tif"))
open("usa_pop_1km.tif", "wb").write(requests.get(pop_url).content)
with rasterio.open("usa_pop_1km.tif") as src:
    pop = src.read(1, window=from_bounds(*served, transform=src.transform)).astype("float64")
    pop[pop == src.nodata] = np.nan
adm = gpd.read_file(requests.get(
    "https://www.geoboundaries.org/api/current/gbOpen/USA/ADM2/").json()["gjDownloadURL"])
county = adm[adm.contains(Point(-115.14, 36.17))].iloc[0]["shapeName"]
print(f"~{np.nansum(pop):,.0f} people in the served metro ({county} County).")

# Cross-check: SWOT wse is a global-geoid datum, not a gauge — trust the *change*, not the
# raw metres; sanity-check the extent against the free JRC 1984–2021 baseline before believing it.

Where the data comes from

This is a multi-agency question, joined together on your own machine. The current level and area come from NASA/CNES SWOT, served by NASA’s PO.DAAC archive (one free Earthdata Login). The decades-long JRC Global Surface Water baseline comes from the EU’s Joint Research Centre. The “who depends on it” side comes from WorldPop (University of Southampton) for population and geoBoundaries (William & Mary) for naming the county or district. Those last three are all free, no-login layers stitched in client-side.

Sources

• SWOT Lake Single-Pass (PO.DAAC): https://podaac.jpl.nasa.gov/dataset/SWOT_L2_HR_LakeSP_2.0
• SWOT mission & hydrology: https://swot.jpl.nasa.gov/
• SWOT Prior Lake Database (lake_id): https://hydroweb.next.theia-land.fr/
• JRC Global Surface Water (free baseline): https://global-surface-water.appspot.com/
• WorldPop population (free, no login): https://www.worldpop.org/
• geoBoundaries (free CC-BY admin boundaries): https://www.geoboundaries.org/