Choosing a Projected CRS for Area and Distance

Every projection is a compromise, and the compromise is chosen — usually by whoever set up the first dataset, sometimes by a default in a tool, occasionally on purpose. When a pipeline computes areas for billing, distances for service coverage or buffers for regulatory setbacks, that inherited choice becomes an error budget nobody has quantified. This guide makes the choice explicit: deciding which quantity must be preserved, shortlisting projections that actually cover the extent, measuring the distortion each one introduces against geodesic truth, and recording the decision with its numbers. It puts a figure on the CRS discipline described in Coordinate Reference System Precision Standards.

Prerequisites

• pyproj 3.6+, shapely 2.0+, geopandas 0.14+. pyproj.Geod supplies the geodesic reference values.
• The data’s true extent, not its nominal one. A “county” layer that includes an offshore asset spans a much larger area than expected, and the outlier drives the worst-case distortion.
• A written tolerance. “Area must be correct to within 0.1%” makes this a decision; without it, every projection looks acceptable.
• The list of candidate systems your organisation already uses. Introducing a new CRS has an integration cost that a marginal accuracy gain rarely justifies.

Step-by-Step Procedure

Step 1 — State what must be preserved

Most national grids are conformal (Transverse Mercator or Lambert Conformal Conic) with a scale factor chosen so that scale error stays within a few parts in ten thousand across the country. For most national-scale work that is well inside any practical tolerance for both area and distance, which is why “use the national grid” is usually the right answer and the interesting cases are the ones where it is not.

Step 2 — Shortlist by extent

# crs_choice/step1_candidates.py
from pyproj import CRS
from pyproj.aoi import AreaOfInterest
from pyproj.database import query_utm_crs_info, query_crs_info


def candidates_for(bounds: tuple[float, float, float, float]) -> list[dict]:
    """Projected systems whose declared area of use covers the data extent."""
    west, south, east, north = bounds
    aoi = AreaOfInterest(west, south, east, north)

    utm = [{"code": f"EPSG:{c.code}", "name": c.name, "kind": "utm"}
           for c in query_utm_crs_info(datum_name="WGS 84", area_of_interest=aoi)]

    others = [{"code": f"{c.auth_name}:{c.code}", "name": c.name, "kind": "projected"}
              for c in query_crs_info(pj_types="PROJECTED_CRS", area_of_interest=aoi,
                                      contains=True)][:20]
    return utm + others


def extent_span_km(bounds) -> tuple[float, float]:
    from pyproj import Geod
    g = Geod(ellps="WGS84")
    west, south, east, north = bounds
    ew = g.inv(west, (south + north) / 2, east, (south + north) / 2)[2] / 1000
    ns = g.inv((west + east) / 2, south, (west + east) / 2, north)[2] / 1000
    return round(ew, 1), round(ns, 1)

Verification: if extent_span_km shows an east-west span above roughly 400 km, a single UTM zone will not cover it cleanly — zones are 6 degrees wide and distortion grows toward their edges. That single number eliminates most bad choices immediately.

Step 3 — Measure the distortion against geodesic truth

# crs_choice/step2_measure.py
import geopandas as gpd
import numpy as np
from pyproj import Geod

GEOD = Geod(ellps="WGS84")


def geodesic_area_m2(geom) -> float:
    """True area on the ellipsoid, independent of any projection."""
    if geom.geom_type == "MultiPolygon":
        return sum(geodesic_area_m2(p) for p in geom.geoms)
    lon, lat = geom.exterior.coords.xy
    area, _perimeter = GEOD.polygon_area_perimeter(list(lon), list(lat))
    return abs(area)


def distortion_report(gdf_wgs84: gpd.GeoDataFrame, candidate_epsg: int,
                      sample: int = 400, seed: int = 20260806) -> dict:
    """Projected area and distance error versus geodesic values, across the extent."""
    s = gdf_wgs84.sample(min(sample, len(gdf_wgs84)), random_state=seed)
    truth = s.geometry.map(geodesic_area_m2).to_numpy()
    projected = s.to_crs(candidate_epsg).area.to_numpy()

    with np.errstate(divide="ignore", invalid="ignore"):
        area_err = 100.0 * (projected - truth) / truth

    pts = s.geometry.representative_point()
    lons, lats = pts.x.to_numpy(), pts.y.to_numpy()
    proj_pts = gpd.GeoSeries(pts, crs=4326).to_crs(candidate_epsg)
    px, py = proj_pts.x.to_numpy(), proj_pts.y.to_numpy()

    geo_d, proj_d = [], []
    for i in range(len(lons) - 1):
        geo_d.append(GEOD.inv(lons[i], lats[i], lons[i + 1], lats[i + 1])[2])
        proj_d.append(float(np.hypot(px[i + 1] - px[i], py[i + 1] - py[i])))
    geo_d, proj_d = np.array(geo_d), np.array(proj_d)
    dist_err = 100.0 * (proj_d - geo_d) / np.where(geo_d == 0, np.nan, geo_d)

    return {
        "epsg": candidate_epsg,
        "area_err_median_pct": round(float(np.nanmedian(area_err)), 5),
        "area_err_p95_abs_pct": round(float(np.nanpercentile(np.abs(area_err), 95)), 5),
        "area_err_max_abs_pct": round(float(np.nanmax(np.abs(area_err))), 5),
        "distance_err_p95_abs_pct": round(float(np.nanpercentile(np.abs(dist_err), 95)), 5),
        "sampled_features": int(len(s)),
    }

