⚙️ How-To CRS & projections Files & formats Notebook By · Updated

How to Work with a Custom or Local CRS in Python

On this page

Problem statement

Not every coordinate system has an EPSG code.

A site survey arrives with coordinates measured from a peg the surveyor knocked in at one corner of the site. A historical dataset uses a projection defined by a PROJ string nobody recognises. A client's engineering drawing is on a local grid with an arbitrary origin. A continental analysis needs a projection centred on your study area rather than on someone else's.

import geopandas as gpd

gdf = gpd.read_file("site_survey.shp")
print(gdf.crs)
print(gdf.total_bounds.round(2).tolist())
None
[23.4, -4.13, 1352.55, 1287.3]

Metres, clearly. From where? set_crs(27700) would place the site off the Isles of Scilly. There is no EPSG code for "1,000 m east of a peg the surveyor put in".

Meanwhile a projection you need may exist and have no code at all:

gdf.to_crs("+proj=aea +lat_1=40 +lat_2=65 +lat_0=53 +lon_0=10 +datum=WGS84 +units=m")

That works, has no EPSG number, and cannot be written to every file format. Custom CRS are entirely usable in Python — the practical issues are how to define them, how to make them travel with the data, and when to stop and get the data georeferenced properly.

Quick answer

Three situations that all present as "no EPSG code". three different problems, three answersis it a CRS?what you needa custom projection✓yesa PROJ string or WKTa local grid, tied!not yetcontrol pointsa local grid, untied✕noask the supplierNever infer a CRS from the coordinate magnitudes spatialworkflow.io
Three different situations that all present as "no EPSG code".
Situation What you need
a projection with no EPSG code a PROJ string or WKT — fully usable
a local grid with a known tie point an affine transform to a real CRS
a local grid with no tie point control points, or the data stays local
a CRS you cannot identify do not guess — ask the supplier 
from pyproj import CRS

