Geodesy Helpers

Ground-measurement helpers in pyramids.base.geodesy. Where pyramids.base.crs answers what a coordinate means, these answer how far apart two coordinates are on the figure of the earth, and how much of a CRS's own unit a given ground distance occupies at a given place.

The ellipsoid always comes from the CRS's own datum, never from a hard-coded WGS 84 default — the same discipline Dataset.cell_area uses when it integrates cell areas.

Why the place matters

There is no single metres-per-degree constant. A degree of longitude is about 111 km at the equator and 0 at the pole, and Web Mercator stretches its own metre by 1 / cos(lat). So ground_distance_in_crs takes the location as a required argument rather than assuming one:

from pyramids.base.geodesy import ground_distance_in_crs

ground_distance_in_crs(100_000.0, crs=4326, at=(0.0, 0.0))    # 0.8983 degrees
ground_distance_in_crs(100_000.0, crs=4326, at=(0.0, 60.0))   # 1.7917 degrees

This is what sizing a scale bar needs: cleopatra's add_scale_bar takes a length in the axes' data units, and converting "100 km on the ground" into those units is a question about the Earth, not about pixels.

Coordinates are always geographic degrees

Every function here takes geographic degrees and uses crs only to select the ellipsoid. Projected coordinates are refused, not answered, because PROJ would return nan for an out-of-range latitude and that would propagate silently. Reproject first with reproject_coordinates.

at is always in (x, y) order — easting then northing, longitude then latitude — whatever axis order the CRS itself declares.

Vector counterparts

FeatureCollection wraps these as the geodesic counterparts of the planar operations it inherits from GeoDataFrame, for callers who hold geometries rather than coordinates:

Geodesic	Planar twin
FeatureCollection.geodesic_distance	.distance
FeatureCollection.geodesic_length	.length
FeatureCollection.geodesic_area	.area

See the FeatureCollection reference.

Functions

pyramids.base.geodesy.geodesic_distance(lon1, lat1, lon2, lat2, *, crs=4326, unit='m')

Ground distance between geographic points, on the CRS's own ellipsoid.

Solves the inverse geodesic problem -- the shortest path across the ellipsoid, not a great circle on an assumed sphere, which differs by up to ~0.5 % and by ~20 km on a pole-to-pole span.

All four coordinate arguments are geographic degrees, whatever crs is: crs selects the ellipsoid to measure on (through its datum), not the frame the inputs are expressed in. Pass projected coordinates through :func:pyramids.base.crs.reproject_coordinates first -- handing them in directly is refused, not answered, because pyproj.Geod would return nan for an out-of-range latitude and that would propagate silently.

Scalars in, scalar out; arrays in, array out -- the loop runs inside PROJ, so a transect of thousands of vertices costs one call. The arrays must all be the same length: PROJ pairs them elementwise and does not broadcast, so a scalar cannot be mixed with an array.

Parameters:

Name	Type	Description	Default
lon1	Any	Longitude(s) of the first point, in degrees.	required
lat1	Any	Latitude(s) of the first point, in degrees.	required
lon2	Any	Longitude(s) of the second point, in degrees.	required
lat2	Any	Latitude(s) of the second point, in degrees.	required
crs	Any	The CRS whose datum names the ellipsoid to measure on. Any form :func:pyramids.base.crs.crs_from_user_input accepts. Default 4326 (WGS 84). A projected CRS is accepted and contributes its own datum's ellipsoid.	4326
unit	str	m (default), km, mi or nmi.	'm'

Returns:

Type	Description
float | FloatArray	float | FloatArray: The distance(s), in unit. A float when every
float | FloatArray	coordinate argument is scalar, otherwise an array of the same shape as
float | FloatArray	the inputs, which PROJ pairs elementwise.

Raises:

Type	Description
ValueError	unit is not recognised, a coordinate is not numeric or not finite, or a latitude lies outside +-90. The message names the offending argument. Longitude is not range-checked: CF uses both the -180..180 and 0..360 conventions and Geod wraps it itself.
CRSError	crs cannot be resolved, or its datum names no ellipsoid.

Examples:

• A degree of longitude at the equator:

  >>> from pyramids.base.geodesy import geodesic_distance
  >>> round(geodesic_distance(0.0, 0.0, 1.0, 0.0))
  111319
  

