Spatial Epi Public Health GIS

Choosing an Equal-Area Projection for Rate Mapping

Any health measure with area in its denominator — incidence per square kilometre, facility density, population density used as a covariate — is only as correct as the projection it was computed in. This guide, part of Coordinate Reference Systems for Public Health, covers selecting an equal-area projection for a study extent and quantifying what remains.

Problem Context & Constraints

No map projection preserves area, distance and shape simultaneously. An equal-area projection preserves area exactly, at the cost of distorting shape and, away from its standard parallels, distance. A conformal projection preserves local angles and shape and distorts area. Between them sits a large family of compromise projections that preserve nothing exactly.

For area-normalised measures the choice is forced: use an equal-area projection or the denominator is wrong. The error is not small at regional scale — computing tract areas in Web Mercator inflates them by the square of the scale factor, which at 45° north is a factor of two.

The remaining decision is which equal-area projection, and it is decided by extent and by what else the same coordinates must support. Albers Equal Area Conic suits mid-latitude regions wider than tall; Lambert Azimuthal Equal Area suits compact or circular extents and polar regions; a national grid such as EPSG:5070 for the conterminous United States is an Albers instance with published parameters, which is preferable to a hand-rolled one because it is reproducible by name.

Areal Error Across a 900 km Extent Percentage areal error plotted against distance from the projection centre across a nine hundred kilometre extent. Albers Equal Area and Lambert Azimuthal Equal Area both stay flat at zero across the whole range. UTM, a conformal projection, rises to about two percent at the zone edge and then jumps when the zone changes. Web Mercator rises steeply and reaches forty percent at the northern edge of the extent. Only the first two are usable for a denominator. Two of these are usable as a denominator 0% 20% 40% areal error Albers & Lambert azimuthal equal area UTM — discontinuous at the zone seam Web Mercator 0 450 km 900 km distance from the projection centre The flat line is not an approximation — equal-area projections preserve area exactly

Prerequisites

  • The study extent as a bounding box in geographic coordinates
  • python 3.11, pyproj 3.6.1, geopandas 1.0.1, shapely 2.0.6, numpy 1.26.4
  • A decision about whether the same coordinates must also support distance work, since equal-area projections distort distance away from their standard parallels

The family choice follows from the shape of the extent rather than from preference:

Which Projection Family Suits Which Extent Three extent shapes matched to projection families. A wide mid-latitude extent, such as the conterminous United States, suits Albers Equal Area Conic with two standard parallels. A compact or circular extent, such as a metropolitan region, suits Lambert Azimuthal Equal Area centred on it. A tall narrow extent, such as a single state running north to south, is served by either, with Albers preferred when a published national grid already covers it. Extent shape picks the family Albers conic wide, mid-latitude Lambert azimuthal compact or polar either — prefer a grid tall and narrow A published authority code beats a hand-derived one whenever one covers the extent

Step-by-Step Solution

# Select and validate an equal-area projection for a study extent.
# Pinned: pyproj==3.6.1, geopandas==1.0.1, shapely==2.0.6, numpy==1.26.4
import logging
import numpy as np
import geopandas as gpd
from pyproj import CRS, Geod, Transformer

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("crs.equalarea")

GEOD = Geod(ellps="WGS84")

# Published national grids are preferred over hand-rolled parameters: they are
# reproducible from an authority code and reviewers can look them up.
KNOWN_EQUAL_AREA = {
    "conus": "EPSG:5070",     # NAD83 / Conus Albers
    "alaska": "EPSG:3338",    # NAD83 / Alaska Albers
    "europe": "EPSG:3035",    # ETRS89-extended / LAEA Europe
    "world": "EPSG:6933",     # WGS 84 / NSIDC EASE-Grid 2.0 Global
}


def suggest_albers(bounds) -> CRS:
    """Albers Equal Area with standard parallels at 1/6 and 5/6 of the latitude span.

    This is the classic rule of thumb and it minimises the (already zero) areal
    error's companion distortions across the extent. Use a published grid instead
    whenever one covers the study area."""
    minx, miny, maxx, maxy = bounds
    span = maxy - miny
    lat1, lat2 = miny + span / 6.0, maxy - span / 6.0
    lon0, lat0 = (minx + maxx) / 2.0, (miny + maxy) / 2.0
    crs = CRS.from_proj4(
        f"+proj=aea +lat_1={lat1:.4f} +lat_2={lat2:.4f} +lat_0={lat0:.4f} "
        f"+lon_0={lon0:.4f} +datum=WGS84 +units=m +no_defs")
    log.info("suggested Albers: lat_1=%.3f lat_2=%.3f lat_0=%.3f lon_0=%.3f",
             lat1, lat2, lat0, lon0)
    return crs


