Fixing CRS Mismatches in Watershed Shapefiles

Watershed boundaries distributed by agencies often arrive in a different coordinate reference system than your DEM — WGS84 geographic coordinates against a UTM-projected raster, or NAD27 against NAD83. As part of coordinate reference system alignment — and the broader Hydrology Data Preparation & DEM Processing pipeline — permanently reprojecting shapefiles to match the DEM’s CRS is the single highest-leverage step you can take before any geometric operation. This page covers the exact Python procedure: detecting the mismatch, repairing geometries, reprojecting with an explicit EPSG code, validating alignment against the raster, and exporting a corrected file that downstream tools can use without further guesswork.

Prerequisites

This page assumes you already have a working Python environment meeting the baseline described in the coordinate reference system alignment overview. Beyond that baseline, this technique requires:

• geopandas >= 0.13 (for stable to_crs() datum-grid support via pyproj >= 3.4)
• shapely >= 2.0 (for make_valid() as a module-level function rather than a method)
• proj-data package installed via conda-forge (required for NAD27→NAD83 grid shifts and other datum transformations)
• A raster DEM in a projected CRS — you need its EPSG code as the reprojection target

conda install -c conda-forge geopandas shapely rasterio pyproj proj-data

If proj-data is absent, pyproj silently falls back to approximate transformations for datum shifts, introducing metre-scale errors that are difficult to detect downstream.

Core Technique: Permanent Reprojection with to_crs()

The central mechanism is geopandas.GeoDataFrame.to_crs(), which applies a pyproj coordinate transform to every vertex of every geometry in the frame and returns a new GeoDataFrame carrying the target CRS. Critically, this only produces a permanently corrected file when you subsequently write the result with to_file() — the .shp and .prj on disk remain unchanged until you do.

Three details determine whether the transformation is geometrically correct:

1. Source CRS must be explicitly known. If gdf.crs is None, to_crs() raises a ValueError. You must call set_crs() with the correct EPSG before reprojecting.
2. Geometry validity must be enforced before transformation. Projection math can produce NaN coordinates or collapsed rings when applied to self-intersecting polygons — geometries that are technically invalid under the OGC spec. shapely.make_valid() resolves these first.
3. Axis order follows the EPSG authority by default. With EPSG codes, geopandas delegates axis order to pyproj, which respects the authority definition (lat/lon for geographic, easting/northing for most projected CRS). When using raw PROJ strings, pass always_xy=True to force (x=easting, y=northing) regardless of authority.

Annotated Code Example

import geopandas as gpd
from pyproj import CRS
from shapely.validation import make_valid  # requires shapely >= 2.0
import rasterio
import logging
import os

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")