• The same degree at 60 degrees north is about half as long, which is the whole reason a ground distance cannot be read off a cell size:

  >>> from pyramids.base.geodesy import geodesic_distance
  >>> round(geodesic_distance(0.0, 60.0, 1.0, 60.0))
  55799
  

• Vectorised, and in kilometres:

  >>> from pyramids.base.geodesy import geodesic_distance
  >>> lon1, lat1 = [0.0, 0.0], [0.0, 60.0]
  >>> lon2, lat2 = [1.0, 1.0], [0.0, 60.0]
  >>> d = geodesic_distance(lon1, lat1, lon2, lat2, unit="km")
  >>> [round(float(v), 1) for v in d]
  [111.3, 55.8]
  

• An unrecognised unit is refused, and the message lists the valid ones:

  >>> from pyramids.base.geodesy import geodesic_distance
  >>> geodesic_distance(0.0, 0.0, 1.0, 0.0, unit="furlong")
  Traceback (most recent call last):
      ...
  ValueError: unknown length unit 'furlong'; expected one of km, m, mi, nmi
  

See Also

ground_distance_in_crs: The inverse -- a ground distance expressed in a CRS's own units.

Source code in src/pyramids/base/geodesy.py

def geodesic_distance(
    lon1: Any,
    lat1: Any,
    lon2: Any,
    lat2: Any,
    *,
    crs: Any = 4326,
    unit: str = "m",
) -> float | FloatArray:
    """Ground distance between geographic points, on the CRS's own ellipsoid.

    Solves the inverse geodesic problem -- the shortest path across the
    ellipsoid, not a great circle on an assumed sphere, which differs by up to
    ~0.5 % and by ~20 km on a pole-to-pole span.

    All four coordinate arguments are **geographic degrees**, whatever `crs` is:
    `crs` selects the *ellipsoid* to measure on (through its datum), not the
    frame the inputs are expressed in. Pass projected coordinates through
    :func:`pyramids.base.crs.reproject_coordinates` first -- handing them in
    directly is refused, not answered, because `pyproj.Geod` would return `nan`
    for an out-of-range latitude and that would propagate silently.

    Scalars in, scalar out; arrays in, array out -- the loop runs inside PROJ,
    so a transect of thousands of vertices costs one call. The arrays must all be
    the **same length**: PROJ pairs them elementwise and does not broadcast, so a
    scalar cannot be mixed with an array.

    Args:
        lon1: Longitude(s) of the first point, in degrees.
        lat1: Latitude(s) of the first point, in degrees.
        lon2: Longitude(s) of the second point, in degrees.
        lat2: Latitude(s) of the second point, in degrees.
        crs: The CRS whose datum names the ellipsoid to measure on. Any form
            :func:`pyramids.base.crs.crs_from_user_input` accepts. Default
            `4326` (WGS 84). A projected CRS is accepted and contributes its
            own datum's ellipsoid.
        unit: `m` (default), `km`, `mi` or `nmi`.

    Returns:
        float | FloatArray: The distance(s), in `unit`. A float when every
        coordinate argument is scalar, otherwise an array of the same shape as
        the inputs, which PROJ pairs elementwise.

    Raises:
        ValueError: `unit` is not recognised, a coordinate is not numeric or not
            finite, or a latitude lies outside +-90. The message names the
            offending argument. Longitude is not range-checked: CF uses both the
            -180..180 and 0..360 conventions and `Geod` wraps it itself.
        CRSError: `crs` cannot be resolved, or its datum names no ellipsoid.

    Examples:
        - A degree of longitude at the equator:
            ```python
            >>> from pyramids.base.geodesy import geodesic_distance
            >>> round(geodesic_distance(0.0, 0.0, 1.0, 0.0))
            111319

            ```
        - The same degree at 60 degrees north is about half as long, which is
          the whole reason a ground distance cannot be read off a cell size:
            ```python
            >>> from pyramids.base.geodesy import geodesic_distance
            >>> round(geodesic_distance(0.0, 60.0, 1.0, 60.0))
            55799

            ```
        - Vectorised, and in kilometres:
            ```python
            >>> from pyramids.base.geodesy import geodesic_distance
            >>> lon1, lat1 = [0.0, 0.0], [0.0, 60.0]
            >>> lon2, lat2 = [1.0, 1.0], [0.0, 60.0]
            >>> d = geodesic_distance(lon1, lat1, lon2, lat2, unit="km")
            >>> [round(float(v), 1) for v in d]
            [111.3, 55.8]

            ```
        - An unrecognised unit is refused, and the message lists the valid ones:
            ```python
            >>> from pyramids.base.geodesy import geodesic_distance
            >>> geodesic_distance(0.0, 0.0, 1.0, 0.0, unit="furlong")
            Traceback (most recent call last):
                ...
            ValueError: unknown length unit 'furlong'; expected one of km, m, mi, nmi

            ```

    See Also:
        ground_distance_in_crs: The inverse -- a ground distance expressed in a
            CRS's own units.
    """
    scale = _length_scale(unit)
    _, geod = _resolve_geod(crs)
    x1 = _as_degrees("lon1", lon1)
    y1 = _as_degrees("lat1", lat1, latitude=True)
    x2 = _as_degrees("lon2", lon2)
    y2 = _as_degrees("lat2", lat2, latitude=True)
    # `inv` returns (forward azimuth, back azimuth, distance); only the third is
    # wanted here, and it accepts scalars and arrays alike.
    try:
        _, _, metres = geod.inv(x1, y1, x2, y2)
    except (ProjError, GeodError) as exc:
        # `GeodError` is a sibling of `ProjError`, not a subclass, so both have
        # to be named. Mismatched array lengths arrive here.
        raise ValueError(
            f"could not measure a geodesic distance: {exc} All four coordinate "
            "arguments must be scalars or arrays of the same length."
        ) from exc
    scaled = np.asarray(metres, dtype=float) / scale
    result: float | FloatArray = float(scaled) if scaled.ndim == 0 else scaled
    return result

