From 4614b1f7bc5b014b3d909df0e3240f3696c118b6 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Mon, 29 Aug 2022 13:38:46 +0300 Subject: [PATCH] Type safe angles and distances --- .../src/jvmMain/kotlin/AngleConversion.kt | 4 - demo/maps/src/jvmMain/kotlin/Main.kt | 12 +- .../center/sciprog/maps/compose/MapFeature.kt | 32 +- .../sciprog/maps/compose/MapFeatureBuilder.kt | 6 +- .../center/sciprog/maps/compose/MapView.kt | 25 +- .../sciprog/maps/compose/MapTextFeature.kt | 4 +- .../center/sciprog/maps/compose/MapViewJvm.kt | 26 +- maps-kt-core/build.gradle.kts | 8 + .../sciprog/maps/coordinates/Distance.kt | 5 +- .../sciprog/maps/coordinates/Ellipsoid.kt | 14 - .../sciprog/maps/coordinates/GeoEllipsoid.kt | 72 ++++ .../coordinates/GeodeticMapCoordinates.kt | 32 +- .../center/sciprog/maps/coordinates/GmcBox.kt | 77 ---- .../sciprog/maps/coordinates/GmcCurve.kt | 348 ++++++++++++++++++ .../sciprog/maps/coordinates/GmcPose.kt | 12 + .../sciprog/maps/coordinates/GmcRectangle.kt | 82 +++++ .../sciprog/maps/coordinates/MapViewPoint.kt | 4 +- .../maps/coordinates/MercatorProjection.kt | 29 +- .../maps/coordinates/WebMercatorProjection.kt | 6 +- .../center/sciprog/maps/coordinates/angles.kt | 92 +++++ .../sciprog/maps/coordinates/DistanceTest.kt | 37 ++ 21 files changed, 752 insertions(+), 175 deletions(-) delete mode 100644 demo/maps/src/jvmMain/kotlin/AngleConversion.kt delete mode 100644 maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/Ellipsoid.kt create mode 100644 maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GeoEllipsoid.kt delete mode 100644 maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcBox.kt create mode 100644 maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcCurve.kt create mode 100644 maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcPose.kt create mode 100644 maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcRectangle.kt create mode 100644 maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/angles.kt create mode 100644 maps-kt-core/src/commonTest/kotlin/center/sciprog/maps/coordinates/DistanceTest.kt diff --git a/demo/maps/src/jvmMain/kotlin/AngleConversion.kt b/demo/maps/src/jvmMain/kotlin/AngleConversion.kt deleted file mode 100644 index a7ec5e5..0000000 --- a/demo/maps/src/jvmMain/kotlin/AngleConversion.kt +++ /dev/null @@ -1,4 +0,0 @@ -import kotlin.math.PI - -fun Double.toDegrees() = this * 180 / PI - diff --git a/demo/maps/src/jvmMain/kotlin/Main.kt b/demo/maps/src/jvmMain/kotlin/Main.kt index 9660683..6bf40a6 100644 --- a/demo/maps/src/jvmMain/kotlin/Main.kt +++ b/demo/maps/src/jvmMain/kotlin/Main.kt @@ -23,7 +23,7 @@ import kotlin.math.PI import kotlin.random.Random private fun GeodeticMapCoordinates.toShortString(): String = - "${(latitude * 180.0 / PI).toString().take(6)}:${(longitude * 180.0 / PI).toString().take(6)}" + "${(latitude.degrees.value).toString().take(6)}:${(longitude.degrees.value).toString().take(6)}" @Composable @@ -50,7 +50,7 @@ fun App() { val pointOne = 55.568548 to 37.568604 - var pointTwo by remember { mutableStateOf(55.929444 to 37.518434) } + var pointTwo: Pair by remember { mutableStateOf(55.929444 to 37.518434) } val pointThree = 60.929444 to 37.518434 MapView( mapTileProvider = mapTileProvider, @@ -59,11 +59,11 @@ fun App() { inferViewBoxFromFeatures = true, onViewChange = { centerCoordinates = focus }, onDrag = { start, end -> - if (start.focus.latitude.toDegrees() in (pointTwo.first - 0.05)..(pointTwo.first + 0.05) && - start.focus.longitude.toDegrees() in (pointTwo.second - 0.05)..(pointTwo.second + 0.05) + if (start.focus.latitude.degrees.value in (pointTwo.first - 0.05)..(pointTwo.first + 0.05) && + start.focus.longitude.degrees.value in (pointTwo.second - 0.05)..(pointTwo.second + 0.05) ) { - pointTwo = pointTwo.first + (end.focus.latitude - start.focus.latitude).toDegrees() to - pointTwo.second + (end.focus.longitude - start.focus.longitude).toDegrees() + pointTwo = pointTwo.first + (end.focus.latitude - start.focus.latitude).degrees.value to + pointTwo.second + (end.focus.longitude - start.focus.longitude).degrees.value false// returning false, because when we are dragging circle we don't want to drag map } else true } diff --git a/maps-kt-compose/src/commonMain/kotlin/center/sciprog/maps/compose/MapFeature.kt b/maps-kt-compose/src/commonMain/kotlin/center/sciprog/maps/compose/MapFeature.kt index 202e30b..355284f 100644 --- a/maps-kt-compose/src/commonMain/kotlin/center/sciprog/maps/compose/MapFeature.kt +++ b/maps-kt-compose/src/commonMain/kotlin/center/sciprog/maps/compose/MapFeature.kt @@ -12,15 +12,15 @@ import androidx.compose.ui.unit.DpSize import androidx.compose.ui.unit.IntSize import androidx.compose.ui.unit.dp import center.sciprog.maps.coordinates.GeodeticMapCoordinates -import center.sciprog.maps.coordinates.GmcBox +import center.sciprog.maps.coordinates.GmcRectangle import center.sciprog.maps.coordinates.wrapAll public interface MapFeature { public val zoomRange: IntRange - public fun getBoundingBox(zoom: Int): GmcBox? + public fun getBoundingBox(zoom: Int): GmcRectangle? } -public fun Iterable.computeBoundingBox(zoom: Int): GmcBox? = +public fun Iterable.computeBoundingBox(zoom: Int): GmcRectangle? = mapNotNull { it.getBoundingBox(zoom) }.wrapAll() internal fun Pair.toCoordinates() = GeodeticMapCoordinates.ofDegrees(first, second) @@ -35,7 +35,7 @@ public class MapFeatureSelector( ) : MapFeature { override val zoomRange: IntRange get() = defaultZoomRange - override fun getBoundingBox(zoom: Int): GmcBox? = selector(zoom).getBoundingBox(zoom) + override fun getBoundingBox(zoom: Int): GmcRectangle? = selector(zoom).getBoundingBox(zoom) } public class MapDrawFeature( @@ -43,9 +43,9 @@ public class MapDrawFeature( override val zoomRange: IntRange = defaultZoomRange, public val drawFeature: DrawScope.() -> Unit, ) : MapFeature { - override fun getBoundingBox(zoom: Int): GmcBox { + override fun getBoundingBox(zoom: Int): GmcRectangle { //TODO add box computation - return GmcBox(position, position) + return GmcRectangle(position, position) } } @@ -56,8 +56,8 @@ public class MapPointsFeature( public val color: Color = Color.Red, public val pointMode: PointMode = PointMode.Points ) : MapFeature { - override fun getBoundingBox(zoom: Int): GmcBox { - return GmcBox(points.first(), points.last()) + override fun getBoundingBox(zoom: Int): GmcRectangle { + return GmcRectangle(points.first(), points.last()) } } @@ -67,7 +67,7 @@ public class MapCircleFeature( public val size: Float = 5f, public val color: Color = Color.Red, ) : MapFeature { - override fun getBoundingBox(zoom: Int): GmcBox = GmcBox(center, center) + override fun getBoundingBox(zoom: Int): GmcRectangle = GmcRectangle(center, center) } public class MapRectangleFeature( @@ -76,7 +76,7 @@ public class MapRectangleFeature( public val size: DpSize = DpSize(5.dp, 5.dp), public val color: Color = Color.Red, ) : MapFeature { - override fun getBoundingBox(zoom: Int): GmcBox = GmcBox(center, center) + override fun getBoundingBox(zoom: Int): GmcRectangle = GmcRectangle(center, center) } public class MapLineFeature( @@ -85,17 +85,17 @@ public class MapLineFeature( override val zoomRange: IntRange = defaultZoomRange, public val color: Color = Color.Red, ) : MapFeature { - override fun getBoundingBox(zoom: Int): GmcBox = GmcBox(a, b) + override fun getBoundingBox(zoom: Int): GmcRectangle = GmcRectangle(a, b) } public class MapArcFeature( - public val oval: GmcBox, + public val oval: GmcRectangle, public val startAngle: Float, public val endAngle: Float, override val zoomRange: IntRange = defaultZoomRange, public val color: Color = Color.Red, ) : MapFeature { - override fun getBoundingBox(zoom: Int): GmcBox = oval + override fun getBoundingBox(zoom: Int): GmcRectangle = oval } public class MapBitmapImageFeature( @@ -104,7 +104,7 @@ public class MapBitmapImageFeature( public val size: IntSize = IntSize(15, 15), override val zoomRange: IntRange = defaultZoomRange, ) : MapFeature { - override fun getBoundingBox(zoom: Int): GmcBox = GmcBox(position, position) + override fun getBoundingBox(zoom: Int): GmcRectangle = GmcRectangle(position, position) } public class MapVectorImageFeature( @@ -113,7 +113,7 @@ public class MapVectorImageFeature( public val size: DpSize, override val zoomRange: IntRange = defaultZoomRange, ) : MapFeature { - override fun getBoundingBox(zoom: Int): GmcBox = GmcBox(position, position) + override fun getBoundingBox(zoom: Int): GmcRectangle = GmcRectangle(position, position) } @Composable @@ -131,5 +131,5 @@ public class MapFeatureGroup( public val children: Map, override val zoomRange: IntRange = defaultZoomRange, ) : MapFeature { - override fun getBoundingBox(zoom: Int): GmcBox? = children.values.mapNotNull { it.getBoundingBox(zoom) }.wrapAll() + override fun getBoundingBox(zoom: Int): GmcRectangle? = children.values.mapNotNull { it.getBoundingBox(zoom) }.wrapAll() } \ No newline at end of file diff --git a/maps-kt-compose/src/commonMain/kotlin/center/sciprog/maps/compose/MapFeatureBuilder.kt b/maps-kt-compose/src/commonMain/kotlin/center/sciprog/maps/compose/MapFeatureBuilder.kt index e916f12..f3c8a48 100644 --- a/maps-kt-compose/src/commonMain/kotlin/center/sciprog/maps/compose/MapFeatureBuilder.kt +++ b/maps-kt-compose/src/commonMain/kotlin/center/sciprog/maps/compose/MapFeatureBuilder.kt @@ -11,7 +11,7 @@ import androidx.compose.ui.unit.DpSize import androidx.compose.ui.unit.dp import center.sciprog.maps.coordinates.Distance import center.sciprog.maps.coordinates.GeodeticMapCoordinates -import center.sciprog.maps.coordinates.GmcBox +import center.sciprog.maps.coordinates.GmcRectangle public typealias FeatureId = String @@ -87,7 +87,7 @@ public fun MapFeatureBuilder.line( ) public fun MapFeatureBuilder.arc( - oval: GmcBox, + oval: GmcRectangle, startAngle: Number, endAngle: Number, zoomRange: IntRange = defaultZoomRange, @@ -109,7 +109,7 @@ public fun MapFeatureBuilder.arc( ): FeatureId = addFeature( id, MapArcFeature( - GmcBox.withCenter(center.toCoordinates(), radius, radius), + GmcRectangle.square(center.toCoordinates(), radius, radius), startAngle.toFloat(), endAngle.toFloat(), zoomRange, diff --git a/maps-kt-compose/src/commonMain/kotlin/center/sciprog/maps/compose/MapView.kt b/maps-kt-compose/src/commonMain/kotlin/center/sciprog/maps/compose/MapView.kt index 0dc4fd4..3771850 100644 --- a/maps-kt-compose/src/commonMain/kotlin/center/sciprog/maps/compose/MapView.kt +++ b/maps-kt-compose/src/commonMain/kotlin/center/sciprog/maps/compose/MapView.kt @@ -21,9 +21,9 @@ public data class MapViewConfig( val onClick: MapViewPoint.() -> Unit = {}, val onDrag: (start: MapViewPoint, end: MapViewPoint) -> Boolean = { _, _ -> true }, val onViewChange: MapViewPoint.() -> Unit = {}, - val onSelect: (GmcBox) -> Unit = {}, + val onSelect: (GmcRectangle) -> Unit = {}, val zoomOnSelect: Boolean = true, - val resetViewPoint: Boolean = false + val resetViewPoint: Boolean = false, ) @Composable @@ -55,20 +55,21 @@ public fun MapView( ) } -internal fun GmcBox.computeViewPoint(mapTileProvider: MapTileProvider): (canvasSize: DpSize) -> MapViewPoint = { canvasSize -> - val zoom = log2( - min( - canvasSize.width.value / width, - canvasSize.height.value / height - ) * PI / mapTileProvider.tileSize - ) - MapViewPoint(center, zoom) -} +internal fun GmcRectangle.computeViewPoint(mapTileProvider: MapTileProvider): (canvasSize: DpSize) -> MapViewPoint = + { canvasSize -> + val zoom = log2( + min( + canvasSize.width.value / longitudeDelta.radians.value, + canvasSize.height.value / latitudeDelta.radians.value + ) * PI / mapTileProvider.tileSize + ) + MapViewPoint(center, zoom) + } @Composable public fun MapView( mapTileProvider: MapTileProvider, - box: GmcBox, + box: GmcRectangle, features: Map = emptyMap(), config: MapViewConfig = MapViewConfig(), modifier: Modifier = Modifier.fillMaxSize(), diff --git a/maps-kt-compose/src/jvmMain/kotlin/center/sciprog/maps/compose/MapTextFeature.kt b/maps-kt-compose/src/jvmMain/kotlin/center/sciprog/maps/compose/MapTextFeature.kt index fef40fb..4a5ebce 100644 --- a/maps-kt-compose/src/jvmMain/kotlin/center/sciprog/maps/compose/MapTextFeature.kt +++ b/maps-kt-compose/src/jvmMain/kotlin/center/sciprog/maps/compose/MapTextFeature.kt @@ -2,7 +2,7 @@ package center.sciprog.maps.compose import androidx.compose.ui.graphics.Color import center.sciprog.maps.coordinates.GeodeticMapCoordinates -import center.sciprog.maps.coordinates.GmcBox +import center.sciprog.maps.coordinates.GmcRectangle import org.jetbrains.skia.Font @@ -13,7 +13,7 @@ public class MapTextFeature( public val color: Color, public val fontConfig: Font.() -> Unit, ) : MapFeature { - override fun getBoundingBox(zoom: Int): GmcBox = GmcBox(position, position) + override fun getBoundingBox(zoom: Int): GmcRectangle = GmcRectangle(position, position) } public fun MapFeatureBuilder.text( 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 27f9496..a57e5d3 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 @@ -30,12 +30,12 @@ private fun Color.toPaint(): Paint = Paint().apply { private fun IntRange.intersect(other: IntRange) = max(first, other.first)..min(last, other.last) internal fun MapViewPoint.move(deltaX: Double, deltaY: Double): MapViewPoint { - val newCoordinates = GeodeticMapCoordinates.ofRadians( - (focus.latitude + deltaY / scaleFactor).coerceIn( + val newCoordinates = GeodeticMapCoordinates( + (focus.latitude + (deltaY / scaleFactor).radians).coerceIn( -MercatorProjection.MAXIMUM_LATITUDE, MercatorProjection.MAXIMUM_LATITUDE ), - focus.longitude + deltaX / scaleFactor + focus.longitude + (deltaX / scaleFactor).radians ) return MapViewPoint(newCoordinates, zoom) } @@ -69,8 +69,8 @@ public actual fun MapView( features.values.computeBoundingBox(1)?.let { box -> val zoom = log2( min( - canvasSize.width.value / box.width, - canvasSize.height.value / box.height + canvasSize.width.value / box.longitudeDelta.radians.value, + canvasSize.height.value / box.latitudeDelta.radians.value ) * PI / mapTileProvider.tileSize ) MapViewPoint(box.center, zoom) @@ -127,7 +127,7 @@ public actual fun MapView( } selectRect?.let { rect -> //Use selection override if it is defined - val gmcBox = GmcBox( + val gmcBox = GmcRectangle( rect.topLeft.toDpOffset().toGeodetic(), rect.bottomRight.toDpOffset().toGeodetic() ) @@ -146,8 +146,10 @@ public actual fun MapView( config.onClick(MapViewPoint(dpPos.toGeodetic(), viewPoint.zoom)) drag(change.id) { dragChange -> val dragAmount = dragChange.position - dragChange.previousPosition - val dpStart = - DpOffset(dragChange.previousPosition.x.toDp(), dragChange.previousPosition.y.toDp()) + val dpStart = DpOffset( + dragChange.previousPosition.x.toDp(), + dragChange.previousPosition.y.toDp() + ) val dpEnd = DpOffset(dragChange.position.x.toDp(), dragChange.position.y.toDp()) if (!config.onDrag( MapViewPoint(dpStart.toGeodetic(), viewPoint.zoom), @@ -232,6 +234,7 @@ public actual fun MapView( feature.size, center = feature.center.toOffset() ) + is MapRectangleFeature -> drawRect( feature.color, topLeft = feature.center.toOffset() - Offset( @@ -240,6 +243,7 @@ public actual fun MapView( ), size = feature.size.toSize() ) + is MapLineFeature -> drawLine(feature.color, feature.a.toOffset(), feature.b.toOffset()) is MapArcFeature -> { val topLeft = feature.oval.topLeft.toOffset() @@ -252,6 +256,7 @@ public actual fun MapView( drawPath(path, color = feature.color, style = Stroke()) } + is MapBitmapImageFeature -> drawImage(feature.image, feature.position.toOffset()) is MapVectorImageFeature -> { val offset = feature.position.toOffset() @@ -262,6 +267,7 @@ public actual fun MapView( } } } + is MapTextFeature -> drawIntoCanvas { canvas -> val offset = feature.position.toOffset() canvas.nativeCanvas.drawString( @@ -272,17 +278,20 @@ public actual fun MapView( feature.color.toPaint() ) } + is MapDrawFeature -> { val offset = feature.position.toOffset() translate(offset.x, offset.y) { feature.drawFeature(this) } } + is MapFeatureGroup -> { feature.children.values.forEach { drawFeature(zoom, it) } } + is MapPointsFeature -> { val points = feature.points.map { it.toOffset() } drawPoints( @@ -292,6 +301,7 @@ public actual fun MapView( pointMode = feature.pointMode ) } + else -> { logger.error { "Unrecognized feature type: ${feature::class}" } } diff --git a/maps-kt-core/build.gradle.kts b/maps-kt-core/build.gradle.kts index 5ca05d2..7924005 100644 --- a/maps-kt-core/build.gradle.kts +++ b/maps-kt-core/build.gradle.kts @@ -15,4 +15,12 @@ kotlin { js(IR) { browser() } + + sourceSets{ + commonTest{ + dependencies{ + implementation(kotlin("test")) + } + } + } } \ No newline at end of file diff --git a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/Distance.kt b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/Distance.kt index dbe35c8..5fc1349 100644 --- a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/Distance.kt +++ b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/Distance.kt @@ -3,7 +3,9 @@ package center.sciprog.maps.coordinates import kotlin.jvm.JvmInline @JvmInline -public value class Distance(public val kilometers: Double) +public value class Distance(public val kilometers: Double) : Comparable { + override fun compareTo(other: Distance): Int = this.kilometers.compareTo(other.kilometers) +} public operator fun Distance.div(other: Distance): Double = kilometers / other.kilometers @@ -13,5 +15,4 @@ public operator fun Distance.minus(other: Distance): Distance = Distance(kilomet public operator fun Distance.times(number: Number): Distance = Distance(kilometers * number.toDouble()) public operator fun Distance.div(number: Number): Distance = Distance(kilometers / number.toDouble()) - public val Distance.meters: Double get() = kilometers * 1000 \ No newline at end of file diff --git a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/Ellipsoid.kt b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/Ellipsoid.kt deleted file mode 100644 index d99aab1..0000000 --- a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/Ellipsoid.kt +++ /dev/null @@ -1,14 +0,0 @@ -package center.sciprog.maps.coordinates - -public class Ellipsoid(public val equatorRadius: Distance, public val polarRadius: Distance) { - public companion object { - public val WGS84: Ellipsoid = Ellipsoid( - equatorRadius = Distance(6378.137), - polarRadius = Distance(6356.7523142) - ) - } -} - -public val Ellipsoid.f: Double get() = (equatorRadius.kilometers - polarRadius.kilometers) / equatorRadius.kilometers - -public val Ellipsoid.inverseF: Double get() = equatorRadius.kilometers / (equatorRadius.kilometers - polarRadius.kilometers) \ No newline at end of file diff --git a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GeoEllipsoid.kt b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GeoEllipsoid.kt new file mode 100644 index 0000000..8e07710 --- /dev/null +++ b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GeoEllipsoid.kt @@ -0,0 +1,72 @@ +package center.sciprog.maps.coordinates + +import center.sciprog.maps.coordinates.GeoEllipsoid.Companion.greatCircleAngleBetween +import kotlin.math.* + + +public class GeoEllipsoid(public val equatorRadius: Distance, public val polarRadius: Distance) { + + /** + * Flattening https://en.wikipedia.org/wiki/Flattening + */ + public val f: Double = (equatorRadius.kilometers - polarRadius.kilometers) / equatorRadius.kilometers + + /** + * Inverse flattening + */ + public val inverseF: Double = equatorRadius.kilometers / (equatorRadius.kilometers - polarRadius.kilometers) + + public val eSquared: Double = 2 * f - f * f + + public companion object { + + public val WGS84: GeoEllipsoid = GeoEllipsoid( + equatorRadius = Distance(6378.137), + polarRadius = Distance(6356.752314245) + ) + + public val GRS80: GeoEllipsoid = GeoEllipsoid( + equatorRadius = Distance(6378.137), + polarRadius = Distance(6356.752314140) + ) + + public val sphere: GeoEllipsoid = GeoEllipsoid( + equatorRadius = Distance(6378.137), + polarRadius = Distance(6378.137) + ) + + /** + * https://en.wikipedia.org/wiki/Great-circle_distance + */ + public fun greatCircleAngleBetween(r1: GMC, r2: GMC): Radians = acos( + sin(r1.latitude) * sin(r2.latitude) + cos(r1.latitude) * cos(r2.latitude) * cos(r1.longitude - r2.longitude) + ).radians + } +} + +/** + * A radius of circle normal to the axis of the ellipsoid at given latitude + */ +internal fun GeoEllipsoid.reducedRadius(latitude: Angle): Distance { + val reducedLatitudeTan = (1 - f) * tan(latitude) + return equatorRadius / sqrt(1.0 + reducedLatitudeTan.pow(2)) +} + + +/** + * Compute distance between two map points using giv + * https://en.wikipedia.org/wiki/Geographical_distance#Lambert's_formula_for_long_lines + */ +public fun GeoEllipsoid.lambertDistanceBetween(r1: GMC, r2: GMC): Distance { + val s = greatCircleAngleBetween(r1, r2) + + val b1: Double = (1 - f) * tan(r1.latitude) + val b2: Double = (1 - f) * tan(r2.latitude) + val p = (b1 + b2) / 2 + val q = (b2 - b1) / 2 + + val x = (s.value - sin(s)) * sin(p).pow(2) * cos(q).pow(2) / cos(s / 2).pow(2) + val y = (s.value + sin(s)) * cos(p).pow(2) * sin(q).pow(2) / sin(s / 2).pow(2) + + return equatorRadius * (s.value - f / 2 * (x + y)) +} \ No newline at end of file diff --git a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GeodeticMapCoordinates.kt b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GeodeticMapCoordinates.kt index c7a75eb..d12c310 100644 --- a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GeodeticMapCoordinates.kt +++ b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GeodeticMapCoordinates.kt @@ -5,10 +5,15 @@ import kotlin.math.PI /** * Geodetic coordinated */ -public class GeodeticMapCoordinates private constructor( - public val latitude: Double, - public val longitude: Double, -){ +public class GeodeticMapCoordinates( + public val latitude: Angle, + longitude: Angle, +) { + public val longitude: Radians = longitude.radians.value.rem(PI / 2).radians + + init { + require(latitude.radians.value in (-PI / 2)..(PI / 2)) { "Latitude $latitude is not in (-PI/2)..(PI/2)" } + } override fun equals(other: Any?): Boolean { if (this === other) return true @@ -29,23 +34,24 @@ public class GeodeticMapCoordinates private constructor( } override fun toString(): String { - return "GeodeticCoordinates(latitude=${latitude / PI * 180} deg, longitude=${longitude / PI * 180} deg)" + return "GMC(latitude=${latitude.degrees.value} deg, longitude=${longitude.degrees.value} deg)" } public companion object { - public fun ofRadians(latitude: Double, longitude: Double): GeodeticMapCoordinates { - require(latitude in (-PI/2)..(PI/2)) { "Latitude $latitude is not in (-PI/2)..(PI/2)" } - return GeodeticMapCoordinates(latitude, longitude.rem(PI / 2)) - } + public fun ofRadians(latitude: Double, longitude: Double): GeodeticMapCoordinates = + GeodeticMapCoordinates(latitude.radians, longitude.radians) - public fun ofDegrees(latitude: Double, longitude: Double): GeodeticMapCoordinates { - require(latitude in (-90.0)..(90.0)) { "Latitude $latitude is not in -90..90" } - return GeodeticMapCoordinates(latitude * PI / 180, (longitude.rem(180) * PI / 180)) - } + public fun ofDegrees(latitude: Double, longitude: Double): GeodeticMapCoordinates = + GeodeticMapCoordinates(latitude.degrees.radians, longitude.degrees.radians) } } +/** + * Short name for GeodeticMapCoordinates + */ +public typealias GMC = GeodeticMapCoordinates + //public interface GeoToScreenConversion { // public fun getScreenX(gmc: GeodeticMapCoordinates): Double // public fun getScreenY(gmc: GeodeticMapCoordinates): Double diff --git a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcBox.kt b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcBox.kt deleted file mode 100644 index 4eb1ddc..0000000 --- a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcBox.kt +++ /dev/null @@ -1,77 +0,0 @@ -package center.sciprog.maps.coordinates - -import kotlin.math.abs -import kotlin.math.cos -import kotlin.math.max -import kotlin.math.min - -public data class GmcBox( - public val a: GeodeticMapCoordinates, - public val b: GeodeticMapCoordinates, -) { - public companion object { - public fun withCenter( - center: GeodeticMapCoordinates, - width: Distance, - height: Distance, - ellipsoid: Ellipsoid = Ellipsoid.WGS84, - ): GmcBox { - val r = ellipsoid.equatorRadius * cos(center.latitude) - val a = GeodeticMapCoordinates.ofRadians( - center.latitude - height / ellipsoid.polarRadius / 2, - center.longitude - width / r / 2 - ) - val b = GeodeticMapCoordinates.ofRadians( - center.latitude + height / ellipsoid.polarRadius / 2, - center.longitude + width / r / 2 - ) - return GmcBox(a, b) - } - } -} - -public val GmcBox.center: GeodeticMapCoordinates - get() = GeodeticMapCoordinates.ofRadians( - (a.latitude + b.latitude) / 2, - (a.longitude + b.longitude) / 2 - ) - -/** - * Minimum longitude - */ -public val GmcBox.left: Double get() = min(a.longitude, b.longitude) - -/** - * maximum longitude - */ -public val GmcBox.right: Double get() = max(a.longitude, b.longitude) - -/** - * Maximum latitude - */ -public val GmcBox.top: Double get() = max(a.latitude, b.latitude) - -/** - * Minimum latitude - */ -public val GmcBox.bottom: Double get() = min(a.latitude, b.latitude) - -//TODO take curvature into account -public val GmcBox.width: Double get() = abs(a.longitude - b.longitude) -public val GmcBox.height: Double get() = abs(a.latitude - b.latitude) - -public val GmcBox.topLeft: GeodeticMapCoordinates get() = GeodeticMapCoordinates.ofRadians(top, left) -public val GmcBox.bottomRight: GeodeticMapCoordinates get() = GeodeticMapCoordinates.ofRadians(bottom, right) - -/** - * Compute a minimal bounding box including all given boxes. Return null if collection is empty - */ -public fun Collection.wrapAll(): GmcBox? { - if (isEmpty()) return null - //TODO optimize computation - val minLat = minOf { it.bottom } - val maxLat = maxOf { it.top } - val minLong = minOf { it.left } - val maxLong = maxOf { it.right } - return GmcBox(GeodeticMapCoordinates.ofRadians(minLat, minLong), GeodeticMapCoordinates.ofRadians(maxLat, maxLong)) -} \ No newline at end of file diff --git a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcCurve.kt b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcCurve.kt new file mode 100644 index 0000000..c13f51e --- /dev/null +++ b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcCurve.kt @@ -0,0 +1,348 @@ +package center.sciprog.maps.coordinates + +import center.sciprog.maps.coordinates.Angle.Companion.pi +import center.sciprog.maps.coordinates.Angle.Companion.piDiv2 +import center.sciprog.maps.coordinates.Angle.Companion.zero +import kotlin.math.* + +/** + * A directed straight (geodetic) segment on a spheroid with given start, direction, end point and distance. + * @param forward coordinate of a start point with forward direction + * @param backward coordinate of an end point with backward direction + */ +public class GmcCurve internal constructor( + public val forward: GmcPose, + public val backward: GmcPose, + public val distance: Distance, +) + +public fun GmcCurve.reversed(): GmcCurve = GmcCurve(backward, forward, distance) + +/** + * Compute a curve alongside a meridian + */ +public fun GeoEllipsoid.meridianCurve( + longitude: Angle, + fromLatitude: Angle, + toLatitude: Angle, + step: Radians = 0.015.radians, +): GmcCurve { + require(fromLatitude in (-piDiv2)..(piDiv2)) { "Latitude must be in (-90, 90) degrees range" } + require(toLatitude in (-piDiv2)..(piDiv2)) { "Latitude must be in (-90, 90) degrees range" } + + fun smallDistance(from: Radians, to: Radians): Distance = equatorRadius * + (1 - eSquared) * + (1 - eSquared * sin(from).pow(2)).pow(-1.5) * + abs((from - to).value) + + val up = toLatitude > fromLatitude + + val integrateFrom: Radians + val integrateTo: Radians + + if (up) { + integrateFrom = fromLatitude.radians + integrateTo = toLatitude.radians + } else { + integrateTo = fromLatitude.radians + integrateFrom = toLatitude.radians + } + + var current = integrateFrom + var s = Distance(0.0) + while (current < integrateTo) { + val next = minOf(current + step, integrateTo) + s += smallDistance(current, next) + current = next + } + + return GmcCurve( + forward = GmcPose(GMC(fromLatitude, longitude), if (up) zero else pi), + backward = GmcPose(GMC(toLatitude, longitude), if (up) pi else zero), + distance = s + ) +} + +/** + * Compute a curve alongside a parallel + */ +public fun GeoEllipsoid.parallelCurve(latitude: Angle, fromLongitude: Angle, toLongitude: Angle): GmcCurve { + require(latitude in (-piDiv2)..(piDiv2)) { "Latitude must be in (-90, 90) degrees range" } + val right = toLongitude > fromLongitude + return GmcCurve( + forward = GmcPose(GMC(latitude, fromLongitude), if (right) piDiv2.radians else -piDiv2.radians), + backward = GmcPose(GMC(latitude, toLongitude), if (right) -piDiv2.radians else piDiv2.radians), + distance = reducedRadius(latitude) * abs((fromLongitude - toLongitude).radians.value) + ) +} + + +/** + * Taken from https://github.com/mgavaghan/geodesy + * https://github.com/mgavaghan/geodesy/blob/ab1c6969dc964ff34929911f055b27525909ef3f/src/main/java/org/gavaghan/geodesy/GeodeticCalculator.java#L58 + * + * Calculate the destination and final bearing after traveling a specified + * distance, and a specified starting bearing, for an initial location. This + * is the solution to the direct geodetic problem. + * + * @param start starting location + * @return solution to the direct geodetic problem + */ +public fun GeoEllipsoid.curveInDirection( + start: GmcPose, + distance: Distance, + precision: Double = 1e-6, +): GmcCurve { + val a: Distance = equatorRadius + val b: Distance = polarRadius + val aSquared: Double = a.kilometers.pow(2) + val bSquared: Double = b.kilometers.pow(2) + val phi1 = start.latitude + val alpha1 = start.bearing + val cosAlpha1: Double = cos(alpha1) + val sinAlpha1: Double = sin(alpha1) + val tanU1: Double = (1.0 - f) * tan(phi1) + val cosU1: Double = 1.0 / sqrt(1.0 + tanU1 * tanU1) + val sinU1 = tanU1 * cosU1 + + // eq. 1 + val sigma1: Radians = atan2(tanU1, cosAlpha1).radians + + // eq. 2 + val sinAlpha: Double = cosU1 * sinAlpha1 + val sin2Alpha = sinAlpha * sinAlpha + val cos2Alpha = 1 - sin2Alpha + val uSquared = cos2Alpha * (aSquared - bSquared) / bSquared + + // eq. 3 + val A: Double = 1 + uSquared / 16384 * (4096 + uSquared * (-768 + uSquared * (320 - 175 * uSquared))) + + // eq. 4 + val B: Double = uSquared / 1024 * (256 + uSquared * (-128 + uSquared * (74 - 47 * uSquared))) + + // iterate until there is a negligible change in sigma + val sOverbA: Radians = (distance / b / A).radians + var sigma: Radians = sOverbA + var sinSigma: Double + var prevSigma: Radians = sOverbA + var sigmaM2: Radians + var cosSigmaM2: Double + var cos2SigmaM2: Double + while (!prevSigma.value.isNaN()) { + // eq. 5 + sigmaM2 = sigma1 * 2.0 + sigma + cosSigmaM2 = cos(sigmaM2) + cos2SigmaM2 = cosSigmaM2 * cosSigmaM2 + sinSigma = sin(sigma) +// val cosSigma: Double = cos(sigma) + + // eq. 6 + val deltaSigma = B * sinSigma * + (cosSigmaM2 + B / 4.0 * (cos(sigma) * (-1 + 2 * cos2SigmaM2) - + B / 6.0 * cosSigmaM2 * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM2))) + + // eq. 7 + sigma = sOverbA + deltaSigma.radians + + // break after converging to tolerance + if (abs((sigma - prevSigma).value) < precision) break + prevSigma = sigma + } + sigmaM2 = sigma1 * 2.0 + sigma + cosSigmaM2 = cos(sigmaM2) + cos2SigmaM2 = cosSigmaM2 * cosSigmaM2 + val cosSigma: Double = cos(sigma) + sinSigma = sin(sigma) + + // eq. 8 + val phi2: Radians = atan2( + sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1, + (1.0 - f) * sqrt( + sin2Alpha + (sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1).pow(2) + ) + ).radians + + // eq. 9 + // This fixes the pole crossing defect spotted by Matt Feemster. When a + // path passes a pole and essentially crosses a line of latitude twice - + // once in each direction - the longitude calculation got messed up. + // + // Using atan2 instead of atan fixes the defect. The change is in the + // next 3 lines. + // + // double tanLambda = sinSigma * sinAlpha1 / (cosU1 * cosSigma - sinU1 * + // sinSigma * cosAlpha1); + // double lambda = Math.atan(tanLambda); + val lambda: Double = atan2( + sinSigma * sinAlpha1, + cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1 + ) + + // eq. 10 + val C = f / 16 * cos2Alpha * (4 + f * (4 - 3 * cos2Alpha)) + + // eq. 11 + val L = lambda - (1 - C) * f * sinAlpha * + (sigma.value + C * sinSigma * (cosSigmaM2 + C * cosSigma * (-1 + 2 * cos2SigmaM2))) + + val endPoint = GMC(phi2, L.radians) + + // eq. 12 + + return GmcCurve( + start, + GmcPose( + endPoint, + atan2(sinAlpha, -sinU1 * sinSigma + cosU1 * cosSigma * cosAlpha1).radians + ), + distance + ) +} + +/** + * Taken from https://github.com/mgavaghan/geodesy + * + * Calculate the geodetic curve between two points on a specified reference + * ellipsoid. This is the solution to the inverse geodetic problem. + * + * @receiver reference ellipsoid to use + * @param start starting coordinates + * @param end ending coordinates + * @return solution to the inverse geodetic problem + */ +public fun GeoEllipsoid.curveBetween(start: GMC, end: GMC, precision: Double = 1e-6): GmcCurve { + // + // All equation numbers refer back to Vincenty's publication: + // See http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf + // + + // get constants + val a = equatorRadius + val b = polarRadius + + // get parameters as radians + val phi1 = start.latitude + val lambda1 = start.longitude + val phi2 = end.latitude + val lambda2 = end.longitude + + // calculations + val a2 = a.kilometers * a.kilometers + val b2 = b.kilometers * b.kilometers + val a2b2b2 = (a2 - b2) / b2 + val omega: Radians = lambda2 - lambda1 + val tanphi1: Double = tan(phi1) + val tanU1 = (1.0 - f) * tanphi1 + val U1: Double = atan(tanU1) + val sinU1: Double = sin(U1) + val cosU1: Double = cos(U1) + val tanphi2: Double = tan(phi2) + val tanU2 = (1.0 - f) * tanphi2 + val U2: Double = atan(tanU2) + val sinU2: Double = sin(U2) + val cosU2: Double = cos(U2) + val sinU1sinU2 = sinU1 * sinU2 + val cosU1sinU2 = cosU1 * sinU2 + val sinU1cosU2 = sinU1 * cosU2 + val cosU1cosU2 = cosU1 * cosU2 + + // eq. 13 + var lambda = omega + + // intermediates we'll need to compute 's' + var A = 0.0 + + var sigma = 0.0 + var deltasigma = 0.0 + var lambda0: Radians + var converged = false + for (i in 0..19) { + lambda0 = lambda + val sinlambda: Double = sin(lambda) + val coslambda: Double = cos(lambda) + + // eq. 14 + val sin2sigma = + cosU2 * sinlambda * cosU2 * sinlambda + (cosU1sinU2 - sinU1cosU2 * coslambda) * (cosU1sinU2 - sinU1cosU2 * coslambda) + val sinsigma: Double = sqrt(sin2sigma) + + // eq. 15 + val cossigma = sinU1sinU2 + cosU1cosU2 * coslambda + + // eq. 16 + sigma = atan2(sinsigma, cossigma) + + // eq. 17 Careful! sin2sigma might be almost 0! + val sinalpha = if (sin2sigma == 0.0) 0.0 else cosU1cosU2 * sinlambda / sinsigma + val alpha: Double = asin(sinalpha) + val cosalpha: Double = cos(alpha) + val cos2alpha = cosalpha * cosalpha + + // eq. 18 Careful! cos2alpha might be almost 0! + val cos2sigmam = if (cos2alpha == 0.0) 0.0 else cossigma - 2 * sinU1sinU2 / cos2alpha + val u2 = cos2alpha * a2b2b2 + val cos2sigmam2 = cos2sigmam * cos2sigmam + + // eq. 3 + A = 1.0 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2))) + + // eq. 4 + val B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2))) + + // eq. 6 + deltasigma = + B * sinsigma * (cos2sigmam + B / 4 * (cossigma * (-1 + 2 * cos2sigmam2) - B / 6 * cos2sigmam * (-3 + 4 * sin2sigma) * (-3 + 4 * cos2sigmam2))) + + // eq. 10 + val C = f / 16 * cos2alpha * (4 + f * (4 - 3 * cos2alpha)) + + // eq. 11 (modified) + lambda = omega + ( + (1 - C) * f * sinalpha * + (sigma + C * sinsigma * (cos2sigmam + C * cossigma * (-1 + 2 * cos2sigmam2))) + ).radians + + // see how much improvement we got + val change: Double = abs((lambda - lambda0) / lambda) + if (i > 1 && change < precision) { + converged = true + break + } + } + + // eq. 19 + val s: Distance = b * A * (sigma - deltasigma) + val alpha1: Radians + val alpha2: Radians + + // didn't converge? must be N/S + if (!converged) { + if (phi1 > phi2) { + alpha1 = pi.radians + alpha2 = 0.0.radians + } else if (phi1 < phi2) { + alpha1 = 0.0.radians + alpha2 = pi.radians + } else { + alpha1 = Double.NaN.radians + alpha2 = Double.NaN.radians + } + } else { + // eq. 20 + alpha1 = atan2( + cosU2 * sin(lambda), + cosU1sinU2 - sinU1cosU2 * cos(lambda) + ).radians + + // eq. 21 + alpha2 = atan2( + cosU1 * sin(lambda), + -sinU1cosU2 + cosU1sinU2 * cos(lambda) + ).radians + pi + } + return GmcCurve( + GmcPose(start, alpha1), + GmcPose(end, alpha2), + s + ) +} \ No newline at end of file diff --git a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcPose.kt b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcPose.kt new file mode 100644 index 0000000..e573dc5 --- /dev/null +++ b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcPose.kt @@ -0,0 +1,12 @@ +package center.sciprog.maps.coordinates + +/** + * A coordinate-bearing pair + */ +public data class GmcPose(val coordinates: GeodeticMapCoordinates, val bearing: Angle) { + val latitude: Angle get() = coordinates.latitude + val longitude: Angle get() = coordinates.longitude +} + +public fun GmcPose.reversed(): GmcPose = copy(bearing = (bearing + Angle.pi).normalized()) + diff --git a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcRectangle.kt b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcRectangle.kt new file mode 100644 index 0000000..bb1322c --- /dev/null +++ b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/GmcRectangle.kt @@ -0,0 +1,82 @@ +package center.sciprog.maps.coordinates + +/** + * A section of the map between two parallels and two meridians. The figure represents a square in a Mercator projection. + * Params are two opposing "corners" of quasi-square. + * + * Note that this is a rectangle only on a Mercator projection. + */ +public data class GmcRectangle( + public val a: GeodeticMapCoordinates, + public val b: GeodeticMapCoordinates, + public val ellipsoid: GeoEllipsoid = GeoEllipsoid.WGS84, +) { + public companion object { + + /** + * A quasi-square section. Note that latitudinal distance could be imprecise for large distances + */ + public fun square( + center: GeodeticMapCoordinates, + width: Distance, + height: Distance, + ellipsoid: GeoEllipsoid = GeoEllipsoid.WGS84, + ): GmcRectangle { + val reducedRadius = ellipsoid.reducedRadius(center.latitude) + val a = GeodeticMapCoordinates( + center.latitude - (height / ellipsoid.polarRadius / 2).radians, + center.longitude - (width / reducedRadius / 2).radians + ) + val b = GeodeticMapCoordinates( + center.latitude + (height / ellipsoid.polarRadius / 2).radians, + center.longitude + (width / reducedRadius / 2).radians + ) + return GmcRectangle(a, b, ellipsoid) + } + } +} + +public val GmcRectangle.center: GeodeticMapCoordinates + get() = GeodeticMapCoordinates( + (a.latitude + b.latitude) / 2, + (a.longitude + b.longitude) / 2 + ) + +/** + * Minimum longitude + */ +public val GmcRectangle.left: Angle get() = minOf(a.longitude, b.longitude) + +/** + * maximum longitude + */ +public val GmcRectangle.right: Angle get() = maxOf(a.longitude, b.longitude) + +/** + * Maximum latitude + */ +public val GmcRectangle.top: Angle get() = maxOf(a.latitude, b.latitude) + +/** + * Minimum latitude + */ +public val GmcRectangle.bottom: Angle get() = minOf(a.latitude, b.latitude) + +public val GmcRectangle.longitudeDelta: Angle get() = abs(a.longitude - b.longitude) +public val GmcRectangle.latitudeDelta: Angle get() = abs(a.latitude - b.latitude) + +public val GmcRectangle.topLeft: GeodeticMapCoordinates get() = GeodeticMapCoordinates(top, left) +public val GmcRectangle.bottomRight: GeodeticMapCoordinates get() = GeodeticMapCoordinates(bottom, right) + +/** + * Compute a minimal bounding box including all given boxes. Return null if collection is empty + */ +public fun Collection.wrapAll(): GmcRectangle? { + if (isEmpty()) return null + //TODO optimize computation + val minLat = minOf { it.bottom } + val maxLat = maxOf { it.top } + val minLong = minOf { it.left } + val maxLong = maxOf { it.right } + return GmcRectangle(GeodeticMapCoordinates(minLat, minLong), GeodeticMapCoordinates(maxLat, maxLong)) +} \ No newline at end of file diff --git a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/MapViewPoint.kt b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/MapViewPoint.kt index 022338f..8825c4c 100644 --- a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/MapViewPoint.kt +++ b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/MapViewPoint.kt @@ -13,7 +13,7 @@ public data class MapViewPoint( } public fun MapViewPoint.move(delta: GeodeticMapCoordinates): MapViewPoint { - val newCoordinates = GeodeticMapCoordinates.ofRadians( + val newCoordinates = GeodeticMapCoordinates( (focus.latitude + delta.latitude).coerceIn( -MercatorProjection.MAXIMUM_LATITUDE, MercatorProjection.MAXIMUM_LATITUDE @@ -30,7 +30,7 @@ public fun MapViewPoint.zoom( copy(zoom = (zoom + zoomDelta).coerceIn(2.0, 18.0)) } else { val difScale = (1 - 2.0.pow(-zoomDelta)) - val newCenter = GeodeticMapCoordinates.ofRadians( + val newCenter = GeodeticMapCoordinates( focus.latitude + (invariant.latitude - focus.latitude) * difScale, focus.longitude + (invariant.longitude - focus.longitude) * difScale ) 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 d018993..337a9c3 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 @@ -5,9 +5,12 @@ package center.sciprog.maps.coordinates -import kotlin.math.* +import center.sciprog.maps.coordinates.Angle.Companion.pi +import kotlin.math.atan +import kotlin.math.ln +import kotlin.math.sinh -public data class MercatorCoordinates(val x: Double, val y: Double) +public data class MercatorCoordinates(val x: Distance, val y: Distance) /** * @param baseLongitude the longitude offset in radians @@ -15,21 +18,21 @@ public data class MercatorCoordinates(val x: Double, val y: Double) * @param correctedRadius optional radius correction to account for ellipsoid model */ public open class MercatorProjection( - public val baseLongitude: Double = 0.0, - protected val radius: Double = DEFAULT_EARTH_RADIUS, - private val correctedRadius: ((GeodeticMapCoordinates) -> Double)? = null, + public val baseLongitude: Angle = Angle.zero, + protected val radius: Distance = DEFAULT_EARTH_RADIUS, + private val correctedRadius: ((GeodeticMapCoordinates) -> Distance)? = null, ) { public fun toGeodetic(mc: MercatorCoordinates): GeodeticMapCoordinates { val res = GeodeticMapCoordinates.ofRadians( atan(sinh(mc.y / radius)), - baseLongitude + mc.x / radius, + baseLongitude.radians.value + (mc.x / radius), ) return if (correctedRadius != null) { val r = correctedRadius.invoke(res) GeodeticMapCoordinates.ofRadians( atan(sinh(mc.y / r)), - baseLongitude + mc.x / r, + baseLongitude.radians.value + mc.x / r, ) } else { res @@ -41,15 +44,15 @@ public open class MercatorProjection( */ public fun toMercator(gmc: GeodeticMapCoordinates): MercatorCoordinates { require(abs(gmc.latitude) <= MAXIMUM_LATITUDE) { "Latitude exceeds the maximum latitude for mercator coordinates" } - val r = correctedRadius?.invoke(gmc) ?: radius + val r: Distance = correctedRadius?.invoke(gmc) ?: radius return MercatorCoordinates( - x = r * (gmc.longitude - baseLongitude), - y = r * ln(tan(PI / 4 + gmc.latitude / 2)) + x = r * (gmc.longitude - baseLongitude).radians.value, + y = r * ln(tan(pi / 4 + gmc.latitude / 2)) ) } - public companion object : MercatorProjection(0.0, 6378137.0) { - public const val MAXIMUM_LATITUDE: Double = 85.05113 - public val DEFAULT_EARTH_RADIUS: Double = radius + public companion object : MercatorProjection(Angle.zero, Distance(6378.137)) { + 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/commonMain/kotlin/center/sciprog/maps/coordinates/WebMercatorProjection.kt b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/WebMercatorProjection.kt index 4273080..d29e4c4 100644 --- a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/WebMercatorProjection.kt +++ b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/WebMercatorProjection.kt @@ -9,7 +9,7 @@ import kotlin.math.* public data class WebMercatorCoordinates(val zoom: Int, val x: Double, val y: Double) -public object WebMercatorProjection { +public object WebMercatorProjection { /** * Compute radians to projection coordinates ratio for given [zoom] factor @@ -32,8 +32,8 @@ public object WebMercatorProjection { val scaleFactor = scaleFactor(zoom.toDouble()) return WebMercatorCoordinates( zoom = zoom, - x = scaleFactor * (gmc.longitude + PI), - y = scaleFactor * (PI - ln(tan(PI / 4 + gmc.latitude / 2))) + x = scaleFactor * (gmc.longitude.radians.value + PI), + y = scaleFactor * (PI - ln(tan(PI / 4 + gmc.latitude.radians.value / 2))) ) } diff --git a/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/angles.kt b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/angles.kt new file mode 100644 index 0000000..c0245ba --- /dev/null +++ b/maps-kt-core/src/commonMain/kotlin/center/sciprog/maps/coordinates/angles.kt @@ -0,0 +1,92 @@ +/* + * 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 center.sciprog.maps.coordinates + +import kotlin.jvm.JvmInline +import kotlin.math.PI + +// Taken from KMath dev version, to be used directly in the future + + +public sealed interface Angle : Comparable { + public val radians: Radians + public val degrees: Degrees + + public operator fun plus(other: Angle): Angle + public operator fun minus(other: Angle): Angle + + public operator fun times(other: Number): Angle + public operator fun div(other: Number): Angle + public operator fun div(other: Angle): Double + public operator fun unaryMinus(): Angle + + public companion object { + public val zero: Angle = 0.radians + public val pi: Angle = PI.radians + public val piTimes2: Angle = (2 * PI).radians + public val piDiv2: Angle = (PI / 2).radians + } +} + +/** + * Type safe radians + */ +@JvmInline +public value class Radians(public val value: Double) : Angle { + override val radians: Radians + get() = this + override val degrees: Degrees + get() = Degrees(value * 180 / PI) + + public override fun plus(other: Angle): Radians = Radians(value + other.radians.value) + public override fun minus(other: Angle): Radians = Radians(value - other.radians.value) + + public override fun times(other: Number): Radians = Radians(value * other.toDouble()) + public override fun div(other: Number): Radians = Radians(value / other.toDouble()) + override fun div(other: Angle): Double = value / other.radians.value + public override fun unaryMinus(): Radians = Radians(-value) + + override fun compareTo(other: Angle): Int = value.compareTo(other.radians.value) +} + +public fun sin(angle: Angle): Double = kotlin.math.sin(angle.radians.value) +public fun cos(angle: Angle): Double = kotlin.math.cos(angle.radians.value) +public fun tan(angle: Angle): Double = kotlin.math.tan(angle.radians.value) + +public val Number.radians: Radians get() = Radians(toDouble()) + +/** + * Type safe degrees + */ +@JvmInline +public value class Degrees(public val value: Double) : Angle { + override val radians: Radians + get() = Radians(value * PI / 180) + override val degrees: Degrees + get() = this + + public override fun plus(other: Angle): Degrees = Degrees(value + other.degrees.value) + public override fun minus(other: Angle): Degrees = Degrees(value - other.degrees.value) + + public override fun times(other: Number): Degrees = Degrees(value * other.toDouble()) + public override fun div(other: Number): Degrees = Degrees(value / other.toDouble()) + override fun div(other: Angle): Double = value / other.degrees.value + public override fun unaryMinus(): Degrees = Degrees(-value) + + override fun compareTo(other: Angle): Int = value.compareTo(other.degrees.value) +} + +public val Number.degrees: Degrees get() = Degrees(toDouble()) + +/** + * Normalized angle to (0, 2PI) for radians or (0, 360) for degrees. + */ +public fun Angle.normalized(): Angle = when (this) { + is Degrees -> (value + 180.0).rem(360.0).degrees + is Radians -> (value + PI).rem(PI * 2).radians +} + +public fun abs(angle: Angle): Angle = if (angle < Angle.zero) -angle else angle \ No newline at end of file diff --git a/maps-kt-core/src/commonTest/kotlin/center/sciprog/maps/coordinates/DistanceTest.kt b/maps-kt-core/src/commonTest/kotlin/center/sciprog/maps/coordinates/DistanceTest.kt new file mode 100644 index 0000000..2e5cf03 --- /dev/null +++ b/maps-kt-core/src/commonTest/kotlin/center/sciprog/maps/coordinates/DistanceTest.kt @@ -0,0 +1,37 @@ +package center.sciprog.maps.coordinates + +import kotlin.test.Ignore +import kotlin.test.Test +import kotlin.test.assertEquals + +internal class DistanceTest { + companion object { + val moscow = GMC.ofDegrees(55.76058287719673, 37.60358622841869) + val spb = GMC.ofDegrees(59.926686023580444, 30.36038109122013) + } + + @Test + fun ellipsoidParameters() { + assertEquals(298.257223563, GeoEllipsoid.WGS84.inverseF, 1e-6) + } + + @Test + @Ignore + fun greatCircleDistance() { + assertEquals( + expected = 632.035, + actual = GeoEllipsoid.greatCircleAngleBetween(moscow, spb).value * + GeoEllipsoid.WGS84.equatorRadius.kilometers, + absoluteTolerance = 0.1 + ) + } + + @Test + fun largeDistance() { + val curve = GeoEllipsoid.WGS84.curveBetween(moscow, spb) + val distance = curve.distance + + assertEquals(632.035426877, distance.kilometers, 0.1) + + } +} \ No newline at end of file