diff --git a/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/MercatorAlgebra.kt b/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/MercatorAlgebra.kt new file mode 100644 index 000000000..32b4362ae --- /dev/null +++ b/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/MercatorAlgebra.kt @@ -0,0 +1,97 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.kmath.geometry + +import space.kscience.kmath.operations.Algebra +import kotlin.math.* + +/** + * Geodetic coordinated + */ +public class GeodeticCoordinates private constructor(public val latitude: Double, public val longitude: Double) { + + override fun equals(other: Any?): Boolean { + if (this === other) return true + if (other == null || this::class != other::class) return false + + other as GeodeticCoordinates + + if (latitude != other.latitude) return false + if (longitude != other.longitude) return false + + return true + } + + override fun hashCode(): Int { + var result = latitude.hashCode() + result = 31 * result + longitude.hashCode() + return result + } + + override fun toString(): String { + return "GeodeticCoordinates(latitude=$latitude, longitude=$longitude)" + } + + + public companion object { + public fun ofRadians(latitude: Double, longitude: Double): GeodeticCoordinates { + require(longitude in (-PI)..(PI)) { "Longitude $longitude is not in (-PI)..(PI)" } + return GeodeticCoordinates(latitude, longitude.rem(PI / 2)) + } + + public fun ofDegrees(latitude: Double, longitude: Double): GeodeticCoordinates { + require(latitude in (-90.0)..(90.0)) { "Latitude $latitude is not in -90..90" } + return GeodeticCoordinates(latitude * PI / 180, (longitude.rem(180) * PI / 180)) + } + } +} + +public data class MercatorCoordinates(val x: Double, val y: Double) + +/** + * @param baseLongitude the longitude offset in radians + * @param radius the average radius of the Earth + * @param correctedRadius optional radius correction to account for ellipsoid model + */ +public open class MercatorAlgebra( + public val baseLongitude: Double = 0.0, + private val radius: Double = DEFAULT_EARTH_RADIUS, + private val correctedRadius: ((GeodeticCoordinates) -> Double)? = null, +) : Algebra { + + public fun MercatorCoordinates.toGeodetic(): GeodeticCoordinates { + val res = GeodeticCoordinates.ofRadians( + atan(sinh(y / radius)), + baseLongitude + x / radius, + ) + return if (correctedRadius != null) { + val r = correctedRadius.invoke(res) + GeodeticCoordinates.ofRadians( + atan(sinh(y / r)), + baseLongitude + x / r, + ) + } else { + res + } + } + + /** + * https://en.wikipedia.org/wiki/Web_Mercator_projection#Formulas + */ + public fun GeodeticCoordinates.toMercator(): MercatorCoordinates { + require(abs(latitude) <= MAXIMUM_LATITUDE) { "Latitude exceeds the maximum latitude for mercator coordinates" } + val r = correctedRadius?.invoke(this) ?: radius + return MercatorCoordinates( + x = r * (longitude - baseLongitude), + y = r * ln(tan(PI / 4 + latitude / 2)) + ) + } + + public companion object : MercatorAlgebra() { + public const val MAXIMUM_LATITUDE: Double = 85.05113 + public const val DEFAULT_EARTH_RADIUS: Double = 6378137.0 + } +} \ No newline at end of file diff --git a/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/WebMercatorAlgebra.kt b/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/WebMercatorAlgebra.kt index 959146072..e93649eef 100644 --- a/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/WebMercatorAlgebra.kt +++ b/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/WebMercatorAlgebra.kt @@ -8,51 +8,13 @@ package space.kscience.kmath.geometry import space.kscience.kmath.operations.Algebra import kotlin.math.* -public class GeodeticCoordinates private constructor(public val latitude: Double, public val longitude: Double) { +public data class TileWebMercatorCoordinates(val zoom: Double, val x: Double, val y: Double) - override fun equals(other: Any?): Boolean { - if (this === other) return true - if (other == null || this::class != other::class) return false - - other as GeodeticCoordinates - - if (latitude != other.latitude) return false - if (longitude != other.longitude) return false - - return true - } - - override fun hashCode(): Int { - var result = latitude.hashCode() - result = 31 * result + longitude.hashCode() - return result - } - - override fun toString(): String { - return "GeodeticCoordinates(latitude=$latitude, longitude=$longitude)" - } - - - public companion object { - public fun ofRadians(latitude: Double, longitude: Double): GeodeticCoordinates { - require(longitude in (-PI)..(PI)) { "Longitude $longitude is not in (-PI)..(PI)" } - return GeodeticCoordinates(latitude, longitude.rem(PI / 2)) - } - - public fun ofDegrees(latitude: Double, longitude: Double): GeodeticCoordinates { - require(latitude in (-90.0)..(90.0)) { "Latitude $latitude is not in -90..90" } - return GeodeticCoordinates(latitude * PI / 180, (longitude.rem(180) * PI / 180)) - } - } -} - -public data class WebMercatorCoordinates(val zoom: Double, val x: Double, val y: Double) - -public object WebMercatorAlgebra : Algebra { +public object WebMercatorAlgebra : Algebra { private fun scaleFactor(zoom: Double) = 256.0 / 2 / PI * 2.0.pow(zoom) - public fun WebMercatorCoordinates.toGeodetic(): GeodeticCoordinates { + public fun TileWebMercatorCoordinates.toGeodetic(): GeodeticCoordinates { val scaleFactor = scaleFactor(zoom) val longitude = x / scaleFactor - PI val latitude = (atan(exp(PI - y / scaleFactor)) - PI / 4) * 2 @@ -62,11 +24,11 @@ public object WebMercatorAlgebra : Algebra { /** * https://en.wikipedia.org/wiki/Web_Mercator_projection#Formulas */ - public fun GeodeticCoordinates.toMercator(zoom: Double): WebMercatorCoordinates { - require(abs(latitude) <= 2 * atan(E.pow(PI)) - PI / 2) { "Latitude exceeds the maximum latitude for mercator coordinates" } + public fun GeodeticCoordinates.toMercator(zoom: Double): TileWebMercatorCoordinates { + require(abs(latitude) <= MercatorAlgebra.MAXIMUM_LATITUDE) { "Latitude exceeds the maximum latitude for mercator coordinates" } val scaleFactor = scaleFactor(zoom) - return WebMercatorCoordinates( + return TileWebMercatorCoordinates( zoom = zoom, x = scaleFactor * (longitude + PI), y = scaleFactor * (PI - ln(tan(PI / 4 + latitude / 2))) @@ -81,7 +43,7 @@ public object WebMercatorAlgebra : Algebra { base: GeodeticCoordinates, target: GeodeticCoordinates, zoom: Double? = null, - ): WebMercatorCoordinates { + ): TileWebMercatorCoordinates { val xOffsetUnscaled = target.longitude - base.longitude val yOffsetUnscaled = ln( tan(PI / 4 + target.latitude / 2) / tan(PI / 4 + base.latitude / 2) @@ -89,7 +51,7 @@ public object WebMercatorAlgebra : Algebra { val computedZoom = zoom ?: floor(log2(PI / max(abs(xOffsetUnscaled), abs(yOffsetUnscaled)))) val scaleFactor = scaleFactor(computedZoom) - return WebMercatorCoordinates( + return TileWebMercatorCoordinates( computedZoom, x = scaleFactor * xOffsetUnscaled, y = scaleFactor * yOffsetUnscaled diff --git a/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/MercatorTest.kt b/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/MercatorTest.kt index 81b74ff6e..4f6060de2 100644 --- a/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/MercatorTest.kt +++ b/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/MercatorTest.kt @@ -16,6 +16,16 @@ class MercatorTest { assertEquals(GeodeticCoordinates.ofDegrees(45.0, 20.0), GeodeticCoordinates.ofDegrees(45.0, 200.0)) } + @Test + fun mercatorConversion() = with(MercatorAlgebra) { + val mskCoordinates = GeodeticCoordinates.ofDegrees(55.7558, 37.6173) + val m = mskCoordinates.toMercator() + val r = m.toGeodetic() + + assertEquals(mskCoordinates.longitude, r.longitude,1e-4) + assertEquals(mskCoordinates.latitude, r.latitude,1e-4) + } + @Test fun webMercatorConversion() = with(WebMercatorAlgebra) { val mskCoordinates = GeodeticCoordinates.ofDegrees(55.7558, 37.6173)