Sampling Raster Values at Point Locations in Python

Attaching a surface value to a point — the temperature at a weather station, the elevation at a trailhead, the NDVI at a field trial plot — looks like a one-liner and behaves like one right up until the whole column comes back as -9999. This guide is for anyone joining a point layer to a grid: it sits under Zonal Statistics & Raster Sampling in Spatial Analysis & Advanced Query Techniques, and covers rasterio's sample(), the vectorised xarray & rioxarray equivalent, and when nearest-pixel is too coarse to be honest.

Why This Approach / What Goes Wrong

Point sampling is a lookup, not a search. sample() inverts the raster's affine transform to turn (x, y) into (row, col), floors the result, and reads that one cell. There is no distance metric and no candidate set — the answer is whichever pixel happens to contain the coordinate. That is worth stating plainly because it is a different operation from matching a point to the nearest feature, which is an index-backed search covered in Nearest Neighbor & KD-Tree Search. Confusing the two is why people occasionally reach for a KD-tree to do something a matrix inversion already did.

Because it is a bare lookup, nothing in the call can fail in an obvious way. Three failure modes share one symptom — a column that is uniform, wrong, and exception-free:

Four sample points over a pixel grid and what each one returns A six by five pixel grid represents a raster. Point A sits well inside the cell at row zero column one and returns that cell's value of 31.2. Point B lies exactly on a vertical cell edge; because the row and column indices are floored, the right-hand cell wins. Point C falls on a shaded nodata pixel and returns the sentinel minus 9999 as an ordinary number. Point D lies outside the grid entirely and returns the nodata value, or 0.0 when the file declares none. No exception is raised in any of the four cases. One point, one cell — and three ways to get nothing back lst_utm32n.tif A B C nodata D 6 × 5 cells · shaded cell carries the fill value what src.sample() yields A · inside cell (row 0, col 1) 31.2 · correct B · exactly on a cell edge indices are floored — the right-hand cell wins C · lands on a nodata pixel -9999.0 · an ordinary float D · outside the raster nodata — or 0.0 if the file declares none no exception is raised in any of the four cases Only an explicit bounds test and an explicit nodata test tell C and D apart from a real reading.
The lookup always succeeds; what varies is whether the number it hands back means anything.

Prerequisites

python -m pip install \
  "rasterio>=1.3,<2" "geopandas>=1.0" "rioxarray>=0.15" "xarray>=2024.3" \
  "rasterstats>=0.19" "scipy>=1.11" "numpy>=1.26" "pandas>=2.0"

Step-by-Step Implementation

1. Open both sides and read the raster's metadata while the handle is alive. Capture crs, nodata and bounds into locals — reaching for them after the with block raises Attempt to read from a closed dataset.

import geopandas as gpd
import numpy as np
import rasterio

sensors = gpd.read_file("air_quality_sensors.geojson")      # 240 points, EPSG:4326

with rasterio.open("lst_utm32n.tif") as src:                # land surface temperature
    grid_crs, fill, bounds = src.crs, src.nodata, src.bounds

print(sensors.crs, "->", grid_crs, "| nodata:", fill)
# EPSG:4326 -> EPSG:32632 | nodata: -9999.0

2. Move the points into the raster's CRS — the step everyone skips. Reprojecting a few thousand vertices is exact and cheap; reprojecting the grid resamples every pixel and changes the values you are about to read. Let the raster define the working CRS.

if sensors.crs is None:
    raise ValueError("sensors carry no CRS — set_crs() before sampling anything")
if not sensors.crs.equals(grid_crs):
    sensors = sensors.to_crs(grid_crs)                      # move the points, not the pixels

coords = np.column_stack([sensors.geometry.x, sensors.geometry.y])   # (n, 2) as (x, y)

sample() wants (easting, northing) — x first. GeoPandas always hands you x from .geometry.x, so the ordering is safe here, but building coordinates by hand from a pyproj transformer is not: Transformer.from_crs("EPSG:4326", grid_crs) follows the authority axis order and returns (lat, lon). Pass always_xy=True and the tuple comes back the way sample() expects.

The five stages of a point-sampling pipeline and the skipped reprojection Five stages run left to right: a sensors layer in EPSG:4326, a CRS alignment stage calling to_crs on the raster CRS, a bounds test that drops points outside the scene, the sample call that reads the nearest cell for every band, and a final stage converting the nodata sentinel to NaN. A dashed red path skips the alignment stage and feeds the raw degree coordinates straight into sample, where every point falls outside the raster bounds and the whole column comes back as nodata with no exception. The sample pipeline — stage 2 is the one that gets skipped 1 · sensors EPSG:4326 240 points 2 · align CRS to_crs(grid_crs) degrees → metres 3 · bounds test src.bounds flag off-scene points 4 · sample src.sample(xy) nearest cell, all bands 5 · fill → NaN isclose(v, fill) column ready to join skip stage 2 → every point falls outside src.bounds and the entire column comes back as the fill value, silently Stages 2 and 3 cost microseconds; without them the other three produce a plausible column of nothing.
Alignment and the bounds test are not defensive extras — they are the only two stages that can tell you the result is real.

