sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] branch geoapi-4.0 updated: Complete (for now) GeodeticCalculator with spherical formulas. Ellipsoidal formulas will need to be considered later.
Date Sat, 11 May 2019 16:07:55 GMT
This is an automated email from the ASF dual-hosted git repository.

desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git


The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
     new a77e130  Complete (for now) GeodeticCalculator with spherical formulas. Ellipsoidal
formulas will need to be considered later.
a77e130 is described below

commit a77e130048538868d1e530e19523b34d2e0a585a
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Sat May 11 17:18:28 2019 +0200

    Complete (for now) GeodeticCalculator with spherical formulas. Ellipsoidal formulas will
need to be considered later.
---
 .../org/apache/sis/geometry/CoordinateFormat.java  |   4 +-
 .../internal/referencing/PositionTransformer.java  |   2 +-
 .../apache/sis/internal/referencing/Resources.java |  10 +
 .../sis/internal/referencing/Resources.properties  |   2 +
 .../internal/referencing/Resources_fr.properties   |   2 +
 .../apache/sis/referencing/GeodeticCalculator.java | 413 +++++++++++++++++----
 .../sis/referencing/GeodeticCalculatorTest.java    |  83 +++++
 .../sis/test/suite/ReferencingTestSuite.java       |   3 +-
 .../org/apache/sis/util/resources/Vocabulary.java  |  20 +
 .../sis/util/resources/Vocabulary.properties       |   4 +
 .../sis/util/resources/Vocabulary_fr.properties    |   4 +
 11 files changed, 476 insertions(+), 71 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/geometry/CoordinateFormat.java
b/core/sis-referencing/src/main/java/org/apache/sis/geometry/CoordinateFormat.java
index 7844846..7ae3bdb 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/geometry/CoordinateFormat.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/geometry/CoordinateFormat.java
@@ -62,8 +62,8 @@ import org.apache.sis.io.CompoundFormat;
  * using the following rules:
  *
  * <ul>
- *   <li>Ordinate values in angular units are formated as angles using {@link AngleFormat}.</li>
- *   <li>Ordinate values in temporal units are formated as dates using {@link DateFormat}.</li>
+ *   <li>Ordinate values in angular units are formatted as angles using {@link AngleFormat}.</li>
+ *   <li>Ordinate values in temporal units are formatted as dates using {@link DateFormat}.</li>
  *   <li>Other values are formatted as numbers using {@link NumberFormat} followed
by the unit symbol
  *       formatted by {@link org.apache.sis.measure.UnitFormat}.</li>
  * </ul>
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionTransformer.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionTransformer.java
index 04cbbce..1fd3713 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionTransformer.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionTransformer.java
@@ -225,6 +225,6 @@ public final class PositionTransformer extends GeneralDirectPosition {
         if (inverse == null) {
             inverse = (forward != null) ? forward.inverse() : MathTransforms.identity(getDimension());
         }
-        return inverse.transform(this, null);
+        return inverse.transform(this, new GeneralDirectPosition(getCoordinateReferenceSystem()));
     }
 }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.java
