Type safe angles and distances
This commit is contained in:
parent
11b278fc81
commit
4614b1f7bc
@ -1,4 +0,0 @@
|
||||
import kotlin.math.PI
|
||||
|
||||
fun Double.toDegrees() = this * 180 / PI
|
||||
|
@ -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<Double, Double> 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
|
||||
}
|
||||
|
@ -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<MapFeature>.computeBoundingBox(zoom: Int): GmcBox? =
|
||||
public fun Iterable<MapFeature>.computeBoundingBox(zoom: Int): GmcRectangle? =
|
||||
mapNotNull { it.getBoundingBox(zoom) }.wrapAll()
|
||||
|
||||
internal fun Pair<Double, Double>.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<FeatureId, MapFeature>,
|
||||
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()
|
||||
}
|
@ -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,
|
||||
|
@ -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 ->
|
||||
internal fun GmcRectangle.computeViewPoint(mapTileProvider: MapTileProvider): (canvasSize: DpSize) -> MapViewPoint =
|
||||
{ canvasSize ->
|
||||
val zoom = log2(
|
||||
min(
|
||||
canvasSize.width.value / width,
|
||||
canvasSize.height.value / height
|
||||
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<FeatureId, MapFeature> = emptyMap(),
|
||||
config: MapViewConfig = MapViewConfig(),
|
||||
modifier: Modifier = Modifier.fillMaxSize(),
|
||||
|
@ -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(
|
||||
|
@ -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}" }
|
||||
}
|
||||
|
@ -15,4 +15,12 @@ kotlin {
|
||||
js(IR) {
|
||||
browser()
|
||||
}
|
||||
|
||||
sourceSets{
|
||||
commonTest{
|
||||
dependencies{
|
||||
implementation(kotlin("test"))
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
@ -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<Distance> {
|
||||
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
|
@ -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)
|
@ -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))
|
||||
}
|
@ -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
|
||||
|
@ -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<GmcBox>.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))
|
||||
}
|
@ -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
|
||||
)
|
||||
}
|
@ -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())
|
||||
|
@ -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<GmcRectangle>.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))
|
||||
}
|
@ -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
|
||||
)
|
||||
|
@ -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
|
||||
}
|
||||
}
|
@ -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)))
|
||||
)
|
||||
}
|
||||
|
||||
|
@ -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<Angle> {
|
||||
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
|
@ -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)
|
||||
|
||||
}
|
||||
}
|
Loading…
Reference in New Issue
Block a user