3. Flag the off-grid points before sampling, then sample the rest. Keeping the mask separate preserves the row count of the table: unsampled rows get NaN, never a silent zero.

left, bottom, right, top = bounds
on_grid = (
    (coords[:, 0] >= left) & (coords[:, 0] < right)
    & (coords[:, 1] > bottom) & (coords[:, 1] <= top)
)

values = np.full(len(sensors), np.nan, dtype="float64")
with rasterio.open("lst_utm32n.tif") as src:
    sampled = np.array(
        [record[0] for record in src.sample(coords[on_grid], indexes=1)],
        dtype="float64",
    )

if fill is not None:
    sampled[np.isclose(sampled, fill)] = np.nan     # a NaN fill is already NaN

values[on_grid] = sampled
sensors["lst_c"] = values

sample() is a generator that performs one 1×1 windowed read per coordinate, so it never materialises the scene — the same property that makes it usable against a remote cloud-optimized GeoTIFF. Passing indexes=1 still yields a one-element array per point, hence the record[0].

4. Take every band in the same pass. For a 12-month composite, ask for all bands at once: each point costs one windowed read covering the full stack rather than twelve separate traversals. The band semantics of such files are covered in Reading Multi-Band TIFFs with Rasterio.

import pandas as pd

with rasterio.open("lst_monthly_utm32n.tif") as src:        # 12 bands, one per month
    band_ids = list(range(1, src.count + 1))
    stack = np.vstack(list(src.sample(coords, indexes=band_ids))).astype("float32")
    monthly_fill = src.nodata

if monthly_fill is not None:
    stack[np.isclose(stack, monthly_fill)] = np.nan

monthly = pd.DataFrame(
    stack,                                                   # shape (n_points, 12)
    columns=[f"lst_m{m:02d}" for m in band_ids],
    index=sensors.index,
)
sensors = sensors.join(monthly)

5. Go vectorised when the raster is already a cube. rioxarray opens the file as a labelled array, and xarray's advanced indexing does the whole point set in one call. The rule that makes it pointwise rather than a cross product is that both indexers share a dimension name.

import rioxarray  # noqa: F401 — registers the .rio accessor
import xarray as xr

lst = rioxarray.open_rasterio("lst_utm32n.tif", masked=True).squeeze("band", drop=True)
pts = sensors.to_crs(lst.rio.crs)

xs = xr.DataArray(pts.geometry.x.to_numpy(), dims="point")   # same dim name …
ys = xr.DataArray(pts.geometry.y.to_numpy(), dims="point")   # … on both indexers

half_cell = abs(lst.rio.resolution()[0]) / 2
nearest = lst.sel(x=xs, y=ys, method="nearest", tolerance=half_cell)
sensors["lst_c"] = nearest.to_numpy()                        # shape (n_points,)

masked=True turns the declared nodata into NaN on read, which removes step 3's sentinel arithmetic entirely. tolerance=half_cell is the part worth copying: without it, a point beyond the edge of the grid snaps to the nearest edge cell and returns a plausible value from the wrong place; with it, xarray raises KeyError: not all values found in index 'x'. Loud beats plausible.

Separate dimension names take a cross product, shared names index pointwise On the left, passing plain lists of x and y values to sel produces the outer product: a 240 by 240 result of 57 600 values, of which only the 240 on the diagonal were ever wanted, and memory grows with the square of the point count. On the right, wrapping both coordinate arrays in DataArrays that share the dimension name point produces a single row of 240 values, one per sample location. The dimension name decides the shape of the answer plain lists → cross product lst.sel(x=x_list, y=y_list, ...) 240 × 240 = 57 600 values 240 wanted memory grows with the square of the point count shared dims="point" → pointwise xs = xr.DataArray(x, dims="point") 240 values · one per point identical dim name on both indexers add tolerance= to reject off-grid points Same call, same data — only the dimension names differ between a 240-element column and a 57 600-cell matrix.
Naming both indexer dimensions "point" is what collapses the cross product into one value per location.

6. Interpolate when the staircase is a lie. Nearest-pixel sampling is correct for categorical grids — land cover, soil class, zoning — where a blend of class 20 and class 40 is class 30, which may not exist. For a smooth field it produces a visible step at every cell boundary, and two stations 30 m apart on a 250 m grid report exactly the same temperature. Bilinear interpolation weights the four surrounding cell centres by distance instead.

