Implement a class to load and sample an EGM96 grid

CHANGES.md

@@ -11,6 +11,7 @@
 - Added full support for custom ellipsoids by setting `TilesetOptions::ellipsoid` when creating a tileset.
   - Many methods have been updated with an additional ellipsoid parameter to support this. The WGS84 ellipsoid is used as a default parameter here to ensure API compatibility.
   - `CESIUM_DISABLE_DEFAULT_ELLIPSOID` can be defined to disable the WGS84 default parameter, exposing through errors the places in your code that are still assuming a WGS84 ellipsoid.
+- Added `CesiumGeospatial::EarthGravitationalModel96Grid` class to allow transforming heights on a WGS84 ellipsoid into heights above mean sea level using the EGM96 model.
 
 ### v0.36.0 - 2024-06-03
 

CesiumGeospatial/CMakeLists.txt

@@ -12,6 +12,7 @@ set_target_properties(CesiumGeospatial
     PROPERTIES
         TEST_SOURCES "${CESIUM_GEOSPATIAL_TEST_SOURCES}"
         TEST_HEADERS "${CESIUM_GEOSPATIAL_TEST_HEADERS}"
+        TEST_DATA_DIR ${CMAKE_CURRENT_LIST_DIR}/test/data
 )
 
 set_target_properties(CesiumGeospatial

CesiumGeospatial/include/CesiumGeospatial/EarthGravitationalModel96Grid.h

@@ -0,0 +1,66 @@
+#pragma once
+
+#include "Library.h"
+
+#include <CesiumGeospatial/Cartographic.h>
+
+#include <gsl/span>
+
+#include <cstdint>
+#include <optional>
+#include <string>
+#include <vector>
+
+namespace CesiumGeospatial {
+
+/**
+ * @brief Loads and queries heights from an Earth Gravitational Model 96 (EGM96)
+ * grid.
+ *
+ * EGM96 is a standard geopotential model of the earth's surface, which can be
+ * used to obtain an approximation of the mean sea level (MSL) at any location
+ * on a WGS84 ellipsoid.
+ */
+class CESIUMGEOSPATIAL_API EarthGravitationalModel96Grid final {
+public:
+  /**
+   * @brief Attempts to create a {@link EarthGravitationalModel96Grid} from the given file.
+   *
+   * This method expects a file in the format of the WW15MGH.DAC 15-arcminute
+   * grid.
+   */
+  static std::optional<EarthGravitationalModel96Grid>
+  fromFile(const std::string& filename);
+
+  /**
+   * Attempts to create a {@link EarthGravitationalModel96Grid} from the given buffer.
+   *
+   * This method expects the buffer to contain the contents of the WW15MGH.DAC
+   * 15-arcminute grid.
+   */
+  static std::optional<EarthGravitationalModel96Grid>
+  fromBuffer(const gsl::span<const std::byte>& buffer);
+
+  /**
+   * @brief Samples the height at the given position.
+   *
+   * @param position The position to sample.
+   * @returns The height (in meters) of the EGM96 surface above the WGS84
+   * ellipsoid. A positive value indicates that MSL is above the ellipsoid's
+   * surface, while a negative value indicates that MSL is below the ellipsoid's
+   * surface.
+   */
+  double sampleHeight(const Cartographic& position) const;
+
+private:
+  EarthGravitationalModel96Grid(std::vector<int16_t>&& gridValues);
+
+  /**
+   * Returns the height of the given grid value, in meters.
+   */
+  double getHeightForIndices(const int vertical, const int horizontal) const;
+
+  std::vector<int16_t> _gridValues;
+};
+
+} // namespace CesiumGeospatial

CesiumGeospatial/src/EarthGravitationalModel96Grid.cpp

@@ -0,0 +1,120 @@
+#include <CesiumGeospatial/EarthGravitationalModel96Grid.h>
+#include <CesiumGeospatial/Ellipsoid.h>
+#include <CesiumUtility/Math.h>
+
+#include <algorithm>
+#include <fstream>
+#include <vector>
+
+using namespace CesiumUtility;
+
+namespace CesiumGeospatial {
+
+// WW15MGH.DAC has 721 rows representing the range (90N, 90S) and 1440 columns
+// representing the range (0E, 360E)
+
+// The number of rows in the file
+const int NUM_ROWS = 721;
+// The number of columns in the file
+const int NUM_COLUMNS = 1440;
+
+std::optional<EarthGravitationalModel96Grid>
+CesiumGeospatial::EarthGravitationalModel96Grid::fromFile(
+    const std::string& filename) {
+  std::ifstream file(filename, std::ios::binary | std::ios::ate);
+  if (!file.good()) {
+    return std::nullopt;
+  }
+
+  std::streamsize size = file.tellg();
+  file.seekg(0, std::ios::beg);
+
+  std::vector<std::byte> buffer(static_cast<size_t>(size));
+  file.read(reinterpret_cast<char*>(buffer.data()), size);
+  return fromBuffer(buffer);
+}
+
+std::optional<EarthGravitationalModel96Grid>
+CesiumGeospatial::EarthGravitationalModel96Grid::fromBuffer(
+    const gsl::span<const std::byte>& buffer) {
+  const size_t expectedBytes = NUM_ROWS * NUM_COLUMNS * sizeof(int16_t);
+
+  if (
+      // Not enough data - is this a valid WW15MGH.DAC?
+      buffer.size_bytes() < expectedBytes ||
+      // WW15MGH.DAC should just be full of 2-byte int16_t, so if it's not even
+      // then something has gone wrong (so we should stop now before we overflow
+      // the end of the gridValues buffer!)
+      buffer.size_bytes() % 2 != 0) {
+    return std::nullopt;
+  }
+
+  const size_t size = buffer.size_bytes();
+  const std::byte* pRead = buffer.data();
+  std::vector<int16_t> gridValues;
+  gridValues.resize(size / 2);
+  std::byte* pWrite = reinterpret_cast<std::byte*>(gridValues.data());
+
+  for (size_t i = 0; i < size; i += 2) {
+    // WW15MGH.DAC is in big endian, so we swap the bytes
+    pWrite[i] = pRead[i + 1];
+    pWrite[i + 1] = pRead[i];
+  }
+
+  return EarthGravitationalModel96Grid(std::move(gridValues));
+}
+
+double EarthGravitationalModel96Grid::sampleHeight(
+    const Cartographic& position) const {
+  const double longitude = Math::zeroToTwoPi(position.longitude);
+  const double latitude =
+      Math::clamp(position.latitude, -Math::PiOverTwo, Math::PiOverTwo);
+
+  const double horizontalIndexDecimal = (NUM_COLUMNS * longitude) / Math::TwoPi;
+  const int horizontalIndex = static_cast<int>(horizontalIndexDecimal);
+
+  const double verticalIndexDecimal =
+      ((NUM_ROWS - 1) * (Math::PiOverTwo - latitude)) / Math::OnePi;
+  const int verticalIndex = static_cast<int>(verticalIndexDecimal);
+
+  // Get the normalized position of the coordinates within the grid tile
+  const double xn = horizontalIndexDecimal - horizontalIndex;
+  const double iyn = verticalIndexDecimal - verticalIndex;
+  const double ixn = 1.0 - xn;
+  const double yn = 1.0 - iyn;
+
+  // Get surrounding grid points
+  const double a1 = getHeightForIndices(horizontalIndex, verticalIndex);
+  const double a2 = getHeightForIndices(horizontalIndex + 1, verticalIndex);
+  const double b1 = getHeightForIndices(horizontalIndex, verticalIndex + 1);
+  const double b2 = getHeightForIndices(horizontalIndex + 1, verticalIndex + 1);
+
+  // Each component contributes to the total based on its normalized position in
+  // the grid square
+  const double result =
+      (a1 * ixn * yn + a2 * xn * yn + b1 * ixn * iyn + b2 * xn * iyn);
+
+  return result;
+}
+
+EarthGravitationalModel96Grid::EarthGravitationalModel96Grid(
+    std::vector<int16_t>&& gridValues)
+    : _gridValues(gridValues) {}
+
+double EarthGravitationalModel96Grid::getHeightForIndices(
+    const int horizontal,
+    const int vertical) const {
+
+  int clampedVertical = vertical;
+  if (vertical > NUM_ROWS - 1) {
+    clampedVertical = NUM_ROWS - 1;
+  }
+
+  const size_t index =
+      static_cast<size_t>(clampedVertical * NUM_COLUMNS + horizontal);
+  const double result = _gridValues[index] / 100.0;
+
+  return result;
+}
+
+} // namespace CesiumGeospatial