Showcase — Gaza, satellite change detection (Khan Younis, 2023 → 2024)¶

Before/after satellite imagery is a standard humanitarian damage-assessment tool. This notebook builds a two-year Sentinel-2 time series over Khan Younis (southern Gaza) as one cloud-masked median composite per time window — so every frame is clean without hand-picking clear dates. It reads each series off disk with pyramids' DatasetCollection.read_multiple_files(), renders the scenes in true colour with Dataset.plot(rgb_options=...), and animates the derived NDVI with the collection's own save_animation(). Presented factually; no interpretation beyond what the pixels show.

What this notebook does¶

1. Composite a Sentinel-2 series — one cloud-masked median composite per ~6-week window — from the true-colour bands B4,B3,B2 and the NIR/red pair B8,B4, into two folders, computing an NDVI raster per window.
2. Read each folder as a collection with DatasetCollection.read_multiple_files() (no manual file lists).
3. True-colour progression — every window rendered as an RGB composite with Dataset.plot(rgb_options=...).
4. Animated true colour — the same scenes as a real-colour GIF.
5. NDVI before / after / change maps with Dataset.plot().
6. Animated NDVI across the windows with DatasetCollection.plot().save_animation().
7. NDVI timeline (matplotlib chart).

NDVI = (NIR − Red)/(NIR + Red): high over vegetation, near zero over bare ground.

Each frame is composited with raw ee — Earth Engine's per-pixel cloud mask is not exposed by the EarthLens GEE facade. Needs GEE_SERVICE_ACCOUNT / GEE_SERVICE_KEY and pyramids-gis[viz].

Setup¶

First the imports and a bit of noise-suppression. pyramids provides Dataset and DatasetCollection (reading + plotting); earthlens provides the unified EarthLens entry point.

In [ ]:

import os
import shutil
import tempfile
import urllib.request
import warnings

import ee
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from IPython.display import Image, display
from loguru import logger
from pyramids.dataset import Dataset, DatasetCollection  # pyramids' own IO + plotting

from earthlens.gee.auth import EarthEngineAuth  # service-account -> ee.Initialize

warnings.filterwarnings("ignore")
logger.remove()  # silence earthlens / pyramids / GDAL log chatter
plt.rcParams["figure.dpi"] = 80  # keep embedded outputs lean

import os import shutil import tempfile import urllib.request import warnings import ee import matplotlib.pyplot as plt import numpy as np import pandas as pd from IPython.display import Image, display from loguru import logger from pyramids.dataset import Dataset, DatasetCollection # pyramids' own IO + plotting from earthlens.gee.auth import EarthEngineAuth # service-account -> ee.Initialize warnings.filterwarnings("ignore") logger.remove() # silence earthlens / pyramids / GDAL log chatter plt.rcParams["figure.dpi"] = 80 # keep embedded outputs lean

Credentials¶

The 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:

GEE_SERVICE_ACCOUNT="your-sa@your-project.iam.gserviceaccount.com"
GEE_SERVICE_KEY="/path/to/your/service-account-key.json"

In [ ]:

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
SERVICE_KEY = os.environ["GEE_SERVICE_KEY"]  # path to the service-account JSON key
if any(x is None for x in (SERVICE_ACCOUNT, SERVICE_KEY)):
    raise ValueError("Missing credentials in .env")

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 SERVICE_KEY = os.environ["GEE_SERVICE_KEY"] # path to the service-account JSON key if any(x is None for x in (SERVICE_ACCOUNT, SERVICE_KEY)): raise ValueError("Missing credentials in .env")

1 · Build the cloud-masked Sentinel-2 series¶

Request parameters¶

The AOI is a small lon/lat box over Khan Younis. WINDOWS is seven ~6-week time windows spanning 2023–2024; each one becomes a single cloud-masked median composite. We also create the two scratch folders the true-colour and derived NDVI composites land in.

In [ ]:

AOI = [34.28, 31.32, 34.33, 31.37]  # Khan Younis, southern Gaza (lon/lat box)
COLL = "COPERNICUS/S2_SR_HARMONIZED"
# (label, window-start, window-end) — one cloud-masked median composite per window
WINDOWS = [
    ("2023-01", "2023-01-01", "2023-02-15"),
    ("2023-05", "2023-05-01", "2023-06-15"),
    ("2023-09", "2023-09-01", "2023-10-15"),
    ("2023-12", "2023-12-01", "2024-01-15"),
    ("2024-03", "2024-03-01", "2024-04-15"),
    ("2024-05", "2024-05-01", "2024-06-15"),
    ("2024-09", "2024-09-01", "2024-10-15"),
]
rgb_dir, ndvi_dir = tempfile.mkdtemp(), tempfile.mkdtemp()