def fix_watershed_crs(
    shapefile_path: str,
    target_epsg: int,
    output_path: str,
    alignment_tolerance_m: float = 100.0,
    dem_path: str | None = None,
) -> gpd.GeoDataFrame:
    """
    Detect and permanently reproject a watershed shapefile to a target CRS.

    Parameters
    ----------
    shapefile_path      : Path to the input .shp file.
    target_epsg         : EPSG code of the desired output CRS (e.g. 32614 for UTM 14N).
    output_path         : Where to write the reprojected shapefile.
    alignment_tolerance_m : Maximum permitted bounding-box offset in metres when
                           cross-checking against a DEM. Ignored if dem_path is None.
    dem_path            : Optional path to the reference DEM for alignment validation.
    """
    if not os.path.exists(shapefile_path):
        raise FileNotFoundError(f"Shapefile not found: {shapefile_path}")

    # --- 1. Load with explicit UTF-8 encoding to prevent attribute-table corruption ---
    # Agency shapefiles sometimes ship with Latin-1 encoded .dbf files; UTF-8 is the
    # safest default — if characters are garbled, retry with encoding="latin-1".
    gdf = gpd.read_file(shapefile_path, encoding="utf-8")
    logging.info("Loaded %d features from %s", len(gdf), shapefile_path)

    # --- 2. Guard against undefined CRS (missing or empty .prj) ---
    # A None CRS means geopandas cannot determine the projection. You must call
    # gdf = gdf.set_crs("EPSG:XXXX") before reprojecting. Stop here so the caller
    # knows they need to supply the correct EPSG rather than silently continuing.
    if gdf.crs is None:
        raise ValueError(
            "No CRS found in the .prj file. Assign the source CRS with "
            "gdf.set_crs('EPSG:XXXX', inplace=True) before calling this function."
        )
    logging.info("Source CRS: %s (EPSG:%s)", gdf.crs.name, gdf.crs.to_epsg())

    # --- 3. Repair invalid geometries before any coordinate transformation ---
    # Watershed polygon boundaries from older agency datasets often have
    # self-intersections or bowtie artefacts. make_valid() uses the GEOS
    # BufferWithStyle algorithm to fix these without discarding vertices.
    invalid_mask = ~gdf.is_valid
    if invalid_mask.any():
        n_invalid = int(invalid_mask.sum())
        logging.warning("Repairing %d invalid geometr%s...", n_invalid, "ies" if n_invalid > 1 else "y")
        gdf.loc[invalid_mask, "geometry"] = (
            gdf.loc[invalid_mask, "geometry"].apply(make_valid)
        )
        # Confirm repair succeeded — any remaining invalids abort the workflow
        still_invalid = (~gdf.is_valid).sum()
        if still_invalid:
            raise RuntimeError(f"{still_invalid} geometr(ies) could not be repaired.")

    # --- 4. Build the target CRS object and compare to source ---
    target_crs = CRS.from_epsg(target_epsg)
    if gdf.crs == target_crs:
        logging.info("Source and target CRS already match. Skipping transformation.")
    else:
        logging.info(
            "Reprojecting: %s → EPSG:%d (%s)",
            gdf.crs.to_string(), target_epsg, target_crs.name,
        )
        # to_crs() applies pyproj's full transformation pipeline, including
        # grid-shift datum changes (NAD27→NAD83, etc.) when proj-data is installed.
        gdf = gdf.to_crs(target_crs)

    # --- 5. Sanity-check projected bounds ---
    bounds = gdf.total_bounds  # (minx, miny, maxx, maxy)
    if target_crs.is_geographic:
        # Geographic CRS: coordinates should stay within ±180 / ±90
        if not (-180 <= bounds[0] <= bounds[2] <= 180 and -90 <= bounds[1] <= bounds[3] <= 90):
            logging.warning(
                "Geographic bounds look suspicious: %s — confirm target EPSG.", bounds
            )
    else:
        # Projected CRS: coordinates should NOT be in the degree range (-180 to 180)
        if -180 <= bounds[0] <= 180 and -90 <= bounds[1] <= 90:
            logging.warning(
                "Projected bounds look degree-scaled (%s) — source CRS may be mislabelled.", bounds
            )

    # --- 6. Optional: cross-check alignment against the reference DEM ---
    if dem_path is not None:
        with rasterio.open(dem_path) as src:
            if str(src.crs.to_epsg()) != str(target_epsg):
                logging.warning(
                    "DEM CRS (EPSG:%s) differs from target EPSG:%d",
                    src.crs.to_epsg(), target_epsg,
                )
            dem_b = src.bounds
            x_offset = abs(dem_b.left - bounds[0])
            y_offset = abs(dem_b.bottom - bounds[1])
            logging.info(
                "Alignment check — X offset: %.1f m, Y offset: %.1f m", x_offset, y_offset
            )
            if x_offset > alignment_tolerance_m or y_offset > alignment_tolerance_m:
                raise ValueError(
                    f"Watershed/DEM origin offset ({x_offset:.1f} m, {y_offset:.1f} m) "
                    f"exceeds tolerance ({alignment_tolerance_m} m). "
                    "Verify EPSG codes for both datasets."
                )

    # --- 7. Export — writes updated coordinates and a fresh .prj sidecar ---
    gdf.to_file(output_path, driver="ESRI Shapefile")
    logging.info("Reprojected watershed saved to %s", output_path)
    return gdf

Parameter Reference

Worked Example: Identifying and Fixing a Degree/Metre Mismatch

A common scenario: you clip a 10 m DEM to a watershed boundary and the output is empty, or you receive an obvious-looking clip that places the catchment polygon in the ocean. Running the function with a mismatched pair produces diagnostic log output you can read immediately:

INFO: Loaded 1 features from NHD_catchment_09010001.shp
INFO: Source CRS: WGS 84 (EPSG:4326)
INFO: Reprojecting: EPSG:4326+4979 → EPSG:32614 (WGS 84 / UTM zone 14N)
INFO: Alignment check — X offset: 2.4 m, Y offset: 1.8 m
INFO: Reprojected watershed saved to NHD_catchment_09010001_utm14n.shp