pyramids.base.geodesy.ground_distance_in_crs(distance_m, *, crs, at, azimuth=_EAST)

How many of crs's own units span distance_m of ground at at.

The question a scale bar asks: cleopatra's add_scale_bar takes a length in the axes' data units, so sizing a "100 km" bar means converting 100 km of ground into the display CRS's units -- and that conversion depends on where on the map it is measured. A degree of longitude runs from 111 km at the equator to 0 at the pole, and Web Mercator stretches by 1 / cos(lat), so at is required rather than defaulted.

The ground distance is walked out from at along azimuth by the direct geodesic problem, and the straight-line separation of the two endpoints is measured back in crs. Because a geodesic heading east is not a parallel, the endpoint drifts slightly in latitude; the returned value is the full separation, which is what a bar drawn between those two points spans.

Parameters:

Name	Type	Description	Default
distance_m	float	The ground distance, in metres. Must be finite and positive.	required
crs	Any	The CRS to express the answer in, in any form :func:pyramids.base.crs.crs_from_user_input accepts. The result is in this CRS's units -- degrees for a geographic CRS, and the CRS's own linear unit (not necessarily metres) for a projected one.	required
at	tuple[float, float]	(x, y) in crs, the place the distance is measured at. Always in x, y order -- easting then northing, or longitude then latitude -- whatever axis order crs itself declares, matching always_xy=True and :func:pyramids.base.crs.reproject_coordinates. Swapping the two is silent for a projected CRS, so it is worth getting right; for a geographic one it is caught when the second member exceeds a latitude.	required
azimuth	float	Direction to measure along, in degrees clockwise from north. Default 90.0 (due east), the axis a horizontal scale bar spans. Wraps, so 450.0 is also due east; it must be finite.	_EAST

Returns:

Name	Type	Description
float	float	The span in crs's own units.

Raises:

Type	Description
ValueError	distance_m is not a finite positive real number, azimuth is not a finite real number, at is not a two-member sequence of finite numbers, or at (or the point distance_m away from it along azimuth) falls outside the CRS's usable domain. A bool is refused for either number, and a mapping or string is refused for at, since both would otherwise pass a length check and measure elsewhere.
CRSError	crs cannot be resolved, its datum names no ellipsoid, or it has no geographic counterpart to walk the geodesic in.

Examples:

• 100 km of ground at the equator, in degrees:

  >>> from pyramids.base.geodesy import ground_distance_in_crs
  >>> round(ground_distance_in_crs(100_000.0, crs=4326, at=(0.0, 0.0)), 4)
  0.8983
  

