Measurement - greatCircle (#131)

* Add greatCircle measurement function.

* Fix several detekt issues

turf/src/commonMain/kotlin/io/github/dellisd/spatialk/turf/Measurement.kt

@@ -488,3 +488,167 @@ fun center(feature: Feature): Point {
 fun center(geometry: Geometry): Point {
     return center(Feature(geometry = geometry))
 }
+
+/**
+ * Calculate great circles routes as [LineString]. Raises error when [start] and [end] are antipodes.
+ *
+ * @param start source position
+ * @param end destination position
+ * @param pointCount number of positions on the arc (including [start] and [end])
+ * @param antimeridianOffset from antimeridian in degrees (default long. = +/- 10deg, geometries within 170deg to
+ * -170deg will be split)
+ *
+ */
+@Suppress("CyclomaticComplexMethod")
+@Throws(IllegalArgumentException::class)
+@ExperimentalTurfApi
+fun greatCircle(start: Position, end: Position, pointCount: Int = 100, antimeridianOffset: Double = 10.0): Geometry {
+
+    val deltaLongitude = start.longitude - end.longitude
+    val deltaLatitude = start.latitude - end.latitude
+
+    // check antipodal positions
+    @Suppress("MagicNumber")
+    require(abs(deltaLatitude) != 0.0 && abs(deltaLongitude % 360) - ANTIMERIDIAN_POS != 0.0) {
+        "Input $start and $end are diametrically opposite, thus there is no single route but rather infinite"
+    }
+
+    val distance = distance(start, end, Units.Radians)
+
+    /**
+     * Calculates the intermediate point on a great circle line
+     *         http://www.edwilliams.org/avform.htm#Intermediate
+     */
+    fun intermediateCoordinate(fraction: Double): Position {
+        val lon1 = radians(start.longitude)
+        val lon2 = radians(end.longitude)
+        val lat1 = radians(start.latitude)
+        val lat2 = radians(end.latitude)
+
+        val a = sin((1 - fraction) * distance) / sin(distance)
+        val b = sin(fraction * distance) / sin(distance)
+        val x = a * cos(lat1) * cos(lon1) + b * cos(lat2) * cos(lon2)
+        val y = a * cos(lat1) * sin(lon1) + b * cos(lat2) * sin(lon2)
+        val z = a * sin(lat1) + b * sin(lat2)
+
+        val lat = degrees(atan2(z, sqrt(x.pow(2) + y.pow(2))))
+        val lon = degrees(atan2(y, x))
+        return Position(lon, lat)
+    }
+
+    @Suppress("LongMethod")
+    fun createCoordinatesAntimeridianAttended(
+        plainArc: List<Position>,
+        antimeridianOffset: Double
+    ): List<List<Position>> {
+        val borderEast = ANTIMERIDIAN_POS - antimeridianOffset
+        val borderWest = ANTIMERIDIAN_NEG + antimeridianOffset
+
+        @Suppress("MagicNumber")
+        val diffSpace = 360.0 - antimeridianOffset
+
+        val passesAntimeridian = plainArc.zipWithNext { a, b ->
+            val diff = abs(a.longitude - b.longitude)
+            (diff > diffSpace &&
+                    ((a.longitude > borderEast && b.longitude < borderWest) ||
+                            (b.longitude > borderEast && a.longitude < borderWest)))
+        }.any()
+
+        val maxSmallDiffLong = plainArc.zipWithNext { a, b -> abs(a.longitude - b.longitude) }
+            .filter { it <= diffSpace } // Filter differences less than or equal to diffSpace
+            .maxByOrNull { it }?.toDouble() ?: 0.0
+
+        val poMulti = mutableListOf<List<Position>>()
+        if (passesAntimeridian && maxSmallDiffLong < antimeridianOffset) {
+            var poNewLS = mutableListOf<Position>()
+            plainArc.forEachIndexed { k, currentPosition ->
+                if (k > 0 && abs(currentPosition.longitude - plainArc[k - 1].longitude) > diffSpace) {
+                    val previousPosition = plainArc[k - 1]
+                    var lon1 = previousPosition.longitude
+                    var lat1 = previousPosition.latitude
+                    var lon2 = currentPosition.longitude
+                    var lat2 = currentPosition.latitude
+
+                    @Suppress("ComplexCondition")
+                    if (lon1 in (ANTIMERIDIAN_NEG + 1..<borderWest) &&
+                        lon2 == ANTIMERIDIAN_POS &&
+                        k + 1 < plainArc.size
+                    ) {
+                        poNewLS.add(Position(ANTIMERIDIAN_NEG, currentPosition.latitude))
+                        poNewLS.add(Position(plainArc[k + 1].longitude, plainArc[k + 1].latitude))
+                        return@forEachIndexed
+                    } else if (
+                        lon1 > borderEast &&
+                        lon1 < ANTIMERIDIAN_POS &&
+                        lon2 == ANTIMERIDIAN_POS &&
+                        k + 1 < plainArc.size
+                    ) {
+                        poNewLS.add(Position(ANTIMERIDIAN_POS, currentPosition.latitude))
+                        poNewLS.add(Position(plainArc[k + 1].longitude, plainArc[k + 1].latitude))
+                        return@forEachIndexed
+                    }
+
+                    if (lon1 < borderWest && lon2 > borderEast) {
+                        val tmpX = lon1
+                        lon1 = lon2
+                        lon2 = tmpX
+                        val tmpY = lat1
+                        lat1 = lat2
+                        lat2 = tmpY
+                    }
+                    if (lon1 > borderEast && lon2 < borderWest) {
+                        @Suppress("MagicNumber")
+                        lon2 += 360.0
+                    }
+
+                    if (ANTIMERIDIAN_POS in lon1..lon2 && lon1 < lon2) {
+                        val ratio = (ANTIMERIDIAN_POS - lon1) / (lon2 - lon1)
+                        val lat = ratio * lat2 + (1 - ratio) * lat1
+                        poNewLS.add(
+                            if (previousPosition.longitude > borderEast) Position(ANTIMERIDIAN_POS, lat)
+                            else Position(ANTIMERIDIAN_NEG, lat)
+                        )
+                        poMulti.add(poNewLS.toList())
+                        poNewLS = mutableListOf() // Clear poNewLS instead of replacing it with an empty list
+                        poNewLS.add(
+                            if (previousPosition.longitude > borderEast) Position(ANTIMERIDIAN_NEG, lat)
+                            else Position(ANTIMERIDIAN_POS, lat)
+                        )
+                    } else {
+                        poNewLS = mutableListOf() // Clear poNewLS instead of replacing it with an empty list
+                        poMulti.add(poNewLS.toList())
+                    }
+                }
+                poNewLS.add(currentPosition) // Adding current position to poNewLS
+            }
+            poMulti.add(poNewLS.toList()) // Adding the last remaining poNewLS to poMulti
+        } else {
+            poMulti.add(plainArc)
+        }
+        return poMulti
+    }
+
+    val arc = buildList {
+        add(start)
+        (1 until (pointCount - 1)).forEach { i ->
+            add(
+                intermediateCoordinate((i + 1).toDouble() / (pointCount - 2 + 1))
+            )
+        }
+        add(end)
+    }
+
+    val coordinates = createCoordinatesAntimeridianAttended(arc, antimeridianOffset)
+    return if (coordinates.size == 1) {
+        LineString(
+            coordinates = coordinates[0],
+            bbox = computeBbox(coordinates[0])
+        )
+    } else {
+        MultiLineString(
+            coordinates = coordinates,
+            bbox = computeBbox(coordinates.flatten())
+        )
+    }
+}
+

turf/src/commonMain/kotlin/io/github/dellisd/spatialk/turf/Units.kt

@@ -8,6 +8,9 @@ import kotlin.native.concurrent.SharedImmutable
 @SharedImmutable
 const val EARTH_RADIUS = 6371008.8
 
+internal const val ANTIMERIDIAN_POS = 180.0
+internal const val ANTIMERIDIAN_NEG = -180.0
+
 /**
  * Supported units of measurement in Turf.
  *

turf/src/commonTest/kotlin/io/github/dellisd/spatialk/turf/TurfMeasurementTest.kt

@@ -1,8 +1,10 @@
 @file:Suppress("MagicNumber")
 
 package io.github.dellisd.spatialk.turf
+
 import io.github.dellisd.spatialk.geojson.BoundingBox
 import io.github.dellisd.spatialk.geojson.LineString
+import io.github.dellisd.spatialk.geojson.MultiLineString
 import io.github.dellisd.spatialk.geojson.Feature
 import io.github.dellisd.spatialk.geojson.Point
 import io.github.dellisd.spatialk.geojson.Polygon
@@ -13,6 +15,8 @@ import io.github.dellisd.spatialk.turf.utils.assertDoubleEquals
 import io.github.dellisd.spatialk.turf.utils.readResource
 import kotlin.test.Test
 import kotlin.test.assertEquals
+import kotlin.test.assertFails
+import kotlin.test.assertIs
 
 @ExperimentalTurfApi
 class TurfMeasurementTest {
@@ -141,4 +145,45 @@ class TurfMeasurementTest {
         assertDoubleEquals(-75.71805238723755, centerPoint.coordinates.longitude, 0.0001)
         assertDoubleEquals(45.3811030151199, centerPoint.coordinates.latitude, 0.0001)
     }
+
+    @Test
+    fun testGreatCircleBasic() {
+        val startLon = -122.0
+        val startLat = 48.0
+        val endLon = 28.125
+        val endLat = 43.32517767999296
+
+        val start = Position(startLon, startLat)
+        val end = Position(endLon, endLat)
+
+        val greatCircle = greatCircle(start, end, pointCount = 99)
+
+        assertEquals(99, greatCircle.coordAll().size)
+        assertIs<LineString>(greatCircle)
+    }
+
+    @Test
+    fun testGreatCirclePassesAntimeridian() {
+        val startLon = -122.349358
+        val startLat = 47.620422
+        val endLon = 77.036560
+        val endLat = 38.897957
+
+        val start = Position(startLon, startLat)
+        val end = Position(endLon, endLat)
+        val geometry = greatCircle(start, end, pointCount = 100)
+        assertIs<MultiLineString>(geometry)
+    }
+
+    @Test
+    fun testGreatCircleAntipodals() {
+        val startLat = 47.620422
+
+        val start = Position(-122.349358, startLat)
+        val antipodal = Position(106.33, startLat)
+
+        assertFails {
+            greatCircle(start, antipodal)
+        }
+    }
 }