gazagnaire.org / ocaml-latlon
WGS 84 positions: great-circle distance, boxes and an index
Star 0
Fork 0
Watch
OCaml 88%
Python 4%
4%
Shell 2%
Dune 2%
15 1 0
fuzz
lib
test
.gitignore
.ocamlformat
dune
dune-project
latlon.opam
LICENSE.md
README.md
README.md

latlon #

WGS 84 positions and what you do with them: great-circle distance, indexes over points, segments and areas, and operations on bounding boxes.

Overview #

The library has four parts, kept together because callers use them together.

Latlon.t is a position in decimal degrees on WGS 84. The distance between two positions is the haversine distance in kilometres on the IUGG mean sphere.

Latlon.Index stores positions, each with a payload, and answers distance queries over them: nearest, k_nearest, within_km and inverse-distance interpolation (idw). Latlon.Segments finds the nearest of a set of line segments. Both give exact answers, including at the poles and across the antimeridian, and neither needs a cell size from the caller.

Latlon.Ring is a closed ring of positions, and Latlon.Rings indexes the box around each ring to find which area contains a point. The box lookup narrows the candidates and Ring.contains makes the decision. Its boundary rule means two neighbouring areas can never both claim a point on their shared edge.

Latlon.Box handles axis-aligned boxes. You can cut a box into tiles (for a server that will not return a whole country in one request), take the box around a disc, subtract one rectangle from another, and measure the distance from a box to a point. Subtraction gives you the part of a region that none of a set of covered areas reached.

It needs OCaml 4.14 or later, and nox-rtree for the index.

Install #

$ opam install latlon

Usage #

# let paris = Latlon.v ~lat:48.8566 ~lon:2.3522
val paris : Latlon.t = <abstr>
# let london = Latlon.v ~lat:51.5074 ~lon:(-0.1278)
val london : Latlon.t = <abstr>
# Latlon.distance_km paris london
- : float = 343.556534880883248

A position is abstract because its longitude is normalised, which gives each place a single representation. distance_km is periodic in longitude, so 225° and −135° are the same meridian to it. An index compares raw coordinates, and without normalisation a point stored as 225° would never match a query at −135°.

# Latlon.lon (Latlon.v ~lat:0. ~lon:225.)
- : float = -135.
# Latlon.equal (Latlon.v ~lat:0. ~lon:180.) (Latlon.v ~lat:0. ~lon:(-180.))
- : bool = true

An index carries a payload of your own type:

# let cities =
    [ paris, "Paris"; london, "London";
      Latlon.v ~lat:45.7640 ~lon:4.8357, "Lyon" ]
val cities : (Latlon.t * string) list =
  [(<abstr>, "Paris"); (<abstr>, "London"); (<abstr>, "Lyon")]
# let index = Latlon.Index.v cities
val index : string Latlon.Index.t = <abstr>
# Latlon.Index.nearest index (Latlon.v ~lat:47. ~lon:3.) |> Option.map snd
- : string option = Some "Lyon"
# Latlon.Index.k_nearest index (Latlon.v ~lat:47. ~lon:3.) 2 |> List.map snd
- : string list = ["Lyon"; "Paris"]

A box can be split into tiles, and one box subtracted from another:

# let window = Latlon.Box.v ~south:47. ~west:(-2.) ~north:51. ~east:2.
val window : Latlon.Box.t =
  {Latlon.Box.south = 47.; west = -2.; north = 51.; east = 2.}
# List.length (Latlon.Box.split ~degrees:2. window)
- : int = 4
# Latlon.Box.subtract window (Latlon.Box.v ~south:49. ~west:(-2.) ~north:51. ~east:2.)
- : Latlon.Box.t list =
[{Latlon.Box.south = 47.; west = -2.; north = 49.; east = 2.}]

Rings.containing lists the areas that hold a point. The index first keeps the rings whose box contains the point, then Ring.contains checks each of them exactly, so a point inside a ring's box but outside the ring belongs to no area:

# let ring corners =
    Latlon.Ring.v (List.map (fun (lat, lon) -> Latlon.v ~lat ~lon) corners)
