Choose a Projected CRS for Analysis in PyQGIS

Data increasingly arrives in WGS 84 longitude and latitude: GPS exports, web APIs, OpenStreetMap extracts, global datasets. That is fine for storage and display and wrong for almost every measurement. A buffer of 0.01 is not a distance, an area in square degrees is not an area, and a nearest-neighbour search in degrees favours north–south neighbours over east–west ones everywhere except the equator. Before buffering, measuring, gridding or interpolating, the data needs to be in a projected CRS — and which one depends on where the data is and what you are measuring.

This recipe belongs to Coordinate Reference Systems in PyQGIS. It decides between a national grid, a UTM zone and an equal-area projection, computes the right choice from a layer's extent, measures the distortion to confirm the choice, and wires the reprojection into a workflow so results come back in the original CRS.

Which projected CRS?Start with the layer extent in WGS 84. If a national or regional grid covers it, such as British National Grid or ETRS89 UTM, use that. Otherwise, if the extent fits within one UTM zone, use that zone. If the extent spans several zones and the analysis is about areas or densities, use an equal-area projection such as LAEA Europe or a continental Albers. For distances over a large extent, use a projection centred on the data, such as an azimuthal equidistant centred on the study area.Let the extent and the measurement decidelayer extentin WGS 84inside a nationalgrid's area of use?yes: national gridfits in oneUTM zone?yes: that UTM zonemeasuringarea or distance?area: equal-areaLAEA, Albersdistance:local AEQD

Prerequisites

  • QGIS 3.40 LTR or newer, or the QGIS 4 series.
  • A layer with a correctly assigned CRS. If it has none, fix that first — see handling missing CRS.
  • A clear statement of what the analysis measures: distances, areas, directions, or just topology.

Get the extent in WGS 84

Every choice below starts from where the data is on the globe. Transform the layer's extent to WGS 84 once and derive everything from it.

from qgis.core import (
    QgsProject, QgsCoordinateReferenceSystem, QgsCoordinateTransform,
)

layer = QgsProject.instance().mapLayersByName("survey_points")[0]
wgs84 = QgsCoordinateReferenceSystem("EPSG:4326")

to_wgs = QgsCoordinateTransform(layer.crs(), wgs84, QgsProject.instance())
ext = to_wgs.transformBoundingBox(layer.extent())
centre = ext.center()

print(f"lon {ext.xMinimum():.3f} to {ext.xMaximum():.3f}, "
      f"lat {ext.yMinimum():.3f} to {ext.yMaximum():.3f}")
print(f"centre {centre.x():.3f}, {centre.y():.3f}; width {ext.width():.2f}° of longitude")

Breakdown: transformBoundingBox densifies the box edges before transforming, so a projected extent maps to a correct geographic box even where edges curve. The width in degrees of longitude is the quantity that decides whether a single UTM zone is enough: zones are six degrees wide, and distortion rises towards their edges.

Prefer the official grid where it applies

National and regional grids — British National Grid, the Dutch RD New, the ETRS89 UTM zones used across Europe, US State Plane zones — are designed for their territory, match the coordinates in official data, and are what your colleagues' layers already use. The CRS's area of use tells you whether the data lies inside it.

CANDIDATES = ["EPSG:27700", "EPSG:28992", "EPSG:25832", "EPSG:25833", "EPSG:2154"]

def covers(crs_id, box):
    crs = QgsCoordinateReferenceSystem(crs_id)
    bounds = crs.bounds()
    return crs.isValid() and bounds.contains(box), crs

usable = [(cid, crs.description()) for cid in CANDIDATES
          for ok, crs in [covers(cid, ext)] if ok]
print("official grids covering the data:", usable or "none")

Breakdown: bounds() returns the CRS's area of use in WGS 84 as recorded in the EPSG database, so a containment test is a quick, authoritative filter. A candidate list keeps the check fast and relevant — keep it to the grids your organisation actually uses. When one applies, use it even if a UTM zone would be marginally less distorted: matching the coordinates of every other dataset in the organisation is worth more than a few parts per million.

