Raster and Vector Do Not Line Up in Python: How to Fix It

Problem statement

You plot an elevation raster and a boundary layer on the same axes and they are nowhere near each other:

fig, ax = plt.subplots()
show(src, ax=ax)
boundary.plot(ax=ax, facecolor="none", edgecolor="red")

The raster fills the frame and the boundary is a red dot in the corner. Or the boundary is invisible because it is 5,000 km away and matplotlib has zoomed out to fit both.

Or the failure is quieter. They look aligned, and then:

ValueError: Input shapes do not overlap raster

Or worse, no error at all — the clip works, the zonal statistics run, and every number is half a cell out.

There are exactly five causes, and each has a distinct signature. Reading the numbers tells you which one you have in about thirty seconds.

Quick answer

Print both bounding boxes. The shape of the numbers identifies the cause:

import rasterio, geopandas as gpd

with rasterio.open("elevation.tif") as src:
    print("raster crs   ", src.crs)
    print("raster bounds", tuple(round(v, 4) for v in src.bounds))
    print("raster res   ", src.res)

gdf = gpd.read_file("boundary.gpkg")
print("vector crs   ", gdf.crs)
print("vector bounds", tuple(round(v, 4) for v in gdf.total_bounds))

Triage table matching bounding-box symptoms to the five causes of raster and vector misalignment. — The numbers name the cause. Read them before changing anything.

What the bounds look like	Cause	Fix
one in degrees (`-2.7 … 53.9`), one in metres (`320000 … 480000`)	different CRS	`gdf.to_crs(src.crs)`
both in degrees, but one is `(53.3, -2.7, …)`	swapped lat/lon	fix the vector's axis order
both plausible, both same CRS, offset by ~1 cell or less	half-pixel / grid offset	use `src.xy()`, not manual arithmetic
raster bounds are `(0, 0, width, height)`	missing transform	the raster has no georeferencing
raster CRS is `None`	missing CRS	identify it; do not guess

gdf = gdf.to_crs(src.crs)          # fixes the most common case

Step-by-step solution

1. Different CRS — the 80% case

Degrees next to metres is unambiguous:

raster crs    EPSG:27700
raster bounds (320000.0, 405000.0, 420000.0, 480000.0)
vector crs    EPSG:4326
vector bounds (-2.7013, 53.3382, -1.9106, 53.9004)

Reproject the vector, not the raster:

gdf = gdf.to_crs(src.crs)
print(gdf.total_bounds)     # (355012.6, 383104.9, 407881.2, 445902.3) ← now comparable

Reprojecting a raster rebuilds every cell and resamples every value; reprojecting a vector transforms a few hundred coordinates and changes nothing else. Move the cheap one.

Do this inside the with block, or capture src.crs first — a common slip is calling to_crs(src.crs) after the dataset has closed.

2. Missing CRS on one side

raster crs    None

vector crs    None

There is nothing to convert. to_crs on a CRS-less GeoDataFrame raises ValueError: Cannot transform naive geometries — see cannot transform naive geometries.

Do not guess. Assigning a CRS relabels coordinates without moving them, so a wrong guess produces data that is confidently in the wrong place. Infer it from the numbers:

minx, miny, maxx, maxy = gdf.total_bounds
if -180 <= minx <= 180 and -90 <= miny <= 90:
    print("looks geographic — probably EPSG:4326")
elif 0 < minx < 800_000 and 0 < miny < 1_400_000:
    print("could be British National Grid (EPSG:27700) — check against a known layer")
else:
    print(f"projected, unknown: {gdf.total_bounds}")

Then verify against a layer whose CRS you trust:

known = gpd.read_file("uk_outline_27700.gpkg")
gdf = gdf.set_crs(27700, allow_override=True)          # the hypothesis
print(gdf.geometry.iloc[0].within(known.union_all()))  # True → hypothesis holds

set_crs labels; to_crs converts. Mixing them up is its own alignment bug — see set_crs vs to_crs.

3. Swapped latitude and longitude

vector bounds (53.3382, -2.7013, 53.9004, -1.9106)

X values of 53 and Y values of −2 in a geographic CRS means the axes are the wrong way round. EPSG:4326 formally defines axis order as latitude, longitude, and most software ignores that in favour of x, y — so data crossing a formal boundary sometimes arrives swapped.

from shapely.ops import transform

gdf["geometry"] = gdf.geometry.apply(
    lambda g: transform(lambda x, y, z=None: (y, x), g)
)
gdf = gdf.set_crs(4326, allow_override=True)

The full diagnosis, including the cases where it is the CRS definition rather than the data, is in my points plot in the ocean.

4. The raster has no transform

raster bounds (0.0, 0.0, 4000.0, 3000.0)
raster res    (1.0, 1.0)

