Drawing an AOI and Reading Pixels Interactively

Capture the drawn shape, reproject it into the raster CRS, and read only the pixels inside it:

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

aoi = gpd.GeoDataFrame.from_features(m.draw_features, crs="EPSG:4326")
with rasterio.open("ndvi.tif") as src:
    arr, transform = mask(src, aoi.to_crs(src.crs).geometry, crop=True, filled=False)
print(float(arr.mean()), int(arr.count()))

The to_crs call is the one that cannot be skipped. This page belongs to interactive raster exploration in Jupyter in Visualization, Tiling & Web Delivery.

What Happens Between the Click and the Read

Bounding box read, polygon mask The drawn polygon is reprojected into the raster CRS. Its bounding box defines the window that is actually read from the file, so only those bytes are fetched. Pixels inside the window but outside the polygon are then masked, so the summary statistics describe only the drawn area. Read the box, keep the shape kept masked never read Dashed: the window read from disk. Solid: the pixels that count toward the statistics.

rasterio.mask.mask with crop=True computes the window from the geometry’s bounds, reads only that window, and masks everything outside the polygon. With filled=False the result is a masked array, so outside pixels and nodata pixels are both excluded from mean() and count() automatically.

Environment & Setup

Package Version pin Used for
leafmap >=0.32 Map widget with a draw control
geopandas >=0.14 Holding and reprojecting the drawn geometry
rasterio >=1.3.0 Windowed, masked reads
numpy >=1.23 Summary statistics over the masked array 
pip install "leafmap>=0.32" "geopandas>=0.14" "rasterio>=1.3.0" "numpy>=1.23"

Complete Working Example

import geopandas as gpd
import leafmap
import numpy as np
import rasterio
from ipywidgets import Output
from rasterio.mask import mask

RASTER = "s3://example-bucket/products/ndvi_20260614.tif"
out = Output()

m = leafmap.Map(center=(0.35, 34.75), zoom=12, draw_control=True)
m.add_cog_layer(RASTER.replace("s3://", "https://example-bucket.s3.amazonaws.com/")
                .replace("example-bucket/", "", 1),
                name="NDVI", rescale="-0.2,0.9", colormap_name="rdylgn")


def summarise(geojson_features: list[dict], raster_path: str) -> dict:
    aoi = gpd.GeoDataFrame.from_features(geojson_features, crs="EPSG:4326")
    with rasterio.open(raster_path) as src:
        aoi_proj = aoi.to_crs(src.crs)                      # degrees -> raster CRS
        arr, transform = mask(src, aoi_proj.geometry, crop=True,
                              filled=False, all_touched=False)
        px_area = abs(src.transform.a * src.transform.e)
    band = arr[0]
    valid = band.compressed()
    return {
        "pixels": int(valid.size),
        "area_ha": round(valid.size * px_area / 10_000, 2),
        "mean": round(float(valid.mean()), 3) if valid.size else None,
        "p10": round(float(np.percentile(valid, 10)), 3) if valid.size else None,
        "p90": round(float(np.percentile(valid, 90)), 3) if valid.size else None,
        "nodata_fraction": round(float(np.ma.getmaskarray(band).mean()), 3),
    }


def on_draw(target, action, geo_json):
    if action != "created":
        return
    with out:
        out.clear_output()
        print(summarise([geo_json], RASTER))


m.draw_control.on_draw(on_draw)
display(m, out)

Reporting nodata_fraction alongside the statistics is the detail that keeps a summary honest. A mean computed over an AOI that is 60% cloud-masked describes the 40% that remained, and without the fraction nobody reading the number would know.

Choosing Between Exact and Touched Pixels

Small AOIs are sensitive to the edge rule For a large AOI the choice between counting pixels whose centre is inside and counting every pixel the boundary touches changes the pixel count by a few per cent. For a small AOI a few pixels across, the same choice can change the count by half and move the mean noticeably, because edge pixels are a large share of the total. Edge pixels as a share of the AOI 5 px wide: 64% edge 100 px wide: 4% edge tiny AOI large AOI Below about 20 pixels across, state which rule you used — the answer depends on it.

