diff --git a/turf/src/commonMain/kotlin/io/github/dellisd/spatialk/turf/Measurement.kt b/turf/src/commonMain/kotlin/io/github/dellisd/spatialk/turf/Measurement.kt index 2ce3f8f..817099a 100644 --- a/turf/src/commonMain/kotlin/io/github/dellisd/spatialk/turf/Measurement.kt +++ b/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, + antimeridianOffset: Double + ): List> { + 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>() + if (passesAntimeridian && maxSmallDiffLong < antimeridianOffset) { + var poNewLS = mutableListOf() + 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.. 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()) + ) + } +} + diff --git a/turf/src/commonMain/kotlin/io/github/dellisd/spatialk/turf/Units.kt b/turf/src/commonMain/kotlin/io/github/dellisd/spatialk/turf/Units.kt index 4d2cf7c..1b08349 100644 --- a/turf/src/commonMain/kotlin/io/github/dellisd/spatialk/turf/Units.kt +++ b/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. * diff --git a/turf/src/commonTest/kotlin/io/github/dellisd/spatialk/turf/TurfMeasurementTest.kt b/turf/src/commonTest/kotlin/io/github/dellisd/spatialk/turf/TurfMeasurementTest.kt index d07f785..d17dd88 100644 --- a/turf/src/commonTest/kotlin/io/github/dellisd/spatialk/turf/TurfMeasurementTest.kt +++ b/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(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(geometry) + } + + @Test + fun testGreatCircleAntipodals() { + val startLat = 47.620422 + + val start = Position(-122.349358, startLat) + val antipodal = Position(106.33, startLat) + + assertFails { + greatCircle(start, antipodal) + } + } }