Choosing an Equal-Area CRS for Forest Area Statistics in Python

Every forest report contains a number in hectares, and that number depends entirely on the projection it was measured in. Compute stand area in Web Mercator at 60° north and you will over-report by a factor of four; compute it in a UTM zone 800 km from the central meridian and you will be off by a fraction of a percent that nobody notices but an auditor. This guide is about picking the projection that makes an area figure defensible, as one decision inside Coordinate Reference Systems for Forestry, part of Ecological GIS Data Foundations in Python.

When to use each projection family

No projection preserves area, distance, shape and direction at once. An area statistic needs an equal-area projection, and the choice within that family is about the shape of the region.

The practical rule for forestry: use a national or regional equal-area projection for any reported area, and use UTM for local geometry work where distances and shapes matter more than a hectare in ten thousand. Never compute area in Web Mercator, and never compute it from unprojected degrees.

Minimal reproducible example

Two lines get the area right; the third proves it. geopandas computes planar area in whatever CRS the frame happens to be in, which is exactly why the reprojection has to be explicit.

import geopandas as gpd


def stand_areas_ha(stands_path: str, equal_area_epsg: int = 3035) -> gpd.GeoDataFrame:
    """Stand areas in hectares, computed in a stated equal-area CRS.

    3035 = ETRS89 / LAEA Europe. Use 5070 for CONUS Albers, 3577 for Australian Albers,
    or a custom Albers for a region that none of the standard codes fits well.
    """
    stands = gpd.read_file(stands_path)
    if stands.crs is None:
        raise ValueError("layer has no CRS — refusing to guess")
    ea = stands.to_crs(epsg=equal_area_epsg)
    out = stands.copy()
    out["area_ha"] = ea.area / 10_000.0
    out["area_crs"] = f"EPSG:{equal_area_epsg}"
    return out

Storing area_crs alongside the number is the difference between a figure and a defensible figure. Two reports on the same ownership that disagree by 3% are almost always two projections, and without that column the disagreement takes a day to resolve.

Where no standard code fits — a region straddling zones, or one whose national grid is conformal rather than equal-area — build a custom Albers centred on the data:

import geopandas as gpd
from pyproj import CRS


def local_albers(gdf: gpd.GeoDataFrame) -> CRS:
    """An Albers Equal Area CRS fitted to this layer's own extent."""
    minx, miny, maxx, maxy = gdf.to_crs(4326).total_bounds
    lat0 = (miny + maxy) / 2.0
    lon0 = (minx + maxx) / 2.0
    # Standard parallels at 1/6 and 5/6 of the latitude span minimise distortion.
    span = maxy - miny
    lat1, lat2 = miny + span / 6.0, maxy - span / 6.0
    return CRS.from_proj4(
        f"+proj=aea +lat_1={lat1:.4f} +lat_2={lat2:.4f} "
        f"+lat_0={lat0:.4f} +lon_0={lon0:.4f} "
        "+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
    )

The one-sixth rule for standard parallels is the standard heuristic and holds distortion under 0.1% across a region a few hundred kilometres tall. Record the resulting PROJ string with the output; a custom CRS that exists only in a script is not reproducible.

Parameter reference

Expected output and verification

The check that settles arguments is to compute the same area three ways and compare. Geodesic area — computed on the ellipsoid itself, with no projection — is the reference truth:

import geopandas as gpd
from pyproj import Geod


def compare_area_crs(stands_path: str, candidates: dict[str, int]) -> list[dict]:
    """Compare projected areas against the geodesic (ellipsoidal) truth."""
    stands = gpd.read_file(stands_path)
    geod = Geod(ellps="WGS84")
    truth_ha = sum(
        abs(geod.geometry_area_perimeter(geom)[0]) for geom in stands.to_crs(4326).geometry
    ) / 10_000.0

    rows = [{"crs": "geodesic (truth)", "area_ha": round(truth_ha, 2), "error_pct": 0.0}]
    for name, epsg in candidates.items():
        ha = stands.to_crs(epsg=epsg).area.sum() / 10_000.0
        rows.append({
            "crs": f"{name} (EPSG:{epsg})",
            "area_ha": round(ha, 2),
            "error_pct": round(100.0 * (ha - truth_ha) / truth_ha, 3),
        })
    return rows


for row in compare_area_crs("stands.gpkg", {
    "LAEA Europe": 3035, "UTM 33N": 32633, "Web Mercator": 3857,
}):
    print(row)

Expect the equal-area result within about 0.1% of geodesic, UTM within a few tenths of a percent, and Web Mercator wrong by tens of percent — the exact figure depending on latitude. Turn that into an assertion so a projection mistake fails the pipeline rather than the audit:

def assert_area_crs_sane(gdf, epsg: int, tolerance_pct: float = 0.5) -> None:
    """Refuse to report areas from a CRS that distorts them beyond tolerance."""
    geod = Geod(ellps="WGS84")
    truth = sum(abs(geod.geometry_area_perimeter(g)[0])
                for g in gdf.to_crs(4326).geometry)
    got = gdf.to_crs(epsg=epsg).area.sum()
    err = 100.0 * (got - truth) / truth
    assert abs(err) <= tolerance_pct, (
        f"EPSG:{epsg} distorts area by {err:.2f}% over this extent — "
        "use an equal-area CRS fitted to the region")

Common pitfalls

• Computing area on an unprojected frame. gdf.area on EPSG:4326 returns square degrees. It does not warn, and the number looks like a plausible small figure.
• Reporting from Web Mercator because that is what the basemap used. Display CRS and analysis CRS are independent choices; keep the map in 3857 and the statistics in an equal-area frame.
• Using one UTM zone across a region that spans three. Stands in the far zones carry the largest error and they are exactly the ones nobody checks.
• Mixing datums silently. Reprojecting NAD27 data as though it were NAD83 shifts geometry by up to 200 m. Area survives, position does not — and the overlay against a raster then samples the wrong pixels.
• Not recording the CRS with the number. Every area figure should travel with the EPSG code or PROJ string that produced it.

Frequently Asked Questions

Related

• Coordinate Reference Systems for Forestry — the wider CRS discipline this decision belongs to
• How to Fix CRS Mismatches in GeoPandas — repairing layers that disagree before you measure anything
• Zonal Statistics for Forest Inventory with rasterstats — where per-hectare figures get multiplied up
• Stratified Random Sampling for Forest Plots — stratum areas are weights, so their projection matters

Up: Coordinate Reference Systems for Forestry · Ecological GIS Data Foundations in Python