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. Someone else suggests the first UTM zone they find. A colleague reuses the national grid from their last project, in another country. Every one of those produces numbers in metres, and none of them warns you. Here is Leeds, its 33 wards summed three ways (the setup is under The data used on this page):

import geopandas as gpd

gdf = gpd.read_file("wards.gpkg")
for epsg in (3857, 2157, 27700):
    print(f"EPSG:{epsg:<6} {gdf.to_crs(epsg).geometry.area.sum() / 1e6:>9,.1f} km²")

EPSG:3857     1,581.5 km²
EPSG:2157       554.4 km²
EPSG:27700      551.7 km²

Same data. Web Mercator nearly triples the district. The Irish grid is off by half a percent — too little to notice, more than enough to be wrong. Only one of the three is right, and the rest of this page is how you would know which.

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.

centroid    53.826°N, -1.506°E
extent      0.51° lon × 0.25° lat
            ≈ 34 km × 27 km
utm zone    30N
crosses     0 zone boundary(ies)

Measure the extent, national grid, UTM or a conic

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.826°N, -1.506°E
extent      0.51° lon × 0.25° lat
            ≈ 34 km × 27 km
utm zone    30N
crosses     0 zone boundary(ies)

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:

| Country / region | EPSG | Notes |
|---|---|---|
| Great Britain | 27700 | OSGB36 / British National Grid |
| Northern Ireland | 29903 | Irish Grid |
| Ireland | 2157 | Irish Transverse Mercator |
| France (metropolitan) | 2154 | RGF93 / Lambert-93 |
| Germany | 25832 / 25833 | ETRS89 / UTM 32N, 33N |
| Netherlands | 28992 | Amersfoort / RD New |
| Belgium | 3812 | ETRS89 / Belgian Lambert 2008 |
| Spain | 25830 | ETRS89 / UTM 30N |
| Switzerland | 2056 | CH1903+ / LV95 |
| Contiguous USA | 5070 | NAD83 / Conus Albers (equal-area) |
| Australia | 7855 etc. | GDA2020 / MGA zones |
| New Zealand | 2193 | NZGD2000 / Transverse Mercator |

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.