Bounds that are exactly (0, 0, width, height) with 1.0 resolution mean the identity transform — the file has no georeferencing at all. This happens with rasters written from a NumPy array without a profile, PNG or JPEG inputs, and scanned images.

with rasterio.open("scan.tif") as src:
    print(src.transform)

| 1.00, 0.00, 0.00|
| 0.00, 1.00, 0.00|
| 0.00, 0.00, 1.00|

Note the positive y-scale — another giveaway, since a real north-up raster has a negative one.

If you know the corner coordinates and cell size, build the transform:

from rasterio.transform import from_origin, from_bounds

transform = from_origin(320000, 480000, 25, 25)          # top-left x, top-left y, xres, yres
# or, if you know the full extent:
transform = from_bounds(320000, 405000, 420000, 480000, width=4000, height=3000)

profile = src.profile.copy()
profile.update(transform=transform, crs="EPSG:27700")

If you do not know them, the raster needs georeferencing against control points — rasterio will not invent it, and neither should you.

5. A half-pixel or sub-cell offset

Both inputs share a CRS, the bounds overlap, and everything is out by less than one cell:

raster bounds (320000.0, 405000.0, 420000.0, 480000.0)
vector bounds (320012.5, 405012.5, 419987.5, 479987.5)     ← exactly 12.5 = res/2

`src.bounds` describes cell edges; `src.xy()` returns cell centres. Half a cell apart, by definition.

A difference of exactly res / 2 is diagnostic: someone mixed cell edges with cell centres. The transform's origin is the outer edge of the top-left pixel; src.xy(row, col) returns the pixel centre.

# WRONG — manual arithmetic from the origin gives edges, not centres
x = src.transform.c + col * src.res[0]

# RIGHT
x, y = src.xy(row, col)

The same applies in reverse. Use src.index(x, y) rather than dividing by the resolution yourself.

The other sub-cell cause is two rasters reprojected independently to the same CRS. Each gets its own default grid, offset by an arbitrary fraction of a cell. The fix is to warp one onto the other's exact transform:

from rasterio.warp import reproject, Resampling
import numpy as np

with rasterio.open("template.tif") as tmpl:
    dst_transform, dst_crs, shape = tmpl.transform, tmpl.crs, (tmpl.height, tmpl.width)

with rasterio.open("other.tif") as src:
    out = np.empty((src.count, *shape), dtype=src.dtypes[0])
    reproject(source=src.read(), destination=out,
              src_transform=src.transform, src_crs=src.crs, src_nodata=src.nodata,
              dst_transform=dst_transform, dst_crs=dst_crs, dst_nodata=src.nodata,
              resampling=Resampling.bilinear)

Code examples

Example 1: one function that diagnoses all five causes

import rasterio, geopandas as gpd
from shapely.geometry import box

def diagnose(raster_path, vector_path):
    with rasterio.open(raster_path) as src:
        r_crs, r_bounds, r_res = src.crs, tuple(src.bounds), src.res
        identity = src.transform.a == 1 and src.transform.e == 1 and src.transform.c == 0
    gdf = gpd.read_file(vector_path)
    v_crs, v_bounds = gdf.crs, tuple(gdf.total_bounds)

    print(f"raster  {r_crs}  res {r_res}")
    print(f"        {tuple(round(v, 3) for v in r_bounds)}")
    print(f"vector  {v_crs}")
    print(f"        {tuple(round(v, 3) for v in v_bounds)}")
    print()

    if identity:
        return "raster has the identity transform — no georeferencing"
    if r_crs is None:
        return "raster has no CRS — identify it before aligning"
    if v_crs is None:
        return "vector has no CRS — identify it, do not guess"
    if v_crs != r_crs:
        return f"different CRS — call gdf.to_crs({r_crs})"

    # same CRS from here on
    if v_bounds[0] > v_bounds[1] and abs(v_bounds[1]) <= 90 < abs(v_bounds[0]):
        return "vector axes look swapped (lat in x)"
    if not box(*r_bounds).intersects(box(*v_bounds)):
        return "same CRS but disjoint extents — different areas entirely"

    dx = min(abs(v_bounds[0] - r_bounds[0]), abs(v_bounds[2] - r_bounds[2]))
    if 0 < dx <= r_res[0]:
        half = abs(dx - r_res[0] / 2) < r_res[0] * 0.05
        return ("half-pixel offset — edges vs centres" if half
                else f"sub-cell offset of {dx:.3f} — grids differ")
    return "extents overlap and CRS match — alignment looks fine"

print(diagnose("elevation.tif", "boundary.gpkg"))

raster  EPSG:27700  res (25.0, 25.0)
        (320000.0, 405000.0, 420000.0, 480000.0)
vector  EPSG:4326
        (-2.701, 53.338, -1.911, 53.9)

different CRS — call gdf.to_crs(EPSG:27700)

