Calculate area of geometry in meters with Python/Shapely

calculate_area.py

import pyproj

from shapely import wkt, geometry, ops
from functools import partial


def get_utm_zone(point: geometry.Point) -> int:
    """
    Helper function to get the UTM Zone EPSG of the input point.

    Parameters:
        point: shapely.geometry.Point

    Returns:
        int: EPSG of the UTM zone of the input point
    """
    prefix = 32600
    if point.y < 0:
        prefix = 32700

    zone = int(((point.x + 180) / 6) + 1)
    return prefix + zone


def calculate_area_using_utm(input_geometry) -> float:
    """
    Calculate the area of the input geometry in square meters.

    Transforming the geometry into UTM gives us a meter-based projection so the area we get back
    is in square meters.

    UTM is a conformal projection (https://en.wikipedia.org/wiki/Conformal_map_projection), so it
    accurately preserves angles at a small-to-medium scale. This makes it particularly effective
    for maps which need to describe directions. However, it is not necessarily ideal for calculating area.

    Parameters:
        input_geometry: shapely geometry to calculate area for

    Returns:
        float: area of the geometry in meters
    """
    utm_epsg = get_utm_zone(input_geometry.centroid)
    # Make sure you use the correct intitial reference here
    proj = partial(pyproj.transform,
                   pyproj.Proj(init='epsg:4326'),
                   pyproj.Proj(init='epsg:{}'.format(utm_epsg))
                   )
    transformed_geom = ops.transform(proj, input_geometry)
    return transformed_geom.area


def calculate_area_using_albers(input_geometry) -> float:
    """
    The Albers Equal Area projection (https://en.wikipedia.org/wiki/Albers_projection) is more effective
    for area calculation. It accurately represents area while distorting scale and shape, which we don't
    particularly care about for this application.

    Parameters:
        input_geometry: shapely geometry to calculate area for

    Returns:
        float: area of the geometry in meters
    """
    proj = partial(pyproj.transform,
                   pyproj.Proj(init='epsg:4326'),
                   pyproj.Proj(proj='aea', lat1=input_geometry.bounds[1], lat2=input_geometry.bounds[3])
                   )
    transformed_geom = ops.transform(proj, input_geometry)
    return transformed_geom.area


def calculate_area_using_web_mercator(input_geometry) -> float:
    """
    BEWARE this method. Web Mercator is a meter-based projection, but doesn't
    preserve area nearly as well as UTM, which is still not as effective as Albers.

    You will get exceedingly large (and incorrect) values with this method.

    Parameters:
        input_geometry: shapely geometry to calculate area for

    Returns:
        float: area of the geometry in meters
    """
    proj = partial(pyproj.transform,
                   pyproj.Proj(init='epsg:4326'),
                   pyproj.Proj(init='epsg:3857')
                   )
    transformed_geom = ops.transform(proj, input_geometry)
    return transformed_geom.area


if __name__ == '__main__':
    geom = wkt.loads('POLYGON (( -0.5 -0.5, 0.5 -0.5, 0.5 0.5, -0.5 0.5, -0.5 -0.5 ))')
    print("Area of a polygon at the equator")
    print('Area in degrees:{:>30}'.format(geom.area))
    print('Area using Web Mercator:{:>37}'.format(calculate_area_using_web_mercator(geom)))
    print('Area using UTM:{:>46}'.format(calculate_area_using_utm(geom)))
    print('Area using Albers Equal Area:{:>31}'.format(calculate_area_using_albers(geom)))
    print()
    geom = wkt.loads('POLYGON (( -70.5 44.5, -69.5 44.5, -69.5 45.5, -70.5 45.5, -70.5 44.5 ))')
    print("Area of a polygon at mid-latitudes")
    print('Area in degrees:{:>30}'.format(geom.area))
    print('Area using Web Mercator:{:>37}'.format(calculate_area_using_web_mercator(geom)))
    print('Area using UTM:{:>45}'.format(calculate_area_using_utm(geom)))
    print('Area using Albers Equal Area:{:>31}'.format(calculate_area_using_albers(geom)))