```python
from pyproj import CRS

crs = CRS.from_epsg(27700)
print(crs.name)
print(crs.area_of_use.name)
print(crs.area_of_use.bounds)

02

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

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))
print(zone_fit(gpd.read_file("nuts2.gpkg")))

✓ entirely in UTM zone 30
✗ spans UTM zones 20–40 (118.9° of longitude) — use a conic projection instead

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
print(f"{len(europe)} NUTS 2 regions, {areas.sum():,.0f} km²")

+proj=aea +lat_1=-4.3709 +lat_2=63.7036 +lat_0=51.6802 +lon_0=14.7862 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
299 NUTS 2 regions, 5,756,490 km²

The extent here is the whole NUTS dataset, which includes the French overseas regions, the Azores and the Canaries, so the centred cone is a poor fit for any of them. For "Europe" in the everyday sense, filter to the continent first — or use the published equal-area CRS for Europe, EPSG:3035, which was designed for exactly this.

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])
print(f"furthest ward edge from the centre of Leeds: {distances.max() / 1000:.1f} km")

furthest ward edge from the centre of Leeds: 10.4 km

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.99928  worst   0.07%
32630                                        median 0.99942  worst   0.06%
3035                                         median 1.00000  worst   0.00%
3857                                         median 2.86321  worst 187.53%
4326                                         median 0.00000  worst 100.00%

Under a tenth of a percent is excellent, under one percent is fine for most purposes, and Web Mercator's error at Leeds' latitude is a different answer rather than a slightly worse one. EPSG:4326 comes out as 100% wrong because its "areas" are in square degrees, a few millionths of the true number. Running this once settles the question for a given study area.

A CRS chooser with reasons, from Leeds to Europe

The data used on this page

Leeds' 33 electoral wards (December 2024) from the ONS, Eurostat's NUTS 2 regions (2024, 1:20 million) from GISCO, Natural Earth's 1:110 million countries, and HM Land Registry's INSPIRE polygons for Leeds — all saved in WGS 84, as layers usually arrive. Land Registry data © Crown copyright and database right 2026, Open Government Licence v3.0; NUTS boundaries © EuroGeographics.

import geopandas as gpd
import pyogrio
import requests

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("wards.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")
gpd.read_file("https://naciscdn.org/naturalearth/110m/cultural/ne_110m_admin_0_countries.zip",
              columns=["NAME", "ISO_A3"]).to_file("world.gpkg")

# HM Land Registry's download sets a cookie and redirects to itself, so use a session
reply = requests.Session().get("https://use-land-property-data.service.gov.uk/datasets/inspire/download/Leeds.zip",
                               timeout=600)
reply.raise_for_status()
with open("leeds_inspire.zip", "wb") as f:
    f.write(reply.content)
pyogrio.set_gdal_config_options({"GML_DOWNLOAD_SCHEMA": "NO"})   # the file's schema URL is unreachable
gpd.read_file("zip://leeds_inspire.zip!Land_Registry_Cadastral_Parcels.gml",
              columns=["INSPIREID"]).to_crs(4326).to_file("parcels.gpkg")
for name in ("wards", "nuts2", "world", "parcels"):
    info = pyogrio.read_info(f"{name}.gpkg")
    print(f"{name + '.gpkg':<14} {info['features']:>8,} features  {info['crs']}")

wards.gpkg           33 features  EPSG:4326
nuts2.gpkg          299 features  EPSG:4326
world.gpkg          177 features  EPSG:4326
parcels.gpkg    293,729 features  EPSG:4326

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,
}

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 [
    ("Leeds 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()

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

── European NUTS2 ────────────────────────────────
chosen  WGS 84 / Equal Earth Greenwich
        EPSG:8857
why     global extent — Equal Earth keeps areas comparable
units   metre
check   area error up to 0.10% on a 100-feature sample

── global countries ──────────────────────────────
chosen  WGS 84 / Equal Earth Greenwich
        EPSG:8857
why     global extent — Equal Earth keeps areas comparable
units   metre
check   area error up to 0.48% 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. Note what it did with the NUTS regions: it reached for a global projection. The NUTS file includes Réunion, French Guiana and the Azores, so its extent is nearly 119° of longitude — the rule did its job on the extent it was given, and the extent was not "Europe". A chooser is only as good as the extent you hand it. (The conic branch returns a custom PROJ string, which has no EPSG code and reports its name as "unknown" — expected, 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.000      0.067            True
         WGS 84 / UTM zone 30N 32630  metre       0.063      0.031            True
OSGB36 / British National Grid 27700  metre       0.073      0.037            True
IRENET95 / Irish Transverse Me  2157  metre       0.436      0.217           False
                        WGS 84  4326 degree     100.000     99.999            True
      WGS 84 / Pseudo-Mercator  3857  metre     187.528     69.560            True

The Irish grid row is the instructive one: 0.44% area error and in_area_of_use = False. Nothing raises, the units are metres, and the error is small enough that nobody would ever notice it by eye — which is exactly what makes it dangerous, because it is several times the national grid's and entirely avoidable. Only the comparison against geodesic truth and the area-of-use check reveal it, which is why both belong in any chooser.

Sorting by error makes the shortlist obvious. Here EPSG:3035 is the most accurate on area — it is equal-area by design — 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

# in a project this lives in one module, say crs_config.py, and everything imports it
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.07% 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."
    ),
)

from pathlib import Path


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
Path("out").mkdir(exist_ok=True)
for_storage(parcels).to_file("out/parcels.gpkg", driver="GPKG")
print(f"{len(parcels):,} titles, median {parcels['area_m2'].median():,.0f} m², "
      f"stored as {gpd.read_file('out/parcels.gpkg', rows=1).crs}")

293,729 titles, median 238 m², stored as EPSG:4326

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.

A projection is a scoping decision

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.

Centroid zones, area_of_use and feet-based grids

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