Skip to content

Commit 99de9f3

Browse files
nighthawknighthawk
authored
Add NaturalEarth projection (#35)
* Add NaturalEarth projection See https://www.shadedrelief.com/NE_proj/index.html * Also add to enums and use Natural Earth as new default for '.automatic' * Bump Linux Swift versions
1 parent 77c4b71 commit 99de9f3

5 files changed

Lines changed: 162 additions & 3 deletions

File tree

‎.github/workflows/swift.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ jobs:
2020
runs-on: ubuntu-latest
2121
strategy:
2222
matrix:
23-
swift: ["5.10", "5.9", "5.7"]
23+
swift: ["6.0", "5.10", "5.9"]
2424
container:
2525
image: swift:${{ matrix.swift }}
2626
steps:

‎Examples/Cassini/ContentView+Model.swift‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@ extension ContentView {
2222
case mercator
2323
case gallPeters
2424
case equalEarth
25+
case naturalEarth
2526
case orthographic
2627
case azimuthal
2728

@@ -85,6 +86,8 @@ extension ContentView {
8586
projection = Projections.GallPeters(reference: reference)
8687
case .equalEarth:
8788
projection = Projections.EqualEarth(reference: reference)
89+
case .naturalEarth:
90+
projection = Projections.NaturalEarth(reference: reference)
8891
case .orthographic:
8992
projection = Projections.Orthographic(reference: reference)
9093
case .azimuthal:

‎Sources/GeoProjector/ProjectionMode.swift‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@ public enum ProjectionMode: String, Codable, Sendable {
1313
case mercator
1414
case equirectangular
1515
case equalEarth
16+
case naturalEarth
1617
case gallPeters
1718
case azimuthal
1819
case orthographic
@@ -22,7 +23,9 @@ public enum ProjectionMode: String, Codable, Sendable {
2223
switch self {
2324
case .automatic where max(boundingBox.longitudeSpan, boundingBox.latitudeSpan) < 180, .orthographic:
2425
return Projections.Orthographic(reference: boundingBox.center)
25-
case .equalEarth, .automatic:
26+
case .naturalEarth, .automatic:
27+
return Projections.NaturalEarth(reference: boundingBox.center)
28+
case .equalEarth:
2629
return Projections.EqualEarth(reference: boundingBox.center)
2730

2831
case .mercator:
@@ -39,7 +42,9 @@ public enum ProjectionMode: String, Codable, Sendable {
3942

4043
public func resolve(for center: GeoJSON.Position = .init(latitude: 0, longitude: 0)) -> Projection {
4144
switch self {
42-
case .automatic, .equalEarth:
45+
case .automatic, .naturalEarth:
46+
return Projections.NaturalEarth(reference: center)
47+
case .equalEarth:
4348
return Projections.EqualEarth(reference: center)
4449
case .orthographic:
4550
return Projections.Orthographic(reference: center)

‎Sources/GeoProjector/Projections+Pseudocylindrical.swift‎

Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -92,4 +92,75 @@ extension Projections {
9292

9393
}
9494

95+
/// Compromise projection that's optimised to look nice as a small map
96+
///
97+
/// See https://www.shadedrelief.com/NE_proj/index.html
98+
public struct NaturalEarth: Projection {
99+
private static let A0 = 0.8707
100+
private static let A1 = -0.131979
101+
private static let A2 = -0.013791
102+
private static let A3 = 0.003971
103+
private static let A4 = -0.001529
104+
private static let B0 = 1.007226
105+
private static let B1 = 0.015085
106+
private static let B2 = -0.044475
107+
private static let B3 = 0.028874
108+
private static let B4 = -0.005916
109+
private static let MAX_Y = 0.8707 * 0.52 * .pi
110+
111+
public let reference: Point
112+
113+
public let projectionSize: Size
114+
115+
public let mapBounds: MapBounds
116+
117+
public init(reference: Point) {
118+
self.reference = reference
119+
120+
// Calculate projection size based on equations
121+
self.projectionSize = .init(
122+
width: 2 * .pi * Self.A0,
123+
height: 2 * Self.MAX_Y
124+
)
125+
126+
let inputCorners: [Point] = [
127+
.init(x: -1 * .pi, y: .pi / 2),
128+
.init(x: .pi, y: .pi / 2),
129+
.init(x: .pi, y: -1 * .pi / 2),
130+
.init(x: -1 * .pi, y: -1 * .pi / 2),
131+
.init(x: -1 * .pi, y: .pi / 2),
132+
]
133+
134+
let boundPoints = zip(inputCorners.dropLast(), inputCorners.dropFirst())
135+
.reduce(into: [Point]()) { acc, next in
136+
acc.append(Self.project(next.0))
137+
acc.append(contentsOf: Interpolator.interpolate(from: next.0, to: next.1, maxDiff: 0.0025, projector: Self.project(_:)).map(\.1))
138+
acc.append(Self.project(next.1))
139+
}
140+
self.mapBounds = .bezier(boundPoints)
141+
}
142+
143+
public func willWrap(_ point: Point) -> Bool {
144+
Projections.willWrap(point, reference: reference)
145+
}
146+
147+
public func project(_ point: Point) -> Point? {
148+
let adjusted = Projections.adjust(point, reference: reference)
149+
return Self.project(adjusted)
150+
}
151+
152+
private static func project(_ point: Point) -> Point {
153+
let phi = point.y
154+
let lam = point.x
155+
156+
let phi2 = phi * phi
157+
let phi4 = phi2 * phi2
158+
159+
let x = lam * (A0 + phi2 * (A1 + phi2 * (A2 + phi4 * phi2 * (A3 + phi2 * A4))))
160+
let y = phi * (B0 + phi2 * (B1 + phi4 * (B2 + B3 * phi2 + B4 * phi4)))
161+
162+
return .init(x: x, y: y)
163+
}
164+
}
165+
95166
}

‎Tests/GeoProjectorTests/NaturalEarthTests.swift‎

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
//
2+
// NaturalEarthTests.swift
3+
// GeoProjector
4+
//
5+
// Created by Adrian Schönig on 3/5/2025.
6+
//
7+
8+
#if canImport(Testing)
9+
import Testing
10+
11+
@testable import GeoProjector
12+
13+
struct NaturalEarthTests {
14+
15+
@Test func matchesReferencePoints() async throws {
16+
// Reference coordinates projected with the polynomial Natural Earth projection (lon, lat, X, Y):
17+
// The radius of the sphere is 6371008.7714
18+
let reference = """
19+
0.0 0.0 0.0 0.0
20+
0.0 22.5 0.0 2525419.569383768
21+
0.0 45.0 0.0 5052537.389973222
22+
0.0 67.5 0.0 7400065.6562573705
23+
0.0 90.0 0.0 9062062.394736718
24+
45.0 0.0 4356790.016612169 0.0
25+
45.0 22.5 4253309.544984069 2525419.569383768
26+
45.0 45.0 3924521.5829515466 5052537.389973222
27+
45.0 67.5 3354937.47115583 7400065.6562573705
28+
45.0 90.0 2397978.2448443635 9062062.394736718
29+
90.0 0.0 8713580.033224339 0.0
30+
90.0 22.5 8506619.089968137 2525419.569383768
31+
90.0 45.0 7849043.165903093 5052537.389973222
32+
90.0 67.5 6709874.94231166 7400065.6562573705
33+
90.0 90.0 4795956.489688727 9062062.394736718
34+
135.0 0.0 1.3070370049836507E7 0.0
35+
135.0 22.5 1.2759928634952208E7 2525419.569383768
36+
135.0 45.0 1.177356474885464E7 5052537.389973222
37+
135.0 67.5 1.0064812413467491E7 7400065.6562573705
38+
135.0 90.0 7193934.734533091 9062062.394736718
39+
180.0 0.0 1.7427160066448677E7 0.0
40+
180.0 22.5 1.7013238179936275E7 2525419.569383768
41+
180.0 45.0 1.5698086331806187E7 5052537.389973222
42+
180.0 67.5 1.341974988462332E7 7400065.6562573705
43+
180.0 90.0 9591912.979377454 9062062.394736718
44+
"""
45+
46+
// Parse reference data
47+
let lines = reference.split(separator: "\n").map { $0.trimmingCharacters(in: .whitespaces) }
48+
let referencePoints = lines.compactMap { line -> (lon: Double, lat: Double, x: Double, y: Double)? in
49+
let components = line.split(separator: " ").compactMap { Double($0.trimmingCharacters(in: .whitespaces)) }
50+
guard components.count == 4 else { return nil }
51+
return (lon: components[0], lat: components[1], x: components[2], y: components[3])
52+
}
53+
54+
// Earth radius used in reference data
55+
let earthRadius = 6371008.7714
56+
57+
// Create projection
58+
let projection = Projections.NaturalEarth(reference: .init(x: 0, y: 0))
59+
60+
// Test reference points
61+
for point in referencePoints {
62+
// Convert degrees to radians for input
63+
let lonRad = point.lon * .pi / 180.0
64+
let latRad = point.lat * .pi / 180.0
65+
66+
// Project the point
67+
let projected = try #require(projection.project(.init(x: lonRad, y: latRad)), "Failed to project point at lon: \(point.lon), lat: \(point.lat)")
68+
69+
// Scale to match reference earth radius
70+
let scaledX = projected.x * earthRadius
71+
let scaledY = projected.y * earthRadius
72+
73+
// Compare with reference data within tolerance
74+
#expect(abs(scaledX - point.x) < 0.1, "X coordinate doesn't match for lon: \(point.lon), lat: \(point.lat). Expected: \(point.x), got: \(scaledX)")
75+
#expect(abs(scaledY - point.y) < 0.1, "Y coordinate doesn't match for lon: \(point.lon), lat: \(point.lat). Expected: \(point.y), got: \(scaledY)")
76+
}
77+
}
78+
}
79+
80+
#endif

0 commit comments

Comments
 (0)