AOI = [34.28, 31.32, 34.33, 31.37] # Khan Younis, southern Gaza (lon/lat box) COLL = "COPERNICUS/S2_SR_HARMONIZED" # (label, window-start, window-end) — one cloud-masked median composite per window WINDOWS = [ ("2023-01", "2023-01-01", "2023-02-15"), ("2023-05", "2023-05-01", "2023-06-15"), ("2023-09", "2023-09-01", "2023-10-15"), ("2023-12", "2023-12-01", "2024-01-15"), ("2024-03", "2024-03-01", "2024-04-15"), ("2024-05", "2024-05-01", "2024-06-15"), ("2024-09", "2024-09-01", "2024-10-15"), ] rgb_dir, ndvi_dir = tempfile.mkdtemp(), tempfile.mkdtemp()

Authenticate Earth Engine, then the cloud-masked composite recipe¶

Per-pixel cloud masking is not exposed by the EarthLens GEE facade, so we initialise Earth Engine through earthlens's EarthEngineAuth (the same service-account path the facade uses) and then composite with raw ee. composite_tif() builds one cloud-masked median composite over a date window in five steps:

1. Filter the S2_SR_HARMONIZED collection to the AOI and the date window.
2. Drop the worst scenes up front by metadata — CLOUDY_PIXEL_PERCENTAGE < 60.
3. Per-pixel cloud mask each remaining scene via the Scene-Classification (SCL) band — removing cloud-shadow (3), medium/high cloud (8/9) and cirrus (10).
4. Median-composite the masked stack: median() ignores masked pixels, so each output pixel is the median of its clear observations across the window — transient clouds simply drop out.
5. Clip to the AOI and download the composite as a GeoTIFF via getDownloadURL.

ndvi_from_nir() then turns the downloaded NIR/red pair into an NDVI raster on disk, reusing the source geotransform/EPSG so the result stays georeferenced.

In [ ]:

EarthEngineAuth.initialize(SERVICE_ACCOUNT, SERVICE_KEY)  # one-time ee.Initialize
GEOM = ee.Geometry.Rectangle(AOI)  # [W, S, E, N]


def mask_s2_clouds(img):
    """Drop cloud / shadow / cirrus pixels via the Sentinel-2 SCL band.

    SCL classes removed: 3 (cloud shadow), 8 (cloud, medium prob),
    9 (cloud, high prob), 10 (thin cirrus).
    """
    scl = img.select("SCL")
    clear = scl.neq(3).And(scl.neq(8)).And(scl.neq(9)).And(scl.neq(10))
    return img.updateMask(clear)


def composite_tif(start, end, bands, scale=10):
    """Cloud-masked median S2 composite over [start, end) for the AOI; return the GeoTIFF path."""
    collection = (
        ee.ImageCollection(COLL)
        .filterBounds(GEOM)
        .filterDate(start, end)
        .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 60))
        .map(mask_s2_clouds)
        .select(bands)
    )
    composite = collection.median().clip(GEOM)
    url = composite.getDownloadURL(
        {"region": GEOM, "scale": scale, "format": "GEO_TIFF"}
    )
    out = os.path.join(tempfile.mkdtemp(), "composite.tif")
    urllib.request.urlretrieve(url, out)
    return out


def ndvi_from_nir(nir_path):
    """Compute NDVI = (NIR - Red)/(NIR + Red) from a B8,B4 raster; return (ndvi_array, source)."""
    src = Dataset.read_file(nir_path)
    a = np.asarray(src.read_array(), dtype="float32")
    ndvi = (a[0] - a[1]) / (a[0] + a[1] + 1e-6)
    return ndvi, src

