How to Extract Raster Values at Point Locations with Rasterio

On this page

Problem statement

Sampling a raster at points — the elevation at Snowdon's summit, the land cover at a monitoring site — returns the value of the one cell under each point, provided the points are in the raster's CRS. It is the step behind:

get elevation values from a DEM at survey points
extract land cover classes for monitoring sites
read temperature values at station coordinates
join raster cell values to point features for analysis or export

In Python, the standard way to do this with Rasterio is to sample the raster at point coordinates. The main requirement is that the point coordinates must be in the same CRS as the raster.

Quick answer

Sampling returns the value of the single raster cell beneath each point.

To extract raster values with Rasterio:

open the raster with rasterio.open()
make sure point coordinates are in the raster CRS
pass point coordinates to src.sample()
store the returned values in a list, table, or GeoDataFrame

Basic pattern:

import rasterio

points = [(260986, 354378), (264182, 357793)]      # Snowdon and Glyder Fawr, British National Grid

with rasterio.open("data/elevation.tif") as src:
    values = list(src.sample(points))

print(values)

[array([1066.5554], dtype=float32), array([947.44055], dtype=float32)]

For a single-band raster, each result contains one value. For a multi-band raster, each point returns one value per band.

Sampling raster values at points, step by step

Move a point a metre across a pixel boundary and the sampled value jumps.

The data used on this page

elevation.tif is the Copernicus GLO-30 digital surface model for Snowdonia, four 1° tiles merged and reprojected to British National Grid at 30 m — a projected DEM, as most national ones are. sample_points.geojson holds the 260 named peaks and hills GeoNames lists in the same area, in longitude and latitude, with their published heights. multiband.tif is a four-band Sentinel-2 L2A window around Snowdon.

import io
import os
import zipfile

import geopandas as gpd
import numpy as np
import pandas as pd
import rasterio
import requests
from rasterio.merge import merge
from rasterio.warp import Resampling, calculate_default_transform, reproject
from rasterio.windows import from_bounds

os.environ.update(AWS_NO_SIGN_REQUEST="YES", GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR")
os.makedirs("data", exist_ok=True)
os.makedirs("output", exist_ok=True)
BOX = (-4.3, 52.9, -3.7, 53.25)

if not os.path.exists("data/elevation.tif"):
    tile = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_{0}_DEM/"
            "Copernicus_DSM_COG_10_{0}_DEM.tif")
    sources = [rasterio.open(tile.format(t)) for t in
               ("N52_00_W005_00", "N52_00_W004_00", "N53_00_W005_00", "N53_00_W004_00")]
    mosaic, transform = merge(sources, bounds=BOX)
    dst_transform, width, height = calculate_default_transform(
        "EPSG:4326", "EPSG:27700", mosaic.shape[2], mosaic.shape[1], *BOX, resolution=30)
    dem = np.full((height, width), -9999, "float32")
    reproject(mosaic[0], dem, src_transform=transform, src_crs="EPSG:4326", dst_transform=dst_transform,
              dst_crs="EPSG:27700", resampling=Resampling.bilinear, dst_nodata=-9999)
    with rasterio.open("data/elevation.tif", "w", driver="GTiff", width=width, height=height, count=1,
                       dtype="float32", crs="EPSG:27700", transform=dst_transform, nodata=-9999,
                       compress="deflate", tiled=True) as dst:
        dst.write(dem, 1)

    names = ["id", "name", "ascii", "alt", "lat", "lon", "fclass", "fcode", "cc", "cc2", "a1", "a2",
             "a3", "a4", "pop", "elev", "dem", "tz", "mod"]
    gb = pd.read_csv(zipfile.ZipFile(io.BytesIO(requests.get(
        "https://download.geonames.org/export/dump/GB.zip", timeout=120).content)).open("GB.txt"),
        sep="\t", header=None, names=names, quoting=3, low_memory=False)
    peaks = gb[gb["fcode"].isin(["MT", "PK", "HLL"]) & gb["lon"].between(BOX[0], BOX[2])
               & gb["lat"].between(BOX[1], BOX[3])]
    gpd.GeoDataFrame(peaks[["name", "elev"]], crs=4326,
                     geometry=gpd.points_from_xy(peaks["lon"], peaks["lat"])).to_file(
        "data/sample_points.geojson")

    item = requests.get("https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/items/"
                        "S2A_30UVD_20250922_0_L2A", timeout=60).json()
    window_bounds = (422880, 5875440, 432880, 5885440)            # EPSG:32630
    bands = []
    for key in ("blue", "green", "red", "nir"):
        with rasterio.open(item["assets"][key]["href"]) as src:
            window = from_bounds(*window_bounds, transform=src.transform)
            bands.append(src.read(1, window=window))
            profile = dict(crs=src.crs, transform=src.window_transform(window))
    with rasterio.open("data/multiband.tif", "w", driver="GTiff", width=bands[0].shape[1],
                       height=bands[0].shape[0], count=4, dtype="uint16", nodata=0,
                       compress="deflate", **profile) as dst:
        dst.write(np.stack(bands))

