From 25e9eed88dbaa6499b5ec08004081fa1fc29f476 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sat, 10 Sep 2022 18:04:42 +0300 Subject: [PATCH] Correct mercator projection. Use epsg naming --- .../center/sciprog/maps/compose/MapViewJvm.kt | 6 +- .../maps/coordinates/MercatorProjection.kt | 62 +++++++++++-------- .../sciprog/maps/coordinates/MercatorTest.kt | 18 ++++++ 3 files changed, 60 insertions(+), 26 deletions(-) create mode 100644 maps-kt-core/src/commonTest/kotlin/center/sciprog/maps/coordinates/MercatorTest.kt diff --git a/maps-kt-compose/src/jvmMain/kotlin/center/sciprog/maps/compose/MapViewJvm.kt b/maps-kt-compose/src/jvmMain/kotlin/center/sciprog/maps/compose/MapViewJvm.kt index b8951fe..0416684 100644 --- a/maps-kt-compose/src/jvmMain/kotlin/center/sciprog/maps/compose/MapViewJvm.kt +++ b/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()) diff --git a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/MercatorProjection.kt b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/MercatorProjection.kt index 03d21f5..b411e84 100644 --- a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/MercatorProjection.kt +++ b/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{ +public interface MapProjection { 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 { + public val ellipsoid: GeoEllipsoid = GeoEllipsoid.sphere, +) : MapProjection { 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 } } \ No newline at end of file diff --git a/maps-kt-core/src/commonTest/kotlin/center/sciprog/maps/coordinates/MercatorTest.kt b/maps-kt-core/src/commonTest/kotlin/center/sciprog/maps/coordinates/MercatorTest.kt new file mode 100644 index 0000000..a4018f6 --- /dev/null +++ b/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) + } +} \ No newline at end of file