Interactive Raster Exploration in Jupyter

For most remote sensing teams the notebook is where the data is actually looked at. Whether a cloud mask is behaving, whether a composite has a seam, whether a classification is plausible — these are answered by panning around an image and clicking on pixels, not by reading a summary statistic. This topic, part of Visualization, Tiling & Web Delivery, covers doing that properly rather than by trial and error.

One decision shapes everything else: whether the map shows an array you loaded or tiles a server rendered. Loading arrays is simple and bounded by memory, so it works for a chip and fails for a scene. Pointing at a tile endpoint has no size limit, streams only what is visible, and — crucially — shows precisely what everyone else will see.

---

Prerequisites

pip install "leafmap>=0.32" "ipyleaflet>=0.18" "rasterio>=1.3.0" "geopandas>=0.14" "matplotlib>=3.8"

A tile endpoint is assumed for the layer work — either a hosted service or a local one, as set up in serving raster tiles with TiTiler.

---

Step-by-Step Workflow

Step 1 — Put the raster on the map as tiles

import leafmap


def add_cog(m: leafmap.Map, url: str, *, name: str,
            rescale: tuple[float, float] = (-0.2, 0.9),
            colormap: str = "rdylgn", tiler: str = "https://tiles.example.org") -> str:
    """Add a COG to the map through a tile endpoint and return the template."""
    template = (
        f"{tiler}/cog/tiles/WebMercatorQuad///.png"
        f"?url={url}&rescale={rescale[0]},{rescale[1]}&colormap_name={colormap}"
    )
    m.add_tile_layer(template, name=name, attribution="Sentinel-2")
    return template


m = leafmap.Map(center=(0.35, 34.75), zoom=11)
add_cog(m, "s3://example-bucket/products/ndvi_20260614.tif", name="NDVI 14 June")
m

Returning the template rather than only mutating the map is a small habit worth keeping: the string is the shareable artifact, and having it in a variable means it can be printed, logged into a notebook cell, or handed to a colleague verbatim.

Step 2 — Draw an area of interest and get it back in Python

import geopandas as gpd
import rasterio
from rasterio.mask import mask as rio_mask


def read_drawn_aoi(m: leafmap.Map, raster_path: str):
    """Read the pixels under the last shape drawn on the map."""
    if not m.draw_features:
        raise ValueError("draw a polygon on the map first")

    aoi = gpd.GeoDataFrame.from_features(m.draw_features, crs="EPSG:4326")
    with rasterio.open(raster_path) as src:
        aoi_proj = aoi.to_crs(src.crs)            # always, before anything else
        arr, transform = rio_mask(src, aoi_proj.geometry, crop=True, filled=False)
    return arr, transform, aoi_proj

rio_mask with crop=True does the window arithmetic and the geometry masking in one call, reading only the pixels inside the shape’s bounding box. For irregular shapes the details and the edge cases are in how to clip rasters to irregular polygon boundaries.

Step 3 — Read pixel values on click

Clicking a pixel and seeing its values across bands is the fastest debugging tool in remote sensing. It finds unscaled bands, shifted masks and missing dates in seconds.

import rasterio
from ipyleaflet import Popup
from ipywidgets import HTML


def attach_pixel_probe(m: leafmap.Map, raster_path: str, band_names: list[str]) -> None:
    src = rasterio.open(raster_path)

    def on_click(**kwargs):
        if kwargs.get("type") != "click":
            return
        lat, lon = kwargs["coordinates"]
        # Transform the clicked degrees into the raster's CRS, then to indices
        xs, ys = rasterio.warp.transform("EPSG:4326", src.crs, [lon], [lat])
        row, col = src.index(xs[0], ys[0])
        if not (0 <= row < src.height and 0 <= col < src.width):
            return
        win = rasterio.windows.Window(col, row, 1, 1)
        values = src.read(window=win).ravel()

        rows = "".join(f"<tr><td>{n}</td><td>{v}</td></tr>"
                       for n, v in zip(band_names, values))
        popup = Popup(location=(lat, lon),
                      child=HTML(f"<table>{rows}</table>"), close_button=True)
        m.add_layer(popup)

    m.on_interaction(on_click)

