Extracting Pixel Values at Point Locations

To read raster values at a set of coordinates, transform the points into the raster’s CRS and pass them to dataset.sample:

import geopandas as gpd
import rasterio

points = gpd.read_file("plots.gpkg")

with rasterio.open("ndvi_20230615.tif") as src:
    pts = points.to_crs(src.crs)
    coords = [(geom.x, geom.y) for geom in pts.geometry]
    values = [v[0] for v in src.sample(coords)]      # one value per point, band 1

points["ndvi"] = values

Point sampling is the degenerate case of the workflow in Zonal Statistics and Vector–Raster Integration — one pixel instead of a polygon’s worth — and it has its own set of traps.

Why This Arises in Remote Sensing Workflows

Point sampling is how ground truth meets imagery. Field plots, flux towers, soil samples, validation points from a stratified design, calibration targets — all of them are coordinates that need a value from every scene in a time series. The output feeds model training, accuracy assessment and calibration, which means an error here propagates into a published accuracy figure rather than into a picture.

Three properties make it deceptively easy to get wrong. Sampling never fails loudly: a point outside the raster returns a number, not an exception. There is no averaging to hide a small misalignment, so a half-pixel offset changes the value completely rather than slightly. And the pixel a point falls in is a step function of position, so two points a metre apart across a boundary can differ by more than the signal being studied.

Why point sampling has no tolerance for small errors A survey point sits four metres from a field boundary on a ten metre grid. GPS uncertainty of five metres is enough to move it into the neighbouring pixel, which belongs to a different field with a different NDVI. Unlike a zonal mean, a point sample offers no averaging to absorb the error. A 10 m grid, a 5 m position uncertainty field A — NDVI 0.74 field B — NDVI 0.31 recorded position dashed circle = 5 m GPS uncertainty, which spans the boundary The sampled value is 0.74 or 0.31 — a difference far larger than any effect the study is measuring.

The practical response is not to abandon point sampling but to record enough context to judge it: the pixel actually read, the distance from the point to the nearest pixel boundary, and optionally a small neighbourhood summary that shows whether the surroundings are homogeneous.

Environment & Setup

Package Version Why
rasterio ≥1.3.0 sample, index, xy, windowed reads
geopandas ≥0.13 Point I/O and CRS transformation
pandas ≥2.0 Tidy output
numpy ≥1.23 Neighbourhood statistics 
pip install "rasterio>=1.3.0" "geopandas>=0.13" "pandas>=2.0"

Complete Working Example

This function samples every point in one raster, flags out-of-bounds points and nodata hits explicitly, and reports how far each point sits from its pixel’s edge.

Fields that make a sampled table defensible A sampled value on its own cannot be assessed. Recording whether the point was inside the scene, whether the pixel was masked, how far the point sat from a pixel boundary and the local variability lets a reviewer decide which samples to trust. What to record alongside every sampled value field why it matters typical use in_bounds outside-scene reads return fill exclude, never treat as zero value is NaN masked ≠ missing scene separate cloud from no-coverage edge_distance_m near a boundary, GPS error flips the pixel drop points under half a pixel local std (3×3) high spread = heterogeneous surroundings exclude from calibration sets Four extra columns turn a list of numbers into a table someone can audit. 
import geopandas as gpd
import numpy as np
import pandas as pd
import rasterio


def sample_points(
    raster_path: str,
    points: gpd.GeoDataFrame,
    *,
    id_field: str = "plot_id",
    band: int = 1,
    date: str | None = None,
) -> pd.DataFrame:
    """Sample one raster at every point, keeping the misses visible."""
    with rasterio.open(raster_path) as src:
        pts = points.to_crs(src.crs) if points.crs != src.crs else points
        xs = pts.geometry.x.to_numpy()
        ys = pts.geometry.y.to_numpy()

        left, bottom, right, top = src.bounds
        inside = (xs >= left) & (xs < right) & (ys > bottom) & (ys <= top)

        values = np.full(len(pts), np.nan, dtype="float64")
        if inside.any():
            coords = list(zip(xs[inside], ys[inside]))
            sampled = np.array([v[band - 1] for v in src.sample(coords, indexes=[band])],
                               dtype="float64")
            if src.nodata is not None:
                sampled[sampled == src.nodata] = np.nan   # a fill value is not a measurement
            values[inside] = sampled

        # Distance from each point to the nearest pixel edge, in CRS units
        px, py = abs(src.transform.a), abs(src.transform.e)
        edge_dist = np.minimum(
            np.minimum((xs - left) % px, px - (xs - left) % px),
            np.minimum((top - ys) % py, py - (top - ys) % py),
        )

    return pd.DataFrame({
        id_field: pts[id_field].to_numpy(),
        "date": date,
        "value": values,
        "in_bounds": inside,
        "edge_distance_m": np.round(edge_dist, 2),
    })


