True-colour approach 2 — a Landsat wet-year vs 2026 baseline¶

Sentinel-2 only reaches back to ~2017, so a 2025-vs-2026 pair has little contrast. Landsat 8/9 (surface reflectance, 30 m) reaches back to 2013, so we can anchor 2026 against a genuinely wet year. July 2021 brought record floods to the Rhine/Meuse basins — a natural high-water bookend. We put a Landsat true-colour composite of summer 2021 next to summer 2026 over the Danube at Budapest.

Landsat's 16-day revisit means a single sensor over a short window catches cloud, so we widen the windows to a full summer and merge Landsat 8 + Landsat 9 — picking, per pixel, the clearer (darker) of the two composites — to beat the cloud.

Needs GEE_SERVICE_ACCOUNT / GEE_SERVICE_KEY and pyramids-gis[viz].

Setup¶

pyramids reads and plots the GeoTIFFs (Dataset / DatasetCollection); earthlens provides the EarthLens entry point. Earth Engine credentials come from a repo-root .env.

In [ ]:

import os
import shutil
import tempfile
import warnings
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
from loguru import logger
from pyramids.dataset import Dataset

from earthlens.core import EarthLens

warnings.filterwarnings("ignore")
logger.remove()
plt.rcParams["figure.dpi"] = 80

from dotenv import find_dotenv, load_dotenv

load_dotenv(find_dotenv(usecwd=True))
SERVICE_ACCOUNT = os.environ["GEE_SERVICE_ACCOUNT"]
SERVICE_KEY = os.environ["GEE_SERVICE_KEY"]

import os import shutil import tempfile import warnings from pathlib import Path import matplotlib.pyplot as plt import numpy as np from loguru import logger from pyramids.dataset import Dataset from earthlens.core import EarthLens warnings.filterwarnings("ignore") logger.remove() plt.rcParams["figure.dpi"] = 80 from dotenv import find_dotenv, load_dotenv load_dotenv(find_dotenv(usecwd=True)) SERVICE_ACCOUNT = os.environ["GEE_SERVICE_ACCOUNT"] SERVICE_KEY = os.environ["GEE_SERVICE_KEY"]

The reach, the two years, and the two sensors¶

Landsat Collection-2 Level-2 surface reflectance (SR_B4/SR_B3/SR_B2 = red/green/blue). We take a cloud-median summer composite from both Landsat 8 and Landsat 9 for a wet year (2021) and for 2026 over the Budapest reach, at Landsat's native 30 m. (2026 stops in early August — the latest data available.)

In [ ]:

AOI = [18.98, 47.40, 19.14, 47.58]
LANDSATS = ["LANDSAT/LC08/C02/T1_L2", "LANDSAT/LC09/C02/T1_L2"]
SCALE = 30
YEARS = {
    "2021 (wet)": ("2021-06-01", "2021-08-31"),
    "2026 (drought)": ("2026-06-01", "2026-08-05"),
}
OUT = Path("out") / "tc2_landsat_wetyear"
OUT.mkdir(parents=True, exist_ok=True)

AOI = [18.98, 47.40, 19.14, 47.58] LANDSATS = ["LANDSAT/LC08/C02/T1_L2", "LANDSAT/LC09/C02/T1_L2"] SCALE = 30 YEARS = { "2021 (wet)": ("2021-06-01", "2021-08-31"), "2026 (drought)": ("2026-06-01", "2026-08-05"), } OUT = Path("out") / "tc2_landsat_wetyear" OUT.mkdir(parents=True, exist_ok=True)

Fetch a summer composite per (year, sensor)¶

A median Landsat composite of the three visible bands for each year from each of Landsat 8 and 9, cached under out/.

In [ ]:

paths = {}
for tag, (start, end) in YEARS.items():
    for ds in LANDSATS:
        sensor = ds.split("/")[1]  # LC08 / LC09
        p = OUT / f"landsat_{tag[:4]}_{sensor}.tif"
        if not p.exists():
            job = EarthLens(
                data_source="gee",
                cadence="raw",
                path=tempfile.mkdtemp(),
                dataset=ds,
                variables=["SR_B4", "SR_B3", "SR_B2"],
                aoi=AOI,
                start=start,
                end=end,
                scale=SCALE,
                reducer="median",
            )
            job.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY)
            out = job.download(progress_bar=False)
            if not out:
                continue
            shutil.copy(str(out[0]), p)
        paths[(tag, sensor)] = p
{f"{t}/{s}": v.name for (t, s), v in paths.items()}

