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

On this page

Problem statement

You plot an elevation raster and a boundary layer on the same axes and they are nowhere near each other. Here the raster is 30 m Copernicus elevation for the southern Pennines and the boundary is the Calderdale district (the setup is under The data used on this page):

import matplotlib
matplotlib.use("Agg")
import geopandas as gpd
import matplotlib.pyplot as plt
import rasterio
from rasterio.plot import show

src = rasterio.open("elevation.tif")
boundary = gpd.read_file("boundary.gpkg")
fig, ax = plt.subplots()
show(src, ax=ax)
boundary.plot(ax=ax, facecolor="none", edgecolor="red")
print(f"axes span x {ax.get_xlim()[0]:,.1f} to {ax.get_xlim()[1]:,.1f}; "
      f"the whole district is {boundary.total_bounds[2] - boundary.total_bounds[0]:.2f} units wide")

axes span x -23,361.9 to 467,192.1; the whole district is 0.45 units wide

The raster fills a frame hundreds of kilometres wide, and the district — a fraction of one unit across — is an invisible dot at the origin.

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:

with rasterio.open("elevation.tif") as src:
    print("raster crs   ", src.crs)
    print("raster bounds", tuple(round(v, 1) 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(float(v), 4) for v in gdf.total_bounds))

raster crs    EPSG:27700
raster bounds (332972.1, 344725.2, 467192.1, 456475.2)
raster res    (30.0, 30.0)
vector crs    EPSG:4326
vector bounds (-2.1733, 53.6153, -1.7272, 53.8256)

The numbers name the cause. Read them before changing anything.

What the bounds look like	Cause	Fix
one in degrees (`-2.2 … 53.8`), one in metres (six digits)	different CRS	`gdf.to_crs(src.crs)`
both in degrees, but one is `(53.6, -2.2, …)`	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

gdf = gdf.to_crs(src.crs)          # fixes the most common case
print("vector bounds", tuple(round(float(v)) for v in gdf.total_bounds))

vector bounds (388663, 413197, 418114, 436596)

Different CRS, missing CRS, swapped axes and more

1. Different CRS — the common case

Degrees next to metres is unambiguous, as the quick answer's printout shows. Reproject the vector, not the raster. 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

A copy of the boundary saved without its CRS — the state a shapefile arrives in when its .prj is lost:

naive = gpd.read_file("boundary_no_crs.gpkg")
print("vector crs   ", naive.crs)
try:
    naive.to_crs(27700)
except ValueError as err:
    print(f"ValueError: {err}")

vector crs    None
ValueError: Cannot transform naive geometries.  Please set a crs on the object first.

There is nothing to convert. 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 = naive.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: {naive.total_bounds}")

looks geographic — probably EPSG:4326

Then verify against a layer whose CRS you trust — here the ONS outline of Great Britain:

known = gpd.read_file("gb_outline_27700.gpkg")
guess = naive.set_crs(4326).to_crs(27700)                 # the hypothesis
print("inside Great Britain:", bool(guess.geometry.iloc[0].within(known.union_all().buffer(1000))))
wrong = naive.set_crs(27700)                              # a wrong guess, for contrast
print("as if it were metres:", bool(wrong.geometry.iloc[0].within(known.union_all())))

inside Great Britain: True
as if it were metres: False

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

A copy with its axes swapped, as data crossing a strict EPSG:4326 boundary sometimes arrives:

swapped = gpd.read_file("boundary_swapped.gpkg")
print("vector bounds", tuple(round(float(v), 4) for v in swapped.total_bounds))

vector bounds (53.6153, -2.1733, 53.8256, -1.7272)

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

fixed = swapped.set_geometry(swapped.geometry.apply(lambda g: transform(lambda x, y, z=None: (y, x), g)))
print("fixed bounds ", tuple(round(float(v), 4) for v in fixed.total_bounds))

print("matches the original:", fixed.geometry.iloc[0].equals_exact(boundary.geometry.iloc[0], 1e-9))

