zig-geo

Positions on the Earth for Zig 0.17, and the questions programs ask about them:

• Is this place inside that area?
• How far apart are these two places, and in which direction does one lie from the other?
• What box contains all of these?

The library has no file format, no projection, no allocation and no dependencies. zig-geojson reads and writes GeoJSON in terms of its Position. Anything else that has coordinates can use the functions directly, because they take slices.

The API documentation is generated from the doc comments, which carry most of the explanation. The workflow publishes it from main to https://jeff.jcollie.page/zig-geo/.

Where this lives

The repository is in three places, all with the same history. The Forgejo instance at https://git.jcollie.dev/jeff/zig-geo is the one with a web view:

$ git clone https://git.jcollie.dev/jeff/zig-geo.git

It is mirrored on Tangled at https://tangled.org/jcollie.dev/zig-geo. It is also on the Radicle network, with the repository identifier

rad:zWWHzy9XC7RCpFsr2ykcfe4NMGX6

rad clone rad:zWWHzy9XC7RCpFsr2ykcfe4NMGX6 fetches it from any node that seeds it. Any of the three is the whole project.

Adding it to a project

$ zig fetch --save git+https://git.jcollie.dev/jeff/zig-geo

The repository is zig-geo, but the package is named geo, so that is the name zig fetch --save records it under in build.zig.zon, and the name the build script asks for:

const geo = b.dependency("geo", .{ .target = target, .optimize = optimize });
exe.root_module.addImport("geo", geo.module("geo"));

main follows Zig 0.17. For Zig 0.16, take the zig-0.16 branch or the v0.1.0 tag:

$ zig fetch --save git+https://git.jcollie.dev/jeff/zig-geo#v0.1.0

Using it

const geo = @import("geo");

const des_moines: geo.Position = .{ .longitude = -93.625, .latitude = 41.5868 };
const chicago: geo.Position = .{ .longitude = -87.6298, .latitude = 41.8781 };

// On the WGS 84 ellipsoid, to the nanometre.
const wgs84 = geo.Ellipsoid.wgs84;
const metres = wgs84.distance(des_moines, chicago);
const heading = wgs84.initialBearing(des_moines, chicago); // degrees from north
const halfway = wgs84.destination(des_moines, heading, metres / 2);

// On a sphere, which is quicker and within about half a per cent.
const roughly = geo.sphere.distance(des_moines, chicago);

// Coordinates on NAD27, measured on its ellipsoid.
const nad27_metres = geo.Ellipsoid.epsg.clarke_1866.distance(des_moines_nad27, chicago_nad27);

// A polygon: an exterior ring, then any holes.
const warned: []const []const geo.Position = &.{exterior};
if (geo.polygon.contains(warned, des_moines)) { ... }

// Cheap rejection before the exact test.
const box = geo.BoundingBox.of(exterior).?;
if (!box.contains(des_moines)) { ... }

Namespace	What it answers
Position	a place: longitude and latitude, on WGS 84 unless the caller says otherwise, and an optional altitude
BoundingBox	the box around some positions; containment and intersection; merging two boxes, taking the short way across the antimeridian
ring	whether a ring is closed, which way it runs, whether it contains a point
polygon	containment with holes, multipolygons, RFC 7946's right-hand rule
Ellipsoid	distance, initial bearing and destination along the ellipsoid's geodesics; every EPSG ellipsoid; their radii, and the spheres made from them
Geodesic	the same in bare latitudes and longitudes, with the arrival azimuth: Karney's algorithm
sphere	great-circle distance, initial bearing, and the destination given both, on a sphere

Models of the Earth

Each question is answered in the model it belongs to:

• Areas use the plane of longitude and latitude. ring, polygon and BoundingBox treat longitude and latitude as x and y. That is how GeoJSON defines a polygon's edges: straight lines in coordinate space (RFC 7946 §3.1.1). An area drawn in GeoJSON is therefore exactly the area these functions test against. A ring that crosses the antimeridian has to be cut there, as §3.1.9 asks.
• Distances on the ellipsoid use its geodesics. Ellipsoid.distance, initialBearing and destination follow the shortest path along the ellipsoid itself. They use Karney's algorithm, ported from his own C in Geodesic, which is accurate to about 15 nm. Unlike Vincenty's method, it converges for every pair of points, including nearly antipodal ones.
• Distances on a sphere are the quick approximation. sphere measures on a sphere, by default of the Earth's mean radius, GRS 80's R₁ = 6 371 008.7714 m. That is within 0.56 % of the ellipsoid everywhere: plenty for "how far away is the storm", but not surveying. The worst case is a short north–south distance at the equator, where the ellipsoid curves more tightly than the sphere.

