Coastal & Marine Spatial Analysis Pipelines

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.

Assign-Don't-Guess Flow for a Missing CRS A dataset with no CRS enters at the left. First decision: is crs None? If no, the file passes. If yes, inspect coordinate ranges. A branch shows two range signatures — values around tens of degrees imply a geographic CRS such as EPSG:4326, values around hundreds of thousands imply a projected UTM grid. If the ranges resolve to a single unambiguous candidate, the flow proceeds to assign the CRS in place (no reprojection) and write a valid ESRI WKT prj sidecar, ending in a verified repaired file. If the ranges are ambiguous between adjacent zones or datums, the flow raises an error rather than guessing. read file .shp/.tif/.xyz crs None? no pass through yes inspect coord ranges deg vs metres one candidate? no ✗ raise ambiguous — never guess yes assign (no reproject) write .prj ESRI WKT verify

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][-180, 180] for longitude and [−90,90][-90, 90] for latitude; projected UTM eastings run roughly 100,000100{,}000–900,000900{,}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")
Where the CRS Actually Lives, Per Format Four format panels compared on one axis: whether the coordinate reference system is stored inside the data file or beside it. The shapefile panel shows the multi-file set .shp, .shx, .dbf and .prj, with the .prj highlighted as a detachable plain-text sidecar that an ordinary copy or upload can leave behind. The GeoPackage panel shows a single container file with an internal gpkg_spatial_ref_sys table. The GeoTIFF panel shows a single file with embedded GeoKey tags. The CSV panel shows columns of numbers with no place to record a CRS at all, dependent on an external agreement. A caption notes that only the first can lose its CRS by accident, and only the last never had one. Shapefile survey.shp — geometry survey.shx — index survey.dbf — attributes survey.prj — the CRS detachable · a copy can lose it GeoPackage survey.gpkg features table gpkg_spatial_ref_sys one file, one copy inseparable · travels with the data GeoTIFF bathy.tif GeoKey tags pixel grid + geotransform embedded · header-resident Plain CSV x, y, depth 454102.3, 4293201.2, -12.4 454103.3, 4293201.9, -12.5 nowhere to put a CRS external convention only Only the shapefile can lose a CRS it once had, and only the CSV never had one — the two cases look identical to a reader and need different repairs. Inferring from coordinate ranges tells you which case you are in. The repair below writes a .prj so the shapefile joins the second column: self-describing, not dependent on a note in someone’s email.

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.

What Coordinate Magnitudes Can and Cannot Tell You Three horizontal magnitude bands drawn on a shared logarithmic-feeling scale of absolute coordinate value. The first band, roughly 0 to 180 with decimals, is labelled geographic degrees, EPSG:4326 or similar. The second band, roughly 100,000 to 900,000 easting paired with 0 to 10,000,000 northing, is labelled UTM metres, with a note that the easting range is the same in every one of the sixty zones. The third band, in the millions of survey feet, is labelled State Plane. A fourth row highlights the overlap region where a value could belong to more than one system, marked "ambiguous — do not guess", with a pointer to reading the survey report or the adjacent tiles instead. 10⁰ 10² 10⁴ 10⁶ 10⁸ absolute value of the stored coordinate geographic degrees |lon| ≤ 180, |lat| ≤ 90 EPSG:4326 · 4269 · fractional UTM metres E 100k–900k · N 0–10M zone unknown from range alone State Plane feet millions of survey feet overlap Magnitude narrows the field; it never names the system. A UTM easting looks the same in all sixty zones: the range proves the projection family, never the zone or the datum. Read those from the survey report — and refuse to write a .prj when the evidence is ambiguous.

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.

Up: Geodetic Datums & EPSG Management