EarthEngineAuth.initialize(SERVICE_ACCOUNT, SERVICE_KEY) # one-time ee.Initialize GEOM = ee.Geometry.Rectangle(AOI) # [W, S, E, N] def mask_s2_clouds(img): """Drop cloud / shadow / cirrus pixels via the Sentinel-2 SCL band. SCL classes removed: 3 (cloud shadow), 8 (cloud, medium prob), 9 (cloud, high prob), 10 (thin cirrus). """ scl = img.select("SCL") clear = scl.neq(3).And(scl.neq(8)).And(scl.neq(9)).And(scl.neq(10)) return img.updateMask(clear) def composite_tif(start, end, bands, scale=10): """Cloud-masked median S2 composite over [start, end) for the AOI; return the GeoTIFF path.""" collection = ( ee.ImageCollection(COLL) .filterBounds(GEOM) .filterDate(start, end) .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 60)) .map(mask_s2_clouds) .select(bands) ) composite = collection.median().clip(GEOM) url = composite.getDownloadURL( {"region": GEOM, "scale": scale, "format": "GEO_TIFF"} ) out = os.path.join(tempfile.mkdtemp(), "composite.tif") urllib.request.urlretrieve(url, out) return out def ndvi_from_nir(nir_path): """Compute NDVI = (NIR - Red)/(NIR + Red) from a B8,B4 raster; return (ndvi_array, source).""" src = Dataset.read_file(nir_path) a = np.asarray(src.read_array(), dtype="float32") ndvi = (a[0] - a[1]) / (a[0] + a[1] + 1e-6) return ndvi, src

Build every window¶

For each window we composite the true-colour bands (B4,B3,B2) and the NIR/red pair (B8,B4), copy the true-colour composite into rgb_dir, derive the NDVI raster, and persist it into ndvi_dir.

In [ ]:

dates = []
for label, start, end in WINDOWS:
    rgb = composite_tif(start, end, ["B4", "B3", "B2"])  # true-colour composite
    nir = composite_tif(start, end, ["B8", "B4"])  # NIR/red pair for NDVI
    tag = label.replace("-", ".")
    shutil.copy(rgb, os.path.join(rgb_dir, f"rgb_{tag}.tif"))
    ndvi, src = ndvi_from_nir(nir)
    out = Dataset.create_from_array(arr=ndvi[None, :, :], geo=src.geotransform, epsg=src.epsg)
    out.to_file(os.path.join(ndvi_dir, f"ndvi_{tag}.tif"))  # persist derived NDVI
    dates.append(label)

print(f"{len(dates)} composites:", dates)

dates = [] for label, start, end in WINDOWS: rgb = composite_tif(start, end, ["B4", "B3", "B2"]) # true-colour composite nir = composite_tif(start, end, ["B8", "B4"]) # NIR/red pair for NDVI tag = label.replace("-", ".") shutil.copy(rgb, os.path.join(rgb_dir, f"rgb_{tag}.tif")) ndvi, src = ndvi_from_nir(nir) out = Dataset.create_from_array(arr=ndvi[None, :, :], geo=src.geotransform, epsg=src.epsg) out.to_file(os.path.join(ndvi_dir, f"ndvi_{tag}.tif")) # persist derived NDVI dates.append(label) print(f"{len(dates)} composites:", dates)

Read each folder as a collection¶

DatasetCollection.read_multiple_files() reads a whole directory in one call — no manual file lists. We build one collection for the true-colour scenes and one for the derived NDVI rasters.

In [ ]:

dc_rgb = DatasetCollection.read_multiple_files(rgb_dir, date=False) if dates else None
dc_ndvi = DatasetCollection.read_multiple_files(ndvi_dir, date=False) if dates else None

dc_rgb = DatasetCollection.read_multiple_files(rgb_dir, date=False) if dates else None dc_ndvi = DatasetCollection.read_multiple_files(ndvi_dir, date=False) if dates else None

2 · True-colour progression — Dataset.plot(rgb_options=...)¶

Each date is rendered as a real-colour RGB composite (B4,B3,B2), 2–98th-percentile stretched.

In [ ]:

if dc_rgb:
    n = dc_rgb.time_length
    fig, ax = plt.subplots(1, n, figsize=(2.4 * n, 2.7))
    for i, d in enumerate(dates):
        dc_rgb.iloc(i).plot(
            rgb_options={"rgb": [0, 1, 2], "percentile": 2}, ax=ax[i], title=d
        )
        ax[i].set_xticks([])
        ax[i].set_yticks([])
    fig.suptitle("Khan Younis - 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")

