Measuring Area Accurately on OSM Polygons

Area is the measurement people get wrong most often and notice least, because every wrong answer is a plausible number. A forest measured in Web Mercator at sixty degrees north is four times too large, and nothing about the figure looks unusual.

Prerequisites

Conceptual minimum

OSM stores geographic coordinates, and degrees are not a unit of distance. A square degree near the equator covers roughly twelve thousand square kilometres; at sixty degrees north it covers about six thousand, because lines of longitude converge. Computing area from degree coordinates therefore produces a number whose unit varies with latitude.

Web Mercator is worse, not better. It is a conformal projection, preserving angles at the cost of area, and its area distortion grows as the square of the secant of the latitude. At sixty degrees it inflates area by a factor of four; at seventy-five, by about fifteen. Because Web Mercator is the projection everything renders in, it is the one people reach for, and it is the least suitable for measurement of any in common use.

Two approaches give correct answers.

An equal-area projection — a local one such as a Lambert Azimuthal Equal Area centred on the feature, or a regional standard — preserves area by construction. Accuracy is excellent for features within a few hundred kilometres of the projection centre.

A geodesic computation measures on the ellipsoid directly, with no projection at all. It is correct at any size and any location, marginally slower, and the right default when features are scattered.

How much Web Mercator inflates area, by latitude Five latitudes with the factor by which a Web Mercator area calculation overstates the true area. At the equator the factor is one, so the result is correct. At thirty degrees it is about one point three. At forty five degrees it is two. At sixty degrees it is four. At seventy five degrees it is about fifteen. A note observes that the inflation grows as the square of the secant of the latitude, so it accelerates sharply towards the poles. Web Mercator area inflation, by latitude 0 degrees correct 30 degrees 1.3x too large 45 degrees 2x too large 60 degrees 4x too large 75 degrees 15x too large Inflation grows as the square of the secant of latitude, so a dataset spanning Europe has errors varying by over a factor of two.
A single correction factor cannot fix this, because the error differs across any dataset large enough to matter.

Runnable solution

python 
from __future__ import annotations

import logging
import math
from dataclasses import dataclass

from pyproj import Geod, Transformer
from shapely.geometry import Polygon
from shapely.geometry.base import BaseGeometry
from shapely.ops import transform

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
logger = logging.getLogger("osm.crs.area")

GEOD = Geod(ellps="WGS84")
LOCAL_LIMIT_KM = 500.0      # beyond this, a local projection loses accuracy


@dataclass(frozen=True)
class AreaResult:
    square_metres: float
    method: str
    note: str = ""


def geodesic_area(geom: BaseGeometry) -> AreaResult:
    """Measure on the ellipsoid. Correct at any size, anywhere."""
    if geom.geom_type == "MultiPolygon":
        total = sum(geodesic_area(p).square_metres for p in geom.geoms)
        return AreaResult(total, "geodesic", "summed over parts")
    if geom.geom_type != "Polygon":
        raise TypeError(f"area is undefined for {geom.geom_type}")

    lons, lats = zip(*geom.exterior.coords)
    area, _perimeter = GEOD.polygon_area_perimeter(lons, lats)
    total = abs(area)
    for ring in geom.interiors:
        hole_lons, hole_lats = zip(*ring.coords)
        hole_area, _ = GEOD.polygon_area_perimeter(hole_lons, hole_lats)
        total -= abs(hole_area)
    return AreaResult(total, "geodesic")


def equal_area_projection_area(geom: BaseGeometry) -> AreaResult:
    """Project onto an equal-area plane centred on the feature, then measure."""
    centre = geom.centroid
    # Lambert Azimuthal Equal Area centred on the feature itself: distortion
    # is zero at the centre and grows slowly with distance from it.
    proj = (f"+proj=laea +lat_0={centre.y} +lon_0={centre.x} "
            f"+x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs")
    transformer = Transformer.from_crs("EPSG:4326", proj, always_xy=True)
    projected = transform(transformer.transform, geom)
    return AreaResult(projected.area, "equal-area projection")


def extent_km(geom: BaseGeometry) -> float:
    minx, miny, maxx, maxy = geom.bounds
    _, _, diagonal = GEOD.inv(minx, miny, maxx, maxy)
    return diagonal / 1000.0


def area_of(geom: BaseGeometry) -> AreaResult:
    """Pick a method from the feature's extent, and say which was used."""
    if geom.is_empty:
        return AreaResult(0.0, "none", "empty geometry")
    if not geom.is_valid:
        # An invalid ring has no well-defined area; do not guess one.
        raise ValueError("geometry is invalid; repair before measuring")

    span = extent_km(geom)
    if span > LOCAL_LIMIT_KM:
        result = geodesic_area(geom)
        return AreaResult(result.square_metres, result.method,
                          f"extent {span:.0f} km exceeds the local-projection limit")
    return equal_area_projection_area(geom)


def mercator_error_factor(latitude_deg: float) -> float:
    """How much a Web Mercator area calculation would overstate, here."""
    return 1.0 / (math.cos(math.radians(latitude_deg)) ** 2)


if __name__ == "__main__":
    square = Polygon([(19.9, 50.0), (20.0, 50.0), (20.0, 50.1), (19.9, 50.1)])
    for result in (area_of(square), geodesic_area(square)):
        logger.info("%-26s %12.1f m2 %s", result.method,
                    result.square_metres, result.note)
    logger.info("Web Mercator here would overstate by %.2fx",
                mercator_error_factor(50.0))

