Compare commits

...

5 Commits

4 changed files with 216 additions and 1 deletions

View File

@ -73,4 +73,4 @@ ksciencePublish {
sonatype() sonatype()
} }
apiValidation.nonPublicMarkers.add("space.kscience.kmath.misc.UnstableKMathAPI") apiValidation.nonPublicMarkers.add("space.kscience.kmath.misc.UnstableKMathAPI")

View File

@ -0,0 +1,97 @@
/*
* 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 space.kscience.kmath.geometry
import space.kscience.kmath.operations.Algebra
import kotlin.math.*
/**
* Geodetic coordinated
*/
public class GeodeticCoordinates private constructor(public val latitude: Double, public val longitude: Double) {
override fun equals(other: Any?): Boolean {
if (this === other) return true
if (other == null || this::class != other::class) return false
other as GeodeticCoordinates
if (latitude != other.latitude) return false
if (longitude != other.longitude) return false
return true
}
override fun hashCode(): Int {
var result = latitude.hashCode()
result = 31 * result + longitude.hashCode()
return result
}
override fun toString(): String {
return "GeodeticCoordinates(latitude=$latitude, longitude=$longitude)"
}
public companion object {
public fun ofRadians(latitude: Double, longitude: Double): GeodeticCoordinates {
require(longitude in (-PI)..(PI)) { "Longitude $longitude is not in (-PI)..(PI)" }
return GeodeticCoordinates(latitude, longitude.rem(PI / 2))
}
public fun ofDegrees(latitude: Double, longitude: Double): GeodeticCoordinates {
require(latitude in (-90.0)..(90.0)) { "Latitude $latitude is not in -90..90" }
return GeodeticCoordinates(latitude * PI / 180, (longitude.rem(180) * PI / 180))
}
}
}
public data class MercatorCoordinates(val x: Double, val y: Double)
/**
* @param baseLongitude the longitude offset in radians
* @param radius the average radius of the Earth
* @param correctedRadius optional radius correction to account for ellipsoid model
*/
public open class MercatorAlgebra(
public val baseLongitude: Double = 0.0,
protected val radius: Double = DEFAULT_EARTH_RADIUS,
private val correctedRadius: ((GeodeticCoordinates) -> Double)? = null,
) : Algebra<TileWebMercatorCoordinates> {
public fun MercatorCoordinates.toGeodetic(): GeodeticCoordinates {
val res = GeodeticCoordinates.ofRadians(
atan(sinh(y / radius)),
baseLongitude + x / radius,
)
return if (correctedRadius != null) {
val r = correctedRadius.invoke(res)
GeodeticCoordinates.ofRadians(
atan(sinh(y / r)),
baseLongitude + x / r,
)
} else {
res
}
}
/**
* https://en.wikipedia.org/wiki/Web_Mercator_projection#Formulas
*/
public fun GeodeticCoordinates.toMercator(): MercatorCoordinates {
require(abs(latitude) <= MAXIMUM_LATITUDE) { "Latitude exceeds the maximum latitude for mercator coordinates" }
val r = correctedRadius?.invoke(this) ?: radius
return MercatorCoordinates(
x = r * (longitude - baseLongitude),
y = r * ln(tan(PI / 4 + latitude / 2))
)
}
public companion object : MercatorAlgebra(0.0, 6378137.0) {
public const val MAXIMUM_LATITUDE: Double = 85.05113
public val DEFAULT_EARTH_RADIUS: Double = radius
}
}

View File