Holding the dataset open across clicks is deliberate: reopening per click adds a round trip to object storage and makes the interaction feel sluggish. Close it when the notebook is done with it, or use a context manager around the whole exploration session.

Step 4 — Compare two dates side by side

import leafmap

before = add_cog_template("s3://example-bucket/products/ndvi_20250614.tif",
                          rescale=(-0.2, 0.9))
after = add_cog_template("s3://example-bucket/products/ndvi_20260614.tif",
                         rescale=(-0.2, 0.9))       # the SAME stretch, deliberately

m = leafmap.Map(center=(0.35, 34.75), zoom=12)
m.split_map(left_layer=before, right_layer=after)
m

Sharing the stretch between the two panes is the entire point. Two panes rendered on their own percentiles will differ in brightness for reasons that have nothing to do with the ground, and a viewer will read that difference as change. The arithmetic behind a real change product is in change detection and differencing workflows.

---

Probing a Time Series at a Point

import matplotlib.pyplot as plt
import rasterio
import rasterio.warp


def probe_series(paths: dict[str, str], lon: float, lat: float):
    """Value at one location across a dated stack of single-band rasters."""
    dates, values = [], []
    for date, path in sorted(paths.items()):
        with rasterio.open(path) as src:
            xs, ys = rasterio.warp.transform("EPSG:4326", src.crs, [lon], [lat])
            row, col = src.index(xs[0], ys[0])
            if not (0 <= row < src.height and 0 <= col < src.width):
                continue
            win = rasterio.windows.Window(col, row, 1, 1)
            arr = src.read(1, window=win, masked=True)
        dates.append(date)
        values.append(float(arr.filled(float("nan")).ravel()[0]))

    fig, ax = plt.subplots(figsize=(7, 3))
    ax.plot(dates, values, marker="o")
    ax.set_ylabel("NDVI")
    ax.tick_params(axis="x", rotation=45)
    return fig, dict(zip(dates, values))

This is a single-pixel read per date, so it is fast even against remote files — a few kilobytes each. The gaps matter as much as the values: a location whose series is missing every date around the seasonal peak will have unreliable percentile features, which is exactly the diagnosis that leads back to building temporal feature stacks from image time series.

---

Keeping Notebook Maps Honest

Three habits keep interactive exploration from misleading the person doing it.

The first is to view at native resolution before drawing conclusions about detail. A map at zoom 10 is showing an overview level, so single-pixel artefacts are invisible and small features are smoothed away. Zoom to the level where one screen pixel is one data pixel before judging whether a mask is tight or a boundary is clean.

The second is to be explicit about which layer is on top. Stacking a classification over a true-colour composite is the right way to judge plausibility, but an opaque top layer hides the evidence; set the opacity to about 0.6 so both are legible, and toggle rather than guess.

The third is to distrust anything that looks perfect at low zoom. Overviews are resampled, and for a class raster resampled with the wrong method they can show a composition the data does not have — the failure described in writing prediction rasters as COGs. If a map looks cleaner zoomed out than zoomed in, suspect the pyramid before congratulating the model.

---

Layer Stacks That Answer a Question

An exploration session usually needs more than one layer, and the order they are stacked in decides what can be judged.

The stack above answers one specific question — does this mask follow what the imagery shows — and it is worth assembling deliberately rather than adding layers as they occur to you. A dark basemap under a dark classification, or an opaque overlay above the only evidence, produces a map that cannot answer anything.

Toggling is underrated as a technique. Flicking a mask layer on and off repeatedly makes a misalignment of even one or two pixels obvious, where a static side-by-side comparison hides it entirely. The same trick applied to two dates is what makes the split view in step four so much more informative than two separate figures.

For a classification, add a legend as a small HTML widget built from the same class table the raster’s colour map came from, so the map and the file agree. Deriving it from the file rather than typing it into the notebook removes the possibility that they disagree, which they eventually will.

---

Parameter Reference

---

From Exploration to a Repeatable Check

The habit worth building is turning a useful exploration into something that runs without a human. Most interactive sessions end with a judgement — “the mask is a pixel too tight along field edges”, “these three scenes have the wrong scale factor” — and that judgement is usually expressible as an assertion.

