True-colour approach 4 — mapping newly exposed riverbed at Novi Sad¶
The idea: freshly exposed sand is bright in true colour, open water is dark, so a pixel that goes from dark (water) to bright (sand) between spring and July is newly exposed riverbed. But brightness alone is fooled by farmland (green April fields → bright harvested July fields read as "exposed" too).
The fix is to bound the search to the river channel: we use NDWI (a one-line water index from the green and NIR bands) to outline the spring water extent, then map, within that channel, the pixels that have since dried out — displayed on the July true-colour image. NDWI only defines "where the river was"; the drying itself is what we map.
Needs
GEE_SERVICE_ACCOUNT/GEE_SERVICE_KEYandpyramids-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.
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 and the two composites¶
The Danube at Novi Sad (braided, sandbar-rich). We take an early-season ("wet") and a July ("dry")
median composite with the visible bands plus NIR (B8) so we can compute NDWI, at 15 m.
AOI = [19.78, 45.20, 19.95, 45.30] # Danube at Novi Sad
S2 = "COPERNICUS/S2_SR_HARMONIZED"
SCALE = 15 # 15 m keeps the export under GEE's 50 MB cap
COMPOSITES = {"wet": ("2026-04-01", "2026-05-15"), "dry": ("2026-07-01", "2026-07-24")}
OUT = Path("out") / "tc4_novisad_exposed"
OUT.mkdir(parents=True, exist_ok=True)
Fetch the wet and dry composites (RGB + NIR, cached)¶
paths = {}
for tag, (start, end) in COMPOSITES.items():
p = OUT / f"novisad_{tag}.tif"
if not p.exists():
job = EarthLens(
data_source="gee",
cadence="raw",
path=tempfile.mkdtemp(),
dataset=S2,
variables=["B4", "B3", "B2", "B8"],
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] = p
{k: v.name for k, v in paths.items()}
Derive the exposed-riverbed mask (NDWI-bounded)¶
NDWI = (green - NIR) / (green + NIR); it is positive over open water. The spring channel is where NDWI was positive in the wet composite; exposed riverbed is that channel where NDWI has since gone negative (no longer water). Bounding by the spring channel is what keeps farmland out of the result.
wet = np.asarray(Dataset.read_file(str(paths["wet"])).read_array(), dtype="float32")
dry = np.asarray(Dataset.read_file(str(paths["dry"])).read_array(), dtype="float32")
def ndwi(a): # bands: 0=B4 red, 1=B3 green, 2=B2 blue, 3=B8 NIR
return (a[1] - a[3]) / (a[1] + a[3] + 1e-6)
channel = ndwi(wet) > 0 # water in spring
still_water = ndwi(dry) > 0 # water in July
exposed = channel & ~still_water # spring channel that has since dried
print(
"spring channel px:",
int(channel.sum()),
"| exposed px:",
int(exposed.sum()),
"| exposed fraction of channel:",
round(float(exposed.sum() / max(channel.sum(), 1)), 3),
)
Map it: July true colour with the exposed bed highlighted¶
The July composite in true colour, with the newly exposed riverbed overlaid in orange — confined to the channel, so farmland no longer contaminates it.
rgb = np.clip(
np.stack([dry[0], dry[1], dry[2]], axis=-1) / np.nanpercentile(dry[:3], 98), 0, 1
)
overlay = np.zeros((*exposed.shape, 4), dtype="float32")
overlay[exposed] = (1.0, 0.5, 0.0, 0.9) # orange where exposed
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].imshow(rgb)
axes[0].set_title("Danube at Novi Sad - July true colour")
axes[1].imshow(rgb)
axes[1].imshow(overlay)
axes[1].set_title("newly exposed riverbed (orange)")
for ax in axes:
ax.set_xticks([])
ax.set_yticks([])
plt.tight_layout()
plt.show()
Notes¶
- The honest lesson: pure-RGB brightness can't separate riverbed from harvested farmland — both go dark-to-bright. Bounding the search to the spring water channel (one NDWI line) fixes that.
- The sediment story is still shown in true colour; NDWI only supplies the channel outline.