fixed bounds  (-2.1733, 53.6153, -1.7272, 53.8256)
matches the original: 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

scan.tif is a piece of the elevation raster written from a bare NumPy array, with no profile — the same state as a PNG, a JPEG or a scanned map:

with rasterio.open("scan.tif") as scan:
    print("bounds", tuple(scan.bounds), "res", scan.res, "crs", scan.crs)
    print(scan.transform)

bounds (0.0, 800.0, 1000.0, 0.0) res (1.0, 1.0) crs None
| 1.00, 0.00, 0.00|
| 0.00, 1.00, 0.00|
| 0.00, 0.00, 1.00|

Bounds that are exactly (0, height, width, 0) — pixel counts — with 1.0 resolution mean the identity transform: no georeferencing at all. Note the positive y-scale in the transform, 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(388000, 437000, 30, 30)          # top-left x, top-left y, xres, yres
# or, if you know the full extent:
transform = from_bounds(388000, 413000, 418000, 437000, width=1000, height=800)

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. The Copernicus tiles themselves show the pattern — each one is named for a whole degree:

with rasterio.open("tiles/Copernicus_DSM_COG_10_N53_00_W002_00_DEM.tif") as tile:
    print("tile bounds ", tuple(round(v, 6) for v in tile.bounds))
    print("half a cell ", round(tile.res[0] / 2, 6), "×", round(tile.res[1] / 2, 6))

    print("centre of the top-left cell", tuple(round(float(v), 6) for v in tile.xy(0, 0)))

tile bounds  (-2.000208, 53.000139, -1.000208, 54.000139)
half a cell  0.000208 × 0.000139
centre of the top-left cell (-2.0, 54.0)

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

A difference of exactly res / 2 is diagnostic: cell edges and cell centres have been mixed. These tiles put a cell centre on each whole degree, so their edges — which is what bounds reports — sit half a cell outside it. 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:

from rasterio.warp import calculate_default_transform

for name in ("N53_00_W003_00", "N53_00_W002_00"):
    with rasterio.open(f"tiles/Copernicus_DSM_COG_10_{name}_DEM.tif") as tile:
        t, w, h = calculate_default_transform(tile.crs, 27700, tile.width, tile.height, *tile.bounds, resolution=30)
    print(f"{name} reprojected alone: origin ({t.c:,.2f}, {t.f:,.2f}), offset from a 30 m grid "
          f"({t.c % 30:.2f}, {t.f % 30:.2f}) m")

N53_00_W003_00 reprojected alone: origin (332,972.05, 456,475.20), offset from a 30 m grid (2.05, 25.20) m
N53_00_W002_00 reprojected alone: origin (400,083.25, 456,474.33), offset from a 30 m grid (3.25, 24.33) m

Two neighbouring tiles, the same CRS, the same 30 m cells — and grids that disagree by a fraction of a cell, so every comparison between them resamples. 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)

One function that diagnoses all five causes

The data used on this page

elevation.tif — the Copernicus 30 m digital elevation model (GLO-30), two 1° tiles covering 53–54°N, 3–1°W, merged and reprojected to British National Grid at 30 m.
boundary.gpkg — the Calderdale district from the ONS, in WGS84; study_area.gpkg the same in British National Grid; boundary_no_crs.gpkg and boundary_swapped.gpkg copies with the CRS stripped and the axes swapped, to demonstrate causes 2 and 3.
gb_outline_27700.gpkg — England, Scotland and Wales from the ONS, in British National Grid.
scan.tif — 800 × 1,000 cells of the elevation raster written with no transform or CRS.
tiles/ — the two Copernicus tiles as downloaded, plus elevation.tif again, for the folder audit.

import io
import shutil
from pathlib import Path

import geopandas as gpd
import numpy as np
import rasterio
import requests
from rasterio.merge import merge
from rasterio.warp import Resampling, calculate_default_transform, reproject
from shapely.ops import transform as swap_axes

