Make the projection system generic (phases 1–5 + follow-ups) by trasch · Pull Request #250 · Outdooractive/gis-tools

trasch · 2026-09-15T10:48:33Z

What it delivers

1. ProjectionKind + capability flags (#237)

ProjectionKind (.geographic / .planar / .geocentric / .undefined) plus convenience flags, wraparoundExtent and crsLength(fromMeters:). Algorithms and antimeridian handling dispatch on the semantic category instead of EPSG codes; the duplicated 4326/3857 antimeridian helper implementations collapsed into shared parameterized ones.

2. ProjectionDefinition registry (#238)

Internal ProjectionDefinition strategy protocol with forward/inverse relative to the EPSG:4326 pivot; ProjectionRegistry maps projections to definitions. noSRID keeps its verbatim-copy semantics.

3. New projections

EPSG:3395 (World Mercator, ellipsoidal Snyder forward + iterative inverse) and EPSG:32662 (Plate Carree) — Add EPSG:3395 and EPSG:32662 projections #239
All 120 UTM zones (EPSG:32601–32660 / 32701–32760), SRID-range init?(srid:) — Add UTM zone projections (EPSG:32601–32660, 32701–32760) #241
Datum CRSs — Ship datum support: ETRS89, NAD27, OSGB 1936 and British National Grid (#243) #249: EPSG:4258 (ETRS89, identity), 4267 (NAD27), 4277 (OSGB 1936), 27700 (British National Grid) on Airy 1830, all verified against pyproj ground-truth fixtures computed before porting:
- public Datum / Ellipsoid metadata model (Projection.datum)
- public HelmertTransformation (7-param, position-vector, EPSG "NAD27 to WGS 84 (4)" ~10 m / "OSGB 1936 to WGS 84 (6)" EPSG:1314 ~2 m), exact rotation matrices, arbitrary source ellipsoids, both directions with explicit result labeling
- internal TransverseMercatorMath(ellipsoid:) — Snyder TM parameterized with a latitude of origin; UtmDefinition delegates to it (UTM suite proves behavior preserved)
- grid shifts (NADCON/OSTN15, sub-meter) deliberately out of scope → tracked as Grid shift file support (NTv2 / NADCON / OSTN15) for sub-meter datum accuracy #248 (load-don't-bundle)

4. Registry-driven WKT matching + coder/tile tests (#240)

wktMatchers: [[String]] alternative fragment sets on each definition; init?(wkt:) probes in registry order (specific before generic). Coders' per-projection nested switches collapsed to generic resolution. Ordering pins: Pseudo-Mercator (3857) and British National Grid (27700) before the generic "Mercator" (3395) fragments — a BNG .prj contains Transverse_Mercator.

UTM .prj files identified via regex with Central-Meridian cross-validation (Identify UTM zones in WKT strings via regex #244); UTM-zone-token strings short-circuit to UTM-or-nothing

5. Struct-based projection + add-only custom registration (#247)

This is a breaking change (public enum → struct), designed to keep source impact minimal:

125 same-named static let constants — dot-shorthand (x == .epsg4326) compiles unchanged; Codable keeps the SRID wire format
Projection.register(CustomProjection) -> Bool — add-only (no unregister/disable by design); public struct CustomProjection value descriptor with @Sendable forward/inverse closures relative to the EPSG:4326 (WGS84) pivot — the documented seam where datum transformations live (Support CRSs on other datums (NAD27, OSGB, CH1903+, ...) #243)
public struct ProjectionExtent replaces the tuple
Performance: hot paths (projected(to:), kind, clamped(), isValid, wraparoundExtent) read the definition captured in the Projection value — zero registry access, no locking. Release benchmarks: full path ~1.7× direct-math baseline (was 35× before definition capture; a parallel test run melted down at 924 s per-run from per-access Mutex reads, which motivated the capture)

6. Batch coordinate conversion (commit `376a561`)

Array<Coordinate3D>.projected(to:) + in-place project(to:): geometries and the WKB/WKT/TWKB coders now batch-convert whole coordinate lists with the conversion setup (TM parameterization, origin constants, Helmert rotation matrices) hoisted once per batch instead of per coordinate. Math is identical — the batch machinery only extracts constants through an internal prepared representation, proven output-identical by the frozen fixture suite across all CRS families. Release benchmarks: −15% (4326→3857), −9% (4326→UTM→4326); BNG lane flat (two-step chain dominates). Semantics: uniform-projection batches only (library invariant, asserted); noSRID values verbatim in both directions.

7. Validation API + UTM convenience (#245)

Coordinate3D.isValid / BoundingBox.isValid against the definition's extent; Projection.utmZone/utmHemisphere, Projection(utmZone:hemisphere:), Projection(utmZone:)-style factories.

8. Swiss and Irish CRSs (commits `6c9b6f7`)

New built-in CRSs, all with EPSG Helmert datum chains and prepared batch transforms:

EPSG:2056 (CH1903+/LV95) / EPSG:21781 (CH1903/LV03) — Swiss oblique Mercator on Bessel 1841 ("CH1903+ to WGS 84 (1)" / EPSG:1676, ~1 m). SwissObliqueMercatorMath is a faithful port of PROJ's somerc.cpp
EPSG:29902 (TM65/Irish Grid) / EPSG:29903 (TM75/Irish Grid) — TM on Modified Airy ("TM65 to WGS 84 (2)" / EPSG:1641, ~1 m)
EPSG:2157 (IRENET95/ITM) — TM on GRS80, ≈ WGS84 identity
Datum model gains Ellipsoid.modifiedAiry/.bessel1841 and Datum.ch1903plus/.ch1903/.tm65/.tm75/.irenet95

9. UTM zone selection (commit `6c9b6f7`)

UtmDefinition.zone(for:) + Projection.utmZone(for:) return the zone for a WGS84 coordinate, implementing the EPSG banding exceptions: Norway (56–64°N, 3–12°E → zone "32V") and Svalbard (72–84°N → 3°-wide zones 31/33/35/37 and the 9°-wide zone 32). Antimeridian-normalized; nil outside −80°–84° (UPS territory).

10. Pan-European and national CRSs (commit `4d2a332`)

Three new math implementations — faithful ports of PROJ 9.8.1, each validated against pyproj ground truth to sub-millimeters before landing:

LambertConformalConicMath (Snyder 15-x + GeographicLib tauf Newton inverse)
LambertAzimuthalEqualAreaMath (all four laea domains incl. the authalic latitude series)
ObliqueStereographicMath (the Gauss-conformal reduction + sterea)

New CRS definitions:

EPSG:25831–25837 ETRS89/UTM belt — TM on exact GRS80 (values differ from the WGS84 zone blocks by decimeters through the flattening)
EPSG:3035 ETRS89-LAEA Europe and EPSG:3034 ETRS89-LCC Europe (the EU-wide statistical/ensemble mappings)
EPSG:2154 RGF93 / Lambert-93, the official French grid
EPSG:28992 Amersfoort / RD New — Helmert "Amersfoort to WGS 84 (4)" + oblique stereographic on Bessel 1841 (Datum.amersfoort)

11. Karney transverse Mercator for UTM (commit `ce7498f`)

Replaces the Snyder formulas in UtmDefinition/Etrs89UtmDefinition with Karney's transverse Mercator (the Krüger series extended to 6th order, arXiv:1002.1417; implemented from the paper, no vendored code) under the unchanged ProjectionDefinition contract — conformal-sphere mapping, Clenshaw-summated Krüger series, Newton iteration for the conformal-latitude inversion. The Snyder TransverseMercatorMath stays for the other TM-based CRSs (BNG, Irish Grids). Composes with future TM-family CRSs (Gauss-Krüger, State Plane TM).

Verification (independent Python implementation first, validated before porting):

Against Karney's published test set TMcoords.dat (287k points, 80-digit exact reference): worst forward 5.7 nm / reverse 5.0 nm within 3900 km of the central meridian — consistent with the paper's 5 nm bound
Against PROJ 9.8 (GeographicLib-based): worst 7.5 nm over 1440 random points across all 120 WGS84 zones; worst 5.6 nm over the GRS80 ETRS89 belt
Round-trip sweep over all zones, latitudes incl. ±89.999° and the poles: worst 6.4 nm (the issue required < 1 µm)
Snyder baseline divergence quantified: ~0.96 mm in-zone (Snyder's meridian-arc truncation), 11 mm at the ±6° zone borders, ~1 km at 45° off the central meridian — the degradation that motivated the swap

Behavior change to acknowledge (sub-millimeter, per #246): Snyder and Karney disagree at the sub-millimeter level in-zone; both are valid UTM within the EPSG-defined accuracy. The UtmAlgorithmTests reference values were re-derived under Karney with the tolerance tightened from 1 mm to 1 µm (e.g. zone 19N easting 331792.115 → 331792.114806 m; zone 20S northing shifts 0.6 mm). Downstream users re-projecting Snyder-era UTM values should expect millimeter-scale shifts.

Performance (release, 10k coordinates, median): UTM round trip registry path 2.71 → 7.6–8.7 ms (~2.8–3.2×; above the 1.5–2× estimate in the issue — the extra ~500 ns/round trip buys nanometer accuracy and guaranteed convergence), batch path 2.47 → 7.4–7.8 ms. Unaffected paths (3857/4978/27700) unchanged within noise.

12. Projection matrix cross-validation (commit `8cb37323`)

New ProjectionMatrixTests: every sensible pair of built-in projections is cross-validated against reference data pre-computed with pyproj/PROJ 9.8.1. The generator (scripts/generate_projection_matrix.py, committed) uses explicit parameterized proj4 strings mirroring the library's parameters — including the Helmert datum transformations in position-vector form — never EPSG codes, so PROJ cannot silently select a better transformation than the library implements. Data (3612 rows in TestData/ProjectionMatrix/matrix.csv) and test code are separate; 459 points span global coverage, regional centers/corners and all 120 UTM zones; 26,158 ordered pairs assert x/y values, the projection label, round trips and batch-API equivalence. Tolerances are measured ceilings (audited: exact-CRS pairs ≤ 39 µm, UTM ≤ 1e-7 m). The known deviations were both fixed by follow-up commits: #252 (laea legacy authalic series) by commit d6fdc0e — removing the EPSG:3035 budget — and #251 (Helmert z ordering) by commit b482f28 — the Helmert height passthrough brought the 176 EPSG:4978 ↔ datum-CRS pairs into the matrix (26,334 pairs). EPSG:32662 remains excluded (the library intentionally stores degrees as x/y).

13. LAEA authalic latitude upgraded to the auxlat series (commit `d6fdc0e`)

The laea inverse now converts the authalic latitude through the Karney auxlat series (PROJ 9.8's pj_auxlat_coeffs, order 6 in n, Clenshaw summation) instead of the legacy 3-term Snyder series: full double precision for |f| ≤ 1/150 vs the legacy ~1 mm loss at large distances. Validated in Python against PROJ before porting (Cairo: 1.2e-8 → 2e-14 deg); measured over the matrix points: inverse worst 0.6 µm, round trip 0.5 µm. The EPSG:3035 pair budget is removed and the fixture regenerated; a new domain-wide round-trip test pins the series at 1 µm / 1e-9 deg. Closes #252. (LambertConformalConicMath checked per the issue: it uses the conformal-latitude Newton iteration, no shared pattern.)

14. Helmert height passthrough (commit `b482f28`)

HelmertTransformation now preserves the input height in both directions (the 3D transform still runs on the ECEF position; only the output z is the input height), matching PROJ's +towgs84 nominal-height semantics — validated in Python against PROJ before porting. This fixes the pivot inconsistency: 4326 → datum CRS → 4978 equals the direct 4326 → 4978 instead of carrying the datum's vertical offset (~26 m for Amersfoort). All 176 EPSG:4978 ↔ datum-CRS pairs are now verified in the projection matrix (26,334 pairs; measured worst 48.9 mm, within the documented 0.1 m ordering-difference budget); the exclusion machinery is gone. The Helmert round trip loosens to ~1e-8 deg — the same looseness PROJ's own +towgs84 exhibits. Closes #251.

15. German and Austrian Gauss-Krüger zone belts (commit `1f31f9d`)

Two new zone-belt definitions on Bessel 1841 with the Karney transverse Mercator: EPSG:31466–31469 (DHDN / Gauss-Krüger zones 2–5, the German cadastral grid, x₀ = zone·1M + 500k) and EPSG:31255–31259 (MGI / Austria Gauss-Krüger, all five zones: the x₀=0 Central/East variants plus M28/M31/M34 with the Ferro-based numbering). New Datum.dhdn ("DHDN to WGS 84 (2)", ~3 m) and Datum.mgi ("MGI to WGS 84 (2)", ~1.5 m). PROJ fixtures computed before finalizing (34 forward values across all 9 zones); one implementation bug caught en route (the first pass omitted the Helmert — 52–89 m off). Documented 1.9 mm residual on MGI: PROJ applies datum shifts as a geodetic-domain approximation, the library as the exact XYZ affine (invisible for DHDN's 0.2″ rotations); covered by a 2 mm datum budget. New GaussKruegerWktIdentification (like the UTM short-circuit) handles GDAL/ESRI names, the Ferro M-tokens, and parameter-based zone resolution with CM/False_Easting cross-validation. The projection matrix now covers all 9 GK CRSs (3,626 rows, 26,700 pairs); 9 new GaussKruegerTests.

17. North American NAD83 CRSs (commit `06a58ae`)

New Datum.nad83 (identity over WGS84, like ETRS89; GRS80) and six CRS families: EPSG:4269 (geodetic, identity), EPSG:26901–26960 (NAD83/UTM belt, Karney TM on GRS80), EPSG:5070 (Conus Albers — the USGS/EPA/NLCD standard) and EPSG:3005 (BC Albers) via a new AlbersEqualAreaMath (Snyder 14-1…14-14), and EPSG:3347/EPSG:3978 (Canadian Lambert conic grids) via the existing LCC math. The authalic-latitude conversion is extracted into a shared AuthalicLatitude helper used by both laea and albers — the Albers inverse matches PROJ 9.8 exactly (auxlat series) and avoids the near-pole Newton stall a naive q-iteration exhibits. Python-first verification caught two transcription bugs before porting; the final 3,000-point sweep measures forward 33 nm / round trips 3.7 µm. WKT identification handles NAD83 datum tokens and routes NAD83 UTM zone tokens to the 269xx belt. The projection matrix now covers 4,409 rows / 41,290 pairs; 6 new Nad83CrsTests.

After merge

Will close Public API for third-party / custom projections #242, Support CRSs on other datums (NAD27, OSGB, CH1903+, ...) #243, Replace Snyder transverse Mercator in UTM with Karney's algorithm #246, Helmert datum transformations shift the height: 4326 -> datum CRS -> 4978 differs from the direct 4326 -> 4978 #251 and LAEA (EPSG:3035) inverse loses ~1-2 mm at large distances: upgrade the authalic-latitude inversion to the Karney auxlat series #252 on merge.
Follow-ups that stay open after this merge (not implemented by this PR):
- Grid shift file support (NTv2 / NADCON / OSTN15) for sub-meter datum accuracy #248 (grid shift files, sub-meter datum accuracy) — deliberately out of scope here
The stacked PRs Introduce ProjectionKind and capability-based projection dispatch #237–Ship datum support: ETRS89, NAD27, OSGB 1936 and British National Grid (#243) #249 are closed as superseded by this PR

Phase 1 of making the projection system generic: replace per-EPSG-code branches with dispatch on a semantic projection category, so that new projections registered in a later phase automatically get sensible default behavior. New API: - ProjectionKind enum (geographic / planar / geocentric / undefined) - Projection.kind, isGeographic, isPlanar, isGeocentric, hasSRID - Projection.wraparoundExtent (±180° for 4326, ±originShift for 3857) - Projection.crsLength(fromMeters:) for meter-to-CRS-unit conversion Refactored to kind/capability dispatch (behavior-neutral): - Geodesic ops (distance, bearing, destination, midpoint, rhumb *) - AntimeridianCutting: duplicated 4326/3857 helpers collapsed into implementations parameterized by wraparoundExtent - BoundingBox date-line splits, padding, expansion, size, normalize - Antimeridian-span checks, unit conversions, grid cell sizing across ~20 algorithm files Also adds a .swiftformat configuration file. Per-CRS switch matrices (Coordinate3D.projected, coders, MapTile) are intentionally unchanged; they become registry data in phase 2.

Phase 2 of making the projection system generic: extract the inline transform math into a ProjectionDefinition strategy protocol and route all conversions through a registry, keeping EPSG:4326 as the pivot. - ProjectionDefinition protocol: forward/inverse relative to the EPSG:4326 pivot, plus validExtent and worldBoundingBox metadata - Definitions for EPSG:4326 (identity), EPSG:3857 (Web Mercator), EPSG:4978 (ECEF, with the geodetic/ECEF helpers moved here) and noSRID (null math) - ProjectionRegistry maps a Projection to its definition - Coordinate3D.projected(to:) now dispatches through the registry; clamped() uses the definition's validExtent; latitudeProjected/ longitudeProjected derive from the full projection - BoundingBox.expanded(byDistance:) and Random's world box use the registry metadata instead of per-CRS switches noSRID keeps its historical semantics, now documented and pinned by tests: verbatim copy to .noSRID, verbatim relabel to EPSG:4326 and EPSG:3857 (planar-meter convention), interpreted as EPSG:4326 into EPSG:4978. New tests: round trips across all projection pairs, pivot routing, noSRID semantics, per-axis helper consistency, altitude/m preservation and clamped() extents.

Phase 3 of making the projection system generic: prove the registry architecture by adding two new WGS84 projections. - EPSG:3395 (WGS 84 / World Mercator): ellipsoidal Mercator with Snyder forward formulas and an iterative inverse (the conformal latitude iteration). Wraps at ±originShift like EPSG:3857, but y values differ away from the equator. Valid y extent is the EPSG-registered bound of ±20_048_966.104014604. - EPSG:32662 (WGS 84 / Plate Carree): identity x/y in degrees, wraps at ±180°. Both are registered with SRID lookup, description, WKT patterns ("Pseudo-Mercator" vs. "Mercator" now disambiguated between 3857 and 3395), planar ProjectionKind and wraparound extents. All algorithms, BoundingBox, MapTile and the pixel conversions pick them up automatically through the ProjectionKind/registry abstractions - no per-algorithm changes needed. The WKB/WKT/TWKB coordinate constructors collapsed their per- projection nested switches into one generic expression (behavior- identical for all existing cases); new projections now decode and encode correctly with embedded SRIDs. Tests: conversion values computed independently with the Snyder formulas, round trips across all five projections, antimeridian cutting, distance/boundingBox/contains, normalize/clamp, destination/bearing via the EPSG:4326 pivot, and WKB round trips in the new projections.

Phase 4 of making the projection system generic: move WKT CRS identification into the projection registry and round out test coverage for the coder and MapTile paths. - ProjectionDefinition gains wktMatchers: [[String]] - each inner list is an alternative fragment set, every fragment of one set must be contained for a match. This preserves the previous A && (B || C) semantics (e.g. "WGS 84" or "WGS_1984") exactly. - ProjectionRegistry probes definitions in registry order (specific patterns first: EPSG:3857 "Pseudo-Mercator" before EPSG:3395 "Mercator"); Projection.init?(wkt:) now delegates to it. noSRID registers no patterns and never matches. - MapTile tests: tile bounding boxes in EPSG:3395/EPSG:32662 verified against manual reprojection of the EPSG:4326 bounds, plus centerCoordinate round trips (the generic pixel-projection path already handles the new projections). - Coder tests: WKT round trips with SRID embedded in the string, TWKB decode fixtures for the new projections (decode-only API), and ProjectionWktMatcherTests covering the 3857/3395 disambiguation, both geographic fragment alternatives and unrecognized strings (GEOGCS["ETRS89"] -> nil). - Test rename: NewProjectionAlgorithmTests.swift split into projection-named files (Epsg3395Tests.swift, Epsg32662Tests.swift) so the names stay meaningful as more projections are added. Behavior note: WKTCoder.decode with an explicit sourceSrid does not consume an embedded "SRID=...;" prefix - decoders of strings produced by WKTCoder.encode must use sourceSrid: nil so the SRID is detected from the string.

Phase 5 of making the projection system generic. Math: standard Snyder transverse Mercator formulas on the WGS84 ellipsoid, verified independently in Python before porting - the meridian arc against Simpson integration of the meridian radius (<0.5mm divergence) and round trips over 252 points covering all zones and latitudes to 89 degrees (max error 1.06e-08 degrees). - UtmDefinition.swift: UtmHemisphere and UtmDefinition implementing both directions of the projection with a per-zone central meridian, the 0.9996 scale factor, the 500_000 false easting and per-hemisphere false northing. Sector extents follow the EPSG convention (easting 100_000-900_000, northing 0-10_000_000); the world bounding box is the zone sector and no WKT pattern is registered (UTM .prj files would need per-zone parsing). - Projection gains all 120 zone cases (northern 32601-32660, southern 32701-32760) and init?(srid:) accepts the contiguous SRID ranges. wraparoundExtent is nil (zone coordinates never wrap), kind is planar; description now defaults to the EPSG label so it scales with the SRID count. - ProjectionRegistry resolves definitions via an eagerly built dictionary; registration coverage is asserted at runtime by the new UtmDefinitionTests sweep over all UTM SRIDs. Tests: SRID-to-zone/hemisphere mapping for all 120 SRIDs, registry coverage sweep, central meridian math per zone, known forward values (zone-edge equator eastings, scaled meridian arcs, both hemispheres), hemisphere symmetry (equal eastings, exact 10_000_000 northing offset), round trips, boundingBox/contains/distance, antimeridian no-ops in date-line zones, normalize/clamp, WKB round trips and random generation within the zone sector.

UTM zone `.prj` files could not be recognized: the fragment-set WKT matching cannot extract the zone number, and a naive fallback would mismatch - UTM `.prj` strings mention "Transverse_Mercator", which the generic EPSG:3395 fragments would match. - File-private UtmWktIdentification namespace with two Swift regex literals: a zone token (UTM / "Universal Transverse Mercator", spaces or underscores as separators, case-insensitive) capturing zone number and hemisphere, and a Central_Meridian parameter for cross-validation against the zone's expected central meridian. - The hemisphere token is mandatory; zones are validated to 1...60; an unparseable or contradicted match yields nil. - A string carrying a UTM zone token short-circuits the matcher: it is resolved to a UTM projection or rejected outright, never falling through to unrelated name fragments. - Tests: realistic ESRI-style `.prj` strings (northern and southern variants), zone/hemisphere token variants, date-line edge zones, central meridian agreement and contradiction cases, and rejected invalid tokens (missing hemisphere, zones outside 1...60). - Clears four redundant #require warnings in the MapTile tests (non-failable boundingBox getters).

Three follow-ups to the projection-genericity effort: Validation API: - Coordinate3D.isValid and BoundingBox.isValid report whether the coordinate/box lies within the valid extent of its projection, resolved through the projection registry. Always true for unbounded projections (EPSG:4978, noSRID). UTM convenience API: - UtmHemisphere is now public (it appears in the public API below) - Projection.utmZone / Projection.utmHemisphere accessors return the zone number or hemisphere of a UTM zone projection, nil for non-UTM projections - Projection(utmZone:hemisphere:) creates a zone projection from a zone number (1...60) and hemisphere Benchmarks: - ProjectionBenchmarks suite (skipped in CI like the other suites): compares the pre-registry inline formulas as a baseline against the registry path, with correctness cross-checks. An EPSG:4326 pivot short-circuit in Coordinate3D.projected(to:) skips the registry lookup and identity copy in the two most common directions (behavior-identical). - Decomposition tests isolate lookup, fixed-definition projection and full-path costs (release: ~5 ns baseline, ~30 ns lookup, ~125 ns full path per coordinate). The overhead only matters for bulk reprojection; algorithm hot paths reproject per feature, and the analysis is documented in the suite docs along with the options if a batch reprojection use case appears. Tests: ValidityTests (extent boundaries for 4326/3857/UTM sectors, unbounded projections, invalid-then-clamp) and UtmDefinitionTests.utmConvenienceApi (zone sweep, invalid zones, non-UTM accessors).

Implements issue #242 in two steps on this branch. Phase A - struct-based Projection: - Projection becomes a public struct with SRID identity: custom ==/hash and a Codable implementation that encodes the raw SRID identically to the previous enum wire format. All 125 built-in projections keep their constant names, so dot-shorthand usage (x == .epsg4326) keeps compiling unchanged. - The generated case lists are gone: kind, capability flags, wraparoundExtent and crsLength resolve through the registry. - init?(srid:) resolves through an SRID alias table (the historical EPSG:3857 aliases) plus a registry check; init?(wkt:) and the UTM token short-circuit are unchanged; descriptions unchanged. Phase B - add-only custom registration: - public CustomProjection struct: a value-type descriptor with @sendable forward/inverse closures relative to the EPSG:4326 (WGS84) pivot, plus kind, wraparound extent, valid extent (public ProjectionExtent struct), world bounding box and WKT fragment sets. The pivot contract makes datum-capable custom definitions possible. - public Projection.register(CustomProjection) is add-only: rejected for SRID <= 0, duplicates and built-in shadowing; there is no unregister/deactivate path by design. - ProjectionExtent becomes the public extent type used by all definition metadata. - The registry keeps a Mutex-guarded add-only snapshot: built-ins seeded once, custom definitions appended, never removed or replaced. Hot paths NEVER read the registry - Projection values capture their definition at construction - so projected(to:), kind, clamped(), isValid and wraparoundExtent resolve via the stored definition with zero locking. The registry is only used for lookups (init?(srid:), Codable), WKT matching and registration. Performance (release benchmarks, before -> after): - full projected(to:) path: ~178 ns/coordinate (35x baseline) -> ~1.7x of the direct-math baseline - No hot-path registry access remains; a parallel 2488-test run previously melted down (924s wall) on per-access Mutex reads, which motivated the definition capture. Tests (+6): custom registration/rejection rules, end-to-end round trips, registry-driven capabilities on a custom definition, WKT matching, Codable semantics (register-before-decode) and a concurrent registration storm. Note for test authors: Swift Testing runs tests in parallel against the process-global registry - reserve SRID blocks and use unique WKT fragments. Also adds a Projections section with a custom-projection example to the README.

Implements #243 for the WGS84-pivot architecture established in #247, shipping datum-capable built-in projections and the public transformation helper for user-defined datum CRSs. Phase 1 - Datum metadata model: - public Ellipsoid (wgs84/grs80/clarke1866/airy1830 with derived eccentricities and axes) and public Datum (wgs84/etrs89/nad27/ osgb1936) - Projection.datum accessor resolved from the definition; ProjectionDefinition gains a datum requirement (default WGS84); CustomProjection gains a datum field. Phase 2 - verified math infrastructure: - public HelmertTransformation: 7-parameter Helmert with position vector convention, exact rotation matrices, geocentric legs on arbitrary source ellipsoids, both directions with explicit result projection labeling. Well-known sets for NAD27 ("NAD27 to WGS 84 (4)", ~10 m) and OSGB36 ("OSGB 1936 to WGS 84 (6)", EPSG:1314, ~2 m) sourced from the EPSG registry and verified against pyproj/PROJ reference values. - internal TransverseMercatorMath(ellipsoid:) generalizing the Snyder TM formulas with a latitude-of-origin; UtmDefinition delegates to it (guarded by the existing UTM known-value tests, behavior provably unchanged). Phase 3 - built-in datum CRSs: - EPSG:4267 NAD27 geodetic, EPSG:4277 OSGB36 geodetic, EPSG:4258 ETRS89 geodetic (identity shift; OS labels ETRS89 data as WGS84 at meter accuracy) and EPSG:27700 British National Grid (helmert + TM on Airy 1830: lat0 49, lon0 -2, k 0.9996012717, FE 400000, FN -100000) - all registered in the registry seed with WKT matcher fragments (BNG probes before the generic Mercator fragments), datum metadata, test-appropriate unbounded extents (BNG keeps its southwest negative-northing region valid). Fixtures: The reference values were computed independently with pyproj/PROJ before porting, at 9-decimal precision; round-trip errors measured at ~3.3e-9 deg (NAD27) and ~1.4e-8 deg (OSGB36). Tests (+12): datum metadata, Helmert both directions + round trips, per-CRS reference fixtures, route-pivot round trips, British National Grid values with the southwest sector, algorithm smoke tests, WKB round trips and WKT matching for the new SRIDs. ETRS89-bearing WKT strings now match EPSG:4258 instead of returning nil. Also extends the README Projections section with the datum CRSs, the accuracy tier note and attribution, and links the grid shift support issue (#248).

Organizational restructure: all projection-related files move from Sources/GISTools/GeoJson/ into a dedicated Projections/ folder, with git-tracked renames. This is also the time to split the formerly 620-line ProjectionDefinition.swift into dedicated files: - ProjectionDefinition.swift now contains only the protocol plus the shared default-datum extension - new ProjectionRegistry.swift holds the mutex-guarded registry - each built-in projection gets its own dedicated file: Epsg4326Definition.swift, Epsg3857Definition.swift, Epsg4978Definition.swift, Epsg3395Definition.swift, Epsg32662Definition.swift and NoSridDefinition.swift, alongside the already-separate UtmDefinition and the four datum CRS definitions The README's Projections section becomes a top-level item between GeoJSON and SwiftData, with a new "Implemented projections" table (EPSG, response, coordinate units, transformation, per-projection source links), a cleaned-up custom registration example and an updated features bullet that lists the full CRS set instead of the old two-projection description. No functional changes: builds with zero warnings, 2500 tests in 183 suites green.

Adds Array<Coordinate3D>.projected(to:) and the in-place project(to:) variant: batches of coordinates sharing their projection convert in one call, with the conversion setup (TM parameterization and origin constants, Helmert rotation matrices) built once per batch instead of per coordinate. The math is identical to the single-coordinate path - the batch machinery only extracts constants into an internal hoisted representation (BatchPreparedTransforms plus prepared steps on HelmertTransformation, delegated through the same transform internals). The full frozen fixture suite proves equivalence. Semantics: all coordinates in a batch must share their projection (the library-wide invariant, asserted); noSRID values are copied verbatim in both directions; same-projection batches are verbatim copies. Migration (behavior-identical output, verified by the whole suite): - geometry projected(to:) implementations (LineString/MultiPoint as single batches; Polygon/MultiLineString/MultiPolygon as ring-level batches) - WKB decoder: coordinate lists decode raw and batch convert once per geometry (decodeCoordinate keeps the single-coordinate Point path; decodeRawCoordinate added) - WKT decoder: scans source-labeled coordinates, batch converts at the end - TWKB: coordinate sequences batch convert once; Point unchanged Also fixes a mislabeling bug found by the batch equivalence test: Etrs89Definition built its identity results with the EPSG:4326-labeled initializer, so EPSG:4258 coordinates carried an EPSG:4326 projection label (values were correct, labels were not; the per-value datum tests could not catch it). ETRS89 now carries its own SRID label in both directions. Benchmark comparison (release build, ms per 10k coordinates): - 4326 to 3857: 0.549 per-coordinate vs 0.469 batch (-15%) - 4326 to UTM to 4326: 2.717 vs 2.461 (-9%) - 4326 to British National Grid: batch 5.595 (hoists both the Helmert matrices and the TM parameterization; no visible win on that lane - two-step chain dominates) Tests (+9 in ArrayProjectionTests): edge cases, same-projection verbatim copies, batch-equals-per-coordinate equivalence across all CRS families (including the datum CRSs and ECEF), round trips, noSRID verbatim semantics, the mutating project(to:) variant. Note for test authors: the batch API asserts on mixed-projection input (library-wide invariant); uniform batches only. The three batch cases join ProjectionBenchmarks (CI-skipped).

New CRS definitions (each with an EPSG Helmert datum chain and prepared batch transforms): - EPSG:2056 CH1903+/LV95 and EPSG:21781 CH1903/LV03 via the new SwissObliqueMercatorMath, a faithful port of PROJ's somerc stage - EPSG:29902 TM65/Irish Grid and EPSG:29903 TM75/Irish Grid - EPSG:2157 IRENET95/Irish Transverse Mercator Datum model: Ellipsoid.modifiedAiry/.bessel1841 and Datum.ch1903plus/ .ch1903/.tm65/.tm75/.irenet95. All reference values pinned against pyproj ground truth before porting (somerc uses PROJ's kR = k0*a/c, not Snyder's R'). Zone selection: UtmDefinition.zone(for:) and Projection.utmZone(for:) pick the UTM zone for a coordinate including the EPSG Norway (32V) and Svalbard (31X/32X/33X/35X/37X) banding exceptions. README: table/bullet entries for the new CRSs, UTM zone-selection snippet, datum accuracy note extension, reference-link cleanup (deduplicated definitions, pruned dead refs, ascending order) at 2529 tests green

All *Definition.swift files (17 files: the EPSG bases, noSRID, UTM, NAD27, OSGB, Swiss and Irish CRSs and the ProjectionDefinition protocol) move to Sources/GISTools/Projections/Definitions; the projection and math infrastructure stays at Projections/. README definition links updated.

…ew, ETRS89/UTM belt New math (faithful ports of PROJ 9.8.1, validated against pyproj to sub-mm): LambertConformalConicMath (Snyder 15-x with the GeographicLib tauf Newton inverse), LambertAzimuthalEqualAreaMath (all four domains with the authalic series) and ObliqueStereographicMath (the Gauss-conformal reduction, including PROJ's cos^4 phi0 in the Gaussian C constant). New CRS definitions: - EPSG:25831-25837 ETRS89/UTM zone belt (TM on exact GRS80, Datum.etrs89) - EPSG:3035 ETRS89-LAEA Europe and EPSG:3034 ETRS89-LCC Europe - EPSG:2154 RGF93/Lambert-93 - EPSG:28992 Amersfoort/RD New: Helmert "Amersfoort to WGS 84 (4)" + oblique stereographic on Bessel 1841 (Datum.amersfoort) Projection.utmZone/utmHemisphere resolve both UTM belts, and WKT strings naming ETRS89 next to a zone token (ETRS89 / UTM zone 32N, ESRI ETRS_1989_UTM_Zone_32N) resolve into EPSG:25831-25837 when defined. README: table rows for the five new CRSs, datum-accuracy note extension (EPSG:4833, RD grids), reference links. Fixtures use in-zone pyproj values matching the library's stated transverse Mercator accuracy tier; the Snyder-vs-PROJ modern-TM differential at far-off-zone extents is documented and tracked in #246. 2537 tests green

Implements Karney's transverse Mercator (the Krueger series extended to 6th order, arXiv:1002.1417) for the UTM zones: conformal-sphere mapping, Clenshaw-summated Krueger series and a Newton iteration for the conformal latitude inversion. Wired into UtmDefinition (WGS84 zones) and Etrs89UtmDefinition (GRS80 belt) under the unchanged ProjectionDefinition contract; the Snyder TransverseMercatorMath stays for the other TM-based CRSs. Verification (independent Python implementation first, per the established discipline): - validated against Karney's published test set TMcoords.dat (287k points, 80-digit exact reference): worst forward 5.7 nm / reverse 5.0 nm within 3900 km of the central meridian, consistent with the paper's 5 nm bound - cross-checked against PROJ 9.8 across all 120 WGS84 zones (worst 7.5 nm over 1440 random points) and the ETRS89 belt (5.6 nm) - round-trip sweep over all zones and latitudes incl. the poles: worst 6.4 nm; Snyder degrades to ~1 km at 45 degrees off the central meridian while Karney stays sub-micrometer throughout the convergence region (~3900 km for the 5 nm bound) Behavior change (documented, sub-millimeter): Snyder and Karney agree within ~1 mm in-zone (both valid UTM); the UtmAlgorithmTests reference values were re-derived under Karney and the tolerance tightened from 1 mm to 1 um (e.g. zone 19N easting 331792.115 -> 331792.114806 m, zone 20S northing shifts 0.6 mm). Benchmarks (release, 10k coordinates, median): UTM round trip registry path 2.71 -> 7.6-8.7 ms (~2.8-3.2x, the extra ~500 ns/round trip buys nanometer accuracy and guaranteed convergence), batch path 2.47 -> 7.4-7.8 ms; unrelated paths unchanged within noise. Closes #246

trasch · 2026-09-15T18:43:22Z

Follow-up commit pushed: Replace Snyder transverse Mercator in UTM with Karney's algorithm (ce7498f), implementing #246 — details in section 12 of the PR description. Key points:

Karney 6th-order Krüger series implemented from the paper (arXiv:1002.1417), Clenshaw summation, Newton inverse; drop-in under the unchanged ProjectionDefinition contract for UtmDefinition (WGS84) and Etrs89UtmDefinition (GRS80)
Validated against Karney's 80-digit test set (worst 5.7 nm within 3900 km of the CM) and PROJ 9.8 across all 120 zones (worst 7.5 nm); round trips ≤ 6.4 nm over all zones/latitudes incl. poles
Sub-millimeter behavior change: in-zone reference values shift ≤ ~1 mm (both algorithms valid UTM); test tolerances tightened from 1 mm to 1 µm
Performance: UTM round trip ~2.8–3.2× (2.71 → 7.6–8.7 ms/10k, release); other projections unchanged

Cross-validates every sensible pair of built-in projections against reference values pre-computed with pyproj/PROJ 9.8.1: - scripts/generate_projection_matrix.py: the fixture generator. Uses explicit parameterized proj4 strings mirroring the library's parameters (including the Helmert datum transformations in position-vector form) rather than EPSG codes, so PROJ cannot silently select a better transformation (e.g. grid shifts) than the library implements. --check verifies the committed files match regeneration. - Tests/GISToolsTests/TestData/ProjectionMatrix/matrix.csv: the reference data (3612 rows, one per test point and projection; 459 points spanning global coverage, regional CRS centers/corners and all 120 UTM zones with central-meridian, zone-border and high-latitude rows). Data and test code are kept separate. - Tests/GISToolsTests/Algorithms/ProjectionMatrixTests.swift: the test logic, loading the CSV at runtime via the #filePath-relative TestData mechanism. Three suites: matrixConversions (26158 ordered pairs, x/y values plus the projection label and m passthrough), matrixRoundTrips (A -> B -> A) and batchEquivalence (the hoisted batch path on a deterministic sparse subset, ECEF sources carrying their z). Tolerances are measured ceilings, not round numbers: base rows use 1e-8 degrees (geographic) and 1 mm (planar/geocentric); audit-measured worst errors are ~1e-7 m for UTM, ~39 um for the exact-CRS pairs, 0.69-1.6 mm for 2157/3035, and the datum-pair errors up to 77 mm (4267/4277). Per-CRS budgets document the two known deviations: - The Helmert-vs-PROJ ordering difference (library applies the 3D translation in the geocentric domain, PROJ's +towgs84 offsets lat/lon directly) - few centimeters to 77 mm, root-caused and tracked as #251. EPSG:4978 pairs with the geodetic datum CRSs are excluded from the matrix until the Helmert z handling is fixed: the datum's vertical offset currently rides along as an altitude through the EPSG:4326 pivot (4326 -> datum -> 4978 differs from the direct 4326 -> 4978 by ~26 m for Amersfoort). - EPSG:3035's transcribed legacy authalic series loses ~1-2 mm at large distances from the projection origin (PROJ 9.8 uses the Karney auxlat series); 2 mm budget, upgrade tracked as #252. EPSG:32662 is excluded: the library intentionally stores degrees as x/y (see Epsg32662Definition), diverging from PROJ's meter-valued +proj=eqc; its behavior is pinned by Epsg32662Tests. 2549 tests in 190 suites, zero warnings.

trasch · 2026-09-16T07:50:48Z

Follow-up commit pushed: Add projection matrix cross-validation against PROJ reference data (8cb37323) — details in section 13 of the PR description.

Key points:

New ProjectionMatrixTests + committed generator script (scripts/generate_projection_matrix.py, --check verifies freshness): 26,158 ordered pairs across 459 points and all 129 built-in CRSs (minus documented exclusions), data kept separate in TestData/ProjectionMatrix/matrix.csv
PROJ references use explicit proj4 strings mirroring the library's parameters (incl. position-vector +towgs84), never EPSG codes — so PROJ can't silently pick grid shifts the library doesn't implement
Tolerances are measured ceilings, audited: exact-CRS pairs ≤ 39 µm, UTM ≤ 1e-7 m; per-CRS budgets only where the library knowingly deviates from PROJ, each tracked as an issue (Helmert datum transformations shift the height: 4326 -> datum CRS -> 4978 differs from the direct 4326 -> 4978 #251 Helmert z ordering / 4978-datum exclusions, LAEA (EPSG:3035) inverse loses ~1-2 mm at large distances: upgrade the authalic-latitude inversion to the Karney auxlat series #252 laea authalic-series precision)

Replaces the legacy 3-term Snyder authalic series (pj_authset/ pj_authlat) in LambertAzimuthalEqualAreaMath with PROJ 9.8's auxlat series (pj_auxlat_coeffs for AUTHALIC -> GEOGRAPHIC, order 6 in the third flattening, Karney "On auxiliary latitudes", arXiv:2212.05818), evaluated with Clenshaw summation. The series reaches full double precision for |f| <= 1/150, versus the legacy series' ~1 mm loss at large distances from the projection origin. Verification (independent Python implementation first, per the established discipline): the series was transcribed into Python and validated against PROJ 9.8 before porting — the Cairo laea inverse (the original finding of the projection matrix, ~3500 km from the EPSG:3035 origin) goes from 1.2e-8 deg (1.3 mm) to 2e-14 deg off. Measured over the matrix points: laea inverse worst 0.6 um, round trip worst 0.5 um. Closes #252. The EPSG:3035 pair budget is removed from the projection matrix (all 3035-sourced pairs now held to the plain 1 mm row tolerance, measured headroom ~1600x); the fixture is regenerated. A new laeaAuxlatSeriesRoundTrips test pins the series across the full EPSG:3035 domain (origin, the failing Cairo distance, the Canary worst case, far corners, equator, southern hemisphere) at 1 um forward / 1e-9 deg round-trip tolerance. The LambertConformalConicMath inverse was checked per the issue: it uses the conformal-latitude tauf Newton iteration, not the authalic series — no change needed. 2550 tests in 190 suites, zero warnings.

trasch · 2026-09-16T08:11:06Z

Follow-up commit pushed: Upgrade the laea authalic latitude to the Karney auxlat series (d6fdc0e), closing #252 — details in section 14 of the PR description.

Key points:

LambertAzimuthalEqualAreaMath inverse now uses the Karney auxlat series (order 6 in n, Clenshaw summation) instead of the legacy 3-term Snyder authalic series: Cairo's laea inverse error goes from 1.3 mm to 2e-14 deg; measured worst over the matrix points is 0.6 µm (inverse) / 0.5 µm (round trip)
Validated in Python against PROJ 9.8 before porting, per the established discipline
The EPSG:3035 pair budget is removed from the projection matrix (fixture regenerated); new laeaAuxlatSeriesRoundTrips test pins the full EPSG:3035 domain at 1 µm / 1e-9 deg
LambertConformalConicMath checked per the issue: conformal-latitude Newton iteration, no shared pattern, no change needed

The Helmert transformation now passes the input height through unchanged in both directions (the 3D transform still runs on the ECEF position, so x/y remain datum-correct at the given height; only the output z is the input height instead of the datum-shifted one). This matches PROJ's +towgs84 nominal-height semantics, validated in Python against PROJ before porting (datum -> WGS84 with h=100 reproduces the original lat/lon/height to ~0.3 mm). Previously the datum's vertical offset rode along as an altitude through the EPSG:4326 pivot: the library was internally asymmetric (the forward path passed the input height through, the inverse applied the full 3D shift), and 4326 -> datum CRS -> 4978 differed from the direct 4326 -> 4978 by the datum's vertical offset (~26 m for Amersfoort on Bessel 1841, issue #251's finding from the projection matrix). Consequences: - All 176 EPSG:4978 <-> geodetic-datum-CRS pairs are now verified in the projection matrix (26334 pairs, up from 26158); the exclusion machinery is removed from the generator and the emitted test. Measured worst errors: 1.9 mm (4978 -> Irish, within its budget) and 48.9 mm (datum -> 4978 via NAD27, within the documented 0.1 m budget for the geocentric-domain Helmert ordering difference). - The Helmert round trip is no longer exactly closed (~1e-8 deg, the same looseness PROJ's own +towgs84 round trips exhibit); HelmertTransformationTests reflects that with an explanatory comment. The z round trip itself stays exact. - The type and method documentation spell out the horizontal semantics; the generator's datum-budget comments note the measured residuals. Closes #251. 2550 tests in 190 suites, zero warnings.

trasch · 2026-09-16T08:33:32Z

Follow-up commit pushed: Preserve the height in the Helmert datum transformations (b482f28), closing #251 — details in section 15 of the PR description.

Key points:

HelmertTransformation passes the input height through in both directions (3D transform still runs on ECEF; output z = input height), matching PROJ's +towgs84 nominal-height semantics — validated in Python against PROJ before porting
Fixes the pivot inconsistency: 4326 → datum CRS → 4978 now equals the direct 4326 → 4978 (was ~26 m off for Amersfoort)
All 176 EPSG:4978 ↔ datum-CRS pairs verified in the projection matrix (26,334 pairs; measured worst 48.9 mm, within the documented 0.1 m Helmert-ordering budget); exclusion machinery removed
Helmert round trip loosens to ~1e-8 deg — identical to PROJ's own +towgs84 round-trip looseness; z round trip stays exact

Both matrix findings (#251, #252) are now fixed; the matrix's remaining datum budgets document the residual Helmert ordering difference only.

New built-in CRSs on the Bessel 1841 ellipsoid with the classical transverse Mercator math (Karney series, consistent with the UTM zones): - EPSG:31466-31469, DHDN / Gauss-Krueger zones 2-5 (the German cadastral grid): k=1, x0 = zone*1M + 500k (the zone-prefix convention), Datum.dhdn with the EPSG "DHDN to WGS 84 (2)" Helmert (598.1, 73.7, 418.2, 0.202", 0.045", -2.455", 6.7 ppm, stated accuracy ~3 m) - EPSG:31255-31259, MGI / Austria Gauss-Krueger (all five zones: the x0=0 Central/East variants plus the M28/M31/M34 three-zone belt with the Ferro-based numbering): y0 = -5M, Datum.mgi with the EPSG "MGI to WGS 84 (2)" Helmert (577.326, 90.129, 463.919, 5.137", 1.474", 5.297", 2.4232 ppm, ~1.5 m) Verification per the established discipline: pyproj/PROJ 9.8.1 fixtures computed before finalizing (34 forward values across all 9 zones plus round trips). One implementation bug caught en route (the first pass omitted the Helmert, 52-89 m off; fixed by wiring the datum chain like the other datum CRSs). Documented 1.9 mm residual on MGI: PROJ applies datum shifts as a geodetic-domain approximation while the library uses the exact XYZ affine - invisible for DHDN's 0.2" rotations, measurable for MGI's 5" ones; covered by the 2 mm datum budget with a comment. Integration: - New GaussKruegerWktIdentification (like the UTM short-circuit): handles GDAL ("DHDN / Gauss-Krueger zone 3") and ESRI ("Gauss_Kruger_DHDN_3") names, the Ferro M-tokens (M28/M31/M34), and parameter-based zone resolution with Central_Meridian / False_Easting cross-validation (the x0=0 variants share CMs with the modern zones, so the false easting disambiguates) - Projection matrix: all 9 GK CRSs added (3626 rows, 26700 pairs) - New GaussKruegerTests (9 tests): SRID mapping, registry coverage, both belts' fixtures, zone-prefix convention, false-northing convention, datum metadata, WKT matching, batch equivalence - README: feature list and two table rows + reference links

trasch · 2026-09-16T10:28:43Z

Follow-up commit pushed: Add German and Austrian Gauss-Krueger zone belts (1f31f9d) — details in section 16 of the PR description.

Key points:

EPSG:31466–31469 (DHDN / Gauss-Krüger zones 2–5) and EPSG:31255–31259 (MGI / Austria GK, all five zones incl. the x₀=0 Central/East variants and M28/M31/M34) on Bessel 1841, Karney transverse Mercator, with new Datum.dhdn / Datum.mgi Helmert chains
PROJ 9.8.1 fixtures computed before finalizing; one implementation bug caught en route (omitted Helmert, 52–89 m)
Documented 1.9 mm residual on MGI (PROJ's geodetic-domain approximation vs the library's exact XYZ affine; invisible for DHDN's 0.2″ rotations) — covered by a 2 mm datum budget
New GaussKruegerWktIdentification (UTM-style short-circuit): GDAL/ESRI names, Ferro M-tokens, parameter-based zone resolution with CM/False_Easting cross-validation
Projection matrix extended to all 9 GK CRSs (3,626 rows, 26,700 pairs); 9 new GaussKruegerTests

New Datum.nad83 (identity over WGS84, like ETRS89; the NAD83/WGS84 difference is continental plate motion with no simple Helmert) on GRS80, with six new CRS definitions: - EPSG:4269 NAD83 geodetic (identity) - EPSG:26901-26960 NAD83/UTM zones 1N-60N (Karney TM on GRS80, mirroring the ETRS89 belt) - EPSG:5070 NAD83 / Conus Albers (the CONUS equal-area standard, USGS/EPA/NLCD) and EPSG:3005 NAD83 / BC Albers via a new AlbersEqualAreaMath (Snyder 14-1..14-14) - EPSG:3347 NAD83 / Statistics Canada Lambert and EPSG:3978 NAD83 / Canada Atlas Lambert via the existing LambertConformalConicMath The authalic-latitude conversion is extracted from LambertAzimuthalEqualAreaMath into a shared internal helper (AuthalicLatitude) used by both laea and albers: the Albers inverse uses the Karney auxlat series, matching PROJ 9.8's inverse exactly and avoiding the Newton-stall failure mode near the poles that a naive q-iteration exhibits. Verification per the established discipline: an independent Python Albers implementation was validated against PROJ before porting - two transcription bugs caught en route (a ratio-vs-log mixup in Snyder's q, and missing longitude normalization past +/-180). Final sweep over 3000 points on EPSG:5070/3005: forward worst 33 nm vs PROJ, round trips 3.7 um. Swift-vs-PROJ spot check to 2e-12 m. Integration: - WKT identification: NAD83 datum tokens (NAD83, NAD_1983, ESRI GCS_North_American_1983) resolve to EPSG:4269; the UTM zone short-circuit routes NAD83 zone tokens to the EPSG:269xx belt - Projection matrix: all new CRSs included (4409 rows, 41290 pairs) - New Nad83CrsTests (6 tests): SRID mapping, fixtures across all CRSs, Albers round trips, datum metadata, WKT matching, batch equivalence - README: feature list, six table rows, reference links 2565 tests in 192 suites, zero warnings.

trasch · 2026-09-16T10:48:47Z

Follow-up commit pushed: Add North American NAD83 CRSs: UTM belt, Albers and Canadian grids (06a58ae) — details in section 17 of the PR description.

Key points:

Datum.nad83 (identity over WGS84, like ETRS89, GRS80) + 6 CRS families: EPSG:4269 (geodetic), EPSG:26901–26960 (NAD83/UTM belt), EPSG:5070 (Conus Albers) / EPSG:3005 (BC Albers) via the new AlbersEqualAreaMath, EPSG:3347/3978 (Canadian Lambert) via the existing LCC math
AuthalicLatitude extracted as a shared helper (laea + albers): the Albers inverse uses the Karney auxlat series, matching PROJ 9.8 exactly and avoiding the near-pole Newton stall
Python-first verification: two transcription bugs caught before porting (ratio-vs-log in Snyder's q, missing longitude normalization past ±180°); measured worst 33 nm forward / 3.7 µm round trips over a 3,000-point sweep
WKT: NAD83 datum tokens → EPSG:4269; NAD83 UTM zone tokens → the 269xx belt
Projection matrix extended to 4,409 rows / 41,290 pairs; 6 new Nad83CrsTests

trasch added 11 commits September 14, 2026 13:53

NoSRID coordinates are copied verbatim

5e639b0

Updated AGENTS.md

5f0cd25

trasch self-assigned this Sep 15, 2026

trasch added the enhancement New feature or request label Sep 15, 2026

trasch added 6 commits September 15, 2026 11:09

trasch added 4 commits September 16, 2026 11:18

Some cleanup

024ae64

Moved some test files around

a2d77af

Moved another two files

65f2278

trasch merged commit d2e00dc into main Sep 16, 2026
1 check passed

trasch deleted the generic-projections branch September 16, 2026 10:57

Navigation Menu

Repository navigation

Make the projection system generic (phases 1–5 + follow-ups) - #250

Make the projection system generic (phases 1–5 + follow-ups)#250
trasch merged 25 commits into
mainfrom
generic-projections

trasch commented Sep 15, 2026 •

edited

Loading

Uh oh!

trasch commented Sep 15, 2026

Uh oh!

trasch commented Sep 16, 2026

Uh oh!

trasch commented Sep 16, 2026

Uh oh!

trasch commented Sep 16, 2026

Uh oh!

trasch commented Sep 16, 2026

Uh oh!

trasch commented Sep 16, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Conversation

trasch commented Sep 15, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

What it delivers

1. ProjectionKind + capability flags (#237)

2. ProjectionDefinition registry (#238)

3. New projections

4. Registry-driven WKT matching + coder/tile tests (#240)

5. Struct-based projection + add-only custom registration (#247)

6. Batch coordinate conversion (commit 376a561)

7. Validation API + UTM convenience (#245)

8. Swiss and Irish CRSs (commits 6c9b6f7)

9. UTM zone selection (commit 6c9b6f7)

10. Pan-European and national CRSs (commit 4d2a332)

11. Karney transverse Mercator for UTM (commit ce7498f)

12. Projection matrix cross-validation (commit 8cb37323)

13. LAEA authalic latitude upgraded to the auxlat series (commit d6fdc0e)

14. Helmert height passthrough (commit b482f28)

15. German and Austrian Gauss-Krüger zone belts (commit 1f31f9d)

17. North American NAD83 CRSs (commit 06a58ae)

After merge

Uh oh!

trasch commented Sep 15, 2026

Uh oh!

trasch commented Sep 16, 2026

Uh oh!

trasch commented Sep 16, 2026

Uh oh!

trasch commented Sep 16, 2026

Uh oh!

trasch commented Sep 16, 2026

Uh oh!

trasch commented Sep 16, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

trasch commented Sep 15, 2026 •

edited

Loading

6. Batch coordinate conversion (commit `376a561`)

8. Swiss and Irish CRSs (commits `6c9b6f7`)

9. UTM zone selection (commit `6c9b6f7`)

10. Pan-European and national CRSs (commit `4d2a332`)

11. Karney transverse Mercator for UTM (commit `ce7498f`)

12. Projection matrix cross-validation (commit `8cb37323`)

13. LAEA authalic latitude upgraded to the auxlat series (commit `d6fdc0e`)

14. Helmert height passthrough (commit `b482f28`)

15. German and Austrian Gauss-Krüger zone belts (commit `1f31f9d`)

17. North American NAD83 CRSs (commit `06a58ae`)