if __name__ == "__main__":
    plots = gpd.read_file("plots.gpkg")
    df = sample_points("ndvi_20230615.tif", plots, date="2023-06-15")
    print(df.head())
    print("outside the scene:", int((~df["in_bounds"]).sum()),
          "· nodata hits:", int(df["value"].isna().sum() - (~df["in_bounds"]).sum()))

The in_bounds mask matters more than it looks: without it, points outside the raster receive the nodata value or zero, which is indistinguishable from a genuine measurement of zero. Separating “outside the scene” from “inside but masked” is the difference between a defensible validation table and a misleading one — the same distinction drawn in Extracting nodata and dtype from a GeoTIFF.

Variant Patterns

1. A neighbourhood instead of a single pixel

For validation work, a 3×3 or 5×5 summary around the point is usually more defensible than the single pixel, because it exposes heterogeneity rather than hiding it.

import numpy as np
import rasterio
from rasterio.windows import Window


def sample_window(src, x: float, y: float, size: int = 3, band: int = 1) -> dict:
    """Summarise a size×size neighbourhood centred on a coordinate."""
    row, col = src.index(x, y)
    half = size // 2
    win = Window(col - half, row - half, size, size).intersection(
        Window(0, 0, src.width, src.height))
    arr = src.read(band, window=win, masked=True)
    if arr.count() == 0:
        return {"n": 0}
    return {
        "n": int(arr.count()),
        "centre": float(src.read(band, window=Window(col, row, 1, 1), masked=True).filled(np.nan)[0, 0]),
        "mean": float(arr.mean()),
        "std": float(arr.std()),          # high std = the point sits on a boundary
    }

A large standard deviation across a small neighbourhood is a warning that the point’s value depends on sub-pixel positioning, and it gives a principled reason to exclude that point from a calibration set.

2. Many points across many dates

Where the time goes in a multi-date sampling run Opening a remote raster costs requests and a header parse; sampling costs one small read per point. Opening once per raster and sampling all points inside that context amortises the open across every point, while opening per point multiplies it by the point count. 500 points × 180 dates open per point open open open … 90,000 opens ≈ 4 hours against object storage open per raster open sample all 500 points ≈ 6 minutes, 180 opens GDAL groups the reads for coordinates passed together, so batching also reduces the request count. Sampling is network-bound, so threads across rasters parallelise well; processes mostly add memory. 
from concurrent.futures import ThreadPoolExecutor

import pandas as pd


def sample_series(dated_rasters, points, workers: int = 8) -> pd.DataFrame:
    """Sample the same points across many dates; one open per raster, threads across rasters."""
    def one(item):
        date, path = item
        return sample_points(path, points, date=date)

    with ThreadPoolExecutor(max_workers=workers) as pool:
        frames = list(pool.map(one, dated_rasters))
    return pd.concat(frames, ignore_index=True)

3. Sampling a multi-band stack in one pass

sample accepts indexes and returns a tuple per point, so a whole spectral signature comes back in one read rather than one read per band.

with rasterio.open("stack_20230615.tif") as src:
    bands = [1, 2, 3, 4]                    # e.g. blue, green, red, nir
    rows = list(src.sample(coords, indexes=bands))

signatures = pd.DataFrame(rows, columns=["blue", "green", "red", "nir"])
signatures["ndvi"] = (signatures["nir"] - signatures["red"]) / (signatures["nir"] + signatures["red"])

Computing the index from the sampled values reproduces exactly what a full-raster index would give at that pixel, provided the scaling was applied first — the ordering argument from Calculating NDVI Directly from xarray DataArrays.

Common Errors

Every sampled value is identical

The points are in a different CRS, so they all land in the same corner pixel — or outside the raster entirely, returning the fill value. Compare a point’s coordinates with src.bounds before sampling.

Values look plausible but are shifted by one pixel

src.index was called with the arguments reversed, or a row/column pair was used where an x/y pair was expected. index(x, y) takes map coordinates and returns (row, col); the inverse is xy(row, col).

Sampling a remote raster is far slower than expected

Each sample call on an open dataset triggers a ranged read, and points scattered across the scene touch many tiles. Sort the coordinates by tile before sampling, or read one window covering a cluster of points, as described in Reading a COG over S3 Without Downloading.

Frequently Asked Questions

Q: Does sample interpolate between pixels? No. It returns the value of the pixel containing the coordinate, with no interpolation. If you need a smoothly varying value, read a small window around the point and interpolate yourself.

Q: What does sample return for a point outside the raster? The dataset’s nodata value, or zero when no nodata is declared — which is indistinguishable from a real measurement. Check the point against the raster bounds before sampling rather than trying to detect it afterwards.

Q: Is sampling one point at a time slow? Sampling itself is fast; opening the file is not. Open once per raster and pass all coordinates in a single sample call, which lets GDAL group the reads it needs.