# a custom projection, centred where you need it
crs = CRS.from_proj4(
    "+proj=aea +lat_1=40 +lat_2=65 +lat_0=53 +lon_0=10 "
    "+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
gdf = gdf.to_crs(crs)

# a local engineering grid, tied to a known point
from shapely.affinity import affine_transform
gdf["geometry"] = gdf.geometry.apply(
    lambda g: affine_transform(g, [1, 0, 0, 1, 383_618.5, 398_050.4]))
gdf = gdf.set_crs(27700)

The distinction that matters: a custom projection is a real CRS with a definition; a local grid is coordinates with no relationship to the earth until you supply one.

What you get — from 1. Define a custom projection with PROJ or WKT

unknown
True
metre
None — no EPSG code, as expected

PROJ strings, data-fitted projections, engineering grids

1. Define a custom projection with PROJ or WKT

from pyproj import CRS

# an equal-area conic centred on Europe
europe = CRS.from_proj4(
    "+proj=aea +lat_1=40 +lat_2=65 +lat_0=53 +lon_0=10 "
    "+x_0=0 +y_0=0 +ellps=GRS80 +datum=WGS84 +units=m +no_defs")

print(europe.name)
print(europe.is_projected)
print(europe.axis_info[0].unit_name)
print(europe.to_epsg(), "— no EPSG code, as expected")
unknown
True
metre
None — no EPSG code, as expected

The common projection parameters:

Parameter Meaning
+proj= the projection: aea, lcc, tmerc, utm, laea, aeqd, merc
+lat_0=, +lon_0= the origin latitude and central meridian
+lat_1=, +lat_2= standard parallels, for conic projections
+x_0=, +y_0= false easting and northing, to keep coordinates positive
+k_0= scale factor at the central meridian
+datum= or +ellps= the reference ellipsoid
+units= m, ft, us-ft
+towgs84= a crude 3- or 7-parameter datum shift

WKT2 is the modern, more expressive form, and it is what recent formats store:

print(europe.to_wkt(pretty=True)[:400])
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Albers Equal Area",

Give it a name so it is legible later:

named = CRS.from_wkt(europe.to_wkt().replace('PROJCRS["unknown"',
                                             'PROJCRS["Europe AEA (project)"', 1))
print(named.name)
Europe AEA (project)

2. Build a custom projection from your data

The most common legitimate use: a projection centred on the study area, which beats any standard code for an unusual extent.

import geopandas as gpd
from pyproj import CRS

def centred_crs(gdf, kind="laea"):
    """A projection centred on this layer. kind: laea, aea, lcc, aeqd, tmerc."""
    g = gdf.to_crs(4326)
    minx, miny, maxx, maxy = g.total_bounds
    c = g.union_all().centroid
    params = [f"+proj={kind}", f"+lat_0={c.y:.6f}", f"+lon_0={c.x:.6f}"]
    if kind in ("aea", "lcc"):
        params += [f"+lat_1={miny + (maxy - miny) / 6:.6f}",
                   f"+lat_2={maxy - (maxy - miny) / 6:.6f}"]
    params += ["+x_0=0", "+y_0=0", "+datum=WGS84", "+units=m", "+no_defs"]
    return CRS.from_proj4(" ".join(params))

study_area = gpd.read_file("study_area.gpkg")        # Leeds' 33 wards
crs = centred_crs(study_area, kind="laea")
print(crs.to_proj4())
+proj=laea +lat_0=53.825779 +lon_0=-1.506029 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs

laea (Lambert Azimuthal Equal Area) preserves area exactly around the centre point. aeqd preserves distance from the centre. aea and lcc suit wide east–west extents. Which to pick follows from what the analysis depends on, as in how to choose the right projected CRS.

Verify it against geodesic truth rather than trusting it:

import numpy as np
from pyproj import Geod

geod = Geod(ellps="WGS84")
sample = study_area.to_crs(4326)
truth = np.array([abs(geod.geometry_area_perimeter(g)[0]) for g in sample.geometry])
got = sample.to_crs(crs).geometry.area.to_numpy()
print(f"worst area error over the 33 wards: {100 * np.max(np.abs(got / truth - 1)):.4f}%")
worst area error over the 33 wards: 0.0000%

3. Handle a local engineering grid

A local site grid anchored to a national grid by a tie point. a local grid becomes a CRS only when something ties itlocal survey grid0, 01353, 1287metres from a pegtied to EPSG:27700429286, 432626plus rotation and scaleThree or more control points let you solve for it rather than assume spatialworkflow.io
A local grid becomes a CRS only when something ties it to a known location.

Site surveys often use a grid whose origin is a peg somewhere on site. Those coordinates have no relationship to the earth until you supply the tie.

With one tie point and no rotation, a translation is enough:

from shapely.affinity import affine_transform
import geopandas as gpd

# the site peg, in a real CRS — from the surveyor's notes
SITE_ORIGIN_BNG = PEG                          # EPSG:27700

def tie_local_grid(gdf, origin, target_crs, *, scale=1.0, rotation_deg=0.0):
    """Move a local grid onto a real CRS. Returns a georeferenced GeoDataFrame."""
    import math
    r = math.radians(rotation_deg)
    a =  scale * math.cos(r)
    b = -scale * math.sin(r)
    d =  scale * math.sin(r)
    e =  scale * math.cos(r)
    matrix = [a, b, d, e, origin[0], origin[1]]      # [a, b, d, e, xoff, yoff]
    out = gdf.copy()
    out["geometry"] = gdf.geometry.apply(lambda g: affine_transform(g, matrix))
    return out.set_crs(target_crs)

site = gpd.read_file("site_survey.shp")          # crs None, local metres
for rotation in (0.0, 1.7284):
    placed = tie_local_grid(site, SITE_ORIGIN_BNG, 27700, rotation_deg=rotation)
    miss = shapely.distance(placed.geometry.centroid.values, buildings.geometry.centroid.values)
    print(f"rotation {rotation:>6.4f}°   worst building {miss.max():6.2f} m from where it really is")
rotation 0.0000°   worst building  51.42 m from where it really is
rotation 1.7284°   worst building   0.01 m from where it really is

Shapely's matrix order is [a, b, d, e, xoff, yoff], giving x' = a·x + b·y + xoff and y' = d·x + e·y + yoff. It is easy to transpose by accident, so test with a known point. And note what the translation alone did: the site grid is turned 1.73° from National Grid north, and ignoring that put the far side of a site barely 1.4 km across 51 m out. The error grows with distance from the peg: 1.73° is 3 cm per metre.

With two or more tie points, solve for the transform rather than assuming it:

import numpy as np

def solve_affine(local_pts, real_pts):
    """Least-squares affine transform from >=3 matched point pairs."""
    local = np.asarray(local_pts, dtype=float)
    real = np.asarray(real_pts, dtype=float)
    if len(local) < 3:
        raise ValueError("an affine transform needs at least 3 point pairs")

    A = np.zeros((2 * len(local), 6))
    A[0::2, 0], A[0::2, 1], A[0::2, 4] = local[:, 0], local[:, 1], 1
    A[1::2, 2], A[1::2, 3], A[1::2, 5] = local[:, 0], local[:, 1], 1
    b = real.ravel()

    params, residuals, rank, _ = np.linalg.lstsq(A, b, rcond=None)
    a, bb, d, e, xoff, yoff = params

    predicted = (A @ params).reshape(-1, 2)
    err = np.linalg.norm(predicted - real, axis=1)
    print(f"{len(local)} control points, RMS residual {np.sqrt((err**2).mean()):.3f} m, "
          f"worst {err.max():.3f} m")
    return [a, bb, d, e, xoff, yoff]

matrix = solve_affine(local_pts=[c[0] for c in CONTROL], real_pts=[c[1] for c in CONTROL])
4 control points, RMS residual 0.003 m, worst 0.003 m

The RMS residual is the check. Millimetres, as here, mean the control points are consistent and the transform fits — these points were rounded to the centimetre, as a survey would record them, and that rounding is all the residual there is. A metre would mean the points are wrong, mis-matched, or the relationship is not affine.

4. Know when to stop and get it georeferenced

Sometimes the honest answer is that the data cannot be placed:

def can_georeference(gdf, control_points=None):
    if gdf.crs is not None:
        return "already has a CRS"
    if control_points and len(control_points) >= 3:
        return "affine transform from control points"
    if control_points and len(control_points) >= 1:
        return "translation only — assumes no rotation or scale"
    return ("cannot georeference: no CRS and no control points. "
            "Ask the supplier for the grid definition or a tie point.")

Do not attempt to infer a CRS from the coordinate values. Numbers in the range 0–5,000 could be metres from a site peg, metres from a national grid origin in another country, feet, or arbitrary drawing units. Assigning a CRS is a factual claim, and a wrong one produces data that is internally consistent and confidently in the wrong place — the failure described in set_crs vs to_crs.

A layer with no CRS is honest. A mislabelled one is not.

5. Make the definition travel with the data

Custom CRS support varies by format, and this is the practical constraint:

import geopandas as gpd
from pathlib import Path

def check_crs_round_trip(gdf, tmp_dir="."):
    """Which formats preserve this CRS?"""
    results = []
    for driver, suffix in [("GPKG", ".gpkg"), ("GeoJSON", ".geojson"),
                           ("ESRI Shapefile", ".shp"), ("FlatGeobuf", ".fgb"),
                           ("Parquet", ".parquet")]:
        path = Path(tmp_dir) / f"crs_test{suffix}"
        try:
            if driver == "Parquet":
                gdf.to_parquet(path)
                back = gpd.read_parquet(path)
            else:
                gdf.to_file(path, driver=driver)
                back = gpd.read_file(path)
            same = back.crs is not None and back.crs.equals(gdf.crs)
            results.append((driver, "✓ preserved" if same else
                            f"✗ became {back.crs.name if back.crs else 'None'}"))
        except Exception as exc:
            results.append((driver, f"✗ {type(exc).__name__}: {str(exc)[:44]}"))
    for driver, note in results:
        print(f"  {driver:<16} {note}")
    return results

custom = study_area.to_crs(centred_crs(study_area))
results = check_crs_round_trip(custom)
  GPKG             ✓ preserved
  GeoJSON          ✗ became WGS 84
  ESRI Shapefile   ✓ preserved
  FlatGeobuf       ✓ preserved
  Parquet          ✓ preserved

GeoJSON is the exception and it is by specification: RFC 7946 mandates WGS 84, so GDAL reprojects on write. Any custom CRS is lost — a genuine incompatibility, not a bug.

For anything else, keep the definition next to the data as well:

import json
from pathlib import Path

def save_with_crs(gdf, path, *, driver="GPKG"):
    path = Path(path)
    gdf.to_file(path, driver=driver)
    path.with_suffix(".crs.json").write_text(json.dumps({
        "wkt": gdf.crs.to_wkt(),
        "proj4": gdf.crs.to_proj4(),
        "epsg": gdf.crs.to_epsg(),
        "note": "custom projection — see project documentation",
    }, indent=2))
    return path

A sidecar is redundant when the format stores the CRS and invaluable when someone converts the file to GeoJSON and loses it.

A local CRS registry and a site survey georeferenced

Three project-local coordinate systems ordered by how far from the origin each stays usable. pick the projection by what you are measuringaeqd · true distance from origindistances to one sitetmerc · scale factor 1.0 at centrea site, under ~20 kmaea · equal area between two parallelsa continentDefine each once in a registry, with a note saying what it is for and where it stops being true spatialworkflow.io
The registry below holds all three, each with the note that says when not to use it.

The data used on this page

The site survey is real building outlines in an invented site grid, made so its truth is known. The setup takes the OpenStreetMap buildings around Leeds station from Geofabrik's West Yorkshire extract of 1 January 2026 and moves them into a local grid: origin PEG at the south-west corner, grid north turned 1.7284° from National Grid north, coordinates rounded to the centimetre. It writes that as site_survey.shp with no CRS — and GeoPandas then writes no .prj at all, which is how most CRS-less shapefiles come to exist. CONTROL is four building corners with both their site-grid and National Grid coordinates, as a surveyor would supply them. study_area.gpkg is Leeds' 33 wards from the ONS, and nuts2.gpkg Eurostat's NUTS 2 regions from GISCO. OpenStreetMap data © OpenStreetMap contributors, ODbL; NUTS © EuroGeographics.

import math
import os

import geopandas as gpd
import numpy as np
import requests
import shapely
from shapely.affinity import affine_transform

OSM = "https://download.geofabrik.de/europe/united-kingdom/england/west-yorkshire-260101-free.shp.zip"
with open("west-yorkshire.zip", "wb") as f:
    f.write(requests.get(OSM, timeout=600).content)
buildings = (gpd.read_file("zip://west-yorkshire.zip!gis_osm_buildings_a_free_1.shp",
                           bbox=(-1.556, 53.790, -1.538, 53.800))
                .to_crs(27700)[["osm_id", "geometry"]].reset_index(drop=True))

PEG = tuple(float(v) for v in np.floor(buildings.total_bounds[:2]))   # the site origin, National Grid metres
turn = math.radians(1.7284)                               # site grid north, from National Grid north
c, s = math.cos(turn), math.sin(turn)
to_site = [c, s, -s, c, -(c * PEG[0] + s * PEG[1]), -(-s * PEG[0] + c * PEG[1])]
survey = buildings.copy()
survey["geometry"] = shapely.set_precision(
    buildings.geometry.apply(lambda g: affine_transform(g, to_site)).values, 0.01)
survey.set_crs(None, allow_override=True).to_file("site_survey.shp")    # no CRS, so no .prj is written
assert not os.path.exists("site_survey.prj")

# four control points: the building corners nearest each corner of the site
corners = shapely.get_coordinates(buildings.geometry.values)
minx, miny, maxx, maxy = buildings.total_bounds
CONTROL = []
for cx, cy in [(minx, miny), (maxx, miny), (maxx, maxy), (minx, maxy)]:
    x, y = corners[np.argmin(np.hypot(corners[:, 0] - cx, corners[:, 1] - cy))]
    site_xy = (to_site[0] * x + to_site[1] * y + to_site[4], to_site[2] * x + to_site[3] * y + to_site[5])
    CONTROL.append((tuple(round(v, 2) for v in site_xy), (round(x, 2), round(y, 2))))

ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
gpd.read_file(f"{ONS}/Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query"
              "?where=LAD24NM%3D%27Leeds%27&outFields=WD24CD,WD24NM&outSR=4326&f=geojson").to_file("study_area.gpkg")
gpd.read_file("https://gisco-services.ec.europa.eu/distribution/v2/nuts/geojson/"
              "NUTS_RG_20M_2024_4326_LEVL_2.geojson").to_file("nuts2.gpkg")
print(f"{len(survey)} surveyed buildings, peg at {PEG}, {len(CONTROL)} control points")
555 surveyed buildings, peg at (429286.0, 432626.0), 4 control points

Example 1: a project-local CRS registry

# in a project this is one module, say crs_registry.py, that everything imports
from pyproj import CRS


class CustomCRS:
    """Named custom coordinate systems for this project, defined once."""

    SITE_LOCAL = CRS.from_proj4(
        "+proj=tmerc +lat_0=53.7952 +lon_0=-1.5474 +k_0=1 "
        "+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

    EUROPE_EQUAL_AREA = CRS.from_proj4(
        "+proj=aea +lat_1=40 +lat_2=65 +lat_0=53 +lon_0=10 "
        "+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

    STUDY_AREA_DISTANCE = CRS.from_proj4(
        "+proj=aeqd +lat_0=53.7952 +lon_0=-1.5474 "
        "+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

    NOTES = {
        "SITE_LOCAL": ("Transverse Mercator on the site centre. Scale factor 1.0 so "
                       "distances are true near the origin. Not for anything beyond "
                       "about 20 km."),
        "EUROPE_EQUAL_AREA": ("Albers, standard parallels 40N and 65N. Equal-area by "
                              "design; shapes sheared at the edges. Use for area "
                              "comparison, not for maps."),
        "STUDY_AREA_DISTANCE": ("Azimuthal equidistant on the study centre. Distances "
                                "from that point are exact; all others are not."),
    }

    @classmethod
    def describe(cls, name):
        crs = getattr(cls, name)
        print(f"{name}")
        print(f"  proj4  {crs.to_proj4()}")
        print(f"  units  {crs.axis_info[0].unit_name}")
        print(f"  note   {cls.NOTES[name]}")
        return crs
CustomCRS.describe("EUROPE_EQUAL_AREA")
areas = gpd.read_file("nuts2.gpkg").to_crs(CustomCRS.EUROPE_EQUAL_AREA).geometry.area
print(f"  {len(areas)} NUTS 2 regions, {areas.sum() / 1e6:,.0f} km²")
EUROPE_EQUAL_AREA
  proj4  +proj=aea +lat_0=53 +lon_0=10 +lat_1=40 +lat_2=65 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs
  units  metre
  note   Albers, standard parallels 40N and 65N. Equal-area by design; shapes sheared at the edges. Use for area comparison, not for maps.
  299 NUTS 2 regions, 5,756,473 km²

The NOTES are the point. A PROJ string scattered through a codebase is unmaintainable; one that says what it is for, what it preserves and where it stops being valid means the next person does not have to reverse-engineer the intent. STUDY_AREA_DISTANCE in particular is easy to misuse — it gives exact distances from one point only, and without the note someone will eventually use it for a general distance matrix.

Example 2: georeferencing a site survey end to end

import numpy as np
import geopandas as gpd
from shapely.affinity import affine_transform
from shapely.geometry import Point

def georeference(gdf, control, target_crs, *, max_residual_m=0.5, verbose=True):
    """Place a local-grid layer on a real CRS using matched control points.

    control: list of ((local_x, local_y), (real_x, real_y)) pairs.
    """
    if gdf.crs is not None:
        raise ValueError(f"the layer already has a CRS ({gdf.crs.to_string()})")
    if len(control) < 3:
        raise ValueError(f"need at least 3 control points, got {len(control)}")

    local = np.array([c[0] for c in control], dtype=float)
    real = np.array([c[1] for c in control], dtype=float)

    A = np.zeros((2 * len(local), 6))
    A[0::2, 0], A[0::2, 1], A[0::2, 4] = local[:, 0], local[:, 1], 1
    A[1::2, 2], A[1::2, 3], A[1::2, 5] = local[:, 0], local[:, 1], 1
    params, *_ = np.linalg.lstsq(A, real.ravel(), rcond=None)

    predicted = (A @ params).reshape(-1, 2)
    residuals = np.linalg.norm(predicted - real, axis=1)
    rms = float(np.sqrt((residuals ** 2).mean()))

    a, b, d, e, xoff, yoff = params
    scale_x = float(np.hypot(a, d))
    scale_y = float(np.hypot(b, e))
    rotation = float(np.degrees(np.arctan2(d, a)))

    if verbose:
        print(f"control points  {len(control)}")
        print(f"scale           {scale_x:.6f} x, {scale_y:.6f} y")
        print(f"rotation        {rotation:+.4f}°")
        print(f"translation     {xoff:,.2f}, {yoff:,.2f}")
        print(f"RMS residual    {rms:.3f} m   worst {residuals.max():.3f} m")
        for i, r in enumerate(residuals):
            flag = "  ⚠" if r > max_residual_m else ""
            print(f"  point {i}: {r:.3f} m{flag}")

    if rms > max_residual_m:
        raise ValueError(
            f"RMS residual {rms:.3f} m exceeds {max_residual_m} m — the control "
            f"points are inconsistent, mis-matched, or the relationship is not affine")

    out = gdf.copy()
    out["geometry"] = gdf.geometry.apply(
        lambda g: affine_transform(g, [a, b, d, e, xoff, yoff]))
    return out.set_crs(target_crs)

survey = gpd.read_file("site_survey.shp")
placed = georeference(survey, control=CONTROL, target_crs=27700)
placed.to_file("site_bng.gpkg", driver="GPKG")
miss = shapely.distance(placed.geometry.centroid.values, buildings.geometry.centroid.values)
print(f"worst building after georeferencing: {miss.max() * 1000:.1f} mm from its OSM position")
control points  4
scale           1.000001 x, 0.999995 y
rotation        +1.7286°
translation     429,286.00, 432,626.00
RMS residual    0.003 m   worst 0.003 m
  point 0: 0.003 m
  point 1: 0.003 m
  point 2: 0.002 m
  point 3: 0.003 m
worst building after georeferencing: 12.2 mm from its OSM position

Three diagnostics make this trustworthy rather than merely functional.

Scale of 1.0 confirms the local grid uses the same units as the target. A scale of 0.3048 would say the survey is in feet; 1.0004 might mean it was measured on a different projection's scale factor.

Rotation of 1.73° is the difference between the site's grid north and the national grid's — a real and expected quantity that the transform recovers rather than assumes. Here the true value is known, because the site grid was made by turning real building outlines by exactly 1.7284°, and the transform recovers it from four corners.

Per-point residuals identify a bad control point. These are all millimetres; one at 4 m would be mis-matched, and the whole solution would be dragged toward it.

Raising when RMS exceeds the tolerance is deliberate: a poor fit means the assumption of an affine relationship is wrong, and proceeding produces data that is plausibly placed and wrong throughout.

Example 3: comparing a custom projection against the alternatives

import numpy as np
import pandas as pd
import geopandas as gpd
from pyproj import CRS, Geod

def evaluate_crs_options(gdf, extra=(), samples=200, seed=0):
    """Compare candidate CRS, including custom ones, on this layer."""
    g = gdf.to_crs(4326)
    minx, miny, maxx, maxy = g.total_bounds
    c = g.union_all().centroid
    geod = Geod(ellps="WGS84")

    candidates = {
        "UTM (estimated)": gdf.estimate_utm_crs(),
        "Web Mercator": CRS.from_epsg(3857),
        "custom LAEA (centred)": CRS.from_proj4(
            f"+proj=laea +lat_0={c.y:.6f} +lon_0={c.x:.6f} "
            f"+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"),
        "custom AEA (centred)": CRS.from_proj4(
            f"+proj=aea +lat_1={miny + (maxy-miny)/6:.4f} "
            f"+lat_2={maxy - (maxy-miny)/6:.4f} +lat_0={c.y:.6f} "
            f"+lon_0={c.x:.6f} +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"),
        **{name: CRS.from_user_input(v) for name, v in extra},
    }

    rng = np.random.default_rng(seed)
    idx = rng.choice(len(g), min(samples, len(g)), replace=False)
    sample = g.iloc[idx]
    true_area = np.array([abs(geod.geometry_area_perimeter(x)[0])
                          for x in sample.geometry])
    true_len = np.array([abs(geod.geometry_area_perimeter(x)[1])
                         for x in sample.geometry])

    rows = []
    for name, crs in candidates.items():
        proj = sample.to_crs(crs)
        with np.errstate(divide="ignore", invalid="ignore"):
            a = proj.geometry.area.to_numpy() / np.where(true_area == 0, np.nan, true_area)
            l = proj.geometry.length.to_numpy() / np.where(true_len == 0, np.nan, true_len)
        rows.append({
            "crs": name,
            "epsg": crs.to_epsg() or "—",
            "area_err_%": round(100 * np.nanmax(np.abs(a - 1)), 3),
            "len_err_%": round(100 * np.nanmax(np.abs(l - 1)), 3),
        })
    df = pd.DataFrame(rows).sort_values("area_err_%")
    print(df.to_string(index=False))
    return df

evaluate_crs_options(gpd.read_file("nuts2.gpkg"),
                     extra=[("EPSG:3035", 3035)])
                  crs  epsg  area_err_%  len_err_%
            EPSG:3035  3035       0.037     34.486
custom LAEA (centred)     —       0.040     34.502
 custom AEA (centred)     —       0.091     35.083
      UTM (estimated) 32630     208.955     75.798
         Web Mercator  3857     515.021    147.813

The equal-area projections — the two custom ones and EPSG:3035 — are all within a tenth of a percent on area, which is what equal-area means, and all distort length badly at the edges: a third out on the worst region, one of the overseas ones thousands of kilometres from any of their centres. EPSG:3035 is a Lambert Azimuthal Equal Area centred on Europe. UTM, excellent for a city, falls apart across a continent — and more so here, because the NUTS file reaches the French overseas regions — because most of the data lies far outside the estimated zone.

The practical conclusion here is to use EPSG:3035, which is identical in accuracy and has a code every tool understands. A custom CRS is worth defining when no standard code fits your extent — not when one already does.

A CRS is a definition, not a number

A custom projection compared with a local grid. only one of these is a coordinate reference systemCustom projection✓has a datum and an ellipsoid✓every place on earth has coordinates✓converts to any other CRS✓PROJ handles it like any EPSG code✓no code is only an interop costLocal grid✕measured from an arbitrary peg✕no relation to the earth✕cannot convert to anything✕needs control points first✕assigning an EPSG code is a fictionA layer honestly marked "no CRS" beats a mislabelled one spatialworkflow.io
One is a definition every point on earth has coordinates in. The other is numbers from a peg.

A coordinate reference system is a definition, not a number. The EPSG code is a lookup key into a registry of definitions, and the definition is what PROJ actually uses. So a CRS with no code is not deficient — it is simply not in the registry, and PROJ handles it identically.

That is why a PROJ string or WKT works everywhere in Python. to_crs("+proj=aea …") performs exactly the same transformation pipeline as to_crs(3035); the only difference is where the parameters came from. The practical costs of having no code are about interoperability: some formats cannot store an arbitrary definition, some tools display only codes, and a code is easier to communicate than a sixty-character string.

The distinction between a custom projection and a local grid is the important one, and they are often conflated.

A custom projection is a genuine CRS: it has a datum, an ellipsoid, and a mathematical mapping from the earth to a plane. Every point on the earth has coordinates in it, and the transformation to any other CRS is well defined. Defining one is normal engineering — projections centred on a study area routinely beat standard codes for unusual extents.

A custom projection and a local engineering grid, and why only one is a CRS. only one of these is a coordinate reference systema custom projectiona local gridis it a CRSyesnohas a datumyesnoreaches the earthby formulaonly through tie pointsto_crs worksyesno — fit an affine firstthe relationship is usually affine: translation, rotation, uniform scale, sometimes a reflectionCoordinates near zero with no tie points cannot be placed on the earth at all spatialworkflow.io
Confusing the two is why a site survey sometimes ends up "in" EPSG:4326 near Null Island.

A local grid is not a CRS at all. It is a set of coordinates measured from an arbitrary origin, with no recorded relationship to the earth. Assigning it an EPSG code does not make it right; the numbers were never in that system. Georeferencing it means discovering the relationship, which requires control points — locations whose coordinates are known in both systems.

That relationship is usually affine, meaning it composes translation, rotation, uniform scale and possibly shear. Three matched point pairs determine six parameters exactly; more than three over-determine it, which is better, because the least-squares residuals then tell you whether the assumption holds. A residual of a few centimetres means the fit is good. A residual of metres means something is wrong — a mis-matched point, a unit mismatch, or a relationship that is not affine because the local survey used a different projection with its own curvature.

Sometimes the honest answer is that it cannot be done. Coordinates of (1000.4, 2408.1) with no tie point and no documentation are unplaceable, and no amount of inference changes that. The temptation is to assign a plausible CRS, which produces data that is internally consistent and confidently mislocated — worse than data honestly marked as having no CRS, because the error is now invisible.

Finally, the interoperability constraint. GeoPackage, FlatGeobuf, shapefile and GeoParquet all store WKT and preserve a custom CRS. GeoJSON does not, and cannot: RFC 7946 mandates WGS 84, so GDAL reprojects on write. That is a specification decision rather than a limitation, and it means GeoJSON is the wrong container for anything in a custom projection. Keeping a sidecar JSON with the WKT and PROJ string is cheap insurance against a conversion elsewhere in the chain quietly discarding it.

Missing EPSG codes, GeoJSON and affine matrices

  • crs.to_epsg() returns None for a custom CRS. That is expected, not an error.
  • GeoJSON cannot store a custom CRS. RFC 7946 mandates WGS 84 and GDAL reprojects on write.
  • Shapely's affine matrix is [a, b, d, e, xoff, yoff] — easy to transpose. Test with a known point.
  • Three control points determine an affine transform exactly, so residuals are zero and tell you nothing. Use four or more.
  • A scale far from 1.0 in a fitted transform usually means a unit mismatch — feet versus metres.
  • +towgs84= is a crude datum shift. Prefer a proper transformation with a grid file where accuracy matters.
  • +proj=aeqd gives exact distances from its centre only. It is not a general-purpose CRS.
  • PROJ strings are lossy relative to WKT2 — they cannot express every datum or transformation nuance.
  • CRS.equals() compares definitions, so two equivalent strings written differently may compare unequal.
  • Name your custom CRS in the WKT. "unknown" in a file six months later is a small archaeology project.

FAQ

Can I use a CRS with no EPSG code?

Yes. PROJ works from the definition, not the code, so a PROJ string or WKT behaves identically. The only costs are interoperability and communication.

How do I define a custom projection?

CRS.from_proj4("+proj=aea +lat_1=… +lat_2=… +lat_0=… +lon_0=… +datum=WGS84 +units=m +no_defs"), or CRS.from_wkt(...) for the more expressive form.

My survey has coordinates starting at zero. What CRS is it?

None. It is a local grid, and it has no relationship to the earth until you supply control points whose coordinates are known in both systems.

How many control points do I need?

At least three for an affine transform, and four or more so the residuals mean something — three fit exactly and tell you nothing about accuracy.

Why did my custom CRS disappear after saving?

You saved to GeoJSON, which mandates WGS 84 by specification. Use GeoPackage, FlatGeobuf or GeoParquet, all of which store WKT.

Should I use a custom projection or a standard one?

A standard one whenever it fits — every tool understands the code. Define a custom projection when no standard CRS suits your extent, and verify it against geodesic truth.

What if I cannot identify a layer's CRS?

Leave it unset and ask the supplier. A layer honestly marked as having no CRS is usable once identified; a mislabelled one is wrong in a way nothing will catch.