• The same 100 km at 60 degrees north needs about twice as many degrees:

  >>> from pyramids.base.geodesy import ground_distance_in_crs
  >>> round(ground_distance_in_crs(100_000.0, crs=4326, at=(0.0, 60.0)), 4)
  1.7917
  

• Web Mercator is in metres, but not ground metres away from the equator -- at 60 north its unit is stretched by about two:

  >>> from pyramids.base.geodesy import ground_distance_in_crs
  >>> y60 = 8399737.89
  >>> round(ground_distance_in_crs(100_000.0, crs=3857, at=(0.0, y60)) / 1000)
  199
  

• A non-positive ground distance is refused rather than returning zero:

  >>> from pyramids.base.geodesy import ground_distance_in_crs
  >>> ground_distance_in_crs(-1.0, crs=4326, at=(0.0, 0.0))
  Traceback (most recent call last):
      ...
  ValueError: distance_m must be positive, got -1.0
  

See Also

geodesic_distance: The forward direction -- ground metres between two geographic points.

Source code in src/pyramids/base/geodesy.py

def ground_distance_in_crs(
    distance_m: float,
    *,
    crs: Any,
    at: tuple[float, float],
    azimuth: float = _EAST,
) -> float:
    """How many of `crs`'s own units span `distance_m` of ground at `at`.

    The question a scale bar asks: cleopatra's `add_scale_bar` takes a length in
    the axes' data units, so sizing a "100 km" bar means converting 100 km of
    ground into the display CRS's units -- and that conversion depends on where
    on the map it is measured. A degree of longitude runs from 111 km at the
    equator to 0 at the pole, and Web Mercator stretches by `1 / cos(lat)`, so
    `at` is required rather than defaulted.

    The ground distance is walked out from `at` along `azimuth` by the direct
    geodesic problem, and the straight-line separation of the two endpoints is
    measured back in `crs`. Because a geodesic heading east is not a parallel,
    the endpoint drifts slightly in latitude; the returned value is the full
    separation, which is what a bar drawn between those two points spans.

    Args:
        distance_m: The ground distance, in **metres**. Must be finite and
            positive.
        crs: The CRS to express the answer in, in any form
            :func:`pyramids.base.crs.crs_from_user_input` accepts. The result is
            in *this* CRS's units -- degrees for a geographic CRS, and the CRS's
            own linear unit (not necessarily metres) for a projected one.
        at: `(x, y)` in `crs`, the place the distance is measured at. Always in
            **x, y order** -- easting then northing, or longitude then latitude
            -- whatever axis order `crs` itself declares, matching
            `always_xy=True` and :func:`pyramids.base.crs.reproject_coordinates`.
            Swapping the two is silent for a projected CRS, so it is worth
            getting right; for a geographic one it is caught when the second
            member exceeds a latitude.
        azimuth: Direction to measure along, in degrees clockwise from north.
            Default `90.0` (due east), the axis a horizontal scale bar spans.
            Wraps, so `450.0` is also due east; it must be finite.

    Returns:
        float: The span in `crs`'s own units.

    Raises:
        ValueError: `distance_m` is not a finite positive real number, `azimuth`
            is not a finite real number, `at` is not a two-member sequence of
            finite numbers, or `at` (or the point `distance_m` away from it along
            `azimuth`) falls outside the CRS's usable domain. A bool is refused
            for either number, and a mapping or string is refused for `at`, since
            both would otherwise pass a length check and measure elsewhere.
        CRSError: `crs` cannot be resolved, its datum names no ellipsoid, or it
            has no geographic counterpart to walk the geodesic in.

    Examples:
        - 100 km of ground at the equator, in degrees:
            ```python
            >>> from pyramids.base.geodesy import ground_distance_in_crs
            >>> round(ground_distance_in_crs(100_000.0, crs=4326, at=(0.0, 0.0)), 4)
            0.8983

            ```
        - The same 100 km at 60 degrees north needs about twice as many degrees:
            ```python
            >>> from pyramids.base.geodesy import ground_distance_in_crs
            >>> round(ground_distance_in_crs(100_000.0, crs=4326, at=(0.0, 60.0)), 4)
            1.7917

            ```
        - Web Mercator is in metres, but not *ground* metres away from the
          equator -- at 60 north its unit is stretched by about two:
            ```python
            >>> from pyramids.base.geodesy import ground_distance_in_crs
            >>> y60 = 8399737.89
            >>> round(ground_distance_in_crs(100_000.0, crs=3857, at=(0.0, y60)) / 1000)
            199

            ```
        - A non-positive ground distance is refused rather than returning zero:
            ```python
            >>> from pyramids.base.geodesy import ground_distance_in_crs
            >>> ground_distance_in_crs(-1.0, crs=4326, at=(0.0, 0.0))
            Traceback (most recent call last):
                ...
            ValueError: distance_m must be positive, got -1.0

            ```

    See Also:
        geodesic_distance: The forward direction -- ground metres between two
            geographic points.
    """
    metres = _as_finite_scalar("distance_m", distance_m)
    if metres <= 0.0:
        raise ValueError(f"distance_m must be positive, got {distance_m!r}")
    bearing = _as_finite_scalar("azimuth", azimuth)
    x, y = _as_point("at", at)
    target, geod = _resolve_geod(crs)
    geodetic, to_degrees = _geodetic_frame(target)
    if target.is_geographic and abs(y * to_degrees) > 90.0:
        # Only detectable on a geographic CRS, and only for the half of the
        # mistake that lands out of range -- but that half includes the common
        # one, a projected northing handed to a geographic CRS.
        raise ValueError(
            f"at must be (x, y) -- longitude then latitude -- for the geographic "
            f"CRS {target.name!r}, but its second member is {y}, which is not a "
            "latitude. Swap the pair, or pass coordinates in the CRS you named."
        )
    try:
        to_lonlat = Transformer.from_crs(target, geodetic, always_xy=True)
        native_lon, native_lat = to_lonlat.transform(x, y)
        # Into degrees before `Geod` sees them, and back out afterwards: the
        # transform answers in the geodetic CRS's own angular unit, which is
        # grads for the NTF (Paris) family. See `_geodetic_frame`.
        lon, lat = native_lon * to_degrees, native_lat * to_degrees
        # The geodesic is walked in lon/lat and the endpoint brought back, so the
        # answer is measured in `crs` rather than assumed proportional to it.
        lon_end, lat_end, _ = geod.fwd(lon, lat, bearing, metres)
        to_target = Transformer.from_crs(geodetic, target, always_xy=True)
        x_end, y_end = to_target.transform(lon_end / to_degrees, lat_end / to_degrees)
    except (ProjError, GeodError) as exc:
        raise ValueError(
            f"could not measure {metres} m at {(x, y)!r} in {target.name!r}: {exc}"
        ) from exc
    span = float(np.hypot(x_end - x, y_end - y))
    if not np.isfinite(span):
        # PROJ reports an un-invertible point as `inf` rather than raising, so a
        # location outside the CRS's domain arrives here, not in the except.
        raise ValueError(
            f"the point {(x, y)!r} (or the point {metres} m from it along "
            f"azimuth {bearing}) falls outside the usable domain of "
            f"{target.name!r}, so the span is undefined"
        )
    return span