with rasterio.open("data/elevation.tif") as src:
    print("elevation.tif:", src.width, "x", src.height, src.crs, src.res)
print("sample_points.geojson:", len(gpd.read_file("data/sample_points.geojson")), "peaks")
with rasterio.open("data/multiband.tif") as src:
    print("multiband.tif:", src.count, "bands", src.width, "x", src.height, src.crs)

elevation.tif: 1376 x 1335 EPSG:27700 (30.0, 30.0)
sample_points.geojson: 260 peaks
multiband.tif: 4 bands 1000 x 1000 EPSG:32630

Check the raster and point data before sampling

Confirm the raster CRS

Before you sample, inspect the raster metadata. Point coordinates must match the raster CRS.

import rasterio

with rasterio.open("data/elevation.tif") as src:
    print("CRS:", src.crs)
    print("Band count:", src.count)
    print("Bounds:", src.bounds)
    print("NoData:", src.nodata)

CRS: EPSG:27700
Band count: 1
Bounds: BoundingBox(left=245393.66384992778, bottom=334972.9179518375, right=286673.6638499278, top=375022.9179518375)
NoData: -9999.0

If your raster is in a projected CRS such as UTM, do not sample it using longitude/latitude coordinates.

Confirm the point coordinate format

Point locations can come from:

a Python list of (x, y) tuples
a CSV with coordinate columns
a GeoJSON file
a shapefile or GeoPackage

Be clear about whether the coordinates are:

longitude/latitude in EPSG:4326
projected coordinates such as meters in UTM or Web Mercator

Decide what output you need

Typical outputs are:

one value per point for a single-band raster
one value per band for a multi-band raster
a new attribute column in a GeoDataFrame
a CSV for non-spatial reporting

Extract raster values from a list of coordinates with Rasterio

Open the raster file

This example shows a single-band GeoTIFF.

import rasterio

raster_path = "data/elevation.tif"

with rasterio.open(raster_path) as src:
    print(src.crs)
    print(src.count)
    print(src.nodata)

EPSG:27700
1
-9999.0

Build the coordinate list

Create a list of point coordinates in the raster CRS.

points = [
    (260986, 354378),     # Snowdon, 1,085 m
    (264182, 357793),     # Glyder Fawr, 999 m
    (258000, 360000),     # the Llanberis pass
]

Sample values with `src.sample()`

Use src.sample() to get the raster value at each point.

import rasterio

raster_path = "data/elevation.tif"
points = [
    (260986, 354378),
    (264182, 357793),
    (258000, 360000),
]

with rasterio.open(raster_path) as src:
    samples = list(src.sample(points))

for point, value in zip(points, samples):
    print(f"Point {point} -> {value[0]:.1f}")

Point (260986, 354378) -> 1066.6
Point (264182, 357793) -> 947.4
Point (258000, 360000) -> 107.2

For a single-band raster, each sampled result is an array with one element, so value[0] gives the pixel value.

Extract raster values for point features in a GeoPandas GeoDataFrame

Read point data with GeoPandas

If your points are stored in GeoJSON or a shapefile, load them with GeoPandas.

import geopandas as gpd

points_gdf = gpd.read_file("data/sample_points.geojson")
print(points_gdf.crs)
print(points_gdf.head())

EPSG:4326
                  name    elev                   geometry
0              Yr Aran   747.0  POINT (-4.08318 53.04277)
1          Moel Tryfan   429.0  POINT (-4.21811 53.08226)
2  Snowdon - Yr Wyddfa  1085.0   POINT (-4.07623 53.0685)
3         Moel Seisiog   532.0      POINT (-3.7 53.08333)
4         Moelwyn Mawr   770.0  POINT (-4.00014 52.98425)