smooth = (
    rioxarray.open_rasterio("lst_utm32n.tif", masked=True)
    .squeeze("band", drop=True)
    .sortby("y")                                             # ascending axis for scipy
    .interp(x=xs, y=ys, method="linear")                     # bilinear, pointwise
)
sensors["lst_bilinear"] = smooth.to_numpy()

The one-line alternative, if you would rather not open a cube, is rasterstats.point_query(sensors, "lst_utm32n.tif", interpolate="bilinear", nodata=-9999), which returns a plain list of values with None where a point is off the grid — bilinear is its default, so pass interpolate="nearest" explicitly for categorical data.

Bilinear weighting of four cell centres versus a nearest-cell transect On the left, a sample point sits inside the square formed by four neighbouring cell centres holding the values 12.4, 18.6, 13.1 and 20.2. Its distance to each centre gives weights of 0.22, 0.45, 0.11 and 0.22, so nearest returns 18.6 while bilinear returns 17.0. On the right, a transect across the same smooth surface shows the nearest-cell result as a staircase that jumps a whole cell value at every boundary, and the bilinear result as a continuous line through the cell centres. Nearest cell vs bilinear on a smooth surface the four surrounding cell centres w = 0.22 12.4 w = 0.45 18.6 13.1 w = 0.11 20.2 w = 0.22 the sample point falls in the shaded top-right cell nearest → 18.6 bilinear → 17.0 a transect across the same surface value nearest — jumps a whole cell at every edge bilinear — continuous through the centres Bilinear needs four valid neighbours — one nodata cell in the window yields NaN, which is the honest answer.
Weights come from the point's distance to each of the four centres; the nearest-cell answer is whichever centre happens to be closest, whatever the other three say.

Verification

Two independent checks earn their keep: the row count must survive, and sample() must agree with a hand-rolled index() plus windowed read on a point you know is inside the grid.

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

with rasterio.open("lst_utm32n.tif") as src:
    x0, y0 = coords[on_grid][0]
    row, col = src.index(x0, y0)                      # floor of the inverse transform
    manual = src.read(1, window=Window(col, row, 1, 1))[0, 0]

assert len(sensors) == len(coords), "rows were dropped — sampling must not filter the table"
assert np.isclose(manual, sensors.loc[on_grid, "lst_c"].iloc[0]), "sample() != index()+read()"
assert not np.isclose(sensors["lst_c"].dropna(), -9999.0).any(), "a fill value survived"
assert sensors["lst_c"].notna().any(), "every value is NaN — CRS or extent mismatch"

print(f"{int((~on_grid).sum())} of {len(sensors)} sensors fell outside the scene")
# 3 of 240 sensors fell outside the scene
print(sensors[["sensor_id", "lst_c", "lst_bilinear"]].head(3).to_string(index=False))
# sensor_id  lst_c  lst_bilinear
#     A-011  31.20         30.94
#     A-012  29.80         30.05
#     A-013    NaN           NaN

The last assertion is the cheapest CRS alarm you can install: an all-NaN column after a bounds mask means the points and the pixels never occupied the same numeric space.

Edge Cases & Debugging

Frequently Asked Questions

Does rasterio.sample() interpolate at all? No. It returns the value of the single cell containing the coordinate, with no smoothing and no weighting, which is exactly right for categorical grids and coarse for continuous ones. Bilinear requires either rasterstats.point_query(..., interpolate="bilinear") or DataArray.interp(method="linear"); there is no interpolation flag on sample() itself.

What happens to points that fall outside the raster? Nothing visible. rasterio yields the dataset's nodata value for out-of-bounds coordinates and substitutes 0 when the file declares none, so an off-scene point is indistinguishable from a real measurement unless you test the coordinate against src.bounds yourself. The xarray path has the better ergonomics here: method="nearest" with a tolerance of half a cell raises a KeyError rather than snapping to the edge.

Is sample() fast enough for a million points? It is a Python generator issuing one windowed read per coordinate, so it handles tens of thousands of points comfortably and becomes the bottleneck somewhere in the hundreds of thousands — especially over HTTP. Past that, switch strategy: read one window covering the whole point cloud, convert coordinates to array indices in NumPy, and index the block directly, or use the vectorised xarray selection from step 5, which is a single indexing operation over a chunked cube. The trade-offs between the two engines are laid out in xarray vs rasterio for Time-Series Rasters.

Should I ever reproject the raster to match the points instead? Only when several rasters must be forced onto one grid anyway. Reprojecting points is an exact coordinate transform; reprojecting a grid resamples every pixel and changes the values you are about to sample — and if the raster is categorical, anything but nearest-neighbour resampling invents classes that were never observed, as covered in Resampling & Overviews When Writing COGs.

Related