The Grand Ethiopian Renaissance Dam filling, 2013 -> 2025¶
The GERD reservoir on the Blue Nile filled over a decade. We watch it from Landsat 8 (2013 onward), one dry-season (Nov-Feb) median composite per year — the dry season is the clear season over the Ethiopian highlands, so each yearly composite is crisp and cloud-free, and the year-to-year jump tells the filling story without the cloud strobing that a month-by-month series over this rainy reach produces.
In [ ]:
Copied!
import os
import shutil
import tempfile
import warnings
import matplotlib.pyplot as plt
import numpy as np
from cleopatra.glyphs.gridded.array_glyph import FrameLabel
from cleopatra.styling.colors import DATA_STYLES, resolve_colormap
from IPython.display import Image, display
from loguru import logger
from pyramids.dataset import Dataset
from pyramids.dataset.collection import DatasetCollection
from earthlens.core import EarthLens
warnings.filterwarnings("ignore")
logger.remove()
plt.rcParams["figure.dpi"] = 80
import os
import shutil
import tempfile
import warnings
import matplotlib.pyplot as plt
import numpy as np
from cleopatra.glyphs.gridded.array_glyph import FrameLabel
from cleopatra.styling.colors import DATA_STYLES, resolve_colormap
from IPython.display import Image, display
from loguru import logger
from pyramids.dataset import Dataset
from pyramids.dataset.collection import DatasetCollection
from earthlens.core import EarthLens
warnings.filterwarnings("ignore")
logger.remove()
plt.rcParams["figure.dpi"] = 80
Set the area of interest, the persistent output/cache directories, and the ffmpeg square-pad
filter. Raw downloads and derived frames are cached under out/gerd/ (never temp), so a re-run
only fetches what is missing.
In [ ]:
Copied!
from dotenv import find_dotenv, load_dotenv
load_dotenv(
find_dotenv(usecwd=True)
) # repo-root .env -> GEE_SERVICE_ACCOUNT / GEE_SERVICE_KEY
SERVICE_ACCOUNT = os.environ.get("GEE_SERVICE_ACCOUNT")
SERVICE_KEY = os.environ.get("GEE_SERVICE_KEY")
AOI = [34.95, 11.05, 35.25, 11.32]
OUT = os.path.join("out", "gerd")
raw_dir = os.path.join(OUT, "raw_yearly") # persistent per-year raw-download cache
rgb_dir = os.path.join(OUT, "rgb_yearly") # derived reflectance frames (persistent)
ndwi_dir = os.path.join(OUT, "ndwi_yearly")
os.makedirs(raw_dir, exist_ok=True)
for _d in (rgb_dir, ndwi_dir):
shutil.rmtree(_d, ignore_errors=True)
os.makedirs(_d, exist_ok=True)
VF = (
"scale=1200:1200:force_original_aspect_ratio=decrease,"
"pad=1200:1200:(ow-iw)/2:(oh-ih)/2:white"
)
from dotenv import find_dotenv, load_dotenv
load_dotenv(
find_dotenv(usecwd=True)
) # repo-root .env -> GEE_SERVICE_ACCOUNT / GEE_SERVICE_KEY
SERVICE_ACCOUNT = os.environ.get("GEE_SERVICE_ACCOUNT")
SERVICE_KEY = os.environ.get("GEE_SERVICE_KEY")
AOI = [34.95, 11.05, 35.25, 11.32]
OUT = os.path.join("out", "gerd")
raw_dir = os.path.join(OUT, "raw_yearly") # persistent per-year raw-download cache
rgb_dir = os.path.join(OUT, "rgb_yearly") # derived reflectance frames (persistent)
ndwi_dir = os.path.join(OUT, "ndwi_yearly")
os.makedirs(raw_dir, exist_ok=True)
for _d in (rgb_dir, ndwi_dir):
shutil.rmtree(_d, ignore_errors=True)
os.makedirs(_d, exist_ok=True)
VF = (
"scale=1200:1200:force_original_aspect_ratio=decrease,"
"pad=1200:1200:(ow-iw)/2:(oh-ih)/2:white"
)
Fetch one dry-season median per year¶
For each year we ask Google Earth Engine (via earthlens) for the median of every Landsat-8
Collection-2 L2 scene in that year's dry season (Nov -> Feb). Collection-2 DN is converted to true
surface reflectance (DN * 0.0000275 - 0.2) so the colour is vivid, and NDWI (Green-NIR)/(Green+NIR)
is derived for the water index.
In [ ]:
Copied!
frames_written = []
for y in range(2013, 2026): # 2013 dry season ... 2025/26 dry season (all in the past)
raw_path = os.path.join(raw_dir, f"gerd_{y}.tif")
if not os.path.exists(raw_path):
scene = EarthLens(
data_source="gee",
cadence="raw",
path=tempfile.mkdtemp(),
dataset="LANDSAT/LC08/C02/T1_L2",
variables=["SR_B4", "SR_B3", "SR_B2", "SR_B5"], # R, G, B, NIR
aoi=AOI,
start=f"{y}-11-01",
end=f"{y + 1}-02-28",
scale=30,
reducer="median",
export_via="url",
)
scene.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY)
try:
out = scene.download(progress_bar=False)
except Exception:
continue
if not out:
continue
shutil.copy(str(out[0]), raw_path)
# Collection-2 L2 surface reflectance = DN * 0.0000275 - 0.2 -> vivid true colour.
src = Dataset.read_file(raw_path)
a = np.clip(
np.asarray(src.read_array(), dtype="float32") * 0.0000275 - 0.2, 0.0, None
)
Dataset.create_from_array(arr=a[:3], geo=src.geotransform, epsg=src.epsg).to_file(
os.path.join(rgb_dir, f"gerd_{y}.tif")
)
ndwi = (a[1] - a[3]) / (a[1] + a[3] + 1e-6) # (Green - NIR)/(Green + NIR)
Dataset.create_from_array(
arr=ndwi[None, :, :], geo=src.geotransform, epsg=src.epsg
).to_file(os.path.join(ndwi_dir, f"gerd_{y}.tif"))
frames_written.append(y)
print(
f"{len(frames_written)} yearly composites: {frames_written[0]} -> {frames_written[-1]}"
)
frames_written = []
for y in range(2013, 2026): # 2013 dry season ... 2025/26 dry season (all in the past)
raw_path = os.path.join(raw_dir, f"gerd_{y}.tif")
if not os.path.exists(raw_path):
scene = EarthLens(
data_source="gee",
cadence="raw",
path=tempfile.mkdtemp(),
dataset="LANDSAT/LC08/C02/T1_L2",
variables=["SR_B4", "SR_B3", "SR_B2", "SR_B5"], # R, G, B, NIR
aoi=AOI,
start=f"{y}-11-01",
end=f"{y + 1}-02-28",
scale=30,
reducer="median",
export_via="url",
)
scene.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY)
try:
out = scene.download(progress_bar=False)
except Exception:
continue
if not out:
continue
shutil.copy(str(out[0]), raw_path)
# Collection-2 L2 surface reflectance = DN * 0.0000275 - 0.2 -> vivid true colour.
src = Dataset.read_file(raw_path)
a = np.clip(
np.asarray(src.read_array(), dtype="float32") * 0.0000275 - 0.2, 0.0, None
)
Dataset.create_from_array(arr=a[:3], geo=src.geotransform, epsg=src.epsg).to_file(
os.path.join(rgb_dir, f"gerd_{y}.tif")
)
ndwi = (a[1] - a[3]) / (a[1] + a[3] + 1e-6) # (Green - NIR)/(Green + NIR)
Dataset.create_from_array(
arr=ndwi[None, :, :], geo=src.geotransform, epsg=src.epsg
).to_file(os.path.join(ndwi_dir, f"gerd_{y}.tif"))
frames_written.append(y)
print(
f"{len(frames_written)} yearly composites: {frames_written[0]} -> {frames_written[-1]}"
)
Load the yearly frames into DatasetCollections with a real (yearly) time axis.
In [ ]:
Copied!
dc_rgb = (
DatasetCollection.from_files(rgb_dir, date_format="%Y", date_regex=r"\d{4}")
if frames_written
else None
)
dc_ndwi = (
DatasetCollection.from_files(ndwi_dir, date_format="%Y", date_regex=r"\d{4}")
if frames_written
else None
)
labels = [str(t)[:4] for t in dc_rgb.time] if dc_rgb else []
print(f"{len(labels)} years: {labels[0]} -> {labels[-1]}" if labels else "no scenes")
dc_rgb = (
DatasetCollection.from_files(rgb_dir, date_format="%Y", date_regex=r"\d{4}")
if frames_written
else None
)
dc_ndwi = (
DatasetCollection.from_files(ndwi_dir, date_format="%Y", date_regex=r"\d{4}")
if frames_written
else None
)
labels = [str(t)[:4] for t in dc_rgb.time] if dc_rgb else []
print(f"{len(labels)} years: {labels[0]} -> {labels[-1]}" if labels else "no scenes")
Representative years — true colour¶
In [ ]:
Copied!
if dc_rgb:
step = max(1, dc_rgb.time_length // 6)
idx = list(range(0, dc_rgb.time_length, step))[:6]
fig, ax = plt.subplots(1, len(idx), figsize=(2.4 * len(idx), 2.7))
ax = np.atleast_1d(ax)
for j, i in enumerate(idx):
a = np.asarray(
dc_rgb.iloc(i).read_array(), dtype="float32"
) # (3, H, W) reflectance
ax[j].imshow(np.clip(np.transpose(a, (1, 2, 0)) / 0.24, 0.0, 1.0))
ax[j].set_title(labels[i])
ax[j].set_xticks([])
ax[j].set_yticks([])
fig.suptitle("GERD reservoir - Landsat 8 true colour (selected years)", y=1.03)
plt.tight_layout()
plt.show()
else:
print("no scenes - set GEE_SERVICE_ACCOUNT / GEE_SERVICE_KEY and rerun")
if dc_rgb:
step = max(1, dc_rgb.time_length // 6)
idx = list(range(0, dc_rgb.time_length, step))[:6]
fig, ax = plt.subplots(1, len(idx), figsize=(2.4 * len(idx), 2.7))
ax = np.atleast_1d(ax)
for j, i in enumerate(idx):
a = np.asarray(
dc_rgb.iloc(i).read_array(), dtype="float32"
) # (3, H, W) reflectance
ax[j].imshow(np.clip(np.transpose(a, (1, 2, 0)) / 0.24, 0.0, 1.0))
ax[j].set_title(labels[i])
ax[j].set_xticks([])
ax[j].set_yticks([])
fig.suptitle("GERD reservoir - Landsat 8 true colour (selected years)", y=1.03)
plt.tight_layout()
plt.show()
else:
print("no scenes - set GEE_SERVICE_ACCOUNT / GEE_SERVICE_KEY and rerun")
True-colour animation, 2013 -> 2025¶
In [ ]:
Copied!
if dc_rgb:
gif = os.path.join(tempfile.mkdtemp(), "true_colour.gif")
mp4 = os.path.join(OUT, "true_colour_yearly.mp4")
# full_bleed: fills the square figure edge-to-edge, no title/margins (like the Aral timelapse).
glyph = dc_rgb.plot(
rgb_options={"rgb": [0, 1, 2], "surface_reflectance": 0.24},
figsize=(8, 8),
full_bleed=True,
)
glyph.animate(
labels,
interval=300,
frame_label=FrameLabel(location=[20, 72], color="white", size=30),
full_bleed=True,
)
glyph.fig.set_dpi(150) # figsize*dpi -> 1200 px
glyph.save_animation(gif, fps=2)
glyph.save_animation(mp4, fps=2, dpi=150, extra_args=["-vf", "scale=1200:1200"])
plt.close("all")
display(Image(filename=gif))
else:
print("no scenes to animate")
if dc_rgb:
gif = os.path.join(tempfile.mkdtemp(), "true_colour.gif")
mp4 = os.path.join(OUT, "true_colour_yearly.mp4")
# full_bleed: fills the square figure edge-to-edge, no title/margins (like the Aral timelapse).
glyph = dc_rgb.plot(
rgb_options={"rgb": [0, 1, 2], "surface_reflectance": 0.24},
figsize=(8, 8),
full_bleed=True,
)
glyph.animate(
labels,
interval=300,
frame_label=FrameLabel(location=[20, 72], color="white", size=30),
full_bleed=True,
)
glyph.fig.set_dpi(150) # figsize*dpi -> 1200 px
glyph.save_animation(gif, fps=2)
glyph.save_animation(mp4, fps=2, dpi=150, extra_args=["-vf", "scale=1200:1200"])
plt.close("all")
display(Image(filename=gif))
else:
print("no scenes to animate")
NDWI animation (water index), 2013 -> 2025¶
In [ ]:
Copied!
if dc_ndwi:
gif = os.path.join(tempfile.mkdtemp(), "ndwi.gif")
mp4 = os.path.join(OUT, "ndwi_yearly.mp4")
water_cmap = resolve_colormap(
next(iter(DATA_STYLES["river_discharge"].values()))["cmap"]
)
glyph = dc_ndwi.plot(band=0, cmap=water_cmap, figsize=(8, 8), full_bleed=True)
glyph.animate(
labels,
interval=300,
frame_label=FrameLabel(location=[20, 72], color="white", size=30),
full_bleed=True,
)
glyph.fig.set_dpi(150)
glyph.save_animation(gif, fps=2)
glyph.save_animation(mp4, fps=2, dpi=150, extra_args=["-vf", "scale=1200:1200"])
plt.close("all")
display(Image(filename=gif))
else:
print("no scenes - set GEE credentials and rerun")
if dc_ndwi:
gif = os.path.join(tempfile.mkdtemp(), "ndwi.gif")
mp4 = os.path.join(OUT, "ndwi_yearly.mp4")
water_cmap = resolve_colormap(
next(iter(DATA_STYLES["river_discharge"].values()))["cmap"]
)
glyph = dc_ndwi.plot(band=0, cmap=water_cmap, figsize=(8, 8), full_bleed=True)
glyph.animate(
labels,
interval=300,
frame_label=FrameLabel(location=[20, 72], color="white", size=30),
full_bleed=True,
)
glyph.fig.set_dpi(150)
glyph.save_animation(gif, fps=2)
glyph.save_animation(mp4, fps=2, dpi=150, extra_args=["-vf", "scale=1200:1200"])
plt.close("all")
display(Image(filename=gif))
else:
print("no scenes - set GEE credentials and rerun")
Reservoir water-surface area per year¶
In [ ]:
Copied!
if dc_ndwi:
px_km2 = (30 * 30) / 1e6
areas = [
float(
np.sum(np.asarray(dc_ndwi.iloc(i).read_array(), dtype="float32") > 0)
* px_km2
)
for i in range(dc_ndwi.time_length)
]
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(list(dc_ndwi.time), areas, marker="o", color="tab:blue")
ax.set(
ylabel="water surface area (km^2)",
title="GERD reservoir water-surface area, 2013 -> 2025",
)
ax.grid(alpha=0.3)
fig.autofmt_xdate()
plt.show()
if dc_ndwi:
px_km2 = (30 * 30) / 1e6
areas = [
float(
np.sum(np.asarray(dc_ndwi.iloc(i).read_array(), dtype="float32") > 0)
* px_km2
)
for i in range(dc_ndwi.time_length)
]
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(list(dc_ndwi.time), areas, marker="o", color="tab:blue")
ax.set(
ylabel="water surface area (km^2)",
title="GERD reservoir water-surface area, 2013 -> 2025",
)
ax.grid(alpha=0.3)
fig.autofmt_xdate()
plt.show()