val ring : (float * float) list -> Latlon.Ring.t = <fun>
# let areas =
    Latlon.Rings.v
      [ ring [ 48., 2.; 49., 2.; 49., 3.; 48., 3. ], "north";
        ring [ 47., 2.; 48., 2.; 48., 3.; 47., 3. ], "south" ]
val areas : string Latlon.Rings.t = <abstr>
# Latlon.Rings.containing areas (Latlon.v ~lat:48.5 ~lon:2.5)
- : string list = ["north"]
# Latlon.Rings.containing areas (Latlon.v ~lat:46. ~lon:2.5)
- : string list = []

The two areas share the parallel at 48°, and a point on it belongs to exactly one of them. A ring owns the southern and western parts of its boundary and leaves the northern and eastern parts to its neighbours, so adjacent areas split their shared edge between them.

# Latlon.Rings.containing areas (Latlon.v ~lat:48. ~lon:2.5)
- : string list = ["north"]

Box.distance_km uses the same haversine as Latlon.distance_km, so comparing it with a radius gives the same result as scanning the points one by one:

# Latlon.Box.distance_km window (Latlon.v ~lat:48. ~lon:(-1.))
- : float = 0.
# Latlon.Box.distance_km window (Latlon.v ~lat:48. ~lon:3.) > 50.
- : bool = true

Coordinates #

Positions are decimal degrees on WGS 84: latitude in [-90, 90], longitude in [-180, 180], in that order everywhere in this library. GeoJSON (RFC 7946 Section 3.1.1) writes them the other way round, longitude first, and its bbox member is [west, south, east, north], so a Latlon.Box.t needs reordering on the way in and out.

Distances are great-circle arcs on a sphere of the IUGG mean radius (6371.0088 km). That is a spherical approximation of the WGS 84 ellipsoid, off by up to about 0.5%; use an ellipsoidal geodesic if that matters.

A position normalises its longitude, but a box is kept as given. The plain box operations do rectangle arithmetic on degrees, so a box with west > east is treated as empty and never wraps across the antimeridian. The distance queries do wrap. Box.around, Index and Segments handle the antimeridian and the poles exactly, because a nearest neighbour search that stopped at ±180° would return wrong answers near it.

Exact nearest-neighbour queries #

An R-tree (nox-rtree) answers box queries. A distance query guesses a radius, asks the R-tree for everything in the boxes that cover the disc of that radius, and doubles the radius until the answer is complete. Since Box.around covers the whole disc, any point outside those boxes is further away than the radius. Once the k-th neighbour found lies within the radius, no nearer point can remain.

This only works if Box.around really contains the disc, so it uses the exact spherical bound asin (sin σ / cos φ), where σ is the disc's angular radius (its radius in kilometres over the Earth's mean radius R = 6371.0088 km) and φ the latitude of its centre. That is the half-width in longitude of the smallest box around the spherical cap, reached where a meridian is tangent to the cap; once the cap reaches a pole (sin σ ≥ cos φ) the box becomes a full band of longitude. A local degrees-per-kilometre scale, σ / cos φ, falls short of it by about R σ³ tan² φ / 6 on the ground: 12 m for a 100 km disc at 60° latitude, 1.6 km for a 500 km one. The bound is checked against S2's Cap.get_rect_bound in test/interop/python.

  • geocaml/geo provides geospatial primitives over Owl, and is the reason this package is not called geo. It is planar: its distance is Euclidean and its only haversine is the initial bearing. It has no index, and its rectangles have no subtraction.
  • rtree is the R-tree this indexes on, forked as nox-rtree. It answers window queries; the nearest-neighbour layer on top is here.
  • geojson is an RFC 7946 codec, which reads the coordinates this library measures.
  • ocaml-proj provides PROJ bindings, for moving between coordinate reference systems. This library stays in WGS 84 degrees and never reprojects.
  • h3 provides bindings to Uber's hexagonal index, which does handle the whole sphere uniformly, at the cost of a C dependency.

License #

ISC. See LICENSE.md.