The geodesic port

Geodesic is a line-by-line port of geodesic.c, Karney's C implementation, as PROJ 9.9.0 ships it. It keeps the C's names and comments so the two can be read side by side. It covers the inverse and direct problems; area, reduced length, geodesic scale and arc mode are left out. Its tests are the 20 reference geodesics from Karney's published test set and PROJ's regression cases for the hard ones: nearly antipodal points, the poles, prolate ellipsoids and NaN inputs.

It is also checked against the C itself on every push. reference/ is a project of its own that fetches Karney's C, geographiclib-c v2.2, and compiles it beside the port. Its geodesic.c and geodesic.h are byte for byte the files the port was made from. Keeping it separate means nothing that only builds the library has to fetch C or compile it.

$ cd reference
$ nix develop ..#default -c zig build run -- --iterations 2000000   # what CI runs
$ nix develop ..#default -c zig build run -- --seconds 600
$ nix develop ..#default -c zig build run -- --seed 12345           # exactly again

Each case goes to both implementations, over every EPSG ellipsoid, a sphere, and three shapes far from the Earth's. The generated inputs are weighted toward the hard cases: nearly antipodal points, very short lines, the poles and the equator, whole degrees, Karney's regression cases nudged by an ulp or two, and coordinates out of range. The two must agree:

• Distances to 2e-8 m. The worst seen in 200 million cases is 1.1e-8 m, three ulps of 20 000 km.
• Azimuths to 1e-8°. Where they differ by more, the C has to confirm the port's azimuth by following it. At the start of a line it must be another shortest path, which happens for lines a few nanometres long and for antipodal points. At the end of a line it must lead back to the start, which near a pole it can do while differing a good deal.
• Arrivals to 1e-7 m.
• Out-of-range input must give NaN from both. In-range input must always get an answer from the port.

They are not bit for bit the same, and can't be: Zig's sin, cos and atan2 round differently from the C library's. A failure prints its seed and a command that replays that one case exactly.

A departure from the C

The reference check found a case where Karney's C gives NaN for valid input. It happens on a prolate ellipsoid, for a short, nearly east–west line between latitudes one ulp apart. Rounding can put the two latitudes' cosines the wrong way round, and a square root then gets a negative argument. In random testing with glibc, about one such line in 9 000 fails. geodesic.c already guards the case where the two come out equal, and the port extends that guard to the reversed case, where the latitudes are also equal to within rounding. The reference check confirms the port's answer by having the C follow it. Upstream's main branch still has this. No Earth ellipsoid is prolate, so Ellipsoid users are never affected.

The sphere's formulas

Two of the sphere's formulas are not the usual ones, and fuzzing is the reason:

• distance uses the spherical case of Vincenty's formula, an arctangent of the angle's sine and cosine, which is well conditioned at every distance. The first version used the haversine formula. That loses precision for nearly antipodal points, and the fuzzer found a detour through a third point that measured shorter than the direct route. The two arguments are also put in a fixed order, so distance(a, b) equals distance(b, a) to the last bit.
• destination works with unit vectors in three dimensions, as in Gade's n-vector representation. The textbook formula finds the new latitude as an arcsine, which loses half its digits near the poles. Setting off for the North Pole from 48 km away, it arrived 9.5 cm short.

Other ellipsoids

Every Ellipsoid measures on itself with distance, initialBearing and destination. It can also stand in a sphere for itself: a sphere.Sphere is just a radius, and Ellipsoid.sphere makes one from one of its radii. The free functions in sphere measure on sphere.mean, which is the same as Ellipsoid.grs80.sphere(.mean).