A session that found unscaled bands becomes a check on the value range of every incoming scene. A session that found a shifted mask becomes an assertion that the mask and the imagery share a transform. A session that found a bad date becomes a quicklook contact sheet reviewed after every ingest. None of those take long to write, and each one means the same defect never needs to be found by hand again.

The notebook is then doing what it is good at — forming a hypothesis quickly from something visible — while the pipeline does what it is good at, which is checking that hypothesis on everything, forever. Keeping the two separate also keeps the notebook honest: a notebook that quietly becomes production code acquires state nobody can reproduce, whereas one that ends in a written assertion leaves the assertion behind and can be thrown away.

---

Verification & Testing

Interactive work resists automated testing, but the pieces underneath it do not.

import geopandas as gpd
import rasterio

aoi = gpd.GeoDataFrame.from_features(m.draw_features, crs="EPSG:4326")
assert not aoi.empty, "nothing drawn"

with rasterio.open(raster_path) as src:
    aoi_proj = aoi.to_crs(src.crs)
    minx, miny, maxx, maxy = aoi_proj.total_bounds
    b = src.bounds
    assert minx < b.right and maxx > b.left, "AOI does not overlap the raster in x"
    assert miny < b.top and maxy > b.bottom, "AOI does not overlap the raster in y"

The overlap assertion is worth having because the failure it catches — a drawn shape that misses the scene entirely, usually after a reprojection was skipped — otherwise surfaces as an empty array that looks like a data problem rather than a coordinate problem.

For the probe, check one known location against an independent read. Open the file separately, read the same window, and confirm the value matches what the popup showed; a mismatch means the click-to-index transform is wrong, which affects every subsequent inspection.

---

Troubleshooting

The map pane is blank

Either the tile endpoint is unreachable from the notebook, or the layer is outside the current view. Open the tile template in a browser tab: if it returns a PNG, the problem is the map; if it errors, the problem is the endpoint.

The drawn polygon reads the wrong pixels

The geometry was used in degrees against a projected raster. Reproject with to_crs(src.crs) before any read, every time.

Clicking returns values that look like nodata everywhere

The dataset was opened without masked=True, so the fill value is being reported as a real number, or the click transform is landing outside the scene. Check both with a location you know the answer for.

The notebook becomes unresponsive after adding a layer

An array was loaded into the browser rather than a tile layer. Switch to a tile endpoint; a full scene as an in-browser image is tens of megabytes of base64.

Exported HTML shows an empty map

The tile endpoint is private or local. Point the layer at a reachable URL before exporting, or the recipient sees a working map with no data in it.

---

Frequently Asked Questions

Q: Should I load the array into the notebook or point at a tile service? Point at a tile service for anything larger than a small window. A tile layer streams only what is on screen and shows exactly what a stakeholder will see, while loading an array is bounded by memory and shows a decimated version that can hide the defect you are looking for.

Q: Why does a drawn polygon come back in the wrong coordinates? Map widgets return geometry in WGS84 degrees regardless of the layer’s projection. Reproject it into the raster CRS before using it for a window or a mask, or every read will land somewhere else entirely.

Q: Is a notebook map safe to share with stakeholders? Only if the tile endpoint it points at is reachable by them. A notebook exported to HTML keeps the map and the layer URL, so anyone with access to the endpoint sees the live layer, and anyone without it sees an empty pane.

Q: Can I use the same code outside Jupyter? The reading and reprojection helpers are ordinary Python and run anywhere. Only the widget layer is notebook-specific, and it has direct equivalents in a small web page built on the same tile template — which is often the better destination once an exploration becomes something other people need repeatedly.

---

Related

• Displaying COGs on a leafmap Web Map — layer setup, basemaps and opacity in detail.
• Drawing an AOI and Reading Pixels Interactively — the draw control, reprojection and masked reads.
• Comparing Two Dates with a Split Map — shared stretches and swipe comparison.
• Inspecting Time Series at a Clicked Pixel — probing a dated stack and reading the gaps.
• Serving Raster Tiles with TiTiler — the endpoint every layer here consumes.