How to Extract Raster Values at Point Locations with Rasterio

Problem statement

A common GIS task is to extract a raster value for each point location. For example:

  • 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

To extract raster values with Rasterio:

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

Basic pattern:

import rasterio

points = [(500000, 4200000), (500100, 4200100)]

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

print(values)

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

Step-by-step solution

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)

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)

Build the coordinate list

Create a list of point coordinates in the raster CRS.

points = [
    (500000, 4200000),
    (500250, 4200250),
    (500500, 4200500),
]

Sample values with src.sample()

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

import rasterio

raster_path = "data/elevation.tif"
points = [
    (500000, 4200000),
    (500250, 4200250),
    (500500, 4200500),
]

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]}")

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())

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.head())

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]

print(points_gdf.head())

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)

Explanation

How Rasterio sampling works

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

Edge cases / notes

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 = [
    (500000, 4200000),
    (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))

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.

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.

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.