Ellipsoid.epsg holds every Earth ellipsoid in the EPSG Geodetic Parameter Dataset that EPSG hasn't deprecated, 50 of them. Each is named after EPSG's own name, so Clarke 1880 (RGS) is epsg.clarke_1880_rgs. epsg.all lists them with their codes. Ellipsoid.fromEpsg(7008) finds one by code, and Ellipsoid.fromName("clarke 1866") finds one by EPSG's name, ignoring case, through the comptime map Ellipsoid.by_name. Ellipsoid.wgs84 and Ellipsoid.grs80 are shortcuts to the two most used.

For coordinates on	Use
WGS 84, GPS, GeoJSON	Ellipsoid.wgs84
NAD83, ETRS89, GDA2020, NZGD2000, ITRF	Ellipsoid.grs80
NAD27	epsg.clarke_1866
ED50	epsg.international_1924
OSGB36	epsg.airy_1830
DHDN, MGI, Tokyo	epsg.bessel_1841
Pulkovo 1942, Beijing 1954	epsg.krassowsky_1940

The radii are:

• R₁, mean: (2a + b) / 3.
• R₂, authalic: the sphere with the ellipsoid's surface area.
• R₃, volumetric: the sphere with the ellipsoid's volume.
• a, equatorial, and b, polar.

For a sphere, the choice of radius matters less than it looks. For one ellipsoid, R₁, R₂ and R₃ agree to within 8 m, about one part in a million. Only a and b are really different: they are a third of a per cent apart. Across all of EPSG's ellipsoids, the mean radii span 2 040 m, from Zach 1812 to International 1924, or one part in three thousand. WGS 84 and GRS 80 differ by less than a tenth of a millimetre.

Choosing an ellipsoid doesn't convert coordinates between datums. A datum is an ellipsoid and where it sits against the Earth. The same place can have NAD27 and NAD83 coordinates up to about a hundred metres apart, and no choice of radius changes that. Picking the matching ellipsoid only measures in the frame the coordinates were given in. Converting between datums is a job for PROJ.

Where epsg comes from

src/epsg_ellipsoids.zig is generated from proj.db, the copy of the EPSG dataset that ships with PROJ. The devshell provides PROJ and Python and sets $PROJ_DB to the database. When PROJ is updated, regenerate the file and commit it:

$ nix develop -c python3 tools/epsg_ellipsoids.py -o src/epsg_ellipsoids.zig

The file's header records which PROJ and which EPSG release it came from.

The EPSG dataset belongs to the International Association of Oil & Gas Producers (IOGP). The generated file is licensed under EPSG's terms of use, which are in LICENSES/LicenseRef-EPSG.txt, and not under MIT. The terms allow the data to be copied and redistributed, but not for profit. They require IOGP's ownership to be acknowledged and the terms to be passed on to anyone who receives the data. They allow exactly two changes to a parameter: converting between a/b and a/1/f, and converting units with EPSG's own conversion factors. So each ellipsoid is written the way EPSG defines it, and a length given in Clarke's feet or Indian feet is multiplied by EPSG's factor for that unit.

Building and testing

Everything runs inside the flake's devshell:

$ nix develop
$ zig build test --summary all        # unit tests, and the fuzz targets on their seeds
$ zig build check                     # also compiles the fuzz driver and the docs server
$ zig fmt --check --exclude zig-pkg --exclude reference/zig-pkg .
$ reuse lint
$ nix build                           # the package, built and tested offline

Fuzzing

tests/fuzz.zig decodes each input as places: pairs of u32 scaled onto the valid ranges. The targets check these properties:

• ellipsoid: on WGS 84 and International 1924, the same properties as the sphere's, and the mean sphere's distance within 0.56 % of the ellipsoid's.
• sphere: on the default sphere and on the smallest and largest of the others, distances are finite, bounded, symmetric and obey the triangle inequality; bearings fall in [0, 360); and destination undoes distance and initialBearing.
• ring: containment doesn't depend on which vertex a ring starts at, or on whether the ring is written closed.
• bbox: a box contains the points it was made from, and a merge of two boxes contains everything either one did.

Zig's own fuzzer runs them with coverage feedback:

$ zig build fuzz --fuzz                                # until interrupted
$ zig build fuzz --fuzz=1M                             # a bounded run

