Inspecting Time Series at a Clicked Pixel

On click, find the scenes covering the point, read a single pixel from each, and plot the valid values against date:

def on_click(**kw):
    if kw.get("type") == "click":
        lat, lon = kw["coordinates"]
        series = read_point_series(lon, lat)       # one tiny read per scene
        plot_series(series)

m.on_interaction(on_click)

Clicking through a handful of locations finds seasonal oddities, bad dates and mis-scaled scenes faster than any aggregate. This page belongs to interactive raster exploration in Jupyter in Visualization, Tiling & Web Delivery.

Why a Point Read Is Cheap

A single-pixel read from a Cloud-Optimized GeoTIFF costs one header fetch and one internal tile — tens of kilobytes regardless of the scene’s size. That is what makes a whole time series interactive.

Fifty scenes, fifty tiny reads A clicked location is transformed into each scene's CRS and pixel indices. For every scene covering it, one range request fetches the header if not cached and a second fetches the single internal tile containing that pixel. Fifty scenes cost about two megabytes in total, compared with tens of gigabytes to download them. One point through a year of scenes click 34.75 E, 0.35 N scene 1 · 40 KB scene 2 · 40 KB scene 3 · 40 KB … 47 more total: about 2 MB versus about 40 GB to download the same fifty scenes in full Threads overlap the latency, so fifty reads finish in under a second. The same property is what makes chip extraction from remote archives practical.

Environment & Setup

Package Version pin Used for
pystac-client >=0.7 Finding every scene that covers the point
rasterio >=1.3.0 Single-window remote reads
matplotlib >=3.8 Plotting the series
leafmap >=0.32 The map and its click events 
pip install "pystac-client>=0.7" "rasterio>=1.3.0" "matplotlib>=3.8" "leafmap>=0.32"

Complete Working Example

from concurrent.futures import ThreadPoolExecutor

import matplotlib.pyplot as plt
import numpy as np
import pystac_client
import rasterio
import rasterio.warp
from rasterio.windows import Window

API = "https://earth-search.aws.element84.com/v1"
SCL_CLEAR = {4, 5, 6, 11}          # vegetation, bare, water, snow


