How to Work with a Custom or Local CRS in Python #

Problem statement #

Not every coordinate system has an EPSG code.

A site survey arrives with coordinates like (1000.4, 2408.1) measured from a corner peg. 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.

gdf = gpd.read_file("site_survey.shp")
print(gdf.crs)                       # None
print(gdf.total_bounds)              # [0.0, 0.0, 4820.1, 3104.6]

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 #

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.

Step-by-step solution #

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)                              # unknown
print(europe.is_projected)                      # True
print(europe.axis_info[0].unit_name)            # metre
print(europe.to_epsg())                         # None — expected

The common projection parameters:

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,

Give it a name so it is legible later:

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

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))

crs = centred_crs(gpd.read_file("study_area.gpkg"), kind="laea")
print(crs.to_proj4())

+proj=laea +lat_0=53.480800 +lon_0=-2.242600 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

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 = gdf.to_crs(4326).sample(100, random_state=0)
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: {100 * np.max(np.abs(got / truth - 1)):.3f}%")

3. Handle a local engineering grid #

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 origin, in a real CRS
SITE_ORIGIN_BNG = (383_618.5, 398_050.4)      # 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
georeferenced = tie_local_grid(site, SITE_ORIGIN_BNG, 27700)
print(georeferenced.total_bounds)

[383618.5 398050.4 388438.6 401155.0]

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.

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=[(0, 0), (4820.1, 0), (4820.1, 3104.6), (0, 3104.6)],
    real_pts=[(383618.5, 398050.4), (388436.2, 398195.8),
              (388531.4, 401296.9), (383713.7, 401151.5)],
)

4 control points, RMS residual 0.084 m, worst 0.142 m

The RMS residual is the check. Under a few centimetres means the control points are consistent and the transform fits; a metre means 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="/tmp"):
    """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 = gpd.read_file("study_area.gpkg").to_crs(centred_crs(study_area))
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.

Code examples #

Example 1: a project-local CRS registry #

# crs_registry.py
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.4808 +lon_0=-2.2426 +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.4808 +lon_0=-2.2426 "
        "+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. Areas exact to "
                              "0.02% across the study extent; 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

import geopandas as gpd
from crs_registry import CustomCRS

CustomCRS.describe("EUROPE_EQUAL_AREA")
areas = gpd.read_file("nuts2.gpkg").to_crs(CustomCRS.EUROPE_EQUAL_AREA).geometry.area

EUROPE_EQUAL_AREA
  proj4  +proj=aea +lat_1=40 +lat_2=65 +lat_0=53 +lon_0=10 +x_0=0 +y_0=0 …
  units  metre
  note   Albers, standard parallels 40N and 65N. Areas exact to 0.02% across the
         study extent; shapes sheared at the edges. Use for area comparison, not
         for maps.

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=[
    ((0.00,    0.00),    (383618.50, 398050.40)),
    ((4820.10, 0.00),    (388436.20, 398195.80)),
    ((4820.10, 3104.60), (388531.40, 401296.90)),
    ((0.00,    3104.60), (383713.70, 401151.50)),
], target_crs=27700)
placed.to_file("site_bng.gpkg", driver="GPKG")

control points  4
scale           1.000012 x, 1.000008 y
rotation        +1.7284°
translation     383,618.50, 398,050.40
RMS residual    0.084 m   worst 0.142 m
  point 0: 0.041 m
  point 1: 0.092 m
  point 2: 0.142 m
  point 3: 0.048 m

Three diagnostics make this trustworthy rather than merely functional.

Scale near 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.

Per-point residuals identify a bad control point. Point 2 at 14 cm against 4 cm elsewhere suggests it deserves a second look; if it were 4 m, it 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_%
 custom AEA (centred)      —       0.018      0.204
custom LAEA (centred)      —       0.024      0.312
            EPSG:3035   3035       0.021      0.208
      UTM (estimated)  32632       8.412      4.118
         Web Mercator   3857     161.240     61.204

The custom projections match EPSG:3035 to within a hundredth of a percent, which is expected — EPSG:3035 is a Lambert Azimuthal Equal Area centred on Europe, with parameters chosen for roughly this extent. And UTM, excellent for a 60 km study area, is 8.4% off across a continent because most of the data lies 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.

Explanation #

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 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.

Edge cases or notes #

• 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.

Internal links #

• Coordinate reference systems explained for Python GIS — the underlying concepts
• EPSG codes explained: how to choose the right CRS in Python — the registry these sit outside
• How to choose the right projected CRS for your study area — when a standard code fits
• How to set, assign and convert CRS in GeoPandas — labelling versus converting
• CRS not found error in GeoPandas — when a code cannot be resolved
• Projected vs geographic CRS — what a projected definition provides
• How to reproject between datums correctly — the datum half of a definition
• GIS vector file formats compared — which formats keep a custom CRS

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.