pyramids.base.geodesy.geodesic_geometry_length(geometry, *, crs=4326, unit='m')

Length of a geometry's lines along the ellipsoid.

The geodesic counterpart of shapely's planar .length, which on a geographic CRS measures in degrees and is meaningless as a ground distance.

geometry must already be in geographic coordinates; crs only selects the ellipsoid. A polygon reports its perimeter, a point reports 0.0, and a multi-part geometry sums its parts -- all of which is pyproj.Geod's own behaviour, passed through rather than reinterpreted.

Parameters:

Name	Type	Description	Default
geometry	Any	Any shapely geometry, in geographic coordinates.	required
crs	Any	The CRS whose datum names the ellipsoid. Default 4326.	4326
unit	str	m (default), km, mi or nmi.	'm'

Returns:

Name	Type	Description
float	float	The length in unit.

Raises:

Type	Description
ValueError	unit is not recognised, geometry is not a shapely geometry, or its coordinates are not finite geographic degrees -- a projected geometry is refused rather than answered with nan.
CRSError	crs cannot be resolved, or its datum names no ellipsoid.

Examples:

• A one-degree line at the equator is ~111 km long:

  >>> from shapely.geometry import LineString
  >>> from pyramids.base.geodesy import geodesic_geometry_length
  >>> round(geodesic_geometry_length(LineString([(0, 0), (1, 0)])))
  111319
  