Six-degree strips, and what happens at the edgeFour vertical strips labelled zones 29, 30, 31 and 32. A small green dataset sits inside zone 30 and uses EPSG 32630. A wide orange dataset spans zones 30 and 31; using zone 30 for all of it leaves the eastern part up to 1 in 1,000 in scale error. The formula zone equals floor of longitude plus 180 over 6 plus 1 is shown, with EPSG 326 plus zone in the northern hemisphere and 327 plus zone in the southern.One zone is fine; two zones is a compromisezone 29zone 30zone 31zone 32zone from longitudefloor((lon+180)/6)+1north: EPSG:326zzsouth: EPSG:327zzgreen: inside zone 30orange: spans 30 and 31scale error grows towards a zone's edges

Compute the UTM zone

When no official grid covers the data and it spans less than a zone, WGS 84 UTM is a safe, universal choice. The zone and hemisphere follow directly from the centre of the extent.

import math

def utm_crs_for(box):
    lon, lat = box.center().x(), box.center().y()
    zone = int(math.floor((lon + 180) / 6)) + 1
    zone = min(max(zone, 1), 60)
    epsg = (32600 if lat >= 0 else 32700) + zone
    crs = QgsCoordinateReferenceSystem(f"EPSG:{epsg}")
    zones_spanned = (int(math.floor((box.xMaximum() + 180) / 6))
                     - int(math.floor((box.xMinimum() + 180) / 6)) + 1)
    return crs, zones_spanned

utm, spanned = utm_crs_for(ext)
print(utm.authid(), utm.description(), "| zones spanned:", spanned)
if spanned > 1:
    print("extent crosses a zone boundary — consider an equal-area or local projection")

Breakdown: The formula maps longitudes to zones 1 to 60, and the clamp handles the antimeridian edge case of exactly 180°. EPSG codes 32601–32660 are the northern hemisphere zones and 32701–32760 the southern; the hemisphere matters because southern zones use a false northing of 10,000 km so coordinates stay positive. Norway and Svalbard have irregular zones that this simple formula ignores; for data there, use the national ETRS89 zones instead. Counting the zones spanned is the flag for the next section.

Area statistics and large extents

UTM preserves shapes and local distances well, but not area across large regions. When the output is an area — hectares of forest per province, land cover shares across a continent — use an equal-area projection so a square kilometre in the north counts the same as one in the south. When the output is distance over a large area, use a projection centred on the data.

from qgis.core import QgsDistanceArea

EQUAL_AREA = {
    "europe": "EPSG:3035",
    "north_america": "ESRI:102008",
    "world": "EPSG:6933",
}

laea = QgsCoordinateReferenceSystem(EQUAL_AREA["europe"])

