From 21688139ebf92b21ead7a46e8b63d718778ec233 Mon Sep 17 00:00:00 2001
From: gprat <guylaine.prat@csgroup.eu>
Date: Wed, 7 Dec 2022 16:37:10 +0100
Subject: [PATCH] Add a tutorial about how to init a tile updater with real DEM
(SRTM)
Fixes #393
---
src/site/markdown/tutorials/tile-updater.md | 732 ++++++++++++++++++++
src/site/site.xml | 3 +-
2 files changed, 734 insertions(+), 1 deletion(-)
diff --git a/src/site/markdown/tutorials/tile-updater.md b/src/site/markdown/tutorials/tile-updater.md
index 1f450887..d19bf873 100644
--- a/src/site/markdown/tutorials/tile-updater.md
+++ b/src/site/markdown/tutorials/tile-updater.md
@@ -12,3 +12,735 @@
limitations under the License.
-->
+# Example with the real DEM SRTM
+
+This tutorial shows how to use the tile updater with a real DEM (for instance SRTM here).
+
+Same case as the tutorial [Direct location with a DEM](./direct-location-with-DEM.html).
+
+**WARNING**
+
+This tutorial will not be able to compile under [Rugged gitlab CI](https:gitlab.orekit.org/orekit/rugged/-/pipelines) as GDAL is not authorized in Rugged official releases.
+Don't commit under official branches (master, develop, release, ...) of Rugged, a pom.xml with the dependency to GDAL.
+
+## Needs
+
+### Install GDAL
+
+Search the web for how to install GDAL according to your OS.
+
+For instance for Linux/Ubuntu:
+
+The first time you want to install GDAL, you need to add a specific repository:
+
+For GDAL 3.x, use the following repository:
+
+ sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
+
+For GDAL 2.x, use the following repository:
+
+ sudo add-apt-repository ppa:ubuntugis/ppa
+
+Then
+
+ sudo apt update
+ sudo apt install gdal-bin gdal-data libgdal-java
+
+### Add GDAL in pom.xml file
+
+In `<properties>` part:
+
+ <!-- GDAL version -->
+ <rugged.gdal.version>3.0.0</rugged.gdal.version>
+ <!-- GDAL native library path -->
+ <rugged.gdal.native.library.path>${env.GDAL_PATH}</rugged.gdal.native.library.path>
+
+and in `<dependencies>` part:
+
+ <dependency>
+ <groupId>org.gdal</groupId>
+ <artifactId>gdal</artifactId>
+ <version>${rugged.gdal.version}</version>
+ <type>jar</type>
+ <optional>false</optional>
+ </dependency>
+
+### Configure the environment variable GDAL_PATH
+
+To access GDAL libraries (.so files), one must configure the environment variable GDAL_PATH as defined in pom.xml (see above).
+
+For instance (for Linux) according to the GDAL version and the linux distribution:
+
+ export GDAL_PATH=/COTS/gdal-3.0.0/swig/java/
+ or
+
+ export GDAL_PATH=/usr/lib/jni
+
+### Get the following SRTM tile
+
+This example needs one tile under a specific directory: `17/05/srtm_17_05.tif`
+
+For instance from : [SRTM Tile Grabber](http://dwtkns.com/srtm/)
+
+ srtm_17_05.zip
+
+and unzip it under a directory like : `mypath/dem-data/SRTM/17/05/`
+or changed the method: `getRasterFilePath(latitude, longitude)`
+
+### Get the GEOID data: egm96_15.gtx
+
+This example needs the geoid file: `egm96_15.gtx`
+
+For instance from the [GeographicLib](https:geographiclib.sourceforge.io/C++/doc/geoid.html)
+
+### Update the path in variables demRootDir and geoidFilePath
+
+Into the following code ...
+
+## Code example
+
+This is a whole example using the SRTM tiles.
+
+ public class DirectLocationWithSRTMdem {
+
+ // Root dir for the SRTM tiles
+ final static String demRootDir = "/mypath/dem-data/SRTM";
+
+ // File path for the GEOID
+ final static String geoidFilePath = "/mypath/dem-data/GEOID/egm96_15.gtx";
+
+ // Geoid elevations
+ static GEOIDelevationReader geoidElevationReader;
+ static GdalTransformTools geoidGTTools;
+ static double[][] geoidElevations;
+
+
+ public static void main(String[] args) {
+
+ try {
+
+ // Initialize Orekit, assuming an orekit-data folder is in user home directory
+ File home = new File(System.getProperty("user.home"));
+ File orekitData = new File(home, "orekit-data");
+ DataContext.getDefault().getDataProvidersManager().addProvider(new DirectoryCrawler(orekitData));
+
+
+ // Initialize the DEM datas
+ // ========================
+ // Initialize GDAL
+ org.gdal.gdal.gdal.AllRegister();
+
+ // Get the GEOID data once for all
+ geoidElevationReader = new GEOIDelevationReader(geoidFilePath);
+ geoidGTTools = new GdalTransformTools(geoidElevationReader.getGeoidGeoTransformInfo(),
+ geoidElevationReader.getGeoidLonSize(), geoidElevationReader.getGeoidLatSize());
+ geoidElevations = geoidElevationReader.getElevationsWholeEarth();
+
+ // Initialize the SRTM tile updater
+
+ // SRTM elevations are given over the GEOID.
+ // Rugged needs elevations over the ellipsoid so the Geoid must be added to raw SRTM elevations.
+ SRTMelevationUpdater srtmUpdater = new SRTMelevationUpdater(demRootDir);
+
+
+ // Sensor's definition
+ // ===================
+ // Line of sight
+ // -------------
+ // The raw viewing direction of pixel i with respect to the instrument is defined by the vector:
+ List<Vector3D> rawDirs = new ArrayList<Vector3D>();
+ for (int i = 0; i < 2000; i++) {
+ // 20° field of view, 2000 pixels
+ rawDirs.add(new Vector3D(0d, i*FastMath.toRadians(20)/2000d, 1d));
+ }
+
+ // The instrument is oriented 10° off nadir around the X-axis, we need to rotate the viewing
+ // direction to obtain the line of sight in the satellite frame
+ LOSBuilder losBuilder = new LOSBuilder(rawDirs);
+ losBuilder.addTransform(new FixedRotation("10-degrees-rotation", Vector3D.PLUS_I, FastMath.toRadians(10)));
+
+ TimeDependentLOS lineOfSight = losBuilder.build();
+
+ // Datation model
+ // --------------
+ // We use Orekit for handling time and dates, and Rugged for defining the datation model:
+ TimeScale gps = TimeScalesFactory.getGPS();
+ AbsoluteDate absDate = new AbsoluteDate("2009-12-11T16:59:30.0", gps);
+ LinearLineDatation lineDatation = new LinearLineDatation(absDate, 1d, 20);
+
+ // Line sensor
+ // -----------
+ // With the LOS and the datation now defined, we can initialize a line sensor object in Rugged:
+ LineSensor lineSensor = new LineSensor("mySensor", lineDatation, Vector3D.ZERO, lineOfSight);
+
+
+ // Satellite position, velocity and attitude
+ // =========================================
+
+ // Reference frames
+ // ----------------
+ // In our application, we simply need to know the name of the frames we are working with. Positions and
+ // velocities are given in the ITRF terrestrial frame, while the quaternions are given in EME2000
+ // inertial frame.
+ Frame eme2000 = FramesFactory.getEME2000();
+ boolean simpleEOP = true; // we don't want to compute tiny tidal effects at millimeter level
+ Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, simpleEOP);
+
+ // Satellite attitude
+ // ------------------
+ ArrayList<TimeStampedAngularCoordinates> satelliteQList = new ArrayList<TimeStampedAngularCoordinates>();
+
+ addSatelliteQ(gps, satelliteQList, "2009-12-11T16:58:42.592937", -0.340236d, 0.333952d, -0.844012d, -0.245684d);
+ addSatelliteQ(gps, satelliteQList, "2009-12-11T16:59:06.592937", -0.354773d, 0.329336d, -0.837871d, -0.252281d);
+ addSatelliteQ(gps, satelliteQList, "2009-12-11T16:59:30.592937", -0.369237d, 0.324612d, -0.831445d, -0.258824d);
+ addSatelliteQ(gps, satelliteQList, "2009-12-11T16:59:54.592937", -0.3836d, 0.319792d, -0.824743d, -0.265299d);
+ addSatelliteQ(gps, satelliteQList, "2009-12-11T17:00:18.592937", -0.397834d, 0.314883d, -0.817777d, -0.271695d);
+ addSatelliteQ(gps, satelliteQList, "2009-12-11T17:00:42.592937", -0.411912d, 0.309895d, -0.810561d, -0.278001d);
+ addSatelliteQ(gps, satelliteQList, "2009-12-11T17:01:06.592937", -0.42581d, 0.304838d, -0.803111d, -0.284206d);
+ addSatelliteQ(gps, satelliteQList, "2009-12-11T17:01:30.592937", -0.439505d, 0.299722d, -0.795442d, -0.290301d);
+ addSatelliteQ(gps, satelliteQList, "2009-12-11T17:01:54.592937", -0.452976d, 0.294556d, -0.787571d, -0.296279d);
+ addSatelliteQ(gps, satelliteQList, "2009-12-11T17:02:18.592937", -0.466207d, 0.28935d, -0.779516d, -0.302131d);
+
+ // Positions and velocities
+ // ------------------------
+ ArrayList<TimeStampedPVCoordinates> satellitePVList = new ArrayList<TimeStampedPVCoordinates>();
+ addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:58:42.592937", -726361.466d, -5411878.485d, 4637549.599d, -2463.635d, -4447.634d, -5576.736d);
+ addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:04.192937", -779538.267d, -5506500.533d, 4515934.894d, -2459.848d, -4312.676d, -5683.906d);
+ addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:25.792937", -832615.368d, -5598184.195d, 4392036.13d, -2454.395d, -4175.564d, -5788.201d);
+ addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:47.392937", -885556.748d, -5686883.696d, 4265915.971d, -2447.273d, -4036.368d, -5889.568d);
+ addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:00:08.992937", -938326.32d, -5772554.875d, 4137638.207d, -2438.478d, -3895.166d, -5987.957d);
+ addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:00:30.592937", -990887.942d, -5855155.21d, 4007267.717d, -2428.011d, -3752.034d, -6083.317d);
+ addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:00:52.192937", -1043205.448d, -5934643.836d, 3874870.441d, -2415.868d, -3607.05d, -6175.6d);
+ addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:01:13.792937", -1095242.669d, -6010981.571d, 3740513.34d, -2402.051d, -3460.291d, -6264.76d);
+ addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:01:35.392937", -1146963.457d, -6084130.93d, 3604264.372d, -2386.561d, -3311.835d, -6350.751d);
+ addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:01:56.992937", -1198331.706d, -6154056.146d, 3466192.446d, -2369.401d, -3161.764d, -6433.531d);
+ addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:02:18.592937", -1249311.381d, -6220723.191d, 3326367.397d, -2350.574d, -3010.159d, -6513.056d);
+
+
+ // Rugged initialization
+ // ---------------------
+ Rugged rugged = new RuggedBuilder().
+ setAlgorithm(AlgorithmId.DUVENHAGE).
+ setDigitalElevationModel(srtmUpdater, SRTMelevationUpdater.NUMBER_TILES_CACHED, false).
+ setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF).
+ setTimeSpan(absDate, absDate.shiftedBy(60.0), 0.01, 5 / lineSensor.getRate(0)).
+ setTrajectory(InertialFrameId.EME2000,
+ satellitePVList, 4, CartesianDerivativesFilter.USE_P,
+ satelliteQList, 4, AngularDerivativesFilter.USE_R).
+ addLineSensor(lineSensor).
+ build();
+
+ // Direct location of a single sensor pixel (first line, first pixel)
+ // ------------------------------------------------------------------
+ Vector3D position = lineSensor.getPosition(); // This returns a zero vector since we set the relative position of the sensor w.r.T the satellite to 0.
+ AbsoluteDate firstLineDate = lineSensor.getDate(0);
+ Vector3D los = lineSensor.getLOS(firstLineDate, 0);
+ GeodeticPoint upLeftPoint = rugged.directLocation(firstLineDate, position, los);
+ System.out.format(Locale.US, "upper left point: φ = %8.3f °, λ = %8.3f °, h = %8.3f m (elevation with SRTM + Geoid))%n",
+ FastMath.toDegrees(upLeftPoint.getLatitude()),
+ FastMath.toDegrees(upLeftPoint.getLongitude()),
+ upLeftPoint.getAltitude());
+
+ // DEM SRTM with Geoid: (value of GEOID checked with QGIS around -30 meters)
+ // upper left point: φ = 37.585 °, λ = -96.950 °, h = 344.024 m
+
+ // DEM SRTM without Geoid (not to be used as elevation must be over ellipsoid for Rugged):
+ // upper left point: φ = 37.585 °, λ = -96.950 °, h = 373.675 m
+
+
+
+ } catch (OrekitException oe) {
+ System.err.println(oe.getLocalizedMessage());
+ System.exit(1);
+ } catch (RuggedException re) {
+ System.err.println(re.getLocalizedMessage());
+ System.exit(1);
+ }
+
+ }
+
+
+ private static void addSatellitePV(TimeScale gps, Frame eme2000, Frame itrf,
+ ArrayList<TimeStampedPVCoordinates> satellitePVList,
+ String absDate,
+ double px, double py, double pz, double vx, double vy, double vz) {
+
+ AbsoluteDate ephemerisDate = new AbsoluteDate(absDate, gps);
+ Vector3D position = new Vector3D(px, py, pz); // in ITRF, unit: m
+ Vector3D velocity = new Vector3D(vx, vy, vz); // in ITRF, unit: m/s
+ PVCoordinates pvITRF = new PVCoordinates(position, velocity);
+ Transform transform = itrf.getTransformTo(eme2000, ephemerisDate);
+ PVCoordinates pvEME2000 = transform.transformPVCoordinates(pvITRF);
+ satellitePVList.add(new TimeStampedPVCoordinates(ephemerisDate, pvEME2000.getPosition(), pvEME2000.getVelocity(), Vector3D.ZERO));
+
+ }
+
+ private static void addSatelliteQ(TimeScale gps, ArrayList<TimeStampedAngularCoordinates> satelliteQList, String absDate,
+ double q0, double q1, double q2, double q3) {
+ AbsoluteDate attitudeDate = new AbsoluteDate(absDate, gps);
+ Rotation rotation = new Rotation(q0, q1, q2, q3, true); // q0 is the scalar term
+ TimeStampedAngularCoordinates pair =
+ new TimeStampedAngularCoordinates(attitudeDate, rotation, Vector3D.ZERO, Vector3D.ZERO);
+ satelliteQList.add(pair);
+ }
+
+
+ /**
+ * To read SRTM tiles (get from http://dwtkns.com/srtm/)
+ */
+ private static class SRTMelevationUpdater implements TileUpdater {
+
+ /** Nb cached tiles send to init rugged. */
+ public static final int NUMBER_TILES_CACHED = 8;
+
+ /** Last lat value used to ask to update a tile. */
+ private double lastLat = Double.NaN;
+
+ /** Last lon value used to ask to update a tile. */
+ private double lastLon = Double.NaN;
+
+ /** Nb times the same lat/lon value is used to ask to update a tile (used to avoid infinite loop). */
+ private int nbCall = 0;
+
+ /** Nb time rugged ask exactly same DEM tile: means that infinite loop. */
+ private static final int NB_TIME_ASK_SAME_TILE = 1000;
+
+ /** Raster root directory, should contains raster files. */
+ private String srtmRootDirectory;
+
+
+ /**
+ * Constructor.
+ */
+ public SRTMelevationUpdater(final String SRTMrootDirectory) {
+
+ this.srtmRootDirectory = SRTMrootDirectory;
+ checkRasterDirectory(SRTMrootDirectory);
+ }
+
+ /** Update given tile using SRTM elevation.
+ * @param latitude latitude that must be covered by the tile (rad)
+ * @param longitude longitude that must be covered by the tile (rad)
+ * @param tile to update
+ */
+ @Override
+ public void updateTile(final double latitude, final double longitude, final UpdatableTile tile) {
+
+ // Check if latitude and longitude already known
+ if (latitude == this.lastLat && longitude == this.lastLon) {
+ this.nbCall++;
+ } else {
+ this.lastLat = latitude;
+ this.lastLon = longitude;
+ this.nbCall = 0;
+ }
+ if (this.nbCall > NB_TIME_ASK_SAME_TILE) {
+ String str = String.format("infinite loop for %3.8f long %3.8f lat ", longitude, latitude);
+ DummyLocalizable message = new DummyLocalizable(str);
+ throw new RuggedException(message);
+ }
+
+ // Compute the SRTM tile name according to the latitude and longitude
+ String rasterFileName = getRasterFilePath(latitude, longitude);
+
+ File rasterFile = new File(rasterFileName);
+ if (!rasterFile.exists()) {
+ // No SRTM tile found. We may be on water.
+ // For instance one can read the geoid ...
+ System.out.format("WARNING: No DEM tile found for latitude (deg) = %3.8f and longitude (deg) = %3.8f." +
+ " SRTM file %s wasn't found. %n",
+ FastMath.toDegrees(latitude), FastMath.toDegrees(longitude), rasterFileName);
+ return;
+ }
+
+ // Open the SRTM tile
+ org.gdal.gdal.Dataset srtmDataset = org.gdal.gdal.gdal.Open(rasterFileName, org.gdal.gdalconst.gdalconst.GA_ReadOnly);
+
+ // Get information about the tile
+ int rasterLonSize = srtmDataset.GetRasterXSize();
+ int rasterLatSize = srtmDataset.GetRasterYSize();
+ double[] srtmGeoTransformInfo = srtmDataset.GetGeoTransform();
+
+ double minLat = Double.NaN;
+ double minLon = Double.NaN;
+ double lonStep = Double.NaN;
+ double latStep = Double.NaN;
+
+ if (srtmGeoTransformInfo.length < 6) { // Check that the geoTransformInfo is correct
+ String str = "GDALGeoTransform has < 6 elements";
+ DummyLocalizable message = new DummyLocalizable(str);
+ throw new RuggedException(message);
+ } else {
+ // Get the latitudes and longitudes datas (in degrees)
+ lonStep = FastMath.abs(srtmGeoTransformInfo[1]);
+ latStep = FastMath.abs(srtmGeoTransformInfo[5]);
+
+ minLon = srtmGeoTransformInfo[0] + 0.5 * lonStep;
+ minLat = srtmGeoTransformInfo[3] - (rasterLatSize - 0.5) * latStep;
+ }
+
+ // Get the raster band from the tile
+ org.gdal.gdal.Band band = srtmDataset.GetRasterBand(1);
+
+ // Define Tile Geometry
+ // TBN: angles needed in radians
+ tile.setGeometry(FastMath.toRadians(minLat), FastMath.toRadians(minLon), FastMath.toRadians(latStep), FastMath.toRadians(lonStep), rasterLatSize, rasterLonSize);
+
+ // Loop over raster values
+ double[] data = new double[rasterLatSize * rasterLonSize];
+
+ // To use geo transform information from GDAL
+ GdalTransformTools srtmGTTools = new GdalTransformTools(srtmGeoTransformInfo, rasterLonSize, rasterLatSize);
+
+ // We read all raster at once to limit the numbers of Band.ReadRaster calls
+ band.ReadRaster(0, 0, rasterLonSize, rasterLatSize, data);
+
+ // Get the no data value from the raster
+ Double[] noDataValue = new Double[1];
+ band.GetNoDataValue(noDataValue);
+
+ // test if the no data value exists
+ Boolean noDataValueExist = false;
+ if (noDataValue[0] != null) noDataValueExist = true;
+
+ // from bottom left to upper right corner
+ for (int iLat = rasterLatSize - 1; iLat >= 0; iLat--) {
+ for (int jLon = 0; jLon < rasterLonSize; jLon++) {
+
+ double elevationOverGeoid = 0.0;
+
+ // Get elevation over GEOID (for SRTM)
+ elevationOverGeoid = data[iLat * rasterLonSize + jLon];
+
+ if (noDataValueExist && (elevationOverGeoid == noDataValue[0])) {
+ elevationOverGeoid = 0.0;
+ }
+
+ // The elevation value we send to rugged must be computed against ellipsoid
+ // => for SRTM , we must add geoid value
+ double lonDeg = srtmGTTools.getXFromPixelLine(jLon, iLat);
+ double latDeg = srtmGTTools.getYFromPixelLine(jLon, iLat);
+
+ // Get the geoid elevation at latDeg and lonDeg
+ int geoidIndexLatitude = (int) FastMath.floor(geoidGTTools.getLineNumFromLatLon(latDeg, lonDeg));
+ int geoidIndexLongitude = (int) FastMath.floor(geoidGTTools.getPixelNumFromLatLon(latDeg, lonDeg));
+
+ double elevationOverEllipsoid = elevationOverGeoid + geoidElevations[geoidIndexLatitude][geoidIndexLongitude];
+
+ // Set elevation over the ellipsoid
+ tile.setElevation(rasterLatSize - 1 - iLat, jLon, elevationOverEllipsoid);
+ }
+ }
+ band.delete();
+ band = null;
+ srtmDataset.delete();
+ srtmDataset = null;
+ }
+
+ /** Get the file path to the SRTM tile.
+ * Need the SRTM tif file named, for instance: srtm_17_05.tif under the directory rootDirectory/17/05/
+ * @param latitude latitude (rad)
+ * @param longitude longitude (rad)
+ * @return the SRTM tile path
+ */
+ private String getRasterFilePath(final double latitude, final double longitude) {
+
+ double latDeg = FastMath.toDegrees(latitude);
+ // Assure that the longitude belongs to [-180, + 180]
+ double lonDeg = FastMath.toDegrees(MathUtils.normalizeAngle(longitude, 0.0));
+
+ // Compute parent dir with longitude value according to SRTM naming convention
+ int parentDirValue = (int) (1 + (lonDeg + 180.) / 5);
+ String parentDir = String.format("%02d", parentDirValue);
+
+ // Compute sub dir with latitude value according to SRTM naming convention
+ int subDirValue = (int) (1 + (60. - latDeg) / 5);
+ String subDir = String.format("%02d", subDirValue);
+
+ String filePath = this.srtmRootDirectory + File.separator + parentDir + File.separator + subDir + File.separator +
+ "srtm_" + parentDir + "_" + subDir + ".tif";
+
+ return filePath;
+ }
+
+ /** Chech the root directory for the SRTM files and the presence of .tif
+ * @param directory
+ * @return true if everything OK
+ */
+ private boolean checkRasterDirectory(final String directory) {
+
+ if (directory == null) {
+
+ String str = "Directory not defined";
+ DummyLocalizable message = new DummyLocalizable(str);
+ throw new RuggedException(message);
+
+ } else {
+ try {
+ Path dir = FileSystems.getDefault().getPath(directory);
+ DirectoryStream<Path> stream = Files.newDirectoryStream(dir);
+ boolean found = false;
+ for (Path path : stream) {
+ if (!found) {
+ File currentFile = path.toFile();
+ if (currentFile.isDirectory()) {
+ found = checkRasterDirectory(currentFile.getAbsolutePath());
+ } else {
+ String filePath = currentFile.getAbsolutePath();
+ if (filePath.matches(".*.tif")) {
+ found = true;
+ }
+ }
+ if (found) {
+ stream.close();
+ return true;
+ }
+ }
+ }
+ stream.close();
+
+ String str = "SRTM raster file (.tif) not found in " + directory;
+ DummyLocalizable message = new DummyLocalizable(str);
+ throw new RuggedException(message);
+
+ } catch (IOException e) {
+ String str = "Directory for SRTM tiles not found " + directory;
+ DummyLocalizable message = new DummyLocalizable(str);
+ throw new RuggedException(message);
+ }
+ }
+ }
+ }
+
+ /**
+ * To read the GEOID
+ */
+ private static class GEOIDelevationReader {
+
+ /** Geoid dataset */
+ private org.gdal.gdal.Dataset geoidDataset = null;
+
+ /** Geoid elevations for the whole Earth (read from Geoid dataset).
+ * Filled in with GDAL convention : [0][0] = North West */
+ private double[][] elevationsWholeEarth = null;
+
+ /** the Geo transform info */
+ private double[] geoidGeoTransformInfo;
+
+ /** Size in longitude of the Geoid (number of columns/pixels). */
+ private int geoidLonSize;
+
+ /** Size in latitude of the Geoid (number of lines). */
+ private int geoidLatSize;
+
+ /**
+ * Constructor
+ * @param geoidFilePath path to geoid file
+ */
+ public GEOIDelevationReader(String geoidFilePath) {
+
+ // Open the Geoid dataset
+ openGeoid(geoidFilePath);
+
+ // Analyze and read the raster: initialize this.elevationsWholeEarth
+ this.geoidGeoTransformInfo = readRaster();
+ }
+
+ /**
+ * Analyse and read the Geoid raster. Initialize elevationsWholeEarth (fill in with GDAL convention : [0][0] = North West)
+ * @return the geoTransformInfo (read from the geoidDataset)
+ */
+ private double[] readRaster() {
+
+ // Analyse the Geoid dataset
+ this.geoidLonSize = geoidDataset.GetRasterXSize();
+ this.geoidLatSize = geoidDataset.GetRasterYSize();
+
+ // Get the Geo transform info: lon origin (deg), lon step (deg), 0, lat origin (deg), 0, lat step (deg)
+ double[] geoTransformInfo = geoidDataset.GetGeoTransform();
+ if (geoTransformInfo.length < 6) { // should not occur; if a problem occurs : the default transform value is set to (0,1,0,0,0,1)
+ String str = "GDALGeoTransform does not contains 6 elements";
+ DummyLocalizable message = new DummyLocalizable(str);
+ throw new RuggedException(message);
+ }
+
+ // Read the raster
+ // ===============
+ // To store the raster data
+ double[] data = new double[this.geoidLatSize * this.geoidLonSize];
+ // Get the raster
+ Band band = geoidDataset.GetRasterBand(1);
+ // We read all the raster once for all
+ band.ReadRaster(0, 0, this.geoidLonSize, this.geoidLatSize, data);
+ // Get the no data value from the raster
+ Double[] noDataValue = new Double[1];
+ band.GetNoDataValue(noDataValue);
+ // test if the no data value exists
+ Boolean noDataValueExist = false;
+ if (noDataValue[0] != null) noDataValueExist = true;
+
+ // Read the full raster all in once ...
+ this.elevationsWholeEarth = new double[this.geoidLatSize][this.geoidLonSize];
+
+ // from bottom left to upper right corner (GDAL convention)
+ for (int iLat = this.geoidLatSize - 1; iLat >= 0; iLat--) {
+ for (int jLon = 0; jLon < this.geoidLonSize; jLon++) {
+ // if jlon = 0 and ilat = rasterLatSize (721) : crash
+ // if iLat = 720 and jlon = rasterLonSize (1440) : crash
+
+ // Read the data with GDAL convention :
+ // iLat = 0 at North; iLat = rasterLatSize - 1 at South;
+ // jLon = 0 at West; jLon = rasterLonSize - 1 at East
+ double elevationOverEllipsoid = data[iLat * this.geoidLonSize + jLon];
+ if (noDataValueExist && (elevationOverEllipsoid == noDataValue[0])) {
+ elevationOverEllipsoid = 0.0;
+ }
+ // Set elevation for current latitude and longitude
+ // IMPORTANT : keep the GDAL convention [0][0] = North West
+ this.elevationsWholeEarth[iLat][jLon] = elevationOverEllipsoid;
+
+ } // end loop jLon
+ } // end loop iLat
+
+
+ // Close Dataset and Band once the raster is read ...
+ band.delete();
+ band = null;
+ geoidDataset.delete();
+ geoidDataset = null;
+
+ // Return the Geo Transform info
+ return geoTransformInfo;
+ }
+
+
+ public double[][] getElevationsWholeEarth() {
+ return elevationsWholeEarth;
+ }
+
+ public double[] getGeoidGeoTransformInfo() {
+ return geoidGeoTransformInfo;
+ }
+
+ public int getGeoidLonSize() {
+ return geoidLonSize;
+ }
+
+ public int getGeoidLatSize() {
+ return geoidLatSize;
+ }
+
+ /** Open the geoid dataset
+ * @param geoidFilePath the geoid file path
+ * @return the geoid dataset
+ */
+ public org.gdal.gdal.Dataset openGeoid(String geoidFilePath) throws RuggedException {
+
+ this.geoidDataset = (org.gdal.gdal.Dataset) org.gdal.gdal.gdal.Open(geoidFilePath, org.gdal.gdalconst.gdalconst.GA_ReadOnly);
+
+ if (this.geoidDataset == null) {
+ throw new RuggedException(RuggedMessages.UNKNOWN_TILE);
+
+ }
+ return this.geoidDataset;
+ }
+ }
+
+ /**
+ * Tool using GeoTransform info returned from GDAL.
+ */
+ static private class GdalTransformTools {
+
+ /** GeoTransform info from Gdal. */
+ private double[] gtInfo;
+
+ /** x size (longitude in deg). */
+ private double xSize;
+
+ /** y size (latitude in deg). */
+ private double ySize;
+
+ /**
+ * Constructor.
+ * @param gtInfo geo transform info from gdal
+ * @param xSize x size (longitude in deg)
+ * @param ySize y size (latitude in deg)
+ */
+ public GdalTransformTools(final double[] gtInfo, final double xSize, final double ySize) {
+ // We suppose gtInfo[5] and gtInfo[1] different from 0
+ // no tests are performed about gtInfo values ... as org.gdal.gdal.Dataset.GetGeoTransform will give
+ // a default transform value set to (0,1,0,0,0,1) if problems
+ this.gtInfo = gtInfo;
+ this.xSize = xSize;
+ this.ySize = ySize;
+ }
+
+ /**
+ * Get X (longitude) coord from pixel and line.
+ * @param pixelNum pixel number
+ * @param lineNum line number
+ * @return X coord (longitude in degrees) from pixel and line. For Gdal = Upper left corner of the cell.
+ */
+ public double getXFromPixelLine(final double pixelNum, final double lineNum) {
+ return gtInfo[0] + gtInfo[1] * pixelNum + gtInfo[2] * lineNum;
+ }
+
+ /**
+ * Get Y (latitude) coord from pixel and line.
+ * @param pixelNum pixel number
+ * @param lineNum line number
+ * @return Y coord (latitude in degrees) from pixel and line. For Gdal = Upper left corner of the cell.
+ */
+ public double getYFromPixelLine(final double pixelNum, final double lineNum) {
+ return gtInfo[3] + gtInfo[4] * pixelNum + gtInfo[5] * lineNum;
+ }
+
+ /**
+ * Get pixel number from latitude and longitude.
+ * @param lat latitude (deg)
+ * @param lon longitude (deg)
+ * @return pixel number (pixel number = index in longitude). For Gdal = Upper left corner of the cell.
+ */
+ public double getPixelNumFromLatLon(final double lat, final double lon) {
+ // We suppose gtInfo[5] and gtInfo[1] different from 0
+ double tmp = gtInfo[2] / gtInfo[5];
+ double value = (lon - tmp * lat - gtInfo[0] + tmp * gtInfo[3]) / (gtInfo[1] - tmp * gtInfo[4]);
+
+ // if value = xSize, don't convert value to get the inverse value from getLatMax and getLonMax
+ if (value >= xSize) {
+ value = value % xSize;
+ }
+ return value;
+ }
+
+ /**
+ * Get line number from latitude and longitude.
+ * @param lat latitude (deg)
+ * @param lon longitude (deg)
+ * @return line number (line number = index in latitude). For Gdal = Upper left corner of the cell.
+ */
+ public double getLineNumFromLatLon(final double lat, final double lon) {
+ // We suppose gtInfo[1] and gtInfo[5] different from 0
+ double tmp = gtInfo[4] / gtInfo[1];
+ double value = (lat - tmp * lon - gtInfo[3] + tmp * gtInfo[0]) / (gtInfo[5] - tmp * gtInfo[2]);
+ if (value >= ySize) { // if value = xSize, don't convert value to get the inverse value from getLatMax and getLonMax
+ value = value % ySize;
+ }
+ return value;
+ }
+ }
+ }
+
+
+**To be noticed:**
+
+if you are using a version of Rugged ≤ 3.0, you need to replace the call:
+
+ setDigitalElevationModel(srtmUpdater, SRTMelevationUpdater.NUMBER_TILES_CACHED, false)
+
+by
+
+ setDigitalElevationModel(srtmUpdater, SRTMelevationUpdater.NUMBER_TILES_CACHED)
+
+and use a DEM with overlapping tiles as explained in [Direct location with a DEM](./direct-location-with-DEM.html).
+
\ No newline at end of file
diff --git a/src/site/site.xml b/src/site/site.xml
index c29146d5..f17381ca 100644
--- a/src/site/site.xml
+++ b/src/site/site.xml
@@ -40,12 +40,13 @@
<item name="Overview" href="/design/overview.html" />
<item name="Technical choices" href="/design/technical-choices.html" />
<item name="Digital Elevation Model" href="/design/digital-elevation-model.html" />
- <item name="Design of the major functions" href="/design/design-major-functions.html" />
+ <item name="Design of the major functions" href="/design/design-major-functions.html" />
</menu>
<menu name="Tutorials">
<item name="Direct location" href="/tutorials/direct-location.html" />
<item name="Direct location with DEM" href="/tutorials/direct-location-with-DEM.html" />
<item name="Inverse location" href="/tutorials/inverse-location.html" />
+ <item name="Tile updater example" href="/tutorials/tile-updater.html" />
<item name="Matlab examples" href="/tutorials/matlab-example.html" />
</menu>
<menu name="Development">
--
GitLab