Fixing Missing .prj Sidecar Files #

A shapefile lands in the inbox with its .shp, .shx, and .dbf present but no .prj; a GeoTIFF has no CRS tag; an XYZ text dump has never carried a projection in its life. Open any of them and the tooling either reports crs=None or silently assumes a default, and from that moment every overlay, buffer, and reprojection mis-registers. This is one of the most common intake failures in Geodetic Datums & EPSG Management, and the fix is subtle: you must assign the CRS the coordinates are already in, never reproject them, and you must refuse to guess when the evidence is ambiguous. This page walks the full repair — inference, assignment, and writing a valid ESRI WKT sidecar — with production code that raises rather than corrupt data silently.

Root cause: where the CRS lives, and why “missing” means “undefined” #

A shapefile is not one file; it is a set. The geometry sits in .shp, the index in .shx, attributes in .dbf, and the coordinate reference system lives entirely in the .prj sidecar as an ESRI-flavour WKT string. Lose the .prj and the coordinates are just unlabelled numbers — the format has nowhere else to record what they mean. GeoTIFF can embed CRS in its GeoKeys, but a careless export can omit them, and a plain XYZ text file has no place for a CRS at all. In every case the correct interpretation is not “assume WGS84”; it is undefined, and any tool that silently supplies a default is manufacturing information that was never there.

The failure reproduces immediately:

import geopandas as gpd

gdf = gpd.read_file("data/jurisdiction_no_prj.shp")
print(gdf.crs)     # -> None
print(gdf.total_bounds)
# -> [ 412000. 4516000.  468000. 4562000.]   <-- metres: clearly a UTM grid, not degrees

The bounds betray the truth the missing .prj withheld: six-figure coordinates are metres in a projected grid, not the degrees a EPSG:4326 default would imply. Assigning WGS84 here — as some tools do silently — would place the data thousands of kilometres away. This is a cousin of the axis-order trap covered in Debugging Axis-Order Inversion in pyproj: the numbers are fine, the interpretation is missing, and once the wrong interpretation is baked in, downstream tools cannot tell. Resolving it correctly is exactly the diagnostic discipline that also drives Fixing CRS Mismatch in QGIS and GeoPandas.

Step-by-step fix with production code #

1. Infer the true CRS from coordinate ranges — and refuse to guess #

Use the coordinate magnitudes plus the known survey region to narrow the candidate. Degrees fall in $[-180, 180]$ for longitude and $[-90, 90]$ for latitude; projected UTM eastings run roughly $100{,}000$–$900{,}000$ m. That distinction is reliable; picking the exact UTM zone or datum from ranges alone is not, so the inference raises when the evidence is ambiguous instead of silently choosing.

import logging

import numpy as np
import pyproj

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s | %(levelname)s | %(message)s",
    datefmt="%Y-%m-%dT%H:%M:%S",
)
logger = logging.getLogger(__name__)


def infer_crs_candidate(
    minx: float,
    miny: float,
    maxx: float,
    maxy: float,
    *,
    expected_utm_epsg: int | None = None,
) -> pyproj.CRS:
    """Infer the CRS class from coordinate ranges; require a hint when projected.

    Degrees resolve to EPSG:4326. Projected metres require an explicit
    `expected_utm_epsg` from survey metadata — the ranges alone cannot
    distinguish adjacent zones or datums, so we raise rather than guess.
    """
    looks_geographic = (
        -180.0 <= minx <= 180.0 and -180.0 <= maxx <= 180.0
        and -90.0 <= miny <= 90.0 and -90.0 <= maxy <= 90.0
    )
    if looks_geographic:
        logger.info("Ranges are within degree bounds; inferring EPSG:4326")
        return pyproj.CRS.from_epsg(4326)

    looks_projected_m = (10_000 <= abs(maxx) and 10_000 <= abs(maxy))
    if looks_projected_m:
        if expected_utm_epsg is None:
            raise ValueError(
                f"Coordinates {(minx, miny, maxx, maxy)} are projected metres, "
                "but the exact CRS is ambiguous. Supply expected_utm_epsg from "
                "the survey metadata; never guess the zone from ranges alone."
            )
        logger.info("Assigning caller-supplied projected CRS EPSG:%d", expected_utm_epsg)
        return pyproj.CRS.from_epsg(expected_utm_epsg)

    raise ValueError(f"Ranges {(minx, miny, maxx, maxy)} match no known CRS class")

2. Assign the CRS in place — never reproject #

set_crs() (vector) and write_crs() (raster) label coordinates that are already in that system. Calling to_crs() here would move the geometry as though it were already correctly tagged, doubling the error. The distinction is the whole game.

import geopandas as gpd


def assign_vector_crs(path: str, expected_utm_epsg: int | None = None) -> gpd.GeoDataFrame:
    """Attach the inferred CRS to an unlabelled vector layer without moving it."""
    gdf = gpd.read_file(path, engine="pyogrio")
    if gdf.crs is not None:
        logger.info("Layer already has CRS %s; nothing to assign", gdf.crs.to_string())
        return gdf

    minx, miny, maxx, maxy = gdf.total_bounds
    crs = infer_crs_candidate(minx, miny, maxx, maxy, expected_utm_epsg=expected_utm_epsg)
    # allow_override=False guards against silently relabelling a set CRS.
    gdf = gdf.set_crs(crs, allow_override=False)
    logger.info("Assigned CRS EPSG:%s to %s", crs.to_epsg(), path)
    return gdf