def areal_error(gdf: gpd.GeoDataFrame, crs) -> dict:
    """Projected polygon area against geodesic area on the ellipsoid.

    The geodesic area is the truth. An equal-area projection should agree with it
    to within floating-point noise; anything above about 0.01% means the CRS is
    not equal-area, whatever its name says."""
    proj = gdf.to_crs(crs)
    geo = gdf.to_crs("EPSG:4326")
    truth = np.array([abs(GEOD.geometry_area_perimeter(g)[0]) for g in geo.geometry])
    got = proj.geometry.area.to_numpy()
    rel = np.abs(got - truth) / truth
    res = {"median_pct": float(100 * np.median(rel)), "max_pct": float(100 * rel.max())}
    log.info("areal error vs geodesic: median %.5f%%, max %.5f%%",
             res["median_pct"], res["max_pct"])
    return res


def assert_equal_area(gdf: gpd.GeoDataFrame, crs, tol_pct: float = 0.01) -> None:
    """Gate: refuse to proceed if the CRS does not actually preserve area."""
    e = areal_error(gdf, crs)
    if e["max_pct"] > tol_pct:
        raise ValueError(
            f"CRS {CRS.from_user_input(crs).to_string()[:60]} has {e['max_pct']:.3f}% "
            f"maximum areal error over this extent; it is not equal-area here")
    log.info("equal-area gate passed")

The gate is the useful part. A CRS is often assumed to be equal-area because it was chosen for that reason years earlier, and testing against geodesic area takes milliseconds and settles the question.

Validation & Edge Cases

1. Test against geodesic area, not against another projection. Comparing two projections tells you they disagree, not which is right. pyproj’s geodesic area on the ellipsoid is the reference:

INFO suggested Albers: lat_1=34.833 lat_2=41.167 lat_0=38.000 lon_0=-79.500
INFO areal error vs geodesic: median 0.00002%, max 0.00009%
INFO equal-area gate passed
INFO areal error vs geodesic: median 26.41%, max 38.77%   [EPSG:3857 — rejected]

2. Do not use an equal-area projection for distance work without checking. Albers distorts distance away from its standard parallels, by a few tenths of a percent within a state and more across a continent. Where the same pipeline needs both, either use two projections explicitly or compute distances geodesically.

3. Prefer a published authority code. EPSG:5070 is reproducible from four characters; a proj4 string with hand-computed parallels is reproducible only if it is stored verbatim, and it will differ between analysts who each applied the rule of thumb to slightly different bounds.

4. Re-derive the parameters if the extent changes materially. A projection tuned for one state is fine for a neighbouring one and poor for the far side of the country, and the standard-parallel rule assumes the extent it was computed from.

5. Watch datum, not only projection. An equal-area projection on the wrong datum still misplaces everything by up to a few hundred metres; the two decisions are independent and both belong in the authority code.

6. Do not confuse an equal-area projection with an equidistant one. The names are similar and the properties are unrelated. An analysis that buffers case points by five kilometres in an Albers projection is measuring a distance in a CRS that does not preserve distance, and the error grows away from the standard parallels. Buffers, nearest-neighbour distances and drive-time snapping all belong in a conformal or equidistant projection, or should be computed geodesically.

7. Watch for area calculations hidden inside other operations. Population density covariates, dasymetric weights, polygon-overlap fractions in an areal interpolation and cell areas in a raster zonal statistic all consume area, and none of them announces it. Auditing a pipeline for area dependence usually turns up two or three operations nobody had classified that way, and each of them inherits whatever CRS was current at the time.

8. Set the analysis CRS once, at ingestion. Reprojecting repeatedly through a pipeline accumulates both floating-point drift and opportunities to forget, and a single documented reprojection at the boundary is easier to audit than six scattered ones.

9. Re-run the gate whenever the boundary file changes. A new vintage can extend the study extent beyond the range the projection was tuned for, and the areal-error check is the cheapest way to find out before the denominators shift.

Auditing a pipeline for area dependence is more productive than it sounds, because most of the dependence is implicit:

Operations That Silently Consume Area Five pipeline operations that depend on area without saying so: a population density covariate, dasymetric weighting, areal interpolation overlap fractions, raster zonal statistics, and any rate expressed per square kilometre. Each inherits whatever coordinate reference is current when it runs, and none of them raises an error in a projection that does not preserve area. Five places area enters without announcing itself population density as a covariate dasymetric weighting habitable land area overlap fractions areal interpolation zonal statistics raster cell area rates per km² the obvious one Only the last one is usually classified as an area calculation audit a pipeline for the other four and expect to find at least two nobody had listed All five inherit whatever CRS happened to be current when they ran

Compliance Notes

  • Record the CRS as an authority code wherever one exists, and store the full definition where one does not.
  • Log the areal-error gate result with the run, since it is the evidence that the denominator is sound.
  • State the projection in the map’s metadata, because a reader recomputing an area from the published geometry needs to know what to reproject to.
  • Keep the display projection separate from the analysis projection. A map may be drawn in Web Mercator for a web viewer as long as every number on it was computed elsewhere, and that separation should be explicit in the pipeline rather than implicit.