index 5c89e58..4af6b0b 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.java
@@ -68,6 +68,11 @@ public final class Resources extends IndexedResourceBundle {
         public static final short AmbiguousEllipsoid_1 = 1;
 
         /**
+         * The azimuth and distance have not been specified.
+         */
+        public static final short AzimuthAndDistanceNotSet = 87;
+
+        /**
          * Can not create objects of type ‘{0}’ from combined URI.
          */
         public static final short CanNotCombineUriAsType_1 = 79;
@@ -197,6 +202,11 @@ public final class Resources extends IndexedResourceBundle {
         public static final short EllipsoidalHeightNotAllowed_1 = 77;
 
         /**
+         * The destination coordinates have not been specified.
+         */
+        public static final short EndPointNotSet = 88;
+
+        /**
          * There is no local registry for version {1} of “{0}” authority. Fallback on
default version
          * for objects creation.
          */
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.properties
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.properties
index 17d7b7b..7d78b45 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.properties
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.properties
@@ -42,6 +42,7 @@ NonConformCRS_3                   = The given \u201c{0}\u201d description
does n
 #
 # Error messages (to be used in exceptions)
 #
+AzimuthAndDistanceNotSet          = The azimuth and distance have not been specified.
 CanNotCombineUriAsType_1          = Can not create objects of type \u2018{0}\u2019 from combined
URI.
 CanNotComputeDerivative           = Can not compute the coordinate operation derivative.
 CanNotConcatenateTransforms_2     = Can not concatenate transforms \u201c{0}\u201d and \u201c{1}\u201d.
@@ -61,6 +62,7 @@ CanNotUseGeodeticParameters_2     = Can not use the {0} geodetic parameters:
{1}
 ColinearAxisDirections_2          = Axis directions {0} and {1} are colinear.
 CoordinateOperationNotFound_2     = Coordinate conversion of transformation from system \u201c{0}\u201d
to \u201c{1}\u201d has not been found.
 DatumOriginShallBeDate            = Origin of temporal datum shall be a date.
+EndPointNotSet                    = The destination coordinates have not been specified.
 DuplicatedParameterName_4         = Name or alias for parameter \u201c{0}\u201d at index
{1} conflict with name \u201c{2}\u201d at index {3}.
 DuplicatedSpatialComponents_1     = Compound coordinate reference systems can not contain
two {0,choice,1#horizontal|2#vertical} components.
 EllipsoidalHeightNotAllowed_1     = Compound coordinate reference systems should not contain
ellipsoidal height. Use a three-dimensional {0,choice,0#geographic|1#projected} system instead.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources_fr.properties
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources_fr.properties
index a4feee5..7f624ba 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources_fr.properties
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources_fr.properties
@@ -47,6 +47,7 @@ NonConformCRS_3                   = La description donn\u00e9e pour \u00ab\u202f
 #
 # Error messages (to be used in exceptions)
 #
+AzimuthAndDistanceNotSet          = Le cap et la distance n\u2019ont pas \u00e9t\u00e9 d\u00e9finis.
 CanNotCombineUriAsType_1          = Ne peut pas cr\u00e9er d\u2019objets de type \u2018{0}\u2019
\u00e0 partir d\u2019un URI combin\u00e9.
 CanNotComputeDerivative           = La d\u00e9riv\u00e9 de l\u2019op\u00e9ration sur les
coordonn\u00e9es ne peut pas \u00eatre calcul\u00e9e.
 CanNotConcatenateTransforms_2     = Les transformations \u00ab\u202f{0}\u202f\u00bb et \u00ab\u202f{1}\u202f\u00bb
ne peuvent pas \u00eatre combin\u00e9es.
@@ -66,6 +67,7 @@ CanNotUseGeodeticParameters_2     = Ne peut pas utiliser les param\u00e8tres
g\u
 ColinearAxisDirections_2          = Les directions d\u2019axes {0} et {1} sont colin\u00e9aires.
 CoordinateOperationNotFound_2     = La conversion ou transformation des coordonn\u00e9es
du syst\u00e8me \u00ab\u202f{0}\u202f\u00bb vers \u00ab\u202f{1}\u202f\u00bb n\u2019a pas
\u00e9t\u00e9 trouv\u00e9e.
 DatumOriginShallBeDate            = L\u2019origine d\u2019un r\u00e9f\u00e9rentiel temporel
doit \u00eatre une date.
+EndPointNotSet                    = Les coordonn\u00e9es de la destination n\u2019ont pas
\u00e9t\u00e9 d\u00e9finies.
 DuplicatedParameterName_4         = Le nom ou un alias pour le param\u00e8tre \u00ab\u202f{0}\u202f\u00bb
\u00e0 l\u2019index {1} duplique le nom \u00ab\u202f{2}\u202f\u00bb \u00e0 l\u2019index {3}.
 DuplicatedSpatialComponents_1     = Un syst\u00e8me de r\u00e9f\u00e9rence des coordonn\u00e9es
ne peut pas contenir deux composantes {0,choice,1#horizontales|2#verticales}.
 EllipsoidalHeightNotAllowed_1     = Un syst\u00e8me de r\u00e9f\u00e9rence des coordonn\u00e9es
ne devrait pas contenir une hauteur ellipso\u00efdale. Utilisez plut\u00f4t un syst\u00e8me
{0,choice,0#g\u00e9ographique|1#projet\u00e9} \u00e0 trois dimensions.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
index 884c4a9..6ad9e73 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
@@ -16,6 +16,12 @@
  */
 package org.apache.sis.referencing;
 
+import java.util.Locale;
+import java.io.IOException;
+import java.io.UncheckedIOException;
+import javax.measure.Unit;
+import javax.measure.quantity.Length;
+
 import org.opengis.referencing.datum.Ellipsoid;
 import org.opengis.referencing.operation.TransformException;
 import org.opengis.referencing.crs.GeographicCRS;
@@ -23,10 +29,15 @@ import org.opengis.referencing.crs.CoordinateReferenceSystem;
 import org.opengis.geometry.coordinate.Position;
 import org.opengis.geometry.DirectPosition;
 
+import org.apache.sis.io.TableAppender;
+import org.apache.sis.measure.Angle;
 import org.apache.sis.measure.Latitude;
+import org.apache.sis.geometry.CoordinateFormat;
 import org.apache.sis.internal.referencing.PositionTransformer;
 import org.apache.sis.internal.referencing.ReferencingUtilities;
+import org.apache.sis.internal.referencing.Resources;
 import org.apache.sis.internal.referencing.Formulas;
+import org.apache.sis.util.resources.Vocabulary;
 import org.apache.sis.util.resources.Errors;
 import org.apache.sis.util.ArgumentChecks;
 
@@ -34,23 +45,18 @@ import static java.lang.Math.*;
 
 
 /**
- * Performs geodetic calculations on a sphere or an ellipsoid.
- * This class calculates the following properties:
- *
- * <ul>
- *   <li>Distance and azimuth between two points.</li>
- *   <li>Point located at a given distance and azimuth from another point.</li>
- * </ul>
- *
- * The calculation uses the following information:
+ * Performs geodetic calculations on a sphere or an ellipsoid. This class computes the distance
between two points,
+ * or conversely the point located at a given distance from another point when navigating
in a given direction.
+ * The distance depends on the path (or track) on Earth surface connecting the two points.
+ * The track can be great circles (shortest path between two points) or rhumb lines (path
with constant heading).
  *
+ * <p>This class uses the following information:</p>
  * <ul>
- *   <li>The {@linkplain #setStartingPosition(Position) starting position}, which is
always considered valid.
+ *   <li>The {@linkplain #setStartPoint(Position) start point}, which is always considered
valid.
  *     It is initially set at (0,0) and can only be changed to another valid value.</li>
- *   <li>One of the following
- *     (the latest specified property overrides the other property and determines what will
be calculated):
+ *   <li>One of the followings (the latest specified property overrides other properties
and determines what will be calculated):
  *     <ul>
- *       <li>the {@linkplain #setDestinationPosition(Position) destination position},
or</li>
+ *       <li>the {@linkplain #setEndPoint(Position) end point}, or</li>
  *       <li>an {@linkplain #setDirection(double, double) azimuth and distance}.</li>
  *     </ul>
  *   </li>
@@ -60,6 +66,9 @@ import static java.lang.Math.*;
  * then a distinct instance of {@code GeodeticCalculator} needs to be created for each thread.
  *
  * @version 1.0
+ *
+ * @see <a href="https://en.wikipedia.org/wiki/Great-circle_navigation">Great-circle
navigation on Wikipedia</a>
+ *
  * @since 1.0
  * @module
  */
@@ -93,34 +102,42 @@ public class GeodeticCalculator {
     private final double radius;
 
     /**
-     * The (<var>latitude</var>, <var>longitude</var>) coordinates
of the first point <strong>in radians</strong>.
-     * This point is set by {@link #setStartingPosition(double, double)}.
+     * The (<var>latitude</var>, <var>longitude</var>) coordinates
of the start point <strong>in radians</strong>.
+     * This point is set by {@link #setStartPoint(double, double)}.
      */
     private double φ1, λ1;
 
     /**
-     * The (<var>latitude</var>, <var>longitude</var>) coordinates
of the destination point <strong>in radians</strong>.
-     * This point is set by {@link #setDestinationPosition(double, double)}.
+     * The (<var>latitude</var>, <var>longitude</var>) coordinates
of the end point <strong>in radians</strong>.
+     * This point is set by {@link #setEndPoint(double, double)}.
      */
     private double φ2, λ2;
 
     /**
-     * The distance and azimuth from the starting point ({@link #φ1},{@link #λ1}) to the
destination point ({@link #φ2},{@link #λ2}).
+     * The azimuth at start point and end point, in radians between -π and +π.
+     * 0° point toward North and values are increasing clockwise.
+     */
+    private double α1, α2;
+
+    /**
+     * The distance from the starting point ({@link #φ1},{@link #λ1}) to the end point
({@link #φ2},{@link #λ2}).
      * The distance is in the same units than ellipsoid axes and the azimuth is in radians.
+     *
+     * @see #getDistanceUnit()
      */
-    private double distance, azimuth;
+    private double distance;
 
     /**
-     * Tells if the destination point is valid.
+     * Tells if the end point is valid.
      * This is {@code false} if {@link #φ2} and {@link #λ2} need to be computed.
      */
-    private boolean isDestinationValid;
+    private boolean isEndPointValid;
 
     /**
-     * Tells if the azimuth and the distance are valid.
-     * This is {@code false} if {@link #distance} and {@link #azimuth} need to be computed.
+     * Tells if the azimuths and the distance are valid.
+     * If {@code false} then {@link #distance}, {@link #α1} and {@link #α2} need to be
computed.
      */
-    private boolean isDirectionValid;
+    private boolean isCourseValid;
 
     /**
      * Constructs a new geodetic calculator expecting coordinates in the supplied CRS.
@@ -168,13 +185,11 @@ public class GeodeticCalculator {
     }
 
     /**
-     * Returns the default CRS of {@link Position} instances. For every method expecting
a {@code Position} argument,
-     * if the {@linkplain DirectPosition#getCoordinateReferenceSystem() CRS of that position}
is unspecified,
-     * then the CRS given by this method is assumed.
-     * Conversely every method returning a {@code Position} value will use this CRS.
-     *
-     * <p>This is the CRS specified at construction time.
-     * This CRS is not necessarily geographic; it may be projected or geocentric.</p>
+     * Returns the Coordinate Reference System (CRS) in which {@code Position}s are represented,
unless otherwise specified.
+     * This is the CRS of all {@link Position} instances returned by methods in this class.
This is also the default CRS
+     * assumed by methods receiving a {@link Position} argument when the given position does
not specify its CRS.
+     * This default CRS is specified at construction time.
+     * It is not necessarily geographic; it may be projected or geocentric.
      *
      * @return the default CRS for {@link Position} instances.
      */
@@ -184,7 +199,8 @@ public class GeodeticCalculator {
 
     /**
      * Returns the coordinate reference system for all methods expecting (φ,λ) as {@code
double} values.
-     * This CRS always use (<var>latitude</var>, <var>longitude</var>)
axis order in degrees.
+     * This CRS always has (<var>latitude</var>, <var>longitude</var>)
axes, in that order and in degrees.
+     * The CRS may contain an additional axis for ellipsoidal height.
      *
      * @return the coordinate reference system of (φ,λ) coordinates.
      */
@@ -193,76 +209,339 @@ public class GeodeticCalculator {
     }
 
     /**
-     * Sets the starting point as geographic coordinates. The azimuth and geodesic distance
values
-     * will be updated as a side effect of this call. They will be recomputed the next time
that
-     * {@link #getAzimuth()} or {@link #getGeodesicDistance()} are invoked.
+     * Sets the starting point as geographic (<var>latitude</var>, <var>longitude</var>)
coordinates.
+     * The {@linkplain #getStartingAzimuth() starting} and {@linkplain #getEndingAzimuth()
ending azimuths},
+     * the {@linkplain #getGeodesicDistance() geodesic distance} and the {@linkplain #getEndPoint()
end point}
+     * are discarded by this method call; some of them will need to be specified again.
      *
-     * @param  latitude   the latitude in decimal degrees between -90 and  +90°
-     * @param  longitude  the longitude in decimal degrees.
+     * @param  latitude   the latitude in degrees between {@value Latitude#MIN_VALUE}° and
{@value Latitude#MAX_VALUE}°.
+     * @param  longitude  the longitude in degrees.
      * @throws IllegalArgumentException if the latitude is out of bounds.
+     *
+     * @see #setEndPoint(double, double)
      */
-    public void setStartingPosition(final double latitude, final double longitude) {
+    public void setStartPoint(final double latitude, final double longitude) {
         ArgumentChecks.ensureBetween("latitude", Latitude.MIN_VALUE, Latitude.MAX_VALUE,
latitude);
         ArgumentChecks.ensureFinite("longitude", longitude);
         φ1 = toRadians(latitude);
         λ1 = toRadians(longitude);
-        isDestinationValid = false;
-        isDirectionValid   = false;
+        isEndPointValid = false;
+        isCourseValid = false;
     }
 
     /**
-     * Sets the starting point in any CRS. The coordinates will be transformed to geographic
coordinates and
-     * given to {@link #setStartingPosition(double, double)}. If the given point is not associated
to a CRS,
-     * then the CRS specified at construction time is assumed.
+     * Sets the starting point as coordinates in arbitrary reference system. This method
transforms the given
+     * coordinates to geographic coordinates, then delegates to {@link #setStartPoint(double,
double)}.
+     * If the given point is not associated to a Coordinate Reference System (CRS), then
this method assumes
+     * the CRS specified at construction time.
      *
-     * @param  position  the starting point in any coordinate reference system.
+     * @param  point  the starting point in any coordinate reference system.
      * @throws TransformException if the coordinates can not be transformed.
+     *
+     * @see #setEndPoint(Position)
      */
-    public void setStartingPosition(final Position position) throws TransformException {
-        final DirectPosition p = userToGeodetic.transform(position.getDirectPosition());
-        setStartingPosition(p.getOrdinate(0), p.getOrdinate(1));
+    public void setStartPoint(final Position point) throws TransformException {
+        final DirectPosition p = userToGeodetic.transform(point.getDirectPosition());
+        setStartPoint(p.getOrdinate(0), p.getOrdinate(1));
     }
 
     /**
      * Returns the starting point in the CRS specified at construction time.
+     * This method returns the last point given to a {@code setStartPoint(…)} method,
+     * transformed to the {@linkplain #getPositionCRS() position CRS}.
+     *
+     * @return the starting point represented in the CRS specified at construction time.
+     * @throws TransformException if the coordinates can not be transformed to {@linkplain
#getPositionCRS() position CRS}.
      *
-     * @return the starting point in user CRS.
-     * @throws TransformException if the coordinates can not be transformed to user CRS.
+     * @see #getEndPoint()
      */
-    public DirectPosition getStartingPosition() throws TransformException {
-        userToGeodetic.setOrdinate(0, toDegrees(φ1));
-        userToGeodetic.setOrdinate(1, toDegrees(λ1));
-        return userToGeodetic.inverseTransform();
+    public DirectPosition getStartPoint() throws TransformException {
+        return geographic(φ1, λ1).inverseTransform();
+    }
+
+    /**
+     * Sets {@link #userToGeodetic} to the given coordinates.
+     * All coordinates in dimension 2 and above (typically the ellipsoidal height) are set
to zero.
+     *
+     * @param  φ  the latitude value to set, in radians.
+     * @param  λ  the longitude value to set, in radians.
+     * @return {@link #userToGeodetic} for convenience.
+     */
+    private PositionTransformer geographic(final double φ, final double λ) {
+        userToGeodetic.setOrdinate(0, toDegrees(φ));
+        userToGeodetic.setOrdinate(1, toDegrees(λ));
+        for (int i=userToGeodetic.getDimension(); --i >= 2;) {
+            userToGeodetic.setOrdinate(i, 0);                   // Set height to ellipsoid
surface.
+        }
+        return userToGeodetic;
     }
 
     /**
-     * Sets the destination point as geographic coordinates. The azimuth and geodesic distance
values
-     * will be updated as a side effect of this call. They will be recomputed the next time
that
-     * {@link #getAzimuth()} or {@link #getGeodesicDistance()} are invoked.
+     * Sets the destination as geographic (<var>latitude</var>, <var>longitude</var>)
coordinates.
+     * The azimuth and geodesic distance values will be updated as an effect of this call.
+     * They will be recomputed next time that {@link #getStartingAzimuth()}, {@link #getEndingAzimuth()}
+     * or {@link #getGeodesicDistance()} is invoked.
      *
-     * @param  latitude   the latitude in decimal degrees between -90 and  +90°
-     * @param  longitude  the longitude in decimal degrees.
+     * @param  latitude   the latitude in degrees between {@value Latitude#MIN_VALUE}° and
{@value Latitude#MAX_VALUE}°.
+     * @param  longitude  the longitude in degrees.
      * @throws IllegalArgumentException if the latitude is out of bounds.
+     *
+     * @see #setStartPoint(double, double)
      */
-    public void setDestinationPosition(final double latitude, final double longitude) {
+    public void setEndPoint(final double latitude, final double longitude) {
         ArgumentChecks.ensureBetween("latitude", Latitude.MIN_VALUE, Latitude.MAX_VALUE,
latitude);
         ArgumentChecks.ensureFinite("longitude", longitude);
         φ2 = toRadians(latitude);
         λ2 = toRadians(longitude);
-        isDestinationValid = false;
-        isDirectionValid   = false;
+        isEndPointValid = true;
+        isCourseValid = false;
     }
 
     /**
-     * Sets the destination point in any CRS. The coordinates will be transformed to geographic
coordinates and
-     * given to {@link #setDestinationPosition(double, double)}. If the given point is not
associated to a CRS,
-     * then the CRS specified at construction time is assumed.
+     * Sets the destination as coordinates in arbitrary reference system. This method transforms
the given
+     * coordinates to geographic coordinates, then delegates to {@link #setEndPoint(double,
double)}.
+     * If the given point is not associated to a Coordinate Reference System (CRS), then
this method assumes
+     * the CRS specified at construction time.
      *
-     * @param  position  the destination point in any coordinate reference system.
+     * @param  position  the destination (end point) in any coordinate reference system.
      * @throws TransformException if the coordinates can not be transformed.
+     *
+     * @see #setStartPoint(Position)
      */
-    public void setDestinationPosition(final Position position) throws TransformException
{
+    public void setEndPoint(final Position position) throws TransformException {
         final DirectPosition p = userToGeodetic.transform(position.getDirectPosition());
-        setDestinationPosition(p.getOrdinate(0), p.getOrdinate(1));
+        setEndPoint(p.getOrdinate(0), p.getOrdinate(1));
+    }
+
+    /**
+     * Returns or computes the destination in the CRS specified at construction time. This
method returns the
+     * point specified in the last call to a {@link #setEndPoint(double, double) setEndPoint(…)}
method, unless
+     * {@link #setDirection(double, double) setDirection(…)} has been invoked more recently.
In the later case,
+     * the end point will be computed from the {@linkplain #getStartPoint() starting point}
and the specified
+     * azimuth and distance.
+     *
+     * @return the destination (end point) represented in the CRS specified at construction
time.
+     * @throws TransformException if the coordinates can not be transformed to {@linkplain
#getPositionCRS() position CRS}.
+     * @throws IllegalStateException if the destination point, azimuth or distance have not
been set.
+     *
+     * @see #getStartPoint()
+     */
+    public DirectPosition getEndPoint() throws TransformException {
+        if (!isEndPointValid) {
+            computeEndPoint();
+        }
+        return geographic(φ2, λ2).inverseTransform();
+    }
+
+    /**
+     * Sets the azimuth and the distance from the {@linkplain #getStartPoint() starting point}.
+     * The direction is relative to geographic North, with values increasing clockwise.
+     * The distance is usually a positive number, but negative numbers are accepted
+     * as courses in direction opposite to the given azimuth.
+     *
+     * <p>The {@linkplain #getEndPoint() end point} will be updated as an effect of
this method call;
+     * it will be recomputed the next time that {@link #getEndPoint()} is invoked.</p>
+     *
+     * @param  azimuth   the starting azimuth in degrees, with 0° toward north and values
increasing clockwise.
+     * @param  distance  the geodesic distance in unit of measurement given by {@link #getDistanceUnit()}.
+     *
+     * @see #getStartingAzimuth()
+     * @see #getGeodesicDistance()
+     */
+    public void setDirection(double azimuth, double distance) {
+        ArgumentChecks.ensureFinite("azimuth",  azimuth);
+        ArgumentChecks.ensureFinite("distance", distance);
+        if (distance < 0) {
+            distance = -distance;
+            azimuth += 180;
+        }
+        α1              = toRadians(IEEEremainder(azimuth, 360));
+        this.distance   = distance;
+        isEndPointValid = false;
+        isCourseValid   = true;
+    }
+
+    /**
+     * Returns the angular heading (relative to geographic North) at the starting point.
+     * This method returns the azimuth normalized to [-180 … +180]° range given in last
call to
+     * <code>{@linkplain #setDirection(double, double) setDirection}(azimuth, …)</code>
method,
+     * unless the {@link #setEndPoint(double, double) setEndPoint(…)} method has been invoked
more recently.
+     * In the later case, the azimuth will be computed from the {@linkplain #getStartPoint
start point}.
+     *
+     * @return the azimuth in degrees from -180° to +180°. 0° is toward North and values
are increasing clockwise.
+     * @throws IllegalStateException if the destination point, azimuth or distance have not
been set.
+     */
+    public double getStartingAzimuth() {
+        if (!isCourseValid) {
+            computeCourse();
+        }
+        return toDegrees(α1);
+    }
+
+    /**
+     * Returns the angular heading (relative to geographic North) at the ending point.
+     *
+     * @return the azimuth in degrees from -180° to +180°. 0° is toward North and values
are increasing clockwise.
+     * @throws IllegalStateException if the destination point, azimuth or distance have not
been set.
+     */
+    public double getEndingAzimuth() {
+        if (!isCourseValid) {
+            computeCourse();
+        } else if (!isEndPointValid) {
+            computeEndPoint();
+        }
+        return toDegrees(α2);
+    }
+
+    /**
+     * Returns the shortest distance from start point to end point.
+     * This is sometime called "great circle" or "orthodromic" distance.
+     * This method returns the absolute distance value set by the last call to
+     * {@link #setDirection(double,double) setDirection(…, distance)} method,
+     * unless the {@link #setEndPoint(double, double) setEndPoint(…)} method has been invoked
more recently.
+     * In the later case, the distance will be computed from the {@linkplain #getStartPoint
start point} to the end point.
+     *
+     * @return The shortest distance in the unit of measurement given by {@link #getDistanceUnit()}.
+     * @throws IllegalStateException if the destination point has not been set.
+     *
+     * @see #getDistanceUnit()
+     */
+    public double getGeodesicDistance() {
+        if (!isCourseValid) {
+            computeCourse();
+        }
+        return distance;
+    }
+
+    /**
+     * Returns the unit of measurement of all distance measurements.
+     * This is the {@linkplain Ellipsoid#getAxisUnit() ellipsoid axis unit}.
+     *
+     * @return the unit of measurement of all distance measurements.
+     *
+     * @see #getGeodesicDistance()
+     */
+    public Unit<Length> getDistanceUnit() {
+        return ellipsoid.getAxisUnit();
+    }
+
+    /**
+     * Computes the end point from the start point, the azimuth and the geodesic distance.
+     * This method should be invoked if the end point or ending azimuth is requested while
+     * {@link #isEndPointValid} is {@code false}.
+     *
+     * <p>The default implementation computes {@link #φ2}, {@link #λ2} and {@link
#α2} using
+     * spherical formulas. Subclasses should override if they can provide ellipsoidal formulas.</p>
+     *
+     * @throws IllegalStateException if the azimuth and the distance have not been set.
+     */
+    void computeEndPoint() {
+        if (!isCourseValid) {
+            throw new IllegalStateException(Resources.format(Resources.Keys.AzimuthAndDistanceNotSet));
+        }
+        final double Δσ    = distance / radius;
+        final double sinΔσ = sin(Δσ);
+        final double cosΔσ = cos(Δσ);
+        final double sinφ1 = sin(φ1);
+        final double cosφ1 = cos(φ1);
+        final double sinα1 = sin(α1);
+        final double cosα1 = cos(α1);
+
+        final double sinΔσ_cosα1 = sinΔσ * cosα1;
+        final double Δλy = sinΔσ * sinα1;
+        final double Δλx = cosΔσ*cosφ1 - sinφ1*sinΔσ_cosα1;
+        final double Δλ  = atan2(Δλy, Δλx);
+
+        final double sinφ2 = sinφ1*cosΔσ + cosφ1*sinΔσ_cosα1;       // First estimation
of φ2.
+        φ2 = atan(sinφ2 / hypot(Δλx, Δλy));                         // Improve accuracy
close to poles.
+        λ2 = IEEEremainder(λ1 + Δλ, 2*PI);
+        α2 = atan2(sinα1, cosΔσ*cosα1 - sinφ1/cosφ1 * sinΔσ);
+        isEndPointValid = true;
+    }
+
+    /**
+     * Computes the geodetic distance and azimuths from the start point and end point.
+     * This method should be invoked if the distance or an azimuth is requested while
+     * {@link #isCourseValid} is {@code false}.
+     *
+     * @throws IllegalStateException if the distance or azimuth has not been set.
+     */
+    private void computeCourse() {
+        if (!isEndPointValid) {
+            throw new IllegalStateException(Resources.format(Resources.Keys.EndPointNotSet));
+        }
+        final double Δλ    = λ2 - λ1;           // No need to reduce to −π … +π
range.
+        final double sinΔλ = sin(Δλ);
+        final double cosΔλ = cos(Δλ);
+        final double sinφ1 = sin(φ1);
+        final double cosφ1 = cos(φ1);
+        final double sinφ2 = sin(φ2);
+        final double cosφ2 = cos(φ2);
+
+        final double cosφ1_sinφ2 = cosφ1 * sinφ2;
+        final double cosφ2_sinφ1 = cosφ2 * sinφ1;
+        final double α1y = cosφ2 * sinΔλ;
+        final double α1x = cosφ1_sinφ2 - cosφ2_sinφ1*cosΔλ;
+
+        α1 = atan2(α1y, α1x);
+        α2 = atan2(cosφ1*sinΔλ, cosφ1_sinφ2*cosΔλ - cosφ2_sinφ1);
+        /*
+         * Δσ = acos(sinφ₁⋅sinφ₂ + cosφ₁⋅cosφ₂⋅cosΔλ) is a first estimation
inaccurate for small distances.
+         * Δσ = atan2(…) computes the same value but with better accuracy.
+         */
+        double Δσ = sinφ1*sinφ2 + cosφ1*cosφ2*cosΔλ;        // Actually Δσ = acos(…).
+        Δσ = atan2(hypot(α1x, α1y), Δσ);
+        distance = radius * Δσ;
+        isCourseValid = true;
+    }
+
+    /**
+     * Returns a string representation of start point, end point, azimuths and distance.
+     *
+     * @return a string representation of this calculator state.
+     */
+    @Override
+    public String toString() {
+        final Locale     locale    = Locale.getDefault();
+        final Vocabulary resources = Vocabulary.getResources(locale);
+        final StringBuilder buffer = new StringBuilder();
+        final String lineSeparator = System.lineSeparator();
+        final CoordinateReferenceSystem crs = userToGeodetic.getCoordinateReferenceSystem();
+        final boolean isGeographic = crs.equals(userToGeodetic.defaultCRS);
+        try {
+            resources.appendLabel(Vocabulary.Keys.CoordinateRefSys, buffer);
+            buffer.append(' ').append(crs.getName().getCode()).append(lineSeparator);
+            final TableAppender table = new TableAppender(buffer, " │ ");
+            table.appendHorizontalSeparator();
+            table.nextColumn(); if (isGeographic) table.append(resources.getString(Vocabulary.Keys.Latitude));
+            table.nextColumn(); if (isGeographic) table.append(resources.getString(Vocabulary.Keys.Longitude));
+            for (int i=crs.getCoordinateSystem().getDimension(); --i >= 2;) {
+                table.nextColumn();     // Insert space for additional coordinates, e.g.
ellipsoidal height.
+            }
+            table.nextColumn();
+            table.append(resources.getString(Vocabulary.Keys.Azimuth)).nextLine();
+            final CoordinateFormat cf = new CoordinateFormat(locale, null);
+            cf.setSeparator("\t");      // For distributing coordinate values on different
columns.
+            boolean endPoint = false;
+            do {
+                table.append(resources.getString(endPoint ? Vocabulary.Keys.EndPoint : Vocabulary.Keys.StartPoint))
+                     .nextColumn();
+                try {
+                    cf.format(endPoint ? getEndPoint() : getStartPoint(), table);
+                    table.nextColumn();
+                    table.append(new Angle(endPoint ? getEndingAzimuth() : getStartingAzimuth()).toString());
+                } catch (IllegalStateException | TransformException e) {
+                    // Ignore.
+                }
+                table.nextLine();
+            } while ((endPoint = !endPoint) == true);
+            table.appendHorizontalSeparator();
+            table.flush();
+            resources.appendLabel(Vocabulary.Keys.GeodesicDistance, buffer);
+        } catch (IOException e) {
+            throw new UncheckedIOException(e);      // Should never happen since we are writting
in a StringBuilder.
+        }
+        buffer.append(String.format(locale, " %f %s", distance, getDistanceUnit()));
+        return buffer.toString();
     }
 }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
new file mode 100644
index 0000000..cafb080
--- /dev/null
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
@@ -0,0 +1,83 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF 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 org.apache.sis.referencing;
+
+import org.opengis.geometry.DirectPosition;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.measure.Units;
+import org.apache.sis.test.TestCase;
+import org.junit.Test;
+
+import static org.junit.Assert.*;
+
+
+/**
+ * Tests {@link GeodeticCalculator}.
+ *
+ * @version 1.0
+ * @since 1.0
+ * @module
+ */
+public final strictfp class GeodeticCalculatorTest extends TestCase {
+    /**
+     * Tests spherical formulas with the example given in Wikipedia.
+     * This computes the great circle route from Valparaíso (33°N 71.6W) to Shanghai (31.4°N
121.8°E).
+     *
+     * @throws TransformException if an error occurred while transforming coordinates.
+     */
+    @Test
+    public void testWikipediaExample() throws TransformException {
+        final GeodeticCalculator c = new GeodeticCalculator(CommonCRS.SPHERE.geographic());
+        c.setStartPoint(-33.0, -71.6);          // Valparaíso
+        c.setEndPoint  ( 31.4, 121.8);          // Shanghai
+        /*
+         * Wikipedia example gives:
+         *
+         *     Δλ = −166.6°
+         *     α₁ = −94.41°
+         *     α₂ = −78.42°
+         *     Δσ = 168.56°   →    taking R = 6371 km, the distance is 18743 km.
+         */
+        assertEquals(Units.METRE,         c.getDistanceUnit());
+        assertEquals("α₁",        -94.41, c.getStartingAzimuth(), 0.005);
+        assertEquals("α₂",        -78.42, c.getEndingAzimuth(),   0.005);
+        assertEquals("distance",   18743, c.getGeodesicDistance() / 1000, 0.5);
+        assertPositionEquals(31.4, 121.8, c.getEndPoint(), 1E-12);                  // Should
be the specified value.
+        /*
+         * Keep start point unchanged, but set above azimuth and distance.
+         * Verify that we get the Shanghai coordinates.
+         */
+        c.setDirection(-94.41, 18743000);
+        assertEquals("α₁",        -94.41, c.getStartingAzimuth(), 1E-12);           //
Should be the specified value.
+        assertEquals("α₂",        -78.42, c.getEndingAzimuth(),   0.01);
+        assertEquals("distance",   18743, c.getGeodesicDistance() / 1000, STRICT);  // Should
be the specified value.
+        assertPositionEquals(31.4, 121.8, c.getEndPoint(), 0.01);
+    }
+
+    /**
+     * Verifies that the given point is equals to the given latitude and longitude.
+     *
+     * @param φ  the expected latitude value, in degrees.
+     * @param λ  the expected longitude value, in degrees.
+     * @param p  the position to verify.
+     * @param ε  the tolerance.
+     */
+    private static void assertPositionEquals(final double φ, final double λ, final DirectPosition
p, final double ε) {
+        assertEquals("φ", φ, p.getOrdinate(0), ε);
+        assertEquals("λ", λ, p.getOrdinate(1), ε);
+    }
+}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
index 6aed0e3..dec780f 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
@@ -219,7 +219,7 @@ import org.junit.BeforeClass;
     org.apache.sis.referencing.factory.GIGS2008.class,
     org.apache.sis.referencing.factory.GIGS2009.class,
 
-    // Following tests use indirectly EPSG factory.
+    // Following tests may use indirectly EPSG factory.
     org.apache.sis.referencing.CommonCRSTest.class,
     org.apache.sis.referencing.factory.CommonAuthorityFactoryTest.class,
     org.apache.sis.referencing.factory.AuthorityFactoryProxyTest.class,
@@ -234,6 +234,7 @@ import org.junit.BeforeClass;
     org.apache.sis.referencing.AuthorityFactoriesTest.class,
     org.apache.sis.referencing.cs.CodesTest.class,
     org.apache.sis.referencing.CRSTest.class,
+    org.apache.sis.referencing.GeodeticCalculatorTest.class,
     org.apache.sis.internal.referencing.DefinitionVerifierTest.class,
     org.apache.sis.internal.referencing.CoordinateOperationsTest.class,
 
diff --git a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.java
b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.java
index 6c548eb..ce4fc0d 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.java
@@ -107,6 +107,11 @@ public final class Vocabulary extends IndexedResourceBundle {
         public static final short AxisChanges = 8;
 
         /**
+         * Azimuth
+         */
+        public static final short Azimuth = 168;
+
+        /**
          * Band {0}
          */
         public static final short Band_1 = 161;
@@ -337,6 +342,11 @@ public final class Vocabulary extends IndexedResourceBundle {
         public static final short EndDate = 134;
 
         /**
+         * End point
+         */
+        public static final short EndPoint = 169;
+
+        /**
          * {0} entr{0,choice,0#y|2#ies}
          */
         public static final short EntryCount_1 = 121;
@@ -382,6 +392,11 @@ public final class Vocabulary extends IndexedResourceBundle {
         public static final short GeocentricRadius = 44;
 
         /**
+         * Geodesic distance
+         */
+        public static final short GeodesicDistance = 170;
+
+        /**
          * Geodetic dataset
          */
         public static final short GeodeticDataset = 45;
@@ -747,6 +762,11 @@ public final class Vocabulary extends IndexedResourceBundle {
         public static final short StartDate = 139;
 
         /**
+         * Start point
+         */
+        public static final short StartPoint = 171;
+
+        /**
          * Subset of {0}
          */
         public static final short SubsetOf_1 = 95;
diff --git a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.properties
b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.properties
index 2933b17..4110944 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.properties
+++ b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.properties
@@ -24,6 +24,7 @@ AngularMinutes          = Minutes
 AngularSeconds          = Seconds
 Attributes              = Attributes
 AxisChanges             = Axis changes
+Azimuth                 = Azimuth
 Band_1                  = Band {0}
 BarometricAltitude      = Barometric altitude
 Cardinality             = Cardinality
@@ -70,6 +71,7 @@ Ellipsoid               = Ellipsoid
 EllipsoidChange         = Ellipsoid change
 EllipsoidalHeight       = Ellipsoidal height
 EndDate                 = End date
+EndPoint                = End point
 EntryCount_1            = {0} entr{0,choice,0#y|2#ies}
 Envelope                = Envelope
 Errors                  = Errors
@@ -79,6 +81,7 @@ FillValue               = Fill value
 Geocentric              = Geocentric
 GeocentricRadius        = Geocentric radius
 GeocentricConversion    = Geocentric conversion
+GeodesicDistance        = Geodesic distance
 GeodeticDataset         = Geodetic dataset
 GeographicIdentifier    = Geographic identifier
 GridExtent              = Grid extent
@@ -152,6 +155,7 @@ Source                  = Source
 SouthBound              = South bound
 StandardDeviation       = Standard deviation
 StartDate               = Start date
+StartPoint              = Start point
 SubsetOf_1              = Subset of {0}
 SupersededBy_1          = Superseded by {0}.
 TemporaryFiles          = Temporary files
diff --git a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary_fr.properties
b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary_fr.properties
index 4c5bebe..6015151 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary_fr.properties
+++ b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary_fr.properties
@@ -31,6 +31,7 @@ AngularMinutes          = Minutes
 AngularSeconds          = Secondes
 Attributes              = Attributs
 AxisChanges             = Changements d\u2019axes
+Azimuth                 = Azimut
 Band_1                  = Bande {0}
 BarometricAltitude      = Altitude barom\u00e9trique
 Cardinality             = Cardinalit\u00e9
@@ -78,6 +79,7 @@ EllipsoidChange         = Changement d\u2019ellipso\u00efde
 EllipsoidalHeight       = Hauteur ellipso\u00efdale
 EntryCount_1            = {0} entr\u00e9e{0,choice,0#|2#s}
 EndDate                 = Date de fin
+EndPoint                = Point d\u2019arriv\u00e9
 Envelope                = Enveloppe
 Errors                  = Erreurs
 Exit                    = Quitter
@@ -86,6 +88,7 @@ FillValue               = Valeur de remplissage
 Geocentric              = G\u00e9ocentrique
 GeocentricRadius        = Rayon g\u00e9ocentrique
 GeocentricConversion    = Conversion g\u00e9ocentrique
+GeodesicDistance        = Distance g\u00e9od\u00e9sique
 GeodeticDataset         = Base de donn\u00e9es g\u00e9od\u00e9sique
 GeographicIdentifier    = Identifiant g\u00e9ographique
 GridExtent              = \u00c9tendue de la grille
@@ -159,6 +162,7 @@ Source                  = Source
 SouthBound              = Limite sud
 StandardDeviation       = \u00c9cart type
 StartDate               = Date de d\u00e9part
+StartPoint              = Point de d\u00e9part
 SubsetOf_1              = Sous-ensemble de {0}
 SupersededBy_1          = Remplac\u00e9 par {0}.
 TemporaryFiles          = Fichiers temporaires


Mime
View raw message