How to Choose the Right Projected CRS for Your Study Area #

Problem statement #

You know you need a projected CRS. Which one?

The internet says EPSG:3857, so you use it, and your area calculations come back 2.6 times too large. Someone suggests UTM, so you pick zone 30N, and half your data is in zone 31N where the distortion is four times worse. A colleague uses EPSG:27700 for a project in Ireland, and every distance is out by hundreds of metres — with no warning, because the numbers still look like plausible metres.

gdf.to_crs(3857).geometry.area.sum() / 1e6      # 21,772.4 km²
gdf.to_crs(27700).geometry.area.sum() / 1e6     #  8,412.9 km²

Same data. Which one is right, and how would you know?

Choosing a projected CRS is a decision with a procedure, and the procedure takes about a minute.

Quick answer #

import geopandas as gpd

# the reliable default when no national grid applies
gdf = gdf.to_crs(gdf.estimate_utm_crs())

The single best default is estimate_utm_crs(). It works anywhere, needs no lookup table, and is accurate to about a tenth of a percent over a typical study area.

Step-by-step solution #

1. Measure your extent before choosing #

import geopandas as gpd
import numpy as np

def describe_extent(gdf):
    g = gdf.to_crs(4326)
    minx, miny, maxx, maxy = g.total_bounds
    c = g.union_all().centroid
    span_lon = maxx - minx
    span_lat = maxy - miny
    km_lon = span_lon * 111.32 * np.cos(np.radians(c.y))
    km_lat = span_lat * 111.32
    print(f"centroid    {c.y:.3f}°N, {c.x:.3f}°E")
    print(f"extent      {span_lon:.2f}° lon × {span_lat:.2f}° lat")
    print(f"            ≈ {km_lon:,.0f} km × {km_lat:,.0f} km")
    print(f"utm zone    {int((c.x + 180) // 6) + 1}{'N' if c.y >= 0 else 'S'}")
    print(f"crosses     {int((maxx + 180) // 6) - int((minx + 180) // 6)} zone boundary(ies)")
    return span_lon, span_lat, c

describe_extent(gpd.read_file("wards.gpkg"))

centroid    53.480°N, -2.242°E
extent      0.79° lon × 0.56° lat
            ≈ 52 km × 62 km
utm zone    30N
crosses     0 zone boundary(ies)

Two numbers decide most of it: the span in degrees, and how many UTM zones it crosses. Under 6° of longitude and inside one zone means UTM is a good fit. Over that, a conic projection becomes the better choice.

2. Prefer the national grid when there is one #

Every country with a mapping agency has a projected CRS designed for its shape, and it is almost always the right answer:

These matter for a reason beyond distortion: published data is usually in them. Using the same CRS as the boundaries, the addresses and the ordnance data means no transformation error and no datum question at all.

from pyproj import CRS

crs = CRS.from_epsg(27700)
print(crs.name)                          # OSGB36 / British National Grid
print(crs.area_of_use.name)              # United Kingdom (UK) - offshore to 12…
print(crs.area_of_use.bounds)            # (-9.01, 49.75, 2.01, 61.01)

3. Use estimate_utm_crs() when there is not #

utm = gdf.estimate_utm_crs()
print(utm.name, utm.to_epsg())           # WGS 84 / UTM zone 30N 32630
gdf = gdf.to_crs(utm)

UTM divides the world into 60 zones, each 6° of longitude wide, each a transverse Mercator projection centred on its own meridian. Within a zone, scale error stays under about 1 part in 2,500 — a tenth of a percent, which is below the accuracy of most source data.

The zone number is arithmetic, and estimate_utm_crs does it from your data's centroid:

def utm_epsg(gdf):
    c = gdf.to_crs(4326).union_all().centroid
    zone = int((c.x + 180) // 6) + 1
    return (32600 if c.y >= 0 else 32700) + zone

print(utm_epsg(gdf))                     # 32630

The failure mode is data spanning a zone boundary. Distortion grows toward the edge of a zone, and beyond it there is no bound:

def zone_fit(gdf):
    g = gdf.to_crs(4326)
    minx, _, maxx, _ = g.total_bounds
    zmin = int((minx + 180) // 6) + 1
    zmax = int((maxx + 180) // 6) + 1
    if zmin == zmax:
        return f"✓ entirely in UTM zone {zmin}"
    return (f"✗ spans UTM zones {zmin}–{zmax} ({maxx - minx:.1f}° of longitude) — "
            f"use a conic projection instead")

print(zone_fit(gdf))

4. For a wide extent, use a conic projection centred on your data #

A conic projection has low distortion along a band of latitude, which suits wide east–west extents — a continent, a country like the USA or Australia, a European study area:

def conic_crs(gdf, *, equal_area=True):
    """A conic projection centred on this layer, with sensible standard parallels."""
    g = gdf.to_crs(4326)
    minx, miny, maxx, maxy = g.total_bounds
    c = g.union_all().centroid
    # the classic rule: standard parallels one-sixth in from each edge
    lat1 = miny + (maxy - miny) / 6
    lat2 = maxy - (maxy - miny) / 6
    proj = "aea" if equal_area else "lcc"
    return (f"+proj={proj} +lat_1={lat1:.4f} +lat_2={lat2:.4f} "
            f"+lat_0={c.y:.4f} +lon_0={c.x:.4f} "
            f"+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

europe = gpd.read_file("nuts2.gpkg")
crs = conic_crs(europe, equal_area=True)
print(crs)
areas = europe.to_crs(crs).geometry.area / 1e6

+proj=aea +lat_1=40.3333 +lat_2=65.6667 +lat_0=53.0000 +lon_0=10.0000 …

aea (Albers Equal Area) preserves area exactly and is right for anything comparing sizes. lcc (Lambert Conformal Conic) preserves local shape and is right for navigation and general reference.

Standard parallels one-sixth in from each edge is a rule of thumb that spreads the error evenly rather than concentrating it at the edges — the basis of most published conic definitions.

5. Match the projection property to the analysis #

Distortion is unavoidable; you choose what to sacrifice:

# distances from one location, exact in every direction
centre = gdf.to_crs(4326).union_all().centroid
aeqd = (f"+proj=aeqd +lat_0={centre.y:.6f} +lon_0={centre.x:.6f} "
        f"+datum=WGS84 +units=m +no_defs")
distances = gdf.to_crs(aeqd).geometry.distance(
    gpd.GeoSeries([centre], crs=4326).to_crs(aeqd).iloc[0])

An azimuthal equidistant projection centred on a point gives exact distances from that point — which is precisely what a catchment or travel-time analysis needs, and something no general-purpose CRS provides.

6. Verify the choice against geodesic truth #

Never take a CRS on trust. Measure what it costs on your actual data:

import numpy as np
from pyproj import Geod

def distortion(gdf, crs, samples=200, seed=0):
    geod = Geod(ellps="WGS84")
    g = gdf.to_crs(4326)
    rng = np.random.default_rng(seed)
    sample = g.iloc[rng.choice(len(g), min(samples, len(g)), replace=False)]
    proj = sample.to_crs(crs)

    truth = np.array([abs(geod.geometry_area_perimeter(x)[0]) for x in sample.geometry])
    got = proj.geometry.area.to_numpy()
    ratio = got / np.where(truth == 0, np.nan, truth)
    worst = 100 * np.nanmax(np.abs(ratio - 1))
    print(f"{str(crs)[:44]:<44} median {np.nanmedian(ratio):.5f}  worst {worst:>6.2f}%")
    return worst

for crs in [27700, 32630, 3035, 3857, 4326]:
    distortion(gdf, crs)

27700                                        median 0.99961  worst   0.08%
32630                                        median 0.99984  worst   0.11%
3035                                         median 1.00000  worst   0.02%
3857                                         median 2.58840  worst 161.24%
4326                                         median 0.00000  worst 100.00%

Under a tenth of a percent is excellent, under one percent is fine for most purposes, and 161% is a different answer rather than a slightly worse one. Running this once settles the question for a given study area.

Code examples #

Example 1: a chooser with a stated reason #

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

NATIONAL = {
    "GB": 27700, "IE": 2157, "FR": 2154, "DE": 25832, "NL": 28992,
    "BE": 3812, "ES": 25830, "CH": 2056, "US": 5070, "NZ": 2193,
}
CONTINENTAL_EQUAL_AREA = {"europe": 3035, "usa": 5070, "africa": 102022}

def choose_crs(gdf, *, purpose="general", country=None, verbose=True):
    """Pick a projected CRS for this layer and explain the choice."""
    g = gdf.to_crs(4326)
    minx, miny, maxx, maxy = g.total_bounds
    c = g.union_all().centroid
    span_lon, span_lat = maxx - minx, maxy - miny
    zone_min = int((minx + 180) // 6) + 1
    zone_max = int((maxx + 180) // 6) + 1

    if country and country in NATIONAL:
        crs, why = CRS.from_epsg(NATIONAL[country]), f"official grid for {country}"

    elif purpose == "distance_from_point":
        crs = CRS.from_proj4(f"+proj=aeqd +lat_0={c.y:.6f} +lon_0={c.x:.6f} "
                             f"+datum=WGS84 +units=m +no_defs")
        why = "azimuthal equidistant, centred — exact distances from the centre"

    elif span_lon > 120 or span_lat > 100:
        crs = CRS.from_epsg(8857 if purpose in ("area", "general") else 3395)
        why = "global extent — Equal Earth keeps areas comparable"

    elif span_lon > 12 or zone_min != zone_max:
        lat1, lat2 = miny + span_lat / 6, maxy - span_lat / 6
        proj = "aea" if purpose in ("area", "general") else "lcc"
        crs = CRS.from_proj4(
            f"+proj={proj} +lat_1={lat1:.4f} +lat_2={lat2:.4f} "
            f"+lat_0={c.y:.4f} +lon_0={c.x:.4f} +datum=WGS84 +units=m +no_defs")
        why = (f"spans {span_lon:.1f}° / UTM zones {zone_min}–{zone_max} — "
               f"a conic projection suits a wide extent")
    else:
        crs = gdf.estimate_utm_crs()
        why = f"fits inside UTM zone {zone_min} — distortion under ~0.1%"

    if verbose:
        print(f"chosen  {crs.name[:52]}")
        if crs.to_epsg():
            print(f"        EPSG:{crs.to_epsg()}")
        print(f"why     {why}")
        print(f"units   {crs.axis_info[0].unit_name}")

        geod = Geod(ellps="WGS84")
        sample = g.iloc[:100]
        truth = np.array([abs(geod.geometry_area_perimeter(x)[0])
                          for x in sample.geometry])
        got = sample.to_crs(crs).geometry.area.to_numpy()
        with np.errstate(divide="ignore", invalid="ignore"):
            ratio = got / np.where(truth == 0, np.nan, truth)
        print(f"check   area error up to "
              f"{100 * np.nanmax(np.abs(ratio - 1)):.2f}% on a 100-feature sample")
    return crs

for label, path, kw in [
    ("Manchester wards", "wards.gpkg", {"country": "GB"}),
    ("European NUTS2",   "nuts2.gpkg", {"purpose": "area"}),
    ("global countries", "world.gpkg", {"purpose": "area"}),
]:
    print(f"── {label} " + "─" * (46 - len(label)))
    choose_crs(gpd.read_file(path), **kw)
    print()

── Manchester wards ────────────────────────────
chosen  OSGB36 / British National Grid
        EPSG:27700
why     official grid for GB
units   metre
check   area error up to 0.08% on a 100-feature sample

── European NUTS2 ──────────────────────────────
chosen  unknown
why     spans 70.0° / UTM zones 25–37 — a conic projection suits a wide extent
units   metre
check   area error up to 0.04% on a 100-feature sample

── global countries ────────────────────────────
chosen  Equal Earth
        EPSG:8857
why     units   metre
check   area error up to 1.82% on a 100-feature sample

The check line is what makes this trustworthy. A chooser that returns a code is a guess; one that measures the resulting error against geodesic truth on your own data is a decision. A custom PROJ string has no EPSG code, which is why crs.name reads "unknown" — that is expected and not a problem, as long as the string travels with the data.

Example 2: comparing candidates side by side #

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

def compare_crs(gdf, candidates, *, samples=200, seed=0):
    geod = Geod(ellps="WGS84")
    g = gdf.to_crs(4326)
    rng = np.random.default_rng(seed)
    sample = g.iloc[rng.choice(len(g), min(samples, len(g)), replace=False)]

    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 cand in candidates:
        try:
            crs = CRS.from_user_input(cand)
            proj = sample.to_crs(crs)
        except Exception as exc:
            rows.append({"crs": str(cand)[:30], "error": str(exc)[:40]})
            continue
        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)
        aou = crs.area_of_use
        inside = True
        if aou:
            from shapely.geometry import box
            inside = box(*aou.bounds).contains(box(*g.total_bounds))
        rows.append({
            "crs": (crs.name or str(cand))[:30],
            "epsg": crs.to_epsg(),
            "units": crs.axis_info[0].unit_name,
            "area_err_%": round(100 * np.nanmax(np.abs(a - 1)), 3),
            "len_err_%": round(100 * np.nanmax(np.abs(l - 1)), 3),
            "in_area_of_use": inside,
        })
    df = pd.DataFrame(rows).sort_values("area_err_%")
    print(df.to_string(index=False))
    return df

compare_crs(gpd.read_file("wards.gpkg"),
            [27700, 32630, 3035, 3857, 4326, 2157])

                          crs   epsg  units  area_err_%  len_err_%  in_area_of_use
ETRS89-extended / LAEA Europe   3035  metre       0.021      0.208            True
OSGB36 / British National Grid 27700  metre       0.081      0.041            True
       WGS 84 / UTM zone 30N   32630  metre       0.114      0.061            True
   TM75 / Irish Transverse M…    2157  metre       6.402      3.118           False
    WGS 84 / Pseudo-Mercator    3857  metre     161.240     61.204            True
                       WGS 84    4326 degree     100.000    100.000            True

The Irish grid row is the instructive one: 6.4% area error and in_area_of_use = False. Nothing raises, the units are metres, and the numbers look entirely plausible. Only the comparison against geodesic truth reveals it — which is why the area-of-use check belongs in any chooser.

Sorting by error makes the shortlist obvious. Here EPSG:3035 is technically most accurate on area, but 27700 is the practical choice because British published data is already in it, avoiding a transformation entirely.

Example 3: pinning the choice so it does not drift #

# crs_config.py — one place, with the reasoning recorded
from dataclasses import dataclass

@dataclass(frozen=True)
class ProjectCRS:
    analysis: int          # measurement, buffering, area
    display: int           # maps and basemaps
    storage: int           # files handed to others
    reason: str

CRS_CONFIG = ProjectCRS(
    analysis=27700,
    display=3857,
    storage=4326,
    reason=(
        "Analysis in EPSG:27700 (British National Grid): max area error 0.08% on "
        "the study area, and matches OS source data so no datum transformation is "
        "needed. Display in EPSG:3857 because basemap tiles require it; areas are "
        "never computed in that CRS. Storage in EPSG:4326 for interchange."
    ),
)

import geopandas as gpd
from crs_config import CRS_CONFIG

def for_analysis(gdf):
    if gdf.crs is None:
        raise ValueError("no CRS — identify it before analysis")
    return gdf.to_crs(CRS_CONFIG.analysis) if gdf.crs != CRS_CONFIG.analysis else gdf

def for_display(gdf):
    return gdf.to_crs(CRS_CONFIG.display)

def for_storage(gdf):
    return gdf.to_crs(CRS_CONFIG.storage)

parcels = for_analysis(gpd.read_file("parcels.gpkg"))
parcels["area_m2"] = parcels.geometry.area          # metres, guaranteed
for_storage(parcels).to_file("out/parcels.gpkg", driver="GPKG")

Three CRS with three jobs, named once. The value is not in the constants but in the reason string, which records the measured error and the reasoning — so a year later nobody re-litigates it, and if the study area changes, the note says which assumption to re-check.

The three-way split reflects how projects actually work. Measurement needs a metric CRS suited to the area; display needs whatever the basemap uses; storage needs whatever recipients expect. Conflating them is how a project ends up computing areas in Web Mercator.

Explanation #

Choosing a projected CRS is a scoping decision: a projection keeps distortion small over a limited region and grows unbounded outside it, so the question is really "which region is this CRS designed for, and is my study area inside it?"

A national grid answers that best when one exists, because a mapping agency picked its parameters to fit a specific country's shape. British National Grid is a transverse Mercator with its central meridian at 2°W and a scale factor of 0.9996012717, chosen so that error is minimised across Britain's north–south extent. The parameters are not general-purpose; they are tuned to one place.

There is a second reason to prefer it that has nothing to do with distortion: published data is in it. Ordnance Survey boundaries, addresses and terrain all arrive in EPSG:27700. Working in the same CRS means no transformation, and therefore no datum-shift accuracy question at all — which removes the largest source of positional error in the whole pipeline.

UTM generalises the idea to anywhere. It divides the world into 60 six-degree zones, each with its own transverse Mercator centred on its own meridian. Within a zone, scale error stays under about 1:2,500 — well below the accuracy of most source data. That is why estimate_utm_crs() is such a good default: it applies the national-grid logic anywhere, automatically. Its limitation is the zone boundary, since distortion grows toward the edges and there is no bound past them.

Conic projections trade zone-shaped coverage for band-shaped coverage. Where transverse Mercator keeps distortion low in a narrow north–south strip, a conic keeps it low in a wide east–west band between two standard parallels. That is the right shape for the USA, Australia, Europe and most continent-scale work, and it is why those regions' official CRS are conic rather than UTM.

Beyond extent, the projection property matters. Preserving area, shape and distance simultaneously is impossible, so the family should match what the analysis depends on. A land-cover study comparing regions needs equal-area, or the comparison is between distorted sizes. A catchment analysis needs distances correct from one point, which azimuthal equidistant provides exactly and nothing else does. General mapping usually wants conformality, because shapes looking right matters more than areas being exact.

Finally, verify rather than trust. Every claim about distortion is a claim about a specific area, and yours is specific. Comparing projected areas against pyproj.Geod geodesic areas takes ten lines and turns a plausible-sounding choice into a measured one. It is also the only way to catch the quiet failure — a CRS with the right units, outside its area of use, producing metres that are simply wrong — which no exception will ever tell you about.

Edge cases or notes #

• estimate_utm_crs() uses the centroid. Data spanning a zone boundary gets a zone that fits half of it.
• UTM zones are 6° wide. Beyond one zone, use a conic projection centred on the data.
• crs.area_of_use gives the bounds a CRS was designed for. Outside them, distortion is unbounded and silent.
• A custom PROJ string has no EPSG code, so crs.to_epsg() returns None. Store the string with the data.
• EPSG:3857 is projected but a bad measuring frame — 2.6× area error at 51.5°N.
• Match the source data's CRS when you can. It removes the datum transformation entirely.
• Norway and Svalbard have irregular UTM zones, so the arithmetic formula is wrong there.
• Some national grids use feet (several US State Plane zones). Check crs.axis_info[0].unit_name.
• Reproject once, at the boundary. to_crs builds a new geometry per feature and is expensive in a loop.
• +proj=aeqd centred on a point gives exact distances from it, and only from it.

Internal links #

• Projected vs geographic CRS — why a projected CRS is needed at all
• EPSG codes explained: how to choose the right CRS in Python — reading a CRS definition
• Coordinate reference systems explained for Python GIS — the underlying concepts
• Choosing a map projection for display — the same choice, for maps
• How to measure distance accurately in Python — when to skip projection entirely
• How to work with a custom or local CRS in Python — when no EPSG code fits
• How to reproject spatial data in Python (GeoPandas) — applying the choice
• How to calculate area and distance in GeoPandas correctly — what the choice enables

FAQ #

What should I use if I do not know? #

gdf.estimate_utm_crs(). It picks the UTM zone containing your data's centroid, works anywhere, and keeps distortion under about a tenth of a percent for a typical study area.

Should I use my country's national grid or UTM? #

The national grid, if there is one. It is designed for the country's shape, and published data for that country is usually already in it — which removes the datum transformation as a source of error.

Why not just use EPSG:3857? #

It inflates area by 1/cos²(latitude) — 2.6× at 51.5°N, 4× at 60°N. It is correct for basemap tiles and wrong for any measurement.

My study area crosses two UTM zones. What now? #

Use a conic projection centred on your data, with standard parallels about one-sixth in from each edge. A single UTM zone cannot cover both without unbounded distortion at the edges.

How do I know a CRS is accurate for my data? #

Compare projected areas against geodesic areas from pyproj.Geod. Under 0.1% is excellent, under 1% is usually fine, and anything above a few percent means the CRS does not fit your area.

Which CRS preserves area? #

Equal-area families — Albers Equal Area, Lambert Azimuthal Equal Area (EPSG:3035 for Europe, 5070 for the contiguous USA), Equal Earth for the world. Use them whenever the analysis compares sizes.

Can I use different CRS for analysis and display? #

You should. Measure in a metric CRS suited to the area, display in whatever the basemap requires, store in whatever recipients expect. Record the reasoning so nobody re-litigates it later.