3. Write a valid .prj in ESRI WKT #

Shapefile consumers — especially ArcGIS — expect the sidecar in ESRI WKT (WKT1 ESRI variant), not OGC WKT2. pyproj serialises either; use the ESRI dialect for .prj so older tooling parses it.

import pathlib


def write_prj_sidecar(shp_path: str, crs: pyproj.CRS) -> pathlib.Path:
    """Write a .prj sidecar in ESRI WKT next to a shapefile."""
    prj_path = pathlib.Path(shp_path).with_suffix(".prj")
    esri_wkt = crs.to_wkt(version="WKT1_ESRI")
    if not esri_wkt:
        raise ValueError(f"Could not serialise ESRI WKT for CRS {crs.to_epsg()}")
    prj_path.write_text(esri_wkt, encoding="utf-8")
    logger.info("Wrote %s (EPSG:%s, ESRI WKT)", prj_path, crs.to_epsg())
    return prj_path

For a GeoTIFF missing its CRS tag, set it in place with rasterio — again assigning, not warping:

import rasterio


def assign_raster_crs(path: str, epsg: int) -> None:
    """Write a CRS into a GeoTIFF's header without resampling pixels."""
    with rasterio.open(path, "r+") as src:
        if src.crs is not None:
            logger.info("%s already has CRS %s; leaving untouched", path, src.crs)
            return
        src.crs = rasterio.crs.CRS.from_epsg(epsg)
        logger.info("Stamped EPSG:%d onto %s (no pixel resampling)", epsg, path)

4. Document provenance #

Record how the CRS was determined — survey metadata, provider email, range inference plus a named analyst — in a sidecar log or the dataset’s metadata. An assigned CRS with no provenance is a future mystery; the parent Geodetic Datums & EPSG Management workflow treats this audit trail as mandatory.

Verification and acceptance test #

Reopen the repaired file, assert the CRS is now set, and confirm the coordinates fall inside the expected region for that CRS. A .prj that writes cleanly but places the data in the wrong hemisphere is not a fix.

import pytest


REGION_BOUNDS_UTM18N = (400_000, 4_200_000, 560_000, 4_600_000)  # Delaware area, metres


def test_repaired_layer_has_crs_and_is_in_region():
    gdf = assign_vector_crs("data/jurisdiction_no_prj.shp", expected_utm_epsg=32618)
    write_prj_sidecar("data/jurisdiction_no_prj.shp", gdf.crs)

    reopened = gpd.read_file("data/jurisdiction_no_prj.shp")
    assert reopened.crs is not None
    assert reopened.crs.to_epsg() == 32618

    minx, miny, maxx, maxy = reopened.total_bounds
    rx0, ry0, rx1, ry1 = REGION_BOUNDS_UTM18N
    assert rx0 <= minx and maxx <= rx1, f"x bounds {minx}-{maxx} outside region"
    assert ry0 <= miny and maxy <= ry1, f"y bounds {miny}-{maxy} outside region"


def test_ambiguous_projected_data_raises():
    # Projected metres with no metadata hint must not be silently assigned.
    with pytest.raises(ValueError, match="ambiguous"):
        infer_crs_candidate(412000, 4516000, 468000, 4562000)

A quick CLI cross-check confirms GDAL now reads the CRS straight from the file:

gdalsrsinfo -o epsg data/jurisdiction_no_prj.shp   # -> EPSG:32618

Edge cases and gotchas #

• ESRI WKT vs OGC WKT2. ArcGIS-era tooling expects the ESRI WKT1 dialect in .prj; handing it OGC WKT2 can produce an “unknown CRS” in some readers. Serialise with to_wkt(version="WKT1_ESRI") for the sidecar, but store WKT2 or the EPSG code in your own metadata where round-tripping fidelity matters.
• .prj present but wrong. A sidecar that exists but declares the wrong system is more dangerous than a missing one, because tools trust it. Run the coordinate-range sanity check even when a .prj is present, and override with allow_override=True only after confirming the declared tag is wrong.
• EPSG code vs proj string. Prefer assigning an EPSG code over a bare PROJ.4 string: proj strings lose the datum’s precise realisation and can silently drop the towgs84 shift, reintroducing a metre-scale error. If a source only gives you a proj string, resolve it to the nearest authoritative EPSG code before writing the sidecar.
• Ranges ambiguous between zones. Coordinates near a UTM zone boundary, or a survey that straddles two zones, can look valid in either. The ranges cannot disambiguate — this is precisely when infer_crs_candidate raises and demands the zone from survey metadata rather than picking one.

Related #

• Debugging Axis-Order Inversion in pyproj — the sibling case where the CRS is present but its lat/lon order is swapped.
• Fixing CRS Mismatch in QGIS and GeoPandas — diagnosing the on-the-fly rendering that hides a missing or wrong CRS.
• Geodetic Datums & EPSG Management — the parent workflow: datum identity, EPSG resolution, and provenance requirements.

Up: Geodetic Datums & EPSG Management