Ordering matters here. Identity transform and missing CRS are checked first, because they make every later comparison meaningless — a bounds comparison against an ungeoreferenced raster produces a confident, wrong answer about offsets.

Example 2: a plot that makes the misalignment obvious

Numbers diagnose; a picture confirms. Force both onto the same axes with explicit extents:

import matplotlib.pyplot as plt
import rasterio
from rasterio.plot import show
import geopandas as gpd

def overlay(raster_path, gdf, ax=None):
    with rasterio.open(raster_path) as src:
        if gdf.crs != src.crs:
            print(f"reprojecting vector {gdf.crs} → {src.crs}")
            gdf = gdf.to_crs(src.crs)
        fig, ax = plt.subplots(figsize=(9, 9)) if ax is None else (ax.figure, ax)
        show(src, ax=ax, cmap="terrain")
        gdf.boundary.plot(ax=ax, color="red", linewidth=1)

        rb, vb = src.bounds, gdf.total_bounds
        ax.set_xlim(min(rb.left, vb[0]), max(rb.right, vb[2]))
        ax.set_ylim(min(rb.bottom, vb[1]), max(rb.top, vb[3]))
        ax.set_title(f"{src.crs} · raster {src.width}×{src.height} @ {src.res[0]}")
    return ax

overlay("elevation.tif", gpd.read_file("boundary.gpkg"))

The set_xlim/set_ylim over the union of both extents is the point. Without it, matplotlib autoscales to whatever was drawn last, and a vector 5,000 km away can look like it is sitting neatly on the raster because the axes silently rescaled between the two draws. See why is my GeoPandas plot blank for the related failure.

Example 3: a folder-level alignment audit

Before a batch run, check the whole set rather than discovering the odd one out at file 180:

from pathlib import Path
from collections import Counter
import rasterio, geopandas as gpd

def audit_alignment(raster_dir, vector_path):
    gdf = gpd.read_file(vector_path)
    rows = []
    for p in sorted(Path(raster_dir).glob("*.tif")):
        with rasterio.open(p) as src:
            ok_crs = src.crs is not None and src.crs == gdf.crs
            identity = src.transform.a == 1 and src.transform.e == 1
            overlaps = (
                False if src.crs is None or identity
                else rasterio.coords.disjoint_bounds(src.bounds, gdf.to_crs(src.crs).total_bounds) is False
            )
            rows.append({"file": p.name, "crs": str(src.crs), "res": src.res[0],
                         "same_crs": ok_crs, "identity": identity, "overlaps": overlaps})

    print(f"vector CRS: {gdf.crs}")
    print(f"raster CRS: {dict(Counter(r['crs'] for r in rows))}")
    print(f"resolutions: {sorted({r['res'] for r in rows})}")
    for r in rows:
        if not r["overlaps"]:
            reason = ("no georeferencing" if r["identity"]
                      else "no CRS" if r["crs"] == "None"
                      else "different CRS" if not r["same_crs"] else "disjoint extent")
            print(f"  ✗ {r['file']:<24} {reason}")
    return rows

audit_alignment("tiles/", "study_area.gpkg")

vector CRS: EPSG:27700
raster CRS: {'EPSG:27700': 238, 'EPSG:4326': 2}
resolutions: [25.0]
  ✗ SJ99SE.tif               different CRS
  ✗ SK00SW.tif               different CRS

Two files out of 240 — exactly the kind of thing that produces a batch job that runs for six hours and then reports two mystery failures. Catching it in the audit costs seconds. This is the same pre-flight discipline as validating pipeline inputs automatically.

Explanation

Stack showing the raster answer as transform plus CRS and the vector answer as coordinates plus CRS, with the parts that must agree. — Two datasets agree only when every layer of the answer agrees. Each layer has its own failure.

Misalignment always reduces to one question: do these two datasets agree about what the coordinate numbers mean?

A raster's answer is stored in two places. The transform maps array indices to coordinate values, and the CRS says which coordinate system those values belong to. A vector's answer is simpler — it stores coordinates directly and a CRS to interpret them. Alignment requires both parts to match, and the failure modes correspond exactly to the parts that can be missing or wrong.

Wrong CRS means both sides have coordinates, but the numbers refer to different systems. This is loud when the units differ, because 53.9 next to 480000 is impossible to miss. It is dangerously quiet between two projected systems with similar ranges — UTM zone 30N and zone 31N produce eastings in the same numeric range, and data in the wrong zone lands a few hundred kilometres off with no obvious tell.

Missing CRS means one side has coordinates that mean nothing. The instinct is to assign the CRS of the other side, which is a hypothesis stated as a fact. Sometimes it is right. When it is wrong, everything downstream is wrong in a way that no later check will catch, because both datasets now agree — incorrectly.

Missing transform means the raster is an image, not a map. Its cells have positions in the array but not on the earth. This cannot be fixed by any amount of reprojection; it needs georeferencing, which is a different operation requiring control points.