Verification: run the report for three or four candidates and compare. The numbers are usually decisive: a national grid across a national extent produces area errors in the fourth decimal place, a correctly chosen UTM zone similar, an adjacent UTM zone an order of magnitude worse, and Web Mercator errors of tens of percent.

Step 4 — Compare against the tolerance

# crs_choice/step3_decide.py
def assess(report: dict, area_tolerance_pct: float, distance_tolerance_pct: float) -> dict:
    ok_area = report["area_err_p95_abs_pct"] <= area_tolerance_pct
    ok_dist = report["distance_err_p95_abs_pct"] <= distance_tolerance_pct
    return {
        **report,
        "area_within_tolerance": ok_area,
        "distance_within_tolerance": ok_dist,
        "verdict": ("suitable" if ok_area and ok_dist else
                    "area only" if ok_area else
                    "distance only" if ok_dist else "unsuitable"),
        "headroom_pct": round(area_tolerance_pct - report["area_err_p95_abs_pct"], 5),
    }

Verification: headroom_pct is the number worth keeping. A choice that only just satisfies the tolerance will fail as soon as the extent grows, and extents always grow.

Step 5 — Record the decision where the data lives

# crs_choice/step4_record.py
import json
from datetime import date


def record(chosen_epsg: int, assessment: dict, rationale: str, extent: tuple, path: str):
    doc = {
        "decided_on": date.today().isoformat(),
        "crs": f"EPSG:{chosen_epsg}",
        "rationale": rationale,
        "extent_wgs84": list(extent),
        "measured_error": {
            "area_p95_abs_pct": assessment["area_err_p95_abs_pct"],
            "distance_p95_abs_pct": assessment["distance_err_p95_abs_pct"],
            "sampled_features": assessment["sampled_features"],
        },
        "review_when": "the data extent changes materially, or the tolerance is tightened",
    }
    with open(path, "w", encoding="utf-8") as fh:
        json.dump(doc, fh, indent=2)
    return doc

Verification: the record should let a successor answer “why this projection” in one read. It also feeds the conformance metadata described in Aligning Spatial Metadata with ISO 19157, where positional accuracy statements need exactly this provenance.

Interpreting Results

Splitting the concerns is legitimate and underused. Nothing prevents storing data in the national grid and computing areas by reprojecting to an equal-area system for that one calculation — provided the practice is documented, because two different area figures for the same parcel is exactly the kind of inconsistency that erodes trust in a dataset.

Gotchas & Edge Cases

GeoSeries.area gives no warning in geographic coordinates. Computing area on an EPSG:4326 layer returns square degrees, silently. Assert the CRS is projected before any metric operation — that check belongs in the rule set, not in a reviewer’s memory.

A projection’s “area of use” is advisory, not enforced. pyproj will happily transform coordinates far outside a system’s declared extent, producing plausible numbers with large errors. Check the extent against the area of use explicitly.

Scale factor is not uniform within a zone. A Transverse Mercator grid has a scale factor below one at the central meridian and above one at the zone edges. Distances measured near the central meridian are slightly short and near the edges slightly long — usually a few parts in ten thousand, occasionally material for engineering work.

Equal-area projections distort shape, and it shows. Buffers become slightly elliptical and angles are not preserved. For point-in-polygon and area work this is irrelevant; for anything involving bearings it is not.

Reprojection is not free at scale. Transforming a large layer for a single calculation costs real time. Where area is computed frequently, store a pre-computed area column in the correct system rather than reprojecting per query.

Vertical measurements are a separate problem. Choosing a horizontal projection says nothing about heights, which need their own datum decision — see the vertical datum discussion in Raster and Elevation Data Quality Checks.

When to Escalate

• No candidate meets the tolerance — that is a genuine finding, and the answer is usually to compute on the ellipsoid directly with geodesic functions rather than to accept a projection that cannot deliver.
• Two teams using different projections for the same layer will publish different areas for the same parcel. Settle it at the governance level, not per project, following Defining Spatial Data Quality Policies.
• A legally significant area figure — tax, ownership, subsidy — should have its projection choice reviewed and recorded formally. The error is small; the exposure is not.
• Changing the CRS of an established dataset affects every stored coordinate, every cached area and every downstream join. Treat it as a migration with a plan, not as a configuration change.

Frequently Asked Questions

Related

• Coordinate Reference System Precision Standards — the precision model this projection choice feeds
• Reprojecting Mixed CRS Datasets with pyproj — executing the transformation once the target is chosen
• Setting Decimal Precision for Survey Boundaries — how many decimals the chosen unit needs

Back to Coordinate Reference System Precision Standards