Correct mercator projection. Use epsg naming

maps-kt-compose/src/jvmMain/kotlin/center/sciprog/maps/compose/MapViewJvm.kt

@@ -251,7 +251,11 @@ public actual fun MapView(
                     val bottomRight = feature.oval.bottomRight.toOffset()
 
                     val path = Path().apply {
-                        addArcRad(Rect(topLeft, bottomRight), feature.startAngle, feature.endAngle - feature.startAngle)
+                        addArcRad(
+                            Rect(topLeft, bottomRight),
+                            feature.startAngle.radians.value.toFloat(),
+                            (feature.endAngle - feature.startAngle).radians.value.toFloat()
+                        )
                     }
 
                     drawPath(path, color = feature.color, style = Stroke())

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

@@ -6,45 +6,46 @@
 package center.sciprog.maps.coordinates
 
 import center.sciprog.maps.coordinates.Angle.Companion.pi
-import kotlin.math.atan
-import kotlin.math.ln
-import kotlin.math.sinh
+import kotlin.math.*
 
 public data class ProjectionCoordinates(val x: Distance, val y: Distance)
 
 /**
  * @param T the type of projection coordinates
  */
-public interface MapProjection<T: Any>{
+public interface MapProjection<T : Any> {
     public fun toGeodetic(pc: T): GeodeticMapCoordinates
     public fun toProjection(gmc: GeodeticMapCoordinates): T
+
+    public companion object{
+        public val epsg3857: MercatorProjection = MercatorProjection()
+    }
 }
 
 
 /**
  * @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
+ * @param ellipsoid - a [GeoEllipsoid] to be used for conversion
  */
 public open class MercatorProjection(
     public val baseLongitude: Angle = Angle.zero,
-    protected val radius: Distance = DEFAULT_EARTH_RADIUS,
-    private val correctedRadius: ((GeodeticMapCoordinates) -> Distance)? = null,
-): MapProjection<ProjectionCoordinates> {
+    public val ellipsoid: GeoEllipsoid = GeoEllipsoid.sphere,
+) : MapProjection<ProjectionCoordinates> {
 
     override fun toGeodetic(pc: ProjectionCoordinates): GeodeticMapCoordinates {
         val res = GeodeticMapCoordinates.ofRadians(
-            atan(sinh(pc.y / radius)),
-            baseLongitude.radians.value + (pc.x / radius),
+            atan(sinh(pc.y / ellipsoid.equatorRadius)),
+            baseLongitude.radians.value + (pc.x / ellipsoid.equatorRadius),
         )
-        return if (correctedRadius != null) {
-            val r = correctedRadius.invoke(res)
-            GeodeticMapCoordinates.ofRadians(
-                atan(sinh(pc.y / r)),
-                baseLongitude.radians.value + pc.x / r,
-            )
-        } else {
+
+        return if (ellipsoid === GeoEllipsoid.sphere) {
             res
+        } else {
+            TODO("Elliptical mercator projection not implemented")
+//            GeodeticMapCoordinates(
+//                atan(sinh(pc.y / ellipsoid.polarRadius)).radians,
+//                res.longitude,
+//            )
         }
     }
 
@@ -53,15 +54,26 @@ public open class MercatorProjection(
      */
     override fun toProjection(gmc: GeodeticMapCoordinates): ProjectionCoordinates {
         require(abs(gmc.latitude) <= MAXIMUM_LATITUDE) { "Latitude exceeds the maximum latitude for mercator coordinates" }
-        val r: Distance = correctedRadius?.invoke(gmc) ?: radius
-        return ProjectionCoordinates(
-            x = r * (gmc.longitude - baseLongitude).radians.value,
-            y = r * ln(tan(pi / 4 + gmc.latitude / 2))
-        )
+        val e = sqrt(ellipsoid.eSquared)
+
+        return if (ellipsoid === GeoEllipsoid.sphere) {
+            ProjectionCoordinates(
+                x = ellipsoid.equatorRadius * (gmc.longitude - baseLongitude).radians.value,
+                y = ellipsoid.equatorRadius * ln(tan(pi / 4 + gmc.latitude / 2))
+            )
+        } else {
+            val sinPhi = sin(gmc.latitude)
+            ProjectionCoordinates(
+                x = ellipsoid.equatorRadius * (gmc.longitude - baseLongitude).radians.value,
+                y = ellipsoid.equatorRadius * ln(
+                    tan(pi / 4 + gmc.latitude / 2) * ((1 - e * sinPhi) / (1 + e * sinPhi)).pow(e / 2)
+                )
+            )
+        }
+
     }
 
-    public companion object : MercatorProjection(Angle.zero, Distance(6378.137)) {
+    public companion object {
         public val MAXIMUM_LATITUDE: Angle = 85.05113.degrees
-        public val DEFAULT_EARTH_RADIUS: Distance = radius
     }
 }

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

@@ -0,0 +1,18 @@
+package center.sciprog.maps.coordinates
+
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class MercatorTest {
+    @Test
+    fun forwardBackward(){
+        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.value, backwards.latitude.degrees.value, 0.001)
+        assertEquals(moscow.longitude.degrees.value, backwards.longitude.degrees.value, 0.001)
+    }
+}