# Operations

The `geo.` operations measure distances, lengths, directions and areas,
and say whether one shape touches or holds another:

| Operation | Gives | Fails | SQL form |
|---|---|---|---|
| `geo.distance(from, to)` | `Length`: the geodesic between two `Point`s on the WGS84 ellipsoid, rounded half to even to the micrometre (`Length`'s scale) | | `ST_Distance(a, b)` |
| `geo.length(line)` | `Length`: a `LineString`'s or `MultiLineString`'s geodesics summed, each edge measured as by `distance` | | `ST_Length(x)` |
| `geo.bearing(from, to)` | `Bearing`: the geodesic's direction as it leaves `from`, degrees clockwise from true north, at least 0 and below 360, to 9 places | `same_point` | `degrees(ST_Azimuth(a, b))` |
| `geo.area(shape)` | `Area`: a `Polygon`'s or `MultiPolygon`'s area on the ellipsoid, holes taken away, to `Area`'s scale | | `ST_Area(x)` |
| `geo.intersects(a, b)` | `Bool`: the two share at least one point, edges included | | `ST_Intersects(a, b)` |
| `geo.contains(a, b)` | `Bool`: no point of `b` lies outside `a`; `a` may also be a `BoundingBox` or a `Circle` | | `ST_Covers(a, b)` |
| `geo.within(a, b)` | `Bool`: `geo.contains(b, a)` | | `ST_CoveredBy(a, b)` |
| `geo.within(a, b, distance:)` | `Bool`: `geo.distance(a, b) <= distance`, both `Point`s, `distance` a `Length` never below 0 | | `ST_DWithin(a, b, d)` |
| `geo.centroid(shape)` | `Point`: the centre on the sphere, weighted by area, length or count by the shape's kind; may lie outside a concave polygon | weights that cancel (two antipodal points) are a defect | `ST_Centroid(x, false)` |
| `geo.bounding_box(shape)` | `BoundingBox`: the smallest box of parallels and meridians holding the shape, an edge's highest or lowest point counted; across the antimeridian when that box is narrower | | none |

Each SQL form is PostGIS's own function for the operation, which
[aggregate queries](counts-sums-and-groups.md#the-aggregates) can use. A
`where:` condition doesn't compile to its form: it searches a little
wider and leaves the answer to Jadpo's code, as [Queries](geospatial-queries.md)
explains. `bounding_box` and `centroid` have no use in a query.

Here a bakery measures its van's trips, its round and its patches:

```jadpo
type SameSpot = Failure {
    kind: FailureKind.invalid_value
}

type VanLeg = Object {
    trip_length: Length
    bearing_out: Bearing
}

type RoundLength = Length {}

type PatchArea = Area {}

/// How far a bakery's van drives to a customer, and which way it sets off.
function first_leg(bakery: Point, customer: Point) fails SameSpot -> VanLeg {
    var setting_off: Bearing = attempt geo.bearing(bakery, customer) {
        same_point: SameSpot  // a pickup at the counter
    }
    return VanLeg {
        trip_length: geo.distance(bakery, customer)
        bearing_out: setting_off
    }
}

/// How far a van drives on its round, stop to stop.
function round_length(van_round: LineString) -> RoundLength {
    return RoundLength(geo.length(van_round))
}

/// How much ground a delivery patch covers, its holes left out.
function patch_area(patch: Polygon) -> PatchArea {
    return PatchArea(geo.area(patch))
}
```

A `Point` literal sets every field, its `altitude` as `none` when it has
none ([Objects](objects.md)). The tests check the
figures above:

```jadpo
test "London to Paris is measured on the ellipsoid" {
    var big_ben: Point = Point {
        latitude: 51.5072
        longitude: -0.1276
        altitude: none
    }
    var eiffel_tower: Point = Point {
        latitude: 48.8566
        longitude: 2.3522
        altitude: none
    }
    assert geo.distance(big_ben, eiffel_tower) == 343896.891267metre
    assert geo.distance(eiffel_tower, big_ben) == geo.distance(big_ben, eiffel_tower)
    var leg: VanLeg = attempt first_leg(big_ben, eiffel_tower)
    assert leg.bearing_out == 148.046085477
    assert first_leg(big_ben, big_ben) fails SameSpot
}

test "a round is its legs summed" {
    var van_round: LineString = LineString {
        points: [
            Point {
                latitude: 51.5072
                longitude: -0.1276
                altitude: none
            },
            Point {
                latitude: 48.8566
                longitude: 2.3522
                altitude: none
            },
            Point {
                latitude: 51.5072
                longitude: -0.1276
                altitude: none
            }
        ]
    }
    assert round_length(van_round) == 687793.782534metre
}
```

## Distance and length

`geo.distance` agrees with GeographicLib to well under a micrometre
before it rounds: it is within 75 femtometres of every distance in
GeographicLib's published test set
([Precision and determinism](precision-and-determinism.md)). PostGIS's
doubles are good to 15 nanometres. London (51.5072° N, 0.1276° W) to
Paris (48.8566° N, 2.3522° E) is 343,896.891267 m. On a sphere it would
be 343,530.338 m, 367 m short (both verified with PostGIS).

`geo.length` adds up each edge's distance as `geo.distance` rounds it, so
a line of two points is exactly their distance. A `MultiLineString`'s
lines are added up the same way. No distance on Earth needs more than 38
digits, and a sum that does is a defect, as it is for `+`.

## Bearing

`geo.bearing` gives a direction to 9 decimal places, as an `Angle` keeps
radians to 9. Two equal points have no direction, so a call maps the
failure `same_point`, as a local time that doesn't exist is mapped
([Zones and local times](zones-and-local-times.md)).

Two points are the same point when they have the same latitude and the
same longitude, -180 being 180, or are the same pole at any longitude.
Two points a picometre apart do have a direction. Whether two points are
equal is known only when the program runs, so `same_point` is settled
only by `attempt`, as decoding is. `same_point` exists only in code:
`ST_Azimuth` gives `NULL` for two equal points (verified), and a query
never needs a bearing.

A bearing that rounds to 360 is 0. At a pole the bearing is
GeographicLib's, measured from the meridian of the longitude the pole is
given with.

## Area

`geo.area` uses Karney's method for polygons, each edge a geodesic.
PostGIS's `ST_Area` gives 12,308,778,361.47 m² for the 1° square at the
equator (verified); Jadpo gives 12,308,778,361.469453 m², and each
geodesic's area matches GeographicLib's test set to the square
millimetre.

It rounds once. Each ring's area, the smaller of the two sides the ring
bounds, is worked out without rounding. The holes are taken away and a
`MultiPolygon`'s polygons added, and then the total is rounded half to
even to the square millimetre. A hole always lies inside its outer ring
(the `rings_valid` rule under
[Polygons are valid when built](predicates-and-shapes.md#polygons-are-valid-when-built)), so holes
never outweigh it.

## Centroid

`geo.centroid` gives a point to put a shape's label at. It isn't a
measure. It is worked out on a sphere, as PostGIS's
`ST_Centroid(x, false)` is: each point becomes a unit vector, the vectors
are weighted and added up, and the sum is turned back into latitude and
longitude.

It gives PostGIS's answer, quirks included:

- each point of a `MultiPoint` weighs one;
- each edge of a line gives both of its ends the edge's angle as their
  weight;
- a polygon is cut into triangles, from the first corner of its first
  outer ring to every edge of every ring. Each triangle's centre, its
  corners added up, is weighted by the triangle's area, always counted
  positive.

So a hole adds weight rather than taking it away. A concave polygon's
centroid moves with the corner its ring starts at, as PostGIS's does, and
may lie outside the polygon. Verified: an L-shaped polygon's centroid is
3.214° N, 3.236° E from one starting corner and 4.151° N, 4.159° E from
another.

A point's centroid is the point itself, without its altitude. When the
weights cancel, as for two opposite points, the call is a defect, as
PostGIS gives no number there either. The result is rounded to the
nearest double, as every coordinate is.

## `contains`, `within` and `intersects`

`geo.contains(a, b)` is true when no point of `b` lies outside `a`, and
`geo.within(a, b)` is `geo.contains(b, a)`. `geo.intersects(a, b)` is
true when the two share at least one point, edges included. A
`MultiPolygon` is the union of its parts, which may overlap.

The SQL form, `ST_Covers` on `geography`, compares two polygons ring by
ring, so it parts from that meaning in four cases:

- it doesn't see a hole of `a` that `b` covers: it says a square with a
  hole covers a square over the hole;
- it counts rings lying along each other as crossing: it says a square
  doesn't cover the same square with a hole;
- with an `a` of several parts, it wants all of `b` in one part, so two
  points each in a different part aren't covered.

`geo.contains` keeps to "no point of `b` lies outside `a`" in each of
these cases, in code and in a query ([Queries](geospatial-queries.md)).

## A box as the container

With a `BoundingBox` as `a`, `geo.contains` compares latitudes and
longitudes directly, because a box's edges are parallels and meridians,
not great circles. In code a box can hold any geometry; an edge's highest
or lowest point, and the longitudes between its ends, count. In a query
the box filter works on `Point` fields only ([Queries](geospatial-queries.md)).

## A circle as the container

With a `Circle`, `geo.contains` is `geo.within(distance:)` of the
circle's centre, the radius being the distance. `within(distance:)` measures between points only, so a
circle holds points only. `geo.contains(circle, shape)` and
`geo.within(shape, circle)` with anything but a `Point` are
`type.mismatch`, and so is `geo.within(a, b, distance:)` unless both are
points. The message cites the decision on point-only distance.

## `within` a distance

`geo.within(a, b, distance: d)` is true when the two points are at most
`d` apart. The distance must be proved never below 0 while compiling, as
a count is, or the call is `type.argument_negative`. A length written
out, or a value of a type with `min: 0metre`, is proved.

## Bounding box

`geo.bounding_box` gives the smallest box of parallels and meridians that
holds the shape. A side worked out from an edge is rounded outwards to a
double, so a shape is always inside its own box. The box crosses the
antimeridian when that box is narrower; when two boxes are equally
narrow, it takes the one that doesn't cross.

A shape that holds a pole reaches every longitude: a point at a pole, a
line through one, or a polygon around or touching one. Its box runs from
-180 to 180, and only such a box holds it, whatever longitude a point at
the pole is written with, since every meridian meets there.

## An edge between opposite points

There is no one great circle between two opposite points of the globe,
so a line's edge between them is just its two ends.

## Operations Jadpo doesn't have

`buffer` isn't an operation ([What geospatial leaves out, and what to use instead](geospatial-left-out.md)).
Neither is the distance between a point and a line or polygon. It comes
when one of the standard library's
[scenario applications](scenario-applications.md)
needs it.