Compute the area and perimeter of a polygon on the WGS84 ellipsoid using geodesic calculations. This gives accurate results for polygons of any size, including those spanning large portions of the globe.

polygon_area(x, id = NULL, polyline = FALSE)

Arguments

x

A two-column matrix or data frame of coordinates (longitude, latitude) in decimal degrees defining polygon vertices, or a list with longitude and latitude components.

id

Optional integer vector identifying separate polygons. Points with the same id are treated as vertices of the same polygon. If NULL (default), all points are treated as a single polygon.

polyline

Logical. If FALSE (default), compute area and perimeter of a closed polygon. If TRUE, compute only the length of a polyline (area will be meaningless).

Value

  • For a single polygon (id = NULL): A list with components:

    • area: Signed area in square meters. Positive for counter-clockwise polygons, negative for clockwise.

    • perimeter: Perimeter in meters.

    • n: Number of vertices.

  • For multiple polygons (id specified): A data frame with columns:

    • id: Polygon identifier

    • area: Signed area in square meters

    • perimeter: Perimeter in meters

    • n: Number of vertices

Details

The polygon area is computed using the geodesic method which accounts for the ellipsoidal shape of the Earth. This is more accurate than spherical approximations, especially for large polygons.

The area is signed: counter-clockwise polygons have positive area, clockwise polygons have negative area. The absolute value gives the actual area.

For very large polygons (more than half the Earth's surface), the sign convention may seem counterintuitive - the "inside" is the smaller region.

The computation uses the WGS84 ellipsoid (the same as GPS).

Examples

# Triangle: London - New York - Rio de Janeiro
pts <- cbind(
  lon = c(0, -74, -43),
  lat = c(52, 41, -23)
)
polygon_area(pts)
#> $area
#> [1] 2.653936e+13
#> 
#> $perimeter
#> [1] 22634340
#> 
#> $n
#> [1] 3
#> 

# Multiple polygons using id
pts <- cbind(
  lon = c(0, -74, -43, 100, 110, 105),
  lat = c(52, 41, -23, 10, 10, 20)
)
polygon_area(pts, id = c(1, 1, 1, 2, 2, 2))
#>   id         area perimeter n
#> 1  1 2.653936e+13  22634340 3
#> 2  2 6.061691e+11   3556143 3

# Polyline length (not a closed polygon)
polygon_area(pts[1:3, ], polyline = TRUE)
#> $area
#> [1] 4.638237e-310
#> 
#> $perimeter
#> [1] 13332388
#> 
#> $n
#> [1] 3
#> 

# Area of Australia (approximate boundary)
australia <- cbind(
  lon = c(113, 153, 153, 142, 129, 113),
  lat = c(-26, -26, -10, -10, -15, -26)
)
result <- polygon_area(australia)
# Area in square kilometers
abs(result$area) / 1e6
#> [1] 5459214