Add NaturalEarth projection

.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:

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:

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)

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)
+    }
+  }
+
 }

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