def _read(href: str, lon: float, lat: float, size: int = 3) -> float:
    with rasterio.open(href) as src:
        xs, ys = rasterio.warp.transform("EPSG:4326", src.crs, [lon], [lat])
        row, col = src.index(xs[0], ys[0])
        win = Window(col - size // 2, row - size // 2, size, size)
        arr = src.read(1, window=win, boundless=True, masked=True).astype("float32")
    return float(arr.mean()) if arr.count() else np.nan


def point_series(lon: float, lat: float, start: str, end: str) -> dict[str, tuple]:
    items = list(pystac_client.Client.open(API).search(
        collections=["sentinel-2-l2a"],
        intersects={"type": "Point", "coordinates": [lon, lat]},
        datetime=f"{start}/{end}",
    ).items())

    def one(item):
        red = _read(item.assets["red"].href, lon, lat)
        nir = _read(item.assets["nir"].href, lon, lat)
        scl = _read(item.assets["scl"].href, lon, lat, size=1)
        ndvi = (nir - red) / (nir + red) if (nir + red) else np.nan
        return item.datetime.date().isoformat(), (ndvi, int(scl) in SCL_CLEAR)

    with ThreadPoolExecutor(max_workers=8) as pool:
        return dict(sorted(pool.map(one, items)))


def plot_series(series: dict[str, tuple], title: str = "") -> None:
    dates = np.array(list(series), dtype="datetime64[D]")
    vals = np.array([v for v, _ in series.values()])
    clear = np.array([c for _, c in series.values()])

    fig, ax = plt.subplots(figsize=(8, 3))
    ax.plot(dates[clear], vals[clear], "-o", color="#15803d", label="clear")
    ax.plot(dates[~clear], vals[~clear], "o", mfc="none", color="#64748b",
            label="masked (cloud, shadow)")
    ax.set_ylabel("NDVI")
    ax.set_ylim(-0.2, 1.0)
    ax.legend(loc="lower right", fontsize=8)
    ax.set_title(title, fontsize=9, loc="left")
    fig.autofmt_xdate()

Two design decisions carry the method. Reading a 3×3 neighbourhood and averaging it smooths the pixel-level registration jitter between acquisitions, which otherwise makes a stable field look noisy. And masked observations are plotted rather than dropped — hollow markers show where the gaps are, which is diagnostic information in its own right. The masking itself follows masking clouds with the Sentinel-2 SCL band.

Reading What the Series Shows

Four things a clicked series reveals A seasonal NDVI curve rises and falls normally. One clear-flagged observation dips sharply, which is residual cloud the mask missed. Several masked observations cluster around the peak, so the peak itself is unobserved. Late in the series the value drops abruptly and stays low, which is genuine change such as a harvest or clearing. What to look for in one pixel's year missed cloud peak unobserved abrupt, sustained drop real change normal green-up

The isolated dip is the commonest finding and the most useful: a single clear-flagged observation far below its neighbours is almost always residual cloud or shadow that the mask missed, and seeing it tells you the mask needs dilation, as discussed in dilating cloud masks to catch thin cirrus.

A cluster of hollow markers around the peak matters for anything downstream that summarises the season. Percentile features computed from that pixel will describe the shoulders of the curve, not its top, which is exactly the issue that motivates carrying an observation count as a feature in building temporal feature stacks from image time series.

The abrupt, sustained drop is the signal change detection looks for. Seeing it at one pixel is anecdote; seeing it at a dozen clicked pixels across a field is a reason to trust — or question — whatever the change map says about that field.

Verification

Confirm the click lands on the right pixel The value returned by the click handler for one date is compared with an independent read of the same scene at the same coordinate. Agreement confirms the coordinate transform and the index lookup; disagreement means every point in the series is being read from the wrong place. One independent read validates the whole chain click handler 2026-06-14: NDVI 0.712 independent read 2026-06-14: NDVI 0.712 match Do this once per session; a mismatch invalidates every series you have looked at. 
import rasterio
import rasterio.warp

lon, lat = 34.7512, 0.3504
series = point_series(lon, lat, "2026-06-01", "2026-06-30")
date, (ndvi, _) = next(iter(series.items()))

# Independent path: open one scene directly and compute at the same place
with rasterio.open(red_href) as r, rasterio.open(nir_href) as n:
    xs, ys = rasterio.warp.transform("EPSG:4326", r.crs, [lon], [lat])
    row, col = r.index(xs[0], ys[0])
    rv = r.read(1, window=((row - 1, row + 2), (col - 1, col + 2))).mean()
    nv = n.read(1, window=((row - 1, row + 2), (col - 1, col + 2))).mean()
assert abs((nv - rv) / (nv + rv) - ndvi) < 1e-4

The comparison is cheap and it validates the whole chain from click to value. A mismatch almost always traces to a skipped coordinate transform or to reading the red asset from one item and the near-infrared from another.

Common Errors

Every value is NaN

The point falls outside the scenes returned, usually because the search used a bbox around the wrong coordinate order. Use intersects with a GeoJSON point in longitude, latitude order.

The series is extremely noisy

Single-pixel reads are picking up registration jitter. Read and average a 3×3 neighbourhood.

Values are in the thousands

The assets are scaled integers. Divide by the scale factor before computing the index, as described in handling nodata and scale factors in band math.

The plot takes a long time to appear

Reads are serial. Use a small thread pool; latency, not bandwidth, is the limit. Eight workers is plenty — beyond that the catalog search, not the reads, dominates.

Frequently Asked Questions

Q: Is reading one pixel from many remote files slow? Each read fetches the file header and one internal tile — a few tens of kilobytes. Fifty scenes take a few seconds serially and well under a second with a small thread pool, which is fast enough for interactive use.

Q: Should I plot masked observations? Plot them differently rather than hiding them. Showing cloud-masked dates as hollow markers makes the gaps visible, and a series with no observations near its peak is a series whose seasonal statistics cannot be trusted.

Q: Why read a small neighbourhood rather than one pixel? Because a single pixel carries geolocation jitter between acquisitions of up to a pixel or so. Averaging a three by three neighbourhood smooths that jitter and gives a series that reflects the surface rather than the registration error.

Q: Can this run against a local stack instead of a catalog? Yes, and it is simpler: open the stacked file once and read the same window across the time dimension. The catalog version exists because most archives are not stacked, and a point read against remote scenes avoids ever building the stack.