Step-by-step walkthrough

  1. Refuse to measure invalid geometry. A self-intersecting ring has no well-defined area, and the number a library returns for one is an artefact of the algorithm rather than a measurement.
  2. Subtract holes explicitly in the geodesic path. The ellipsoid area function measures a single ring, so interior rings must be computed and subtracted rather than assumed handled.
  3. Take absolute values. The sign of a geodesic ring area depends on winding direction, and OSM ring orientation is not guaranteed.
  4. Centre the projection on the feature. An equal-area projection has zero distortion at its centre, so centring it per feature keeps accuracy high without needing a regional standard.
  5. Switch on extent, not on location. A feature spanning more than a few hundred kilometres accumulates error in any local projection, and the geodesic path has no such limit.
  6. Report the method. A stored area whose method is unknown cannot be compared against one computed differently, and the two can differ by a fraction of a percent, which matters when the numbers are summed.
  7. Keep the Mercator factor available. Being able to say by how much a wrong calculation would have been wrong is what convinces somebody that the existing figures need recomputing.
Four ways to compute area, and what each is good for A grid of four methods against their correctness and their appropriate use. Computing in degrees gives a number whose unit varies with latitude and is never correct. Computing in Web Mercator gives a number inflated by a latitude-dependent factor and is never correct for measurement. A local equal-area projection is accurate for features within a few hundred kilometres of its centre and is fast. A geodesic computation on the ellipsoid is accurate at any size and location and is the right default when features are scattered. Four methods, two of which are never right Correct? Use it for Degrees never nothing Web Mercator never rendering only Local equal-area within ~500 km most features Geodesic always large or scattered Web Mercator appears in the table only because it is what most pipelines are already using without having chosen to.
The bottom two rows agree to well within a percent on ordinary features, so either is a defensible default.
The decision each area computation makes, in order Four steps. The validate step refuses invalid geometry, since a self-intersecting ring has no well-defined area and any number returned is an artefact. The size step measures the feature's diagonal extent, which is what decides whether a local projection remains accurate. The method step picks a feature-centred equal-area projection for ordinary features and a geodesic computation for large or scattered ones. The record step stores which method produced the figure, so areas computed at different times remain comparable. Validate, size, choose, record validate refuse invalid rings no area is defined size diagonal extent decides the method choose projection or geodesic both are correct record store the method keeps figures comparable Recording the method sounds pedantic until two datasets computed differently are summed and the total is defensible to nobody.
The first step is the one most often skipped, and it is the only one that can turn a wrong number into an error.

Verification

  • A known area matches. Take a feature whose area is published independently and confirm agreement within a fraction of a percent.
  • The two correct methods agree. For a moderate feature, the local projection and the geodesic computation should differ negligibly.
  • Latitude does not change the answer. Move a test polygon of fixed geodesic size to several latitudes; the computed area should stay constant.
  • Holes are subtracted. A polygon with a hole must measure less than the same outline without it.
  • Invalid geometry raises. Feed a bowtie and confirm an error rather than a number.

Common errors and fixes

Symptom Root cause One-line fix
Areas grow with latitude Computed in Web Mercator Use an equal-area projection or a geodesic method
Areas in the thousandths Computed in degrees Project or measure geodetically before reporting
Holes not subtracted Only the exterior ring measured Compute and subtract each interior ring
Negative areas Ring winding taken as meaningful Take absolute values; OSM winding is not guaranteed
Large features slightly off Local projection used beyond its range Switch to geodesic above a few hundred kilometres
Stored areas incomparable Method not recorded Record which method produced each figure
Nonsense area on a bowtie Invalid geometry measured anyway Validate before measuring and refuse invalid input

Specification reference

Web Mercator preserves angles and inflates area by a factor equal to the square of the secant of the latitude, which reaches four at sixty degrees. Equal-area projections such as Lambert Azimuthal Equal Area preserve area exactly at the cost of shape, and geodesic polygon area computation on a reference ellipsoid is independent of any projection. See the pyproj Geod documentation for the ellipsoidal computation and Coordinate Reference Systems in OSM for the projection background.

Frequently Asked Questions

Why is Web Mercator so bad for area?

Because it is conformal: it preserves local angles, and the price of that is that scale varies with latitude. Area, being scale squared, is inflated by the square of the secant of the latitude — a factor of two at forty-five degrees and four at sixty. It is also the projection everything renders in, so it is the one already loaded in most pipelines, which is exactly why the mistake is so common.

Should I use an equal-area projection or a geodesic computation?

Either is correct for ordinary features, and they agree to well within a percent. A local equal-area projection centred on the feature is slightly faster and entirely adequate up to a few hundred kilometres of extent. A geodesic computation has no size limit and needs no per-feature projection setup, which makes it the simpler default when features are scattered across a continent.

Can I just apply a correction factor to Mercator areas?

Only for a single point, which is rarely what you have. The factor depends on latitude, so a polygon spanning several degrees has a different error at its north and south edges, and a dataset spanning a country has errors differing by tens of percent between its extremes. A single correction is exact nowhere except at the latitude it was derived for.

Does ring winding direction matter?

For the geodesic computation, yes, in that it determines the sign of the result — which is why absolute values are taken. OSM does not guarantee a winding direction on multipolygon rings, and the assembly discussed in Handling Multipolygon Members with No Role derives containment geometrically rather than from orientation, so the sign carries no information worth preserving.

Up one level: Coordinate Reference Systems in OSM.