An X offset of 2.4 m is sub-pixel for a 10 m DEM — acceptable. If you see offsets in the thousands-of-metres range, the source EPSG in the .prj file is almost certainly wrong. Cross-check against the agency metadata (FGDC or ISO 19115 XML) or the original download portal.

After reprojection, verify area conservation: compute gdf.geometry.area (in square metres for projected CRS) both before and after reprojection. For catchments originally in geographic CRS, a 10–30 % area discrepancy in degree units relative to projected-CRS metres is expected and correct — if the areas match exactly, the transformation may have been skipped silently.

Before running any DEM pit-filling algorithms or clipping operations, this validation ensures the catchment boundary and raster share the same spatial frame of reference.

Gotchas and Edge Cases

• .prj file present but incorrect. USGS National Hydrography Dataset shapefiles and some NWM boundary products historically shipped with placeholder or incomplete .prj tags. gdf.crs.to_epsg() returns a code that looks right but is off by one UTM zone or uses the wrong datum. Always cross-reference with the download metadata before assuming the .prj is authoritative.

• NAD27 datum shifts require proj-data. If you are working with pre-1980 USGS quadrangle boundaries, the source datum may be NAD27. Without the proj-data conda package, pyproj applies a simplified approximate conversion that can be off by 20–100 m in the contiguous US — large enough to shift stream channels off-network after spatial resolution resampling.

• Multi-polygon features from watershed union operations. If the shapefile was generated by dissolving sub-basins, some features may be MULTIPOLYGON or GEOMETRYCOLLECTION objects. make_valid() does not change geometry type, but to_crs() handles all geometry types natively. Only operations that assume POLYGON inputs (some older richdem wrappers) will fail — inspect with gdf.geom_type.value_counts() first.

• Silent success on already-projected files. If the source EPSG matches the target EPSG exactly, to_crs() returns the original frame unchanged with no log message unless you add the equality check shown in the code above. Without that check, a pipeline that runs twice on the same file will not report an error — but you also get no confirmation that the transform ran, making reproducibility audits harder.

• Vertical datum is independent. Fixing the horizontal CRS does not touch the elevation values in a DEM. If your hydrological model requires NAVD88 orthometric heights and your LiDAR DEM is in NAVD88 but your boundary comes from a GCS-referenced dataset, the horizontal alignment fix here is complete. If the DEM is in ellipsoidal heights (common from some Copernicus products), a separate geoid conversion step is required before computing slope or running flow-direction algorithms — this is outside the scope of shapefile reprojection.

FAQ

Why does on-the-fly projection in QGIS not fix my shapefile’s CRS?

Desktop GIS tools reproject geometry on screen without writing new coordinates to disk. The .shp and .prj files remain unchanged, so Python pipelines that read those files directly still encounter the original CRS. On-the-fly reprojection is useful for visual inspection but must never be treated as a substitute for a permanent to_crs() + to_file() write cycle.

What happens if I run DEM pit-filling on a mismatched shapefile boundary?

The clipping mask will be misaligned by tens to hundreds of metres, leaving real depressions inside the boundary unfilled and potentially filling valid terrain outside it. This corrupts all downstream flow accumulation and watershed delineation. Completing DEM pit-filling algorithms before CRS alignment is one of the most common ordering mistakes in production pipelines.

How do I handle a shapefile with no .prj file?

geopandas sets gdf.crs to None. You must assign the correct CRS manually with gdf = gdf.set_crs("EPSG:XXXX") before calling to_crs(). Consult the data provider’s metadata — FGDC XML, ISO 19115 records, or the download portal’s layer properties — to identify the original projection. If metadata is absent, check whether coordinate values are in degree range (±180 / ±90) for geographic CRS, or large integer values for projected CRS.

Related Topics

• Coordinate Reference System Alignment — parent overview covering CRS detection, datum handling, and raster/vector alignment strategies
• Hydrology Data Preparation & DEM Processing — full DEM conditioning pipeline from raw data acquisition through flow analysis
• Resampling DEMs Without Losing Hydrologic Connectivity — the next conditioning step after CRS alignment is confirmed
• DEM Pit Filling Algorithms — sink removal that must follow CRS alignment to avoid clipping artefacts
• Best Practices for Filling Sinks in High-Resolution LiDAR Data — complementary LiDAR pit-fill workflow after boundary reprojection is complete