From 8303092cc36af2b946b5fb60ceeb5533baeb4614 Mon Sep 17 00:00:00 2001 From: Thomas Rasch Date: Thu, 25 Jul 2024 14:40:58 +0200 Subject: [PATCH] Added GISTool.convertToDegrees(fromMeters:atLatitude:) (#59) Also added GISTool.convertToCoordinate(fromPixelX:pixelY:atZoom:tileSideLength:projection:) and GISTool.metersPerPixel(atZoom:latitude:tileSideLength:) (they were in MapTile before) --- Sources/GISTools/Algorithms/Conversions.swift | 64 ++++++++++++++++++- Sources/GISTools/GeoJson/BoundingBox.swift | 15 ++--- Sources/GISTools/Other/MapTile.swift | 37 +++++------ .../Algorithms/ConversionTests.swift | 40 ++++++++++++ Tests/GISToolsTests/Other/MapTileTests.swift | 27 +------- 5 files changed, 127 insertions(+), 56 deletions(-) create mode 100644 Tests/GISToolsTests/Algorithms/ConversionTests.swift diff --git a/Sources/GISTools/Algorithms/Conversions.swift b/Sources/GISTools/Algorithms/Conversions.swift index 7af20fc..10c75b6 100644 --- a/Sources/GISTools/Algorithms/Conversions.swift +++ b/Sources/GISTools/Algorithms/Conversions.swift @@ -5,6 +5,8 @@ import Foundation // Ported from https://github.com/Turfjs/turf/tree/master/packages/turf-helpers +// MARK: Measurements + extension GISTool { /// Unit of measurement. @@ -12,6 +14,7 @@ extension GISTool { case acres case centimeters case centimetres + /// Latitude case degrees case feet case inches @@ -31,7 +34,7 @@ extension GISTool { public static func factor(for unit: Unit) -> Double? { switch unit { case .centimeters, .centimetres: return earthRadius * 100.0 - case .degrees: return earthRadius / 111_325.0 + case .degrees: return earthRadius / (GISTool.earthCircumference / 360.0) case .feet: return earthRadius * 3.28084 case .inches: return earthRadius * 39.370 case .kilometers, .kilometres: return earthRadius / 1000.0 @@ -98,3 +101,62 @@ extension GISTool { } } + +// MARK: - Pixels + +extension GISTool { + + /// Converts pixel coordinates in a given zoom level to a coordinate. + public static func convertToCoordinate( + fromPixelX pixelX: Double, + pixelY: Double, + atZoom zoom: Int, + tileSideLength: Double = GISTool.tileSideLength, + projection: Projection = .epsg4326) + -> Coordinate3D + { + let resolution = metersPerPixel(atZoom: zoom, tileSideLength: tileSideLength) + + let coordinateXY = Coordinate3D( + x: pixelX * resolution - GISTool.originShift, + y: pixelY * resolution - GISTool.originShift) + + if projection == .epsg4326 { + return coordinateXY.projected(to: projection) + } + + return coordinateXY + } + + /// Resolution (meters/pixel) for a given zoom level (measured at `latitude`, defaults to the equator). + public static func metersPerPixel( + atZoom zoom: Int, + latitude: CLLocationDegrees = 0.0, // equator + tileSideLength: Double = GISTool.tileSideLength) + -> Double + { + (cos(latitude * Double.pi / 180.0) * 2.0 * Double.pi * GISTool.equatorialRadius / tileSideLength) / pow(2.0, Double(zoom)) + } + +} + +// MARK: - Meters/latitude + +extension GISTool { + + public static func convertToDegrees( + fromMeters meters: Double, + atLatitude latitude: CLLocationDegrees) + -> (latitudeDegrees: CLLocationDegrees, longitudeDegrees: CLLocationDegrees) + { + // Length of one minute at this latitude + let oneDegreeLatitudeDistance: CLLocationDistance = GISTool.earthCircumference / 360.0 // ~111 km + let oneDegreeLongitudeDistance: CLLocationDistance = cos(latitude * Double.pi / 180.0) * oneDegreeLatitudeDistance + + let longitudeDistance: Double = (meters / oneDegreeLongitudeDistance) + let latitudeDistance: Double = (meters / oneDegreeLatitudeDistance) + + return (latitudeDistance, longitudeDistance) + } + +} diff --git a/Sources/GISTools/GeoJson/BoundingBox.swift b/Sources/GISTools/GeoJson/BoundingBox.swift index 0843943..92855ba 100644 --- a/Sources/GISTools/GeoJson/BoundingBox.swift +++ b/Sources/GISTools/GeoJson/BoundingBox.swift @@ -77,17 +77,12 @@ public struct BoundingBox: northEast.longitude += padding case .epsg4326: - // Length of one minute at this latitude - let oneDegreeLatitudeDistance: CLLocationDistance = GISTool.earthCircumference / 360.0 // ~111 km - let oneDegreeLongitudeDistance: CLLocationDistance = cos(southWest.latitude * Double.pi / 180.0) * oneDegreeLatitudeDistance + let latLongDegrees = GISTool.convertToDegrees(fromMeters: padding, atLatitude: southWest.latitude) - let longitudeDistance: Double = (padding / oneDegreeLongitudeDistance) - let latitudeDistance: Double = (padding / oneDegreeLatitudeDistance) - - southWest.latitude -= latitudeDistance - northEast.latitude += latitudeDistance - southWest.longitude -= longitudeDistance - northEast.longitude += longitudeDistance + southWest.latitude -= latLongDegrees.latitudeDegrees + northEast.latitude += latLongDegrees.latitudeDegrees + southWest.longitude -= latLongDegrees.longitudeDegrees + northEast.longitude += latLongDegrees.longitudeDegrees case .noSRID: break // Don't know what to do -> ignore diff --git a/Sources/GISTools/Other/MapTile.swift b/Sources/GISTools/Other/MapTile.swift index 9bdff78..ddfee69 100644 --- a/Sources/GISTools/Other/MapTile.swift +++ b/Sources/GISTools/Other/MapTile.swift @@ -70,8 +70,8 @@ public struct MapTile: CustomStringConvertible, Sendable { let pixelX: Double = (Double(x) + 0.5) * GISTool.tileSideLength let pixelY: Double = (Double(y) + 0.5) * GISTool.tileSideLength - return MapTile.pixelCoordinate( - pixelX: pixelX, + return GISTool.convertToCoordinate( + fromPixelX: pixelX, pixelY: pixelY, atZoom: z, tileSideLength: GISTool.tileSideLength, @@ -89,14 +89,14 @@ public struct MapTile: CustomStringConvertible, Sendable { // Flip y let y = (1 << z) - 1 - y - let southWest = MapTile.pixelCoordinate( - pixelX: Double(x) * GISTool.tileSideLength, + let southWest = GISTool.convertToCoordinate( + fromPixelX: Double(x) * GISTool.tileSideLength, pixelY: Double(y) * GISTool.tileSideLength, atZoom: z, tileSideLength: GISTool.tileSideLength, projection: projection) - let northEast = MapTile.pixelCoordinate( - pixelX: Double(x + 1) * GISTool.tileSideLength, + let northEast = GISTool.convertToCoordinate( + fromPixelX: Double(x + 1) * GISTool.tileSideLength, pixelY: Double(y + 1) * GISTool.tileSideLength, atZoom: z, tileSideLength: GISTool.tileSideLength, @@ -161,7 +161,8 @@ public struct MapTile: CustomStringConvertible, Sendable { // MARK: - Conversion pixel to meters - /// Converts pixel coordinates in a given zoom level to EPSG:3857. + /// Converts pixel coordinates in a given zoom level to a coordinate. + @available(*, deprecated, renamed: "GISTool.convertToCoordinate", message: "This method has been moved to the GISTool namespace") public static func pixelCoordinate( pixelX: Double, pixelY: Double, @@ -170,34 +171,30 @@ public struct MapTile: CustomStringConvertible, Sendable { projection: Projection = .epsg4326) -> Coordinate3D { - let resolution = metersPerPixel(at: zoom, tileSideLength: tileSideLength) - - let coordinateXY = Coordinate3D( - x: pixelX * resolution - GISTool.originShift, - y: pixelY * resolution - GISTool.originShift) - - if projection == .epsg4326 { - return coordinateXY.projected(to: projection) - } - - return coordinateXY + GISTool.convertToCoordinate( + fromPixelX: pixelX, + pixelY: pixelY, + atZoom: zoom, + tileSideLength: tileSideLength, + projection: projection) } // MARK: - Meters per pixel /// Resolution (meters/pixel) for a given zoom level (measured at `latitude`, defaults to the equator). + @available(*, deprecated, renamed: "GISTool.metersPerPixel", message: "This method has been moved to the GISTool namespace") public static func metersPerPixel( at zoom: Int, latitude: Double = 0.0, // equator tileSideLength: Double = GISTool.tileSideLength) -> Double { - (cos(latitude * Double.pi / 180.0) * 2.0 * Double.pi * GISTool.equatorialRadius / tileSideLength) / pow(2.0, Double(zoom)) + GISTool.metersPerPixel(atZoom: zoom, latitude: latitude, tileSideLength: tileSideLength) } /// Resolution (meters/pixel) for a given zoom level measured at the tile center. public var metersPerPixel: Double { - MapTile.metersPerPixel(at: z, latitude: centerCoordinate().latitude) + GISTool.metersPerPixel(atZoom: z, latitude: centerCoordinate().latitude) } // MARK: - Private diff --git a/Tests/GISToolsTests/Algorithms/ConversionTests.swift b/Tests/GISToolsTests/Algorithms/ConversionTests.swift new file mode 100644 index 0000000..9ce12f8 --- /dev/null +++ b/Tests/GISToolsTests/Algorithms/ConversionTests.swift @@ -0,0 +1,40 @@ +@testable import GISTools +import XCTest + +final class ConversionTests: XCTestCase { + + func testMetersPerPixelAtEquator() { + var mppAtZoomLevels: [Double] = Array(repeating: 0.0, count: 21) + mppAtZoomLevels[0] = 156_543.03392804096 + + for zoom in 1...20 { + mppAtZoomLevels[zoom] = mppAtZoomLevels[zoom - 1] / 2.0 + + XCTAssertEqual(GISTool.metersPerPixel(atZoom: zoom), + mppAtZoomLevels[zoom], + accuracy: 0.00001) + } + } + + func testMetersPerPixelAt45() { + var mppAtZoomLevels: [Double] = Array(repeating: 0.0, count: 21) + mppAtZoomLevels[0] = 110_692.6408380335 + + for zoom in 1...20 { + mppAtZoomLevels[zoom] = mppAtZoomLevels[zoom - 1] / 2.0 + + XCTAssertEqual(GISTool.metersPerPixel(atZoom: zoom, latitude: 45.0), + mppAtZoomLevels[zoom], + accuracy: 0.00001) + } + } + + func testMetersAtLatitude() throws { + let meters = 10000.0 + let degreesLatitude1 = try XCTUnwrap(GISTool.convert(length: meters, from: .meters, to: .degrees)) + let degreesLatitude2 = GISTool.convertToDegrees(fromMeters: meters, atLatitude: 0.0).latitudeDegrees + + XCTAssertEqual(degreesLatitude1, degreesLatitude2, accuracy: 0.00000001) + } + +} diff --git a/Tests/GISToolsTests/Other/MapTileTests.swift b/Tests/GISToolsTests/Other/MapTileTests.swift index c465d4b..f878784 100644 --- a/Tests/GISToolsTests/Other/MapTileTests.swift +++ b/Tests/GISToolsTests/Other/MapTileTests.swift @@ -161,33 +161,10 @@ final class MapTileTests: XCTestCase { } func testMetersPerPixelAtEquator() { - var mppAtZoomLevels: [Double] = Array(repeating: 0.0, count: 21) - mppAtZoomLevels[0] = 156_543.03392804096 - - for zoom in 1...20 { - mppAtZoomLevels[zoom] = mppAtZoomLevels[zoom - 1] / 2.0 - - XCTAssertEqual(MapTile.metersPerPixel(at: zoom), - mppAtZoomLevels[zoom], - accuracy: 0.00001) - } - - // Alternative let worldTile = MapTile(x: 0, y: 0, z: 0) - XCTAssertEqual(worldTile.metersPerPixel, mppAtZoomLevels[0], accuracy: 0.00001) - } - - func testMetersPerPixelAt45() { - var mppAtZoomLevels: [Double] = Array(repeating: 0.0, count: 21) - mppAtZoomLevels[0] = 110_692.6408380335 + let mppAtZoom0 = 156_543.03392804096 - for zoom in 1...20 { - mppAtZoomLevels[zoom] = mppAtZoomLevels[zoom - 1] / 2.0 - - XCTAssertEqual(MapTile.metersPerPixel(at: zoom, latitude: 45.0), - mppAtZoomLevels[zoom], - accuracy: 0.00001) - } + XCTAssertEqual(worldTile.metersPerPixel, mppAtZoom0, accuracy: 0.00001) } func testQquadkey() {