From 41c0671ebc1bc7505327ad0cde6c3387ebdcddae Mon Sep 17 00:00:00 2001 From: "Thomas Rasch (OpenCode)" Date: Mon, 14 Sep 2026 19:37:36 +0000 Subject: [PATCH] Add all 120 UTM zone projections (EPSG:32601-32660, 32701-32760) Phase 5 of making the projection system generic. Math: standard Snyder transverse Mercator formulas on the WGS84 ellipsoid, verified independently in Python before porting - the meridian arc against Simpson integration of the meridian radius (<0.5mm divergence) and round trips over 252 points covering all zones and latitudes to 89 degrees (max error 1.06e-08 degrees). - UtmDefinition.swift: UtmHemisphere and UtmDefinition implementing both directions of the projection with a per-zone central meridian, the 0.9996 scale factor, the 500_000 false easting and per-hemisphere false northing. Sector extents follow the EPSG convention (easting 100_000-900_000, northing 0-10_000_000); the world bounding box is the zone sector and no WKT pattern is registered (UTM .prj files would need per-zone parsing). - Projection gains all 120 zone cases (northern 32601-32660, southern 32701-32760) and init?(srid:) accepts the contiguous SRID ranges. wraparoundExtent is nil (zone coordinates never wrap), kind is planar; description now defaults to the EPSG label so it scales with the SRID count. - ProjectionRegistry resolves definitions via an eagerly built dictionary; registration coverage is asserted at runtime by the new UtmDefinitionTests sweep over all UTM SRIDs. Tests: SRID-to-zone/hemisphere mapping for all 120 SRIDs, registry coverage sweep, central meridian math per zone, known forward values (zone-edge equator eastings, scaled meridian arcs, both hemispheres), hemisphere symmetry (equal eastings, exact 10_000_000 northing offset), round trips, boundingBox/contains/distance, antimeridian no-ops in date-line zones, normalize/clamp, WKB round trips and random generation within the zone sector. --- Sources/GISTools/GeoJson/Projection.swift | 274 +++++++++++++++++- .../GeoJson/ProjectionDefinition.swift | 81 ++++-- Sources/GISTools/GeoJson/UtmDefinition.swift | 210 ++++++++++++++ .../Algorithms/UtmAlgorithmTests.swift | 194 +++++++++++++ .../GeoJson/UtmDefinitionTests.swift | 82 ++++++ 5 files changed, 807 insertions(+), 34 deletions(-) create mode 100644 Sources/GISTools/GeoJson/UtmDefinition.swift create mode 100644 Tests/GISToolsTests/Algorithms/UtmAlgorithmTests.swift create mode 100644 Tests/GISToolsTests/GeoJson/UtmDefinitionTests.swift diff --git a/Sources/GISTools/GeoJson/Projection.swift b/Sources/GISTools/GeoJson/Projection.swift index d33a6b8..8336c25 100644 --- a/Sources/GISTools/GeoJson/Projection.swift +++ b/Sources/GISTools/GeoJson/Projection.swift @@ -21,6 +21,247 @@ public enum Projection: /// EPSG:32662 - WGS 84 / Plate Carree, equirectangular /// (https://epsg.io/32662). case epsg32662 = 32662 + // UTM zones (northern EPSG:32601-32660, southern EPSG:32701-32760). + /// EPSG:32601 - UTM zone 1N (https://epsg.io/32601). + case epsg32601 = 32601 + /// EPSG:32602 - UTM zone 2N (https://epsg.io/32602). + case epsg32602 = 32602 + /// EPSG:32603 - UTM zone 3N (https://epsg.io/32603). + case epsg32603 = 32603 + /// EPSG:32604 - UTM zone 4N (https://epsg.io/32604). + case epsg32604 = 32604 + /// EPSG:32605 - UTM zone 5N (https://epsg.io/32605). + case epsg32605 = 32605 + /// EPSG:32606 - UTM zone 6N (https://epsg.io/32606). + case epsg32606 = 32606 + /// EPSG:32607 - UTM zone 7N (https://epsg.io/32607). + case epsg32607 = 32607 + /// EPSG:32608 - UTM zone 8N (https://epsg.io/32608). + case epsg32608 = 32608 + /// EPSG:32609 - UTM zone 9N (https://epsg.io/32609). + case epsg32609 = 32609 + /// EPSG:32610 - UTM zone 10N (https://epsg.io/32610). + case epsg32610 = 32610 + /// EPSG:32611 - UTM zone 11N (https://epsg.io/32611). + case epsg32611 = 32611 + /// EPSG:32612 - UTM zone 12N (https://epsg.io/32612). + case epsg32612 = 32612 + /// EPSG:32613 - UTM zone 13N (https://epsg.io/32613). + case epsg32613 = 32613 + /// EPSG:32614 - UTM zone 14N (https://epsg.io/32614). + case epsg32614 = 32614 + /// EPSG:32615 - UTM zone 15N (https://epsg.io/32615). + case epsg32615 = 32615 + /// EPSG:32616 - UTM zone 16N (https://epsg.io/32616). + case epsg32616 = 32616 + /// EPSG:32617 - UTM zone 17N (https://epsg.io/32617). + case epsg32617 = 32617 + /// EPSG:32618 - UTM zone 18N (https://epsg.io/32618). + case epsg32618 = 32618 + /// EPSG:32619 - UTM zone 19N (https://epsg.io/32619). + case epsg32619 = 32619 + /// EPSG:32620 - UTM zone 20N (https://epsg.io/32620). + case epsg32620 = 32620 + /// EPSG:32621 - UTM zone 21N (https://epsg.io/32621). + case epsg32621 = 32621 + /// EPSG:32622 - UTM zone 22N (https://epsg.io/32622). + case epsg32622 = 32622 + /// EPSG:32623 - UTM zone 23N (https://epsg.io/32623). + case epsg32623 = 32623 + /// EPSG:32624 - UTM zone 24N (https://epsg.io/32624). + case epsg32624 = 32624 + /// EPSG:32625 - UTM zone 25N (https://epsg.io/32625). + case epsg32625 = 32625 + /// EPSG:32626 - UTM zone 26N (https://epsg.io/32626). + case epsg32626 = 32626 + /// EPSG:32627 - UTM zone 27N (https://epsg.io/32627). + case epsg32627 = 32627 + /// EPSG:32628 - UTM zone 28N (https://epsg.io/32628). + case epsg32628 = 32628 + /// EPSG:32629 - UTM zone 29N (https://epsg.io/32629). + case epsg32629 = 32629 + /// EPSG:32630 - UTM zone 30N (https://epsg.io/32630). + case epsg32630 = 32630 + /// EPSG:32631 - UTM zone 31N (https://epsg.io/32631). + case epsg32631 = 32631 + /// EPSG:32632 - UTM zone 32N (https://epsg.io/32632). + case epsg32632 = 32632 + /// EPSG:32633 - UTM zone 33N (https://epsg.io/32633). + case epsg32633 = 32633 + /// EPSG:32634 - UTM zone 34N (https://epsg.io/32634). + case epsg32634 = 32634 + /// EPSG:32635 - UTM zone 35N (https://epsg.io/32635). + case epsg32635 = 32635 + /// EPSG:32636 - UTM zone 36N (https://epsg.io/32636). + case epsg32636 = 32636 + /// EPSG:32637 - UTM zone 37N (https://epsg.io/32637). + case epsg32637 = 32637 + /// EPSG:32638 - UTM zone 38N (https://epsg.io/32638). + case epsg32638 = 32638 + /// EPSG:32639 - UTM zone 39N (https://epsg.io/32639). + case epsg32639 = 32639 + /// EPSG:32640 - UTM zone 40N (https://epsg.io/32640). + case epsg32640 = 32640 + /// EPSG:32641 - UTM zone 41N (https://epsg.io/32641). + case epsg32641 = 32641 + /// EPSG:32642 - UTM zone 42N (https://epsg.io/32642). + case epsg32642 = 32642 + /// EPSG:32643 - UTM zone 43N (https://epsg.io/32643). + case epsg32643 = 32643 + /// EPSG:32644 - UTM zone 44N (https://epsg.io/32644). + case epsg32644 = 32644 + /// EPSG:32645 - UTM zone 45N (https://epsg.io/32645). + case epsg32645 = 32645 + /// EPSG:32646 - UTM zone 46N (https://epsg.io/32646). + case epsg32646 = 32646 + /// EPSG:32647 - UTM zone 47N (https://epsg.io/32647). + case epsg32647 = 32647 + /// EPSG:32648 - UTM zone 48N (https://epsg.io/32648). + case epsg32648 = 32648 + /// EPSG:32649 - UTM zone 49N (https://epsg.io/32649). + case epsg32649 = 32649 + /// EPSG:32650 - UTM zone 50N (https://epsg.io/32650). + case epsg32650 = 32650 + /// EPSG:32651 - UTM zone 51N (https://epsg.io/32651). + case epsg32651 = 32651 + /// EPSG:32652 - UTM zone 52N (https://epsg.io/32652). + case epsg32652 = 32652 + /// EPSG:32653 - UTM zone 53N (https://epsg.io/32653). + case epsg32653 = 32653 + /// EPSG:32654 - UTM zone 54N (https://epsg.io/32654). + case epsg32654 = 32654 + /// EPSG:32655 - UTM zone 55N (https://epsg.io/32655). + case epsg32655 = 32655 + /// EPSG:32656 - UTM zone 56N (https://epsg.io/32656). + case epsg32656 = 32656 + /// EPSG:32657 - UTM zone 57N (https://epsg.io/32657). + case epsg32657 = 32657 + /// EPSG:32658 - UTM zone 58N (https://epsg.io/32658). + case epsg32658 = 32658 + /// EPSG:32659 - UTM zone 59N (https://epsg.io/32659). + case epsg32659 = 32659 + /// EPSG:32660 - UTM zone 60N (https://epsg.io/32660). + case epsg32660 = 32660 + /// EPSG:32701 - UTM zone 1S (https://epsg.io/32701). + case epsg32701 = 32701 + /// EPSG:32702 - UTM zone 2S (https://epsg.io/32702). + case epsg32702 = 32702 + /// EPSG:32703 - UTM zone 3S (https://epsg.io/32703). + case epsg32703 = 32703 + /// EPSG:32704 - UTM zone 4S (https://epsg.io/32704). + case epsg32704 = 32704 + /// EPSG:32705 - UTM zone 5S (https://epsg.io/32705). + case epsg32705 = 32705 + /// EPSG:32706 - UTM zone 6S (https://epsg.io/32706). + case epsg32706 = 32706 + /// EPSG:32707 - UTM zone 7S (https://epsg.io/32707). + case epsg32707 = 32707 + /// EPSG:32708 - UTM zone 8S (https://epsg.io/32708). + case epsg32708 = 32708 + /// EPSG:32709 - UTM zone 9S (https://epsg.io/32709). + case epsg32709 = 32709 + /// EPSG:32710 - UTM zone 10S (https://epsg.io/32710). + case epsg32710 = 32710 + /// EPSG:32711 - UTM zone 11S (https://epsg.io/32711). + case epsg32711 = 32711 + /// EPSG:32712 - UTM zone 12S (https://epsg.io/32712). + case epsg32712 = 32712 + /// EPSG:32713 - UTM zone 13S (https://epsg.io/32713). + case epsg32713 = 32713 + /// EPSG:32714 - UTM zone 14S (https://epsg.io/32714). + case epsg32714 = 32714 + /// EPSG:32715 - UTM zone 15S (https://epsg.io/32715). + case epsg32715 = 32715 + /// EPSG:32716 - UTM zone 16S (https://epsg.io/32716). + case epsg32716 = 32716 + /// EPSG:32717 - UTM zone 17S (https://epsg.io/32717). + case epsg32717 = 32717 + /// EPSG:32718 - UTM zone 18S (https://epsg.io/32718). + case epsg32718 = 32718 + /// EPSG:32719 - UTM zone 19S (https://epsg.io/32719). + case epsg32719 = 32719 + /// EPSG:32720 - UTM zone 20S (https://epsg.io/32720). + case epsg32720 = 32720 + /// EPSG:32721 - UTM zone 21S (https://epsg.io/32721). + case epsg32721 = 32721 + /// EPSG:32722 - UTM zone 22S (https://epsg.io/32722). + case epsg32722 = 32722 + /// EPSG:32723 - UTM zone 23S (https://epsg.io/32723). + case epsg32723 = 32723 + /// EPSG:32724 - UTM zone 24S (https://epsg.io/32724). + case epsg32724 = 32724 + /// EPSG:32725 - UTM zone 25S (https://epsg.io/32725). + case epsg32725 = 32725 + /// EPSG:32726 - UTM zone 26S (https://epsg.io/32726). + case epsg32726 = 32726 + /// EPSG:32727 - UTM zone 27S (https://epsg.io/32727). + case epsg32727 = 32727 + /// EPSG:32728 - UTM zone 28S (https://epsg.io/32728). + case epsg32728 = 32728 + /// EPSG:32729 - UTM zone 29S (https://epsg.io/32729). + case epsg32729 = 32729 + /// EPSG:32730 - UTM zone 30S (https://epsg.io/32730). + case epsg32730 = 32730 + /// EPSG:32731 - UTM zone 31S (https://epsg.io/32731). + case epsg32731 = 32731 + /// EPSG:32732 - UTM zone 32S (https://epsg.io/32732). + case epsg32732 = 32732 + /// EPSG:32733 - UTM zone 33S (https://epsg.io/32733). + case epsg32733 = 32733 + /// EPSG:32734 - UTM zone 34S (https://epsg.io/32734). + case epsg32734 = 32734 + /// EPSG:32735 - UTM zone 35S (https://epsg.io/32735). + case epsg32735 = 32735 + /// EPSG:32736 - UTM zone 36S (https://epsg.io/32736). + case epsg32736 = 32736 + /// EPSG:32737 - UTM zone 37S (https://epsg.io/32737). + case epsg32737 = 32737 + /// EPSG:32738 - UTM zone 38S (https://epsg.io/32738). + case epsg32738 = 32738 + /// EPSG:32739 - UTM zone 39S (https://epsg.io/32739). + case epsg32739 = 32739 + /// EPSG:32740 - UTM zone 40S (https://epsg.io/32740). + case epsg32740 = 32740 + /// EPSG:32741 - UTM zone 41S (https://epsg.io/32741). + case epsg32741 = 32741 + /// EPSG:32742 - UTM zone 42S (https://epsg.io/32742). + case epsg32742 = 32742 + /// EPSG:32743 - UTM zone 43S (https://epsg.io/32743). + case epsg32743 = 32743 + /// EPSG:32744 - UTM zone 44S (https://epsg.io/32744). + case epsg32744 = 32744 + /// EPSG:32745 - UTM zone 45S (https://epsg.io/32745). + case epsg32745 = 32745 + /// EPSG:32746 - UTM zone 46S (https://epsg.io/32746). + case epsg32746 = 32746 + /// EPSG:32747 - UTM zone 47S (https://epsg.io/32747). + case epsg32747 = 32747 + /// EPSG:32748 - UTM zone 48S (https://epsg.io/32748). + case epsg32748 = 32748 + /// EPSG:32749 - UTM zone 49S (https://epsg.io/32749). + case epsg32749 = 32749 + /// EPSG:32750 - UTM zone 50S (https://epsg.io/32750). + case epsg32750 = 32750 + /// EPSG:32751 - UTM zone 51S (https://epsg.io/32751). + case epsg32751 = 32751 + /// EPSG:32752 - UTM zone 52S (https://epsg.io/32752). + case epsg32752 = 32752 + /// EPSG:32753 - UTM zone 53S (https://epsg.io/32753). + case epsg32753 = 32753 + /// EPSG:32754 - UTM zone 54S (https://epsg.io/32754). + case epsg32754 = 32754 + /// EPSG:32755 - UTM zone 55S (https://epsg.io/32755). + case epsg32755 = 32755 + /// EPSG:32756 - UTM zone 56S (https://epsg.io/32756). + case epsg32756 = 32756 + /// EPSG:32757 - UTM zone 57S (https://epsg.io/32757). + case epsg32757 = 32757 + /// EPSG:32758 - UTM zone 58S (https://epsg.io/32758). + case epsg32758 = 32758 + /// EPSG:32759 - UTM zone 59S (https://epsg.io/32759). + case epsg32759 = 32759 + /// EPSG:32760 - UTM zone 60S (https://epsg.io/32760). + case epsg32760 = 32760 /// Initialize a Projection with a SRID number. /// @@ -36,6 +277,9 @@ public enum Projection: case 4978: self = .epsg4978 case 3395: self = .epsg3395 case 32662: self = .epsg32662 + // UTM zones 1-60, northern and southern. + case 32_601 ... 32_660, 32_701 ... 32_760: + self = Projection(rawValue: srid) ?? .noSRID default: return nil } } @@ -68,11 +312,7 @@ public enum Projection: public var description: String { switch self { case .noSRID: return "No SRID" - case .epsg3857: return "EPSG:3857" - case .epsg4326: return "EPSG:4326" - case .epsg4978: return "EPSG:4978" - case .epsg3395: return "EPSG:3395" - case .epsg32662: return "EPSG:32662" + default: return "EPSG:" + String(rawValue) } } @@ -89,6 +329,17 @@ extension Projection { case .epsg3857, .epsg3395, .epsg32662: .planar case .epsg4326: .geographic case .epsg4978: .geocentric + case .epsg32601, .epsg32602, .epsg32603, .epsg32604, .epsg32605, .epsg32606, .epsg32607, .epsg32608, .epsg32609, .epsg32610, .epsg32611, .epsg32612, + .epsg32613, .epsg32614, .epsg32615, .epsg32616, .epsg32617, .epsg32618, .epsg32619, .epsg32620, .epsg32621, .epsg32622, .epsg32623, .epsg32624, + .epsg32625, .epsg32626, .epsg32627, .epsg32628, .epsg32629, .epsg32630, .epsg32631, .epsg32632, .epsg32633, .epsg32634, .epsg32635, .epsg32636, + .epsg32637, .epsg32638, .epsg32639, .epsg32640, .epsg32641, .epsg32642, .epsg32643, .epsg32644, .epsg32645, .epsg32646, .epsg32647, .epsg32648, + .epsg32649, .epsg32650, .epsg32651, .epsg32652, .epsg32653, .epsg32654, .epsg32655, .epsg32656, .epsg32657, .epsg32658, .epsg32659, .epsg32660, + .epsg32701, .epsg32702, .epsg32703, .epsg32704, .epsg32705, .epsg32706, .epsg32707, .epsg32708, .epsg32709, .epsg32710, .epsg32711, .epsg32712, + .epsg32713, .epsg32714, .epsg32715, .epsg32716, .epsg32717, .epsg32718, .epsg32719, .epsg32720, .epsg32721, .epsg32722, .epsg32723, .epsg32724, + .epsg32725, .epsg32726, .epsg32727, .epsg32728, .epsg32729, .epsg32730, .epsg32731, .epsg32732, .epsg32733, .epsg32734, .epsg32735, .epsg32736, + .epsg32737, .epsg32738, .epsg32739, .epsg32740, .epsg32741, .epsg32742, .epsg32743, .epsg32744, .epsg32745, .epsg32746, .epsg32747, .epsg32748, + .epsg32749, .epsg32750, .epsg32751, .epsg32752, .epsg32753, .epsg32754, .epsg32755, .epsg32756, .epsg32757, .epsg32758, .epsg32759, .epsg32760: + .planar } } @@ -123,6 +374,17 @@ extension Projection { GISTool.originShift case .epsg4326, .epsg32662: 180.0 + case .epsg32601, .epsg32602, .epsg32603, .epsg32604, .epsg32605, .epsg32606, .epsg32607, .epsg32608, .epsg32609, .epsg32610, .epsg32611, .epsg32612, + .epsg32613, .epsg32614, .epsg32615, .epsg32616, .epsg32617, .epsg32618, .epsg32619, .epsg32620, .epsg32621, .epsg32622, .epsg32623, .epsg32624, + .epsg32625, .epsg32626, .epsg32627, .epsg32628, .epsg32629, .epsg32630, .epsg32631, .epsg32632, .epsg32633, .epsg32634, .epsg32635, .epsg32636, + .epsg32637, .epsg32638, .epsg32639, .epsg32640, .epsg32641, .epsg32642, .epsg32643, .epsg32644, .epsg32645, .epsg32646, .epsg32647, .epsg32648, + .epsg32649, .epsg32650, .epsg32651, .epsg32652, .epsg32653, .epsg32654, .epsg32655, .epsg32656, .epsg32657, .epsg32658, .epsg32659, .epsg32660, + .epsg32701, .epsg32702, .epsg32703, .epsg32704, .epsg32705, .epsg32706, .epsg32707, .epsg32708, .epsg32709, .epsg32710, .epsg32711, .epsg32712, + .epsg32713, .epsg32714, .epsg32715, .epsg32716, .epsg32717, .epsg32718, .epsg32719, .epsg32720, .epsg32721, .epsg32722, .epsg32723, .epsg32724, + .epsg32725, .epsg32726, .epsg32727, .epsg32728, .epsg32729, .epsg32730, .epsg32731, .epsg32732, .epsg32733, .epsg32734, .epsg32735, .epsg32736, + .epsg32737, .epsg32738, .epsg32739, .epsg32740, .epsg32741, .epsg32742, .epsg32743, .epsg32744, .epsg32745, .epsg32746, .epsg32747, .epsg32748, + .epsg32749, .epsg32750, .epsg32751, .epsg32752, .epsg32753, .epsg32754, .epsg32755, .epsg32756, .epsg32757, .epsg32758, .epsg32759, .epsg32760: + nil } } @@ -138,4 +400,4 @@ extension Projection { isGeographic ? meters / 111_325.0 : meters } -} +} \ No newline at end of file diff --git a/Sources/GISTools/GeoJson/ProjectionDefinition.swift b/Sources/GISTools/GeoJson/ProjectionDefinition.swift index a4006fe..7cfbf47 100644 --- a/Sources/GISTools/GeoJson/ProjectionDefinition.swift +++ b/Sources/GISTools/GeoJson/ProjectionDefinition.swift @@ -50,32 +50,66 @@ protocol ProjectionDefinition: Sendable { /// Registry providing the ``ProjectionDefinition`` for each ``Projection``. enum ProjectionRegistry { - // MARK: - Lookup + // MARK: - Definitions - /// The definitions in match order. More specific WKT patterns must - /// come before more generic ones. - private static let allDefinitions: [any ProjectionDefinition] = [ - epsg3857Definition, - epsg3395Definition, - epsg32662Definition, - epsg4978Definition, - epsg4326Definition, - noSridDefinition, - ] + private static let noSridDefinition = NoSridDefinition() + private static let epsg3857Definition = Epsg3857Definition() + private static let epsg4326Definition = Epsg4326Definition() + private static let epsg4978Definition = Epsg4978Definition() + private static let epsg3395Definition = Epsg3395Definition() + private static let epsg32662Definition = Epsg32662Definition() + + /// All definitions keyed by projection, built eagerly so lookups are + /// constant time and every `Projection` case is covered. + /// + /// `UtmTests` sweeps every UTM SRID and asserts that lookups resolve to + /// definitions whose ``ProjectionDefinition/projection`` matches. + private static let definitions: [Projection: any ProjectionDefinition] = { + var map: [Projection: any ProjectionDefinition] = [ + .noSRID: noSridDefinition, + .epsg3857: epsg3857Definition, + .epsg4326: epsg4326Definition, + .epsg4978: epsg4978Definition, + .epsg3395: epsg3395Definition, + .epsg32662: epsg32662Definition, + ] + + for srid in 32_601 ... 32_760 { + guard let projection = Projection(srid: srid), + let utmDefinition = UtmDefinition.definition(for: projection) + else { continue } + map[projection] = utmDefinition + } + + return map + }() + + /// The definitions with WKT patterns, in match order. More specific + /// patterns must come before more generic ones (e.g. EPSG:3857's + /// "Pseudo-Mercator" before EPSG:3395's "Mercator"). UTM zones register + /// no WKT patterns. + private static var wktDefinitions: [any ProjectionDefinition] { + [ + epsg3857Definition, + epsg3395Definition, + epsg32662Definition, + epsg4978Definition, + epsg4326Definition, + ] + } + + // MARK: - Lookup /// The definition implementing the math for a projection. /// /// - Parameter projection: The projection to look up /// - Returns: The definition for the projection static func definition(for projection: Projection) -> any ProjectionDefinition { - switch projection { - case .noSRID: noSridDefinition - case .epsg3857: epsg3857Definition - case .epsg4326: epsg4326Definition - case .epsg4978: epsg4978Definition - case .epsg3395: epsg3395Definition - case .epsg32662: epsg32662Definition + guard let definition = definitions[projection] else { + assertionFailure("No projection definition registered for \(projection)") + return noSridDefinition } + return definition } /// The first definition with a WKT alternative fully contained in the string. @@ -83,22 +117,13 @@ enum ProjectionRegistry { /// - Parameter wkt: A WKT projection string /// - Returns: The matching definition, or `nil` if the string is not recognised static func definition(matchingWkt wkt: String) -> (any ProjectionDefinition)? { - allDefinitions.first { definition in + wktDefinitions.first { definition in definition.wktMatchers.contains { alternative in alternative.allSatisfy { wkt.contains($0) } } } } - // MARK: - Definitions - - private static let noSridDefinition = NoSridDefinition() - private static let epsg3857Definition = Epsg3857Definition() - private static let epsg4326Definition = Epsg4326Definition() - private static let epsg4978Definition = Epsg4978Definition() - private static let epsg3395Definition = Epsg3395Definition() - private static let epsg32662Definition = Epsg32662Definition() - } // MARK: - EPSG:4326 diff --git a/Sources/GISTools/GeoJson/UtmDefinition.swift b/Sources/GISTools/GeoJson/UtmDefinition.swift new file mode 100644 index 0000000..bedc493 --- /dev/null +++ b/Sources/GISTools/GeoJson/UtmDefinition.swift @@ -0,0 +1,210 @@ +import Foundation + +/// The hemisphere of a UTM zone. +enum UtmHemisphere: Sendable { + + /// Northern zone (EPSG:32601–32660, false northing 0). + case north + + /// Southern zone (EPSG:32701–32760, false northing 10_000_000 m). + case south + + /// The false northing applied to y values of the zone. + var falseNorthing: Double { + switch self { + case .north: 0.0 + case .south: 10_000_000.0 + } + } + +} + +/// Transverse Mercator math for a single UTM zone. +/// +/// Uses the standard Snyder transverse Mercator formulas (*Map Projections - +/// A Working Manual*, USGS PP 1395) on the WGS84 ellipsoid with the UTM +/// central scale factor of 0.9996. Coordinates are easting/northing in +/// meters with a false easting of 500_000 m and a false northing of +/// 0 m (northern zones) or 10_000_000 m (southern zones). +/// +/// The formulas lose accuracy farther than about ±3–4° from the zone's +/// central meridian, which is inherent to the transverse Mercator series +/// and matches the use expected of UTM coordinates. +struct UtmDefinition: ProjectionDefinition { + + /// The projection case of the zone (EPSG:326xx/327xx). + let projection: Projection + + /// The UTM zone number, 1–60. + let zone: Int + + /// The hemisphere the zone covers. + let hemisphere: UtmHemisphere + + /// The UTM central scale factor. + private static let k0 = 0.9996 + + /// The false easting applied to x values of every zone. + private static let falseEasting = 500_000.0 + + /// The longitude of the zone's central meridian, in radians. + var centralMeridian: Double { + Double((zone - 1) * 6 - 180 + 3) * .pi / 180.0 + } + + var kind: ProjectionKind { .planar } + + /// The EPSG-defined sector extent of the zone. + var validExtent: (minX: Double, minY: Double, maxX: Double, maxY: Double)? { + ( + minX: 100_000.0, + minY: 0.0, + maxX: 900_000.0, + maxY: 10_000_000.0 + ) + } + + /// The sector of the zone in its own coordinates (the full extent of + /// the projection: a UTM zone covers no "world" extent). + var worldBoundingBox: BoundingBox? { + guard let extent = validExtent else { return nil } + return BoundingBox( + southWest: Coordinate3D(x: extent.minX, y: extent.minY, projection: projection), + northEast: Coordinate3D(x: extent.maxX, y: extent.maxY, projection: projection)) + } + + /// UTM zones have no unique WKT fragment. + var wktMatchers: [[String]] { + [] + } + + // MARK: - Lookup + + /// Returns the UTM definition for a projection, if the projection is a + /// UTM zone SRID. + /// + /// - Parameter projection: The projection to look up + /// - Returns: The definition for the zone, or `nil` if the projection is not UTM + static func definition(for projection: Projection) -> UtmDefinition? { + let srid = projection.rawValue + + let isNorthern = srid >= 32_601 && srid <= 32_660 + let isSouthern = srid >= 32_701 && srid <= 32_760 + guard isNorthern || isSouthern else { return nil } + + let hemisphere: UtmHemisphere = isNorthern ? .north : .south + let zone = isNorthern ? srid - 32_600 : srid - 32_700 + guard zone >= 1, zone <= 60 else { return nil } + + return UtmDefinition(projection: projection, zone: zone, hemisphere: hemisphere) + } + + // MARK: - Conversions + + func forward(_ coordinate: Coordinate3D) -> Coordinate3D { + let a = GISTool.equatorialRadius + let e2 = GISTool.wgs84EccentricitySquared + let e1sq = e2 / (1.0 - e2) + + let phi = coordinate.latitude.degreesToRadians + let deltaLambda = (coordinate.longitude.degreesToRadians) - centralMeridian + + let sinPhi = sin(phi) + let cosPhi = cos(phi) + let tanPhi = tan(phi) + + let radiusSquared = 1.0 - e2 * sinPhi * sinPhi + let nRadius = a / sqrt(radiusSquared) + let tRadius = tanPhi * tanPhi + let cTerm = e1sq * cosPhi * cosPhi + let aTerm = deltaLambda * cosPhi + + let meridianArc = Self.meridianArc(phi) + + let easting = Self.falseEasting + + Self.k0 * nRadius * ( + aTerm + + (1.0 - tRadius + cTerm) * pow(aTerm, 3) / 6.0 + + (5.0 - 18.0 * tRadius + tRadius * tRadius + 72.0 * cTerm - 58.0 * e1sq) * pow(aTerm, 5) / 120.0 + ) + + let northing = Self.k0 * ( + meridianArc + + nRadius * tanPhi * ( + aTerm * aTerm / 2.0 + + (5.0 - tRadius + 9.0 * cTerm + 4.0 * cTerm * cTerm) * pow(aTerm, 4) / 24.0 + + (61.0 - 58.0 * tRadius + tRadius * tRadius + 600.0 * cTerm - 330.0 * e1sq) * pow(aTerm, 6) / 720.0 + ) + ) + hemisphere.falseNorthing + + return Coordinate3D( + x: easting, + y: northing, + z: coordinate.altitude, + m: coordinate.m, + projection: projection) + } + + func inverse(_ coordinate: Coordinate3D) -> Coordinate3D { + let a = GISTool.equatorialRadius + let e2 = GISTool.wgs84EccentricitySquared + let e1sq = e2 / (1.0 - e2) + + let x = coordinate.longitude - Self.falseEasting + let y = coordinate.latitude - hemisphere.falseNorthing + + let meridianArc = y / Self.k0 + let mu = meridianArc / (a * (1.0 - e2 / 4.0 - 3.0 * e2 * e2 / 64.0 - 5.0 * pow(e2, 3) / 256.0)) + + let e1 = (1.0 - sqrt(1.0 - e2)) / (1.0 + sqrt(1.0 - e2)) + + let phi1 = mu + + (3.0 * e1 / 2.0 - 27.0 * pow(e1, 3) / 32.0) * sin(2.0 * mu) + + (21.0 * e1 * e1 / 16.0 - 55.0 * pow(e1, 4) / 32.0) * sin(4.0 * mu) + + (151.0 * pow(e1, 3) / 96.0) * sin(6.0 * mu) + + (1097.0 * pow(e1, 4) / 512.0) * sin(8.0 * mu) + + let radiusSquared1 = 1.0 - e2 * sin(phi1) * sin(phi1) + let nRadius1 = a / sqrt(radiusSquared1) + let tRadius1 = tan(phi1) * tan(phi1) + let cTerm1 = e1sq * cos(phi1) * cos(phi1) + let meridianRadius1 = a * (1.0 - e2) / pow(radiusSquared1, 1.5) + let dTerm = x / (nRadius1 * Self.k0) + + let latitude = phi1 + - (nRadius1 * tan(phi1) / meridianRadius1) * ( + dTerm * dTerm / 2.0 + - (5.0 + 3.0 * tRadius1 + 10.0 * cTerm1 - 4.0 * cTerm1 * cTerm1 - 9.0 * e1sq) * pow(dTerm, 4) / 24.0 + + (61.0 + 90.0 * tRadius1 + 298.0 * cTerm1 + 45.0 * tRadius1 * tRadius1 - 252.0 * e1sq - 3.0 * cTerm1 * cTerm1) * pow(dTerm, 6) / 720.0 + ) + + let lambdaTerm = (dTerm + - (1.0 + 2.0 * tRadius1 + cTerm1) * pow(dTerm, 3) / 6.0 + + (5.0 - 2.0 * cTerm1 + 28.0 * tRadius1 - 3.0 * cTerm1 * cTerm1 + 8.0 * e1sq + 24.0 * tRadius1 * tRadius1) * pow(dTerm, 5) / 120.0) + / cos(phi1) + let longitude = centralMeridian + lambdaTerm + + return Coordinate3D( + latitude: latitude * 180.0 / .pi, + longitude: longitude * 180.0 / .pi, + altitude: coordinate.altitude, + m: coordinate.m) + } + + // MARK: - Snippet helpers + + /// Snyder meridian arc length (formulas 3-21) from the equator to `phi`, + /// in meters. + private static func meridianArc(_ phi: Double) -> Double { + let a = GISTool.equatorialRadius + let e2 = GISTool.wgs84EccentricitySquared + + return a * ( + (1.0 - e2 / 4.0 - 3.0 * e2 * e2 / 64.0 - 5.0 * pow(e2, 3) / 256.0) * phi + - (3.0 * e2 / 8.0 + 3.0 * e2 * e2 / 32.0 + 45.0 * pow(e2, 3) / 1024.0) * sin(2.0 * phi) + + (15.0 * e2 * e2 / 256.0 + 45.0 * pow(e2, 3) / 1024.0) * sin(4.0 * phi) + - (35.0 * pow(e2, 3) / 3072.0) * sin(6.0 * phi) + ) + } + +} diff --git a/Tests/GISToolsTests/Algorithms/UtmAlgorithmTests.swift b/Tests/GISToolsTests/Algorithms/UtmAlgorithmTests.swift new file mode 100644 index 0000000..43bd8c9 --- /dev/null +++ b/Tests/GISToolsTests/Algorithms/UtmAlgorithmTests.swift @@ -0,0 +1,194 @@ +#if canImport(CoreLocation) +import CoreLocation +#endif +import Foundation +@testable import GISTools +import Testing + +/// Algorithm round trips and smoke tests in UTM zone projections. +struct UtmAlgorithmTests { + + // MARK: - Known reference values + + // Validates forward values against independently computed references: + // the Snyder formulas were verified in Python against numerically + // integrated meridian arcs before being ported. + @Test + func knownReferenceValues() async throws { + // Zone 31N at 50N/3E (central meridian): x is exactly the false + // easting, y is the scaled meridian arc to 50 degrees (verified + // against numeric integration of the meridian radius). + let zone31 = Coordinate3D(latitude: 50.0, longitude: 3.0).projected(to: .epsg32631) + #expect(abs(zone31.longitude - 500_000.0) < 0.000001) + #expect(abs(zone31.latitude - 5_538_630.703) < 0.001) + + // Zone edges at the equator are the known ±333_978.56 meter offsets + // around the false easting (zones 1 and 60 both touch the date line). + let zoneEdgeWest = Coordinate3D(latitude: 0.0, longitude: -180.0).projected(to: .epsg32601) + #expect(abs(zoneEdgeWest.longitude - 166_021.443) < 0.001) + #expect(abs(zoneEdgeWest.latitude) < 0.000001) + + let zoneEdgeEast = Coordinate3D(latitude: 0.0, longitude: 174.0).projected(to: .epsg32660) + #expect(abs(zoneEdgeEast.longitude - 166_021.443) < 0.001) + + // Zone 19N at 41N/71W. + let zone19 = Coordinate3D(latitude: 41.0, longitude: -71.0).projected(to: .epsg32619) + #expect(abs(zone19.longitude - 331_792.115) < 0.001) + #expect(abs(zone19.latitude - 4_540_683.529) < 0.001) + + // Zone 20S at 45.5S/63W: northing is below 10_000_000 (false northing). + let zone20south = Coordinate3D(latitude: -45.5, longitude: -63.0).projected(to: .epsg32720) + #expect(abs(zone20south.longitude - 500_000.0) < 0.000001) + #expect(abs(zone20south.latitude - 4_961_503.495) < 0.001) + } + + // Validates round trips between EPSG:4326 and both hemispheres of the + // same zone: eastings are identical and northings differ by exactly the + // 10_000_000 false-northing offset. + @Test + func hemisphereSymmetry() async throws { + let equator = Coordinate3D(latitude: 0.0, longitude: -71.0) + let north = equator.projected(to: .epsg32619) + let south = equator.projected(to: .epsg32719) + + #expect(abs(north.longitude - south.longitude) < 0.0000000001) + #expect(abs(south.latitude - (north.latitude + 10_000_000.0)) < 0.000001) + + #expect(abs(south.projected(to: .epsg4326).longitude - -71.0) < 0.0000000001) + #expect(abs(south.projected(to: .epsg4326).latitude - 0.0) < 0.0000000001) + } + + // MARK: - Round trips + + // Round-trips coordinates through each point's own zone across all + // longitudes and near-polar latitudes. + @Test + func roundTrips() async throws { + let latitudes: [Double] = [0.0, 41.0, -45.5, 80.0, -80.0, 89.0, -89.0] + let longitudes: [Double] = [-180.0, -177.0, -150.0, -71.0, 0.0, 3.0, 45.0, 150.0, 174.0, 179.0] + + for longitude in longitudes { + let zone = Int((longitude + 180.0).rounded(.down) / 6.0) + 1 + for latitude in latitudes { + let srid: Int = latitude >= 0 ? 32_600 + zone : 32_700 + zone + let projection: Projection = try #require(Projection(srid: srid)) + let original = Coordinate3D(latitude: latitude, longitude: longitude) + + let projected = original.projected(to: projection) + #expect(abs(projected.longitude - 500_000.0) < 350_000.0) // inside the zone + + let back = projected.projected(to: .epsg4326) + #expect(abs(back.latitude - latitude) < 0.00001) + #expect(abs(back.longitude - longitude) < 0.00001) + } + } + } + + // MARK: - Algorithm smoke tests + + // Validates bounding boxes and containment in a UTM zone. + @Test + func boundingBoxAndContains() async throws { + let coordinates: [[Coordinate3D]] = [[ + Coordinate3D(x: 300_000.0, y: 4_000_000.0, projection: .epsg32619), + Coordinate3D(x: 400_000.0, y: 4_000_000.0, projection: .epsg32619), + Coordinate3D(x: 400_000.0, y: 4_100_000.0, projection: .epsg32619), + Coordinate3D(x: 300_000.0, y: 4_100_000.0, projection: .epsg32619), + Coordinate3D(x: 300_000.0, y: 4_000_000.0, projection: .epsg32619), + ]] + let polygon = try #require(Polygon(coordinates, calculateBoundingBox: true)) + #expect(polygon.projection == .epsg32619) + + let boundingBox = try #require(polygon.boundingBox) + #expect(boundingBox.projection == .epsg32619) + #expect(boundingBox.southWest.x == 300_000.0) + #expect(boundingBox.northEast.x == 400_000.0) + + #expect(polygon.contains(Coordinate3D(x: 350_000.0, y: 4_050_000.0, projection: .epsg32619))) + #expect(!polygon.contains(Coordinate3D(x: 450_000.0, y: 4_050_000.0, projection: .epsg32619))) + } + + // Validates Euclidean distance (planar kind) in a UTM zone. + @Test + func distance() async throws { + let origin = Coordinate3D(x: 500_000.0, y: 4_000_000.0, projection: .epsg32619) + let point = Coordinate3D(x: 500_000.0, y: 4_300_000.0, projection: .epsg32619) + + #expect(abs(origin.distance(from: point) - 300_000.0) < 0.000001) + } + + // Validates that the antimeridian machinery is a no-op for UTM zones: + // coordinates in a zone never wrap, so nothing crosses or cuts. + @Test + func antimeridianIsNoop() async throws { + #expect(Projection.epsg32660.wraparoundExtent == nil) + + // A geographic line across the date line, expressed in zone 60/1 + // coordinates: no wrap is possible within a zone. + let lineCoordinates: [Coordinate3D] = [ + Coordinate3D(x: 700_000.0, y: 0.0, projection: .epsg32660), + Coordinate3D(x: 100_000.0, y: 10_000.0, projection: .epsg32660), + ] + let line = try #require(LineString(lineCoordinates)) + + #expect(!line.crossesAntimeridian) + let cut = line.cutAtAntimeridian() + #expect(cut.features.count == 1) + } + + // Validates normalize/clamp in a zone sector. + @Test + func normalizeAndClamp() async throws { + // Zone coordinates never wrap (no wraparound extent). + let normalized = Coordinate3D(x: 950_000.0, y: 5_000_000.0, projection: .epsg32619).normalized() + #expect(normalized.longitude == 950_000.0) + + let clamped = Coordinate3D(x: 50_000.0, y: 20_000_000.0, projection: .epsg32620).clamped() + #expect(clamped.longitude == 100_000.0) + #expect(clamped.latitude == 10_000_000.0) + } + + // Validates geodesic operations via the EPSG:4326 pivot. + @Test + func destinationAndBearing() async throws { + // 100_000 meters north of a zone point. + // `destination()` uses a spherical earth model, so the grid northern + // distance differs from the elipsoidal meridian arc by a few hundred + // meters at mid latitudes; the tolerance accounts for that. + let start = Coordinate3D(x: 500_000.0, y: 4_000_000.0, projection: .epsg32619) + let destination = start.destination(distance: 100_000.0, bearing: 0.0) + #expect(destination.projection == .epsg32619) + #expect(abs(destination.longitude - 500_000.0) < 1.0) + #expect(abs(destination.latitude - 4_100_000.0) < 500.0) + #expect(destination.latitude > start.latitude) + } + + // MARK: - Coders + + // Validates a WKB round trip with an embedded UTM zone SRID. + @Test + func wkbRoundTrip() async throws { + let point = Point(Coordinate3D(x: 331_792.115, y: 4_540_683.529, projection: .epsg32619)) + + let wkb = try #require(WKBCoder.encode(geometry: point, targetProjection: .epsg32619)) + let decoded = try #require(WKBCoder.decode( + wkb: wkb, + sourceSrid: nil, + targetProjection: .epsg32619) as? Point) + #expect(decoded.projection == .epsg32619) + #expect(abs(decoded.coordinate.x - 331_792.115) < 0.000001) + #expect(abs(decoded.coordinate.y - 4_540_683.529) < 0.000001) + } + + // Validates random coordinate generation in a zone's world box. + @Test + func random() async throws { + let coordinate = BoundingBox.randomCoordinate(projection: .epsg32619) + #expect(coordinate.projection == .epsg32619) + #expect(coordinate.longitude >= 100_000.0) + #expect(coordinate.longitude <= 900_000.0) + #expect(coordinate.latitude >= 0.0) + #expect(coordinate.latitude <= 10_000_000.0) + } + +} diff --git a/Tests/GISToolsTests/GeoJson/UtmDefinitionTests.swift b/Tests/GISToolsTests/GeoJson/UtmDefinitionTests.swift new file mode 100644 index 0000000..a108f1c --- /dev/null +++ b/Tests/GISToolsTests/GeoJson/UtmDefinitionTests.swift @@ -0,0 +1,82 @@ +import Foundation +@testable import GISTools +import Testing + +/// Tests for the UTM zone SRIDs and their definitions. +struct UtmDefinitionTests { + + /// Validates that all 120 UTM SRIDs map to their zone projection and + /// resolve to a registered definition with the expected metadata. + @Test + func sridMapping() async throws { + for zone in 1 ... 60 { + let northern = try #require(Projection(srid: 32_600 + zone)) + #expect(northern.rawValue == 32_600 + zone) + let northernDefinition = UtmDefinition.definition(for: northern) + #expect(northernDefinition?.zone == zone) + #expect(northernDefinition?.hemisphere == .north) + #expect(northernDefinition?.projection == northern) + + let southern = try #require(Projection(srid: 32_700 + zone)) + let southernDefinition = UtmDefinition.definition(for: southern) + #expect(southernDefinition?.zone == zone) + #expect(southernDefinition?.hemisphere == .south) + #expect(southernDefinition?.projection == southern) + } + + // Non-UTM SRIDs have no UTM definition. + #expect(UtmDefinition.definition(for: .epsg4326) == nil) + #expect(UtmDefinition.definition(for: .epsg3857) == nil) + #expect(UtmDefinition.definition(for: .noSRID) == nil) + } + + /// Validates that every projection case resolves to a registered + /// definition (registry coverage). + @Test + func registryCoverage() async throws { + for srid in [0, 3857, 4326, 4978, 3395, 32662] { + let projection = try #require(Projection(srid: srid)) + #expect(ProjectionRegistry.definition(for: projection).projection == projection) + } + for srid in 32_601 ... 32_660 { + let projection = try #require(Projection(srid: srid)) + #expect(ProjectionRegistry.definition(for: projection).projection == projection) + } + for srid in 32_701 ... 32_760 { + let projection = try #require(Projection(srid: srid)) + #expect(ProjectionRegistry.definition(for: projection).projection == projection) + } + } + + /// Validates zone metadata. + @Test + func zoneMetadata() async throws { + for zone in 1 ... 60 { + let projection = try #require(Projection(srid: 32_600 + zone)) + let definition = try #require(UtmDefinition.definition(for: projection)) + + // Central meridian: zone 1 is at -177 degrees, zones step by 6 degrees. + let expectedCentralLongitude = Double((zone - 1) * 6 - 180 + 3) + #expect(abs(definition.centralMeridian * 180.0 / .pi - expectedCentralLongitude) < 0.0000000001) + } + + // Sector extent per hemisphere convention (false northing caps y). + let northern = try #require(UtmDefinition.definition(for: .epsg32619)) + #expect(northern.hemisphere.falseNorthing == 0.0) + #expect(northern.validExtent?.minY == 0.0) + #expect(northern.validExtent?.maxY == 10_000_000.0) + #expect(northern.validExtent?.minX == 100_000.0) + #expect(northern.validExtent?.maxX == 900_000.0) + + // The world bounding box covers the zone sector; UTM zones do not wrap. + let worldBox = try #require(northern.worldBoundingBox) + #expect(worldBox.projection == .epsg32619) + #expect(worldBox.southWest.x == 100_000.0) + #expect(worldBox.northEast.x == 900_000.0) + + // Description strings. + #expect(Projection.epsg32619.description == "EPSG:32619") + #expect(Projection.epsg32719.description == "EPSG:32719") + } + +}