Skip to content

Commit

Permalink
Added GISTool.convertToDegrees(fromMeters:atLatitude:)
Browse files Browse the repository at this point in the history
  • Loading branch information
trasch committed Jul 25, 2024
1 parent 004d356 commit e432f64
Show file tree
Hide file tree
Showing 5 changed files with 118 additions and 52 deletions.
64 changes: 63 additions & 1 deletion Sources/GISTools/Algorithms/Conversions.swift
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,16 @@ import Foundation

// Ported from https://github.com/Turfjs/turf/tree/master/packages/turf-helpers

// MARK: Measurements

extension GISTool {

/// Unit of measurement.
public enum Unit {
case acres
case centimeters
case centimetres
/// Latitude
case degrees
case feet
case inches
Expand All @@ -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
Expand Down Expand Up @@ -98,3 +101,62 @@ extension GISTool {
}

}

// MARK: - Pixels

extension GISTool {

/// Converts pixel coordinates in a given zoom level to EPSG:3857.
public static func pixelCoordinate(
pixelX: Double,
pixelY: Double,
atZoom zoom: Int,
tileSideLength: Double = GISTool.tileSideLength,
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
}

/// Resolution (meters/pixel) for a given zoom level (measured at `latitude`, defaults to the equator).
public static func metersPerPixel(
at 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)
}

}
15 changes: 5 additions & 10 deletions Sources/GISTools/GeoJson/BoundingBox.swift
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
24 changes: 8 additions & 16 deletions Sources/GISTools/Other/MapTile.swift
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ 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(
return GISTool.pixelCoordinate(
pixelX: pixelX,
pixelY: pixelY,
atZoom: z,
Expand All @@ -89,13 +89,13 @@ public struct MapTile: CustomStringConvertible, Sendable {
// Flip y
let y = (1 << z) - 1 - y

let southWest = MapTile.pixelCoordinate(
let southWest = GISTool.pixelCoordinate(
pixelX: Double(x) * GISTool.tileSideLength,
pixelY: Double(y) * GISTool.tileSideLength,
atZoom: z,
tileSideLength: GISTool.tileSideLength,
projection: projection)
let northEast = MapTile.pixelCoordinate(
let northEast = GISTool.pixelCoordinate(
pixelX: Double(x + 1) * GISTool.tileSideLength,
pixelY: Double(y + 1) * GISTool.tileSideLength,
atZoom: z,
Expand Down Expand Up @@ -162,6 +162,7 @@ public struct MapTile: CustomStringConvertible, Sendable {
// MARK: - Conversion pixel to meters

/// Converts pixel coordinates in a given zoom level to EPSG:3857.
@available(*, deprecated, renamed: "GISTool.pixelCoordinate", message: "This method has been moved to the GISTool namespace")
public static func pixelCoordinate(
pixelX: Double,
pixelY: Double,
Expand All @@ -170,34 +171,25 @@ 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.pixelCoordinate(pixelX: 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(at: 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(at: z, latitude: centerCoordinate().latitude)
}

// MARK: - Private
Expand Down
40 changes: 40 additions & 0 deletions Tests/GISToolsTests/Algorithms/ConversionTests.swift
Original file line number Diff line number Diff line change
@@ -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(at: 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(at: 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)
}

}
27 changes: 2 additions & 25 deletions Tests/GISToolsTests/Other/MapTileTests.swift
Original file line number Diff line number Diff line change
Expand Up @@ -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() {
Expand Down

0 comments on commit e432f64

Please sign in to comment.