• A polygon reports its perimeter, in kilometres here:

  >>> from shapely.geometry import Polygon
  >>> from pyramids.base.geodesy import geodesic_geometry_length
  >>> square = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
  >>> round(geodesic_geometry_length(square, unit="km"), 1)
  443.8
  

• A point has no length:

  >>> from shapely.geometry import Point
  >>> from pyramids.base.geodesy import geodesic_geometry_length
  >>> geodesic_geometry_length(Point(12.5, 41.9))
  0.0
  

See Also

geodesic_geometry_area: The area counterpart.

Source code in src/pyramids/base/geodesy.py

def geodesic_geometry_length(
    geometry: Any,
    *,
    crs: Any = 4326,
    unit: str = "m",
) -> float:
    """Length of a geometry's lines along the ellipsoid.

    The geodesic counterpart of shapely's planar `.length`, which on a
    geographic CRS measures in degrees and is meaningless as a ground distance.

    `geometry` must already be in **geographic** coordinates; `crs` only selects
    the ellipsoid. A polygon reports its perimeter, a point reports `0.0`, and a
    multi-part geometry sums its parts -- all of which is `pyproj.Geod`'s own
    behaviour, passed through rather than reinterpreted.

    Args:
        geometry: Any shapely geometry, in geographic coordinates.
        crs: The CRS whose datum names the ellipsoid. Default `4326`.
        unit: `m` (default), `km`, `mi` or `nmi`.

    Returns:
        float: The length in `unit`.

    Raises:
        ValueError: `unit` is not recognised, `geometry` is not a shapely
            geometry, or its coordinates are not finite geographic degrees -- a
            projected geometry is refused rather than answered with `nan`.
        CRSError: `crs` cannot be resolved, or its datum names no ellipsoid.

    Examples:
        - A one-degree line at the equator is ~111 km long:
            ```python
            >>> from shapely.geometry import LineString
            >>> from pyramids.base.geodesy import geodesic_geometry_length
            >>> round(geodesic_geometry_length(LineString([(0, 0), (1, 0)])))
            111319

            ```
        - A polygon reports its perimeter, in kilometres here:
            ```python
            >>> from shapely.geometry import Polygon
            >>> from pyramids.base.geodesy import geodesic_geometry_length
            >>> square = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
            >>> round(geodesic_geometry_length(square, unit="km"), 1)
            443.8

            ```
        - A point has no length:
            ```python
            >>> from shapely.geometry import Point
            >>> from pyramids.base.geodesy import geodesic_geometry_length
            >>> geodesic_geometry_length(Point(12.5, 41.9))
            0.0

            ```

    See Also:
        geodesic_geometry_area: The area counterpart.
    """
    scale = _length_scale(unit)
    _, geod = _resolve_geod(crs)
    _geometry_in_degrees("geometry", geometry)
    try:
        metres = geod.geometry_length(geometry)
    except (ProjError, GeodError) as exc:
        raise ValueError(f"could not measure {geometry!r}: {exc}") from exc
    # The bounds check cannot see a nan vertex in the interior of a line --
    # shapely's `bounds` skips it -- so the result is checked as well. `Geod`
    # answers out-of-domain input with nan rather than raising.
    if not math.isfinite(metres):
        raise ValueError(
            f"measuring {geometry!r} produced a non-finite length; its "
            "coordinates are not usable geographic degrees"
        )
    return float(metres) / scale

pyramids.base.geodesy.geodesic_geometry_area(geometry, *, crs=4326, unit='m2')

Area of a geometry on the ellipsoid.

The geodesic counterpart of shapely's planar .area, which on a geographic CRS reports square degrees -- a quantity that varies with latitude and so cannot be compared between rows of a grid.

geometry must already be in geographic coordinates; crs only selects the ellipsoid. The result is always non-negative: pyproj.Geod signs the area by ring orientation (negative for a clockwise ring), which describes winding rather than size, so the magnitude is returned. A line or point reports 0.0.

Parameters:

Name	Type	Description	Default
geometry	Any	Any shapely geometry, in geographic coordinates.	required
crs	Any	The CRS whose datum names the ellipsoid. Default 4326.	4326
unit	str	m2 (default), km2 or ha.	'm2'

Returns:

Name	Type	Description
float	float	The area in unit, never negative.