TILE = "https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_{0}_DEM/Copernicus_DSM_COG_10_{0}_DEM.tif"
Path("tiles").mkdir(exist_ok=True)
for t in ("N53_00_W003_00", "N53_00_W002_00"):
    Path(f"tiles/Copernicus_DSM_COG_10_{t}_DEM.tif").write_bytes(requests.get(TILE.format(t), timeout=600).content)
tiles = [rasterio.open(p) for p in sorted(Path("tiles").glob("*.tif"))]
mosaic, mosaic_t = merge(tiles)
t, w, h = calculate_default_transform(tiles[0].crs, "EPSG:27700", mosaic.shape[2], mosaic.shape[1],
                                      *rasterio.transform.array_bounds(mosaic.shape[1], mosaic.shape[2], mosaic_t),
                                      resolution=30)
dem = np.full((h, w), -9999, dtype="float32")
reproject(mosaic[0], dem, src_transform=mosaic_t, src_crs=tiles[0].crs, dst_transform=t, dst_crs="EPSG:27700",
          resampling=Resampling.bilinear, dst_nodata=-9999)
with rasterio.open("elevation.tif", "w", driver="GTiff", height=h, width=w, count=1, dtype="float32",
                   crs="EPSG:27700", transform=t, nodata=-9999) as dst:
    dst.write(dem, 1)
with rasterio.open("scan.tif", "w", driver="GTiff", height=800, width=1000, count=1, dtype="float32") as dst:
    dst.write(dem[1500:2300, 1800:2800], 1)                  # no transform, no CRS
shutil.copy("elevation.tif", "tiles/elevation_bng.tif")

ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"


def ons(layer, where, fields, sr):
    r = requests.get(f"{ONS}/{layer}/FeatureServer/0/query", timeout=120,
                     params={"where": where, "outFields": fields, "outSR": sr, "f": "geojson"})
    return gpd.read_file(io.BytesIO(r.content)).set_crs(sr, allow_override=True)


calderdale = ons("Local_Authority_Districts_DEC_2025_Boundaries_UK_BFC", "LAD25NM = 'Calderdale'", "LAD25NM", 4326)
calderdale.to_file("boundary.gpkg")
calderdale.to_crs(27700).to_file("study_area.gpkg")
calderdale.set_crs(None, allow_override=True).to_file("boundary_no_crs.gpkg")
calderdale.set_geometry(calderdale.geometry.apply(lambda g: swap_axes(lambda x, y, z=None: (y, x), g))).to_file(
    "boundary_swapped.gpkg")
ons("Countries_December_2024_Boundaries_UK_BUC", "CTRY24NM <> 'Northern Ireland'", "CTRY24NM", 27700).to_file(
    "gb_outline_27700.gpkg")
print(f"elevation.tif {w:,} × {h:,}; tiles/: {sorted(p.name for p in Path('tiles').glob('*.tif'))}")

elevation.tif 4,474 × 3,725; tiles/: ['Copernicus_DSM_COG_10_N53_00_W002_00_DEM.tif', 'Copernicus_DSM_COG_10_N53_00_W003_00_DEM.tif', 'elevation_bng.tif']

Example 1: one function that diagnoses all five causes

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(float(v) for v in gdf.total_bounds)

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

    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"

    def lands_on_raster(bounds):
        footprint = gpd.GeoSeries([box(*bounds)], crs=v_crs).to_crs(r_crs).iloc[0]
        return box(*r_bounds).intersects(footprint)

    if v_crs.is_geographic and not lands_on_raster(v_bounds):
        swapped = (v_bounds[1], v_bounds[0], v_bounds[3], v_bounds[2])
        if lands_on_raster(swapped):
            return "vector axes look swapped (lat in x) — swapping them lands on the raster"
    if v_crs != r_crs:
        return f"different CRS — call gdf.to_crs({r_crs})"
    if not box(*r_bounds).intersects(box(*v_bounds)):
        return "same CRS but disjoint extents — different areas entirely"
    return "extents overlap and CRS match — alignment looks fine"


