From 081304f11007add35ca0e04fb6ad0ea3c023d4d4 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Thu, 16 Jun 2022 20:46:43 +0300
Subject: [PATCH] Web mercator conversion

---
 .../kmath/geometry/WebMercatorAlgebra.kt      | 67 +++++++++++++++++++
 .../kscience/kmath/geometry/MercatorTest.kt   | 24 +++++++
 2 files changed, 91 insertions(+)
 create mode 100644 kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/WebMercatorAlgebra.kt
 create mode 100644 kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/MercatorTest.kt

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
new file mode 100644
index 000000000..f4ea339b8
--- /dev/null
+++ b/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/WebMercatorAlgebra.kt
@@ -0,0 +1,67 @@
+/*
+ * 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.*
+
+public class GeodeticCoordinates private constructor(public val latitude: Double, public val longitude: Double) {
+    init {
+        require(longitude in (-PI)..(PI)) { "Longitude $longitude is not in (-PI)..(PI)" }
+        require(latitude in (-PI / 2)..(PI / 2)) { "Latitude $latitude is not in (-PI/2)..(PI/2)" }
+    }
+
+    public companion object {
+        public fun ofRadians(longitude: Double, latitude: Double): GeodeticCoordinates =
+            GeodeticCoordinates(latitude, longitude)
+
+        public fun ofDegrees(longitude: Double, latitude: Double): GeodeticCoordinates =
+            GeodeticCoordinates(latitude * PI / 180, longitude * PI / 180)
+    }
+}
+
+public data class WebMercatorCoordinates(val zoom: Double, val x: Double, val y: Double)
+
+public object WebMercatorAlgebra : Algebra<WebMercatorCoordinates> {
+    fun WebMercatorCoordinates.toGeodetic(): GeodeticCoordinates = TODO()
+
+    /**
+     * 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" }
+
+        val scaleFactor = 256.0 / 2 / PI * 2.0.pow(zoom)
+        return WebMercatorCoordinates(
+            zoom = zoom,
+            x = scaleFactor * (longitude + PI),
+            y = scaleFactor * (PI - ln(tan(PI / 4 + latitude / 2)))
+        )
+    }
+
+    /**
+     * Compute and offset of [target] coordinate relative to [base] coordinate. If [zoom] is null, then optimal zoom
+     * will be computed to put the resulting x and y coordinates between -127.0 and 128.0
+     */
+    public fun computeOffset(
+        base: GeodeticCoordinates,
+        target: GeodeticCoordinates,
+        zoom: Double? = null,
+    ): WebMercatorCoordinates {
+        val xOffsetUnscaled = target.longitude - base.longitude
+        val yOffsetUnscaled = ln(
+            tan(PI / 4 + target.latitude / 2) / tan(PI / 4 + base.latitude / 2)
+        )
+
+        val computedZoom = zoom ?: floor( log2(PI / max(abs(xOffsetUnscaled), abs(yOffsetUnscaled))))
+        val scaleFactor = 256.0 / 2 / PI * 2.0.pow(computedZoom)
+        return WebMercatorCoordinates(
+            computedZoom,
+            x = scaleFactor * xOffsetUnscaled,
+            y = scaleFactor * yOffsetUnscaled
+        )
+    }
+}
\ No newline at end of file
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
new file mode 100644
index 000000000..5514e5521
--- /dev/null
+++ b/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/MercatorTest.kt
@@ -0,0 +1,24 @@
+/*
+ * 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 kotlin.test.Test
+import kotlin.test.assertTrue
+
+class MercatorTest {
+
+    @Test
+    fun mercatorOffset(){
+        val mskCoordinates = GeodeticCoordinates.ofDegrees(55.7558, 37.6173)
+        val spbCoordinates = GeodeticCoordinates.ofDegrees(59.9311, 30.3609)
+
+        val offset =  WebMercatorAlgebra.computeOffset(mskCoordinates, spbCoordinates)
+
+        assertTrue { offset.x in -127.0..128.0 }
+        assertTrue { offset.y in -127.0..128.0 }
+        assertTrue { offset.zoom > 1.0 }
+    }
+}
\ No newline at end of file