diff --git a/.github/workflows/swift.yml b/.github/workflows/swift.yml index c03dc88..38157cb 100644 --- a/.github/workflows/swift.yml +++ b/.github/workflows/swift.yml @@ -20,7 +20,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - swift: ["5.10", "5.9", "5.7"] + swift: ["6.0", "5.10", "5.9"] container: image: swift:${{ matrix.swift }} steps: diff --git a/Examples/Cassini/ContentView+Model.swift b/Examples/Cassini/ContentView+Model.swift index ef6111f..d053a72 100644 --- a/Examples/Cassini/ContentView+Model.swift +++ b/Examples/Cassini/ContentView+Model.swift @@ -22,6 +22,7 @@ extension ContentView { case mercator case gallPeters case equalEarth + case naturalEarth case orthographic case azimuthal @@ -85,6 +86,8 @@ extension ContentView { projection = Projections.GallPeters(reference: reference) case .equalEarth: projection = Projections.EqualEarth(reference: reference) + case .naturalEarth: + projection = Projections.NaturalEarth(reference: reference) case .orthographic: projection = Projections.Orthographic(reference: reference) case .azimuthal: diff --git a/Sources/GeoProjector/ProjectionMode.swift b/Sources/GeoProjector/ProjectionMode.swift index d8749e6..ef39ae0 100644 --- a/Sources/GeoProjector/ProjectionMode.swift +++ b/Sources/GeoProjector/ProjectionMode.swift @@ -13,6 +13,7 @@ public enum ProjectionMode: String, Codable, Sendable { case mercator case equirectangular case equalEarth + case naturalEarth case gallPeters case azimuthal case orthographic @@ -22,7 +23,9 @@ public enum ProjectionMode: String, Codable, Sendable { switch self { case .automatic where max(boundingBox.longitudeSpan, boundingBox.latitudeSpan) < 180, .orthographic: return Projections.Orthographic(reference: boundingBox.center) - case .equalEarth, .automatic: + case .naturalEarth, .automatic: + return Projections.NaturalEarth(reference: boundingBox.center) + case .equalEarth: return Projections.EqualEarth(reference: boundingBox.center) case .mercator: @@ -39,7 +42,9 @@ public enum ProjectionMode: String, Codable, Sendable { public func resolve(for center: GeoJSON.Position = .init(latitude: 0, longitude: 0)) -> Projection { switch self { - case .automatic, .equalEarth: + case .automatic, .naturalEarth: + return Projections.NaturalEarth(reference: center) + case .equalEarth: return Projections.EqualEarth(reference: center) case .orthographic: return Projections.Orthographic(reference: center) diff --git a/Sources/GeoProjector/Projections+Pseudocylindrical.swift b/Sources/GeoProjector/Projections+Pseudocylindrical.swift index a449020..012bee0 100644 --- a/Sources/GeoProjector/Projections+Pseudocylindrical.swift +++ b/Sources/GeoProjector/Projections+Pseudocylindrical.swift @@ -92,4 +92,75 @@ extension Projections { } + /// Compromise projection that's optimised to look nice as a small map + /// + /// See https://www.shadedrelief.com/NE_proj/index.html + public struct NaturalEarth: Projection { + private static let A0 = 0.8707 + private static let A1 = -0.131979 + private static let A2 = -0.013791 + private static let A3 = 0.003971 + private static let A4 = -0.001529 + private static let B0 = 1.007226 + private static let B1 = 0.015085 + private static let B2 = -0.044475 + private static let B3 = 0.028874 + private static let B4 = -0.005916 + private static let MAX_Y = 0.8707 * 0.52 * .pi + + public let reference: Point + + public let projectionSize: Size + + public let mapBounds: MapBounds + + public init(reference: Point) { + self.reference = reference + + // Calculate projection size based on equations + self.projectionSize = .init( + width: 2 * .pi * Self.A0, + height: 2 * Self.MAX_Y + ) + + let inputCorners: [Point] = [ + .init(x: -1 * .pi, y: .pi / 2), + .init(x: .pi, y: .pi / 2), + .init(x: .pi, y: -1 * .pi / 2), + .init(x: -1 * .pi, y: -1 * .pi / 2), + .init(x: -1 * .pi, y: .pi / 2), + ] + + let boundPoints = zip(inputCorners.dropLast(), inputCorners.dropFirst()) + .reduce(into: [Point]()) { acc, next in + acc.append(Self.project(next.0)) + acc.append(contentsOf: Interpolator.interpolate(from: next.0, to: next.1, maxDiff: 0.0025, projector: Self.project(_:)).map(\.1)) + acc.append(Self.project(next.1)) + } + self.mapBounds = .bezier(boundPoints) + } + + public func willWrap(_ point: Point) -> Bool { + Projections.willWrap(point, reference: reference) + } + + public func project(_ point: Point) -> Point? { + let adjusted = Projections.adjust(point, reference: reference) + return Self.project(adjusted) + } + + private static func project(_ point: Point) -> Point { + let phi = point.y + let lam = point.x + + let phi2 = phi * phi + let phi4 = phi2 * phi2 + + let x = lam * (A0 + phi2 * (A1 + phi2 * (A2 + phi4 * phi2 * (A3 + phi2 * A4)))) + let y = phi * (B0 + phi2 * (B1 + phi4 * (B2 + B3 * phi2 + B4 * phi4))) + + return .init(x: x, y: y) + } + } + } diff --git a/Tests/GeoProjectorTests/NaturalEarthTests.swift b/Tests/GeoProjectorTests/NaturalEarthTests.swift new file mode 100644 index 0000000..17b9999 --- /dev/null +++ b/Tests/GeoProjectorTests/NaturalEarthTests.swift @@ -0,0 +1,80 @@ +// +// NaturalEarthTests.swift +// GeoProjector +// +// Created by Adrian Schönig on 3/5/2025. +// + +#if canImport(Testing) +import Testing + +@testable import GeoProjector + +struct NaturalEarthTests { + + @Test func matchesReferencePoints() async throws { + // Reference coordinates projected with the polynomial Natural Earth projection (lon, lat, X, Y): + // The radius of the sphere is 6371008.7714 + let reference = """ + 0.0 0.0 0.0 0.0 + 0.0 22.5 0.0 2525419.569383768 + 0.0 45.0 0.0 5052537.389973222 + 0.0 67.5 0.0 7400065.6562573705 + 0.0 90.0 0.0 9062062.394736718 + 45.0 0.0 4356790.016612169 0.0 + 45.0 22.5 4253309.544984069 2525419.569383768 + 45.0 45.0 3924521.5829515466 5052537.389973222 + 45.0 67.5 3354937.47115583 7400065.6562573705 + 45.0 90.0 2397978.2448443635 9062062.394736718 + 90.0 0.0 8713580.033224339 0.0 + 90.0 22.5 8506619.089968137 2525419.569383768 + 90.0 45.0 7849043.165903093 5052537.389973222 + 90.0 67.5 6709874.94231166 7400065.6562573705 + 90.0 90.0 4795956.489688727 9062062.394736718 + 135.0 0.0 1.3070370049836507E7 0.0 + 135.0 22.5 1.2759928634952208E7 2525419.569383768 + 135.0 45.0 1.177356474885464E7 5052537.389973222 + 135.0 67.5 1.0064812413467491E7 7400065.6562573705 + 135.0 90.0 7193934.734533091 9062062.394736718 + 180.0 0.0 1.7427160066448677E7 0.0 + 180.0 22.5 1.7013238179936275E7 2525419.569383768 + 180.0 45.0 1.5698086331806187E7 5052537.389973222 + 180.0 67.5 1.341974988462332E7 7400065.6562573705 + 180.0 90.0 9591912.979377454 9062062.394736718 + """ + + // Parse reference data + let lines = reference.split(separator: "\n").map { $0.trimmingCharacters(in: .whitespaces) } + let referencePoints = lines.compactMap { line -> (lon: Double, lat: Double, x: Double, y: Double)? in + let components = line.split(separator: " ").compactMap { Double($0.trimmingCharacters(in: .whitespaces)) } + guard components.count == 4 else { return nil } + return (lon: components[0], lat: components[1], x: components[2], y: components[3]) + } + + // Earth radius used in reference data + let earthRadius = 6371008.7714 + + // Create projection + let projection = Projections.NaturalEarth(reference: .init(x: 0, y: 0)) + + // Test reference points + for point in referencePoints { + // Convert degrees to radians for input + let lonRad = point.lon * .pi / 180.0 + let latRad = point.lat * .pi / 180.0 + + // Project the point + let projected = try #require(projection.project(.init(x: lonRad, y: latRad)), "Failed to project point at lon: \(point.lon), lat: \(point.lat)") + + // Scale to match reference earth radius + let scaledX = projected.x * earthRadius + let scaledY = projected.y * earthRadius + + // Compare with reference data within tolerance + #expect(abs(scaledX - point.x) < 0.1, "X coordinate doesn't match for lon: \(point.lon), lat: \(point.lat). Expected: \(point.x), got: \(scaledX)") + #expect(abs(scaledY - point.y) < 0.1, "Y coordinate doesn't match for lon: \(point.lon), lat: \(point.lat). Expected: \(point.y), got: \(scaledY)") + } + } +} + +#endif \ No newline at end of file