Sub-cell offsets are the interesting case, because both datasets are correct and still disagree. A raster's grid is arbitrary — it was chosen by whoever produced it. Two rasters of the same place, at the same resolution, in the same CRS, will not share a grid unless someone made them. calculate_default_transform picks an origin from the bounding box, so independent reprojections of two tiles produce two different grids. The arrays even have the same shape, so NumPy happily multiplies them cell by cell, comparing slightly different ground in every cell. The only fix is to warp onto a shared explicit transform.

The edge-versus-centre distinction is the last half-cell. A pixel is an area, and it has both an outer boundary and a centre. src.bounds and transform.c/transform.f describe edges; src.xy() describes centres. Both conventions are correct, and code that mixes them is out by exactly res / 2 — which is why that number is such a reliable diagnostic. Use the accessor methods and the question never arises.

Edge cases or notes

src.crs == gdf.crs can be False for equivalent CRS with different WKT. Compare with pyproj.CRS.from_user_input(a).equals(b) when the codes look right but the check fails.
EPSG:4326 versus EPSG:4979 (3-D) compare unequal despite the same horizontal datum. Drop Z or normalise the CRS.
rasterio reads bounds from the transform, so a wrong transform gives wrong bounds with no error.
A negative res[1] does not appear — src.res returns absolute values. Check src.transform.e for the sign.
Reprojecting the raster to fix alignment is usually the wrong call. It resamples every value; move the vector instead.
Web Mercator (EPSG:3857) bounds run to ±20,037,508 — a common source of "the numbers look wrong" when overlaying with a basemap. See how to add a basemap with contextily.
Shapefiles store the CRS in a sidecar .prj. A missing .prj is the most common cause of a CRS-less vector.
gdalinfo and ogrinfo -al -so print the same information from the command line and are faster for a quick check.
rasterio.coords.disjoint_bounds(a, b) is a ready-made overlap test that handles the ordering conventions for you.

Internal links

How to fix CRS mismatch in GeoPandas — the vector-to-vector version of the same problem
The raster data model explained — transform, CRS, and why both are needed
How to reproject a raster in Python with Rasterio — when the raster really must move
How to set, assign and convert CRS in GeoPandas — labelling versus converting
My points plot in the ocean: fixing swapped latitude and longitude — the axis-order case in depth
Cannot transform naive geometries to CRS — what a missing CRS raises
How to clip a raster to a polygon in Python — where "shapes do not overlap" comes from
Coordinate precision and floating point in GIS explained — why sub-cell comparisons are fragile

FAQ

Why do I get "Input shapes do not overlap raster"?

The polygon is in a different CRS from the raster, almost always. Print both bounding boxes: degrees next to metres is the giveaway. Fix with gdf.to_crs(src.crs).

Should I reproject the raster or the vector?

The vector. Reprojecting a raster resamples every cell and changes the values; reprojecting a vector transforms a few hundred coordinates and changes nothing else.

My raster has no CRS. Can I assign the vector's?

Only if you have evidence it is correct. Assigning a CRS relabels coordinates without moving them, so a wrong guess produces data that is confidently in the wrong place. Test the hypothesis against a layer you trust.

Everything is offset by exactly half a cell. Why?

Cell edges are being compared with cell centres. src.bounds and the transform's origin describe edges; src.xy() returns centres. Use the accessors instead of manual arithmetic.

Two rasters are in the same CRS but still do not align.

They have different grids. Reprojecting each independently gives each its own origin. Warp one onto the other's exact transform.

My raster bounds are (0, 0, 4000, 3000). What does that mean?

The identity transform — the file has no georeferencing. It cannot be reprojected into place; it needs georeferencing against known control points.

How do I check alignment for a whole folder at once?

Compare CRS, resolution and overlap per file, as in Example 3. A two-file exception in a 240-file set is exactly the failure that surfaces six hours into a batch run.

Raster and Vector Do Not Line Up in Python: How to Fix It #

Problem statement #

Quick answer #

Step-by-step solution #

1. Different CRS — the 80% case #

2. Missing CRS on one side #

3. Swapped latitude and longitude #

4. The raster has no transform #

5. A half-pixel or sub-cell offset #

Code examples #

Example 1: one function that diagnoses all five causes #

Example 2: a plot that makes the misalignment obvious #

Example 3: a folder-level alignment audit #

Explanation #

Edge cases or notes #

Internal links #

FAQ #

Why do I get "Input shapes do not overlap raster"? #

Should I reproject the raster or the vector? #

My raster has no CRS. Can I assign the vector's? #

Everything is offset by exactly half a cell. Why? #

Two rasters are in the same CRS but still do not align. #

My raster bounds are (0, 0, 4000, 3000). What does that mean? #

How do I check alignment for a whole folder at once? #