Ship datum support: ETRS89, NAD27, OSGB 1936 and British National Grid (#243)

AGENTS.md

@@ -173,7 +173,14 @@ When a function **call** is split across multiple lines, place each argument on
 ## General instructions
 
 - DO NOT take any shortcuts while implementing an algorithm. Correctness is the highest priority
-- DO NOT commit changes unless the user tells you to do so, ALWAYS let the user review your changes
+
+## NEVER COMMIT OR PUSH WITHOUT THE USER'S EXPLICIT REVIEW — THIS IS THE HIGHEST PRIORITY RULE IN THIS FILE
+
+- **NEVER commit without the user having explicitly reviewed the changes.** Not once, not "just this once", not even when the user asks for a PR implicitly or you believe they'd want it. Zero exceptions without a direct instruction.
+- **NEVER push** to any remote without the user's explicit, per-action consent — this includes branch moves, force-pushes, and pushes done as part of "creating a PR".
+- **ALWAYS leave changes uncommitted** in the working tree for review. Let the user run `git add`/`git commit`/`git push` themselves, or explicitly grant permission.
+- **When in ANY doubt, ask first.** A wrongly published change cannot be unpublished; a delayed commit costs nothing. This has caused real problems before (public review-bypass, lost user control); treat every `git commit`, `git push`, and `gc` remote action as requiring an unambiguous "yes, commit/push" from the user.
+- An explicit user request to commit applies **only to that request** — never carry permission forward into follow-up turns; each new task starts back at "working tree only".
 - DO NOT create free functions (un-namespaced top-level functions). Always use a `private enum` namespace or extensions on existing types.
 - Code MUST compile cleanly, with no warnings
 - New algorithms and bug fixes MUST include tests for all projections (EPSG:4326, EPSG:3857, EPSG:4978, noSRID)

README.md

@@ -451,7 +451,11 @@ let projectedBack = mercator.projected(to: .epsg4326)
 print(projection.description, projection.srid)
 ```
 
-Densification, buffer, distance etc. automatically take the projection into account. Beyond the built-in projections (EPSG:4326, 3857, 4978, 3395, 32662 and all 120 UTM zones), custom projections can be registered - registration is add-only, applied for the whole process, typically at startup:
+Densification, buffer, distance etc. automatically take the projection into account. Beyond the built-in projections (EPSG:4326, 3857, 4978, 3395, 32662, 4258/ETRS89, 4267/NAD27, 4277/OSGB 1936, 27700/British National Grid and all 120 UTM zones), custom projections can be registered - registration is add-only, applied for the whole process, typically at startup:
+```swift
+let coordinate = Coordinate3D(latitude: 41.0, longitude: -71.0).projected(to: .epsg4267)  // NAD27
+let datum = Projection.epsg27700.datum  // Datum.osgb1936, Airy 1830
+```
 
 ```swift
 let custom = CustomProjection(
@@ -473,6 +477,11 @@ let customProjection = try Projection(srid: 900_001)
 
 All algorithms dispatch on the projection *kind* (geographic/planar/geocentric) and honor the definition's capabilities (wraparound extents, valid ranges, world bounding boxes), so custom projections work across the whole library like built-in ones.
 
+### Datum accuracy note
+The datum-capable built-in projections (NAD27, OSGB 1936, British National Grid) use the published EPSG Helmert transformations ("NAD27 to WGS 84 (4)", ~10 m; "OSGB 1936 to WGS 84 (6)" / EPSG:1314, ~2 m). Sub-meter-centimeter accuracy for NAD27 (NADCON) or Great Britain (OSTN15, EPSG:7709) requires grid shift files, which the library deliberately does not bundle (see [#248]). Datum transformations are parameterized in HELMERT via `HelmertTransformation` for use in your own `CustomProjection` definitions.
+
+Attribution: EPSG parameter values are based on the EPSG dataset (https://epsg.org) used under its terms; OS transform parameters reference the Ordnance Survey *Guide to Coordinate Systems in Great Britain*.
+
 ## BoundingBox
 [Implementation][18] / [BoundingBox test cases][19]
 
@@ -1702,3 +1711,5 @@ Thomas Rasch, Outdooractive
 
 [image-1]:	https://img.shields.io/endpoint?url=https%3A%2F%2Fswiftpackageindex.com%2Fapi%2Fpackages%2FOutdooractive%2Fgis-tools%2Fbadge%3Ftype%3Dswift-versions
 [image-2]:	https://img.shields.io/endpoint?url=https%3A%2F%2Fswiftpackageindex.com%2Fapi%2Fpackages%2FOutdooractive%2Fgis-tools%2Fbadge%3Ftype%3Dplatforms
+
+[#248]:	https://github.com/Outdooractive/gis-tools/issues/248

Sources/GISTools/Extensions/DoubleExtensions.swift

@@ -18,6 +18,13 @@ extension Double {
         self.remainder(dividingBy: 2.0 * .pi) * 180.0 / .pi
     }
 
+    /// Converts an angle in arc seconds to radians.
+    ///
+    /// - Returns: The angle in radians.
+    public var arcSecondsToRadians: Self {
+        self / 3_600.0 * .pi / 180.0
+    }
+
 }
 
 extension Double {

Sources/GISTools/GeoJson/CustomProjection.swift

@@ -50,6 +50,11 @@ public struct CustomProjection: Sendable {
     /// algorithm dispatch.
     public var kind: ProjectionKind
 
+    /// The geodetic datum of the projection's coordinate frame. For
+    /// non-WGS84 datums, ``forward``/``inverse`` perform the datum
+    /// transformation relative to the pivot. Defaults to ``Datum/wgs84``.
+    public var datum: Datum
+
     /// The absolute value beyond which the horizontal axis wraps around
     /// (e.g. the antimeridian for geographic coordinate systems).
     /// `nil` disables wrapping.
@@ -82,6 +87,7 @@ public struct CustomProjection: Sendable {
     public init(
         srid: Int,
         kind: ProjectionKind,
+        datum: Datum = .wgs84,
         wraparoundExtent: Double? = nil,
         validExtent: ProjectionExtent? = nil,
         worldBoundingBox: BoundingBox? = nil,
@@ -91,6 +97,7 @@ public struct CustomProjection: Sendable {
     ) {
         self.srid = srid
         self.kind = kind
+        self.datum = datum
         self.wraparoundExtent = wraparoundExtent
         self.validExtent = validExtent
         self.worldBoundingBox = worldBoundingBox
@@ -137,6 +144,7 @@ struct RegisteredCustomDefinition: ProjectionDefinition {
 
     var projection: Projection { Projection(uncheckedSrid: custom.srid, definition: self) }
     var kind: ProjectionKind { custom.kind }
+    var datum: Datum { custom.datum }
     var wraparoundExtent: Double? { custom.wraparoundExtent }
     var validExtent: ProjectionExtent? { custom.validExtent }
     var worldBoundingBox: BoundingBox? { custom.worldBoundingBox }

Sources/GISTools/GeoJson/Datum.swift

@@ -0,0 +1,125 @@
+
+/// The ellipsoid of a geodetic datum.
+///
+/// Reference values follow the EPSG registry (see EPSG guidance note 7-2).
+public struct Ellipsoid:
+    Hashable,
+    Sendable,
+    CustomStringConvertible
+{
+
+    /// Semi-major axis, in meters.
+    public let semiMajorAxis: Double
+
+    /// The inverse of the flattening (a / (a - b)).
+    public let inverseFlattening: Double
+
+    public init(
+        semiMajorAxis: Double,
+        inverseFlattening: Double
+    ) {
+        self.semiMajorAxis = semiMajorAxis
+        self.inverseFlattening = inverseFlattening
+    }
+
+    // MARK: - Well-known ellipsoids
+
+    /// WGS84 (EPSG:7030).
+    public static let wgs84 = Ellipsoid(
+        semiMajorAxis: 6_378_137.0,
+        inverseFlattening: 298.257223563)
+
+    /// GRS80 (EPSG:7019), the ETRS89 ellipsoid.
+    public static let grs80 = Ellipsoid(
+        semiMajorAxis: 6_378_137.0,
+        inverseFlattening: 298.257222101)
+
+    /// Clarke 1866 (EPSG:7008), the NAD27 ellipsoid.
+    public static let clarke1866 = Ellipsoid(
+        semiMajorAxis: 6_378_206.4,
+        inverseFlattening: 294.978698213898)
+
+    /// Airy 1830 (EPSG:7001), the OSGB 1936 ellipsoid.
+    public static let airy1830 = Ellipsoid(
+        semiMajorAxis: 6_377_563.396,
+        inverseFlattening: 299.3249646)
+
+    /// A textual description of the receiver.
+    public var description: String {
+        "Ellipsoid(a: \(semiMajorAxis), 1/f: \(inverseFlattening))"
+    }
+
+    // MARK: - Derived values
+
+    /// First eccentricity squared, e².
+    public var eccentricitySquared: Double {
+        let f = 1.0 / inverseFlattening
+        return 2.0 * f - f * f
+    }
+
+    /// Semi-minor axis, in meters.
+    public var semiMinorAxis: Double {
+        semiMajorAxis * (1.0 - 1.0 / inverseFlattening)
+    }
+
+}
+
+/// A geodetic datum, combining an ellipsoid with its realization.
+public struct Datum:
+    Hashable,
+    Sendable,
+    CustomStringConvertible
+{
+
+    /// The foundation all library projections are based on (GPS etc.).
+    public static let wgs84 = Datum(
+        name: "WGS 84",
+        realization: "World Geodetic System 1984",
+        ellipsoid: .wgs84)
+
+    /// European Terrestrial Reference System 1989. Computed positions
+    /// effectively coincide with WGS84 at meter accuracy (centimeter-level
+    /// drift is not represented).
+    public static let etrs89 = Datum(
+        name: "ETRS89",
+        realization: "European Terrestrial Reference System 1989",
+        ellipsoid: .grs80)
+
+    /// The classical North American datum (Clarke 1866), defined through
+    /// triangulation.
+    public static let nad27 = Datum(
+        name: "NAD 1927",
+        realization: "North American Datum of 1927",
+        ellipsoid: .clarke1866)
+
+    /// Ordnance Survey of Great Britain 1936 (Airy 1830).
+    public static let osgb1936 = Datum(
+        name: "OSGB 1936",
+        realization: "Ordnance Survey Great Britain 1936 triangulation",
+        ellipsoid: .airy1830)
+
+    /// A human readable datum name.
+    public let name: String
+
+    /// The realization the datum is based on.
+    public let realization: String
+
+    /// The ellipsoid of the datum.
+    public let ellipsoid: Ellipsoid
+
+    public init(
+        name: String,
+        realization: String,
+        ellipsoid: Ellipsoid
+    ) {
+        self.name = name
+        self.realization = realization
+        self.ellipsoid = ellipsoid
+    }
+
+    /// A textual description of the receiver.
+    public var description: String {
+        "\(name) (\(realization), \(ellipsoid))"
+    }
+
+}

Sources/GISTools/GeoJson/Etrs89Definition.swift

@@ -0,0 +1,54 @@
+#if canImport(CoreLocation)
+import CoreLocation
+#endif
+import Foundation
+
+/// ETRS89 math for EPSG:4258 (ETRS89 geodetic coordinates on the GRS80
+/// ellipsoid).
+///
+/// ETRS89 and WGS84 effectively coincide at meter accuracy: the plate-fixed
+/// ETRS89 frame differs from WGS84 by centimeters that grow over time, far
+/// below the tolerance of this library's algorithms. The transformation is
+/// therefore the identity (same behavior as [Ordnance Survey's
+/// approach](https://docs.os.uk/osngd/getting-started/os-ngd-fundamentals/coordinate-reference-systems)
+/// of providing ETRS89 data labelled as WGS84).
+///
+/// For cm-level accuracy an epoch-dependent transformation would be needed,
+/// which is out of scope for this library.
+struct Etrs89Definition: ProjectionDefinition {
+
+    var projection: Projection { Projection(uncheckedSrid: 4258, definition: self) }
+    var kind: ProjectionKind { .geographic }
+    var datum: Datum { .etrs89 }
+
+    var wraparoundExtent: Double? { 180.0 }
+
+    var validExtent: ProjectionExtent? {
+        ProjectionExtent(minX: -180.0, minY: -90.0, maxX: 180.0, maxY: 90.0)
+    }
+
+    var worldBoundingBox: BoundingBox? {
+        BoundingBox.world
+    }
+
+    var wktMatchers: [[String]] {
+        [
+            ["GEOGCS", "ETRS89"],
+            ["GEOGCS", "ETRS_1989"],
+            ["GEOGCS", "ETRS 1989"],
+        ]
+    }
+
+    func forward(_ coordinate: Coordinate3D) -> Coordinate3D {
+        Coordinate3D(
+            latitude: coordinate.latitude,
+            longitude: coordinate.longitude,
+            altitude: coordinate.altitude,
+            m: coordinate.m)
+    }
+
+    func inverse(_ coordinate: Coordinate3D) -> Coordinate3D {
+        forward(coordinate)
+    }
+
+}