# Comparing Coordinate Projection Methods For Polygon Area Calculations

## Why Area Measurements Change with Coordinate Systems

A map projection is a systematic transformation of the latitudes and longitudes of locations on the surface of a sphere or ellipsoid into locations on a plane. This process inevitably introduces distortions in shape, area, distance and direction as the three-dimensional Earth is transformed into a two-dimensional map.

Different map projections use different mathematical formulas and techniques to project locations from the Earth’s surface to a planar coordinate system. Some projections aim to preserve shape, others aim to preserve area, while some try to find a compromise between multiple types of distortions. Each projection has its specific use cases where it is more appropriate to use.

The change in shape of spatial features when changing coordinate systems consequently also affects calculations of feature areas. The same polygon representing an area on the Earth will have different planar area measurements when that polygon is defined in different projected coordinate systems. This is because the projection process stretches, compresses and deforms the shape of the polygon to varying degrees.

For example, a polygon enclosing North America has an area of about 24,256,000 km2 when calculated geodesically based on latitude and longitude on the ellipsoid. That same polygon has a measured planar area of 23,839,500 km2 in the Equal Area Cylindrical projection, which preserves area at the expense of shape. In the Mercator projection, which preserves shape but greatly distorts area near the poles, the North America polygon only has a calculated planar area of 15,564,000 km2.

## Methods to Calculate Area

There are two main methods used to calculate the area of polygons on maps – planar area and geodesic area.

### Planar Area

Planar area calculates the area of the polygon as it is projected onto a two-dimensional surface, based on the Cartesian coordinates in that projection. It effectively treats the map as a flat plane rather than part of a sphere/ellipsoid. Area formulas like triangulation or the Shoelace formula can be used. Planar area is easier to calculate but does not give an accurate, real-world area measurement for most projections.

### Geodesic Area

Geodesic area takes the actual curved surface of the Earth into account. The area integrals are done over latitude-longitude geodesics rather than on a projected grid. This provides an accurate, real-world area but requires more complex math and computations. A spherical or ellipsoidal Earth model must be chosen. Geodesic area tends to be more computationally expensive to calculate compared to planar area.

In summary:

- Planar area
- Conceptually simpler
- Computationally faster
- Distorts area in most projections

- Geodesic area
- Models Earth’s surface
- Provides real-world area
- More computationally intensive

## Sample Polygon Area Analysis

As an example, we will analyze a polygon enclosing part of Southern Africa and calculate its area using different projections and methods.

### Polygon Definition

The polygon vertices are defined below based on the WGS 84 ellipsoid in latitude and longitude degrees:

- Point 1: -28, 12
- Point 2: -28, 40
- Point 3: -35, 40
- Point 4: -35, 12

### Projections Analyzed

Three projections are selected to demonstrate the range of possible planar areas for this polygon:

- Mercator
- Robinson
- Equal Earth

### Planar Area Results

The planar area is calculated using the Shoelace formula for each projection coordinate space. Results:

- Mercator Projection: 2,104,000 km2
- Robinson Projection: 2,276,000 km2
- Equal Earth: 2,162,000 km2

The different projections stretch the polygon to different degrees resulting in over 8% difference between extrema. The Equal Earth provides a compromise projection for calculating approximate polygon area.

### Geodesic Area Calculation

To calculate the most accurate, real-world area, geodesic methods are used with an oblate spheroid model approximating the Earth’s surface. The geodesic area is calculated to be 2,201,891 km2.

As the geodesic area takes the curvature of the Earth into account, it should be considered as the authoritative measure for the area enclosed by this polygon.

## Recommendations for Selecting Projection and Area Method

### When to Use Planar Area

Planar area is appropriate for simple shape and area comparisons as long as the distortion effects of the projection are small enough to be acceptable. It can give a quick, easy area indication in many cases. Planar area may also be preferred in cases where projecting and flattening complex geodesic polygons would be computationally prohibitive.

### When to Use Geodesic Area