Reproject points to the raster CRS if needed

This is the most important step when point CRS and raster CRS differ.

import geopandas as gpd
import rasterio

points_gdf = gpd.read_file("data/sample_points.geojson")

with rasterio.open("data/elevation.tif") as src:
    raster_crs = src.crs

if points_gdf.crs != raster_crs:
    points_gdf = points_gdf.to_crs(raster_crs)

Filter to valid point geometries before sampling

Before extracting coordinates, remove null geometries and non-point features.

points_gdf = points_gdf[points_gdf.geometry.notnull()]
points_gdf = points_gdf[points_gdf.geometry.geom_type == "Point"]

Add sampled raster values as a new column

Extract the point coordinates from the geometry column, sample the raster, and store the result in a new attribute field.

import geopandas as gpd
import rasterio
import numpy as np

points_gdf = gpd.read_file("data/sample_points.geojson")
points_gdf = points_gdf[points_gdf.geometry.notnull()]
points_gdf = points_gdf[points_gdf.geometry.geom_type == "Point"]

with rasterio.open("data/elevation.tif") as src:
    if points_gdf.crs != src.crs:
        points_gdf = points_gdf.to_crs(src.crs)

    coords = [(geom.x, geom.y) for geom in points_gdf.geometry]
    sampled = list(src.sample(coords))

    values = []
    for item in sampled:
        val = item[0]
        if src.nodata is not None and val == src.nodata:
            values.append(np.nan)
        else:
            values.append(val)

points_gdf["raster_value"] = values

print(points_gdf.sort_values("elev", ascending=False).head().to_string(index=False))

               name   elev                      geometry  raster_value
Snowdon - Yr Wyddfa 1085.0 POINT (260985.655 354378.055)   1066.555420
      Carnedd Ugain 1065.0 POINT (261072.686 355160.336)   1055.745605
   Carnedd Llewelyn 1064.0 POINT (268360.973 364374.029)   1058.298096
     Carnedd Dafydd 1044.0 POINT (266311.961 363067.271)   1040.617920
        Glyder Fawr  999.0  POINT (264182.097 357792.95)    947.440552

elev is the height GeoNames publishes for each summit; raster_value is the DEM cell under it. The next section explains why the second is lower.

This is the standard workflow for extracting raster cell values to point features.

Handle single-band and multi-band rasters

Single-band raster output

With one band, each sampled item contains one value.

single_band_values = [item[0] for item in sampled]

This is the normal approach to read a raster value by coordinate in Python.

Multi-band raster output

For a multi-band raster, each point returns one value per band.

import geopandas as gpd
import rasterio

points_gdf = gpd.read_file("data/sample_points.geojson")
points_gdf = points_gdf[points_gdf.geometry.notnull()]
points_gdf = points_gdf[points_gdf.geometry.geom_type == "Point"]

with rasterio.open("data/multiband.tif") as src:
    if points_gdf.crs != src.crs:
        points_gdf = points_gdf.to_crs(src.crs)

    coords = [(geom.x, geom.y) for geom in points_gdf.geometry]
    sampled = list(src.sample(coords))

if sampled:
    for band_index in range(len(sampled[0])):
        points_gdf[f"band_{band_index + 1}"] = [row[band_index] for row in sampled]

inside = points_gdf[points_gdf["band_4"] > 0]
print(f"{len(inside)} of {len(points_gdf)} peaks inside the image")
print(inside[["name", "band_1", "band_2", "band_3", "band_4"]].head(4).to_string(index=False))

22 of 260 peaks inside the image
               name  band_1  band_2  band_3  band_4
            Yr Aran    1250    1376    1284    3672
Snowdon - Yr Wyddfa    2028    2120    2094    3604
        Glyder Fawr     908    1070    1194    2304
        Glyder Fach    1116    1238    1328    1882

The multiband raster is a 10 km Sentinel-2 window around Snowdon — blue, green, red and near-infrared at 10 m, from 22 September 2025 — so most of the 260 peaks fall outside it and sample as 0, the image's nodata. The bands are still raw digital numbers; divide by 10,000 for reflectance.

This is useful when sampling imagery or stacked environmental layers.

Save the results for later GIS use

