Showcase — the GERD reservoir filling on the Blue Nile (2020 → 2024)¶
The Grand Ethiopian Renaissance Dam impounds the Blue Nile; its reservoir was filled in stages from 2020 to 2024. This notebook fetches one dry-season Sentinel-2 composite per year, reads them with DatasetCollection.read_multiple_files(), shows the fill in true colour with Dataset.plot(rgb_options=...), animates the open water with an NDWI stack, and measures the water-surface area each year.
What this notebook does¶
- Fetch one dry-season Sentinel-2 composite per year (2020–2024) — true colour
B4,B3,B2plusB8(NIR) — and derive an NDWI raster per year. - True-colour progression of the five years with
Dataset.plot(rgb_options=...). - Animated true colour — the same composites as a real-colour GIF.
- Animate the open water (NDWI) with
DatasetCollection.plot().save_animation(). - Measure the reservoir's water-surface area per year from NDWI.
NDWI = (Green − NIR)/(Green + NIR) > 0 marks open water.
Needs
GEE_SERVICE_ACCOUNT/GEE_SERVICE_KEYandpyramids-gis[viz].
Setup¶
First the imports and a bit of noise-suppression. pyramids provides Dataset / DatasetCollection (reading + plotting); earthlens provides the unified EarthLens entry point. tempfile / shutil / os handle the scratch directories the composites and GIFs are written into.
import os
import shutil
import tempfile
import warnings
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import Image, display
from loguru import logger
from pyramids.dataset import Dataset, DatasetCollection # pyramids' own IO + plotting
from earthlens import EarthLens
warnings.filterwarnings("ignore")
logger.remove()
plt.rcParams["figure.dpi"] = 80 # keep embedded outputs lean
Credentials¶
This notebook reads its credentials from a .env file at the repository root
(loaded with python-dotenv; your PATH is left untouched). Create that .env
with at least:
GEE_SERVICE_ACCOUNT="your-sa@your-project.iam.gserviceaccount.com"
GEE_SERVICE_KEY="/path/to/your/service-account-key.json"
from dotenv import find_dotenv, load_dotenv
load_dotenv(find_dotenv(usecwd=True))
SERVICE_ACCOUNT = os.environ[
"GEE_SERVICE_ACCOUNT"
] # Earth Engine service-account email (.env)
SERVICE_KEY = os.environ[
"GEE_SERVICE_KEY"
] # path to the service-account JSON key (.env)
1 · Fetch the yearly composites¶
For each dry season (Nov–Dec) from 2020 to 2024 we ask GEE for one median Sentinel-2 composite over the dam AOI — true-colour B4,B3,B2 plus B8 (NIR). Each composite is download()ed to its own temp directory, then copied into a shared rgb_dir under a year-stamped name so the whole series can be read back as one collection later.
Request parameters and scratch directories¶
AOI is the bounding box over the reservoir and C the Sentinel-2 collection id. rgb_dir holds the true-colour (+NIR) composites and ndwi_dir the derived NDWI rasters; both are temp directories that the per-year frames are copied into.
AOI = [34.95, 11.05, 35.25, 11.32]
C = "COPERNICUS/S2_SR_HARMONIZED"
rgb_dir, ndwi_dir = tempfile.mkdtemp(), tempfile.mkdtemp()
labels = []
Fetch, copy, and derive NDWI per year¶
For each year we build the EarthLens request, then download() it on its own line (no longer chained off the constructor) and take the first written path. The composite is copied into rgb_dir, an NDWI raster is derived from the green and NIR bands and written into ndwi_dir, and the year is recorded in labels.
for y in range(2020, 2025):
scene = EarthLens(
data_source="gee",
cadence="raw",
path=tempfile.mkdtemp(),
dataset=C,
variables=["B4", "B3", "B2", "B8"],
aoi=AOI,
start=f"{y}-11-01",
end=f"{y}-12-28",
scale=30,
reducer="median",
)
scene.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY)
outputs = scene.download(progress_bar=False)
p = str(outputs[0])
if not p:
continue
shutil.copy(p, os.path.join(rgb_dir, f"gerd_{y}.tif")) # true colour (+NIR)
a = np.asarray(Dataset.read_file(p).read_array(), dtype="float32")
ndwi = (a[1] - a[3]) / (a[1] + a[3] + 1e-6) # (Green - NIR)/(Green + NIR)
src = Dataset.read_file(p)
ndwi_ds = Dataset.create_from_array(
arr=ndwi[None, :, :], geo=src.geotransform, epsg=src.epsg
)
ndwi_ds.to_file(os.path.join(ndwi_dir, f"ndwi_{y}.tif")) # persist derived NDWI
labels.append(str(y))
Read the two series back as collections¶
DatasetCollection.read_multiple_files() reads each year-stamped folder into a single time-stacked collection — dc_rgb for the true-colour composites and dc_ndwi for the NDWI rasters — which the plotting and animation cells below consume.
dc_rgb = DatasetCollection.read_multiple_files(rgb_dir, date=False) if labels else None
dc_ndwi = (
DatasetCollection.read_multiple_files(ndwi_dir, date=False) if labels else None
)
print(f"{len(labels)} yearly composites:", labels)
2 · True-colour progression — Dataset.plot(rgb_options=...)¶
Each year rendered as a real-colour composite, side by side; the reservoir grows behind the dam across the five panels.
if dc_rgb:
n = dc_rgb.time_length
fig, ax = plt.subplots(1, n, figsize=(2.4 * n, 2.7))
for i, y in enumerate(labels):
dc_rgb.iloc(i).plot(
rgb_options={"rgb": [0, 1, 2], "percentile": 2}, ax=ax[i], title=y
)
ax[i].set_xticks([])
ax[i].set_yticks([])
fig.suptitle("GERD reservoir - Sentinel-2 true colour", y=1.03)
plt.tight_layout()
plt.show()
else:
print("no scenes - set GEE_SERVICE_ACCOUNT / GEE_SERVICE_KEY and rerun")
3 · Animated true colour — DatasetCollection.plot(rgb_options=...).save_animation()¶
The five dry-season composites as a real-colour animation — the reservoir visibly grows behind the dam. DatasetCollection.plot(rgb_options=...) composites the RGB bands per frame and animates them across the years.
if dc_rgb:
gif = os.path.join(tempfile.mkdtemp(), "true_colour.gif")
dc_rgb.plot(
rgb_options={"rgb": [0, 1, 2], "percentile": 2},
figsize=(5, 5),
title="GERD reservoir - Sentinel-2 true colour",
).save_animation(gif, fps=1.2)
plt.close("all") # drop the static builder frame; keep only the GIF
display(Image(filename=gif))
else:
print("no scenes to animate")
4 · Animated reservoir filling (NDWI) — DatasetCollection.plot().save_animation()¶
NDWI is a single-band water index (high over open water), so the collection animates it directly — the blue pool spreads year on year as the reservoir fills.
if dc_ndwi:
gif = os.path.join(tempfile.mkdtemp(), "ndwi.gif")
dc_ndwi.plot(
band=0, cmap="Blues", figsize=(5, 5), title="GERD reservoir - NDWI 2020 -> 2024"
).save_animation(gif, fps=1)
plt.close("all")
display(Image(filename=gif))
else:
print("no scenes - set GEE credentials and rerun")
5 · Reservoir surface area each year (NDWI)¶
NDWI = (Green − NIR)/(Green + NIR) > 0 marks open water; counting those pixels gives the surface area, which climbs year on year.
Count the water pixels per year¶
For each composite we recompute NDWI from the green and NIR bands and count the pixels above zero, scaling by the 30 m × 30 m pixel area to get square kilometres.
if dc_rgb:
px_km2 = (30 * 30) / 1e6
areas = []
for i in range(dc_rgb.time_length):
a = np.asarray(dc_rgb.iloc(i).read_array(), dtype="float32")
ndwi = (a[1] - a[3]) / (a[1] + a[3] + 1e-6)
areas.append(float(np.sum(ndwi > 0) * px_km2))
print("reservoir area (km^2):", {l: round(v) for l, v in zip(labels, areas)})
Plot area by year¶
The per-year surface areas as a bar chart, each bar annotated with its value in square kilometres.
if dc_rgb:
fig, ax = plt.subplots(figsize=(7, 4))
ax.bar(labels, areas, color="tab:blue")
for i, v in enumerate(areas):
ax.text(i, v, f"{v:.0f}", ha="center", va="bottom")
ax.set(ylabel="water surface area (km^2)", title="GERD reservoir area by year")
plt.show()
6 · Annual fill rate¶
Differencing the yearly surface areas gives how much new water the reservoir gained each year — the pace of the staged filling.
if dc_rgb and len(areas) > 1:
rates = [areas[i] - areas[i - 1] for i in range(1, len(areas))]
yl = [f"{labels[i-1]}->{labels[i]}" for i in range(1, len(labels))]
fig, ax = plt.subplots(figsize=(8, 4))
ax.bar(yl, rates, color="teal")
for i, v in enumerate(rates):
ax.text(i, v, f"+{v:.0f}", ha="center", va="bottom")
ax.set(ylabel="area gained (km^2/yr)", title="GERD reservoir annual fill rate")
plt.show()
print("annual gain (km^2):", dict(zip(yl, [round(r) for r in rates])))
else:
print("no areas to difference")
Recap¶
The yearly series is read from one folder with DatasetCollection.read_multiple_files(); the progression renders in true colour with Dataset.plot(rgb_options=...) and animates as a real-colour GIF, the NDWI stack animates with DatasetCollection.plot().save_animation(), and counting NDWI>0 pixels turns it into a surface-area-per-year curve.
Try it yourself¶
- Point
AOIat another reservoir, or animate a different index. - See the GEE backend reference.