For authoritative area figures, geodesic area should be used wherever possible. It gives accurate, real-world values representing the actual surface area on the Earth’s surface. Geodesic calculations are essential for any geospatial analytics requiring high accuracy analyses over large regions and globally significant calculations.

### Projection Selection Factors

Projection choice depends on the geospatial extent, required accuracy, expected distortions and purpose of the analysis work. Conformal projections like Mercator preserve shape well and are suitable for small scale mapping and analysis for mid-latitude regions. Equal-area projections such as Eckert IV give accurate area representations but distort shape. Compromise projections provide balance. Additional factors like developing appropriate tessellation for topology and data formats should also be considered.

## Example Code for Polygon Area Analysis

Below is Python code that loads polygon vertex coordinates, projects to different coordinate reference systems, calculates planar and geodesic areas, and compares the values.

“`python

import pyproj

from math import radians, cos, sin, asin, sqrt

# Polygon points

points = [(-28, 12), (-28, 40), (-35, 40), (-35, 12)]

# Function to calculate geodesic area

def geodesic_area(points):

# Convert degree coordinates to radians

rad_points = [(radians(lat), radians(lon)) for lat, lon in points]

# WGS84 ellipsoid axes

a = 6378137.0

b = 6356752.3142

area = 0.0

# Apply geodesic formula

for i in range(-1, len(rad_points)-1):

p1 = rad_points[i]

p2 = rad_points[i+1]

x1 = a * cos(p1[0]) * cos(p1[1])

y1 = a * cos(p1[0]) * sin(p1[1])

z1 = b * sin(p1[0])

x2 = a * cos(p2[0]) * cos(p2[1])

y2 = a * cos(p2[0]) * sin(p2[1])

z2 = b * sin(p2[0])

crossprod = x1*y2 – x2*y1

area += sqrt(x2*x2 + y2*y2 + z2*z2) * atan2(crossprod, x1*x2 + y1*y2 + z1*z2)

# Adjust for duplicate endpoints

area = area * 0.5

return abs(area) # Return positive value

# Shoelace algorithm for planar area

def shoelace_area(points):

area = 0.0

for i in range(-1, len(points)-1):

p1 = points[i]

p2 = points[i+1]

area += (p2[0] + p1[0]) * (p2[1] – p1[1])

# Return absolute value

return abs(area / 2)

# Define projections

merc = pyproj.Proj(“+proj=merc +a=6378137 +b=6378137”)

robin = pyproj.Proj(“+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs”)

ee = pyproj.Proj(“+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs”)

# Convert points

merc_pts = [merc(*pt) for pt in points]

robin_pts = [robin(*pt) for pt in points]

ee_pts = [ee(*pt) for pt in points]

# Calculate areas

merc_area = shoelace_area(merc_pts)

robin_area = shoelace_area(robin_pts)

ee_area = shoelace_area(ee_pts)

geod_area = geodesic_area(points)

# Print areas

print(f”Mercator Projection: {merc_area / 1e6:.3f} million km^2″)

print(f”Robinson Projection {robin_area / 1e6:.3f} million km^2″)

print(f”Equal Earth: {ee_area / 1e6:.3f} million km^2″)

print(f”Geodesic: {geod_area / 1e6:.3f} million km^2″)

“`

The projection transformations and planar/geodesic area functions would be defined elsewhere. The main application logic loads the vertex data, defines projections, converts coordinates, calculates areas with different methods, and prints a comparison report.

## Conclusion and Summary

When doing geographic area analysis, it is critical to select an appropriate projection and use the proper planar or geodesic calculations to get accurate results. Planar area is simpler but can be highly distorted depending on projection properties and geographic extent. Geodesic area gives real-world values but requires complex spherical calculations.

Understanding the effects of map projections and topological transformations can help analysts select effective methods and coordinate reference systems for spatial analysis like area computation. Matching projection properties to analysis needs and computational constraints allows maximizing accuracy while controlling precision and performance.