if dc_rgb: n = dc_rgb.time_length fig, ax = plt.subplots(1, n, figsize=(2.4 * n, 2.7)) for i, d in enumerate(dates): dc_rgb.iloc(i).plot( rgb_options={"rgb": [0, 1, 2], "percentile": 2}, ax=ax[i], title=d ) ax[i].set_xticks([]) ax[i].set_yticks([]) fig.suptitle("Khan Younis - 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 · Before / after — true colour¶

The first and last scenes side by side at full size.

In [ ]:

if dc_rgb:
    fig, ax = plt.subplots(1, 2, figsize=(9, 4.5))
    dc_rgb.iloc(0).plot(
        rgb_options={"rgb": [0, 1, 2], "percentile": 2},
        ax=ax[0],
        title=f"{dates[0]} (before)",
    )
    dc_rgb.iloc(dc_rgb.time_length - 1).plot(
        rgb_options={"rgb": [0, 1, 2], "percentile": 2},
        ax=ax[1],
        title=f"{dates[-1]} (after)",
    )
    for a in ax:
        a.set_xticks([])
        a.set_yticks([])
    plt.tight_layout()
    plt.show()

if dc_rgb: fig, ax = plt.subplots(1, 2, figsize=(9, 4.5)) dc_rgb.iloc(0).plot( rgb_options={"rgb": [0, 1, 2], "percentile": 2}, ax=ax[0], title=f"{dates[0]} (before)", ) dc_rgb.iloc(dc_rgb.time_length - 1).plot( rgb_options={"rgb": [0, 1, 2], "percentile": 2}, ax=ax[1], title=f"{dates[-1]} (after)", ) for a in ax: a.set_xticks([]) a.set_yticks([]) plt.tight_layout() plt.show()

4 · Animated true colour — DatasetCollection.plot(rgb_options=...).save_animation()¶

The seven scenes as a real-colour animation: DatasetCollection.plot(rgb_options=...) composites the RGB bands per frame (the same percentile stretch Dataset.plot applies) and animates them across all dates.

In [ ]:

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="Khan Younis - 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")

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="Khan Younis - 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")

5 · NDVI before / after / change — Dataset.plot()¶

NDVI is a single-band index, so it is rendered with a colour map rather than as an RGB composite. We pull the May 2023 and May 2024 NDVI arrays and difference them.

In [ ]:

if dc_ndvi:
    i_b, i_a = dates.index("2023-05"), dates.index("2024-05")
    nd_b = dc_ndvi.iloc(i_b).read_array(band=0)
    nd_a = dc_ndvi.iloc(i_a).read_array(band=0)

if dc_ndvi: i_b, i_a = dates.index("2023-05"), dates.index("2024-05") nd_b = dc_ndvi.iloc(i_b).read_array(band=0) nd_a = dc_ndvi.iloc(i_a).read_array(band=0)

Render the before, after, and change maps side by side: the two dates with a green colour map and the 2024−2023 difference with a diverging one.

In [ ]:

if dc_ndvi:
    fig, ax = plt.subplots(1, 3, figsize=(12, 4))
    dc_ndvi.iloc(i_b).plot(
        band=0, ax=ax[0], cmap="RdYlGn", title="NDVI - May 2023", vmin=-0.2, vmax=0.6
    )
    dc_ndvi.iloc(i_a).plot(
        band=0, ax=ax[1], cmap="RdYlGn", title="NDVI - May 2024", vmin=-0.2, vmax=0.6
    )
    diff = Dataset.create_from_array(
        arr=(nd_a - nd_b)[None, :, :],
        geo=dc_ndvi.iloc(0).geotransform,
        epsg=dc_ndvi.iloc(0).epsg,
    )
    diff.plot(
        band=0,
        ax=ax[2],
        cmap="RdBu",
        title="NDVI change (2024 - 2023)",
        vmin=-0.5,
        vmax=0.5,
    )
    plt.tight_layout()
    plt.show()

