Direct Location with a DEM

Implementing the interface TileUpdater for DEM loading.

In this tutorial, we will not include real DEM data. Instead we are going to create a fake DEM representing a volcano in a form of a perfect cone, similar to the Mayon volcano in the Philippines, except that we will locate it somewhere just below our satellite. This example is already part of Rugged tests cases, the source code is available in the package org.orekit.rugged.raster, file VolcanicConeElevationUpdater.java.

The class VolcanicConeElevationUpdater implements the interface TileUpdater with its method updateTile. The method is in charge of loading a tile. The extent of the tile must be such that it covers at least the ground point with coordinates (latitude, longitude) which are passed as arguments to the method. The tile is an object of type UpdatableTile which has two methods :

• setGeometry(minLatitude, minLongitude, latitudeStep, longitudeStep, latitudeRows, longitudeRows) : initializes the extent of the new tile before loading

• setElevation(latIdx, longIdx, elevation) : fills the tile buffer at indices (latIdx, lonIdx) with value elevation

Here's the source code of the class VolcanicConeElevationUpdater :

import org.hipparchus.util.FastMath;
import org.orekit.rugged.raster.TileUpdater;
import org.orekit.rugged.raster.UpdatableTile;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.utils.Constants;

public class VolcanicConeElevationUpdater implements TileUpdater {

    private GeodeticPoint summit;
    private double        slope;
    private double        base;
    private double        size;
    private int           n;

    public VolcanicConeElevationUpdater(GeodeticPoint summit, double slope, double base,
                                        double size, int n) {
        this.summit = summit;
        this.slope  = slope;
        this.base   = base;
        this.size   = size;
        this.n      = n;
    }

    public void updateTile(double latitude, double longitude, UpdatableTile tile)
        throws RuggedException {
        double step         = size / (n - 1);
        double minLatitude  = size * FastMath.floor(latitude  / size);
        double minLongitude = size * FastMath.floor(longitude / size);
        double sinSlope     = FastMath.sin(slope);
        tile.setGeometry(minLatitude, minLongitude, step, step, n, n);
        for (int i = 0; i < n; ++i) {
            double cellLatitude = minLatitude + i * step;
            for (int j = 0; j < n; ++j) {
                double cellLongitude = minLongitude + j * step;
                double distance       = Constants.WGS84_EARTH_EQUATORIAL_RADIUS *
                                        FastMath.hypot(cellLatitude  - summit.getLatitude(),
                                                       cellLongitude - summit.getLongitude());
                double altitude = FastMath.max(summit.getAltitude() - distance * sinSlope,
                                               base);
                tile.setElevation(i, j, altitude);
            }
        }
    }

}

Important notes on DEM tiles :

• Ground point elevation are obtained by bilinear interpolation between 4 neighbouring cells. There is no specific algorithm for border management in this computation. As a consequence, a point falling on the border of the tile is considered outside. DEM tiles must be overlapping by at least one line/column in all directions. Two options are available:

  • until version 3.0, the DEM tiles should comply to the overlapping obligation.
  • after version 3.0, the DEM tiles could be seamless (i.e. not overlapping). In that case, a specific algorithm creates zipper tiles on the fly, in order to comply to the overlapping obligation. The case of overlapping DEM tiles is still possible. The choice is given through the API setDigitalElevationModel (see below paragraph Initializing Rugged with a DEM).

• In Rugged terminology, the minimum latitude and longitude correspond to the centre of the farthest Southwest cell of the DEM. Be careful if using GDAL to pass the correct information as there is half a pixel shift with respect to the lower left corner coordinates in gdalinfo.

The following diagram illustrates needed DEM tiling with one line/column overlaps between neighbouring tiles :

This diagram tries to represent the meaning of the different parameters in the definition of a tile :