Raises:

Type	Description
ValueError	unit is not recognised, geometry is not a shapely geometry, or its coordinates are not finite geographic degrees -- a projected geometry is refused rather than answered with nan.
CRSError	crs cannot be resolved, or its datum names no ellipsoid.

Examples:

• A one-degree square at the equator covers ~12 309 square kilometres:

  >>> from shapely.geometry import Polygon
  >>> from pyramids.base.geodesy import geodesic_geometry_area
  >>> square = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
  >>> round(geodesic_geometry_area(square, unit="km2"), 1)
  12308.8
  

• Ring orientation does not change the size:

  >>> from shapely.geometry import Polygon
  >>> from pyramids.base.geodesy import geodesic_geometry_area
  >>> clockwise = Polygon([(0, 0), (0, 1), (1, 1), (1, 0)])
  >>> round(geodesic_geometry_area(clockwise, unit="km2"), 1)
  12308.8
  

• A line encloses nothing:

  >>> from shapely.geometry import LineString
  >>> from pyramids.base.geodesy import geodesic_geometry_area
  >>> geodesic_geometry_area(LineString([(0, 0), (1, 0)]))
  0.0
  

See Also

geodesic_geometry_length: The length counterpart.

Source code in src/pyramids/base/geodesy.py

def geodesic_geometry_area(
    geometry: Any,
    *,
    crs: Any = 4326,
    unit: str = "m2",
) -> float:
    """Area of a geometry on the ellipsoid.

    The geodesic counterpart of shapely's planar `.area`, which on a geographic
    CRS reports square degrees -- a quantity that varies with latitude and so
    cannot be compared between rows of a grid.

    `geometry` must already be in **geographic** coordinates; `crs` only selects
    the ellipsoid. The result is always non-negative: `pyproj.Geod` signs the
    area by ring orientation (negative for a clockwise ring), which describes
    winding rather than size, so the magnitude is returned. A line or point
    reports `0.0`.

    Args:
        geometry: Any shapely geometry, in geographic coordinates.
        crs: The CRS whose datum names the ellipsoid. Default `4326`.
        unit: `m2` (default), `km2` or `ha`.

    Returns:
        float: The area in `unit`, never negative.

    Raises:
        ValueError: `unit` is not recognised, `geometry` is not a shapely
            geometry, or its coordinates are not finite geographic degrees -- a
            projected geometry is refused rather than answered with `nan`.
        CRSError: `crs` cannot be resolved, or its datum names no ellipsoid.

    Examples:
        - A one-degree square at the equator covers ~12 309 square kilometres:
            ```python
            >>> from shapely.geometry import Polygon
            >>> from pyramids.base.geodesy import geodesic_geometry_area
            >>> square = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
            >>> round(geodesic_geometry_area(square, unit="km2"), 1)
            12308.8

            ```
        - Ring orientation does not change the size:
            ```python
            >>> from shapely.geometry import Polygon
            >>> from pyramids.base.geodesy import geodesic_geometry_area
            >>> clockwise = Polygon([(0, 0), (0, 1), (1, 1), (1, 0)])
            >>> round(geodesic_geometry_area(clockwise, unit="km2"), 1)
            12308.8

            ```
        - A line encloses nothing:
            ```python
            >>> from shapely.geometry import LineString
            >>> from pyramids.base.geodesy import geodesic_geometry_area
            >>> geodesic_geometry_area(LineString([(0, 0), (1, 0)]))
            0.0

            ```

    See Also:
        geodesic_geometry_length: The length counterpart.
    """
    scale = _area_scale(unit)
    _, geod = _resolve_geod(crs)
    _geometry_in_degrees("geometry", geometry)
    try:
        area, _ = geod.geometry_area_perimeter(geometry)
    except (ProjError, GeodError) as exc:
        raise ValueError(f"could not measure {geometry!r}: {exc}") from exc
    # As in `geodesic_geometry_length`: a nan vertex inside a ring does not reach
    # `bounds`, so the result is checked too.
    if not math.isfinite(area):
        raise ValueError(
            f"measuring {geometry!r} produced a non-finite area; its "
            "coordinates are not usable geographic degrees"
        )
    # `Geod` signs the area by ring orientation; a clockwise ring is negative.
    # That encodes winding, not size, and every caller here wants the size.
    return abs(float(area)) / scale