Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • orekit/rugged
  • sdinot/rugged
  • yzokras/rugged
  • youngcle/rugged-mod
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 3592 additions and 0 deletions
src/test/resources/replay/replay-inverse-loc-02.txt 0 → 100644
View file @ acf902ea
# Rugged library dump file, created on 2015-03-24T07:07:27Z
# all units are SI units (m, m/s, rad ...)
algorithm: DUVENHAGE
DEM tile: t0 latMin -4.014257279586958e-01 latStep 6.817692390602850e-05 latRows 257 lonMin 2.495820830351891e+00 lonStep 6.817692390602850e-05 lonCols 257
DEM cell: t0 latIndex 192 lonIndex 0 elevation -6.124011631478857e+02
DEM cell: t0 latIndex 128 lonIndex 128 elevation 3.151327114467009e+03
DEM cell: t0 latIndex 123 lonIndex 101 elevation 1.838123293450827e+03
DEM cell: t0 latIndex 123 lonIndex 102 elevation 1.870339491485550e+03
DEM cell: t0 latIndex 124 lonIndex 101 elevation 1.832853313250435e+03
DEM cell: t0 latIndex 124 lonIndex 102 elevation 1.866819199623890e+03
sensor: sensorName s0 nbPixels 2000 position 1.500000000000000e+00 0.000000000000000e+00 -2.000000000000000e-01
inverse location: sensorName s0 latitude -3.930305027145593e-01 longitude 2.502739578717502e+00 elevation 1.853045504720662e+03 minLine 0 maxLine 2000 lightTime true aberration true refraction false
sensor mean plane: sensorName s0 minLine 0 maxLine 2000 maxEval 50 accuracy 1.000000000000000e-02 normal 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 cachedResults 5 lineNumber 1.753079999999941e+03 date 2012-01-01T12:30:01.12961999999991Z target -4.733584770718903e+06 3.515842173828561e+06 -2.428300713188566e+06 targetDirection 4.440892098500626e-16 -7.496646388847980e-01 6.618178973144536e-01 -1.357198500050715e-08 -1.198160098465607e-08 -1.357198500050715e-08 lineNumber 1.753079999999936e+03 date 2012-01-01T12:30:01.12961999999990Z target -4.733634454810624e+06 3.515791377291927e+06 -2.428292161157289e+06 targetDirection 1.221245327087672e-15 -7.496479456645485e-01 6.618368058373009e-01 -1.357209409884309e-08 -1.198230643978911e-08 -1.357209409884309e-08 lineNumber 1.753080000000098e+03 date 2012-01-01T12:30:01.12962000000015Z target -4.733684171891571e+06 3.515740571824192e+06 -2.428283613116900e+06 targetDirection 2.664535259100376e-15 -7.496312519673904e-01 6.618557139390744e-01 -1.357220357852602e-08 -1.198301226854942e-08 -1.357220357852602e-08 lineNumber 1.753080000000016e+03 date 2012-01-01T12:30:01.12962000000002Z target -4.733734047015606e+06 3.515689711809538e+06 -2.428275079609979e+06 targetDirection 1.110223024625157e-15 -7.496145577933011e-01 6.618746216197999e-01 -1.357231495356307e-08 -1.198371980774984e-08 -1.357231495356307e-08 lineNumber 1.753080000000124e+03 date 2012-01-01T12:30:01.12962000000019Z target -4.733783920649958e+06 3.515638855451895e+06 -2.428266547191616e+06 targetDirection 1.665334536937735e-15 -7.495978631423109e-01 6.618935288794415e-01 -1.357242629275803e-08 -1.198442735245515e-08 -1.357242629275803e-08
sensor datation: sensorName s0 lineNumber 1.000000000000000e+03 date 2012-01-01T12:30:00.00000000000000Z
span: minDate 2012-01-01T12:29:57.50000000000000Z maxDate 2012-01-01T12:30:02.50000000000000Z tStep 1.000000000000000e-03 tolerance 5.000000000000000e+00 inertialFrame EME2000
transform: index 2500 body r -8.085963389171905e-01 -3.465415132416124e-04 4.896468952533136e-04 -5.883634938068593e-01 Ω -8.740475534355121e-08 1.215132763920863e-09 -7.292109805268457e-05 ΩDot -1.642299174832473e-16 8.973031065833714e-17 1.983408395826415e-19 spacecraft p 1.384770770635060e+04 3.157810312896036e+03 -7.179504513258174e+06 v -3.180198832979256e+01 -6.779040795561070e+00 8.276861119339028e+00 a -9.324121873999545e-01 -8.319728238691058e+00 1.345786832372742e-03 r -6.828948903548651e-01 4.142451171315383e-01 -3.878489660867486e-01 4.600312284675194e-01 Ω -1.009782076972198e-03 1.982726233227751e-04 1.647740165426390e-04 ΩDot -3.649972672343429e-07 2.007836645730816e-07 -1.257082123046451e-06
ellipsoid: ae 6.378137000000000e+06 f 3.352810664747481e-03 frame ITRF_CIO_CONV_2010_SIMPLE_EOP
sensor datation: sensorName s0 lineNumber 1.753080000000073e+03 date 2012-01-01T12:30:01.12962000000011Z
transform: index 3630 body r -8.086205790715331e-01 -3.465213511527127e-04 4.896611176157392e-04 -5.883301787746077e-01 Ω -8.740475552913246e-08 1.215132865316375e-09 -7.292109805268537e-05 ΩDot -1.642303315278493e-16 8.973086775115428e-17 1.983413459366029e-19 spacecraft p 1.381161686138390e+04 3.149035127746873e+03 -7.179495156850004e+06 v -3.198309695268937e+01 -7.869713919215563e+00 8.283718407079078e+00 a -9.313885738644379e-01 -8.324782304826440e+00 1.354189237138481e-03 r -6.826577475060918e-01 4.145470548488056e-01 -3.882264160785311e-01 4.597927673908209e-01 Ω -1.009835304084986e-03 1.983572786405474e-04 1.645625945551068e-04 ΩDot -3.648903522451360e-07 2.009701619867746e-07 -1.257044812963575e-06
sensor rate: sensorName s0 lineNumber 1.753080000000073e+03 rate 6.666666666666666e+02
sensor datation: sensorName s0 lineNumber 1.753079999999983e+03 date 2012-01-01T12:30:01.12961999999997Z
sensor LOS: sensorName s0 date 2012-01-01T12:30:01.12961999999997Z pixelNumber 999 los 0.000000000000000e+00 -7.660363365669057e-01 6.427972705753770e-01
sensor LOS: sensorName s0 date 2012-01-01T12:30:01.12961999999997Z pixelNumber 1000 los 0.000000000000000e+00 -7.660525495492120e-01 6.427779486954669e-01
sensor LOS: sensorName s0 date 2012-01-01T12:30:01.12961999999997Z pixelNumber 0 los 0.000000000000000e+00 -7.495978631423039e-01 6.618935288794506e-01
sensor LOS: sensorName s0 date 2012-01-01T12:30:01.12961999999997Z pixelNumber 1 los 0.000000000000000e+00 -7.496145577932898e-01 6.618746216198131e-01
sensor LOS: sensorName s0 date 2012-01-01T12:30:01.12961999999997Z pixelNumber 4 los 0.000000000000000e+00 -7.496646388847948e-01 6.618178973144580e-01
sensor LOS: sensorName s0 date 2012-01-01T12:30:01.12961999999997Z pixelNumber 5 los 0.000000000000000e+00 -7.496813316281102e-01 6.617989883705653e-01
sensor LOS: sensorName s0 date 2012-01-01T12:30:01.12961999999997Z pixelNumber 6 los 0.000000000000000e+00 -7.496980238944815e-01 6.617800790056386e-01
sensor datation: sensorName s0 lineNumber 1.753079999999893e+03 date 2012-01-01T12:30:01.12961999999984Z
inverse location result: lineNumber 1.753079999999893e+03 pixelNumber 5.000000000428289e+00
inverse location: sensorName s0 latitude -3.930323254721546e-01 longitude 2.502727645691323e+00 elevation 1.847490708944211e+03 minLine 0 maxLine 2000 lightTime true aberration true refraction false
sensor datation: sensorName s0 lineNumber 1.753080000000180e+03 date 2012-01-01T12:30:01.12962000000027Z
sensor rate: sensorName s0 lineNumber 1.753080000000180e+03 rate 6.666666666666666e+02
sensor datation: sensorName s0 lineNumber 1.753080000000150e+03 date 2012-01-01T12:30:01.12962000000022Z
sensor LOS: sensorName s0 date 2012-01-01T12:30:01.12962000000022Z pixelNumber 7 los 0.000000000000000e+00 -7.497147156838981e-01 6.617611692196904e-01
sensor datation: sensorName s0 lineNumber 1.753080000000120e+03 date 2012-01-01T12:30:01.12962000000018Z
inverse location result: lineNumber 1.753080000000120e+03 pixelNumber 6.000000000570598e+00
inverse location: sensorName s0 latitude -3.930341476421113e-01 longitude 2.502715716159622e+00 elevation 1.841953904196975e+03 minLine 0 maxLine 2000 lightTime true aberration true refraction false
sensor datation: sensorName s0 lineNumber 1.753080000000450e+03 date 2012-01-01T12:30:01.12962000000067Z
sensor rate: sensorName s0 lineNumber 1.753080000000450e+03 rate 6.666666666666666e+02
sensor LOS: sensorName s0 date 2012-01-01T12:30:01.12962000000027Z pixelNumber 8 los 0.000000000000000e+00 -7.497314069963494e-01 6.617422590127325e-01
sensor datation: sensorName s0 lineNumber 1.753079999999910e+03 date 2012-01-01T12:30:01.12961999999986Z
inverse location result: lineNumber 1.753079999999910e+03 pixelNumber 7.000000000854453e+00
src/test/resources/replay/replay-inverse-loc-03.txt 0 → 100644
View file @ acf902ea
# Rugged library dump file, created on 2015-07-23T14:25:51Z
# all units are SI units (m, m/s, rad ...)
sensor: sensorName s0 nbPixels 2552 position 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
inverse location: sensorName s0 latitude 1.426113542156674e+00 longitude -8.868185262195961e-01 elevation 3.102209319249422e+01 minLine -23040 maxLine 39851 lightTime false aberration false refraction false
sensor mean plane: sensorName s0 minLine -23040 maxLine 39851 maxEval 50 accuracy 1.000000000000000e-02 normal 9.994827482027009e-01 -2.652758464258975e-02 1.818029972819145e-02 cachedResults 6 lineNumber 2.496503931691230e+04 date 2015-07-07T18:40:12.48786757028467Z target 5.829885691015603e+05 -7.150938152185429e+05 6.289925119642745e+06 targetDirection -1.574993670955663e-02 8.911384689865436e-02 9.958969132317702e-01 -2.069121568391233e-07 2.118875364268554e-08 -2.069121568391233e-07 lineNumber 2.498990311831761e+04 date 2015-07-07T18:40:12.52680428328538Z target 5.827831035018093e+05 -7.152819110623585e+05 6.290040934506435e+06 targetDirection -1.574653262138034e-02 8.923613956787135e-02 9.958860166231007e-01 -2.106025594148216e-07 6.293169679808365e-08 -2.106025594148216e-07 lineNumber 2.501446735659665e+04 date 2015-07-07T18:40:12.56527188043036Z target 5.825914009565986e+05 -7.154869255677856e+05 6.290306437975926e+06 targetDirection -1.574243698611968e-02 8.937839743168323e-02 9.958733241483490e-01 -2.139580511938570e-07 1.003165269577007e-07 -2.139580511938570e-07 lineNumber 2.503906484840334e+04 date 2015-07-07T18:40:12.60379155259964Z target 5.823982602071301e+05 -7.156902096509174e+05 6.290556499572599e+06 targetDirection -1.573844749711850e-02 8.951962023030334e-02 9.958607025403716e-01 -2.164310151145369e-07 1.276529299636964e-07 -2.164310151145369e-07 lineNumber 2.506393183649445e+04 date 2015-07-07T18:40:12.64273325595032Z target 5.821928312311379e+05 -7.158784186376551e+05 6.290673076115117e+06 targetDirection -1.573492577246402e-02 8.964463942750402e-02 9.958495120914844e-01 -2.189036239175299e-07 1.552935576937779e-07 -2.189036239175299e-07 lineNumber 2.508894303441162e+04 date 2015-07-07T18:40:12.68190079188860Z target 5.819807540563018e+05 -7.160584640312950e+05 6.290717393595084e+06 targetDirection -1.573177124337613e-02 8.976098395370098e-02 9.958390819470831e-01 -2.207684005657831e-07 1.758718590383612e-07 -2.207684005657831e-07
sensor datation: sensorName s0 lineNumber 8.405500000000000e+03 date 2015-07-07T18:39:46.55562900000000Z
span: minDate 2015-07-07T18:38:55.00000000000000Z maxDate 2015-07-07T18:40:35.80000000000000Z tStep 1.000000000000000e-01 tolerance 1.000000000000000e+01 inertialFrame EME2000
transform: index 516 body r -2.177937215786530e-01 -7.327669555860111e-04 1.897190121950331e-04 -9.759945296466600e-01 Ω -1.102119044807511e-07 3.608696529657919e-09 -7.292106718403156e-05 ΩDot -1.242024347024099e-16 -2.613635776865916e-17 1.864594523542536e-19 spacecraft p -3.566422630398301e+02 -4.258914869799046e+02 -7.158956735101774e+06 v -7.418089019766342e+01 -3.428198127164766e+02 -1.850797604708845e-01 a 1.095670688921161e+02 1.182390694645536e+03 -7.505797508270007e-02 r -6.234703793985741e-02 7.494015157786920e-01 6.578869048593570e-01 4.117080790195870e-02 Ω 1.016965607321450e-03 -1.766672912355467e-04 2.200251772998514e-04 ΩDot -8.541831902061166e-07 1.677502045483661e-04 4.581593874451730e-05
ellipsoid: ae 6.378137000000000e+06 f 3.352810664747481e-03 frame ITRF_CIO_CONV_2010_SIMPLE_EOP
sensor datation: sensorName s0 lineNumber 2.496503925597668e+04 date 2015-07-07T18:40:12.48786747485948Z
transform: index 775 body r -2.187152841723702e-01 -7.325893313267924e-04 1.904097590204008e-04 -9.757884255956278e-01 Ω -1.102119076976303e-07 3.608695853043322e-09 -7.292106718404156e-05 ΩDot -1.242034967780088e-16 -2.611193104220005e-17 1.864622725616803e-19 spacecraft p 5.893360284424853e+01 2.270010215908405e+02 -7.158967334657327e+06 v -1.256875338933382e+02 -2.315544869961368e+02 -6.143872625348799e-01 a -2.991050511747990e+02 -4.922561897467523e+03 -1.631427367003653e-01 r -7.121951955935480e-02 7.504990196615283e-01 6.562822802914443e-01 3.118605611938691e-02 Ω 1.025791280551210e-03 -1.602743286567663e-04 2.237358366388949e-04 ΩDot -5.506323746359040e-05 -6.847191984170615e-04 5.295499556290598e-05
sensor rate: sensorName s0 lineNumber 2.496503925597668e+04 rate 6.385696040868454e+02
sensor datation: sensorName s0 lineNumber 2.496503931691246e+04 date 2015-07-07T18:40:12.48786757028491Z
sensor LOS: sensorName s0 date 2015-07-07T18:40:12.48786757028491Z pixelNumber 1275 los -1.530469679155502e-02 1.052282317151670e-01 9.943303201181274e-01
sensor LOS: sensorName s0 date 2015-07-07T18:40:12.48786757028491Z pixelNumber 1276 los -1.530505436550880e-02 1.052156640641227e-01 9.943316445464334e-01
sensor LOS: sensorName s0 date 2015-07-07T18:40:12.48786757028491Z pixelNumber 0 los -1.481288944757218e-02 1.212648644244146e-01 9.925096528308137e-01
sensor LOS: sensorName s0 date 2015-07-07T18:40:12.48786757028491Z pixelNumber 1 los -1.481330341735178e-02 1.212522760703962e-01 9.925111846168859e-01
sensor LOS: sensorName s0 date 2015-07-07T18:40:12.48786757028491Z pixelNumber 2550 los -1.572509919997303e-02 8.921598849965143e-02 9.958881607144354e-01
sensor LOS: sensorName s0 date 2015-07-07T18:40:12.48786757028491Z pixelNumber 2551 los -1.572540117814331e-02 8.920343838289327e-02 9.958892801604283e-01
sensor datation: sensorName s0 lineNumber 2.496331776312614e+04 date 2015-07-07T18:40:12.48517161705554Z
sensor LOS: sensorName s0 date 2015-07-07T18:40:12.48517161705554Z pixelNumber 2550 los -1.572509919997303e-02 8.921598849965143e-02 9.958881607144354e-01
sensor LOS: sensorName s0 date 2015-07-07T18:40:12.48517161705554Z pixelNumber 2551 los -1.572540117814331e-02 8.920343838289327e-02 9.958892801604283e-01
inverse location result: lineNumber 2.496331776312614e+04 pixelNumber 2.558141579736633e+03
src/tutorials/java/fr/cs/examples/AtmosphericRefractionExamples.java 0 → 100644
View file @ acf902ea
package fr.cs.examples;
import java.io.File;
import java.net.URISyntaxException;
import java.text.DecimalFormat;
import java.text.DecimalFormatSymbols;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.stat.descriptive.DescriptiveStatistics;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataContext;
import org.orekit.data.DirectoryCrawler;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.Transform;
import org.orekit.rugged.api.AlgorithmId;
import org.orekit.rugged.api.BodyRotatingFrameId;
import org.orekit.rugged.api.EllipsoidId;
import org.orekit.rugged.api.InertialFrameId;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.api.RuggedBuilder;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.LinearLineDatation;
import org.orekit.rugged.linesensor.SensorPixel;
import org.orekit.rugged.los.FixedRotation;
import org.orekit.rugged.los.LOSBuilder;
import org.orekit.rugged.los.TimeDependentLOS;
import org.orekit.rugged.refraction.AtmosphericRefraction;
import org.orekit.rugged.refraction.MultiLayerModel;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.AngularDerivativesFilter;
import org.orekit.utils.CartesianDerivativesFilter;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.TimeStampedAngularCoordinates;
import org.orekit.utils.TimeStampedPVCoordinates;
/** A simple example of how to use the atmospheric refraction correction.
*/
public class AtmosphericRefractionExamples {
public static void main(String[] args) throws URISyntaxException {
// 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));
// 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(0., i*FastMath.toRadians(20.)/2000., 1.));
}
// 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, 1., 20.);
// Line sensor
// -----------
// With the LOS and the datation now define, 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.340236, 0.333952, -0.844012, -0.245684);
addSatelliteQ(gps, satelliteQList, "2009-12-11T16:59:06.592937", -0.354773, 0.329336, -0.837871, -0.252281);
addSatelliteQ(gps, satelliteQList, "2009-12-11T16:59:30.592937", -0.369237, 0.324612, -0.831445, -0.258824);
addSatelliteQ(gps, satelliteQList, "2009-12-11T16:59:54.592937", -0.3836, 0.319792, -0.824743, -0.265299);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:00:18.592937", -0.397834, 0.314883, -0.817777, -0.271695);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:00:42.592937", -0.411912, 0.309895, -0.810561, -0.278001);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:01:06.592937", -0.42581, 0.304838, -0.803111, -0.284206);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:01:30.592937", -0.439505, 0.299722, -0.795442, -0.290301);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:01:54.592937", -0.452976, 0.294556, -0.787571, -0.296279);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:02:18.592937", -0.466207, 0.28935, -0.779516, -0.302131);
// Positions and velocities
// ------------------------
ArrayList<TimeStampedPVCoordinates> satellitePVList = new ArrayList<TimeStampedPVCoordinates>();
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:58:42.592937", -726361.466, -5411878.485, 4637549.599, -2463.635, -4447.634, -5576.736);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:04.192937", -779538.267, -5506500.533, 4515934.894, -2459.848, -4312.676, -5683.906);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:25.792937", -832615.368, -5598184.195, 4392036.13, -2454.395, -4175.564, -5788.201);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:47.392937", -885556.748, -5686883.696, 4265915.971, -2447.273, -4036.368, -5889.568);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:00:08.992937", -938326.32, -5772554.875, 4137638.207, -2438.478, -3895.166, -5987.957);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:00:30.592937", -990887.942, -5855155.21, 4007267.717, -2428.011, -3752.034, -6083.317);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:00:52.192937", -1043205.448, -5934643.836, 3874870.441, -2415.868, -3607.05, -6175.600);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:01:13.792937", -1095242.669, -6010981.571, 3740513.34, -2402.051, -3460.291, -6264.760);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:01:35.392937", -1146963.457, -6084130.93, 3604264.372, -2386.561, -3311.835, -6350.751);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:01:56.992937", -1198331.706, -6154056.146, 3466192.446, -2369.401, -3161.764, -6433.531);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:02:18.592937", -1249311.381, -6220723.191, 3326367.397, -2350.574, -3010.159, -6513.056);
// Rugged initialization
// ---------------------
RuggedBuilder builder = new RuggedBuilder().
setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID).
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 Rugged without atmospheric refraction model
Rugged ruggedWithout = builder.build();
// Defines atmospheric refraction model (with the default multi layers model)
AtmosphericRefraction atmosphericRefraction = new MultiLayerModel(ruggedWithout.getEllipsoid());
// One can use its own atmospheric model that must extend the abstract refraction.AtmosphericRefraction class.
// Set the optional grid steps, for the inverse location computation. Useless for direct location computation.
// This setting is optional as default values exist. But can be useful for instance to improve results by given smaller steps
// or to speed up the computation bigger steps.
int pixelStep = 100;
int lineStep = 100;
atmosphericRefraction.setGridSteps(pixelStep, lineStep);
// To compute the inverse location grid with atmospheric refraction, a default margin is available:
// atmosphericRefraction.getComputationParameters().getDefaultInverseLocMargin();
// but can be modified with:
// atmosphericRefraction.setInverseLocMargin(margin);
// Build Rugged with atmospheric refraction model
builder.setRefractionCorrection(atmosphericRefraction);
Rugged ruggedWith = builder.build();
// Direct location on a line WITHOUT and WITH atmospheric correction
// -----------------------------------------------------------------
double chosenLine = 200.;
// Without atmosphere
GeodeticPoint[] gpWithoutAtmosphericRefractionCorrection = ruggedWithout.directLocation("mySensor", chosenLine);
// With atmosphere
GeodeticPoint[] gpWithAtmosphericRefractionCorrection = ruggedWith.directLocation("mySensor", chosenLine);
double earthRadius = ruggedWithout.getEllipsoid().getEquatorialRadius();
DescriptiveStatistics statDistance = new DescriptiveStatistics();
for (int i = 0; i < gpWithAtmosphericRefractionCorrection.length; i++) {
double currentRadius = earthRadius + (gpWithAtmosphericRefractionCorrection[i].getAltitude()+ gpWithoutAtmosphericRefractionCorrection[i].getAltitude())/2.;
double distance = computeDistanceInMeter(currentRadius, gpWithAtmosphericRefractionCorrection[i], gpWithoutAtmosphericRefractionCorrection[i]);
statDistance.addValue(distance);
}
System.out.format("Distance must be > 0 and < 2m:" +
" Max = " + DFS.format(statDistance.getMax()) +
" Min = " + DFS.format(statDistance.getMin()) +
" Median = " + DFS.format(statDistance.getPercentile(50.)) +
" Mean = " + DFS.format(statDistance.getMean()) +
" Std deviation= " + DFS.format(statDistance.getStandardDeviation()) +
"\n");
// Inverse loc WITH atmospheric correction
// ==========================================================================
final double epsilonPixel = 1.e-3;
final double epsilonLine = 1.1e-2;
int minLine = (int) FastMath.floor(lineSensor.getLine(ruggedWithout.getMinDate()));
int maxLine = (int) FastMath.ceil(lineSensor.getLine(ruggedWithout.getMaxDate()));
for (int i = 0; i < gpWithAtmosphericRefractionCorrection.length; i++) {
// to check if we go back to the initial point when taking the geodetic point with atmospheric correction
GeodeticPoint gpWith = gpWithAtmosphericRefractionCorrection[i];
SensorPixel sensorPixelReverseWith = ruggedWith.inverseLocation("mySensor", gpWith, minLine, maxLine);
if (sensorPixelReverseWith != null) {
// Compare the computed pixel vs expected
if (FastMath.abs(i - sensorPixelReverseWith.getPixelNumber()) > epsilonPixel) {
System.out.format("Problem with pixel " + i + " . Delta pixel= " + DFS.format((i - sensorPixelReverseWith.getPixelNumber())) + "\n");
}
// Compare the computed line vs the expected
if (FastMath.abs(chosenLine - sensorPixelReverseWith.getLineNumber())> epsilonLine) {
System.out.format("Probem with line, for pixel " + i + ". Delta line= " + DFS.format((chosenLine - sensorPixelReverseWith.getLineNumber())) + "\n");
}
} else {
System.out.println("Inverse location failed for pixel " + i + " with atmospheric refraction correction for geodetic point computed with" );
}
} // end loop on pixel i
System.out.format("Inverse location gave same pixel and line as direct location input for" +
" epsilon pixel= " + DFS.format(epsilonPixel) +
" and epsilon line= " + DFS.format(epsilonLine) +
"\n");
}
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);
}
private static double computeDistanceInMeter(double earthRadius, final GeodeticPoint gp1, final GeodeticPoint gp2) {
// get vectors on unit sphere from angular coordinates
final Vector3D p1 = new Vector3D(gp1.getLongitude(), gp1.getLatitude());
final Vector3D p2 = new Vector3D(gp2.getLongitude(), gp2.getLatitude());
return earthRadius * Vector3D.angle(p1, p2);
}
/** Definition of print format. */
private static final DecimalFormatSymbols SYMBOLS = new DecimalFormatSymbols(Locale.US);
/** Defined of number of digit after the point. */
private static final DecimalFormat DFS = new DecimalFormat("#.#####", SYMBOLS);
}
src/tutorials/java/fr/cs/examples/DirectLocation.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You 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.
*/
package fr.cs.examples;
import java.io.File;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataContext;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.Transform;
import org.orekit.rugged.api.AlgorithmId;
import org.orekit.rugged.api.BodyRotatingFrameId;
import org.orekit.rugged.api.EllipsoidId;
import org.orekit.rugged.api.InertialFrameId;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.api.RuggedBuilder;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.LinearLineDatation;
import org.orekit.rugged.los.FixedRotation;
import org.orekit.rugged.los.LOSBuilder;
import org.orekit.rugged.los.TimeDependentLOS;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.AngularDerivativesFilter;
import org.orekit.utils.CartesianDerivativesFilter;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.TimeStampedAngularCoordinates;
import org.orekit.utils.TimeStampedPVCoordinates;
public class DirectLocation {
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));
// 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.IGNORE_DEM_USE_ELLIPSOID).
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%n",
FastMath.toDegrees(upLeftPoint.getLatitude()),
FastMath.toDegrees(upLeftPoint.getLongitude()),
upLeftPoint.getAltitude());
} 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);
}
}
src/tutorials/java/fr/cs/examples/DirectLocationWithDEM.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You 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.
*/
package fr.cs.examples;
import java.io.File;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataContext;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.Transform;
import org.orekit.rugged.api.AlgorithmId;
import org.orekit.rugged.api.BodyRotatingFrameId;
import org.orekit.rugged.api.EllipsoidId;
import org.orekit.rugged.api.InertialFrameId;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.api.RuggedBuilder;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.LinearLineDatation;
import org.orekit.rugged.los.FixedRotation;
import org.orekit.rugged.los.LOSBuilder;
import org.orekit.rugged.los.TimeDependentLOS;
import org.orekit.rugged.raster.TileUpdater;
import org.orekit.rugged.raster.UpdatableTile;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.AngularDerivativesFilter;
import org.orekit.utils.CartesianDerivativesFilter;
import org.orekit.utils.Constants;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.TimeStampedAngularCoordinates;
import org.orekit.utils.TimeStampedPVCoordinates;
public class DirectLocationWithDEM {
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));
// 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(new VolcanicConeElevationUpdater(new GeodeticPoint(FastMath.toRadians(37.58),
FastMath.toRadians(-96.95),
2463.0),
FastMath.toRadians(30.0), 16.0,
FastMath.toRadians(1.0), 1201),
8).
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%n",
FastMath.toDegrees(upLeftPoint.getLatitude()),
FastMath.toDegrees(upLeftPoint.getLongitude()),
upLeftPoint.getAltitude());
} 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);
}
private static 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;
}
@Override
public void updateTile(double latitude, double longitude, UpdatableTile tile) {
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);
}
}
}
}
}
src/tutorials/java/fr/cs/examples/InverseLocation.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You 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.
*/
package fr.cs.examples;
import java.io.File;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.data.DataContext;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.Transform;
import org.orekit.rugged.api.AlgorithmId;
import org.orekit.rugged.api.BodyRotatingFrameId;
import org.orekit.rugged.api.EllipsoidId;
import org.orekit.rugged.api.InertialFrameId;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.api.RuggedBuilder;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.LinearLineDatation;
import org.orekit.rugged.linesensor.SensorPixel;
import org.orekit.rugged.los.FixedRotation;
import org.orekit.rugged.los.LOSBuilder;
import org.orekit.rugged.los.TimeDependentLOS;
import org.orekit.rugged.utils.RoughVisibilityEstimator;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.AngularDerivativesFilter;
import org.orekit.utils.CartesianDerivativesFilter;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.TimeStampedAngularCoordinates;
import org.orekit.utils.TimeStampedPVCoordinates;
public class InverseLocation {
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));
// 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:
String sensorName = "mySensor";
LineSensor lineSensor = new LineSensor(sensorName, 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
// ---------------------
RuggedBuilder ruggedBuilder = new RuggedBuilder();
ruggedBuilder.setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID);
ruggedBuilder.setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF);
ruggedBuilder.setTimeSpan(absDate, absDate.shiftedBy(60.0), 0.01, 5 / lineSensor.getRate(0));
ruggedBuilder.setTrajectory(InertialFrameId.EME2000,
satellitePVList, 4, CartesianDerivativesFilter.USE_P,
satelliteQList, 4, AngularDerivativesFilter.USE_R).
addLineSensor(lineSensor);
Rugged rugged = ruggedBuilder.build();
// Inverse location of a Geodetic Point
// ------------------------------------
// Point defined by its latitude, longitude and altitude
double latitude = FastMath.toRadians(37.585);
double longitude = FastMath.toRadians(-96.949);
double altitude = 0.0d;
// For a GeodeticPoint : angles are given in radians and altitude in meters
GeodeticPoint gp = new GeodeticPoint(latitude, longitude, altitude);
// Search the sensor pixel seeing point
// ....................................
// interval of lines where to search the point
int minLine = 0;
int maxLine = 100;
SensorPixel sensorPixel = rugged.inverseLocation(sensorName, gp, minLine, maxLine);
// we need to test if the sensor pixel is found in the prescribed lines otherwise the sensor pixel is null
// Be careful: no exception is thrown if the pixel is not found
if (sensorPixel != null){
// The pixel was found
System.out.format(Locale.US, "Sensor Pixel found : line = %5.3f, pixel = %5.3f %n", sensorPixel.getLineNumber(), sensorPixel.getPixelNumber());
} else {
// The pixel was not found: set to null
System.out.println("Sensor Pixel is null: point cannot be seen between the prescribed line numbers\n");
}
// Find the date at which the sensor sees the ground point
// .......................................................
AbsoluteDate dateLine = rugged.dateLocation(sensorName, gp, minLine, maxLine);
System.out.println("Date at which the sensor sees the ground point : " + dateLine);
// How to find min and max lines ? with the RoughVisibilityEstimator
// -------------------------------
// Create a RoughVisibilityEstimator for inverse location
OneAxisEllipsoid oneAxisEllipsoid = ruggedBuilder.getEllipsoid();
Frame pvFrame = ruggedBuilder.getInertialFrame();
// Initialize the RoughVisibilityEstimator
RoughVisibilityEstimator roughVisibilityEstimator = new RoughVisibilityEstimator(oneAxisEllipsoid, pvFrame, satellitePVList);
// Compute the approximated line with a rough estimator
AbsoluteDate roughLineDate = roughVisibilityEstimator.estimateVisibility(gp);
double roughLine = lineSensor.getLine(roughLineDate);
// Compute the min / max lines interval using a margin around the roughLine
int sensorMinLine= 0;
int sensorMaxLine = 1000;
int margin = 100;
int minLineRough = (int) FastMath.max(roughLine - margin, sensorMinLine);
int maxLineRough = (int) FastMath.min(roughLine + margin, sensorMaxLine);
SensorPixel sensorPixelRoughLine = rugged.inverseLocation(sensorName, gp, minLineRough, maxLineRough);
System.out.format(Locale.US, "Rough line found = %5.1f; InverseLocation gives (margin of %d around rough line): line = %5.3f, pixel = %5.3f %n", roughLine, margin, sensorPixelRoughLine.getLineNumber(), sensorPixel.getPixelNumber());
} 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);
}
}
src/tutorials/java/fr/cs/examples/refiningPleiades/GroundRefining.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You 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.
*/
package fr.cs.examples.refiningPleiades;
import java.io.File;
import java.util.List;
import java.util.Locale;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataContext;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.orbits.Orbit;
import org.orekit.rugged.adjustment.measurements.SensorToGroundMapping;
import org.orekit.rugged.api.AlgorithmId;
import org.orekit.rugged.api.BodyRotatingFrameId;
import org.orekit.rugged.api.EllipsoidId;
import org.orekit.rugged.api.InertialFrameId;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.api.RuggedBuilder;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.AngularDerivativesFilter;
import org.orekit.utils.CartesianDerivativesFilter;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.TimeStampedAngularCoordinates;
import org.orekit.utils.TimeStampedPVCoordinates;
import fr.cs.examples.refiningPleiades.generators.GroundMeasurementGenerator;
import fr.cs.examples.refiningPleiades.generators.Noise;
import fr.cs.examples.refiningPleiades.metrics.DistanceTools;
import fr.cs.examples.refiningPleiades.models.OrbitModel;
import fr.cs.examples.refiningPleiades.models.PleiadesViewingModel;
/**
* Class for testing refining (Ground Control Points GCP study)
* with or without noisy measurements
* @author Jonathan Guinet
* @author Lucie Labat-Allee
* @author Guylaine Prat
* @see SensorToGroundMapping
* @see GroundMeasurementGenerator
* @since 2.0
*/
public class GroundRefining extends Refining {
/** Pleiades viewing model */
private static PleiadesViewingModel pleiadesViewingModel;
/** Sensor name */
private static String sensorName;
/** Main function
*/
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 refining context
// ---------------------------
GroundRefining refining = new GroundRefining();
// Sensor's definition: create Pleiades viewing model
// --------------------------------------------------
System.out.format("**** Build Pleiades viewing model and orbit definition **** %n");
AbsoluteDate minDate = pleiadesViewingModel.getMinDate();
AbsoluteDate maxDate = pleiadesViewingModel.getMaxDate();
AbsoluteDate refDate = pleiadesViewingModel.getDatationReference();
LineSensor lineSensor = pleiadesViewingModel.getLineSensor();
// Satellite position, velocity and attitude: create an orbit model
// ----------------------------------------------------------------
OrbitModel orbitmodel = new OrbitModel();
BodyShape earth = orbitmodel.createEarth();
NormalizedSphericalHarmonicsProvider gravityField = orbitmodel.createGravityField();
Orbit orbit = orbitmodel.createOrbit(gravityField.getMu(), refDate);
// Nadir's pointing
final double [] rollPoly = {0.0,0.0,0.0};
double[] pitchPoly = {0.0,0.0};
double[] yawPoly = {0.0,0.0,0.0};
orbitmodel.setLOFTransform(rollPoly, pitchPoly, yawPoly, minDate);
// Satellite attitude
List<TimeStampedAngularCoordinates> satelliteQList =
orbitmodel.orbitToQ(orbit, earth, minDate.shiftedBy(-0.0), maxDate.shiftedBy(+0.0), 0.25);
int nbQPoints = 2;
// Position and velocities
PVCoordinates PV = orbit.getPVCoordinates(earth.getBodyFrame());
List<TimeStampedPVCoordinates> satellitePVList =
orbitmodel.orbitToPV(orbit, earth, minDate.shiftedBy(-0.0), maxDate.shiftedBy(+0.0), 0.25);
int nbPVPoints = 8;
// Convert frame and coordinates type in one call
GeodeticPoint gp = earth.transform(PV.getPosition(), earth.getBodyFrame(), orbit.getDate());
System.out.format(Locale.US, "Geodetic Point at date %s : φ = %8.10f °, λ = %8.10f %n",
orbit.getDate().toString(),
FastMath.toDegrees(gp.getLatitude()),
FastMath.toDegrees(gp.getLongitude()));
// Rugged initialization
// ---------------------
System.out.format("\n**** Rugged initialization **** %n");
RuggedBuilder ruggedBuilder = new RuggedBuilder();
ruggedBuilder.addLineSensor(lineSensor);
ruggedBuilder.setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID);
ruggedBuilder.setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF);
ruggedBuilder.setTimeSpan(minDate,maxDate, 0.001, 5.0);
ruggedBuilder.setTrajectory(InertialFrameId.EME2000, satellitePVList,nbPVPoints,
CartesianDerivativesFilter.USE_PV, satelliteQList,
nbQPoints, AngularDerivativesFilter.USE_R);
ruggedBuilder.setName("Rugged_refining");
Rugged rugged = ruggedBuilder.build();
// Compute ground sample distance (GSD)
// ------------------------------------
double [] gsd = refining.computeGSD(rugged, lineSensor);
System.out.format("GSD - X: %2.2f Y: %2.2f **** %n", gsd[0], gsd[1]);
// Initialize disruptions:
// -----------------------
// Introduce rotations around instrument axes (roll and pitch translations, scale factor)
System.out.format("\n**** Add disruptions: roll and pitch rotations, scale factor **** %n");
double rollValue = FastMath.toRadians(-0.01);
double pitchValue = FastMath.toRadians(0.02);
double factorValue = 1.05;
System.out.format("roll: %3.5f \tpitch: %3.5f \tfactor: %3.5f \n",rollValue, pitchValue, factorValue);
// Apply disruptions on physical model
refining.applyDisruptions(rugged, sensorName, rollValue, pitchValue, factorValue);
// Generate measurements (observations) from physical model disrupted
// ------------------------------------------------------------------
int lineSampling = 1000;
int pixelSampling = 1000;
// Noise definition
// distribution: gaussian(0), vector dimension:3
final Noise noise = new Noise(0,3);
// lat mean, long mean, alt mean
final double[] mean = {0.0,0.0,0.0};
// lat std, long std, alt std}
final double[] standardDeviation = {0.0,0.0,0.0};
noise.setMean(mean);
noise.setStandardDeviation(standardDeviation);
// Ground measurements
GroundMeasurementGenerator measurements = refining.generateNoisyPoints(lineSampling, pixelSampling,
rugged, sensorName,
pleiadesViewingModel.getDimension(),
noise);
// Compute ground truth residuals
// ------------------------------
System.out.format("\n**** Ground truth residuals **** %n");
refining.computeMetrics(measurements.getGroundMapping(), rugged, false);
// Initialize physical model without disruptions
// ---------------------------------------------
System.out.format("\n**** Initialize physical model without disruptions: reset Roll/Pitch/Factor **** %n");
refining.resetModel(rugged, sensorName, true);
// Compute initial residuals
// -------------------------
System.out.format("\n**** Initial Residuals **** %n");
refining.computeMetrics(measurements.getGroundMapping(), rugged, false);
// Start optimization
// ------------------
System.out.format("\n**** Start optimization **** %n");
int maxIterations = 100;
double convergenceThreshold = 1.e-11;
refining.optimization(maxIterations, convergenceThreshold, measurements.getObservables(), rugged);
// Check estimated values
// ----------------------
System.out.format("\n**** Check parameters ajustement **** %n");
refining.paramsEstimation(rugged, sensorName, rollValue, pitchValue, factorValue);
// Compute statistics
// ------------------
System.out.format("\n**** Compute Statistics **** %n");
// Residuals computed in meters
refining.computeMetrics(measurements.getGroundMapping(), rugged, false);
// Residuals computed in degrees
refining.computeMetrics(measurements.getGroundMapping(), rugged, true);
} catch (OrekitException oe) {
System.err.println(oe.getLocalizedMessage());
System.exit(1);
} catch (RuggedException re) {
System.err.println(re.getLocalizedMessage());
System.exit(1);
}
}
/** Constructor */
private GroundRefining() {
sensorName = "line";
pleiadesViewingModel = new PleiadesViewingModel(sensorName);
}
/** Estimate ground sample distance (GSD)
* @param LineSensor line sensor
* @return the GSD
*/
private double[] computeGSD(final Rugged rugged, final LineSensor lineSensor) {
// Get number of line
int dimension = pleiadesViewingModel.getDimension();
// Get position
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.
// Get upper left geodetic point
AbsoluteDate firstLineDate = lineSensor.getDate(0);
Vector3D los = lineSensor.getLOS(firstLineDate,0);
GeodeticPoint upLeftPoint = rugged.directLocation(firstLineDate, position, los);
los = lineSensor.getLOS(firstLineDate,dimension-1);
// Get center geodetic point
AbsoluteDate lineDate = lineSensor.getDate(dimension/2);
los = lineSensor.getLOS(lineDate,dimension/2);
// Get upper right geodetic point
int pixelPosition = dimension-1;
los = lineSensor.getLOS(firstLineDate,pixelPosition);
GeodeticPoint upperRight = rugged.directLocation(firstLineDate, position, los);
// Get lower left geodetic point
AbsoluteDate lineDate_y = lineSensor.getDate(dimension-1);
los = lineSensor.getLOS(lineDate_y,0);
GeodeticPoint lowerLeft = rugged.directLocation(lineDate_y, position, los);
double gsdX = DistanceTools.computeDistanceInMeter(upLeftPoint, upperRight)/dimension;
double gsdY = DistanceTools.computeDistanceInMeter(upLeftPoint, lowerLeft)/dimension;
double [] gsd = {gsdX, gsdY};
return gsd;
}
// /**
// * Get the Pleiades viewing model
// * @return the Pleiades viewing model
// */
// public PleiadesViewingModel getPleiadesViewingModel() {
// return pleiadesViewingModel;
// }
//
// /**
// * Set the Pleiades viewing model
// * @param pleiadesViewingModel Pleiades viewing model to set
// */
// public void setPleiadesViewingModel(final PleiadesViewingModel pleiadesViewingModel) {
// this.pleiadesViewingModel = pleiadesViewingModel;
// }
//
// /**
// * Get the orbit model
// * @return the orbit model
// */
// public OrbitModel getOrbitmodel() {
// return orbitmodel;
// }
//
// /**
// * Set the orbit model
// * @param orbitmodel the orbit model to set
// */
// public void setOrbitmodel(final OrbitModel orbitmodel) {
// this.orbitmodel = orbitmodel;
// }
//
// /**
// * Get the sensor name
// * @return the sensor name
// */
// public String getSensorName() {
// return sensorName;
// }
//
// /**
// * Get the Rugged instance
// * @return the rugged instance
// */
// public Rugged getRugged() {
// return rugged;
// }
//
// /**
// * Set the Rugged instance
// * @param rugged the Rugged instance to set
// */
// public void setRugged(final Rugged rugged) {
// this.rugged = rugged;
// }
}
src/tutorials/java/fr/cs/examples/refiningPleiades/InterRefining.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You 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.
*/
package fr.cs.examples.refiningPleiades;
import java.io.File;
import java.util.Arrays;
import java.util.List;
import java.util.Locale;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataContext;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.orbits.Orbit;
import org.orekit.rugged.adjustment.measurements.SensorToGroundMapping;
import org.orekit.rugged.api.AlgorithmId;
import org.orekit.rugged.api.BodyRotatingFrameId;
import org.orekit.rugged.api.EllipsoidId;
import org.orekit.rugged.api.InertialFrameId;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.api.RuggedBuilder;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.AngularDerivativesFilter;
import org.orekit.utils.CartesianDerivativesFilter;
import org.orekit.utils.Constants;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.TimeStampedAngularCoordinates;
import org.orekit.utils.TimeStampedPVCoordinates;
import fr.cs.examples.refiningPleiades.generators.GroundMeasurementGenerator;
import fr.cs.examples.refiningPleiades.generators.InterMeasurementGenerator;
import fr.cs.examples.refiningPleiades.generators.Noise;
import fr.cs.examples.refiningPleiades.metrics.DistanceTools;
import fr.cs.examples.refiningPleiades.models.OrbitModel;
import fr.cs.examples.refiningPleiades.models.PleiadesViewingModel;
/**
* Class for testing refining (tie points study)
* with or without noisy measurements
* @author Jonathan Guinet
* @author Lucie Labat-Allee
* @author Guylaine Prat
* @see SensorToGroundMapping
* @see GroundMeasurementGenerator
* @since 2.0
*/
public class InterRefining extends Refining {
/** Pleiades viewing model A */
private static PleiadesViewingModel pleiadesViewingModelA;
/** Pleiades viewing model B */
private static PleiadesViewingModel pleiadesViewingModelB;
/** Sensor name A corresponding to viewing model A */
private static String sensorNameA;
/** Sensor name A corresponding to viewing model B */
private static String sensorNameB;
/** Main function
*/
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 refining context
// ---------------------------
InterRefining refining = new InterRefining();
// Sensors's definition: creation of 2 Pleiades viewing models
// -----------------------------------------------------------
System.out.format("**** Build Pleiades viewing model A and its orbit definition **** %n");
// 1/- Create First Pleiades Viewing Model
final AbsoluteDate minDateA = pleiadesViewingModelA.getMinDate();
final AbsoluteDate maxDateA = pleiadesViewingModelA.getMaxDate();
final AbsoluteDate refDateA = pleiadesViewingModelA.getDatationReference();
LineSensor lineSensorA = pleiadesViewingModelA.getLineSensor();
System.out.format(Locale.US, "Viewing model A - date min: %s max: %s ref: %s %n", minDateA, maxDateA, refDateA);
// ----Satellite position, velocity and attitude: create orbit model A
OrbitModel orbitmodelA = new OrbitModel();
BodyShape earthA = orbitmodelA.createEarth();
Orbit orbitA = orbitmodelA.createOrbit(Constants.EIGEN5C_EARTH_MU, refDateA);
// ----If no LOF Transform Attitude Provider is Nadir Pointing Yaw Compensation
final double [] rollPoly = {0.0,0.0,0.0};
double[] pitchPoly = {0.025,0.0};
double[] yawPoly = {0.0,0.0,0.0};
orbitmodelA.setLOFTransform(rollPoly, pitchPoly, yawPoly, minDateA);
// ----Satellite attitude
List<TimeStampedAngularCoordinates> satelliteQListA = orbitmodelA.orbitToQ(orbitA, earthA, minDateA.shiftedBy(-0.0), maxDateA.shiftedBy(+0.0), 0.25);
int nbQPoints = 2;
// ----Position and velocities
PVCoordinates PVA = orbitA.getPVCoordinates(earthA.getBodyFrame());
List<TimeStampedPVCoordinates> satellitePVListA = orbitmodelA.orbitToPV(orbitA, earthA, minDateA.shiftedBy(-0.0), maxDateA.shiftedBy(+0.0), 0.25);
int nbPVPoints = 8;
// ----Convert frame and coordinates type in one call
GeodeticPoint gpA = earthA.transform(PVA.getPosition(), earthA.getBodyFrame(), orbitA.getDate());
System.out.format(Locale.US, "(A) Geodetic Point at date %s : φ = %8.10f °, λ = %8.10f %n",
orbitA.getDate().toString(),
FastMath.toDegrees(gpA.getLatitude()),
FastMath.toDegrees(gpA.getLongitude()));
// Rugged A initialization
// ---------------------
System.out.format("**** Rugged A initialization **** %n");
RuggedBuilder ruggedBuilderA = new RuggedBuilder();
ruggedBuilderA.addLineSensor(lineSensorA);
ruggedBuilderA.setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID);
ruggedBuilderA.setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF);
ruggedBuilderA.setTimeSpan(minDateA,maxDateA, 0.001, 5.0);
ruggedBuilderA.setTrajectory(InertialFrameId.EME2000, satellitePVListA, nbPVPoints,
CartesianDerivativesFilter.USE_PV, satelliteQListA,
nbQPoints, AngularDerivativesFilter.USE_R);
ruggedBuilderA.setLightTimeCorrection(false);
ruggedBuilderA.setAberrationOfLightCorrection(false);
ruggedBuilderA.setName("RuggedA");
Rugged ruggedA = ruggedBuilderA.build();
System.out.format("\n**** Build Pleiades viewing model B and its orbit definition **** %n");
// 2/- Create Second Pleiades Viewing Model
final AbsoluteDate minDateB = pleiadesViewingModelB.getMinDate();
final AbsoluteDate maxDateB = pleiadesViewingModelB.getMaxDate();
final AbsoluteDate refDateB = pleiadesViewingModelB.getDatationReference();
LineSensor lineSensorB = pleiadesViewingModelB.getLineSensor();
System.out.format(Locale.US, "Viewing model B - date min: %s max: %s ref: %s %n", minDateB, maxDateB, refDateB);
// ----Satellite position, velocity and attitude: create orbit model B
OrbitModel orbitmodelB = new OrbitModel();
BodyShape earthB = orbitmodelB.createEarth();
Orbit orbitB = orbitmodelB.createOrbit(Constants.EIGEN5C_EARTH_MU, refDateB);
// ----Satellite attitude
List<TimeStampedAngularCoordinates> satelliteQListB = orbitmodelB.orbitToQ(orbitB, earthB, minDateB.shiftedBy(-0.0), maxDateB.shiftedBy(+0.0), 0.25);
// ----Position and velocities
PVCoordinates PVB = orbitB.getPVCoordinates(earthB.getBodyFrame());
List<TimeStampedPVCoordinates> satellitePVListB = orbitmodelB.orbitToPV(orbitB, earthB, minDateB.shiftedBy(-0.0), maxDateB.shiftedBy(+0.0), 0.25);
GeodeticPoint gpB = earthB.transform(PVB.getPosition(), earthB.getBodyFrame(), orbitB.getDate());
System.out.format(Locale.US, "(B) Geodetic Point at date %s : φ = %8.10f °, λ = %8.10f %n",orbitA.getDate().toString(),
FastMath.toDegrees(gpB.getLatitude()),
FastMath.toDegrees(gpB.getLongitude()));
// Rugged B initialization
// ---------------------
System.out.format("**** Rugged B initialization **** %n");
RuggedBuilder ruggedBuilderB = new RuggedBuilder();
ruggedBuilderB.addLineSensor(lineSensorB);
ruggedBuilderB.setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID);
ruggedBuilderB.setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF);
ruggedBuilderB.setTimeSpan(minDateB,maxDateB, 0.001, 5.0);
ruggedBuilderB.setTrajectory(InertialFrameId.EME2000, satellitePVListB, nbPVPoints,
CartesianDerivativesFilter.USE_PV, satelliteQListB,
nbQPoints, AngularDerivativesFilter.USE_R);
ruggedBuilderB.setLightTimeCorrection(false);
ruggedBuilderB.setAberrationOfLightCorrection(false);
ruggedBuilderB.setName("RuggedB");
Rugged ruggedB = ruggedBuilderB.build();
// Compute distance between LOS
// -----------------------------
double distance = refining.computeDistanceBetweenLOS(ruggedA, lineSensorA, ruggedB, lineSensorB);
System.out.format("distance %f meters %n",distance);
// Initialize disruptions:
// -----------------------
// Introduce rotations around instrument axes (roll and pitch translations, scale factor)
System.out.format("\n**** Add disruptions on both acquisitions (A,B): roll and pitch rotations, scale factor **** %n");
double rollValueA = FastMath.toRadians(0.004);
double pitchValueA = FastMath.toRadians(0.0008);
double rollValueB = FastMath.toRadians(-0.004);
double pitchValueB = FastMath.toRadians(0.0008);
double factorValue = 1.0;
System.out.format("Acquisition A roll: %3.5f \tpitch: %3.5f \tfactor: %3.5f \n",rollValueA,pitchValueA,factorValue);
System.out.format("Acquisition B roll: %3.5f \tpitch: %3.5f \tfactor: %3.5f \n",rollValueB,pitchValueB,factorValue);
// Apply disruptions on physical model (acquisition A)
refining.applyDisruptions(ruggedA, sensorNameA, rollValueA, pitchValueA, factorValue);
// Apply disruptions on physical model (acquisition B)
refining.applyDisruptions(ruggedB, sensorNameB, rollValueB, pitchValueB, factorValue);
// Generate measurements (observations) from physical model disrupted
// ------------------------------------------------------------------
int lineSampling = 1000;
int pixelSampling = 1000;
// Noise definition
// distribution: gaussian(0), vector dimension: 2
final Noise noise = new Noise(0,2);
// {pixelA mean, pixelB mean}
final double[] mean = {5.0,0.0};
// {pixelB std, pixelB std}
final double[] standardDeviation = {0.1,0.1};
noise.setMean(mean);
noise.setStandardDeviation(standardDeviation);
System.out.format("\tSensor A mean: %f std %f %n",mean[0],standardDeviation[0]);
System.out.format("\tSensor B mean: %f std %f %n",mean[1],standardDeviation[1]);
// Inter-measurements (between both viewing models)
InterMeasurementGenerator measurements = refining.generateNoisyPoints(lineSampling, pixelSampling,
ruggedA, sensorNameA,
pleiadesViewingModelA.getDimension(),
ruggedB, sensorNameB,
pleiadesViewingModelB.getDimension(),
noise);
// Compute ground truth residuals
// ------------------------------
System.out.format("\n**** Ground truth residuals **** %n");
refining.computeMetrics(measurements.getInterMapping(), ruggedA, ruggedB, false);
// Initialize physical models without disruptions
// -----------------------------------------------
System.out.format("\n**** Initialize physical models (A,B) without disruptions: reset Roll/Pitch/Factor **** %n");
refining.resetModel(ruggedA, sensorNameA, false);
refining.resetModel(ruggedB, sensorNameB, false);
// Compute initial residuals
// -------------------------
System.out.format("\n**** Initial Residuals **** %n");
refining.computeMetrics(measurements.getInterMapping(), ruggedA, ruggedB, false);
// Start optimization
// ------------------
System.out.format("\n**** Start optimization **** %n");
final int maxIterations = 100;
final double convergenceThreshold = 1.e-7;
List<Rugged> ruggedList = Arrays.asList(ruggedA, ruggedB);
refining.optimization(maxIterations, convergenceThreshold,
measurements.getObservables(), ruggedList);
// Check estimated values
// ----------------------
System.out.format("\n**** Check parameters ajustement **** %n");
System.out.format("Acquisition A:%n");
refining.paramsEstimation(ruggedA, sensorNameA, rollValueA, pitchValueA, factorValue);
System.out.format("Acquisition B:%n");
refining.paramsEstimation(ruggedB, sensorNameB, rollValueB, pitchValueB, factorValue);
// Compute statistics
// ------------------
System.out.format("\n**** Compute Statistics **** %n");
refining.computeMetrics(measurements.getInterMapping(), ruggedA, ruggedB, false);
} catch (OrekitException oe) {
System.err.println(oe.getLocalizedMessage());
System.exit(1);
} catch (RuggedException re) {
System.err.println(re.getLocalizedMessage());
System.exit(1);
}
}
/** Constructor */
private InterRefining() {
sensorNameA = "SensorA";
final double incidenceAngleA = -5.0;
final String dateA = "2016-01-01T11:59:50.0";
pleiadesViewingModelA = new PleiadesViewingModel(sensorNameA, incidenceAngleA, dateA);
sensorNameB = "SensorB";
final double incidenceAngleB = 0.0;
final String dateB = "2016-01-01T12:02:50.0";
pleiadesViewingModelB = new PleiadesViewingModel(sensorNameB, incidenceAngleB, dateB);
}
/** Estimate distance between LOS
* @param lineSensorA line sensor A
* @param lineSensorB line sensor B
* @return GSD
*/
private double computeDistanceBetweenLOS(final Rugged ruggedA, final LineSensor lineSensorA,
final Rugged ruggedB, final LineSensor lineSensorB) {
// Get number of line of sensors
int dimensionA = pleiadesViewingModelA.getDimension();
int dimensionB = pleiadesViewingModelB.getDimension();
Vector3D positionA = lineSensorA.getPosition();
// This returns a zero vector since we set the relative position of the sensor w.r.T the satellite to 0.
AbsoluteDate lineDateA = lineSensorA.getDate(dimensionA/2);
Vector3D losA = lineSensorA.getLOS(lineDateA,dimensionA/2);
GeodeticPoint centerPointA = ruggedA.directLocation(lineDateA, positionA, losA);
System.out.format(Locale.US, "\ncenter geodetic position A : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n",
FastMath.toDegrees(centerPointA.getLatitude()),
FastMath.toDegrees(centerPointA.getLongitude()),centerPointA.getAltitude());
Vector3D positionB = lineSensorB.getPosition();
// This returns a zero vector since we set the relative position of the sensor w.r.T the satellite to 0.
AbsoluteDate lineDateB = lineSensorB.getDate(dimensionB/2);
Vector3D losB = lineSensorB.getLOS(lineDateB,dimensionB/2);
GeodeticPoint centerPointB = ruggedB.directLocation(lineDateB, positionB, losB);
System.out.format(Locale.US, "center geodetic position B : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n",
FastMath.toDegrees(centerPointB.getLatitude()),
FastMath.toDegrees(centerPointB.getLongitude()),centerPointB.getAltitude());
lineDateB = lineSensorB.getDate(0);
losB = lineSensorB.getLOS(lineDateB,0);
GeodeticPoint firstPointB = ruggedB.directLocation(lineDateB, positionB, losB);
System.out.format(Locale.US, "first geodetic position B : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n",
FastMath.toDegrees(firstPointB.getLatitude()),
FastMath.toDegrees(firstPointB.getLongitude()),firstPointB.getAltitude());
lineDateB = lineSensorB.getDate(dimensionB-1);
losB = lineSensorB.getLOS(lineDateB,dimensionB-1);
GeodeticPoint lastPointB = ruggedB.directLocation(lineDateB, positionB, losB);
System.out.format(Locale.US, "last geodetic position B : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n",
FastMath.toDegrees(lastPointB.getLatitude()),
FastMath.toDegrees(lastPointB.getLongitude()),lastPointB.getAltitude());
double distance = DistanceTools.computeDistanceInMeter(centerPointA, centerPointB);
return distance;
}
}
src/tutorials/java/fr/cs/examples/refiningPleiades/Refining.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You 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.
*/
package fr.cs.examples.refiningPleiades;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
import java.util.List;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer.Optimum;
import org.orekit.rugged.adjustment.AdjustmentContext;
import org.orekit.rugged.adjustment.measurements.Observables;
import org.orekit.rugged.adjustment.measurements.SensorToGroundMapping;
import org.orekit.rugged.adjustment.measurements.SensorToSensorMapping;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.errors.RuggedMessages;
import fr.cs.examples.refiningPleiades.generators.GroundMeasurementGenerator;
import fr.cs.examples.refiningPleiades.generators.InterMeasurementGenerator;
import fr.cs.examples.refiningPleiades.generators.Noise;
import fr.cs.examples.refiningPleiades.metrics.LocalisationMetrics;
/**
* Class for refining problems common methods
* @author Jonathan Guinet
* @author Lucie Labat-Allee
* @author Guylaine Prat
* @see SensorToGroundMapping
* @see SensorToSensorMapping
* @see GroundMeasurementGenerator
* @see InterMeasurementGenerator
* @since 2.0
*/
public class Refining {
/**
* Constructor
*/
public Refining() {
}
/** Apply disruptions on acquisition for roll, pitch and scale factor
* @param rugged Rugged instance
* @param sensorName line sensor name
* @param rollValue rotation on roll value
* @param pitchValue rotation on pitch value
* @param factorValue scale factor
*/
public void applyDisruptions(final Rugged rugged, final String sensorName,
final double rollValue, final double pitchValue, final double factorValue) {
rugged.
getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(sensorName + rollSuffix)).
findFirst().get().setValue(rollValue);
rugged.
getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(sensorName + pitchSuffix)).
findFirst().get().setValue(pitchValue);
rugged.
getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(factorName)).
findFirst().get().setValue(factorValue);
}
/** Generate measurements without noise (sensor to ground mapping)
* @param lineSampling line sampling
* @param pixelSampling pixel sampling
* @param rugged Rugged instance
* @param sensorName line sensor name
* @param dimension number of line of the sensor
* @return ground measurements generator (sensor to ground mapping)
*/
public GroundMeasurementGenerator generatePoints(final int lineSampling, final int pixelSampling,
final Rugged rugged, final String sensorName,
final int dimension) {
GroundMeasurementGenerator measurements = new GroundMeasurementGenerator(rugged, sensorName, dimension);
System.out.format("\n**** Generate measurements (without noise; sensor to ground mapping) **** %n");
// Generation measurements without noise
measurements.createMeasurement(lineSampling, pixelSampling);
System.out.format("Number of tie points generated: %d %n", measurements.getMeasurementCount());
return measurements;
}
/** Generate measurements without noise (sensor to sensor mapping)
* @param lineSampling line sampling
* @param pixelSampling pixel sampling
* @param ruggedA Rugged instance of acquisition A
* @param sensorNameA line sensor name A
* @param dimensionA dimension for acquisition A
* @param ruggedB Rugged instance of acquisition B
* @param sensorNameB line sensor name B
* @param dimensionB dimension for acquisition B
* @return inter measurements generator (sensor to sensor mapping)
*/
public InterMeasurementGenerator generatePoints(final int lineSampling, final int pixelSampling,
final Rugged ruggedA, final String sensorNameA, final int dimensionA,
final Rugged ruggedB, final String sensorNameB, final int dimensionB) {
// Outliers control
final double outlierValue = 1e+2;
// Earth constraint weight
final double earthConstraintWeight = 0.1;
// Generate measurements with constraints on outliers control and Earth distance
InterMeasurementGenerator measurements = new InterMeasurementGenerator(ruggedA, sensorNameA, dimensionA,
ruggedB, sensorNameB, dimensionB,
outlierValue,
earthConstraintWeight);
System.out.format("\n**** Generate measurements (without noise; sensor to sensor mapping) **** %n");
// Generation measurements without noise
measurements.createMeasurement(lineSampling, pixelSampling);
System.out.format("Number of tie points generated: %d %n", measurements.getMeasurementCount());
return measurements;
}
/** Generate noisy measurements (sensor to ground mapping)
* @param lineSampling line sampling
* @param pixelSampling pixel sampling
* @param rugged Rugged instance
* @param sensorName line sensor name
* @param dimension dimension
* @param noise Noise structure to generate noisy measurements
* @return ground measurements generator (sensor to ground mapping)
*/
public GroundMeasurementGenerator generateNoisyPoints(final int lineSampling, final int pixelSampling,
final Rugged rugged, final String sensorName, final int dimension,
final Noise noise) {
// Generate ground measurements
GroundMeasurementGenerator measurements = new GroundMeasurementGenerator(rugged, sensorName, dimension);
System.out.format("\n**** Generate noisy measurements (sensor to ground mapping) **** %n");
// Generate noisy measurements
measurements.createNoisyMeasurement(lineSampling, pixelSampling, noise);
System.out.format("Number of tie points generated: %d %n", measurements.getMeasurementCount());
return measurements;
}
/** Generate noisy measurements (sensor to sensor mapping)
* @param lineSampling line sampling
* @param pixelSampling pixel sampling
* @param ruggedA Rugged instance of acquisition A
* @param sensorNameA line sensor name A
* @param dimensionA dimension for acquisition A
* @param ruggedB Rugged instance of acquisition B
* @param sensorNameB line sensor name B
* @param dimensionB dimension for acquisition B
* @param noise noise structure to generate noisy measurements
* @return inter-measurements generator (sensor to sensor mapping)
*/
public InterMeasurementGenerator generateNoisyPoints(final int lineSampling, final int pixelSampling,
final Rugged ruggedA, final String sensorNameA, final int dimensionA,
final Rugged ruggedB, final String sensorNameB, final int dimensionB,
final Noise noise) {
// Outliers control
final double outlierValue = 1.e+2;
// Earth constraint weight
final double earthConstraintWeight = 0.1;
// Generate measurements with constraints on Earth distance and outliers control
InterMeasurementGenerator measurements = new InterMeasurementGenerator(ruggedA, sensorNameA, dimensionA,
ruggedB, sensorNameB, dimensionB,
outlierValue,
earthConstraintWeight);
System.out.format("\n**** Generate noisy measurements (sensor to sensor mapping) **** %n");
// Generation noisy measurements
measurements.createNoisyMeasurement(lineSampling, pixelSampling, noise);
System.out.format("Number of tie points generated: %d %n", measurements.getMeasurementCount());
return measurements;
}
/** Compute metrics to evaluate geometric performances in location,
* for Ground Control Points GCP study.
* @param groundMapping sensor to ground mapping
* @param rugged Rugged instance
* @param unit flag to know if distance is computed in meters (false) or with angular (true)
*/
public void computeMetrics(final SensorToGroundMapping groundMapping,
final Rugged rugged, final boolean unit) {
String stUnit = null;
if(unit) {
stUnit="degrees";
} else {
stUnit="meters";
}
LocalisationMetrics residuals = new LocalisationMetrics(groundMapping, rugged, unit);
System.out.format("Max: %3.4e Mean: %3.4e %s%n",residuals.getMaxResidual(),residuals.getMeanResidual(), stUnit);
}
/** Compute metrics to evaluate geometric performances in location,
* for tie points study.
* @param interMapping sensor to sensor mapping
* @param ruggedA Rugged instance A
* @param ruggedB Rugged instance B
* @param unit flag to know if distance is computed in meters (false) or with angular (true)
*/
public void computeMetrics(final SensorToSensorMapping interMapping,
final Rugged ruggedA, final Rugged ruggedB,
final boolean unit) {
String stUnit = null;
if(unit) stUnit="degrees";
else stUnit="meters";
LocalisationMetrics residuals = new LocalisationMetrics(interMapping, ruggedA, ruggedB, unit);
System.out.format("Max: %1.4e Mean: %1.4e %s%n",residuals.getMaxResidual(),residuals.getMeanResidual(), stUnit);
System.out.format("LOS distance Max: %1.4e Mean: %1.4e %s%n",residuals.getLosMaxDistance(),residuals.getLosMeanDistance(), stUnit);
System.out.format("Earth distance Max: %1.4e Mean: %1.4e %s%n",residuals.getEarthMaxDistance(),residuals.getEarthMeanDistance(), stUnit);
}
/** Reset a model
* @param rugged Rugged instance
* @param sensorName line sensor name
* @param isSelected flag to known if factor parameter is selected or not
*/
public void resetModel(final Rugged rugged, final String sensorName, final boolean isSelected) {
rugged.
getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(sensorName + rollSuffix)
|| driver.getName().equals(sensorName + pitchSuffix)).
forEach(driver -> {
driver.setSelected(true);
driver.setValue(0.0);
});
rugged.
getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(factorName)).
forEach(driver -> {
driver.setSelected(isSelected);
// default value: no Z scale factor applied
driver.setValue(1.0);
});
}
/** Start optimization to adjust parameters (Ground Control Points GCP study).
* @param maxIterations iterations max
* @param convergenceThreshold threshold of convergence
* @param measurements ground measurements
* @param rugged Rugged instance
*/
public void optimization(final int maxIterations, final double convergenceThreshold,
final Observables measurements, final Rugged rugged) {
System.out.format("Iterations max: %d\tconvergence threshold: %3.6e \n", maxIterations, convergenceThreshold);
AdjustmentContext adjustmentContext = new AdjustmentContext(Collections.singletonList(rugged), measurements);
Optimum optimum = adjustmentContext.estimateFreeParameters(Collections.singletonList(rugged.getName()), maxIterations, convergenceThreshold);
// Print statistics
System.out.format("Max value: %3.6e %n", optimum.getResiduals().getMaxValue());
System.out.format("Optimization performed in %d iterations \n", optimum.getIterations());
System.out.format("Optimization performed with %d evaluations of model (objective) function \n", optimum.getEvaluations());
System.out.format("RMSE: %f \n", optimum.getRMS());
}
/** Start optimization to adjust parameters (tie points study).
* @param maxIterations iterations max
* @param convergenceThreshold threshold of convergence
* @param measurements measurements
* @param ruggeds Rugged instances A and B
*/
public void optimization(final int maxIterations, final double convergenceThreshold,
final Observables measurements,
final Collection<Rugged> ruggeds) {
System.out.format("Iterations max: %d\tconvergence threshold: %3.6e \n", maxIterations, convergenceThreshold);
if(ruggeds.size()!= 2 ) {
throw new RuggedException(RuggedMessages.UNSUPPORTED_REFINING_CONTEXT,ruggeds.size());
}
AdjustmentContext adjustmentContext = new AdjustmentContext(ruggeds, measurements);
List<String> ruggedNameList = new ArrayList<String>();
for(Rugged rugged : ruggeds) {
ruggedNameList.add(rugged.getName());
}
Optimum optimum = adjustmentContext.estimateFreeParameters(ruggedNameList, maxIterations, convergenceThreshold);
// Print statistics
System.out.format("Max value: %3.6e %n", optimum.getResiduals().getMaxValue());
System.out.format("Optimization performed in %d iterations \n", optimum.getIterations());
System.out.format("Optimization performed with %d evaluations of model (objective) function \n", optimum.getEvaluations());
System.out.format("RMSE: %f \n", optimum.getRMS());
}
/** Check adjusted parameters of an acquisition
* @param rugged Rugged instance
* @param sensorName line sensor name
* @param rollValue rotation on roll value
* @param pitchValue rotation on pitch value
* @param factorValue scale factor
*/
public void paramsEstimation(final Rugged rugged, final String sensorName,
final double rollValue, final double pitchValue, final double factorValue) {
// Estimate Roll
double estimatedRoll = rugged.getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(sensorName + rollSuffix)).
findFirst().get().getValue();
double rollError = (estimatedRoll - rollValue);
System.out.format("Estimated roll: %3.5f\troll error: %3.6e %n", estimatedRoll, rollError);
// Estimate pitch
double estimatedPitch = rugged.getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(sensorName + pitchSuffix)).
findFirst().get().getValue();
double pitchError = (estimatedPitch - pitchValue);
System.out.format("Estimated pitch: %3.5f\tpitch error: %3.6e %n", estimatedPitch, pitchError);
// Estimate scale factor
double estimatedFactor = rugged.getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(factorName)).
findFirst().get().getValue();
double factorError = (estimatedFactor - factorValue);
System.out.format("Estimated factor: %3.5f\tfactor error: %3.6e %n", estimatedFactor, factorError);
}
/**
* Part of the name of parameter drivers
*/
static final String rollSuffix = "_roll";
static final String pitchSuffix = "_pitch";
static final String factorName = "factor";
public static String getRollsuffix() {
return rollSuffix;
}
public static String getPitchsuffix() {
return pitchSuffix;
}
public static String getFactorname() {
return factorName;
}
}
src/tutorials/java/fr/cs/examples/refiningPleiades/generators/GroundMeasurementGenerator.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You 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.
*/
package fr.cs.examples.refiningPleiades.generators;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.random.GaussianRandomGenerator;
import org.hipparchus.random.UncorrelatedRandomVectorGenerator;
import org.hipparchus.random.Well19937a;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.adjustment.measurements.Observables;
import org.orekit.rugged.adjustment.measurements.SensorToGroundMapping;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.SensorPixel;
import org.orekit.time.AbsoluteDate;
/** Ground measurements generator (sensor to ground mapping).
* @author Jonathan Guinet
* @author Lucie Labat-Allee
* @author Guylaine Prat
* @since 2.0
*/
public class GroundMeasurementGenerator implements Measurable {
/** Sensor to ground mapping. */
private SensorToGroundMapping groundMapping;
/** Observables which contains ground mapping. */
private Observables observables;
/** Rugged instance */
private Rugged rugged;
/** Line sensor */
private LineSensor sensor;
/** Number of lines of the sensor */
private int dimension;
/** Number of measurements */
private int measurementCount;
/** Simple constructor.
* @param rugged Rugged instance
* @param sensorName sensor name
* @param dimension number of line of the sensor
*/
public GroundMeasurementGenerator(final Rugged rugged, final String sensorName, final int dimension) {
// Generate reference mapping
this.groundMapping = new SensorToGroundMapping(rugged.getName(), sensorName);
// Create observables for one model
this.observables = new Observables(1);
this.rugged = rugged;
this.sensor = rugged.getLineSensor(groundMapping.getSensorName());
this.dimension = dimension;
this.measurementCount = 0;
}
/** Get the sensor to ground mapping
* @return the sensor to ground mapping
*/
public SensorToGroundMapping getGroundMapping() {
return groundMapping;
}
/** Get the observables which contains the ground mapping
* @return the observables which contains the ground mapping
*/
public Observables getObservables() {
return observables;
}
@Override
public int getMeasurementCount() {
return measurementCount;
}
@Override
public void createMeasurement(final int lineSampling, final int pixelSampling) {
for (double line = 0; line < dimension; line += lineSampling) {
final AbsoluteDate date = sensor.getDate(line);
for (int pixel = 0; pixel < sensor.getNbPixels(); pixel += pixelSampling) {
final GeodeticPoint gp2 = rugged.directLocation(date, sensor.getPosition(),
sensor.getLOS(date, pixel));
groundMapping.addMapping(new SensorPixel(line, pixel), gp2);
// increment the number of measurements
measurementCount++;
}
}
observables.addGroundMapping(groundMapping);
}
@Override
public void createNoisyMeasurement(final int lineSampling, final int pixelSampling, final Noise noise) {
// Estimate latitude and longitude errors (rad)
final Vector3D latLongError = estimateLatLongError();
// Get noise features
final double[] mean = noise.getMean(); // [latitude, longitude, altitude] mean
final double[] std = noise.getStandardDeviation(); // [latitude, longitude, altitude] standard deviation
final double latErrorMean = mean[0] * latLongError.getX();
final double lonErrorMean = mean[1] * latLongError.getY();
final double latErrorStd = std[0] * latLongError.getX();
final double lonErrorStd = std[1] * latLongError.getY();
// Gaussian random generator
// Build a null mean random uncorrelated vector generator with standard deviation corresponding to the estimated error on ground
final double meanGenerator[] = {latErrorMean, lonErrorMean, mean[2]};
final double stdGenerator[] = {latErrorStd, lonErrorStd, std[2]};
// seed has been fixed for tests purpose
final GaussianRandomGenerator rng = new GaussianRandomGenerator(new Well19937a(0xefac03d9be4d24b9l));
final UncorrelatedRandomVectorGenerator rvg = new UncorrelatedRandomVectorGenerator(meanGenerator, stdGenerator, rng);
for (double line = 0; line < dimension; line += lineSampling) {
final AbsoluteDate date = sensor.getDate(line);
for (int pixel = 0; pixel < sensor.getNbPixels(); pixel += pixelSampling) {
// Components of generated vector follow (independent) Gaussian distribution
final Vector3D vecRandom = new Vector3D(rvg.nextVector());
final GeodeticPoint gp2 = rugged.directLocation(date, sensor.getPosition(),
sensor.getLOS(date, pixel));
final GeodeticPoint gpNoisy = new GeodeticPoint(gp2.getLatitude() + vecRandom.getX(),
gp2.getLongitude() + vecRandom.getY(),
gp2.getAltitude() + vecRandom.getZ());
groundMapping.addMapping(new SensorPixel(line, pixel), gpNoisy);
// increment the number of measurements
measurementCount++;
}
}
this.observables.addGroundMapping(groundMapping);
}
/** Compute latitude and longitude errors
* @return the latitude and longitude errors (rad)
*/
private Vector3D estimateLatLongError() {
final int pix = sensor.getNbPixels() / 2;
final int line = (int) FastMath.floor(pix); // assumption : same number of line and pixels;
final AbsoluteDate date = sensor.getDate(line);
final GeodeticPoint gp_pix0 = rugged.directLocation(date, sensor.getPosition(), sensor.getLOS(date, pix));
final AbsoluteDate date1 = sensor.getDate(line + 1);
final GeodeticPoint gp_pix1 = rugged.directLocation(date1, sensor.getPosition(), sensor.getLOS(date1, pix + 1));
final double latErr = FastMath.abs(gp_pix0.getLatitude() - gp_pix1.getLatitude());
final double lonErr = FastMath.abs(gp_pix0.getLongitude() - gp_pix1.getLongitude());
return new Vector3D(latErr, lonErr, 0.0);
}
}
src/tutorials/java/fr/cs/examples/refiningPleiades/generators/InterMeasurementGenerator.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You 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.
*/
package fr.cs.examples.refiningPleiades.generators;
import org.hipparchus.random.GaussianRandomGenerator;
import org.hipparchus.random.UncorrelatedRandomVectorGenerator;
import org.hipparchus.random.Well19937a;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.errors.RuggedMessages;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.SensorPixel;
import org.orekit.rugged.adjustment.measurements.Observables;
import org.orekit.rugged.adjustment.measurements.SensorToSensorMapping;
import org.orekit.rugged.utils.SpacecraftToObservedBody;
import org.orekit.time.AbsoluteDate;
import fr.cs.examples.refiningPleiades.metrics.DistanceTools;
/**
* Inter-measurements generator (sensor to sensor mapping).
* @author Jonathan Guinet
* @author Lucie Labatallee
* @author Guylaine Prat
* @since 2.0
*/
public class InterMeasurementGenerator implements Measurable {
/** Mapping from sensor A to sensor B. */
private SensorToSensorMapping interMapping;
/** Observables which contains sensor to sensor mapping.*/
private Observables observables;
/** Rugged instance corresponding to the viewing model A */
private Rugged ruggedA;
/** Rugged instance corresponding to the viewing model B */
private Rugged ruggedB;
/** Sensor A */
private LineSensor sensorA;
/** Sensor B */
private LineSensor sensorB;
/** Number of measurements */
private int measurementCount;
/** Sensor name B */
private String sensorNameB;
/** Number of line for acquisition A */
private int dimensionA;
/** Number of line for acquisition B */
private int dimensionB;
/** Limit value for outlier points */
private double outlier;
/** Default constructor: measurements generation without outlier points control
* and without Earth distance constraint.
* @param ruggedA Rugged instance corresponding to the viewing model A
* @param sensorNameA sensor name A
* @param dimensionA number of line for acquisition A
* @param ruggedB Rugged instance corresponding to the viewing model B
* @param sensorNameB sensor name B
* @param dimensionB number of line for acquisition B
*/
public InterMeasurementGenerator(final Rugged ruggedA, final String sensorNameA, final int dimensionA,
final Rugged ruggedB, final String sensorNameB, final int dimensionB) {
// Initialize parameters
initParams(ruggedA, sensorNameA, dimensionA, ruggedB, sensorNameB, dimensionB);
// Generate reference mapping, without Earth distance constraint
interMapping = new SensorToSensorMapping(sensorNameA, sensorNameB);
// Create observables for two models
observables = new Observables(2);
}
/** Constructor for measurements generation taking into account outlier points control,
* without Earth distance constraint.
* @param ruggedA Rugged instance corresponding to the viewing model A
* @param sensorNameA sensor name A
* @param dimensionA dimension for acquisition A
* @param ruggedB Rugged instance corresponding to the viewing model B
* @param sensorNameB sensor name B
* @param dimensionB dimension for acquisition B
* @param outlier limit value for outlier points
*/
public InterMeasurementGenerator(final Rugged ruggedA, final String sensorNameA, final int dimensionA,
final Rugged ruggedB, final String sensorNameB, final int dimensionB,
final double outlier) {
this(ruggedA, sensorNameA, dimensionA, ruggedB, sensorNameB, dimensionB);
this.outlier = outlier;
}
/** Constructor for measurements generation taking into account outlier points control,
* and Earth distance constraint.
* @param ruggedA Rugged instance corresponding to the viewing model A
* @param sensorNameA sensor name A
* @param dimensionA dimension for acquisition A
* @param ruggedB Rugged instance corresponding to the viewing model B
* @param sensorNameB sensor name B
* @param dimensionB dimension for acquisition B
* @param outlier limit value for outlier points
* @param earthConstraintWeight weight given to the Earth distance constraint
* with respect to the LOS distance (between 0 and 1).
*/
public InterMeasurementGenerator(final Rugged ruggedA, final String sensorNameA, final int dimensionA,
final Rugged ruggedB, final String sensorNameB, final int dimensionB,
final double outlier, final double earthConstraintWeight) {
// Initialize parameters
initParams(ruggedA, sensorNameA, dimensionA, ruggedB, sensorNameB, dimensionB);
// Generate reference mapping, with Earth distance constraints.
// Weighting will be applied as follow :
// (1-earthConstraintWeight) for losDistance weighting
// earthConstraintWeight for earthDistance weighting
interMapping = new SensorToSensorMapping(sensorNameA, ruggedA.getName(), sensorNameB, ruggedB.getName(), earthConstraintWeight);
// Outlier points control
this.outlier = outlier;
// Observables which contains sensor to sensor mapping
this.observables = new Observables(2);
}
/** Get the mapping from sensor A to sensor B
* @return the mapping from sensor A to sensor B
*/
public SensorToSensorMapping getInterMapping() {
return interMapping;
}
/** Get the observables which contains sensor to sensor mapping
* @return the observables which contains sensor to sensor mapping
*/
public Observables getObservables() {
return observables;
}
@Override
public int getMeasurementCount() {
return measurementCount;
}
@Override
public void createMeasurement(final int lineSampling, final int pixelSampling) {
// Search the sensor pixel seeing point
final int minLine = 0;
final int maxLine = dimensionB - 1;
for (double line = 0; line < dimensionA; line += lineSampling) {
final AbsoluteDate dateA = sensorA.getDate(line);
for (double pixelA = 0; pixelA < sensorA.getNbPixels(); pixelA += pixelSampling) {
final GeodeticPoint gpA = ruggedA.directLocation(dateA, sensorA.getPosition(),
sensorA.getLOS(dateA, pixelA));
final SensorPixel sensorPixelB = ruggedB.inverseLocation(sensorNameB, gpA, minLine, maxLine);
// We need to test if the sensor pixel is found in the prescribed lines
// otherwise the sensor pixel is null
if (sensorPixelB != null) {
final AbsoluteDate dateB = sensorB.getDate(sensorPixelB.getLineNumber());
final double pixelB = sensorPixelB.getPixelNumber();
final SpacecraftToObservedBody scToBodyA = ruggedA.getScToBody();
final GeodeticPoint gpB = ruggedB.directLocation(dateB, sensorB.getPosition(),
sensorB.getLOS(dateB, pixelB));
final double geoDistance = DistanceTools.computeDistanceInMeter(gpA, gpB);
if (geoDistance < this.outlier) {
final SensorPixel realPixelA = new SensorPixel(line, pixelA);
final SensorPixel realPixelB = new SensorPixel(sensorPixelB.getLineNumber(), sensorPixelB.getPixelNumber());
final AbsoluteDate realDateA = sensorA.getDate(realPixelA.getLineNumber());
final AbsoluteDate realDateB = sensorB.getDate(realPixelB.getLineNumber());
final double[] distanceLOSB = ruggedB.distanceBetweenLOS(sensorA, realDateA, realPixelA.getPixelNumber(), scToBodyA,
sensorB, realDateB, realPixelB.getPixelNumber());
final double losDistance = 0.0;
final double earthDistance = distanceLOSB[1];
interMapping.addMapping(realPixelA, realPixelB, losDistance, earthDistance);
// Increment the number of measurements
this.measurementCount++;
}
}
}
}
this.observables.addInterMapping(interMapping);
}
/** Tie point creation from direct 1ocation with Rugged A and inverse location with Rugged B
* @param lineSampling sampling along lines
* @param pixelSampling sampling along columns
* @param noise errors to apply to measure for pixel A and pixel B
*/
public void createNoisyMeasurement(final int lineSampling, final int pixelSampling, final Noise noise) {
// Get noise features (errors)
// [pixelA, pixelB] mean
final double[] mean = noise.getMean();
// [pixelA, pixelB] standard deviation
final double[] std = noise.getStandardDeviation();
// Search the sensor pixel seeing point
final int minLine = 0;
final int maxLine = dimensionB - 1;
final double meanA[] = { mean[0], mean[0] };
final double stdA[] = { std[0], std[0] };
final double meanB[] = { mean[1], mean[1] };
final double stdB[] = { std[1], std[1] };
// Seed has been fixed for tests purpose
final GaussianRandomGenerator rngA = new GaussianRandomGenerator(new Well19937a(0xefac03d9be4d24b9l));
final UncorrelatedRandomVectorGenerator rvgA = new UncorrelatedRandomVectorGenerator(meanA, stdA, rngA);
// Seed has been fixed for tests purpose
final GaussianRandomGenerator rngB = new GaussianRandomGenerator(new Well19937a(0xdf1c03d9be0b34b9l));
final UncorrelatedRandomVectorGenerator rvgB = new UncorrelatedRandomVectorGenerator(meanB, stdB, rngB);
for (double line = 0; line < dimensionA; line += lineSampling) {
final AbsoluteDate dateA = sensorA.getDate(line);
for (double pixelA = 0; pixelA < sensorA.getNbPixels(); pixelA += pixelSampling) {
final GeodeticPoint gpA = ruggedA.directLocation(dateA, sensorA.getPosition(),
sensorA.getLOS(dateA, pixelA));
final SensorPixel sensorPixelB = ruggedB.inverseLocation(sensorNameB, gpA, minLine, maxLine);
// We need to test if the sensor pixel is found in the prescribed lines
// otherwise the sensor pixel is null
if (sensorPixelB != null) {
final AbsoluteDate dateB = sensorB.getDate(sensorPixelB.getLineNumber());
final double pixelB = sensorPixelB.getPixelNumber();
// Get spacecraft to body transform of Rugged instance A
final SpacecraftToObservedBody scToBodyA = ruggedA.getScToBody();
final GeodeticPoint gpB = ruggedB.directLocation(dateB, sensorB.getPosition(),
sensorB.getLOS(dateB, pixelB));
final double geoDistance = DistanceTools.computeDistanceInMeter(gpA, gpB);
// Create the inter mapping if distance is below outlier value
if (geoDistance < outlier) {
final double[] vecRandomA = rvgA.nextVector();
final double[] vecRandomB = rvgB.nextVector();
final SensorPixel realPixelA = new SensorPixel(line + vecRandomA[0], pixelA + vecRandomA[1]);
final SensorPixel realPixelB = new SensorPixel(sensorPixelB.getLineNumber() + vecRandomB[0],
sensorPixelB.getPixelNumber() + vecRandomB[1]);
final AbsoluteDate realDateA = sensorA.getDate(realPixelA.getLineNumber());
final AbsoluteDate realDateB = sensorB.getDate(realPixelB.getLineNumber());
final double[] distanceLOSB = ruggedB.distanceBetweenLOS(sensorA, realDateA, realPixelA.getPixelNumber(), scToBodyA,
sensorB, realDateB, realPixelB.getPixelNumber());
final Double losDistance = 0.0;
final Double earthDistance = distanceLOSB[1];
interMapping.addMapping(realPixelA, realPixelB, losDistance, earthDistance);
// Increment the number of measurements
this.measurementCount++;
} // end test if geoDistance < outlier
} // end test if sensorPixelB != null
} // end loop on pixel of sensorA
} // end loop on line of sensorA
this.observables.addInterMapping(interMapping);
}
/** Default constructor: measurements generation without outlier points control
* and Earth distances constraint.
* @param rA Rugged instance A
* @param sNameA sensor name A
* @param dimA dimension for acquisition A
* @param rB Rugged instance B
* @param sNameB sensor name B
* @param dimB dimension for acquisition B
*/
private void initParams(final Rugged rA, final String sNameA, final int dimA,
final Rugged rB, final String sNameB, final int dimB) {
this.sensorNameB = sNameB;
// Check that sensors's name is different
if (sNameA.contains(sNameB)) {
throw new RuggedException(RuggedMessages.DUPLICATED_PARAMETER_NAME, sNameA);
}
this.ruggedA = rA;
this.ruggedB = rB;
this.sensorA = rA.getLineSensor(sNameA);
this.sensorB = rB.getLineSensor(sNameB);
this.dimensionA = dimA;
this.dimensionB = dimB;
this.measurementCount = 0;
}
}
src/tutorials/java/fr/cs/examples/refiningPleiades/generators/Measurable.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You 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.
*/
package fr.cs.examples.refiningPleiades.generators;
/** For measurements generator.
* @author Lucie Labat-Allee
* @author Guylaine Prat
* @since 2.0
*/
public interface Measurable {
/** Get the number of measurements
* @return the number of measurements
*/
int getMeasurementCount();
/** Create measurements (without noise)
* @param lineSampling line sampling
* @param pixelSampling pixel sampling
*/
void createMeasurement(int lineSampling, int pixelSampling);
/** Create noisy measurements
* @param lineSampling line sampling
* @param pixelSampling pixel sampling
* @param noise the noise to add to the measurements
*/
void createNoisyMeasurement(int lineSampling, int pixelSampling, Noise noise);
}
src/tutorials/java/fr/cs/examples/refiningPleiades/generators/Noise.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You 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.
*/
package fr.cs.examples.refiningPleiades.generators;
/** Class for adding a noise to measurements.
* @author Lucie Labat-Allee
* @author Guylaine Prat
* @since 2.0
*/
public class Noise {
/** Type of distribution. */
private static final int GAUSSIAN = 0;
/** Dimension. */
private int dimension;
/** Mean. */
private double[] mean;
/** Standard deviation. */
private double[] standardDeviation;
/** Distribution. */
private int distribution = GAUSSIAN;
/** Build a new instance.
* @param distribution noise type
* @param dimension noise dimension
*/
public Noise(final int distribution, final int dimension) {
this.mean = new double[dimension];
this.standardDeviation = new double[dimension];
this.dimension = dimension;
this.distribution = distribution;
}
/** Get the mean.
* @return the mean
*/
public double[] getMean() {
return mean.clone();
}
/** Set the mean.
* @param meanValue the mean to set
*/
public void setMean(final double[] meanValue) {
this.mean = meanValue.clone();
}
/** Get the standard deviation.
* @return the standard deviation
*/
public double[] getStandardDeviation() {
return standardDeviation.clone();
}
/** Set the standard deviation.
* @param standardDeviationValue the standard deviation to set
*/
public void setStandardDeviation(final double[] standardDeviationValue) {
this.standardDeviation = standardDeviationValue.clone();
}
/** Get the distribution.
* @return the distribution
*/
public int getDistribution() {
return distribution;
}
/** Get the dimension.
* @return the dimension
*/
public int getDimension() {
return dimension;
}
}
src/tutorials/java/fr/cs/examples/refiningPleiades/metrics/DistanceTools.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You 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.
*/
package fr.cs.examples.refiningPleiades.metrics;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.utils.Constants;
/**
* Class for computing geodetic and anguler distance between two geodetic points.
* @author Jonathan Guinet
* @author Guylaine Prat
* @since 2.0
*/
public class DistanceTools {
/** Private constructor for utility class.
*/
private DistanceTools() {
}
/** Choice the method for computing geodetic distance between two points.
* @param geoPoint1 first geodetic point (rad)
* @param geoPoint2 second geodetic point (rad)
* @param computeAngular if true, distance will be angular, otherwise will be in meters
* @return distance in meters or radians if flag computeAngular is true
*/
public static double computeDistance(final GeodeticPoint geoPoint1, final GeodeticPoint geoPoint2,
final boolean computeAngular) {
if (computeAngular == true) {
return computeDistanceAngular(geoPoint1, geoPoint2);
} else {
return computeDistanceInMeter(geoPoint1, geoPoint2);
}
}
/** Compute a geodetic distance in meters between two geodetic points..
* @param geoPoint1 first geodetic point (rad)
* @param geoPoint2 second geodetic point (rad)
* @return distance in meters
*/
public static double computeDistanceInMeter(final GeodeticPoint geoPoint1, final GeodeticPoint geoPoint2) {
// get vectors on unit sphere from angular coordinates
final Vector3D p1 = new Vector3D(geoPoint1.getLatitude(), geoPoint1.getLongitude()); //
final Vector3D p2 = new Vector3D(geoPoint2.getLatitude(), geoPoint2.getLongitude());
final double distance = Constants.WGS84_EARTH_EQUATORIAL_RADIUS * Vector3D.angle(p1, p2);
return distance;
}
/** Compute an angular distance between two geodetic points.
* @param geoPoint1 first geodetic point (rad)
* @param geoPoint2 second geodetic point (rad)
* @return angular distance in radians
*/
public static double computeDistanceAngular(final GeodeticPoint geoPoint1, final GeodeticPoint geoPoint2) {
final double lonDiff = geoPoint1.getLongitude() - geoPoint2.getLongitude();
final double latDiff = geoPoint1.getLatitude() - geoPoint2.getLatitude();
return FastMath.hypot(lonDiff, latDiff);
}
}
src/tutorials/java/fr/cs/examples/refiningPleiades/metrics/LocalisationMetrics.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You 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.
*/
package fr.cs.examples.refiningPleiades.metrics;
import java.util.Iterator;
import java.util.Map;
import java.util.Set;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.adjustment.measurements.SensorToGroundMapping;
import org.orekit.rugged.adjustment.measurements.SensorToSensorMapping;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.SensorPixel;
import org.orekit.rugged.utils.SpacecraftToObservedBody;
import org.orekit.time.AbsoluteDate;
/**
* Class for testing geometric performances in absolute location.
* Metrics are computed for two scenarios: ground control points and tie points.
* @see SensorToSensorMapping
* @see SensorToGroundMapping
* @author Jonathan Guinet
* @author Lucie Labat-Allee
* @author Guylaine Prat
* @since 2.0
*/
public class LocalisationMetrics {
/** Maximum residual distance. */
private double resMax;
/** Mean residual distance. */
private double resMean;
/** LOS distance max. */
private double losDistanceMax;
/** LOS distance mean. */
private double losDistanceMean;
/** Earth distance max. */
private double earthDistanceMax;
/** Earth distance mean.*/
private double earthDistanceMean;
/** Compute metrics corresponding to the Ground Control Points (GCP) study.
* @param measMapping Mapping of observations/measurements = the ground truth
* @param rugged Rugged instance
* @param computeAngular flag to know if distance is computed in meters (false) or with angular (true)
*/
public LocalisationMetrics(final SensorToGroundMapping measMapping, final Rugged rugged, final boolean computeAngular) {
// Initialization
this.resMax = 0.0;
this.resMean = 0.0;
// Compute metrics - Case of Sensor to Ground mapping (Ground Control Points GCP study)
computeGCPmetrics(measMapping, rugged, computeAngular);
}
/** Compute metrics corresponding to the tie points study.
* @param measMapping Mapping of observations/measurements = the ground truth
* @param ruggedA Rugged instance corresponding to viewing model A
* @param ruggedB Rugged instance corresponding to viewing model B
* @param computeAngular flag to know if distance is computed in meters (false) or with angular (true)
*/
public LocalisationMetrics(final SensorToSensorMapping measMapping, final Rugged ruggedA, final Rugged ruggedB,
final boolean computeAngular) {
// Initialization
this.resMax = 0.0;
this.resMean = 0.0;
this.losDistanceMax = 0.0;
this.losDistanceMean = 0.0;
this.earthDistanceMax = 0.0;
this.earthDistanceMean = 0.0;
// Compute metrics - Case of Sensor to Sensor mapping (tie points study)
computeTiePointsMetrics(measMapping, ruggedA, ruggedB, computeAngular);
}
/**
* Compute metrics: case of ground control points (GCP).
* @param measMapping Mapping of observations/measurements = the ground truth
* @param rugged Rugged instance
* @param computeAngular flag to know if distance is computed in meters (false) or with angular (true)
*/
public void computeGCPmetrics(final SensorToGroundMapping measMapping, final Rugged rugged, final boolean computeAngular) {
// Mapping of observations/measurements = the ground truth
final Set<Map.Entry<SensorPixel, GeodeticPoint>> measurementsMapping;
final LineSensor lineSensor;
// counter for compute mean distance
double count = 0;
// Initialization
measurementsMapping = measMapping.getMapping();
lineSensor = rugged.getLineSensor(measMapping.getSensorName());
// number of measurements
int nbMeas = measurementsMapping.size();
final Iterator<Map.Entry<SensorPixel, GeodeticPoint>> gtIt = measurementsMapping.iterator();
// Browse map of measurements
while (gtIt.hasNext()) {
final Map.Entry<SensorPixel, GeodeticPoint> gtMeas = gtIt.next();
final SensorPixel gtSP = gtMeas.getKey();
final GeodeticPoint gtGP = gtMeas.getValue();
final AbsoluteDate date = lineSensor.getDate(gtSP.getLineNumber());
final GeodeticPoint esGP = rugged.directLocation(date, lineSensor.getPosition(),
lineSensor.getLOS(date,
(int) gtSP.getPixelNumber()));
// Compute distance
double distance = DistanceTools.computeDistance(esGP, gtGP, computeAngular);
count += distance;
if (distance > resMax) {
resMax = distance;
}
}
// Mean of residuals
resMean = count / nbMeas;
}
/**
* Compute metrics: case of tie points.
* @param measMapping Mapping of observations/measurements = the ground truth
* @param ruggedA Rugged instance corresponding to viewing model A
* @param ruggedB Rugged instance corresponding to viewing model B
* @param computeAngular Flag to know if distance is computed in meters (false) or with angular (true)
*/
public void computeTiePointsMetrics(final SensorToSensorMapping measMapping, final Rugged ruggedA, final Rugged ruggedB,
final boolean computeAngular) {
// Mapping of observations/measurements = the ground truth
final Set<Map.Entry<SensorPixel, SensorPixel>> measurementsMapping;
final LineSensor lineSensorA; // line sensor corresponding to rugged A
final LineSensor lineSensorB; // line sensor corresponding to rugged B
double count = 0; // counter for computing remaining mean distance
double losDistanceCount = 0; // counter for computing LOS distance mean
double earthDistanceCount = 0; // counter for computing Earth distance mean
int i = 0; // increment of measurements
// Initialization
measurementsMapping = measMapping.getMapping();
lineSensorA = ruggedA.getLineSensor(measMapping.getSensorNameA());
lineSensorB = ruggedB.getLineSensor(measMapping.getSensorNameB());
int nbMeas = measurementsMapping.size(); // number of measurements
// Browse map of measurements
for (Iterator<Map.Entry<SensorPixel, SensorPixel>> gtIt = measurementsMapping.iterator();
gtIt.hasNext();
i++) {
if (i == measurementsMapping.size()) {
break;
}
final Map.Entry<SensorPixel, SensorPixel> gtMapping = gtIt.next();
final SensorPixel spA = gtMapping.getKey();
final SensorPixel spB = gtMapping.getValue();
final AbsoluteDate dateA = lineSensorA.getDate(spA.getLineNumber());
final AbsoluteDate dateB = lineSensorB.getDate(spB.getLineNumber());
final double pixelA = spA.getPixelNumber();
final double pixelB = spB.getPixelNumber();
final Vector3D losA = lineSensorA.getLOS(dateA, pixelA);
final Vector3D losB = lineSensorB.getLOS(dateB, pixelB);
final GeodeticPoint gpA = ruggedA.directLocation(dateA, lineSensorA.getPosition(), losA);
final GeodeticPoint gpB = ruggedB.directLocation(dateB, lineSensorB.getPosition(), losB);
final SpacecraftToObservedBody scToBodyA = ruggedA.getScToBody();
// Estimated distances (LOS / Earth)
final double[] distances = ruggedB.distanceBetweenLOS(lineSensorA, dateA, pixelA, scToBodyA, lineSensorB, dateB, pixelB);
// LOS distance control
final double estLosDistance = distances[0]; // LOS distance estimation
if (estLosDistance > losDistanceMax) {
losDistanceMax = estLosDistance;
}
losDistanceCount += estLosDistance;
// Earth distance control
final double estEarthDistance = distances[1]; // Earth distance estimation
final double measEarthDistance = measMapping.getBodyDistance(i).doubleValue(); // Earth measurement distance
final double earthDistance = FastMath.abs(estEarthDistance - measEarthDistance);
if (earthDistance > earthDistanceMax) {
earthDistanceMax = earthDistance;
}
earthDistanceCount += earthDistance;
// Compute remaining distance
double distance = DistanceTools.computeDistance(gpB, gpA, computeAngular);
count += distance;
if (distance > resMax) {
resMax = distance;
}
}
resMean = count / nbMeas;
losDistanceMean = losDistanceCount / nbMeas;
earthDistanceMean = earthDistanceCount / nbMeas;
}
/** Get the max residual.
* @return maximum of residuals
*/
public double getMaxResidual() {
return resMax;
}
/** Get the mean of residuals.
* @return mean of residuals
*/
public double getMeanResidual() {
return resMean;
}
/** Get the LOS maximum distance.
* @return LOS maximum distance
*/
public double getLosMaxDistance() {
return losDistanceMax;
}
/** Get the LOS mean distance.
* @return LOS mean distance
*/
public double getLosMeanDistance() {
return losDistanceMean;
}
/** Get the Earth distance maximum residual.
* @return Earth distance max residual
*/
public double getEarthMaxDistance() {
return earthDistanceMax;
}
/** Get the Earth distance mean residual.
* @return Earth distance mean residual
*/
public double getEarthMeanDistance() {
return earthDistanceMean;
}
}
src/tutorials/java/fr/cs/examples/refiningPleiades/models/OrbitModel.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You 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.
*/
package fr.cs.examples.refiningPleiades.models;
import java.util.ArrayList;
import java.util.List;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.RotationConvention;
import org.hipparchus.geometry.euclidean.threed.RotationOrder;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.attitudes.LofOffset;
import org.orekit.attitudes.NadirPointing;
import org.orekit.attitudes.TabulatedLofOffset;
import org.orekit.attitudes.YawCompensation;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.LOFType;
import org.orekit.frames.Transform;
import org.orekit.orbits.CircularOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.Propagator;
import org.orekit.propagation.analytical.KeplerianPropagator;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.utils.AngularDerivativesFilter;
import org.orekit.utils.Constants;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.TimeStampedAngularCoordinates;
import org.orekit.utils.TimeStampedPVCoordinates;
/**
* Orbit Model class to generate positions-velocities and attitude quaternions.
* <p>
* the aim of this class is to simulate the orbit model of an LEO satellite
* </p>
* @author Jonathan Guinet
* @author Guylaine Prat
* @since 2.0
*/
public class OrbitModel {
/** Flag to change Attitude Law (by default Nadir Pointing Yaw Compensation). */
private boolean userDefinedLOFTransform;
/** User defined Roll law. */
private double[] lofTransformRollPoly;
/** User defined Pitch law. */
private double[] lofTransformPitchPoly;
/** User defined Yaw law. */
private double[] lofTransformYawPoly;
/** Reference date. */
private AbsoluteDate refDate;
/** Default constructor.
*/
public OrbitModel() {
userDefinedLOFTransform = false;
}
/** Generate satellite ephemeris.
*/
public static void addSatellitePV(final TimeScale gps, final Frame eme2000, final Frame itrf,
final List<TimeStampedPVCoordinates> satellitePVList,
final String absDate,
final double px, final double py, final double pz,
final double vx, final double vy, final double vz) {
final AbsoluteDate ephemerisDate = new AbsoluteDate(absDate, gps);
final Vector3D position = new Vector3D(px, py, pz);
final Vector3D velocity = new Vector3D(vx, vy, vz);
final PVCoordinates pvITRF = new PVCoordinates(position, velocity);
final Transform transform = itrf.getTransformTo(eme2000, ephemerisDate);
final Vector3D pEME2000 = transform.transformPosition(pvITRF.getPosition());
final Vector3D vEME2000 = transform.transformVector(pvITRF.getVelocity());
satellitePVList.add(new TimeStampedPVCoordinates(ephemerisDate, pEME2000, vEME2000, Vector3D.ZERO));
}
/** Generate satellite attitude.
*/
public void addSatelliteQ(final TimeScale gps,
final List<TimeStampedAngularCoordinates> satelliteQList,
final String absDate,
final double q0, final double q1, final double q2, final double q3) {
final AbsoluteDate attitudeDate = new AbsoluteDate(absDate, gps);
final Rotation rotation = new Rotation(q0, q1, q2, q3, true);
final TimeStampedAngularCoordinates pair =
new TimeStampedAngularCoordinates(attitudeDate, rotation, Vector3D.ZERO, Vector3D.ZERO);
satelliteQList.add(pair);
}
/** Create an Earth.
* @return the Earth as the WGS84 ellipsoid
*/
public BodyShape createEarth() {
return new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING,
FramesFactory.getITRF(IERSConventions.IERS_2010, true));
}
/** Created a gravity field.
* @return normalized spherical harmonics coefficients
*/
public NormalizedSphericalHarmonicsProvider createGravityField() {
return GravityFieldFactory.getNormalizedProvider(12, 12);
}
/** Create an orbit at a chosen date.
* @param mu Earth gravitational constant
* @return the orbit
*/
public Orbit createOrbit(final double mu, final AbsoluteDate date) {
// the following orbital parameters have been computed using
// Orekit tutorial about phasing, using the following configuration:
//
// orbit.date = 2012-01-01T00:00:00.000
// phasing.orbits.number = 143
// phasing.days.number = 10
// sun.synchronous.reference.latitude = 0
// sun.synchronous.reference.ascending = false
// sun.synchronous.mean.solar.time = 10:30:00
// gravity.field.degree = 12
// gravity.field.order = 12
final Frame eme2000 = FramesFactory.getEME2000();
return new CircularOrbit(694000.0 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
-4.029194321683225E-4,
0.0013530362644647786,
FastMath.toRadians(98.2), // Pleiades inclination 98.2 deg
FastMath.toRadians(-86.47 + 180),
FastMath.toRadians(135.9 + 0.3),
PositionAngle.TRUE,
eme2000,
date,
mu);
}
/** Set the Local Orbital Frame characteristics.
*/
public void setLOFTransform(final double[] rollPoly, final double[] pitchPoly,
final double[] yawPoly, final AbsoluteDate date) {
this.userDefinedLOFTransform = true;
lofTransformRollPoly = rollPoly.clone();
lofTransformPitchPoly = pitchPoly.clone();
lofTransformYawPoly = yawPoly.clone();
this.refDate = date;
}
/** Recompute the polynom coefficients with shift.
*/
private double getPoly(final double[] poly, final double shift) {
double val = 0.0;
for (int coef = 0; coef < poly.length; coef++) {
val = val + poly[coef] * FastMath.pow(shift, coef);
}
return val;
}
/** Get the offset.
*/
private Rotation getOffset(final BodyShape earth, final Orbit orbit, final double shift){
final LOFType type = LOFType.VVLH;
final double roll = getPoly(lofTransformRollPoly, shift);
final double pitch = getPoly(lofTransformPitchPoly, shift);
final double yaw = getPoly(lofTransformYawPoly, shift);
final LofOffset law = new LofOffset(orbit.getFrame(), type, RotationOrder.XYZ,
roll, pitch, yaw);
final Rotation offsetAtt = law.
getAttitude(orbit, orbit.getDate().shiftedBy(shift), orbit.getFrame()).
getRotation();
final NadirPointing nadirPointing = new NadirPointing(orbit.getFrame(), earth);
final Rotation nadirAtt = nadirPointing.
getAttitude(orbit, orbit.getDate().getDate().shiftedBy(shift), orbit.getFrame()).
getRotation();
final Rotation offsetProper = offsetAtt.compose(nadirAtt.revert(), RotationConvention.VECTOR_OPERATOR);
return offsetProper;
}
/** Create the attitude provider.
*/
public AttitudeProvider createAttitudeProvider(final BodyShape earth, final Orbit orbit) {
if (userDefinedLOFTransform) {
final LOFType type = LOFType.VVLH;
final List<TimeStampedAngularCoordinates> list = new ArrayList<TimeStampedAngularCoordinates>();
for (double shift = -10.0; shift < +10.0; shift += 1e-2) {
list.add(new TimeStampedAngularCoordinates(refDate
.shiftedBy(shift), getOffset(earth, orbit, shift),
Vector3D.ZERO,
Vector3D.ZERO));
}
final TabulatedLofOffset tabulated = new TabulatedLofOffset(orbit.getFrame(), type, list,
2, AngularDerivativesFilter.USE_R);
return tabulated;
} else {
return new YawCompensation(orbit.getFrame(), new NadirPointing(orbit.getFrame(), earth));
}
}
/** Generate the orbit.
*/
public List<TimeStampedPVCoordinates> orbitToPV(final Orbit orbit, final BodyShape earth,
final AbsoluteDate minDate, final AbsoluteDate maxDate,
final double step) {
final Propagator propagator = new KeplerianPropagator(orbit);
propagator.setAttitudeProvider(createAttitudeProvider(earth, orbit));
propagator.propagate(minDate);
final List<TimeStampedPVCoordinates> list = new ArrayList<TimeStampedPVCoordinates>();
propagator.getMultiplexer().add(step,
currentState ->
list.add(new TimeStampedPVCoordinates(currentState.getDate(),
currentState.getPVCoordinates().getPosition(),
currentState.getPVCoordinates().getVelocity(),
Vector3D.ZERO)));
propagator.propagate(maxDate);
return list;
}
/** Generate the attitude.
*/
public List<TimeStampedAngularCoordinates> orbitToQ(final Orbit orbit, final BodyShape earth,
final AbsoluteDate minDate, final AbsoluteDate maxDate,
final double step) {
final Propagator propagator = new KeplerianPropagator(orbit);
propagator.setAttitudeProvider(createAttitudeProvider(earth, orbit));
propagator.propagate(minDate);
final List<TimeStampedAngularCoordinates> list = new ArrayList<>();
propagator.getMultiplexer().add(step,
currentState ->
list.add(new TimeStampedAngularCoordinates(currentState.getDate(),
currentState.getAttitude().getRotation(),
Vector3D.ZERO,
Vector3D.ZERO)));
propagator.propagate(maxDate);
return list;
}
}
src/tutorials/java/fr/cs/examples/refiningPleiades/models/PleiadesViewingModel.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You 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.
*/
package fr.cs.examples.refiningPleiades.models;
import java.util.ArrayList;
import java.util.List;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.RotationConvention;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.rugged.linesensor.LineDatation;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.LinearLineDatation;
import org.orekit.rugged.los.FixedRotation;
import org.orekit.rugged.los.FixedZHomothety;
import org.orekit.rugged.los.LOSBuilder;
import org.orekit.rugged.los.TimeDependentLOS;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
import fr.cs.examples.refiningPleiades.Refining;
/**
* Pleiades viewing model class definition.
* The aim of this class is to simulate PHR sensor.
* @author Jonathan Guinet
* @author Lucie Labat-Allee
* @author Guylaine Prat
* @since 2.0
*/
public class PleiadesViewingModel {
/** Pleiades parameters. */
private static final double FOV = 1.65; // 20km - alt 694km
private static final int DIMENSION = 40000;
private static final double LINE_PERIOD = 1.e-4;
private double angle;
private LineSensor lineSensor;
private String date;
private String sensorName;
/** Simple constructor.
* <p>
* initialize PleiadesViewingModel with
* sensorName="line", incidenceAngle = 0.0, date = "2016-01-01T12:00:00.0"
* </p>
*/
public PleiadesViewingModel(final String sensorName) {
this(sensorName, 0.0, "2016-01-01T12:00:00.0");
}
/** PleiadesViewingModel constructor.
* @param sensorName sensor name
* @param incidenceAngle incidence angle
* @param referenceDate reference date
*/
public PleiadesViewingModel(final String sensorName, final double incidenceAngle, final String referenceDate) {
this.sensorName = sensorName;
this.date = referenceDate;
this.angle = incidenceAngle;
this.createLineSensor();
}
/** Create raw fixed Line Of sight
*/
public LOSBuilder rawLOS(final Vector3D center, final Vector3D normal, final double halfAperture, final int n) {
final List<Vector3D> list = new ArrayList<Vector3D>(n);
for (int i = 0; i < n; ++i) {
final double alpha = (halfAperture * (2 * i + 1 - n)) / (n - 1);
list.add(new Rotation(normal, alpha, RotationConvention.VECTOR_OPERATOR).applyTo(center));
}
return new LOSBuilder(list);
}
/** Build a LOS provider
*/
public TimeDependentLOS buildLOS() {
final LOSBuilder losBuilder = rawLOS(new Rotation(Vector3D.PLUS_I,
FastMath.toRadians(this.angle),
RotationConvention.VECTOR_OPERATOR).applyTo(Vector3D.PLUS_K),
Vector3D.PLUS_I, FastMath.toRadians(FOV / 2), DIMENSION);
losBuilder.addTransform(new FixedRotation(sensorName + Refining.getRollsuffix(), Vector3D.MINUS_I, 0.00));
losBuilder.addTransform(new FixedRotation(sensorName + Refining.getPitchsuffix(), Vector3D.MINUS_J, 0.00));
// factor is a common parameters shared between all Pleiades models
losBuilder.addTransform(new FixedZHomothety(Refining.getFactorname(), 1.0));
return losBuilder.build();
}
/** Get the reference date.
*/
public AbsoluteDate getDatationReference() {
// We use Orekit for handling time and dates, and Rugged for defining the datation model:
final TimeScale utc = TimeScalesFactory.getUTC();
return new AbsoluteDate(date, utc);
}
/** Get the min date.
*/
public AbsoluteDate getMinDate() {
return lineSensor.getDate(0);
}
/** Get the max date.
*/
public AbsoluteDate getMaxDate() {
return lineSensor.getDate(DIMENSION);
}
/** Get the line sensor.
*/
public LineSensor getLineSensor() {
return lineSensor;
}
/** Get the sensor name.
*/
public String getSensorName() {
return sensorName;
}
/** Get the number of lines of the sensor
* @return the number of lines of the sensor
*/
public int getDimension() {
return DIMENSION;
}
/** Create the line sensor
*/
private void createLineSensor() {
// Offset of the MSI from center of mass of satellite
// one line sensor
// los: swath in the (YZ) plane, looking at 50° roll, 2.6" per pixel
final Vector3D msiOffset = new Vector3D(0, 0, 0);
final TimeDependentLOS lineOfSight = buildLOS();
final double rate = 1. / LINE_PERIOD;
// linear datation model: at reference time we get the middle line, and the rate is one line every 1.5ms
final LineDatation lineDatation = new LinearLineDatation(getDatationReference(), DIMENSION / 2, rate);
lineSensor = new LineSensor(sensorName, lineDatation, msiOffset, lineOfSight);
}
}