Skip to content

Route projections through a ProjectionDefinition registry #238

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
trasch wants to merge 2 commits into from
Closed

Route projections through a ProjectionDefinition registry #238

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4a05232
Route projections through a ProjectionDefinition registry
trasch Sep 14, 2026
5e639b0
NoSRID coordinates are copied verbatim
trasch Sep 14, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 3 additions & 16 deletions Sources/GISTools/Algorithms/Random.swift
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,24 +118,11 @@ extension BoundingBox {

/// The world bounding box projected to the given CRS.
private static func worldBox(projection: Projection) -> BoundingBox {
switch projection {
case .epsg4326:
guard let worldBox = ProjectionRegistry.definition(for: projection).worldBoundingBox else {
// No defined world extent: fall back to the EPSG:4326 world box.
return BoundingBox.world
case .epsg3857:
let s = GISTool.originShift
return BoundingBox(
southWest: Coordinate3D(x: -s, y: -s, projection: .epsg3857),
northEast: Coordinate3D(x: s, y: s, projection: .epsg3857))
case .epsg4978:
let r = GISTool.equatorialRadius
return BoundingBox(
southWest: Coordinate3D(x: -r, y: -r, z: -r, projection: .epsg4978),
northEast: Coordinate3D(x: r, y: r, z: r, projection: .epsg4978))
case .noSRID:
return BoundingBox(
southWest: Coordinate3D(x: -180.0, y: -90.0, projection: .noSRID),
northEast: Coordinate3D(x: 180.0, y: 90.0, projection: .noSRID))
}
return worldBox
}

/// A random position within the world bounding box.
Expand Down
31 changes: 10 additions & 21 deletions Sources/GISTools/GeoJson/BoundingBox.swift
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -295,28 +295,17 @@ public struct BoundingBox:
// that build a single axis-aligned envelope (e.g. PostGIS
// `ST_MakeEnvelope`) or clip features to the box. Clamp the wrapped
// coordinate back to the world boundary instead.
switch projection {
case .epsg3857:
let upper = GISTool.originShift
let lower = -GISTool.originShift
if ne.longitude < southWest.longitude { ne.longitude = upper }
if sw.longitude > northEast.longitude { sw.longitude = lower }
if let wrapExtent = projection.wraparoundExtent {
if ne.longitude < southWest.longitude { ne.longitude = wrapExtent }
if sw.longitude > northEast.longitude { sw.longitude = -wrapExtent }
// Final safety clamp for any remaining out-of-range values.
sw.longitude = min(upper, max(lower, sw.longitude))
ne.longitude = min(upper, max(lower, ne.longitude))
sw.latitude = min(upper, max(lower, sw.latitude))
ne.latitude = min(upper, max(lower, ne.latitude))

case .epsg4326:
if ne.longitude < southWest.longitude { ne.longitude = 180.0 }
if sw.longitude > northEast.longitude { sw.longitude = -180.0 }
sw.longitude = min(180.0, max(-180.0, sw.longitude))
ne.longitude = min(180.0, max(-180.0, ne.longitude))
sw.latitude = min(90.0, max(-90.0, sw.latitude))
ne.latitude = min(90.0, max(-90.0, ne.latitude))

case .epsg4978, .noSRID:
break
sw.longitude = min(wrapExtent, max(-wrapExtent, sw.longitude))
ne.longitude = min(wrapExtent, max(-wrapExtent, ne.longitude))
}

if let extent = ProjectionRegistry.definition(for: projection).validExtent {
sw.latitude = min(extent.maxY, max(extent.minY, sw.latitude))
ne.latitude = min(extent.maxY, max(extent.minY, ne.latitude))
}

return BoundingBox(southWest: sw, northEast: ne)
Expand Down
274 changes: 38 additions & 236 deletions Sources/GISTools/GeoJson/Coordinate3D.swift
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -407,39 +407,27 @@ extension Coordinate3D {
self = self.clamped()
}

/// Clamped to [[-180,-90], [180,90]].
/// Clamped to the valid coordinate range of the receiver's projection.
///
/// For EPSG:4326 and EPSG:3857, clamps the longitude to the valid range.
/// For EPSG:4978 (ECEF) this is a no-op.
/// For projections with a defined extent (e.g. EPSG:4326, EPSG:3857),
/// clamps both axes to that extent. For unbounded projections
/// (e.g. EPSG:4978 (ECEF)) this is a no-op.
///
/// - Returns: A copy clamped to the valid coordinate range
public func clamped() -> Coordinate3D {
switch projection {
case .epsg3857:
guard longitude < -GISTool.originShift || longitude > GISTool.originShift
|| latitude < -GISTool.originShift || latitude > GISTool.originShift
else { return self }

return Coordinate3D(
x: min(GISTool.originShift, max(-GISTool.originShift, longitude)),
y: min(GISTool.originShift, max(-GISTool.originShift, latitude)),
z: altitude,
m: m)

case .epsg4326:
guard longitude < -180.0 || longitude > 180.0
|| latitude < -90.0 || latitude > 90.0
else { return self }
guard let extent = ProjectionRegistry.definition(for: projection).validExtent
else { return self }

return Coordinate3D(
latitude: min(90.0, max(-90.0, latitude)),
longitude: min(180.0, max(-180.0, longitude)),
altitude: altitude,
m: m)
guard longitude < extent.minX || longitude > extent.maxX
|| latitude < extent.minY || latitude > extent.maxY
else { return self }

default:
return self
}
return Coordinate3D(
x: min(extent.maxX, max(extent.minX, longitude)),
y: min(extent.maxY, max(extent.minY, latitude)),
z: altitude,
m: m,
projection: projection)
}

}
Expand All @@ -450,65 +438,32 @@ extension Coordinate3D: Projectable {

/// Reproject this coordinate.
///
/// Conversions are routed through the EPSG:4326 pivot: the receiver is
/// converted to EPSG:4326 first, then into the target projection.
///
/// Coordinates without an SRID are always copied verbatim, regardless of
/// the target projection: their values are meaningless without a CRS, so
/// no transformation math is applied to them.
///
/// - Parameters:
/// - newProjection: The target projection
/// - Returns: A new coordinate in the requested projection
public func projected(to newProjection: Projection) -> Coordinate3D {
guard newProjection != projection else { return self }

switch newProjection {
case .epsg3857:
switch projection {
case .epsg3857:
return self
case .epsg4326, .noSRID:
return Coordinate3D(
x: longitudeProjected(to: newProjection),
y: latitudeProjected(to: newProjection),
z: altitude,
m: m,
projection: newProjection)
case .epsg4978:
return projected(to: .epsg4326).projected(to: .epsg3857)
}

case .epsg4326:
switch projection {
case .epsg4326:
return self
case .epsg3857, .noSRID:
return Coordinate3D(
x: longitudeProjected(to: newProjection),
y: latitudeProjected(to: newProjection),
z: altitude,
m: m,
projection: newProjection)
case .epsg4978:
let (lat, lon, alt) = Coordinate3D.ecefToGeodetic(
x: longitude, y: latitude, z: altitude ?? 0.0)
return Coordinate3D(latitude: lat, longitude: lon, altitude: alt, m: m)
}

case .epsg4978:
switch projection {
case .epsg4978:
return self
case .epsg4326, .noSRID:
let (x, y, z) = Coordinate3D.geodeticToEcef(
latitude: latitude, longitude: longitude, altitude: altitude ?? 0.0)
return Coordinate3D(x: x, y: y, z: z, m: m, projection: .epsg4978)
case .epsg3857:
return projected(to: .epsg4326).projected(to: .epsg4978)
}

case .noSRID:
// Coordinates without an SRID are copied verbatim into any target:
// without a CRS there is nothing to transform from.
if projection == .noSRID || newProjection == .noSRID {
return Coordinate3D(
x: longitude,
y: latitude,
z: altitude,
m: m,
projection: .noSRID)
projection: newProjection)
}

let pivot = ProjectionRegistry.definition(for: projection).inverse(self)
return ProjectionRegistry.definition(for: newProjection).forward(pivot)
}

/// Project this coordinate's latitude.
Expand All @@ -517,50 +472,13 @@ extension Coordinate3D: Projectable {
/// - newProjection: The target projection
/// - Returns: The latitude value in the requested projection
public func latitudeProjected(to newProjection: Projection) -> Double {
switch newProjection {
case .epsg3857:
switch projection {
case .epsg3857, .noSRID:
return latitude
case .epsg4326:
var y: Double = log(tan((90.0 + latitude) * Double.pi / 360.0)) / (Double.pi / 180.0)
y *= GISTool.originShift / 180.0
return y
case .epsg4978:
return projected(to: .epsg4326).latitudeProjected(to: .epsg3857)
}

case .epsg4326:
switch projection {
case .epsg4326, .noSRID:
return latitude
case .epsg3857:
return 180.0 / Double.pi * (2.0 * atan(exp((latitude / GISTool.originShift) * 180.0 * Double.pi / 180.0)) - Double.pi / 2.0)
case .epsg4978:
let (lat, _, _) = Coordinate3D.ecefToGeodetic(
x: longitude,
y: latitude,
z: altitude ?? 0.0)
return lat
}

case .epsg4978:
switch projection {
case .epsg4978, .noSRID:
return latitude
case .epsg4326:
let (_, y, _) = Coordinate3D.geodeticToEcef(
latitude: latitude,
longitude: longitude,
altitude: altitude ?? 0.0)
return y
case .epsg3857:
return projected(to: .epsg4326).latitudeProjected(to: .epsg4978)
}

case .noSRID:
// Coordinates without an SRID are treated as already being in the
// target projection for per-axis lookups.
if projection == .noSRID {
return latitude
}

return projected(to: newProjection).latitude
}

/// Project this coordinate's longitude.
Expand All @@ -569,129 +487,13 @@ extension Coordinate3D: Projectable {
/// - newProjection: The target projection
/// - Returns: The longitude value in the requested projection
public func longitudeProjected(to newProjection: Projection) -> Double {
switch newProjection {
case .epsg3857:
switch projection {
case .epsg3857, .noSRID:
return longitude
case .epsg4326:
return longitude * GISTool.originShift / 180.0
case .epsg4978:
return projected(to: .epsg4326).longitudeProjected(to: .epsg3857)
}

case .epsg4326:
switch projection {
case .epsg4326, .noSRID:
return longitude
case .epsg3857:
return (longitude / GISTool.originShift) * 180.0
case .epsg4978:
let (_, lon, _) = Coordinate3D.ecefToGeodetic(
x: longitude,
y: latitude,
z: altitude ?? 0.0)
return lon
}

case .epsg4978:
switch projection {
case .epsg4978, .noSRID:
return longitude
case .epsg4326:
let (x, _, _) = Coordinate3D.geodeticToEcef(
latitude: latitude,
longitude: longitude,
altitude: altitude ?? 0.0)
return x
case .epsg3857:
return projected(to: .epsg4326).longitudeProjected(to: .epsg4978)
}

case .noSRID:
// Coordinates without an SRID are treated as already being in the
// target projection for per-axis lookups.
if projection == .noSRID {
return longitude
}
}

}

// MARK: - EPSG:4978 (ECEF) Helpers

extension Coordinate3D {

/// Convert geodetic (EPSG:4326) to geocentric (EPSG:4978) coordinates.
///
/// Uses the WGS84 ellipsoid: a = 6,378,137 m, 1/f = 298.257223563.
///
/// - Parameters:
/// - latitude: Latitude in degrees
/// - longitude: Longitude in degrees
/// - altitude: Height above ellipsoid in meters
/// - Returns: ECEF X, Y, Z in meters
fileprivate static func geodeticToEcef(
latitude: Double,
longitude: Double,
altitude: Double
) -> (x: Double, y: Double, z: Double) {
let a = GISTool.equatorialRadius
let e2 = GISTool.wgs84EccentricitySquared

let phi = latitude * .pi / 180.0
let lambda = longitude * .pi / 180.0
let h = altitude

let sinPhi = sin(phi)
let N = a / sqrt(1.0 - e2 * sinPhi * sinPhi)

let x = (N + h) * cos(phi) * cos(lambda)
let y = (N + h) * cos(phi) * sin(lambda)
let z = ((1.0 - e2) * N + h) * sinPhi

return (x, y, z)
}

/// Convert geocentric (EPSG:4978) to geodetic (EPSG:4326) coordinates.
///
/// Uses the iterative Bowring method (typically converges in ~3 iterations).
///
/// - Parameters:
/// - x: ECEF X in meters
/// - y: ECEF Y in meters
/// - z: ECEF Z in meters
/// - Returns: Latitude (degrees), longitude (degrees), altitude (meters)
fileprivate static func ecefToGeodetic(
x: Double,
y: Double,
z: Double
) -> (latitude: Double, longitude: Double, altitude: Double) {
let a = GISTool.equatorialRadius
let e2 = GISTool.wgs84EccentricitySquared

let p = sqrt(x * x + y * y)

guard p > GISTool.intersectionEpsilon else {
let lat = z >= 0.0 ? 90.0 : -90.0
let h = abs(z) - a * (1.0 - e2)
return (lat, 0.0, h)
}

var phi = atan2(z, p * (1.0 - e2))
var h: Double = 0.0
for _ in 0 ..< 10 {
let sinPhi = sin(phi)
let cosPhi = cos(phi)
let N = a / sqrt(1.0 - e2 * sinPhi * sinPhi)
h = p / cosPhi - N
phi = atan2(z * (N + h), p * ((1.0 - e2) * N + h))
}

// Clamp latitude to the valid geodetic range. The iterative solver
// can diverge for points near the geocenter (huge negative altitude),
// producing |lat| > 90.
let lat = Swift.min(90.0, Swift.max(-90.0, phi * 180.0 / .pi))
let lon = atan2(y, x) * 180.0 / .pi

return (lat, lon, h)
return projected(to: newProjection).longitude
}

}
Expand Down
Loading