Skip to content
Snippets Groups Projects
Commit 07b49291 authored by Luc Maisonobe's avatar Luc Maisonobe
Browse files

Updated tutorials from forge wiki.

parent 57247f33
No related branches found
No related tags found
No related merge requests found
<!--- Copyright 2013-2014 CS Systèmes d'Information
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
-->
# Direct Location with a DEM
The aim of this tutorial is to compute a direct location grid by intersection of the line of sight with a DEM (Digital Elevation Model), using Duvenhage's algorithm. This algorithm is the most performant one in Rugged.
## Feeding Rugged with DEM tiles
Rugged does not parse DEM files but takes buffers of elevation data as input. It is up to the calling application to read the DEM and load the data into buffers. Rugged provides a tile mecanism with cache for large DEMs allowing the user to load one tile at a time. This is in line with the format of world coverage DEMs such as SRTM. Rugged offers an interface for updating the DEM tiles in cache with a callback function triggered everytime a coordinate falls outside the current region.
The calling application must implement the callback function for loading the tiles. We recommend to use GDAL (http://www.gdal.org/) to parse the DEM files as it handles most formats (including geoTIFF for Aster DEM or DTED for SRTM). Rugged does not include the parsing of the DEM, by design, and yet we could have used GDAL. We want Rugged to remain a low-level library that does not pull too many third-party libraries.
### 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 volcanoe 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.apache.commons.math3.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.api.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. 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, in a similar way as for the SRTM DEMs.
* 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 proper DEM tiling with one line/column overlaps between neighbouring tiles :
![DEM tiles overlap](../images/DEM-tiles-overlap.png)
This diagram tries to represent the meaning of the different parameters in the definition of a tile :
![tile description](../images/tile-description.png)
## Initializing Rugged with a DEM
The initialization step differs slightly from the first tutorial [[DirectLocation|Direct location]], as we need to pass the information about our TileUpdater.
Instantiate an object derived from TileUpdater :
TileUpdater updater = new VolcanicConeElevationUpdater(summit, slope, base, FastMath.toRadians(1.0), 100);
int nbTiles = 8 ; //number max of tiles in Rugged cache
AlgorithmId algoId = AlgorithmId.DUVENHAGE;
Initialize Rugged with these parameters :
Rugged rugged = new Rugged(updater, nbTiles, algoId,
EllipsoidId.WGS84, InertialFrameId.EME2000, BodyRotatingFrameId.ITRF,
startDate, stopDate, 0.1, 10.0,
satellitePVList, 6, CartesianDerivativesFilter.USE_P,
satelliteQList, 8, AngularDerivativesFilter.USE_R);
## Computing a direct location grid
In a similar way as in the first tutorial [direct location](./direct-location.html), we call Rugged direct location method. This time it is called in a loop so as to generate a full grid on disk.
DataOutputStream dos = new DataOutputStream(new FileOutputStream("demDirectLoc.c1"));
int lineStep = (maxLine - minLine) / nbLineStep;
int pxStep = (maxPx - minPx) / nbPxStep;
ArrayList<GeodeticPoint> pointList = new ArrayList<GeodeticPoint>();
for (int i = 0; i < nbLineStep; i++) {
int currentLine = minLine + i * lineStep;
for (int j = 0; j < nbPxStep; j++) {
int currentPx = minPx + j * pxStep;
// Call to direct localization on current line
Vector3D position = lineSensor.getPosition();
AbsoluteDate currentLineDate = lineSensor.getDate(currentLine);
Vector3D los = lineSensor.getLos(absDate, currentPx);
GeodeticPoint point = rugged.directLocation(currentLineDate, position, los);
// Write latitude
if (point != null) {
dos.writeFloat((float) FastMath.toDegrees(point.getLatitude()));
} else {
dos.writeFloat((float) base);
}
// Store the GeodeticPoint to write longitude and altitude later
pointList.add(point);
}
}
for (GeodeticPoint point : pointList) {
if (point != null) {
dos.writeFloat((float) FastMath.toDegrees(point.getLongitude()));
} else {
dos.writeFloat((float) base);
}
}
for (GeodeticPoint point : pointList) {
if (point != null) {
dos.writeFloat((float) point.getAltitude());
} else {
dos.writeFloat((float) base);
}
}
## Source code
The source code is available in DirectLocationGrid.java
......@@ -21,8 +21,9 @@ list of positions, velocities and attitude quaternions recorded during the acqui
passing all this information to Rugged, we will be able to precisely locate each point of
the image on the Earth. Well, not exactly precise, as this first tutorial does not use a
Digital Elevation Model, but considers the Earth as an ellipsoid. The DEM will be added in
a second tutorial. The objective here is limited to explain how to initialize everything
Rugged needs to know about the sensor and the acquisition.
a second tutorial [Direct location with DEM](./direct-location-with-DEM.html). The objective
here is limited to explain how to initialize everything Rugged needs to know about the sensor
and the acquisition.
## Sensor's definition
......@@ -31,7 +32,6 @@ Rugged needs to know about the sensor and the acquisition.
Let's start by defining the imager. The sensor model is described by its viewing geometry
(i.e. the line of sight of each physical pixel) and by its datation model.
### Line of sight
......@@ -58,8 +58,8 @@ The viewing direction of pixel i with respect to the instrument is defined by th
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
Rotation rotation = new Rotation(axis, FastMath.toRadians(10));
viewDir = rotation.applyTo(viewDir);
Rotation rotation = new Rotation(Vector3D.PLUS_I, FastMath.toRadians(10));
Vector3D los = rotation.applyTo(viewDir);
The viewing direction is the line of sight of one pixel. In the case of a single line sensor,
Rugged expects to receive an array of LOS. The LOS can be time dependent (useful when applying
......@@ -70,7 +70,6 @@ instantiated from the viewing direction vector.
List<TimeDependentLOS> lineOfSight = new ArrayList<TimeDependentLOS>();
lineOfSight.add(new FixedLOS(los));
### Datation model
The datation model gives the line number in the image as a function of time, and vice-versa the
......@@ -115,13 +114,12 @@ auxiliary telemetry.
Note that for simulation purpose, we could also use Orekit to simulate the orbit. It offers very
convenient functions to propagate sun-synchronous orbits with yaw steering compensation (typical
orbits for Earth Observation satellites).
orbits for Earth Observation satellites).
### Reference frames
Rugged expects the positions, velocities and quaternions to be expressed in an inertial frame.
Rugged expects the positions (unit: m), velocities (unit: m/s) and quaternions to be expressed in an inertial frame.
All standard inertial and terrestrial frames are implemented in Orekit, so we just need to specify
which frames we are using and convert if necessary.
......@@ -131,49 +129,54 @@ bulletins and other physical data are provided within the orekit-data folder. Th
to configure Orekit to use this data. More information is given
[here](https://www.orekit.org/forge/projects/orekit/wiki/Configuration).
In our application, we simply need to know the name of the frames we are working with. Position and
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.
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.utils.IERSConventions;
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
The attitude quaternions are grouped in a list of TimeStampedAngularCoordinates objects:
import org.orekit.utils.TimeStampedAngularCoordinates;
List<TimeStampedAngularCoordinates> satelliteQList = new ArrayList<TimeStampedAngularCoordinates>();`
List<TimeStampedAngularCoordinates> satelliteQList = new ArrayList<TimeStampedAngularCoordinates>();
Each attitude sample (quaternion, time) is added to the list
Each attitude sample (quaternion, time) is added to the list,
AbsoluteDate attitudeDate = new AbsoluteDate(gpsDateAsString, TimeScalesFactory.getGPS());
Rotation rotation = new Rotation(q0, q1, q2, q3, true); //q0 is the scalar term
Rotation rotation = new Rotation(q0, q1, q2, q3, true); // q0 is the scalar term
Vector3D rotationRate = Vector3D.ZERO;
TimeStampedAngularCoordinates pair = new TimeStampedAngularCoordinates(attitudeDate, rotation, rotationRate);
satelliteQList.add(pair);`
Vector3D rotationAcceleration = Vector3D.ZERO;
TimeStampedAngularCoordinates pair = new TimeStampedAngularCoordinates(attitudeDate, rotation, rotationRate, rotationAcceleration);
satelliteQList.add(pair);
where, for instance, gpsDateAsString is set to "2009-12-11T10:49:55.899994"
### Position and velocities
Similarly the position and velocities will be set in a list of `TimeStampedPVCoordinates`. Before being
Similarly the positions and velocities will be set in a list of `TimeStampedPVCoordinates`. Before being
added to the list, they must be transformed to EME2000:
import org.orekit.utils.TimeStampedPVCoordinates;
import org.orekit.utils.PVCoordinates;
import org.orekit.frames.Transform;
List<TimeStampedPVCoordinates> satellitePVList = new ArrayList<TimeStampedPVCoordinates>();
AbsoluteDate ephemerisDate = new AbsoluteDate(gpsDateAsString, TimeScalesFactory.getGPS());
Vector3D position = new Vector3D(px, py, pz);
Vector3D velocity = new Vector3D(vx, vy, vz);
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(pITRF);
satellitePVList.add(new TimeStampedPVCoordinates(ephemerisDate, pEME2000, vEME2000));
PVCoordinates pvEME2000 = transform.transformPVCoordinates(pvITRF);
satellitePVList.add(new TimeStampedPVCoordinates(ephemerisDate, pvEME2000.getPosition(), pvEME2000.getVelocity(), Vector3D.ZERO)));
## Rugged initialization
......@@ -190,13 +193,13 @@ Finally we can initialize Rugged. It looks like this:
import org.orekit.utils.IERSConventions;
Rugged rugged = new Rugged (demTileUpdater, nbTiles, demAlgoId,
EllipsoidId.WGS84, InertialFrameId.EME2000, BodyRotatingFrameId.ITRF,
acquisitionStartDate, acquisitionStopDate, timeTolerance,
acquisitionStartDate, acquisitionStopDate, tStep, timeTolerance,
satellitePVList, nbPVPoints, CartesianDerivativesFilter.USE_P,
satelliteQList, nbQPoints, AngularDerivativesFilter.USE_R, 0.1);
satelliteQList, nbQPoints, AngularDerivativesFilter.USE_R);
Argh, that sounds complicated. It is not so difficult since we have already defined most of what is
needed. Let's describe the arguments line by line:
The first 3 are related to the DEM. Since we are ignoring the DEM, they can be set respectively to
`null`, `0` and `AlgorithmId.IGNORE_DEM_USE_ELLIPSOID`.
......@@ -204,7 +207,8 @@ The next 3 arguments defines the ellipsoid and the reference frames: `EllipsoidI
`InertialFrameId.EME2000`, `BodyRotatingFrameId.ITRF`
On the third line, we find the time interval of the acquisition: acquisitionStartDate, acquisitionStopDate,
timeTolerance (margin allowed for extrapolation, in seconds). This is an important information as Rugged
tStep (step at which the pre-computed frames transforms cache will be filled), timeTolerance (margin
allowed for extrapolation, in seconds). This is an important information as Rugged
will pre-compute a lot of stuff at initialization in order to optimize further calls to the direct and
inverse location functions.
......@@ -217,25 +221,18 @@ The sensor models are added after initialization. We can add as many as we want.
rugged.addLineSensor(lineSensor);
## Direct location
Finally everything is set to do some real work. Let's try to locate a point on Earth
Upper left point (first line, first pixel):
import org.orekit.bodies.GeodeticPoint;
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(1);
Vector3D los = lineSensor.getLos(firstLineDate, 0);
GeodeticPoint upLeftPoint = rugged.directLocation(firstLineDate, position, los);
The results are ....
## Source code
The source code is available in ....
The source code is available in DirectLocation.java
src/site/resources/images/DEM-tiles-overlap.png

9.2 KiB

src/site/resources/images/tile-description.png

14.9 KiB

......@@ -26,29 +26,30 @@
</bannerRight>
<body>
<menu name="Rugged">
<item name="Overview" href="/index.html" />
<item name="Getting the sources" href="/sources.html" />
<item name="Building" href="/building.html" />
<item name="Configuration" href="/configuration.html" />
<item name="FAQ" href="/faq.html" />
<item name="License" href="/license.html" />
<item name="Downloads" href="/downloads.html" />
<item name="Changes" href="/changes-report.html" />
<item name="Contact" href="/contact.html" />
<item name="Overview" href="/index.html" />
<item name="Getting the sources" href="/sources.html" />
<item name="Building" href="/building.html" />
<item name="Configuration" href="/configuration.html" />
<item name="FAQ" href="/faq.html" />
<item name="License" href="/license.html" />
<item name="Downloads" href="/downloads.html" />
<item name="Changes" href="/changes-report.html" />
<item name="Contact" href="/contact.html" />
</menu>
<menu name="Design">
<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="Preliminary design" href="/design/preliminary-design.html" />
<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="Preliminary design" href="/design/preliminary-design.html" />
</menu>
<menu name="Tutorials">
<item name="Direct location" href="/tutorials/direct-location.html" />
<item name="Direct location" href="/tutorials/direct-location.html" />
<item name="Direct location with DEM" href="/tutorials/direct-location-with-DEM.html" />
</menu>
<menu name="Development">
<item name="Contributing" href="/contributing.html" />
<item name="Guidelines" href="/guidelines.html" />
<item name="Javadoc" href="/apidocs/index.html" />
<item name="Contributing" href="/contributing.html" />
<item name="Guidelines" href="/guidelines.html" />
<item name="Javadoc" href="/apidocs/index.html" />
</menu>
<menu ref="reports"/>
</body>
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment