From d89497bf40df9139924faba5377a35ca2f9beaa8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Bardon?= Date: Sat, 1 Oct 2022 23:41:53 +0200 Subject: [PATCH] =?UTF-8?q?=F0=9F=9A=A7=20WIP=20Introduce=20geodesic=20typ?= =?UTF-8?q?es?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../xcshareddata/xcschemes/Geodesy.xcscheme | 67 ++++ .../xcshareddata/xcschemes/WGS84.xcscheme | 79 +++++ .../xcschemes/swift-geo-Package.xcscheme | 48 +++ Package.resolved | 8 + Package.swift | 17 +- Sources/Geodesy/Conversions.swift | 159 +++++++++ Sources/Geodesy/Geodesy.swift | 326 ++++++++++++++++++ Sources/WGS84/WGS84.swift | 114 ++++++ Tests/WGS84Tests/WSG84Tests.swift | 30 ++ 9 files changed, 847 insertions(+), 1 deletion(-) create mode 100644 .swiftpm/xcode/xcshareddata/xcschemes/Geodesy.xcscheme create mode 100644 .swiftpm/xcode/xcshareddata/xcschemes/WGS84.xcscheme create mode 100644 Sources/Geodesy/Conversions.swift create mode 100644 Sources/Geodesy/Geodesy.swift create mode 100644 Sources/WGS84/WGS84.swift create mode 100644 Tests/WGS84Tests/WSG84Tests.swift diff --git a/.swiftpm/xcode/xcshareddata/xcschemes/Geodesy.xcscheme b/.swiftpm/xcode/xcshareddata/xcschemes/Geodesy.xcscheme new file mode 100644 index 0000000..beee8bf --- /dev/null +++ b/.swiftpm/xcode/xcshareddata/xcschemes/Geodesy.xcscheme @@ -0,0 +1,67 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/.swiftpm/xcode/xcshareddata/xcschemes/WGS84.xcscheme b/.swiftpm/xcode/xcshareddata/xcschemes/WGS84.xcscheme new file mode 100644 index 0000000..1229093 --- /dev/null +++ b/.swiftpm/xcode/xcshareddata/xcschemes/WGS84.xcscheme @@ -0,0 +1,79 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/.swiftpm/xcode/xcshareddata/xcschemes/swift-geo-Package.xcscheme b/.swiftpm/xcode/xcshareddata/xcschemes/swift-geo-Package.xcscheme index 8e240c7..48aa5bb 100644 --- a/.swiftpm/xcode/xcshareddata/xcschemes/swift-geo-Package.xcscheme +++ b/.swiftpm/xcode/xcshareddata/xcschemes/swift-geo-Package.xcscheme @@ -90,6 +90,34 @@ ReferencedContainer = "container:"> + + + + + + + + + + + + + + + + { + associatedtype OldCRS: CoordinateReferenceSystem + associatedtype NewCRS: CoordinateReferenceSystem + +// static var name: String { get } +// static var code: Int { get } + + static func apply, C2: Coordinate>(on coordinate: C1) -> C2 +} + +public protocol ReversibleConversion: Conversion { + static func unapply, C2: Coordinate>(from coordinate: C2) -> C1 +} + +// MARK: EPSG 9659 (Geographic 3D to 2D conversions) + +public enum EPSG9659: ReversibleConversion +where OldCRS: GeographicCRS & ThreeDimensionsCRS, + NewCRS: GeographicCRS & TwoDimensionsCRS, + NewCRS.Datum.PrimeMeridian == OldCRS.Datum.PrimeMeridian +{ + // FIXME: Static stored properties not supported in generic types +// static let name: String = "Geographic 3D to 2D conversions" +// static let code: Int = 9659 + + /// Geographic 3D to 2D conversion. + /// + /// section 4.1.4 + public static func apply, C2: Coordinate>(on coordinate: C1) -> C2 { + return C2.init(x: .init(coordinate.x), y: .init(coordinate.y)) + } + + /// Geographic 2D to 3D conversion. + /// + /// section 4.1.1 + public static func unapply, C2: Coordinate>(from coordinate: C2) -> C1 { + return C1.init(x: .init(Double(coordinate.x)), y: .init(coordinate.y), z: .init(0.0)) + } +} + +public extension ThreeDimensionsCoordinate { + func transformed(to newCRS: NewCRS.Type) -> Coordinate2DOf + where CRS: GeographicCRS, + NewCRS: GeographicCRS & TwoDimensionsCRS, + NewCRS.Datum.PrimeMeridian == CRS.Datum.PrimeMeridian + { + EPSG9659.apply(on: self) + } +} + +public extension TwoDimensionsCoordinate { + func transformed(to newCRS: NewCRS.Type) -> Coordinate3DOf + where CRS: GeographicCRS, + NewCRS: GeographicCRS & ThreeDimensionsCRS, + NewCRS.Datum.PrimeMeridian == CRS.Datum.PrimeMeridian + { + EPSG9659.unapply(from: self) + } +} + +// MARK: EPSG 9602 (Geographic/geocentric conversions) + +public enum EPSG9602: ReversibleConversion +where OldCRS: GeographicCRS & ThreeDimensionsCRS, + NewCRS: GeocentricCRS & ThreeDimensionsCRS, + NewCRS.Datum.PrimeMeridian == OldCRS.Datum.PrimeMeridian +{ + // FIXME: Static stored properties not supported in generic types +// static let name: String = "Geographic/geocentric conversions" +// static let code: Int = 9602 + + /// Geographic to geocentric conversion. + /// + /// + /// page 101 + public static func apply, C2: Coordinate>(on coordinate: C1) -> C2 { + // φ and λ are respectively the latitude and longitude (related to Greenwich) of the point + let φ = Double(coordinate.x).toRadians + let λ = Double(coordinate.y).toRadians + // Height above the ellipsoid (topographic height plus geoidal height) + let h = Double(coordinate.z) + + let a = NewCRS.Datum.Ellipsoid.semiMajorAxis + let e2 = NewCRS.Datum.Ellipsoid.e2 + + // ν = a / (1 - e^2*sin^2(φ))^0.5 + let ν = a / sqrt(1 - e2 * sin(sin(φ))) + + let x = (ν + h) * cos(φ) * cos(λ) + let y = (ν + h) * cos(φ) * sin(λ) + let z = ((1 - e2) * ν + h) * sin(φ) + + return C2.init(x: .init(x), y: .init(y), z: .init(z)) + } + + /// Geocentric to geographic conversion. + /// + /// + /// page 101 + public static func unapply, C2: Coordinate>(from coordinate: C2) -> C1 { + let x = Double(coordinate.x) + let y = Double(coordinate.y) + let z = Double(coordinate.z) + + let a = NewCRS.Datum.Ellipsoid.semiMajorAxis + let b = NewCRS.Datum.Ellipsoid.semiMinorAxis + let e2 = NewCRS.Datum.Ellipsoid.e2 + let ε = NewCRS.Datum.Ellipsoid.ε + + // p = (X^2 + Y^2)^0.5 + let p = sqrt(pow(x, 2) + pow(y, 2)) + // q = atan2[(Z a) , (p b)] + let q = atan2(z * a, p * b) + + // φ = atan2[(Z + ε b sin^3(q)) , (p - e^2 a cos^3(q))] + let φ = atan2(z + ε * b * sin(sin(sin(q))), p - e2 * a * cos(cos(cos(q)))) + // λ = atan2(Y, X) + let λ = atan2(y, x) + // ν = a / (1 - e^2*sin^2(φ))^0.5 + let ν = a / sqrt(1 - e2 * sin(sin(φ))) + // h = (p / cos φ) - ν + let h = (p / cos(φ)) - ν + + return C1.init(x: .init(φ.fromRadians), y: .init(λ.fromRadians), z: .init(h)) + } +} + +public extension ThreeDimensionsCoordinate { + func transformed(to newCRS: NewCRS.Type) -> Coordinate3DOf + where CRS: GeographicCRS, + NewCRS: GeocentricCRS & ThreeDimensionsCRS, + NewCRS.Datum.PrimeMeridian == CRS.Datum.PrimeMeridian + { + EPSG9602.apply(on: self) + } + + func transformed(to newCRS: NewCRS.Type) -> Coordinate3DOf + where CRS: GeocentricCRS, + NewCRS: GeographicCRS & ThreeDimensionsCRS, + NewCRS.Datum.PrimeMeridian == CRS.Datum.PrimeMeridian + { + EPSG9602.unapply(from: self) + } +} + +fileprivate extension FloatingPoint { + var toRadians: Self { (self / 180) * Self.pi } + var fromRadians: Self { (self / Self.pi) * 180 } +} diff --git a/Sources/Geodesy/Geodesy.swift b/Sources/Geodesy/Geodesy.swift new file mode 100644 index 0000000..599e5d0 --- /dev/null +++ b/Sources/Geodesy/Geodesy.swift @@ -0,0 +1,326 @@ +// +// Geodesy.swift +// SwiftGeo +// +// Created by Rémi Bardon on 01/10/2022. +// Copyright © 2022 Rémi Bardon. All rights reserved. +// + +import Foundation +import Tagged + +// MARK: - Coordinates + +public protocol Coordinate: Hashable, CustomStringConvertible { + associatedtype CRS: CoordinateReferenceSystem + typealias Components = CRS.CoordinateSystem.Values + + var components: Components { get } + + init(components: Components) +} + +public extension Coordinate { + var description: String { "[\(CRS.name)]\(String(describing: self.components))" } +} + +// MARK: 2D Coordinates + +public extension Coordinate where CRS.CoordinateSystem: TwoDimensionsCS { + var x: CRS.CoordinateSystem.Axis1.Value { self.components.0 } + var y: CRS.CoordinateSystem.Axis2.Value { self.components.1 } + + init( + x: CRS.CoordinateSystem.Axis1.Value, + y: CRS.CoordinateSystem.Axis2.Value + ) { + self.init(components: (x, y)) + } +} + +public protocol TwoDimensionsCoordinate: Coordinate where CRS: TwoDimensionsCRS { + typealias X = CRS.CoordinateSystem.Axis1.Value + typealias Y = CRS.CoordinateSystem.Axis2.Value + + var x: X { get } + var y: Y { get } + + init(x: X, y: Y) +} + +public extension TwoDimensionsCoordinate { + var components: Components { (self.x, self.y) } + + init(components: Components) { + self.init(x: components.0, y: components.1) + } +} + +public struct Coordinate2DOf: TwoDimensionsCoordinate where CRS: TwoDimensionsCRS { + public let x: X + public let y: Y + + public init(x: X, y: Y) { + #warning("TODO: Perform validations") + self.x = x + self.y = y + } + + public init(latitude: X, longitude: Y) { + self.init(x: latitude, y: longitude) + } +} + +// MARK: 3D Coordinates + +public extension Coordinate where CRS.CoordinateSystem: ThreeDimensionsCS { + var x: CRS.CoordinateSystem.Axis1.Value { self.components.0 } + var y: CRS.CoordinateSystem.Axis2.Value { self.components.1 } + var z: CRS.CoordinateSystem.Axis3.Value { self.components.2 } + + init( + x: CRS.CoordinateSystem.Axis1.Value, + y: CRS.CoordinateSystem.Axis2.Value, + z: CRS.CoordinateSystem.Axis3.Value + ) { + self.init(components: (x, y, z)) + } +} + +public protocol ThreeDimensionsCoordinate: Coordinate where CRS: ThreeDimensionsCRS { + typealias X = CRS.CoordinateSystem.Axis1.Value + typealias Y = CRS.CoordinateSystem.Axis2.Value + typealias Z = CRS.CoordinateSystem.Axis3.Value + + var x: X { get } + var y: Y { get } + var z: Z { get } + + init(x: X, y: Y, z: Z) +} + +public extension ThreeDimensionsCoordinate { + var components: Components { (self.x, self.y, self.z) } + + init(components: Components) { + self.init(x: components.0, y: components.1, z: components.2) + } +} + +public struct Coordinate3DOf: ThreeDimensionsCoordinate where CRS: ThreeDimensionsCRS { + public let x: X + public let y: Y + public let z: Z + + public init(x: X, y: Y, z: Z) { + #warning("TODO: Perform validations") + self.x = x + self.y = y + self.z = z + } + + public init(latitude: X, longitude: Y, altitude: Z) { + self.init(x: latitude, y: longitude, z: altitude) + } +} + +// MARK: - Coordinate Reference System (CRS/SRS) + +public protocol CoordinateReferenceSystem { + associatedtype Datum: DatumProtocol + associatedtype CoordinateSystem: Geodesy.CoordinateSystem + static var name: String { get } + static var code: Int { get } +} + +public protocol TwoDimensionsCRS: CoordinateReferenceSystem +where CoordinateSystem: TwoDimensionsCS {} +public protocol ThreeDimensionsCRS: CoordinateReferenceSystem +where CoordinateSystem: ThreeDimensionsCS {} + +public protocol GeocentricCRS: CoordinateReferenceSystem, ThreeDimensionsCRS {} +public protocol GeographicCRS: CoordinateReferenceSystem {} + +// MARK: - Datum + +public protocol DatumProtocol { + associatedtype Ellipsoid: ReferenceEllipsoid + associatedtype PrimeMeridian: MeridianProtocol +} + +public protocol DatumEnsemble: DatumProtocol { + associatedtype PrimaryMember: DatumEnsembleMember + where PrimaryMember.Ellipsoid == Self.Ellipsoid, + PrimaryMember.PrimeMeridian == Self.PrimeMeridian + static var name: String { get } + static var code: Int { get } +} +public protocol DatumEnsembleMember: DatumProtocol { + static var name: String { get } + static var code: Int { get } +} +public protocol DynamicGeodeticDatum: DatumEnsembleMember {} + +// MARK: Reference Ellipsoid + +public protocol ReferenceEllipsoid { + static var name: String { get } + static var code: Int { get } + static var semiMajorAxis: Double { get } + static var semiMinorAxis: Double { get } + static var inverseFlattening: Double { get } + static var flattening: Double { get } + /// Eccentricity of the ellipsoid. + static var eccentricity: Double { get } + /// Eccentricity of the ellipsoid, squared. + /// `e^2 = (a^2 -b^2)/a^2 = 2f -f^2` + static var e2: Double { get } + /// `ε = e^2 / (1 - e^2)` + static var ε: Double { get } +} + +public extension ReferenceEllipsoid { + static var flattening: Double { 1 / inverseFlattening } + /// Source: + static var semiMinorAxis: Double { semiMajorAxis * (1 - flattening) } + static var eccentricity: Double { sqrt(e2) } + static var e2: Double { (2 * flattening) - pow(flattening, 2) } + static var ε: Double { e2 / (1 - e2) } +} + +// MARK: Meridian + +public protocol MeridianProtocol { + static var name: String { get } + static var code: Int { get } + static var greenwichLongitude: Double { get } +} + +public enum Greenwich: MeridianProtocol { + public static let name: String = "Greenwich" + public static let code: Int = 8901 + public static let greenwichLongitude: Double = 0 +} + +// MARK: - Coordinate System (CS) + +public protocol CoordinateSystem { + associatedtype Axes + associatedtype Values + static var name: String { get } + static var code: Int { get } +} + +public protocol TwoDimensionsCS: CoordinateSystem +where Axes == (Axis1, Axis2), + Values == (Axis1.Value, Axis2.Value) +{ + associatedtype Axis1: Axis + associatedtype Axis2: Axis + // associatedtype Axes = (Axis1, Axis2) + // associatedtype Values = (Axis1.Value, Axis2.Value) +} + +public protocol ThreeDimensionsCS: CoordinateSystem +where Axes == (Axis1, Axis2, Axis3), + Values == (Axis1.Value, Axis2.Value, Axis3.Value) +{ + associatedtype Axis1: Axis + associatedtype Axis2: Axis + associatedtype Axis3: Axis + // associatedtype Axes = (Axis1, Axis2, Axis3) + // associatedtype Values = (Axis1.Value, Axis2.Value, Axis3.Value) +} + +public enum GeocentricCartesian3DCS: ThreeDimensionsCS { + public typealias Axis1 = GeocentricX + public typealias Axis2 = GeocentricY + public typealias Axis3 = GeocentricZ + public static let name: String = "Cartesian 3D CS (geocentric)" + public static let code: Int = 6500 +} + +public enum Ellipsoidal3DCS: ThreeDimensionsCS { + public typealias Axis1 = GeodeticLatitude + public typealias Axis2 = GeodeticLongitude + public typealias Axis3 = EllipsoidalHeight + public static let name: String = "Ellipsoidal 3D CS" + public static let code: Int = 6423 +} + +public enum Ellipsoidal2DCS: TwoDimensionsCS { + public typealias Axis1 = GeodeticLatitude + public typealias Axis2 = GeodeticLongitude + public static let name: String = "Ellipsoidal 2D CS" + public static let code: Int = 6422 +} + +// MARK: Axes + +public protocol Axis { + associatedtype UnitOfMeasurement: UnitOfMeasurementProtocol + associatedtype Value: BinaryFloatingPoint + static var name: String { get } + static var abbreviation: String { get } +} + +public extension Axis { + typealias Value = Tagged +} + +public enum GeocentricX: Axis { + public typealias UnitOfMeasurement = Metre + public static let name: String = "Geocentric X" + public static let abbreviation: String = "X" +} + +public enum GeocentricY: Axis { + public typealias UnitOfMeasurement = Metre + public static let name: String = "Geocentric Y" + public static let abbreviation: String = "Y" +} + +public enum GeocentricZ: Axis { + public typealias UnitOfMeasurement = Metre + public static let name: String = "Geocentric Z" + public static let abbreviation: String = "Z" +} + +public enum GeodeticLatitude: Axis { + public typealias UnitOfMeasurement = Degree + public static let name: String = "Geodetic latitude" + public static let abbreviation: String = "Lat" +} + +public enum GeodeticLongitude: Axis { + public typealias UnitOfMeasurement = Degree + public static let name: String = "Geodetic longitude" + public static let abbreviation: String = "Lon" +} + +public enum EllipsoidalHeight: Axis { + public typealias UnitOfMeasurement = Metre + public static let name: String = "Ellipsoidal height" + public static let abbreviation: String = "h" +} + +// MARK: - Units of Measurement + +public protocol UnitOfMeasurementProtocol: Hashable, Sendable { + static var name: String { get } + static var code: Int { get } +} + +public protocol Length: UnitOfMeasurementProtocol {} + +public protocol Angle: UnitOfMeasurementProtocol {} + +public enum Metre: Length { + public static let name: String = "Metre" + public static let code: Int = 9001 +} + +public enum Degree: Angle { + public static let name: String = "Degree" + public static let code: Int = 9122 +} diff --git a/Sources/WGS84/WGS84.swift b/Sources/WGS84/WGS84.swift new file mode 100644 index 0000000..ca86977 --- /dev/null +++ b/Sources/WGS84/WGS84.swift @@ -0,0 +1,114 @@ +// +// WGS84.swift +// SwiftGeo +// +// Created by Rémi Bardon on 01/10/2022. +// Copyright © 2022 Rémi Bardon. All rights reserved. +// + +import Geodesy + +// MARK: - Coordinates + +public typealias Coordinate2D = Coordinate2DOf +public typealias Coordinate3D = Coordinate3DOf + +// MARK: - Coordinate Reference System (CRS/SRS) + +/// +public enum WGS84Geographic2DCRS: GeographicCRS, TwoDimensionsCRS { + public typealias Datum = WGS84Ensemble + public typealias CoordinateSystem = Ellipsoidal2DCS + + public static let name: String = "WGS 84 (geographic 2D)" + public static let code: Int = 4326 +} + +public typealias EPSG4326 = WGS84Geographic2DCRS + +/// +public enum WGS84Geographic3DCRS: GeographicCRS, ThreeDimensionsCRS { + public typealias Datum = WGS84Ensemble + public typealias CoordinateSystem = Ellipsoidal3DCS + + public static let name: String = "WGS 84 (geographic 3D)" + public static let code: Int = 4979 +} + +public typealias EPSG4979 = WGS84Geographic3DCRS + +/// +public enum WGS84GeocentricCRS: GeocentricCRS, ThreeDimensionsCRS { + public typealias Datum = WGS84Ensemble + public typealias CoordinateSystem = GeocentricCartesian3DCS + + public static let name: String = "WGS 84 (geocentric)" + public static let code: Int = 4978 +} + +public typealias EPSG4978 = WGS84GeocentricCRS + +// MARK: - Datum + +public enum WGS84Ensemble: DatumEnsemble { + public typealias Ellipsoid = EPSG7030 + public typealias PrimeMeridian = Greenwich + public typealias PrimaryMember = Member7 + public static let name: String = "World Geodetic System 1984 ensemble" + public static let code: Int = 6326 + + public enum Member7: DatumEnsembleMember { + public typealias Ellipsoid = EPSG7030 + public typealias PrimeMeridian = Greenwich + public static let name: String = "World Geodetic System 1984 (G2139)" + public static let code: Int = 1309 + } +} + +/// +/// +public typealias EPSG6326 = WGS84Ensemble + +public typealias EPSG1309 = WGS84Ensemble.Member7 + +// MARK: Reference Ellipsoid + +/// World Geodetic System 1984 ensemble member 7 ellipsoid +/// +public enum EPSG7030: ReferenceEllipsoid { + public static let name: String = "WGS 84" + public static let code: Int = 7030 + public static let semiMajorAxis: Double = 6378137 + public static let inverseFlattening: Double = 298.257_223_563 +} + +// MARK: - Coordinate System (CS) + +public enum GeocentricCartesian3DCS: ThreeDimensionsCS { + public typealias Axis1 = GeocentricX + public typealias Axis2 = GeocentricY + public typealias Axis3 = GeocentricZ + public static let name: String = "Cartesian 3D CS (geocentric)" + public static let code: Int = 6500 +} + +public typealias EPSG6500 = GeocentricCartesian3DCS + +public enum Ellipsoidal3DCS: ThreeDimensionsCS { + public typealias Axis1 = GeodeticLatitude + public typealias Axis2 = GeodeticLongitude + public typealias Axis3 = EllipsoidalHeight + public static let name: String = "Ellipsoidal 3D CS" + public static let code: Int = 6423 +} + +public typealias EPSG6423 = Ellipsoidal3DCS + +public enum Ellipsoidal2DCS: TwoDimensionsCS { + public typealias Axis1 = GeodeticLatitude + public typealias Axis2 = GeodeticLongitude + public static let name: String = "Ellipsoidal 2D CS" + public static let code: Int = 6422 +} + +public typealias EPSG6422 = Ellipsoidal2DCS diff --git a/Tests/WGS84Tests/WSG84Tests.swift b/Tests/WGS84Tests/WSG84Tests.swift new file mode 100644 index 0000000..04458ba --- /dev/null +++ b/Tests/WGS84Tests/WSG84Tests.swift @@ -0,0 +1,30 @@ +// +// WGS84Tests.swift +// SwiftGeo +// +// Created by Rémi Bardon on 01/10/2022. +// Copyright © 2022 Rémi Bardon. All rights reserved. +// + +import Geodesy +import WGS84 +import XCTest + +final class ConversionTests: XCTestCase { + /// Example from . + func testEPSG9602() { + let c1 = Coordinate2D(x: 0, y: 90) + let c2 = c1.transformed(to: WGS84Geographic3DCRS.self) + XCTAssertEqual(c2, Coordinate3D(x: 0, y: 90, z: 0)) + let c3 = c2.transformed(to: WGS84GeocentricCRS.self) + XCTAssertEqual(c3, Coordinate3DOf( + x: 3.905482530786651e-10, + y: .init(WGS84Ensemble.Ellipsoid.semiMajorAxis), + z: 0 + )) + let c4 = c3.transformed(to: WGS84Geographic3DCRS.self) + XCTAssertEqual(c4, Coordinate3D(x: 0, y: 90, z: 0)) + let c5 = c4.transformed(to: WGS84Geographic2DCRS.self) + XCTAssertEqual(c5, c1) + } +}