Polygon Area, in depth: the shoelace derivation, exact integers and Pick's theorem, cancellation at map coordinates, holes, self-intersection and area on the Earth

Deriving the shoelace formula

Take a polygon with vertices v0 to v(n-1) in order, with v(n) = v0. For each edge, form the triangle from the origin to the edge's two endpoints. Its signed area is half the 2D cross product x_i y_(i+1) - x_(i+1) y_i: positive if the origin sees the edge turn counterclockwise, negative if clockwise. Summing over all edges, the triangles outside the polygon are counted once positively and once negatively and cancel, and every point inside is covered once net. So:

2A = sum over i of (x_i * y_(i+1) - x_(i+1) * y_i)

A > 0  ->  vertices are counterclockwise
A < 0  ->  vertices are clockwise

The same result follows from Green's theorem: area is the line integral of (x dy - y dx) / 2 around the boundary, and on a straight edge that integral is exactly the cross-product term. The origin can be any point. That freedom is the key to the numerical fix later. The sign is not a nuisance to discard with abs(); it is the orientation of the ring, and several formats depend on it.

Worked example, checked with Pick's theorem

Take the pentagon v0 (1, 1), v1 (7, 2), v2 (8, 6), v3 (4, 9) and v4 (2, 5).

Reversing the vertex order gives -68. Because every coordinate is an integer, 2A is an integer, and with integer arithmetic the result is exact. Store 2A, not A, if you need exactness: halving is the only step that can introduce a fraction.

Pick's theorem gives an independent check for lattice polygons: A = I + B/2 - 1, where I counts interior lattice points and B counts lattice points on the boundary. An edge with offsets (dx, dy) carries gcd(|dx|, |dy|) boundary points, not counting its start. Four of the edges have gcd 1; the edge from (4,9) to (2,5) has offsets (-2, -4) with gcd 2, because it passes through (3,7). So B = 6. Then I = A - B/2 + 1 = 34 - 3 + 1 = 32. Counting lattice points by brute force found 32 interior and 6 boundary points.

Two implementations

Two implementations cover almost every case: exact integers when coordinates are integral, and shifted floating point otherwise.

import math

def twice_area_int(pts):
    """Exact 2A for integer coordinates (Python ints never overflow)."""
    s = 0
    for (x1, y1), (x2, y2) in zip(pts, pts[1:] + pts[:1]):
        s += x1 * y2 - x2 * y1
    return s

def area_float(pts):
    """Signed area, with every vertex translated so that pts[0] is the origin."""
    x0, y0 = pts[0]
    terms = []
    for i in range(1, len(pts) - 1):
        ax, ay = pts[i][0] - x0, pts[i][1] - y0
        bx, by = pts[i + 1][0] - x0, pts[i + 1][1] - y0
        terms.append(ax * by - bx * ay)
    return math.fsum(terms) / 2

Choosing v0 as the origin makes the first and last terms vanish, so the loop is a fan of triangles from v0. In C, Java or Go, integer coordinates up to about 2 * 10^9 keep each cross-product term inside a signed 64-bit integer, but a sum of many terms can still overflow; use 128-bit accumulation or check the range. math.fsum adds the terms with exact rounding, which matters when terms of mixed sign cancel.

Cancellation at map coordinates

Web Mercator and UTM coordinates are metres measured from a distant origin, so a 10-metre parcel has coordinates around 10^6 to 10^7. The raw shoelace multiplies those large numbers, producing products near 10^14, and then subtracts nearly equal products to recover an answer near 10^2. A double carries about 16 significant digits, so much of the precision is gone before the subtraction happens.

Measured on 1,000 random triangles of up to 50 m by 50 m, placed at Web Mercator offsets of about 7.4 million and 12.2 million metres, against exact rational arithmetic on the same floating-point inputs:

The fix costs two subtractions per vertex and improves accuracy by eleven orders of magnitude. The subtractions of nearby coordinates are exact or nearly so, and the products are then of numbers near 10, not 10^7. Do it always; there is no case where the raw form is better.

Holes, rings and orientation conventions

A polygon with holes is an outer ring plus inner rings, and its area is the outer area minus the hole areas. If rings follow an orientation convention, the signed sum does the subtraction by itself: an outer ring from (0,0) to (10,8) counterclockwise gives +80, a 3 by 3 hole clockwise gives -9, and the sum is 71.

The trap is that formats disagree on the convention:

• GeoJSON (RFC 7946) says exterior rings are counterclockwise and holes clockwise, but also says parsers should not reject polygons that break the rule. Data in the wild violates it routinely, especially data written before 2016.
• Esri shapefiles use the opposite: outer rings clockwise, holes counterclockwise.

Because orientation cannot be trusted, robust code uses roles, not signs: take the absolute area of the ring the format marks as exterior, and subtract the absolute area of each ring marked as a hole. Turf.js's area function does exactly this. Only use the sign itself when your own pipeline produced and normalised the rings. In a multipolygon, sum the per-polygon results; overlapping parts are an invalid input and are counted twice.

