Add ellipsoid mercator backward transformation

maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/Distance.kt

@@ -1,7 +1,9 @@
 package center.sciprog.maps.coordinates
 
+import kotlinx.serialization.Serializable
 import kotlin.jvm.JvmInline
 
+@Serializable
 @JvmInline
 public value class Distance internal constructor(public val kilometers: Double) : Comparable<Distance> {
     override fun compareTo(other: Distance): Int = this.kilometers.compareTo(other.kilometers)

maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GeoEllipsoid.kt

@@ -1,11 +1,12 @@
 package center.sciprog.maps.coordinates
 
+import kotlinx.serialization.Serializable
 import space.kscience.kmath.geometry.Angle
 import space.kscience.kmath.geometry.tan
 import kotlin.math.pow
 import kotlin.math.sqrt
 
-
+@Serializable
 public class GeoEllipsoid(public val equatorRadius: Distance, public val polarRadius: Distance) {
 
     /**
@@ -18,8 +19,13 @@ public class GeoEllipsoid(public val equatorRadius: Distance, public val polarRa
      */
     public val inverseF: Double = equatorRadius.kilometers / (equatorRadius.kilometers - polarRadius.kilometers)
 
+    /**
+     *  Eccentricity squared
+     */
     public val eSquared: Double = 2 * f - f * f
 
+    public val eccentricity: Double = sqrt(eSquared)
+
     public companion object {
 
         public val WGS84: GeoEllipsoid = GeoEllipsoid(

maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GeodeticMapCoordinates.kt

@@ -1,5 +1,6 @@
 package center.sciprog.maps.coordinates
 
+import kotlinx.serialization.Serializable
 import space.kscience.kmath.geometry.Angle
 import space.kscience.kmath.geometry.degrees
 import space.kscience.kmath.geometry.normalized
@@ -10,6 +11,7 @@ import space.kscience.kmath.geometry.radians
  *
  * @param elevation is optional
  */
+@Serializable
 public class GeodeticMapCoordinates(
     public val latitude: Angle,
     public val longitude: Angle,

maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcPose.kt

@@ -1,11 +1,13 @@
 package center.sciprog.maps.coordinates
 
+import kotlinx.serialization.Serializable
 import space.kscience.kmath.geometry.Angle
 import space.kscience.kmath.geometry.normalized
 
 /**
  * A coordinate-bearing pair
  */
+@Serializable
 public data class GmcPose(val coordinates: GeodeticMapCoordinates, val bearing: Angle) {
     val latitude: Angle get() = coordinates.latitude
     val longitude: Angle get() = coordinates.longitude

maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/MercatorProjection.kt

@@ -5,6 +5,7 @@
 
 package center.sciprog.maps.coordinates
 
+import kotlinx.serialization.Serializable
 import space.kscience.kmath.geometry.*
 import kotlin.math.*
 
@@ -27,25 +28,42 @@ public interface MapProjection<T : Any> {
  * @param baseLongitude the longitude offset in radians
  * @param ellipsoid - a [GeoEllipsoid] to be used for conversion
  */
+@Serializable
 public open class MercatorProjection(
     public val baseLongitude: Angle = Angle.zero,
     public val ellipsoid: GeoEllipsoid = GeoEllipsoid.sphere,
 ) : MapProjection<ProjectionCoordinates> {
 
-    override fun toGeodetic(pc: ProjectionCoordinates): GeodeticMapCoordinates {
-        val res = GeodeticMapCoordinates.ofRadians(
-            atan(sinh(pc.y / ellipsoid.equatorRadius)),
-            baseLongitude.radians + (pc.x / ellipsoid.equatorRadius),
-        )
 
+    /**
+     * Taken from https://github.com/geotools/geotools/blob/main/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/Mercator.java#L164
+     */
+    private fun cphi2(ts: Double): Double {
+        val eccnth: Double = 0.5 * ellipsoid.eccentricity
+        var phi: Double = PI / 2 - 2.0 * atan(ts)
+        for (i in 0 until 15) {
+            val con: Double = ellipsoid.eccentricity * sin(phi)
+            val dphi: Double = PI / 2 - 2.0 * atan(ts * ((1 - con) / (1 + con)).pow(eccnth)) - phi
+            phi += dphi
+            if (abs(dphi) <= 1e-10) {
+                return phi
+            }
+        }
+        error("Inverse mercator projection transformation failed to converge")
+    }
+
+
+    override fun toGeodetic(pc: ProjectionCoordinates): GeodeticMapCoordinates {
         return if (ellipsoid === GeoEllipsoid.sphere) {
-            res
+            GeodeticMapCoordinates.ofRadians(
+                atan(sinh(pc.y / ellipsoid.equatorRadius)),
+                baseLongitude.radians + (pc.x / ellipsoid.equatorRadius),
+            )
         } else {
-            TODO("Elliptical mercator projection not implemented")
+            GeodeticMapCoordinates.ofRadians(
-//            GeodeticMapCoordinates(
+                cphi2(exp(-(pc.y / ellipsoid.equatorRadius))),
-//                atan(sinh(pc.y / ellipsoid.polarRadius)).radians,
+                baseLongitude.radians + (pc.x / ellipsoid.equatorRadius)
-//                res.longitude,
+            )
-//            )
         }
     }
 

maps-kt-core/src/commonTest/kotlin/center/sciprog/maps/coordinates/MercatorTest.kt

@@ -6,14 +6,24 @@ import kotlin.test.assertEquals
 
 class MercatorTest {
     @Test
-    fun forwardBackward(){
+    fun sphereForwardBackward(){
         val moscow = Gmc.ofDegrees(55.76058287719673, 37.60358622841869)
         val mercator = MapProjection.epsg3857.toProjection(moscow)
         //https://epsg.io/transform#s_srs=4326&t_srs=3857&x=37.6035862&y=55.7605829
         assertEquals(4186.0120709, mercator.x.kilometers, 1e-4)
         assertEquals(7510.9013658, mercator.y.kilometers, 1e-4)
         val backwards = MapProjection.epsg3857.toGeodetic(mercator)
-        assertEquals(moscow.latitude.degrees, backwards.latitude.degrees, 0.001)
+        assertEquals(moscow.latitude.degrees, backwards.latitude.degrees, 1e-6)
-        assertEquals(moscow.longitude.degrees, backwards.longitude.degrees, 0.001)
+        assertEquals(moscow.longitude.degrees, backwards.longitude.degrees, 1e-6)
+    }
+
+    @Test
+    fun ellipseForwardBackward(){
+        val moscow = Gmc.ofDegrees(55.76058287719673, 37.60358622841869)
+        val projection = MercatorProjection(ellipsoid = GeoEllipsoid.WGS84)
+        val mercator = projection.toProjection(moscow)
+        val backwards = projection.toGeodetic(mercator)
+        assertEquals(moscow.latitude.degrees, backwards.latitude.degrees, 1e-6)
+        assertEquals(moscow.longitude.degrees, backwards.longitude.degrees, 1e-6)
     }
 }