paths = {} for tag, (start, end) in YEARS.items(): for ds in LANDSATS: sensor = ds.split("/")[1] # LC08 / LC09 p = OUT / f"landsat_{tag[:4]}_{sensor}.tif" if not p.exists(): job = EarthLens( data_source="gee", cadence="raw", path=tempfile.mkdtemp(), dataset=ds, variables=["SR_B4", "SR_B3", "SR_B2"], aoi=AOI, start=start, end=end, scale=SCALE, reducer="median", ) job.authenticate(service_account=SERVICE_ACCOUNT, service_key=SERVICE_KEY) out = job.download(progress_bar=False) if not out: continue shutil.copy(str(out[0]), p) paths[(tag, sensor)] = p {f"{t}/{s}": v.name for (t, s), v in paths.items()}

Merge the two sensors per year — keep the clearer pixel¶

Clouds are bright, so for each pixel we keep whichever sensor's composite is darker (less cloud). That gives one cloud-reduced true-colour composite per year.

In [ ]:

def merge_clearer(a, b):
    a = np.asarray(a, dtype="float32")
    b = np.asarray(b, dtype="float32")
    r, c = min(a.shape[1], b.shape[1]), min(a.shape[2], b.shape[2])
    a, b = a[:, :r, :c], b[:, :r, :c]
    use_a = a.mean(0) <= b.mean(0)  # keep the darker (clearer) pixel
    return np.where(use_a[None], a, b)


merged = {}
for tag in YEARS:
    comps = [
        np.asarray(
            Dataset.read_file(str(paths[(tag, s)]).replace("\\", "/")).read_array(),
            dtype="float32",
        )
        for s in ("LC08", "LC09")
        if (tag, s) in paths
    ]
    merged[tag] = comps[0] if len(comps) == 1 else merge_clearer(comps[0], comps[1])
{tag: m.shape for tag, m in merged.items()}

def merge_clearer(a, b): a = np.asarray(a, dtype="float32") b = np.asarray(b, dtype="float32") r, c = min(a.shape[1], b.shape[1]), min(a.shape[2], b.shape[2]) a, b = a[:, :r, :c], b[:, :r, :c] use_a = a.mean(0) <= b.mean(0) # keep the darker (clearer) pixel return np.where(use_a[None], a, b) merged = {} for tag in YEARS: comps = [ np.asarray( Dataset.read_file(str(paths[(tag, s)]).replace("\\", "/")).read_array(), dtype="float32", ) for s in ("LC08", "LC09") if (tag, s) in paths ] merged[tag] = comps[0] if len(comps) == 1 else merge_clearer(comps[0], comps[1]) {tag: m.shape for tag, m in merged.items()}

Wet year vs drought year, side by side¶

Landsat SR values are scaled integers, so a percentile stretch handles the display. The 2021 flood summer against the 2026 drought summer over the same reach.

In [ ]:

fig, axes = plt.subplots(1, 2, figsize=(12, 5.5))
for ax, tag in zip(axes, YEARS):
    a = merged[tag]
    rgb = np.clip(
        np.stack([a[0], a[1], a[2]], axis=-1) / np.nanpercentile(a[:3], 98), 0, 1
    )
    ax.imshow(rgb)
    ax.set_title(f"Danube at Budapest - {tag}")
    ax.set_xticks([])
    ax.set_yticks([])
fig.suptitle(
    "Landsat true colour (LC08+LC09 merged) - wet year vs drought year", y=1.02
)
plt.tight_layout()
plt.show()

fig, axes = plt.subplots(1, 2, figsize=(12, 5.5)) for ax, tag in zip(axes, YEARS): a = merged[tag] rgb = np.clip( np.stack([a[0], a[1], a[2]], axis=-1) / np.nanpercentile(a[:3], 98), 0, 1 ) ax.imshow(rgb) ax.set_title(f"Danube at Budapest - {tag}") ax.set_xticks([]) ax.set_yticks([]) fig.suptitle( "Landsat true colour (LC08+LC09 merged) - wet year vs drought year", y=1.02 ) plt.tight_layout() plt.show()

Notes¶

• Landsat's 30 m is coarser than Sentinel-2's 10 m, but its multi-decade archive is the point — pick any wet reference year (2013, 2016, 2021).
• Merging LC08+LC09 and widening the window are what make the drought-year composite usable; a single sensor over a few weeks is often too cloudy.