centred = QgsCoordinateReferenceSystem.fromProj(
    f"+proj=aeqd +lat_0={centre.y():.4f} +lon_0={centre.x():.4f} "
    "+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

def distortion_report(layer, target):
    da = QgsDistanceArea()
    da.setSourceCrs(layer.crs(), QgsProject.instance().transformContext())
    da.setEllipsoid("EPSG:7030")
    xf = QgsCoordinateTransform(layer.crs(), target, QgsProject.instance())
    ratios = []
    for f in layer.getFeatures():
        g = f.geometry()
        if g.isEmpty() or g.area() == 0:
            continue
        true_area = da.measureArea(g)
        g.transform(xf)
        ratios.append(g.area() / true_area)
    return min(ratios), max(ratios)

parcels = QgsProject.instance().mapLayersByName("land_parcels")[0]
for label, crs in (("UTM", utm), ("LAEA Europe", laea), ("local AEQD", centred)):
    lo, hi = distortion_report(parcels, crs)
    print(f"{label:<12} planar/ellipsoidal area {lo:.5f} – {hi:.5f}")

Breakdown: QgsDistanceArea with an ellipsoid measures areas on the ellipsoid itself, which is the reference any planar area should be compared against. The ratio of planar to ellipsoidal area per feature shows how much a projection inflates or shrinks areas across the dataset: an equal-area projection gives ratios within a few parts per million everywhere, while UTM drifts by a part in a thousand or more far from its central meridian. ESRI: codes are available because PROJ ships the Esri registry alongside EPSG. The azimuthal equidistant projection centred on the data preserves distances from its centre, which suits service-radius and travel-distance work over a region too big for one UTM zone. The planar-versus-ellipsoidal distinction is also discussed in calculating polygon areas.

Reproject for the analysis, not for storageThe source layer in EPSG 4326 is reprojected to a temporary layer in the chosen analysis CRS. Buffers, distances and areas are computed on the temporary layer. Numeric results are joined back to the source by feature id, and any output geometry is reprojected back to the original CRS or delivered in the analysis CRS with the choice documented.Measure in metres, deliver in the source CRSsourceEPSG:4326reprojecttemporary copyanalysebuffers, areasdeliverjoin or reprojectrecord the analysis CRS in the output metadata

Reproject only for the analysis

Reprojecting the stored data is rarely necessary. Reproject a temporary copy, do the metric work there, and bring results back — numbers by joining on id, geometries by reprojecting.

import processing

analysis_crs = next((QgsCoordinateReferenceSystem(cid) for cid, _ in usable), utm)

work = processing.run("native:reprojectlayer", {
    "INPUT": layer, "TARGET_CRS": analysis_crs, "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

buffers = processing.run("native:buffer", {
    "INPUT": work, "DISTANCE": 250, "SEGMENTS": 16, "DISSOLVE": False,
    "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

delivered = processing.run("native:reprojectlayer", {
    "INPUT": buffers, "TARGET_CRS": layer.crs(),
    "OUTPUT": "/data/work/survey_buffers_250m.gpkg",
})["OUTPUT"]
print("buffers computed in", analysis_crs.authid(), "delivered in", layer.crs().authid())

Breakdown: The analysis CRS is the first official grid that covered the data, falling back to the computed UTM zone — the decision logic from this guide in one line. Buffering in metres in a projected CRS and then reprojecting to WGS 84 gives buffers that are correct on the ground and slightly egg-shaped on a lon/lat map, which is the truthful result. Delivering in the source CRS keeps the output compatible with whatever the data owner uses; delivering in the analysis CRS is equally valid if documented. Buffering itself is covered in buffering a geometry.

QGIS version compatibility

QgsCoordinateReferenceSystem.bounds() and transformBoundingBox have been available since QGIS 3.0; fromProj since 3.10. QgsDistanceArea.measureArea has been stable throughout 3.x. Esri codes such as ESRI:102008 are available wherever PROJ 6 or later is bundled, which covers every supported release and the QGIS 4 series.

Troubleshooting

  • Buffers are huge or tiny. The distance was applied in degrees; reproject first.
  • Areas differ between colleagues by a fraction of a percent. Different analysis CRSs; agree one and document it.
  • UTM coordinates are negative or enormous. Wrong hemisphere code, or the data is far outside the zone.
  • bounds() is empty. The CRS has no recorded area of use, typical for custom definitions.
  • An equal-area result looks distorted on the map. Equal-area projections trade shape for area; display in another CRS if needed.

Conclusion

Transform the extent to WGS 84, prefer an official grid whose area of use covers the data, fall back to the computed UTM zone for compact extents, and use an equal-area or data-centred projection when the extent is large or the output is an area. Confirm the choice with a planar-versus-ellipsoidal comparison, reproject a temporary copy for the analysis, and document the CRS you measured in.

Frequently Asked Questions

Is Web Mercator ever suitable for analysis? No. EPSG:3857 inflates scale by a factor of 1/cos(latitude) — about 60% at 50° north — so distances and areas are badly wrong away from the equator.

Can I avoid projecting by measuring on the ellipsoid? For individual lengths and areas, yes, with QgsDistanceArea. Buffers, grids and interpolation still need planar coordinates.

What about data spanning the antimeridian? Use a projection centred on the data, such as a local azimuthal equidistant, so coordinates are continuous across 180°.

Should the project CRS match the analysis CRS? It does not have to. Processing uses layer CRSs; the project CRS only affects display and on-the-fly measurements.