Add support of Gravitational Model offset.

README.md

@@ -25,7 +25,7 @@ repositories {
 }
 
 dependencies {
-    implementation 'earth.worldwind:worldwind:1.6.2'
+    implementation 'earth.worldwind:worldwind:1.6.3'
 }
 ```
 

build.gradle.kts

@@ -16,7 +16,7 @@ buildscript {
 
 allprojects {
     group = "earth.worldwind"
-    version = "1.6.2"
+    version = "1.6.3"
 
     extra.apply {
         set("minSdk", 21)

worldwind/build.gradle.kts

@@ -41,10 +41,11 @@ kotlin {
         val mokoVersion = "0.24.3"
         val ktorVersion = "2.3.12"
         val ormliteVersion = "6.1"
+        val coroutinesVerion = "1.9.0"
         commonMain {
             dependencies {
                 api("org.jetbrains.kotlinx:kotlinx-datetime:0.6.1")
-                api("org.jetbrains.kotlinx:kotlinx-coroutines-core:1.9.0")
+                api("org.jetbrains.kotlinx:kotlinx-coroutines-core:$coroutinesVerion")
                 implementation("org.jetbrains.kotlinx:kotlinx-serialization-json:1.7.3")
                 implementation("io.ktor:ktor-client-core:$ktorVersion")
                 implementation("io.github.pdvrieze.xmlutil:serialization:0.90.3")
@@ -57,6 +58,7 @@ kotlin {
         commonTest {
             dependencies {
                 implementation(kotlin("test"))
+                implementation("org.jetbrains.kotlinx:kotlinx-coroutines-test:$coroutinesVerion")
                 implementation("dev.icerock.moko:resources-test:$mokoVersion")
             }
         }

worldwind/src/androidMain/kotlin/earth/worldwind/globe/geoid/EGM96Geoid.kt

@@ -0,0 +1,24 @@
+package earth.worldwind.globe.geoid
+
+import dev.icerock.moko.resources.AssetResource
+import kotlinx.coroutines.CoroutineScope
+import kotlinx.coroutines.launch
+import java.io.FileNotFoundException
+import java.nio.ByteBuffer
+import java.nio.ByteOrder
+import java.nio.ShortBuffer
+
+private var deltas: ShortBuffer? = null
+actual val isInitialized get() = deltas != null
+
+actual fun loadData(offsetsFile: AssetResource, scope: CoroutineScope) {
+    scope.launch {
+        // NOTE! It is not possible to provide Context here, that's why we are using pure Java resources access
+        val filePath = "assets/${offsetsFile.path}"
+        val bytes = EGM96Geoid.javaClass.classLoader?.getResourceAsStream(filePath)?.use { it.readBytes() }
+            ?: throw FileNotFoundException("Couldn't open resource as stream at: $filePath")
+        deltas = ByteBuffer.wrap(bytes).order(ByteOrder.BIG_ENDIAN).asShortBuffer()
+    }
+}
+
+actual fun getValue(k: Int) = deltas?.get(k) ?: 0

worldwind/src/commonMain/kotlin/earth/worldwind/globe/Globe.kt

@@ -1,7 +1,10 @@
 package earth.worldwind.globe
 
 import earth.worldwind.geom.*
+import earth.worldwind.geom.Angle.Companion.degrees
 import earth.worldwind.globe.elevation.ElevationModel
+import earth.worldwind.globe.geoid.EGM96Geoid
+import earth.worldwind.globe.geoid.Geoid
 import earth.worldwind.globe.projection.GeographicProjection
 import earth.worldwind.globe.projection.Wgs84Projection
 import kotlin.math.PI
@@ -21,7 +24,11 @@ open class Globe(
      * Indicates the geographic projection used by this globe. The projection specifies this globe's Cartesian
      * coordinate system.
      */
-    var projection: GeographicProjection = Wgs84Projection()
+    var projection: GeographicProjection = Wgs84Projection(),
+    /**
+     * Represents Gravitational Model of the globe
+     */
+    var geoid: Geoid = EGM96Geoid()
 ) {
     /**
      * Represents the elevations for an area, often but not necessarily the whole globe.
@@ -184,7 +191,44 @@ open class Globe(
      * @return Elevation in meters in specified location
      */
     fun getElevation(latitude: Angle, longitude: Angle, retrieve: Boolean = false) =
-        elevationModel.getHeight(latitude, longitude, retrieve).toDouble()
+        (elevationModel.getElevation(latitude, longitude, retrieve) + geoid.getOffset(latitude, longitude)).toDouble()
+
+    /**
+     * Gets elevation values for the specified sector with required width and height resolution, including Geoid offset
+     *
+     * @param gridSector specified sector to determine elevation values
+     * @param gridWidth value matrix width
+     * @param gridHeight value matrix height
+     * @param result pre-allocated array for the result. Must be width * height size.
+     */
+    fun getElevationGrid(gridSector: Sector, gridWidth: Int, gridHeight: Int, result: FloatArray) {
+        elevationModel.getElevationGrid(gridSector, gridWidth, gridHeight, result)
+        // Apply Gravitational Model offset
+        val minLat = gridSector.minLatitude.inDegrees
+        val minLon = gridSector.minLongitude.inDegrees
+        val deltaLat = gridSector.deltaLatitude.inDegrees / (gridHeight - 1)
+        val deltaLon = gridSector.deltaLongitude.inDegrees / (gridWidth - 1)
+        var h = 0
+        for (i in 0 until gridHeight) {
+            val lat = minLat + i * deltaLat
+            for (j in 0 until gridWidth) {
+                val lon = minLon + j * deltaLon
+                result[h++] += geoid.getOffset(lat.degrees, lon.degrees)
+            }
+        }
+    }
+
+    /**
+     * Gets elevation limits at specified sector, including Geoid offset
+     *
+     * @param sector specified sector to determine elevation limits
+     * @return pre-allocated array for the result. Must be size of 2.
+     */
+    fun getElevationLimits(sector: Sector, result: FloatArray) {
+        elevationModel.getElevationLimits(sector, result)
+        // Apply Gravitational Model offset
+        for (i in result.indices) result[i] += geoid.getOffset(sector.centroidLatitude, sector.centroidLongitude)
+    }
 
     /**
      * Get absolute position with terrain elevation at specified coordinates

worldwind/src/commonMain/kotlin/earth/worldwind/globe/elevation/ElevationModel.kt

@@ -46,22 +46,22 @@ open class ElevationModel(): Iterable<ElevationCoverage> {
 
     override fun iterator() = coverages.iterator()
 
-    fun getHeight(latitude: Angle, longitude: Angle, retrieve: Boolean): Float {
+    fun getElevation(latitude: Angle, longitude: Angle, retrieve: Boolean): Float {
         // coverages composite from fine to coarse
         for (i in coverages.indices.reversed()) {
-            val height = coverages[i].getHeight(latitude, longitude, retrieve)
+            val height = coverages[i].getElevation(latitude, longitude, retrieve)
             if (height != null) return height
         }
         return 0f
     }
 
-    fun getHeightGrid(gridSector: Sector, gridWidth: Int, gridHeight: Int, result: FloatArray) {
+    fun getElevationGrid(gridSector: Sector, gridWidth: Int, gridHeight: Int, result: FloatArray) {
         // coverages composite from coarse to fine
-        for (i in coverages.indices) coverages[i].getHeightGrid(gridSector, gridWidth, gridHeight, result)
+        for (i in coverages.indices) coverages[i].getElevationGrid(gridSector, gridWidth, gridHeight, result)
     }
 
-    fun getHeightLimits(sector: Sector, result: FloatArray) {
+    fun getElevationLimits(sector: Sector, result: FloatArray) {
         // coverage order is irrelevant
-        for (i in coverages.indices) coverages[i].getHeightLimits(sector, result)
+        for (i in coverages.indices) coverages[i].getElevationLimits(sector, result)
     }
 }

worldwind/src/commonMain/kotlin/earth/worldwind/globe/elevation/coverage/AbstractElevationCoverage.kt

@@ -34,26 +34,26 @@ abstract class AbstractElevationCoverage: ElevationCoverage {
 
     override fun hasUserProperty(key: Any) = userProperties?.containsKey(key) == true
 
-    override fun getHeight(latitude: Angle, longitude: Angle, retrieve: Boolean): Float? {
+    override fun getElevation(latitude: Angle, longitude: Angle, retrieve: Boolean): Float? {
         return if (isEnabled) {
             val key = 31 * latitude.inDegrees.hashCode() + longitude.inDegrees.hashCode()
-            heightCache[key] ?: doGetHeight(latitude, longitude, retrieve)?.also {
+            heightCache[key] ?: doGetElevation(latitude, longitude, retrieve)?.also {
                 heightCache.put(key, it, 1)
             }
         } else null
     }
 
-    override fun getHeightGrid(gridSector: Sector, gridWidth: Int, gridHeight: Int, result: FloatArray) {
-        if (isEnabled) doGetHeightGrid(gridSector, gridWidth, gridHeight, result)
+    override fun getElevationGrid(gridSector: Sector, gridWidth: Int, gridHeight: Int, result: FloatArray) {
+        if (isEnabled) doGetElevationGrid(gridSector, gridWidth, gridHeight, result)
     }
 
-    override fun getHeightLimits(sector: Sector, result: FloatArray) {
-        if (isEnabled) doGetHeightLimits(sector, result)
+    override fun getElevationLimits(sector: Sector, result: FloatArray) {
+        if (isEnabled) doGetElevationLimits(sector, result)
     }
 
-    protected abstract fun doGetHeight(latitude: Angle, longitude: Angle, retrieve: Boolean): Float?
+    protected abstract fun doGetElevation(latitude: Angle, longitude: Angle, retrieve: Boolean): Float?
 
-    protected abstract fun doGetHeightGrid(gridSector: Sector, gridWidth: Int, gridHeight: Int, result: FloatArray)
+    protected abstract fun doGetElevationGrid(gridSector: Sector, gridWidth: Int, gridHeight: Int, result: FloatArray)
 
-    protected abstract fun doGetHeightLimits(sector: Sector, result: FloatArray)
+    protected abstract fun doGetElevationLimits(sector: Sector, result: FloatArray)
 }

worldwind/src/commonMain/kotlin/earth/worldwind/globe/elevation/coverage/AbstractTiledElevationCoverage.kt

@@ -75,7 +75,7 @@ abstract class AbstractTiledElevationCoverage(
         updateTimestamp()
     }
 
-    override fun doGetHeight(latitude: Angle, longitude: Angle, retrieve: Boolean): Float? {
+    override fun doGetElevation(latitude: Angle, longitude: Angle, retrieve: Boolean): Float? {
         if (!tileMatrixSet.sector.contains(latitude, longitude)) return null // no coverage in the specified location
         val targetIdx = tileMatrixSet.entries.size - 1 // retrieve height from last available matrix
         for (idx in targetIdx downTo 0) {
@@ -110,7 +110,7 @@ abstract class AbstractTiledElevationCoverage(
         return null // did not find a tile
     }
 
-    override fun doGetHeightGrid(gridSector: Sector, gridWidth: Int, gridHeight: Int, result: FloatArray) {
+    override fun doGetElevationGrid(gridSector: Sector, gridWidth: Int, gridHeight: Int, result: FloatArray) {
         if (!tileMatrixSet.sector.intersects(gridSector)) return  // no coverage in the specified sector
         val targetPixelSpan = gridSector.deltaLatitude.div(gridHeight)
         val targetIdx = tileMatrixSet.indexOfMatrixNearest(targetPixelSpan)
@@ -126,7 +126,7 @@ abstract class AbstractTiledElevationCoverage(
         }
     }
 
-    override fun doGetHeightLimits(sector: Sector, result: FloatArray) {
+    override fun doGetElevationLimits(sector: Sector, result: FloatArray) {
         if (!tileMatrixSet.sector.intersects(sector)) return  // no coverage in the specified sector
         val targetPixelSpan = sector.deltaLatitude.div(GET_HEIGHT_LIMIT_SAMPLES)
         val targetIdx = tileMatrixSet.indexOfMatrixNearest(targetPixelSpan)

worldwind/src/commonMain/kotlin/earth/worldwind/globe/elevation/coverage/ElevationCoverage.kt

@@ -62,7 +62,7 @@ interface ElevationCoverage {
      * @param retrieve if true, then the value will be retrieved from a remote source for the next frame
      * @return elevation value at specified location
      */
-    fun getHeight(latitude: Angle, longitude: Angle, retrieve: Boolean): Float?
+    fun getElevation(latitude: Angle, longitude: Angle, retrieve: Boolean): Float?
 
     /**
      * Gets elevation values for the specified sector with required width and height resolution
@@ -71,12 +71,12 @@ interface ElevationCoverage {
      * @param gridHeight value matrix height
      * @param result pre-allocated array for the result. Must be width * height size.
      */
-    fun getHeightGrid(gridSector: Sector, gridWidth: Int, gridHeight: Int, result: FloatArray)
+    fun getElevationGrid(gridSector: Sector, gridWidth: Int, gridHeight: Int, result: FloatArray)
 
     /**
      * Gets elevation limits at specified sector
      * @param sector specified sector to determine elevation limits
      * @return pre-allocated array for the result. Must be size of 2.
      */
-    fun getHeightLimits(sector: Sector, result: FloatArray)
+    fun getElevationLimits(sector: Sector, result: FloatArray)
 }

worldwind/src/commonMain/kotlin/earth/worldwind/globe/geoid/EGM96Geoid.kt

@@ -0,0 +1,91 @@
+package earth.worldwind.globe.geoid
+
+import dev.icerock.moko.resources.AssetResource
+import earth.worldwind.MR
+import earth.worldwind.geom.Angle
+import kotlinx.coroutines.CoroutineScope
+import kotlinx.coroutines.DelicateCoroutinesApi
+import kotlinx.coroutines.GlobalScope
+
+expect val isInitialized: Boolean
+expect fun loadData(offsetsFile: AssetResource, scope: CoroutineScope)
+expect fun getValue(k: Int): Short
+
+/**
+ * Computes EGM96 geoid offsets.
+ *
+ * A file with the offset grid must be passed to the constructor. This file must have 721 rows of 1440 2-byte integer
+ * values. Each row corresponding to a latitude, with the first row corresponding to +90 degrees (90 North). The integer
+ * values must be in centimeters.
+ *
+ * @param offsetsFile a resource pointing to a file with the geoid offsets. See the class description above for a
+ * description of the file.
+ */
+@OptIn(DelicateCoroutinesApi::class)
+open class EGM96Geoid(offsetsFile: AssetResource = MR.assets.EGM96_dat, scope: CoroutineScope = GlobalScope) : Geoid {
+    init {
+        loadData(offsetsFile, scope)
+    }
+
+    override fun getOffset(latitude: Angle, longitude: Angle): Float {
+        // Return 0 for all offsets if the file not loaded yet or failed to load.
+        if (!isInitialized) return 0f
+
+        val lat = latitude.inDegrees
+        val lon = if (longitude.inDegrees >= 0.0) longitude.inDegrees else longitude.inDegrees + 360.0
+
+        var topRow = ((90.0 - lat) / INTERVAL).toInt()
+        if (lat <= -90.0) topRow = NUM_ROWS - 2
+        val bottomRow = topRow + 1
+
+        // Note that the number of columns does not repeat the column at 0 longitude, so we must force the right
+        // column to 0 for any longitude that's less than one interval from 360, and force the left column to the
+        // last column of the grid.
+        var leftCol = (lon / INTERVAL).toInt()
+        var rightCol = leftCol + 1
+        if (lon >= 360.0 - INTERVAL) {
+            leftCol = NUM_COLS - 1
+            rightCol = 0
+        }
+
+        val latBottom = 90.0 - bottomRow * INTERVAL
+        val lonLeft = leftCol * INTERVAL
+
+        val ul = getPostOffset(topRow, leftCol)
+        val ll = getPostOffset(bottomRow, leftCol)
+        val lr = getPostOffset(bottomRow, rightCol)
+        val ur = getPostOffset(topRow, rightCol)
+
+        val u = (lon - lonLeft) / INTERVAL
+        val v = (lat - latBottom) / INTERVAL
+
+        val pll = (1.0 - u) * (1.0 - v)
+        val plr = u * (1.0 - v)
+        val pur = u * v
+        val pul = (1.0 - u) * v
+
+        val offset = pll * ll + plr * lr + pur * ur + pul * ul
+
+        return (offset / 100.0).toFloat() // convert centimeters to meters
+    }
+
+    internal fun getPostOffset(row: Int, col: Int) = getValue(row * NUM_COLS + col)
+
+    companion object {
+        // Description of the EGMA96 offsets file:
+        // See: http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/binary/binarygeoid.html
+        //    The total size of the file is 2,076,480 bytes. This file was created
+        //    using an INTEGER*2 data type format and is an unformatted direct access
+        //    file. The data on the file is arranged in records from north to south.
+        //    There are 721 records on the file starting with record 1 at 90 N. The
+        //    last record on the file is at latitude 90 S. For each record, there
+        //    are 1,440 15 arc-minute geoid heights arranged by longitude from west to
+        //    east starting at the Prime Meridian (0 E) and ending 15 arc-minutes west
+        //    of the Prime Meridian (359.75 E). On file, the geoid heights are in units
+        //    of centimeters. While retrieving the Integer*2 values on file, divide by
+        //    100 and this will produce a geoid height in meters.
+        internal const val INTERVAL = 15.0 / 60.0 // 15' angle delta
+        internal const val NUM_ROWS = 721
+        internal const val NUM_COLS = 1440
+    }
+}

worldwind/src/commonMain/kotlin/earth/worldwind/globe/geoid/Geoid.kt

@@ -0,0 +1,17 @@
+package earth.worldwind.globe.geoid
+
+import earth.worldwind.geom.Angle
+
+/**
+ * Representation of the globe Gravitational Model
+ */
+interface Geoid {
+    /**
+     * Calculates Gravitational Model offset at specified location
+     *
+     * @param latitude Input latitude
+     * @param longitude Input longitude
+     * @return Gravitational Model offset
+     */
+    fun getOffset(latitude: Angle, longitude: Angle): Float
+}

worldwind/src/commonMain/kotlin/earth/worldwind/globe/terrain/TerrainTile.kt

@@ -40,7 +40,7 @@ open class TerrainTile(sector: Sector, level: Level, row: Int, column: Int): Til
         val timestamp = rc.elevationModelTimestamp
         if (timestamp != heightTimestamp) {
             heightGrid.fill(0f)
-            globe.elevationModel.getHeightGrid(sector, tileWidth, tileHeight, heightGrid)
+            globe.getElevationGrid(sector, tileWidth, tileHeight, heightGrid)
             // Calculate height vertex buffer from height grid
             for (r in 0 until level.tileHeight) for (c in 0 until level.tileWidth) {
                 heights[(r + 1) * (level.tileWidth + 2) + c + 1] = heightGrid[r * level.tileWidth + c]

worldwind/src/commonMain/kotlin/earth/worldwind/render/AbstractSurfaceRenderable.kt

@@ -41,7 +41,7 @@ abstract class AbstractSurfaceRenderable(sector: Sector, displayName: String? =
         // initialize the heights for elevation model scan
         heightLimits[0] = Float.MAX_VALUE
         heightLimits[1] = -Float.MAX_VALUE
-        globe.elevationModel.getHeightLimits(sector, heightLimits)
+        globe.getElevationLimits(sector, heightLimits)
         // check for valid height limits
         if (heightLimits[0] > heightLimits[1]) heightLimits.fill(0f)
     }

worldwind/src/commonMain/kotlin/earth/worldwind/util/AbstractTile.kt

@@ -111,7 +111,7 @@ abstract class AbstractTile(
         // initialize the heights for elevation model scan
         heightLimits[0] = Float.MAX_VALUE
         heightLimits[1] = -Float.MAX_VALUE
-        globe.elevationModel.getHeightLimits(sector, heightLimits)
+        globe.getElevationLimits(sector, heightLimits)
         // check for valid height limits
         if (heightLimits[0] > heightLimits[1]) heightLimits.fill(0f)
     }