Export points with new raster value fields

After sampling, save the updated GeoDataFrame.

points_gdf.to_file("output/points_with_values.gpkg", driver="GPKG")
points_gdf.to_file("output/points_with_values.geojson", driver="GeoJSON")

Export a non-spatial table

If you only need attributes, write a CSV.

points_gdf.drop(columns="geometry").to_csv("output/points_with_values.csv", index=False)

How sampling works and why CRS matters

How Rasterio sampling works

A CRS mismatch does not raise here either — it just samples empty space.

src.sample() takes an iterable of (x, y) coordinates and returns the value of the raster cell for each coordinate. It is a point sampling method, not a zonal statistic or neighborhood average.

If you need interpolation or summary values from surrounding cells, use a different workflow.

Why CRS mismatches cause wrong values

If your raster is in UTM and your points are in EPSG:4326, the coordinates refer to different places. Sampling will either return values from the wrong cells or values outside the intended area.

For more background, see How to Reproject Spatial Data in Python (GeoPandas).

When sampled values may be missing or wrong

Unexpected values usually come from:

points outside the raster extent
cells marked as nodata
wrong CRS
null or non-point geometries in the input data

Measured on the Snowdonia peaks:

import time

import numpy as np

peaks = gpd.read_file("data/sample_points.geojson")
with rasterio.open("data/elevation.tif") as src:
    bng = peaks.to_crs(src.crs)
    dem = np.array([v[0] for v in src.sample([(g.x, g.y) for g in bng.geometry])])
    difference = dem - peaks["elev"].to_numpy()
    print(f"DEM minus published height: median {np.median(difference):+.1f} m, "
          f"quartiles {np.percentile(difference, 25):+.1f} / {np.percentile(difference, 75):+.1f} m")
    glyder = bng[bng["name"] == "Glyder Fawr"].geometry.iloc[0]
    around = src.read(1, window=rasterio.windows.from_bounds(
        glyder.x - 150, glyder.y - 150, glyder.x + 150, glyder.y + 150, src.transform))
    print(f"Glyder Fawr: cell {dem[bng['name'] == 'Glyder Fawr'][0]:.1f} m, highest cell within "
          f"150 m {around.max():.1f} m, published 999 m, GeoNames position "
          f"{peaks[peaks['name'] == 'Glyder Fawr'].geometry.iloc[0].coords[0]}")

    wrong = np.array([v[0] for v in src.sample([(g.x, g.y) for g in peaks.geometry])])
    print("lon/lat sampled against a BNG raster:", np.unique(wrong))
    print("outside the extent:", [v[0] for v in src.sample([(0, 0)])],
          "| masked=True:", [v.mask[0] for v in src.sample([(0, 0)], masked=True)])

    rng = np.random.default_rng(0)
    xs = rng.uniform(src.bounds.left, src.bounds.right, 100_000)
    ys = rng.uniform(src.bounds.bottom, src.bounds.top, 100_000)
    start = time.perf_counter()
    one_by_one = np.array([v[0] for v in src.sample(zip(xs, ys))])
    t_sample = time.perf_counter() - start
    start = time.perf_counter()
    rows, cols = rasterio.transform.rowcol(src.transform, xs, ys)
    all_at_once = src.read(1)[rows, cols]
    t_array = time.perf_counter() - start
print(f"100,000 points: sample() {t_sample:.2f} s, read + index {t_array:.3f} s, "
      f"same values: {np.array_equal(one_by_one, all_at_once)}")

DEM minus published height: median -8.0 m, quartiles -13.5 / -4.4 m
Glyder Fawr: cell 947.4 m, highest cell within 150 m 989.8 m, published 999 m, GeoNames position (-4.03, 53.1)
lon/lat sampled against a BNG raster: [-9999.]
outside the extent: [np.float32(-9999.0)] | masked=True: [np.True_]
100,000 points: sample() 1.45 s, read + index 0.019 s, same values: True

Four lessons. A DEM cell is an average over 30 m, and a summit is the highest point of it, so sampled peaks come out low — Snowdon's cell is 18 m under the published 1,085 m. A point sample is also only as good as the point: Glyder Fawr's cell is 52 m low because GeoNames places it at a position rounded to hundredths of a degree, on the slope, while the highest cell within 150 m is 10 m under the summit. Points in the wrong CRS do not raise: all 260 land outside the raster and return its nodata value, which is the "every value is the same" symptom. And sample() is a loop over points; when the raster fits in memory, reading the band once and indexing it is around a hundred times faster and returns the same numbers.