@ -0,0 +1,60 @@
/*
* 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 space.kscience.kmath.geometry
import space.kscience.kmath.operations.Algebra
import kotlin.math.*
public data class TileWebMercatorCoordinates(val zoom: Double, val x: Double, val y: Double)
public object WebMercatorAlgebra : Algebra<TileWebMercatorCoordinates> {
private fun scaleFactor(zoom: Double) = 256.0 / 2 / PI * 2.0.pow(zoom)
public fun TileWebMercatorCoordinates.toGeodetic(): GeodeticCoordinates {
val scaleFactor = scaleFactor(zoom)
val longitude = x / scaleFactor - PI
val latitude = (atan(exp(PI - y / scaleFactor)) - PI / 4) * 2
return GeodeticCoordinates.ofRadians(latitude, longitude)
}
/**
* https://en.wikipedia.org/wiki/Web_Mercator_projection#Formulas
*/
public fun GeodeticCoordinates.toMercator(zoom: Double): TileWebMercatorCoordinates {
require(abs(latitude) <= MercatorAlgebra.MAXIMUM_LATITUDE) { "Latitude exceeds the maximum latitude for mercator coordinates" }
val scaleFactor = scaleFactor(zoom)
return TileWebMercatorCoordinates(
zoom = zoom,
x = scaleFactor * (longitude + PI),
y = scaleFactor * (PI - ln(tan(PI / 4 + latitude / 2)))
)
}
/**
* Compute and offset of [target] coordinate relative to [base] coordinate. If [zoom] is null, then optimal zoom
* will be computed to put the resulting x and y coordinates between -127.0 and 128.0
*/
public fun computeOffset(
base: GeodeticCoordinates,
target: GeodeticCoordinates,
zoom: Double? = null,
): TileWebMercatorCoordinates {
val xOffsetUnscaled = target.longitude - base.longitude
val yOffsetUnscaled = ln(
tan(PI / 4 + target.latitude / 2) / tan(PI / 4 + base.latitude / 2)
)
val computedZoom = zoom ?: floor(log2(PI / max(abs(xOffsetUnscaled), abs(yOffsetUnscaled))))
val scaleFactor = scaleFactor(computedZoom)
return TileWebMercatorCoordinates(
computedZoom,
x = scaleFactor * xOffsetUnscaled,
y = scaleFactor * yOffsetUnscaled
)
}
}

View File

@ -0,0 +1,58 @@
/*
* 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 space.kscience.kmath.geometry
import kotlin.test.Test
import kotlin.test.assertEquals
import kotlin.test.assertTrue
class MercatorTest {
@Test
fun boundaries() {
assertEquals(GeodeticCoordinates.ofDegrees(45.0, 20.0), GeodeticCoordinates.ofDegrees(45.0, 200.0))
}
@Test
fun mercatorConversion() = with(MercatorAlgebra) {
val mskCoordinates = GeodeticCoordinates.ofDegrees(55.7558, 37.6173)
val m = mskCoordinates.toMercator()
val r = m.toGeodetic()
assertEquals(mskCoordinates.longitude, r.longitude, 1e-4)
assertEquals(mskCoordinates.latitude, r.latitude, 1e-4)
}
@Test
fun webMercatorConversion() = with(WebMercatorAlgebra) {
val mskCoordinates = GeodeticCoordinates.ofDegrees(55.7558, 37.6173)
val m = mskCoordinates.toMercator(2.0)
val r = m.toGeodetic()
assertEquals(mskCoordinates.longitude, r.longitude, 1e-4)
assertEquals(mskCoordinates.latitude, r.latitude, 1e-4)
}
@Test
fun webMercatorOffset() {
val mskCoordinates = GeodeticCoordinates.ofDegrees(55.7558, 37.6173)
val spbCoordinates = GeodeticCoordinates.ofDegrees(59.9311, 30.3609)
val offset = WebMercatorAlgebra.computeOffset(mskCoordinates, spbCoordinates)
assertTrue { offset.x in -127.0..128.0 }
assertTrue { offset.y in -127.0..128.0 }
assertTrue { offset.zoom > 0.0 }
}
@Test
fun webMercatorAbsolute() {
val mskCoordinates = GeodeticCoordinates.ofDegrees(55.7558, 37.6173)
val wmc = with(WebMercatorAlgebra) { mskCoordinates.toMercator(13.0) }
assertEquals(wmc.x, 1267712.0, 50.0)
assertEquals(wmc.y, 655616.0, 50.0)
}
}