Self-intersecting rings

On a ring that crosses itself the formula still returns a number, but not the area you want. It returns the sum over the plane of the winding number at each point. The figure eight (0,0), (2,2), (2,0), (0,2) has two equal lobes traversed in opposite directions, so 2A is exactly 0. The bow-tie (0,0), (4,4), (4,0), (0,2) has lobes of +4/3 and -16/3, so the formula reports -4 while the covered area is 20/3.

No arithmetic error flags this. Detect self-intersection before computing area, with a sweep-line intersection test or a geometry library's validity check, and either reject the ring or repair it into simple pieces first. The winding number itself is explained in Point in Polygon.

Area on the Earth

Applying the shoelace formula to longitude and latitude in degrees gives square degrees, and a square degree is not a fixed area. On a sphere of radius 6,371 km, a 1-degree by 1-degree cell covers 12,363.7 square kilometres at the equator, 8,666.2 at 45 degrees north and 6,088.4 at 60 degrees. Projecting to Web Mercator first is worse in a different way: the same cells measure 1.000, 2.036 and 4.124 times their true area, because Mercator scales area by the square of the secant of latitude.

For small regions, project into a local equal-area projection and use the planar formula. For anything larger, use a spherical formula. The line-integral formula of Chamberlain and Duquette (JPL, 2007) is short and is the one Turf.js cites:

import math
R = 6371008.8                        # mean Earth radius in metres

def ring_area_sphere(ring):          # ring: [(lon_deg, lat_deg), ...], not closed
    n, s = len(ring), 0.0
    for i in range(n):
        lon_prev = ring[i - 1][0]
        lon_next = ring[(i + 1) % n][0]
        s += math.radians(lon_next - lon_prev) * math.sin(math.radians(ring[i][1]))
    return -s * R * R / 2            # sign follows ring orientation; check yours

On the 1-degree cell at 60 degrees it returns 6,088.42 square kilometres, matching the closed-form value. Rings that cross the antimeridian need longitudes unwrapped first. A sphere still differs from the WGS84 ellipsoid by a few tenths of a percent; for legal, survey or billing areas, use Karney's ellipsoidal geodesic algorithms from GeographicLib, which pyproj also exposes.

Streaming, parallel and incremental area

Every term of the sum depends on only two consecutive vertices, which gives the formula three useful properties for large data.

• Streaming. A ring with millions of vertices, such as a coastline, can be measured in one pass while it is read, keeping only the first vertex, the previous vertex and a running sum. Close the ring with the last-to-first term at the end.
• Parallel. Split the vertex list into chunks that overlap by one vertex, sum each chunk's terms independently and add the partial sums. Every chunk must translate by the same reference point, for example the ring's first vertex or a bounding-box corner, or the partial sums are not compatible.
• Incremental. Moving vertex i changes only the two terms that touch it, so an editor can update the area in constant time per drag: subtract the old two terms and add the new ones. Inserting a vertex between i and i+1 replaces one term with two. Recompute from scratch periodically, because a long chain of float updates accumulates rounding.

The same per-edge terms, weighted by vertex sums, also give the centroid, so a pipeline that needs both should compute them in the same pass.

Operational guidance

• Store the method with the number. Record whether an area is planar, spherical or ellipsoidal, and in what units, next to the value. Mixed methods in one column are a common source of reconciliation bugs.
• Test with relative tolerances. Compare planar float results with the exact integer or rational answer to a relative tolerance such as 1e-12, and spherical results with a reference library to a tolerance you have measured.
• Count rejected geometry. Log how many rings fail validation and why. A sudden rise usually means an upstream export changed its orientation or closing convention.
• Treat near-zero areas as a signal. Slivers and collapsed rings with area below a threshold are usually digitising errors; flag them rather than silently keeping them.

Failure modes

• Raw coordinates near 10^7. Errors of hundredths of a square metre per small polygon, which add up across millions of parcels. Translate first.
• Trusting ring orientation. A clockwise exterior from a shapefile read as GeoJSON turns a hole into added area. Use ring roles.
• Self-intersecting input. The answer is a winding-weighted sum, possibly zero. Validate before measuring.
• Degrees treated as planar. Square degrees times the equatorial scale overstate area by a factor of about 2 at 60 degrees latitude. Use a spherical or ellipsoidal formula.
• Closed rings counted twice. GeoJSON repeats the first vertex at the end. The formula tolerates it, but code that indexes pts[n-1] as a distinct vertex may not.

Trade-offs

What to do next

1. Replace any raw shoelace in your code with the translated version and fsum, and add the pentagon from this article (2A = 68) as a unit test.
2. Add a Pick's theorem property test on random lattice polygons to catch sign and indexing bugs.
3. Decide the orientation convention for each data source, and compute hole areas by ring role rather than sign.
4. Run a validity check before area on any user-supplied geometry.
5. Audit every place that computes area from lat/lon and move it to a spherical or ellipsoidal method.
6. Read Convex Hull, in depth for the same orientation predicate used for hulls.