Glyder Fawr is the outlier because its point is on the slope, not because its cell is wrong.

Points outside the extent, nodata and large batches

Points outside the raster extent

If a point falls outside src.bounds, handle it explicitly before or after sampling. A bounds check is the safest option if you need reliable control over outside points.

import rasterio

points = [
    (260986, 354378),
    (900000, 9999999),
]

with rasterio.open("data/elevation.tif") as src:
    valid_points = []
    for x, y in points:
        if src.bounds.left <= x <= src.bounds.right and src.bounds.bottom <= y <= src.bounds.top:
            valid_points.append((x, y))

    samples = list(src.sample(valid_points))
print(f"{len(valid_points)} of {len(points)} points inside; samples {samples}")

1 of 2 points inside; samples [array([1066.5554], dtype=float32)]

Outside points should be checked explicitly, because how they appear in output depends on your sampling and nodata handling.

Nodata values in the raster

If the raster uses a nodata value, convert it to None or NaN after sampling so your output is easier to analyze.

Invalid geometries

This workflow expects point geometries. If your dataset contains null geometries or non-point features, filter them before calling .x and .y.

points_gdf = points_gdf[points_gdf.geometry.notnull()]
points_gdf = points_gdf[points_gdf.geometry.geom_type == "Point"]

Sampling large point datasets

For large point files, avoid repeated raster opens inside loops. Open the raster once, build the coordinate list once, and sample in one pass.

Raster resolution and pixel alignment

The extracted value depends on raster cell size and grid alignment. Two nearby points can return different values if they fall in different cells.

Internal links

For CRS background and reprojection, see How to Reproject Spatial Data in Python (GeoPandas).

If your points come from a shapefile, read How to Read a Shapefile in Python with GeoPandas.

If you want to save the sampled points as GeoJSON, see How to Export GeoJSON in Python with GeoPandas.

If your sampled values look wrong, first check CRS, geometry type, nodata, and raster bounds.

How to read a USGS ShakeMap and find who was shaken — sampling an intensity grid at exposure units
How to extract climate covariates at occurrence points — the same sampling, for species records

FAQ

How do I extract raster values at point locations with Rasterio?

Open the raster with rasterio.open(), prepare point coordinates in the raster CRS, and pass them to src.sample(). For point files, read them with GeoPandas, reproject if needed, then add the sampled values as new columns.

Do the point coordinates need to match the raster CRS?

Yes. If the raster and points use different CRSs, reproject the points to the raster CRS before sampling.

How do I handle nodata values when sampling a raster?

Check src.nodata and replace matching sampled values with None or numpy.nan. That makes later filtering and statistics easier.

Can Rasterio extract values for multiple bands at the same point?

Yes. For a multi-band raster, each point returns one value per band. Store them in separate fields such as band_1, band_2, and band_3.

What happens if a point falls outside the raster extent?

Check outside points explicitly with src.bounds if you need predictable handling. That is safer than assuming how those locations will appear in sampled output.

How do I sample only one band from a multi-band raster?

Pass indexes= to sample(), for example src.sample(coords, indexes=1) for the first band. Each returned array then holds only that band's value, which avoids reading every band when you need just one.

Can I sample a raster directly from point coordinates without GeoPandas?

Yes. src.sample() accepts any iterable of (x, y) tuples, so a plain Python list or a NumPy array of coordinates works. GeoPandas is only a convenience for reading point files and handling their CRS.

How do I get the row and column of the cell a point falls in?

Use src.index(x, y), which returns the (row, col) of the cell containing that coordinate. This is useful when you want to read or modify the array directly rather than only sampling values.

Why are my sampled elevation values all the same number?

This usually means the points landed outside the raster extent or on nodata cells, so sample() returned the raster's fill or nodata value. Check each point against src.bounds and compare results to src.nodata before trusting them.

Is `src.sample()` slow for hundreds of thousands of points?

It reads each cell on demand: measured, 100,000 points took well over a second, where reading the band once and indexing the array took around a hundredth of that and returned the same values. Open the raster once and pass the full coordinate list in one call; if the raster fits in memory, index the array instead.