The default pixel-centre rule is right for area and statistics because it treats each pixel as belonging to exactly one side of the boundary. all_touched=True is useful when the point is to see every pixel that could be relevant — inspecting a thin feature such as a river bank, say — but it inflates the count and biases the mean toward whatever surrounds the shape. The same trade-off appears in computing zonal statistics with rasterstats, where it matters for the same reasons.

Beyond a Single Band

The same pattern extends naturally to a multi-band stack or a time series, and both are more useful than a single mean.

For a stack, read every band inside the shape and report a small spectral profile — the mean reflectance per band — which is often the fastest way to tell whether an unusual patch is water, bare soil, burned ground or cloud shadow. A profile that is flat and low across the near and shortwave infrared is water; one that rises steeply into the shortwave is bare soil; one with a deep near-infrared trough after a fire is a burn scar. Reading those signatures off a drawn polygon takes seconds and settles questions that a single index cannot.

import numpy as np
import rasterio
from rasterio.mask import mask


def spectral_profile(stack_path: str, aoi_proj) -> dict[str, float]:
    with rasterio.open(stack_path) as src:
        arr, _ = mask(src, aoi_proj.geometry, crop=True, filled=False)
        names = list(src.descriptions)
    return {n: round(float(arr[i].mean()), 4) for i, n in enumerate(names)}

For a time series, apply the same mask to each date and collect the mean, which gives a small per-polygon trajectory rather than a per-pixel one. Averaging over the polygon suppresses pixel noise, so the seasonal curve of a field is far cleaner than any single pixel’s — at the cost of mixing whatever else the polygon contains, which is why the edge-pixel rule above matters more for small shapes. The per-pixel alternative is covered in inspecting time series at a clicked pixel.

Verification

Check overlap before reading An AOI reprojected correctly falls inside the raster bounds and the read proceeds. An AOI used in degrees against a projected raster lands near the coordinate origin, far outside the bounds, and the overlap check fails with a clear message instead of returning an empty array that looks like missing data. Fail on the coordinate error, not on the empty result reprojected AOI inside the bounds degrees read as metres: near the origin 
import geopandas as gpd
import rasterio

aoi = gpd.GeoDataFrame.from_features(m.draw_features, crs="EPSG:4326")
with rasterio.open(RASTER) as src:
    minx, miny, maxx, maxy = aoi.to_crs(src.crs).total_bounds
    b = src.bounds
    overlaps = minx < b.right and maxx > b.left and miny < b.top and maxy > b.bottom
    assert overlaps, "AOI does not overlap the raster — was it reprojected?"

Failing on the overlap check turns a confusing empty result into a clear message about coordinates. Without it, a skipped reprojection produces either a ValueError from the mask function or an all-masked array, both of which look like a data problem rather than the geometry problem they are.

A second check worth doing by eye: draw a small AOI over a feature whose value you know — a lake should have low NDVI, a forest high — and confirm the summary agrees. It takes ten seconds and validates the whole chain from click to statistic.

Common Errors

ValueError: Input shapes do not overlap raster

The AOI was used in degrees against a projected raster. Reproject with to_crs(src.crs) first.

Statistics include obviously wrong values

The raster has no nodata set, so fill values are counted. Pass nodata= to mask, or fix the file’s metadata.

The handler fires but nothing prints

Output from a widget callback goes nowhere unless captured. Route it through an Output widget as the example does.

Reading a large AOI is slow

The read is at full resolution over a large area. Read from an overview with out_shape for exploration, and note that the statistics are approximate.

Frequently Asked Questions

Q: What CRS does a drawn shape come back in? WGS84 longitude and latitude, always, regardless of the basemap or the layers on the map. Reproject it into the raster’s CRS before using it for any read.

Q: How large an area can I read interactively? As much as fits comfortably in memory at full resolution — a few hundred square kilometres at 10 metres for a handful of bands. For larger areas read from an overview with out_shape, and say so in the summary, since statistics from an overview are approximations.

Q: Can the drawn AOI be saved for later? Yes — write the reprojected GeoDataFrame to a GeoPackage. Saving it in the raster’s CRS alongside a WGS84 copy avoids re-deriving it and makes the analysis reproducible without the notebook session.

Q: Can I draw several AOIs and compare them? Yes. Collect each drawn feature with an identifier, run the same summary on each, and assemble the results into a table. Keep the rule for edge pixels identical across them, or small and large AOIs will be compared on different footing.