From b84315fff16ba4c2061d923aabb917bf49b4cd69 Mon Sep 17 00:00:00 2001 From: "Thomas Rasch (OpenCode)" Date: Tue, 15 Sep 2026 10:04:37 +0000 Subject: [PATCH 1/2] Support non-WGS84 datums: ETRS89/NAD27/OSGB1936/BNG built-in CRSs 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). --- README.md | 13 +- .../Extensions/DoubleExtensions.swift | 7 + .../GISTools/GeoJson/CustomProjection.swift | 8 + Sources/GISTools/GeoJson/Datum.swift | 125 +++++++ .../GISTools/GeoJson/Etrs89Definition.swift | 54 +++ .../GeoJson/HelmertTransformation.swift | 321 ++++++++++++++++++ .../GISTools/GeoJson/Nad27Definition.swift | 51 +++ .../GeoJson/Osgb1936BngDefinition.swift | 95 ++++++ .../GISTools/GeoJson/Osgb1936Definition.swift | 48 +++ Sources/GISTools/GeoJson/Projection.swift | 15 + .../GeoJson/ProjectionDefinition.swift | 26 ++ .../GeoJson/TransverseMercatorMath.swift | 164 +++++++++ Sources/GISTools/GeoJson/UtmDefinition.swift | 98 ++---- .../Algorithms/DatumCrsTests.swift | 165 +++++++++ Tests/GISToolsTests/GeoJson/DatumTests.swift | 52 +++ .../GeoJson/HelmertTransformationTests.swift | 65 ++++ .../GeoJson/ProjectionWktMatcherTests.swift | 3 +- 17 files changed, 1235 insertions(+), 75 deletions(-) create mode 100644 Sources/GISTools/GeoJson/Datum.swift create mode 100644 Sources/GISTools/GeoJson/Etrs89Definition.swift create mode 100644 Sources/GISTools/GeoJson/HelmertTransformation.swift create mode 100644 Sources/GISTools/GeoJson/Nad27Definition.swift create mode 100644 Sources/GISTools/GeoJson/Osgb1936BngDefinition.swift create mode 100644 Sources/GISTools/GeoJson/Osgb1936Definition.swift create mode 100644 Sources/GISTools/GeoJson/TransverseMercatorMath.swift create mode 100644 Tests/GISToolsTests/Algorithms/DatumCrsTests.swift create mode 100644 Tests/GISToolsTests/GeoJson/DatumTests.swift create mode 100644 Tests/GISToolsTests/GeoJson/HelmertTransformationTests.swift diff --git a/README.md b/README.md index d6929783..dc5d77f1 100644 --- a/README.md +++ b/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 diff --git a/Sources/GISTools/Extensions/DoubleExtensions.swift b/Sources/GISTools/Extensions/DoubleExtensions.swift index 5e24a9b1..895574aa 100644 --- a/Sources/GISTools/Extensions/DoubleExtensions.swift +++ b/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 { diff --git a/Sources/GISTools/GeoJson/CustomProjection.swift b/Sources/GISTools/GeoJson/CustomProjection.swift index e17c2b0f..2ae06573 100644 --- a/Sources/GISTools/GeoJson/CustomProjection.swift +++ b/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 } diff --git a/Sources/GISTools/GeoJson/Datum.swift b/Sources/GISTools/GeoJson/Datum.swift new file mode 100644 index 00000000..73a2d0e7 --- /dev/null +++ b/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))" + } + +} diff --git a/Sources/GISTools/GeoJson/Etrs89Definition.swift b/Sources/GISTools/GeoJson/Etrs89Definition.swift new file mode 100644 index 00000000..d6d42eb1 --- /dev/null +++ b/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) + } + +} diff --git a/Sources/GISTools/GeoJson/HelmertTransformation.swift b/Sources/GISTools/GeoJson/HelmertTransformation.swift new file mode 100644 index 00000000..6b4cc0f9 --- /dev/null +++ b/Sources/GISTools/GeoJson/HelmertTransformation.swift @@ -0,0 +1,321 @@ +#if canImport(CoreLocation) +import CoreLocation +#endif +import Foundation + +/// A Helmert datum transformation (translation + rotation + scale). +/// +/// Coordinates are transformed between a geodetic datum and WGS84. The +/// parameters follow the **position vector convention** and are listed in +/// the direction *datum → WGS84* (the direction the EPSG registry uses). +/// Both directions are supported: +/// ``transform(datumToWgs84:)`` and ``transform(wgs84ToDatum:)``. +/// +/// The transformation operates on geocentric (ECEF) coordinates, so the +/// datum's ``Datum/ellipsoid`` must be correct: the geocentric legs use it. +/// +/// - Note: Well-known transformations for NAD27 and OSGB 1936 are provided +/// (``nad27``/``osgb1936``), sourced from the EPSG registry. Sub-meter +/// datum accuracy requires grid shift files (NADCON/OSTN15), which the +/// library does not incorporate (see the grid shift support issue). +/// +/// ## Example +/// +/// ```swift +/// let helmert = HelmertTransformation.osgb1936 +/// let coordinate = Coordinate3D(latitude: 51.4778, longitude: -0.0015, altitude: 46.0) +/// let osgb = helmert.transform(wgs84ToDatum: coordinate) // OSGB36 lat/lon +/// ``` +public struct HelmertTransformation: + Sendable, + CustomStringConvertible +{ + + // MARK: - Well-known transformations + + /// **NAD27 to WGS 84 (4)** ("NAD27 to WGS 84 (4)" in the EPSG registry, + /// stated accuracy approximately 10 m): a translation-only + /// transformation. + /// + /// Sub-meter or centimeter-level accuracy for NAD27 requires NADCON + /// grid shifts, which the library does not provide (see the grid shift + /// support issue). + public static let nad27 = HelmertTransformation( + datum: .nad27, + dx: -8.0, + dy: 160.0, + dz: 176.0) + + /// **OSGB 1936 to WGS 84 (6)** (EPSG:1314, the OS-documented + /// approximate Helmert transformation, stated accuracy approximately + /// 2 m; the OS-recommended route for sub-meter accuracy is the OSTN15 + /// grid, which the library does not incorporate: see the grid shift + /// support issue). + public static let osgb1936 = HelmertTransformation( + datum: .osgb1936, + dx: 446.448, + dy: -125.157, + dz: 542.06, + rx: 0.15, + ry: 0.247, + rz: 0.842, + scalePpm: -20.489) + + // MARK: - Parameters + + /// The datum frame the transformation maps *from* (its geodetic + /// coordinates feed the geocentric leg). + public let datum: Datum + + /// Translation along the x/y/z axes, in meters (datum → WGS84). + public let dx: Double + /// Translation along the y axis, in meters (datum → WGS84). + public let dy: Double + /// Translation along the z axis, in meters (datum → WGS84). + public let dz: Double + + /// Rotation about the x axis, in arc seconds (datum → WGS84, + /// position vector convention). + public let rx: Double + /// Rotation about the y axis, in arc seconds. + public let ry: Double + /// Rotation about the z axis, in arc seconds. + public let rz: Double + + /// Scale difference, in parts per million (datum → WGS84). + public let scalePpm: Double + + /// A textual description of the receiver. + public var description: String { + "HelmertTransformation <- \(datum.name) (dx: \(dx), dy: \(dy), dz: \(dz), rx: \(rx), ry: \(ry), rz: \(rz), s: \(scalePpm))" + } + + // MARK: - Initialization + + /// Creates a Helmert transformation *from* a datum frame *towards* + /// WGS84 (position vector convention). + /// + /// - Parameters: + /// - datum: The datum the transformation maps *from* + /// - dx: Translation along the x axis, in meters + /// - dy: Translation along the y axis, in meters + /// - dz: Translation along the z axis, in meters + /// - rx: Rotation about the x axis, in arc seconds + /// - ry: Rotation about the y axis, in arc seconds + /// - rz: Rotation about the z axis, in arc seconds + /// - scalePpm: Scale difference, in parts per million + public init( + datum: Datum, + dx: Double, + dy: Double, + dz: Double, + rx: Double = 0.0, + ry: Double = 0.0, + rz: Double = 0.0, + scalePpm: Double = 0.0 + ) { + self.dx = dx + self.dy = dy + self.dz = dz + self.rx = rx + self.ry = ry + self.rz = rz + self.scalePpm = scalePpm + self.datum = datum + } + + // MARK: - Transformation + + /// Transforms WGS84 coordinates into the datum's geographic coordinate + /// system. + /// + /// - Parameters: + /// - coordinate: A coordinate in EPSG:4326 + /// - projection: The projection label carried by the result. + /// Geographic CRSs of the datum pass their own SRID here to keep + /// the label correct (e.g. EPSG:4267). Defaults to EPSG:4326, + /// marking the values as "geodetic but without a datum-specific + /// SRID". + /// - Returns: A coordinate in the datum's frame (`latitude`/ + /// `longitude` in degrees) + public func transform( + wgs84ToDatum coordinate: Coordinate3D, + projection: Projection? = nil + ) -> Coordinate3D { + Self.transform( + coordinate, + sourceEllipsoid: Datum.wgs84.ellipsoid, + targetEllipsoid: datum.ellipsoid, + rotationIsForward: false, + helmert: self, + targetProjection: projection ?? .epsg4326) + } + + /// Transforms datum-frame coordinates into WGS84. + /// + /// - Parameter coordinate: A coordinate in the datum's frame + /// - Returns: A coordinate in EPSG:4326 + public func transform(datumToWgs84 coordinate: Coordinate3D) -> Coordinate3D { + Self.transform( + coordinate, + sourceEllipsoid: datum.ellipsoid, + targetEllipsoid: Datum.wgs84.ellipsoid, + rotationIsForward: true, + helmert: self, + targetProjection: .epsg4326) + } + + // MARK: - Shared implementation + + /// Shared implementation: `coordinate` is in the frame of + /// `sourceEllipsoid`; `rotationIsForward` selects the datum → WGS84 + /// direction of the stored parameters. + private static func transform( + _ coordinate: Coordinate3D, + sourceEllipsoid: Ellipsoid, + targetEllipsoid: Ellipsoid, + rotationIsForward: Bool, + helmert: HelmertTransformation, + targetProjection: Projection + ) -> Coordinate3D { + let phi = coordinate.latitude.degreesToRadians + let lambda = coordinate.longitude.degreesToRadians + let h = coordinate.altitude ?? 0.0 + + // Geodetic -> geocentric on the source ellipsoid. + let (x1, y1, z1) = geodeticToEcef( + latitude: phi, + longitude: lambda, + height: h, + ellipsoid: sourceEllipsoid) + + // Position vector convention: E = d + s * R * X (datum → WGS84), + // with R the rotation matrix built from (rx, ry, rz). + // The WGS84 → datum direction applies R^T and inverted shift/scale. + var ecef: (x: Double, y: Double, z: Double) + if rotationIsForward { + let scale = 1.0 + helmert.scalePpm / 1_000_000.0 + let (rxx, rxy, rxz, ryx, ryy, ryz, rzx, rzy, rzz) = rotationEntries( + rotationAboutX: helmert.rx.arcSecondsToRadians, + rotationAboutY: helmert.ry.arcSecondsToRadians, + rotationAboutZ: helmert.rz.arcSecondsToRadians, + transposed: false) + ecef.x = helmert.dx + scale * (rxx * x1 + rxy * y1 + rxz * z1) + ecef.y = helmert.dy + scale * (ryx * x1 + ryy * y1 + ryz * z1) + ecef.z = helmert.dz + scale * (rzx * x1 + rzy * y1 + rzz * z1) + } + else { + let scale = 1.0 / (1.0 + helmert.scalePpm / 1_000_000.0) + let (rxx, rxy, rxz, ryx, ryy, ryz, rzx, rzy, rzz) = rotationEntries( + rotationAboutX: helmert.rx.arcSecondsToRadians, + rotationAboutY: helmert.ry.arcSecondsToRadians, + rotationAboutZ: helmert.rz.arcSecondsToRadians, + transposed: true) + ecef.x = scale * (rxx * (x1 - helmert.dx) + rxy * (y1 - helmert.dy) + rxz * (z1 - helmert.dz)) + ecef.y = scale * (ryx * (x1 - helmert.dx) + ryy * (y1 - helmert.dy) + ryz * (z1 - helmert.dz)) + ecef.z = scale * (rzx * (x1 - helmert.dx) + rzy * (y1 - helmert.dy) + rzz * (z1 - helmert.dz)) + } + + // Geocentric -> geodetic on the target ellipsoid. + let (latitude, longitude, altitude) = ecefToGeodetic( + x: ecef.x, + y: ecef.y, + z: ecef.z, + ellipsoid: targetEllipsoid) + + return Coordinate3D( + x: longitude, + y: latitude, + z: altitude, + m: coordinate.m, + projection: targetProjection) + } + + /// The position vector rotation matrix R = Rz * Ry * Rx, built from + /// arc seconds; optionally transposed. + static func rotationEntries( + rotationAboutX rxRadians: Double, + rotationAboutY ryRadians: Double, + rotationAboutZ rzRadians: Double, + transposed: Bool + ) -> (rxx: Double, rxy: Double, rxz: Double, ryx: Double, ryy: Double, ryz: Double, rzx: Double, rzy: Double, rzz: Double) { + let cx = cos(rxRadians) + let sx = sin(rxRadians) + let cy = cos(ryRadians) + let sy = sin(ryRadians) + let cz = cos(rzRadians) + let sz = sin(rzRadians) + + let rxx = cz * cy + let rxy = cz * sy * sx - sz * cx + let rxz = cz * sy * cx + sz * sx + let ryx = sz * cy + let ryy = sz * sy * sx + cz * cx + let ryz = sz * sy * cx - cz * sx + let rzx = -sy + let rzy = cy * sx + let rzz = cy * cx + + if transposed { + return (rxx, ryx, rzx, rxy, ryy, rzy, rxz, ryz, rzz) + } + return (rxx, rxy, rxz, ryx, ryy, ryz, rzx, rzy, rzz) + } + + /// Geodetic to geocentric conversion on an arbitrary ellipsoid. + private static func geodeticToEcef( + latitude phi: Double, + longitude lambda: Double, + height: Double, + ellipsoid: Ellipsoid + ) -> (x: Double, y: Double, z: Double) { + let e2 = ellipsoid.eccentricitySquared + + let sinPhi = sin(phi) + let cosPhi = cos(phi) + let radiusSquared = 1.0 - e2 * sinPhi * sinPhi + let nRadius = ellipsoid.semiMajorAxis / sqrt(radiusSquared) + + let x = (nRadius + height) * cosPhi * cos(lambda) + let y = (nRadius + height) * cosPhi * sin(lambda) + let z = ((1.0 - e2) * nRadius + height) * sinPhi + + return (x, y, z) + } + + /// Geocentric to geodetic conversion on an arbitrary ellipsoid + /// (iterative Bowring method, as used for EPSG:4978). + private static func ecefToGeodetic( + x: Double, + y: Double, + z: Double, + ellipsoid: Ellipsoid + ) -> (latitude: CLLocationDegrees, longitude: CLLocationDegrees, altitude: CLLocationDistance) { + let a = ellipsoid.semiMajorAxis + let e2 = ellipsoid.eccentricitySquared + + let p = sqrt(x * x + y * y) + + guard p > GISTool.intersectionEpsilon else { + let latitude = z >= 0.0 ? 90.0 : -90.0 + let altitude = abs(z) - a * (1.0 - e2) + return (latitude, 0.0, altitude) + } + + var phi = atan2(z, p * (1.0 - e2)) + var height: Double = 0.0 + for _ in 0 ..< 10 { + let sinPhi = sin(phi) + let cosPhi = cos(phi) + let radius = a / sqrt(1.0 - e2 * sinPhi * sinPhi) + height = p / cosPhi - radius + phi = atan2(z * (radius + height), p * ((1.0 - e2) * radius + height)) + } + + let latitude = Swift.min(90.0, Swift.max(-90.0, phi * 180.0 / .pi)) + let longitude = atan2(y, x) * 180.0 / .pi + + return (latitude, longitude, height) + } + +} diff --git a/Sources/GISTools/GeoJson/Nad27Definition.swift b/Sources/GISTools/GeoJson/Nad27Definition.swift new file mode 100644 index 00000000..fe8c4bfd --- /dev/null +++ b/Sources/GISTools/GeoJson/Nad27Definition.swift @@ -0,0 +1,51 @@ +#if canImport(CoreLocation) +import CoreLocation +#endif +import Foundation + +/// NAD27 math for EPSG:4267 (NAD 1927 geodetic coordinates on the +/// Clarke 1866 ellipsoid). +/// +/// Transformations use the EPSG "NAD27 to WGS 84 (4)" Helmert +/// transformation (translation-only, stated accuracy approximately 10 m). +/// Sub-meter accuracy for NAD27 requires NADCON grid shifts, which the +/// library does not incorporate (see the grid shift support issue). +struct Nad27Definition: ProjectionDefinition { + + private static let helmert = HelmertTransformation.nad27 + + var projection: Projection { Projection(uncheckedSrid: 4267, definition: self) } + var kind: ProjectionKind { .geographic } + var datum: Datum { .nad27 } + + var wraparoundExtent: Double? { 180.0 } + + /// The expression is global. Unbounded like EPSG:4326 (the extent of + /// the underlying datum's realization is not enforced). + 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", "NAD27"], + ["GEOGCS", "NAD_1927"], + ["GEOGCS", "NAD 1927"], + ["GEOGCS", "North_American_1927"], + ["GEOGCS", "North American 1927"], + ] + } + + func forward(_ coordinate: Coordinate3D) -> Coordinate3D { + Self.helmert.transform(wgs84ToDatum: coordinate, projection: projection) + } + + func inverse(_ coordinate: Coordinate3D) -> Coordinate3D { + Self.helmert.transform(datumToWgs84: coordinate) + } + +} diff --git a/Sources/GISTools/GeoJson/Osgb1936BngDefinition.swift b/Sources/GISTools/GeoJson/Osgb1936BngDefinition.swift new file mode 100644 index 00000000..74a32d76 --- /dev/null +++ b/Sources/GISTools/GeoJson/Osgb1936BngDefinition.swift @@ -0,0 +1,95 @@ +#if canImport(CoreLocation) +import CoreLocation +#endif +import Foundation + +/// British National Grid math for EPSG:27700 (OSGB 1936 / British National +/// Grid): a modified transverse Mercator projection on the Airy 1830 +/// ellipsoid with the OS-documented approximate Helmert transformation +/// (stated accuracy approximately 2 m; the OS-recommended route for +/// sub-meter accuracy is the OSTN15 grid, which the library does not +/// incorporate: see the grid shift support issue). +/// +/// Reference parameters (EPSG:27700 / OS Guide, appendix A): +/// TM origin at 49°N / 2°W, scale factor 0.9996012717, false easting +/// 400_000 m, false northing −100_000 m. +/// +/// The valid extent is `nil` (unbounded): coordinates south of the false +/// northing (northern Scilly, southwest of the projection origin) carry +/// *negative* northings, so a nominal sector extent would reject or clamp +/// valid values. +struct Osgb1936BngDefinition: ProjectionDefinition { + + private static let helmert = HelmertTransformation.osgb1936 + + /// The projection parameters of the British National Grid. + private static let transverseMercator = TransverseMercatorMath( + ellipsoid: .airy1830, + latitudeOfOrigin: 49.0, + longitudeOfOrigin: -2.0, + scaleFactor: 0.9996012717, + falseEasting: 400_000.0, + falseNorthing: -100_000.0) + + var projection: Projection { Projection(uncheckedSrid: 27700, definition: self) } + var kind: ProjectionKind { .planar } + var datum: Datum { .osgb1936 } + + /// British National Grid coordinates never wrap. + var wraparoundExtent: Double? { + nil + } + + /// Unbounded: the National Grid sector includes negative northings + /// southwest of the projection origin. + var validExtent: ProjectionExtent? { + nil + } + + var worldBoundingBox: BoundingBox? { + nil + } + + var wktMatchers: [[String]] { + [ + ["PROJCS", "British National Grid"], + ["PROJCS", "OSGB_1936_British_National_Grid"], + ["OSGB 1936 / British National Grid"], + ] + } + + func forward(_ coordinate: Coordinate3D) -> Coordinate3D { + let osgb36 = Self.helmert.transform(wgs84ToDatum: coordinate) + let (easting, northing) = Self.transverseMercator.forward( + latitude: osgb36.latitude, + longitude: osgb36.longitude) + + return Coordinate3D( + x: easting, + y: northing, + z: coordinate.altitude, + m: coordinate.m, + projection: projection) + } + + func inverse(_ coordinate: Coordinate3D) -> Coordinate3D { + let (latitude, longitude) = Self.transverseMercator.inverse( + x: coordinate.longitude, + y: coordinate.latitude) + + let osgb36 = Coordinate3D( + latitude: latitude, + longitude: longitude, + altitude: coordinate.altitude, + m: coordinate.m) + + let wgs84 = Self.helmert.transform(datumToWgs84: osgb36) + + return Coordinate3D( + latitude: wgs84.latitude, + longitude: wgs84.longitude, + altitude: wgs84.altitude, + m: coordinate.m) + } + +} diff --git a/Sources/GISTools/GeoJson/Osgb1936Definition.swift b/Sources/GISTools/GeoJson/Osgb1936Definition.swift new file mode 100644 index 00000000..e7803d04 --- /dev/null +++ b/Sources/GISTools/GeoJson/Osgb1936Definition.swift @@ -0,0 +1,48 @@ +#if canImport(CoreLocation) +import CoreLocation +#endif +import Foundation + +/// OSGB 1936 math for EPSG:4277 (OSGB 1936 geodetic coordinates on the +/// Airy 1830 ellipsoid). +/// +/// Transformations use the OS-documented approximate Helmert +/// transformation (EPSG:1314, "OSGB 1936 to WGS 84 (6)", stated accuracy +/// approximately 2 m). The OS-recommended route for sub-meter accuracy is +/// the OSTN15 grid, which the library does not incorporate (see the grid +/// shift support issue). +struct Osgb1936Definition: ProjectionDefinition { + + private static let helmert = HelmertTransformation.osgb1936 + + var projection: Projection { Projection(uncheckedSrid: 4277, definition: self) } + var kind: ProjectionKind { .geographic } + var datum: Datum { .osgb1936 } + + 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", "OSGB 1936"], + ["GEOGCS", "OSGB_1936"], + ["GEOGCS", "OSGB36"], + ] + } + + func forward(_ coordinate: Coordinate3D) -> Coordinate3D { + Self.helmert.transform(wgs84ToDatum: coordinate, projection: projection) + } + + func inverse(_ coordinate: Coordinate3D) -> Coordinate3D { + Self.helmert.transform(datumToWgs84: coordinate) + } + +} diff --git a/Sources/GISTools/GeoJson/Projection.swift b/Sources/GISTools/GeoJson/Projection.swift index 58dd22f7..00f08adf 100644 --- a/Sources/GISTools/GeoJson/Projection.swift +++ b/Sources/GISTools/GeoJson/Projection.swift @@ -32,6 +32,15 @@ public struct Projection: /// EPSG:32662 - WGS 84 / Plate Carree, equirectangular /// (https://epsg.io/32662). public static let epsg32662 = Projection.builtin(srid: 32_662) + /// EPSG:4258 - ETRS89 geodetic degrees (https://epsg.io/4258). + /// Effectively coincides with WGS84 at meter accuracy. + public static let epsg4258 = Projection.builtin(srid: 4258) + /// EPSG:4267 - NAD27 geodetic degrees (https://epsg.io/4267). + public static let epsg4267 = Projection.builtin(srid: 4267) + /// EPSG:4277 - OSGB 1936 geodetic degrees (https://epsg.io/4277). + public static let epsg4277 = Projection.builtin(srid: 4277) + /// EPSG:27700 - OSGB 1936 / British National Grid (https://epsg.io/27700). + public static let epsg27700 = Projection.builtin(srid: 27_700) // UTM zones (northern EPSG:32601-32660, southern EPSG:32701-32760). /// EPSG:32601 - UTM zone 1N (https://epsg.io/32601). public static let epsg32601 = Projection.builtin(srid: 32601) @@ -457,6 +466,12 @@ extension Projection { kind != .undefined } + /// The geodetic datum of the receiver's coordinate frame, resolved from + /// the projection's captured definition. + public var datum: Datum { + definition.datum + } + /// The absolute value beyond which the receiver's horizontal axis wraps /// around (±180° for EPSG:4326, ±`originShift` meters for EPSG:3857), /// or `nil` if the horizontal axis does not wrap. diff --git a/Sources/GISTools/GeoJson/ProjectionDefinition.swift b/Sources/GISTools/GeoJson/ProjectionDefinition.swift index 5a365e3b..958b66d9 100644 --- a/Sources/GISTools/GeoJson/ProjectionDefinition.swift +++ b/Sources/GISTools/GeoJson/ProjectionDefinition.swift @@ -15,6 +15,9 @@ protocol ProjectionDefinition: Sendable { /// The semantic category of the projection. var kind: ProjectionKind { get } + /// The geodetic datum of the projection's coordinate frame. + var datum: Datum { get } + /// The valid coordinate extent of the projection, or `nil` if the /// coordinate space is unbounded (e.g. geocentric or no SRID). var validExtent: ProjectionExtent? { get } @@ -53,6 +56,17 @@ protocol ProjectionDefinition: Sendable { } +extension ProjectionDefinition { + + /// The geodetic datum of the projection's coordinate frame. + /// + /// Default is ``Datum/wgs84``; datum-capable projections override this. + var datum: Datum { + .wgs84 + } + +} + /// Registry providing ``ProjectionDefinition`` implementations for /// registered projections. /// @@ -71,6 +85,10 @@ enum ProjectionRegistry { private static let epsg4978Definition = Epsg4978Definition() private static let epsg3395Definition = Epsg3395Definition() private static let epsg32662Definition = Epsg32662Definition() + private static let nad27Definition = Nad27Definition() + private static let etrs89Definition = Etrs89Definition() + private static let osgb1936Definition = Osgb1936Definition() + private static let osgb1936BngDefinition = Osgb1936BngDefinition() /// All built-in definitions: the 6 base definitions and all UTM zones, /// keyed by canonical SRID. @@ -85,6 +103,10 @@ enum ProjectionRegistry { 4978: epsg4978Definition, 3395: epsg3395Definition, 32_662: epsg32662Definition, + 4258: etrs89Definition, + 4267: nad27Definition, + 4277: osgb1936Definition, + 27700: osgb1936BngDefinition, ] for srid in 32_601 ... 32_660 { @@ -110,8 +132,12 @@ enum ProjectionRegistry { private static let builtinWktDefinitions: [any ProjectionDefinition] = [ epsg3857Definition, + osgb1936BngDefinition, epsg3395Definition, epsg32662Definition, + etrs89Definition, + nad27Definition, + osgb1936Definition, epsg4978Definition, epsg4326Definition, ] diff --git a/Sources/GISTools/GeoJson/TransverseMercatorMath.swift b/Sources/GISTools/GeoJson/TransverseMercatorMath.swift new file mode 100644 index 00000000..0be9fa5b --- /dev/null +++ b/Sources/GISTools/GeoJson/TransverseMercatorMath.swift @@ -0,0 +1,164 @@ +#if canImport(CoreLocation) +import CoreLocation +#endif +import Foundation + +/// Internal transverse Mercator math, parameterized by ellipsoid and +/// projection origin (Snyder formulas, *Map Projections - A Working +/// Manual*, USGS PP 1395). +/// +/// Used by ``UtmDefinition`` (EPSG:32xxx, origin at the equator) and by the +/// datument-capable built-in projected CRSs (EPSG:27700 etc.), and exposed +/// publicly through ``CustomProjection/transverseMercator(datum:helmert:)``. +struct TransverseMercatorMath: Sendable { + + /// The ellipsoid the projection is computed on. + let ellipsoid: Ellipsoid + + /// The latitude of the projection's origin, in degrees. + let latitudeOfOrigin: Double + + /// The longitude of the projection's central meridian, in degrees. + let longitudeOfOrigin: Double + + /// The central scale factor (e.g. 0.9996 for the UTM convention). + let scaleFactor: Double + + /// The false easting applied to x values, in meters. + let falseEasting: Double + + /// The false northing applied to y values. + let falseNorthing: Double + + init( + ellipsoid: Ellipsoid, + latitudeOfOrigin: Double = 0.0, + longitudeOfOrigin: Double, + scaleFactor: Double, + falseEasting: Double, + falseNorthing: Double = 0.0 + ) { + self.ellipsoid = ellipsoid + self.latitudeOfOrigin = latitudeOfOrigin + self.longitudeOfOrigin = longitudeOfOrigin + self.scaleFactor = scaleFactor + self.falseEasting = falseEasting + self.falseNorthing = falseNorthing + } + + // MARK: - Forward + + /// Snyder forward formulas (equations 8-9..8-25): geodetic degrees into + /// grid coordinates in meters. + func forward( + latitude: Double, + longitude: Double + ) -> (x: Double, y: Double) { + let a = ellipsoid.semiMajorAxis + let e2 = ellipsoid.eccentricitySquared + let e1sq = e2 / (1.0 - e2) + + let phi = latitude.degreesToRadians + let deltaLambda = longitude.degreesToRadians - longitudeOfOrigin.degreesToRadians + + let sinPhi = sin(phi) + let cosPhi = cos(phi) + let tanPhi = tan(phi) + + let radiusSquared = 1.0 - e2 * sinPhi * sinPhi + let nRadius = a / sqrt(radiusSquared) + let tRadius = tanPhi * tanPhi + let cTerm = e1sq * cosPhi * cosPhi + let aTerm = deltaLambda * cosPhi + + let meridianArc = Self.meridianArc(phi, ellipsoid: ellipsoid) + + let easting = falseEasting + + scaleFactor * nRadius * ( + aTerm + + (1.0 - tRadius + cTerm) * pow(aTerm, 3) / 6.0 + + (5.0 - 18.0 * tRadius + tRadius * tRadius + 72.0 * cTerm - 58.0 * e1sq) * pow(aTerm, 5) / 120.0 + ) + + let northing = falseNorthing + scaleFactor * ( + meridianArc - Self.meridianArc(latitudeOfOrigin.degreesToRadians, ellipsoid: ellipsoid) + + nRadius * tanPhi * ( + aTerm * aTerm / 2.0 + + (5.0 - tRadius + 9.0 * cTerm + 4.0 * cTerm * cTerm) * pow(aTerm, 4) / 24.0 + + (61.0 - 58.0 * tRadius + tRadius * tRadius + 600.0 * cTerm - 330.0 * e1sq) * pow(aTerm, 6) / 720.0 + ) + ) + + return (easting, northing) + } + + // MARK: - Inverse + + /// Snyder inverse formulas: grid coordinates in meters into geodetic + /// degrees. + func inverse( + x: Double, + y: Double + ) -> (latitude: Double, longitude: Double) { + let a = ellipsoid.semiMajorAxis + let e2 = ellipsoid.eccentricitySquared + let e1sq = e2 / (1.0 - e2) + + let x = x - falseEasting + let y = y - falseNorthing + + let meridianArc = y / scaleFactor + + Self.meridianArc(latitudeOfOrigin.degreesToRadians, ellipsoid: ellipsoid) + let mu = meridianArc / (a * (1.0 - e2 / 4.0 - 3.0 * e2 * e2 / 64.0 - 5.0 * pow(e2, 3) / 256.0)) + + let e1 = (1.0 - sqrt(1.0 - e2)) / (1.0 + sqrt(1.0 - e2)) + + let phi1 = mu + + (3.0 * e1 / 2.0 - 27.0 * pow(e1, 3) / 32.0) * sin(2.0 * mu) + + (21.0 * e1 * e1 / 16.0 - 55.0 * pow(e1, 4) / 32.0) * sin(4.0 * mu) + + (151.0 * pow(e1, 3) / 96.0) * sin(6.0 * mu) + + (1097.0 * pow(e1, 4) / 512.0) * sin(8.0 * mu) + + let radiusSquared1 = 1.0 - e2 * sin(phi1) * sin(phi1) + let nRadius1 = a / sqrt(radiusSquared1) + let tRadius1 = tan(phi1) * tan(phi1) + let cTerm1 = e1sq * cos(phi1) * cos(phi1) + let meridianRadius1 = a * (1.0 - e2) / pow(radiusSquared1, 1.5) + let dTerm = x / (nRadius1 * scaleFactor) + + let latitude = phi1 + - (nRadius1 * tan(phi1) / meridianRadius1) * ( + dTerm * dTerm / 2.0 + - (5.0 + 3.0 * tRadius1 + 10.0 * cTerm1 - 4.0 * cTerm1 * cTerm1 - 9.0 * e1sq) * pow(dTerm, 4) / 24.0 + + (61.0 + 90.0 * tRadius1 + 298.0 * cTerm1 + 45.0 * tRadius1 * tRadius1 - 252.0 * e1sq - 3.0 * cTerm1 * cTerm1) * pow(dTerm, 6) / 720.0 + ) + + let lambdaTerm = (dTerm + - (1.0 + 2.0 * tRadius1 + cTerm1) * pow(dTerm, 3) / 6.0 + + (5.0 - 2.0 * cTerm1 + 28.0 * tRadius1 - 3.0 * cTerm1 * cTerm1 + 8.0 * e1sq + 24.0 * tRadius1 * tRadius1) * pow(dTerm, 5) / 120.0) + / cos(phi1) + let longitude = longitudeOfOrigin.degreesToRadians + lambdaTerm + + return (latitude * 180.0 / .pi, longitude * 180.0 / .pi) + } + + // MARK: - Meridian arc + + /// Snyder meridian arc length (equation 3-21) from the equator to + /// `phi` (radians), in meters. + static func meridianArc( + _ phi: Double, + ellipsoid: Ellipsoid + ) -> Double { + let a = ellipsoid.semiMajorAxis + let e2 = ellipsoid.eccentricitySquared + + return a * ( + (1.0 - e2 / 4.0 - 3.0 * e2 * e2 / 64.0 - 5.0 * pow(e2, 3) / 256.0) * phi + - (3.0 * e2 / 8.0 + 3.0 * e2 * e2 / 32.0 + 45.0 * pow(e2, 3) / 1024.0) * sin(2.0 * phi) + + (15.0 * e2 * e2 / 256.0 + 45.0 * pow(e2, 3) / 1024.0) * sin(4.0 * phi) + - (35.0 * pow(e2, 3) / 3072.0) * sin(6.0 * phi) + ) + } + +} diff --git a/Sources/GISTools/GeoJson/UtmDefinition.swift b/Sources/GISTools/GeoJson/UtmDefinition.swift index 13f7911a..2b38b381 100644 --- a/Sources/GISTools/GeoJson/UtmDefinition.swift +++ b/Sources/GISTools/GeoJson/UtmDefinition.swift @@ -127,41 +127,27 @@ struct UtmDefinition: ProjectionDefinition { // MARK: - Conversions + // The UTM convention applies the Snyder transverse Mercator formulas on + // the WGS84 ellipsoid with the 0.9996 scale factor and the per-zone + // false easting/northing; delegating to the parameterized math keeps + // the UTM zones and other TM-based CRSs (EPSG:27700 etc.) on the same + // implementation. + + /// The Snyder transverse Mercator setup of the zone. + var transverseMercator: TransverseMercatorMath { + TransverseMercatorMath( + ellipsoid: Ellipsoid.wgs84, + latitudeOfOrigin: 0.0, + longitudeOfOrigin: centralMeridian * 180.0 / .pi, + scaleFactor: Self.k0, + falseEasting: Self.falseEasting, + falseNorthing: hemisphere.falseNorthing) + } + func forward(_ coordinate: Coordinate3D) -> Coordinate3D { - let a = GISTool.equatorialRadius - let e2 = GISTool.wgs84EccentricitySquared - let e1sq = e2 / (1.0 - e2) - - let phi = coordinate.latitude.degreesToRadians - let deltaLambda = (coordinate.longitude.degreesToRadians) - centralMeridian - - let sinPhi = sin(phi) - let cosPhi = cos(phi) - let tanPhi = tan(phi) - - let radiusSquared = 1.0 - e2 * sinPhi * sinPhi - let nRadius = a / sqrt(radiusSquared) - let tRadius = tanPhi * tanPhi - let cTerm = e1sq * cosPhi * cosPhi - let aTerm = deltaLambda * cosPhi - - let meridianArc = Self.meridianArc(phi) - - let easting = Self.falseEasting - + Self.k0 * nRadius * ( - aTerm - + (1.0 - tRadius + cTerm) * pow(aTerm, 3) / 6.0 - + (5.0 - 18.0 * tRadius + tRadius * tRadius + 72.0 * cTerm - 58.0 * e1sq) * pow(aTerm, 5) / 120.0 - ) - - let northing = Self.k0 * ( - meridianArc - + nRadius * tanPhi * ( - aTerm * aTerm / 2.0 - + (5.0 - tRadius + 9.0 * cTerm + 4.0 * cTerm * cTerm) * pow(aTerm, 4) / 24.0 - + (61.0 - 58.0 * tRadius + tRadius * tRadius + 600.0 * cTerm - 330.0 * e1sq) * pow(aTerm, 6) / 720.0 - ) - ) + hemisphere.falseNorthing + let (easting, northing) = transverseMercator.forward( + latitude: coordinate.latitude, + longitude: coordinate.longitude) return Coordinate3D( x: easting, @@ -172,47 +158,13 @@ struct UtmDefinition: ProjectionDefinition { } func inverse(_ coordinate: Coordinate3D) -> Coordinate3D { - let a = GISTool.equatorialRadius - let e2 = GISTool.wgs84EccentricitySquared - let e1sq = e2 / (1.0 - e2) - - let x = coordinate.longitude - Self.falseEasting - let y = coordinate.latitude - hemisphere.falseNorthing - - let meridianArc = y / Self.k0 - let mu = meridianArc / (a * (1.0 - e2 / 4.0 - 3.0 * e2 * e2 / 64.0 - 5.0 * pow(e2, 3) / 256.0)) - - let e1 = (1.0 - sqrt(1.0 - e2)) / (1.0 + sqrt(1.0 - e2)) - - let phi1 = mu - + (3.0 * e1 / 2.0 - 27.0 * pow(e1, 3) / 32.0) * sin(2.0 * mu) - + (21.0 * e1 * e1 / 16.0 - 55.0 * pow(e1, 4) / 32.0) * sin(4.0 * mu) - + (151.0 * pow(e1, 3) / 96.0) * sin(6.0 * mu) - + (1097.0 * pow(e1, 4) / 512.0) * sin(8.0 * mu) - - let radiusSquared1 = 1.0 - e2 * sin(phi1) * sin(phi1) - let nRadius1 = a / sqrt(radiusSquared1) - let tRadius1 = tan(phi1) * tan(phi1) - let cTerm1 = e1sq * cos(phi1) * cos(phi1) - let meridianRadius1 = a * (1.0 - e2) / pow(radiusSquared1, 1.5) - let dTerm = x / (nRadius1 * Self.k0) - - let latitude = phi1 - - (nRadius1 * tan(phi1) / meridianRadius1) * ( - dTerm * dTerm / 2.0 - - (5.0 + 3.0 * tRadius1 + 10.0 * cTerm1 - 4.0 * cTerm1 * cTerm1 - 9.0 * e1sq) * pow(dTerm, 4) / 24.0 - + (61.0 + 90.0 * tRadius1 + 298.0 * cTerm1 + 45.0 * tRadius1 * tRadius1 - 252.0 * e1sq - 3.0 * cTerm1 * cTerm1) * pow(dTerm, 6) / 720.0 - ) - - let lambdaTerm = (dTerm - - (1.0 + 2.0 * tRadius1 + cTerm1) * pow(dTerm, 3) / 6.0 - + (5.0 - 2.0 * cTerm1 + 28.0 * tRadius1 - 3.0 * cTerm1 * cTerm1 + 8.0 * e1sq + 24.0 * tRadius1 * tRadius1) * pow(dTerm, 5) / 120.0) - / cos(phi1) - let longitude = centralMeridian + lambdaTerm + let (latitude, longitude) = transverseMercator.inverse( + x: coordinate.longitude, + y: coordinate.latitude) return Coordinate3D( - latitude: latitude * 180.0 / .pi, - longitude: longitude * 180.0 / .pi, + latitude: latitude, + longitude: longitude, altitude: coordinate.altitude, m: coordinate.m) } diff --git a/Tests/GISToolsTests/Algorithms/DatumCrsTests.swift b/Tests/GISToolsTests/Algorithms/DatumCrsTests.swift new file mode 100644 index 00000000..c21ca60d --- /dev/null +++ b/Tests/GISToolsTests/Algorithms/DatumCrsTests.swift @@ -0,0 +1,165 @@ +#if canImport(CoreLocation) +import CoreLocation +#endif +import Foundation +@testable import GISTools +import Testing + +/// Algorithm round trips and smoke tests for EPSG:4267 (NAD27 geodetic). +struct Nad27CrsTests { + + /// Validates forward conversions against reference values computed + /// independently with pyproj/PROJ (via the Helmert-only + /// "NAD27 to WGS 84 (4)" transformation). + @Test + func referenceValues() async throws { + let fixtures: [(lat: Double, lon: Double, expectLat: Double, expectLon: Double)] = [ + (40.0, -100.0, 39.999990516, -99.999582393), + (35.0, -105.0, 34.999924993, -104.999461737), + (45.0, -68.0, 45.000057988, -68.000666074), + (61.0, -149.0, 61.000455150, -148.997389203), + (25.0, -80.0, 24.999595268, -80.000197174), + ] + + for fixture in fixtures { + let projected = Coordinate3D(latitude: fixture.lat, longitude: fixture.lon) + .projected(to: .epsg4267) + #expect(abs(projected.latitude - fixture.expectLat) < 0.000000001) + #expect(abs(projected.longitude - fixture.expectLon) < 0.000000001) + } + } + + /// Validates round trips through all projections with the datum pivot. + @Test + func roundTrips() async throws { + let fixtures: [(lat: Double, lon: Double)] = [ + (40.0, -100.0), (35.0, -105.0), (45.0, -68.0), (61.0, -149.0), (25.0, -80.0), + ] + + for fixture in fixtures { + let original = Coordinate3D(latitude: fixture.lat, longitude: fixture.lon, altitude: 300.0, m: 3.0) + + let nad27 = original.projected(to: .epsg4267) + #expect(nad27.projection == .epsg4267) + #expect(Projection.epsg4267.datum == Datum.nad27) + + let back = nad27.projected(to: .epsg4326) + #expect(abs(back.latitude - fixture.lat) < 0.000002) + #expect(abs(back.longitude - fixture.lon) < 0.000002) + + // Through the WGS84 pivot into the other datum-capable CRSs. + for target in [Projection.epsg4326, .epsg4258] { + let projected = original.projected(to: target) + let homeAgain = projected.projected(to: .epsg4326) + #expect(abs(homeAgain.latitude - fixture.lat) < 0.000002) + #expect(abs(homeAgain.longitude - fixture.lon) < 0.000002) + } + } + + // Round trips through the British National Grid (in-domain points: + // the TM series is only valid on British territory, like any + // transverse Mercator-based CRS beyond its zone). + for fixture in [ + (lat: 52.0, lon: -1.0), (lat: 55.5, lon: -3.0), (lat: 58.5, lon: -5.0), + ] { + let original = Coordinate3D(latitude: fixture.lat, longitude: fixture.lon) + let target = original.projected(to: .epsg27700) + #expect(target.projection == .epsg27700) + let homeAgain = target.projected(to: .epsg4326) + #expect(abs(homeAgain.latitude - fixture.lat) < 0.000002) + #expect(abs(homeAgain.longitude - fixture.lon) < 0.000002) + } + } + + /// Validates BNG (27700) E/N coordinates against reference values + /// computed independently with pyproj/PROJ. + @Test + func britishNationalGrid() async throws { + let fixtures: [(lat: Double, lon: Double, expectE: Double, expectN: Double)] = [ + (52.0, -1.0, 468_748.5557, 233_978.3388), + (55.5, -3.0, 336_926.9327, 623_367.5946), + (49.95, -6.3, 91_648.8496, 14_401.7063), + (58.5, -5.0, 225_281.6028, 960_772.9988), + (50.5, 0.8, 598_665.9800, 70_433.2732), + (59.0, -2.0, 400_099.7506, 1_012_543.3994), + ] + + for fixture in fixtures { + let bng = Coordinate3D(latitude: fixture.lat, longitude: fixture.lon) + .projected(to: .epsg27700) + #expect(abs(bng.longitude - fixture.expectE) < 0.001) + #expect(abs(bng.latitude - fixture.expectN) < 0.001) + } + + // Southwest grid values are valid (EPSG:27700 is unbounded, so + // even out-of-sector coordinates like the far southwest stay valid + // instead of being clamped away). + let scilly = Coordinate3D(latitude: 49.95, longitude: -6.3).projected(to: .epsg27700) + #expect(abs(scilly.longitude - 91_648.8496) < 0.001) + #expect(abs(scilly.latitude - 14_401.7063) < 0.001) + #expect(scilly.isValid) + let back = scilly.projected(to: .epsg4326) + #expect(abs(back.latitude - 49.95) < 0.000002) + #expect(abs(back.longitude - -6.3) < 0.000002) + } + + /// Validates algorithm behavior in the geographic datum CRSs + /// (kind-based dispatch: Haversine distance, antimeridian handling). + @Test + func algorithmSmokeTests() async throws { + // Geographic kind: haversine distance between London/Paris areas. + let london = Coordinate3D(x: -0.12, y: 51.5, projection: .epsg4267) + let paris = Coordinate3D(x: 2.35, y: 48.85, projection: .epsg4267) + let distance = london.distance(from: paris) + #expect(distance > 300_000.0) + #expect(distance < 400_000.0) + + // Planar kind for BNG. + let bng1 = Coordinate3D(x: 468_748.5557, y: 233_978.3388, projection: .epsg27700) + let bng2 = Coordinate3D(x: 468_748.5557, y: 234_978.3388, projection: .epsg27700) + #expect(abs(bng1.distance(from: bng2) - 1_000.0) < 0.000001) + + // Antimeridian handling in the datum geographic CRSs (wrap at ±180°). + #expect(Projection.epsg4267.wraparoundExtent == 180.0) + #expect(Projection.epsg27700.wraparoundExtent == nil) + + let normalized = Coordinate3D( + x: -100.0, + y: 40.0, + projection: .epsg4267).normalized() + #expect(abs(normalized.longitude - -100.0) < 0.00000001) + } + + /// Validates WKB round trips with embedded SRIDs for the new projections. + @Test + func wkbRoundTrip() async throws { + for (projection, x, y) in [ + (Projection.epsg4267, -99.999582393, 39.999990516), + (Projection.epsg4258, 2.35, 48.85), + (Projection.epsg4277, -0.998473603, 51.999560347), + (Projection.epsg27700, 468_748.5557, 233_978.3388), + ] { + let point = Point(Coordinate3D(x: x, y: y, projection: projection)) + let wkb = try #require(WKBCoder.encode(geometry: point, targetProjection: projection)) + let decoded = try #require(WKBCoder.decode( + wkb: wkb, + sourceSrid: nil, + targetProjection: projection) as? Point) + + #expect(decoded.projection == projection) + #expect(abs(decoded.coordinate.x - x) < 0.000001) + #expect(abs(decoded.coordinate.y - y) < 0.000001) + } + } + + /// Validates WKT recognition for the new datum-capable projections. + @Test + func wktMatching() async throws { + #expect(Projection(wkt: #"GEOGCS["NAD 1927"]"#) == .epsg4267) + #expect(Projection(wkt: #"GEOGCS["GCS_North_American_1927"]"#) == .epsg4267) + #expect(Projection(wkt: #"GEOGCS["ETRS 1989"]"#) == .epsg4258) + #expect(Projection(wkt: #"GEOGCS["OSGB_1936"]"#) == .epsg4277) + #expect(Projection(wkt: #"PROJCS["OSGB 1936 / British National Grid",...Transverse_Mercator...]"#) == .epsg27700) + } + +} diff --git a/Tests/GISToolsTests/GeoJson/DatumTests.swift b/Tests/GISToolsTests/GeoJson/DatumTests.swift new file mode 100644 index 00000000..1dee665e --- /dev/null +++ b/Tests/GISToolsTests/GeoJson/DatumTests.swift @@ -0,0 +1,52 @@ +#if canImport(CoreLocation) +import CoreLocation +#endif +import Foundation +@testable import GISTools +import Testing + +/// Tests for the datum metadata model. +struct DatumTests { + + /// Validates the well-known ellipsoid reference values (EPSG registry). + @Test + func ellipsoidReferenceValues() async throws { + #expect(Ellipsoid.wgs84.semiMajorAxis == 6_378_137.0) + #expect(Ellipsoid.wgs84.inverseFlattening == 298.257223563) + #expect(abs(Ellipsoid.wgs84.semiMinorAxis - 6_356_752.31424518) < 0.000001) + + #expect(Ellipsoid.grs80.semiMajorAxis == 6_378_137.0) + #expect(Ellipsoid.grs80.inverseFlattening == 298.257222101) + + #expect(Ellipsoid.clarke1866.semiMajorAxis == 6_378_206.4) + #expect(Ellipsoid.clarke1866.inverseFlattening == 294.978698213898) + + #expect(Ellipsoid.airy1830.semiMajorAxis == 6_377_563.396) + #expect(Ellipsoid.airy1830.inverseFlattening == 299.3249646) + } + + /// Validates the well-known datum constants. + @Test + func datumConstants() async throws { + #expect(Datum.wgs84.ellipsoid == Ellipsoid.wgs84) + #expect(Datum.etrs89.ellipsoid == Ellipsoid.grs80) + #expect(Datum.nad27.ellipsoid == Ellipsoid.clarke1866) + #expect(Datum.osgb1936.ellipsoid == Ellipsoid.airy1830) + } + + /// Validates that the built-in projections report their datums. + @Test + func projectionDatums() async throws { + #expect(Projection.epsg4326.datum == Datum.wgs84) + #expect(Projection.epsg3857.datum == Datum.wgs84) + #expect(Projection.epsg4978.datum == Datum.wgs84) + #expect(Projection.epsg3395.datum == Datum.wgs84) + #expect(Projection.epsg32662.datum == Datum.wgs84) + + #expect(Projection.epsg4258.datum == Datum.etrs89) + #expect(Projection.epsg4267.datum == Datum.nad27) + #expect(Projection.epsg4277.datum == Datum.osgb1936) + #expect(Projection.epsg27700.datum == Datum.osgb1936) + } + +} diff --git a/Tests/GISToolsTests/GeoJson/HelmertTransformationTests.swift b/Tests/GISToolsTests/GeoJson/HelmertTransformationTests.swift new file mode 100644 index 00000000..9d5b0da4 --- /dev/null +++ b/Tests/GISToolsTests/GeoJson/HelmertTransformationTests.swift @@ -0,0 +1,65 @@ +#if canImport(CoreLocation) +import CoreLocation +#endif +import Foundation +@testable import GISTools +import Testing + +/// Tests for the Helmert datum transformation helper. +/// +/// The reference values were computed independently with pyproj/PROJ 9 +/// (EPSG "NAD27 to WGS 84 (4)", accuracy ~10 m, and the OS-documented +/// "OSGB 1936 to WGS 84 (6)" transformation, EPSG:1314, accuracy ~2 m). +struct HelmertTransformationTests { + + /// Validates WGS84 → datum coordinates against the reference values. + @Test + func wgs84ToDatum() async throws { + // Reference: pyproj, Helmert inverse of "NAD27 to WGS 84 (4)". + // (wgs84 lat/lon in, NAD27 lat/lon out) + let helmertNad27 = HelmertTransformation.nad27 + let nad27 = helmertNad27.transform(wgs84ToDatum: Coordinate3D(latitude: 40.0, longitude: -100.0)) + #expect(abs(nad27.latitude - 39.999990516) < 0.00000001) + #expect(abs(nad27.longitude - -99.999582393) < 0.00000001) + + // Reference: pyproj, Helmert inverse of "OSGB 1936 to WGS 84 (6)" + // (zero height) + let helmertOsgb = HelmertTransformation.osgb1936 + let osgb = helmertOsgb.transform(wgs84ToDatum: Coordinate3D(latitude: 52.0, longitude: -1.0)) + #expect(abs(osgb.latitude - 51.999560347) < 0.00000001) + #expect(abs(osgb.longitude - -0.998473603) < 0.00000001) + } + + /// Validates the datum → WGS84 direction against the reference values. + @Test + func datumToWgs84() async throws { + let helmertNad27 = HelmertTransformation.nad27 + let wgs84 = helmertNad27.transform(datumToWgs84: Coordinate3D(latitude: 39.999990516, longitude: -99.999582393)) + #expect(abs(wgs84.latitude - 40.0) < 0.00000001) + #expect(abs(wgs84.longitude - -100.0) < 0.00000001) + + let helmertOsgb = HelmertTransformation.osgb1936 + let wgs84Uk = helmertOsgb.transform(datumToWgs84: Coordinate3D(latitude: 51.999560347, longitude: -0.998473603)) + #expect(abs(wgs84Uk.latitude - 52.0) < 0.00000001) + #expect(abs(wgs84Uk.longitude - -1.0) < 0.0000001) + } + + /// Validates the rotation-free and full 7-parameter matrix versions + /// round trip through both directions. + @Test + func roundTrips() async throws { + // NAD27 has no rotations/scale; OSGB applies the full 7 parameters. + for (helmert, testCoordinate) in [ + (HelmertTransformation.nad27, Coordinate3D(latitude: 40.0, longitude: -100.0, altitude: 300.0)), + (HelmertTransformation.osgb1936, Coordinate3D(latitude: 52.5, longitude: -1.5, altitude: 250.0)), + ] { + let datumFrame = helmert.transform(wgs84ToDatum: testCoordinate) + let back = helmert.transform(datumToWgs84: datumFrame) + #expect(abs(back.latitude - testCoordinate.latitude) < 0.000000001) + #expect(abs(back.longitude - testCoordinate.longitude) < 0.000000001) + #expect(abs((back.altitude ?? 0.0) - (testCoordinate.altitude ?? 0.0)) < 0.00000001) + #expect(back.m == testCoordinate.m) + } + } + +} diff --git a/Tests/GISToolsTests/GeoJson/ProjectionWktMatcherTests.swift b/Tests/GISToolsTests/GeoJson/ProjectionWktMatcherTests.swift index 6f856ed8..437ec3f0 100644 --- a/Tests/GISToolsTests/GeoJson/ProjectionWktMatcherTests.swift +++ b/Tests/GISToolsTests/GeoJson/ProjectionWktMatcherTests.swift @@ -27,7 +27,8 @@ struct ProjectionWktMatcherTests { @Test func unrecognized() async throws { #expect(Projection(wkt: "PROJCS[\"Something Else\"]") == nil) - #expect(Projection(wkt: "GEOGCS[\"ETRS89\"]") == nil) + // ETRS89 is a registered projection with its own WKT fragments. + #expect(Projection(wkt: "GEOGCS[\"ETRS89\"]") == .epsg4258) #expect(Projection(wkt: "") == nil) } From 5f0cd252f488ebff6342e621a3276640c5430600 Mon Sep 17 00:00:00 2001 From: Thomas Rasch Date: Tue, 15 Sep 2026 12:45:49 +0200 Subject: [PATCH 2/2] Updated AGENTS.md --- AGENTS.md | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/AGENTS.md b/AGENTS.md index 2b2b1b92..2e2be353 100644 --- a/AGENTS.md +++ b/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)