for raster, vector in [("elevation.tif", "boundary.gpkg"), ("elevation.tif", "boundary_no_crs.gpkg"),
                       ("elevation.tif", "boundary_swapped.gpkg"), ("scan.tif", "study_area.gpkg"),
                       ("elevation.tif", "study_area.gpkg")]:
    print(f"→ {diagnose(raster, vector)}\n")

raster  EPSG:27700  res (30.0, 30.0)
vector  EPSG:4326  (-2.173, 53.615, -1.727, 53.826)
→ different CRS — call gdf.to_crs(EPSG:27700)

raster  EPSG:27700  res (30.0, 30.0)
vector  None  (-2.173, 53.615, -1.727, 53.826)
→ vector has no CRS — identify it, do not guess

raster  EPSG:27700  res (30.0, 30.0)
vector  EPSG:4326  (53.615, -2.173, 53.826, -1.727)
→ vector axes look swapped (lat in x) — swapping them lands on the raster

raster  None  res (1.0, 1.0)
vector  EPSG:27700  (388662.758, 413196.876, 418113.962, 436596.385)
→ raster has the identity transform — no georeferencing

raster  EPSG:27700  res (30.0, 30.0)
vector  EPSG:27700  (388662.758, 413196.876, 418113.962, 436596.385)
→ extents overlap and CRS match — alignment looks fine

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. The swapped-axes test comes before the CRS comparison for the same reason: a swapped layer is still labelled EPSG:4326, so "reproject it" would be exactly the wrong advice.

The verdicts example 1 printed for the five pairs. The swapped copy is labelled EPSG:4326 like the correct one — only the axis test tells them apart.

Example 2: a plot that makes the misalignment obvious

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

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


ax = overlay("elevation.tif", gpd.read_file("boundary.gpkg"))
ax.figure.savefig("overlay.png", dpi=70)

reprojecting vector EPSG:4326 → EPSG:27700

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 hundreds of kilometres 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 halfway through:

from collections import Counter


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:
            identity = src.transform.a == 1 and src.transform.e == 1
            overlaps = (False if src.crs is None or identity
                        else not rasterio.coords.disjoint_bounds(src.bounds, gdf.to_crs(src.crs).total_bounds))
            rows.append({"file": p.name, "crs": str(src.crs), "res": tuple(round(r, 6) for r in src.res),
                         "same_crs": src.crs == gdf.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["same_crs"] or not r["overlaps"]:
            reason = ("no georeferencing" if r["identity"] else "no CRS" if r["crs"] == "None"
                      else "disjoint extent" if not r["overlaps"] else "different CRS — reproject before stacking")
            print(f"  ✗ {r['file']:<40} {reason}")
    return rows


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

vector CRS: EPSG:27700
raster CRS: {'EPSG:4326': 2, 'EPSG:27700': 1}
resolutions: [(0.000417, 0.000278), (30.0, 30.0)]
  ✗ Copernicus_DSM_COG_10_N53_00_W002_00_DEM.tif different CRS — reproject before stacking
  ✗ Copernicus_DSM_COG_10_N53_00_W003_00_DEM.tif different CRS — reproject before stacking

Two of the three rasters are in a different CRS from the study area and would need reprojecting before any cell-by-cell work with it — and with each other only after they share a grid, as step 5 showed. Catching it in the audit costs seconds, against a batch job that runs for hours and then reports mystery failures. This is the same pre-flight discipline as validating pipeline inputs automatically.

Do the datasets agree on what numbers mean

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.

Each of the five causes is one of those parts missing or wrong, and the steps above read each from the numbers. One case they cannot read is two projected systems with similar ranges: UTM zones 30N and 31N produce eastings in the same numeric range, so data labelled with the wrong zone lands a few hundred kilometres off with nothing in the bounds to give it away.

Equivalent CRSs that compare unequal

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 couple of odd files in a set of hundreds is exactly the failure that surfaces hours into a batch run.