if dc_ndvi: fig, ax = plt.subplots(1, 3, figsize=(12, 4)) dc_ndvi.iloc(i_b).plot( band=0, ax=ax[0], cmap="RdYlGn", title="NDVI - May 2023", vmin=-0.2, vmax=0.6 ) dc_ndvi.iloc(i_a).plot( band=0, ax=ax[1], cmap="RdYlGn", title="NDVI - May 2024", vmin=-0.2, vmax=0.6 ) diff = Dataset.create_from_array( arr=(nd_a - nd_b)[None, :, :], geo=dc_ndvi.iloc(0).geotransform, epsg=dc_ndvi.iloc(0).epsg, ) diff.plot( band=0, ax=ax[2], cmap="RdBu", title="NDVI change (2024 - 2023)", vmin=-0.5, vmax=0.5, ) plt.tight_layout() plt.show()

Summarise the shift: mean NDVI and the vegetated fraction (NDVI > 0.3) before vs after.

In [ ]:

if dc_ndvi:
    veg = lambda x: float(np.mean(x > 0.3))
    print(
        f"mean NDVI 2023 {np.nanmean(nd_b):.3f} -> 2024 {np.nanmean(nd_a):.3f}; "
        f"vegetated {veg(nd_b):.1%} -> {veg(nd_a):.1%}"
    )

if dc_ndvi: veg = lambda x: float(np.mean(x > 0.3)) print( f"mean NDVI 2023 {np.nanmean(nd_b):.3f} -> 2024 {np.nanmean(nd_a):.3f}; " f"vegetated {veg(nd_b):.1%} -> {veg(nd_a):.1%}" )

6 · Animated NDVI — DatasetCollection.plot().save_animation()¶

The collection animates one band across all seven dates; NDVI is the natural single-band quantity to watch. The plot() call and save_animation() are split onto their own lines.

In [ ]:

if dc_ndvi:
    gif = os.path.join(tempfile.mkdtemp(), "ndvi.gif")
    glyph = dc_ndvi.plot(
        band=0,
        cmap="RdYlGn",
        figsize=(5, 5),
        title="Khan Younis - NDVI 2023 -> 2024",
    )
    glyph.save_animation(gif, fps=1.2)
    plt.close("all")  # drop the static builder frame; keep only the GIF
    display(Image(filename=gif))

if dc_ndvi: gif = os.path.join(tempfile.mkdtemp(), "ndvi.gif") glyph = dc_ndvi.plot( band=0, cmap="RdYlGn", figsize=(5, 5), title="Khan Younis - NDVI 2023 -> 2024", ) glyph.save_animation(gif, fps=1.2) plt.close("all") # drop the static builder frame; keep only the GIF display(Image(filename=gif))

7 · NDVI timeline across all seven dates¶

Mean NDVI per date as a line chart, with a dashed reference line at the 0.3 vegetated threshold.

In [ ]:

if dc_ndvi:
    series_v = [
        float(np.nanmean(dc_ndvi.iloc(i).read_array(band=0)))
        for i in range(dc_ndvi.time_length)
    ]
    fig, ax = plt.subplots(figsize=(9, 4))
    ax.plot([pd.Timestamp(d) for d in dates], series_v, "-o", color="tab:green")
    ax.axhline(0.3, color="0.7", ls="--", lw=0.8)
    ax.set(ylabel="mean NDVI", title="Khan Younis - mean vegetation index over time")
    fig.autofmt_xdate()
    plt.show()

if dc_ndvi: series_v = [ float(np.nanmean(dc_ndvi.iloc(i).read_array(band=0))) for i in range(dc_ndvi.time_length) ] fig, ax = plt.subplots(figsize=(9, 4)) ax.plot([pd.Timestamp(d) for d in dates], series_v, "-o", color="tab:green") ax.axhline(0.3, color="0.7", ls="--", lw=0.8) ax.set(ylabel="mean NDVI", title="Khan Younis - mean vegetation index over time") fig.autofmt_xdate() plt.show()

Recap¶

Each frame is a cloud-masked median composite over a ~6-week window — the SCL band masks clouds per pixel and median() fills from the clear observations, so no hand-picked clear dates are needed. Both series are read from disk with DatasetCollection.read_multiple_files(). The progression and before/after panels render in true colour with Dataset.plot(rgb_options=...), the same scenes animate as a real-colour GIF (frames prepared by cleopatra's prepare_array), and the NDVI maps/animation come from Dataset.plot() / DatasetCollection.plot().save_animation().

Try it yourself¶

• Re-point AOI, add windows, tighten the CLOUDY_PIXEL_PERCENTAGE / SCL mask, or animate a different derived index.
• See the GEE backend reference.