The bounded report names only the first fuzz test in the binary, but every target ran. tools/fuzz.zig is a mutation loop of its own, without coverage, that CI runs because it is bounded by a count, seeded, and replays a single input exactly:

$ zig build fuzz-run -- --iterations 60000             # what CI runs
$ zig build fuzz-run -- --seconds 600 --target sphere
$ zig build fuzz-run -- --input fuzz-findings/x.bin --target sphere

API documentation

$ zig build docs         # into zig-out/docs
$ zig build docs-serve   # and read them at http://127.0.0.1:8000/

The docs have to be served rather than opened as files. The viewer is a WebAssembly program that fetches its sources at runtime, and a browser refuses that from a file:// page.

References cited

Butler, H., M. Daly, A. Doyle, S. Gillies, S. Hagen, and T. Schaub. 2016. The GeoJSON Format. RFC 7946. RFC Editor. https://www.rfc-editor.org/info/rfc7946. §3.1.1 is why edges are straight in coordinate space, §3.1.6 is the right-hand rule, §3.1.9 is the cut at the antimeridian, and §5.2 is how a box across it is written.

Franklin, W. Randolph. n.d. "PNPOLY - Point Inclusion in Polygon Test." Accessed September 24, 2026. https://wrfranklin.org/Research/Short_Notes/pnpoly.html. The even–odd crossing test in ring.contains, including the half-open comparison of latitudes that counts a ray through a vertex once.

Gade, Kenneth. 2010. "A Non-Singular Horizontal Position Representation." Journal of Navigation 63 (3): 395–417. https://doi.org/10.1017/s0373463309990415. Positions as unit vectors, which is how sphere.destination avoids the arcsine that loses precision near the poles.

IOGP. n.d. "EPSG Geodetic Parameter Dataset." International Association of Oil & Gas Producers. Accessed September 24, 2026. https://epsg.org/. The source of Ellipsoid.epsg, generated from release v13.102 as PROJ 9.9.0 ships it in proj.db; the datums listed against each ellipsoid come from the same place.

Karney, Charles F. F. 2013. "Algorithms for Geodesics." Journal of Geodesy 87 (1): 43–55. https://doi.org/10.1007/s00190-012-0578-z. The algorithm in Geodesic, which is ported from the author's C implementation, geodesic.c, as PROJ 9.9.0 carries it (MIT licensed). Its tests are the reference geodesics and regression cases of that implementation's geodtest.c.

Karney, Charles F. F. 2022. "Addenda and Errata for Papers on Geodesics." GeographicLib. June 8, 2022. Accessed September 24, 2026. https://geographiclib.sourceforge.io/geod-addenda.html. What geodesic.c does beyond the 2013 paper, and so what Geodesic does: it takes a Newton step in the direct problem when |f| > 0.01, brackets the Newton iteration of the inverse problem, multiplies the A₂ series by 1 + ε instead of dividing it by 1 − ε, and breaks the degeneracy of lines along the equator.

Karney, Charles F. F. 2025. GeographicLib C Library. Version 2.2. C. https://github.com/geographiclib/geographiclib-c. The implementation reference/ checks Geodesic against, whose geodesic.c is the file it was ported from, and whose NaN for some valid inputs on a prolate ellipsoid is the one place the port departs from it.

Moritz, H. 2000. "Geodetic Reference System 1980." Journal of Geodesy 74 (1): 128–33. https://doi.org/10.1007/s001900050278. The source of sphere.earth_mean_radius, the ellipsoid's mean radius R₁, and the GRS 80 values Ellipsoid's tests check against. Its authalic radius R₂, 6 371 007.1810 m, is 0.12 mm higher than the closed form gives when worked to fifty digits.

Vincenty, T. 1975. "Direct and Inverse Solutions of Geodesics on the Ellipsoid with Application of Nested Equations." Survey Review 23 (176): 88–93. https://doi.org/10.1179/sre.1975.23.176.88. The inverse formula whose spherical case is sphere.distance. On the ellipsoid, Geodesic uses Karney's method instead, because Vincenty's iteration fails to converge for nearly antipodal points.

License

MIT, except src/epsg_ellipsoids.zig, which is EPSG data under EPSG's terms of use (see Where epsg comes from). The project is REUSE compliant: every file says what it is licensed under, and reuse lint checks that.