feature(Geoid): add support for geoid heights

docs/config.json

@@ -29,13 +29,15 @@
                 "Extent",
                 "OrientationUtils",
                 "OBB",
-                "DEMUtils"
+                "DEMUtils",
+                "GeoidGrid"
             ],
 
             "Layer": [
                 "Layer",
                 "ColorLayer",
                 "ElevationLayer",
+                "GeoidLayer",
                 "GeometryLayer",
                 "FeatureGeometryLayer",
                 "TiledGeometryLayer",

examples/js/GUI/GuiTools.js

@@ -44,8 +44,10 @@ function GuiTools(domId, view, w) {
         document.body.appendChild(element);
         this.colorGui = this.gui.addFolder('Color Layers');
         this.elevationGui = this.gui.addFolder('Elevation Layers');
+        this.geoidGui = this.gui.addFolder('Geoid Layers');
         this.elevationGui.hide();
         this.colorGui.hide();
+        this.geoidGui.hide();
         this.view = view;
         view.addEventListener('layers-order-changed', (function refreshColorGui() {
             var i;
@@ -64,6 +66,8 @@ GuiTools.prototype.addLayerGUI = function fnAddLayerGUI(layer) {
         this.addImageryLayerGUI(layer);
     } else if (layer.isElevationLayer) {
         this.addElevationLayerGUI(layer);
+    } else if (layer.isGeoidLayer) {
+        this.addGeoidLayerGUI(layer);
     }
 };
 
@@ -107,6 +111,19 @@ GuiTools.prototype.addElevationLayerGUI = function addElevationLayerGUI(layer) {
     }).bind(this));
 };
 
+GuiTools.prototype.addGeoidLayerGUI = function addGeoidLayerGUI(layer) {
+    if (this.geoidGui.hasFolder(layer.id)) { return; }
+    this.geoidGui.show();
+    var folder = this.geoidGui.addFolder(layer.id);
+    folder.add({ frozen: layer.frozen }, 'frozen').onChange(function refreshFrozenGui(value) {
+        layer.frozen = value;
+    });
+    folder.add({ visible: layer.visible }, 'visible').onChange((function updateVisibility(value) {
+        layer.visible = value;
+        this.view.notifyChange(layer);
+    }).bind(this));
+};
+
 GuiTools.prototype.addImageryLayersGUI = function addImageryLayersGUI(layers) {
     var i;
     var seq = itowns.ImageryLayers.getColorLayersIdOrderedBySequence(layers);

src/Converter/convertToTile.js

@@ -81,6 +81,7 @@ export default {
 
             if (parent) {
                 tile.setBBoxZ({ min: parent.obb.z.min, max: parent.obb.z.max });
+                tile.geoidHeight = parent.geoidHeight;
             }
 
             return tile;

src/Core/Geographic/GeoidGrid.js

@@ -0,0 +1,135 @@
+import * as THREE from 'three';
+import Coordinates from 'Core/Geographic/Coordinates';
+import CRS from 'Core/Geographic/Crs';
+
+
+const coord = new Coordinates('EPSG:4326');
+const indexes = new THREE.Vector2();
+
+
+function biLinearInterpolation(indexes, getData) {
+    const j = Math.floor(indexes.x);
+    const i = Math.floor(indexes.y);
+
+    const u = indexes.x - j;
+    const v = indexes.y - i;
+
+    return (1 - u) * (
+        (1 - v) * getData(i, j) + v * getData(i + 1, j)
+    ) + u * (
+        (1 - v) * getData(i, j + 1) + v * (getData(i + 1, j + 1))
+    );
+}
+
+
+/**
+ * An instance of `GeoidGrid` allows accessing some geoid height grid data from geographic instances (like some
+ * `{@link Coordinates}`). The geoid height grid data must contain geoid height values for a set of geographic points
+ * regularly dispatched on a planar surface.
+ *
+ * @property    {Extent}        extent      The geographic extent of the geoid height grid data.
+ * @property    {THREE.Vector2} step        The distance between two consecutive points of the geoid height grid. The
+                                            * `x` value stands for the distance along the West-East direction, and the
+                                            * `y` value stands for the distance along the South-North direction.
+ * @property    {THREE.Vector2} dimensions  The planar dimensions of the geoid height grid data extent.
+ * @property    {THREE.Vector2} dataSize    The number of values in the gridded data along the West-East direction (`x`
+                                            * axis) and the South-North direction (`y` axis).
+ *
+ * @example
+ * // Create a set of gridded data.
+ * const data = [
+ *     [1, 2, 3],
+ *     [2, 3, 4],
+ *     [3, 4, 5],
+ * ];
+ * // This set of data presents the following spatial distribution of geoid heights values :
+ * //
+ * //    Latitudes  ^
+ * //               |
+ * //         41.0  |   3  4  5
+ * //         40.5  |   2  3  4
+ * //         40.0  |   1  2  3
+ * //               |------------->
+ * //                   1  2  3     Longitudes
+ *
+ * // Create a GeoidGrid allowing to access the gridded data.
+ * const geoidGrid = new GeoidGrid(
+ *     new Extent('EPSG:4326', 1, 3, 40, 41),
+ *     new THREE.Vector2(1, 0.5),
+ *     (verticalIndex, horizontalIndex) => data[verticalIndex][horizontalIndex],
+ * );
+ *
+ * // Access a value of geoid height at some geographic coordinates.
+ * // The value is interpolated from the gridded data.
+ * const value = geoidGrid.getHeightAtCoordinates(
+ *     new Coordinates('EPSG:4326', 1.5, 40.25)
+ * );
+ * // This should return 2.0, which is the result from the bi-linear
+ * // interpolation at the center of the `[[1, 2], [2, 3]]` subsection
+ * // of the grid data.
+ */
+class GeoidGrid {
+    /**
+     * @param   {Extent}        extent      The geographic extent of the geoid height grid data.
+     * @param   {THREE.Vector2} step        The distance between two consecutive points of the geoid height grid. The
+                                            * `x` value stands for the distance along the West-East direction, and the
+                                            * `y` value stands for the distance along the South-North direction.
+     * @param   {function}      getData     A method that allows reading a value in the geoid height grid from its
+                                            * vertical and horizontal indexes. The lower an index, the lower the
+                                            * coordinate on the corresponding axis - 0 being the index of the minimal
+                                            * coordinate of the gridded data on a given axis. In other words :
+                                            * - `getData(0, 0)` must return the geoid height value at the SOUTH-WEST
+                                            *   corner of your data extent.
+                                            * - `getData(0, j)` must return a geoid height on the southern limit of your
+                                            *   data extent.
+                                            * - `getData(i, 0)` must return a geoid height on the western limit of your
+                                            *   data extent.
+                                            * - if your gridded data has dimensions (rowNumber, colNumber),
+                                            *   `getData(rowNumber - 1, colNumber - 1)` must return the geoid height at
+                                            *   the NORTH-EAST corner of your data extent.
+     */
+    constructor(
+        extent,
+        step,
+        getData,
+    ) {
+        CRS.isGeographic(extent.crs);
+
+        this.extent = extent;
+        this.step = new THREE.Vector2(step.x, step.y || step.x);
+
+        this.dimensions = this.extent.planarDimensions();
+        this.dataSize = new THREE.Vector2().addVectors(this.step, this.dimensions).divide(this.step).round();
+
+        this.getData = getData;
+    }
+
+    /**
+     * Get the value of the geoid height at given geographic `{@link Coordinates}`. The geoid height value is
+     * bi-linearly interpolated from the gridded data accessed by the `GeoidGrid` instance.
+     *
+     * @param   {Coordinates}   coordinates     Geographic coordinates to get the geoid height value at.
+     *
+     * @returns {number}    The geoid height value at the given `{@link Coordinates}`, bi-interpolated from the gridded
+                            * data accessed by the `GeoidGrid` instance.
+     */
+    getHeightAtCoordinates(coordinates) {
+        coordinates.as(this.extent.crs, coord);
+
+        indexes.set(
+            (this.dataSize.x - 1) * (coord.x - this.extent.west) / (this.dimensions.x),
+            (this.dataSize.y - 1) * (coord.y - this.extent.south) / (this.dimensions.y),
+        );
+
+        // TODO : add management for global GeoidGrid.
+        if (
+            indexes.x < 0 || indexes.x >= this.dataSize.x - 1
+            || indexes.y < 0 || indexes.y >= this.dataSize.y - 1
+        ) { return 0; }
+
+        return biLinearInterpolation(indexes, this.getData);
+    }
+}
+
+
+export default GeoidGrid;

src/Core/TileMesh.js

@@ -43,6 +43,8 @@ class TileMesh extends THREE.Mesh {
         this.isTileMesh = true;
 
         this.domElements = {};
+
+        this.geoidHeight = 0;
     }
 
     /**

src/Layer/GeoidLayer.js

@@ -0,0 +1,93 @@
+import Layer from 'Layer/Layer';
+import LayerUpdateState from 'Layer/LayerUpdateState';
+
+
+/**
+ * `GeoidLayer` is a specific `{@link Layer}` which supports geoid height data. When added to a `{@link View}`, it
+ * vertically translates each of the view's tiles by a proper geoid height value. For a given tile, the geoid height
+ * value used for translation is the geoid height computed at the center of the tile.
+ *
+ * @example
+ * // Create a GeoidLayer from a GTX geoid heights file.
+ * const geoidLayer = new GeoidLayer('geoid', {
+ *     source: new FileSource({
+ *         url: 'url-to-some-GTX-geoid-heights-file.gtx',
+ *         crs: 'EPSG:4326',
+ *         format: 'application/gtx',
+ *     }),
+ * });
+ */
+class GeoidLayer extends Layer {
+    /**
+     * Creates a new instance of `GeoidLayer`.
+     *
+     * @param   {string}    id              An unique identifier for the layer.
+     * @param   {Object}    config          The layer configuration. All elements in it will be merged as is in the
+                                            * layer. For example, if the configuration contains three elements `name,
+                                            * protocol, extent`, these elements will be available using `layer.name` or
+                                            * something else depending on the property name. Only `config.source`
+                                            * parameter is mandatory.
+     * @param   {Object}    config.source   The source of the geoid data displayed by the `GeoidLayer`. It is mandatory
+                                            * that the source data for a `GeoidLayer` be parsed into a
+                                            * `{@link GeoidGrid}`. You can refer to `{@link GTXParser}`,
+                                            * `{@link GDFParser}` and `{@link ISGParser}` to see how three standard
+                                            * geoid height grid file formats are parsed into `{@link GeoidGrid}`.
+     */
+    constructor(id, config = {}) {
+        super(id, config);
+        this.isGeoidLayer = true;
+        this.defineLayerProperty('visible', true);
+    }
+
+    updateNodeZ(node, parent) {
+        node.position.z += (this.visible ? 1 : -1) * (node.geoidHeight - parent.geoidHeight);
+        node.updateMatrix();
+        node.updateMatrixWorld(true);
+    }
+
+    update(context, layer, node, parent) {
+        if (!parent || !node.material) {
+            return;
+        }
+
+        // Don't update tile if its zoom is not within the layer's zoom limits
+        const extentsDestination = node.getExtentsByProjection(layer.crs);
+        const zoom = extentsDestination[0].zoom;
+        if (zoom > layer.zoom.max || zoom < layer.zoom.min) {
+            return;
+        }
+
+        if (node.layerUpdateState[layer.id] === undefined) {
+            node.layerUpdateState[layer.id] = new LayerUpdateState();
+
+            const updateNodeZ = () => this.updateNodeZ(node, parent);
+            layer.addEventListener('visible-property-changed', updateNodeZ);
+            node.addEventListener('dispose', () => {
+                layer.removeEventListener('visible-property-changed', updateNodeZ);
+            });
+        }
+
+        if (
+            layer.frozen
+            || !layer.visible
+            || !node.material.visible
+            || !node.layerUpdateState[layer.id].canTryUpdate()
+        ) {
+            return;
+        }
+
+        node.layerUpdateState[layer.id].newTry();
+
+        return this.getData(node.extent, extentsDestination).then(
+            (result) => {
+                node.geoidHeight = result.getHeightAtCoordinates(node.extent.center());
+                this.updateNodeZ(node, parent);
+
+                node.layerUpdateState[layer.id].noMoreUpdatePossible();
+            },
+        );
+    }
+}
+
+
+export default GeoidLayer;

src/Layer/TiledGeometryLayer.js

@@ -66,6 +66,8 @@ class TiledGeometryLayer extends GeometryLayer {
         super(id, object3d, config);
 
         this.isTiledGeometryLayer = true;
+        // TODO : this should be add in a preprocess method specific to GeoidLayer.
+        this.object3d.geoidHeight = 0;
 
         this.protocol = 'tile';
 

src/Main.js

@@ -6,6 +6,7 @@ export const REVISION = conf.version;
 // Geographic tools
 export { default as Extent } from 'Core/Geographic/Extent';
 export { default as Coordinates } from 'Core/Geographic/Coordinates';
+export { default as GeoidGrid } from 'Core/Geographic/GeoidGrid';
 export { default as CRS } from 'Core/Geographic/Crs';
 
 export { default as Ellipsoid, ellipsoidSizes } from 'Core/Math/Ellipsoid';
@@ -60,6 +61,7 @@ export { default as GlobeLayer } from 'Core/Prefab/Globe/GlobeLayer';
 export { default as PlanarLayer } from 'Core/Prefab/Planar/PlanarLayer';
 export { default as LabelLayer } from 'Layer/LabelLayer';
 export { default as EntwinePointTileLayer } from 'Layer/EntwinePointTileLayer';
+export { default as GeoidLayer } from 'Layer/GeoidLayer';
 
 // Sources provided by default in iTowns
 // A custom source should at least implements Source

src/Process/FeatureProcessing.js

@@ -90,6 +90,7 @@ export default {
                 group.layer = layer;
                 group.matrixWorld.copy(node.matrixWorld).invert();
                 group.matrixWorld.decompose(group.position, group.quaternion, group.scale);
+                group.position.z += node.geoidHeight;
                 node.add(group.add(result));
                 group.updateMatrixWorld(true);
             } else {

test/unit/geoidgrid.js

@@ -0,0 +1,32 @@
+import assert from 'assert';
+import Extent from 'Core/Geographic/Extent';
+import GeoidGrid from 'Core/Geographic/GeoidGrid';
+import Coordinates from 'Core/Geographic/Coordinates';
+
+
+describe('GeoidGrid', function () {
+    let geoidGrid;
+    const data = [
+        [1, 2, 3, 4, 5],
+        [2, 3, 4, 5, 6],
+        [3, 4, 5, 6, 7],
+        [4, 5, 6, 7, 8],
+        [5, 6, 7, 8, 9],
+    ];
+
+    before(function () {
+        geoidGrid = new GeoidGrid(
+            new Extent('EPSG:4326', 1, 5, 1, 5),
+            { x: 1 },
+            (verticalIndex, horizontalIndex) => data[verticalIndex][horizontalIndex],
+        );
+    });
+
+    it('should return the correct geoid height from coordinates', function () {
+        assert.strictEqual(geoidGrid.getHeightAtCoordinates(new Coordinates('EPSG:4326', 1.5, 1.5)), 2);
+    });
+
+    it('should return a null geoid height for coordinates outside of data extent', function () {
+        assert.strictEqual(geoidGrid.getHeightAtCoordinates(new Coordinates('EPSG:4326', 6, 6)), 0);
+    });
+});