Regular mercator conversion

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<TileWebMercatorCoordinates> {
+
+    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
+    }
+}

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 {
+public object WebMercatorAlgebra : Algebra<TileWebMercatorCoordinates> {
-        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<WebMercatorCoordinates> {
 
     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<WebMercatorCoordinates> {
     /**
      * https://en.wikipedia.org/wiki/Web_Mercator_projection#Formulas
      */
-    public fun GeodeticCoordinates.toMercator(zoom: Double): WebMercatorCoordinates {
+    public fun GeodeticCoordinates.toMercator(zoom: Double): TileWebMercatorCoordinates {
-        require(abs(latitude) <= 2 * atan(E.pow(PI)) - PI / 2) { "Latitude exceeds the maximum latitude for mercator coordinates" }
+        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<WebMercatorCoordinates> {
         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<WebMercatorCoordinates> {
 
         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

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)