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: Initial (and incomplete) implementation of loxodromic distance using spherical formulas.
Date Wed, 15 May 2019 19:20:17 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 9a15749  Initial (and incomplete) implementation of loxodromic distance using spherical
formulas.
9a15749 is described below

commit 9a157499e3753835cc67deb5d825bd48a606cb41
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed May 15 21:19:30 2019 +0200

    Initial (and incomplete) implementation of loxodromic distance using spherical formulas.
---
 .../apache/sis/referencing/GeodeticCalculator.java | 125 ++++++++++++++++-----
 .../sis/referencing/GeodeticCalculatorTest.java    |  17 +--
 2 files changed, 103 insertions(+), 39 deletions(-)

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 cbaff00..4d513c3 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
@@ -131,7 +131,7 @@ public class GeodeticCalculator {
     private double α1, α2;
 
     /**
-     * The distance from the starting point ({@link #φ1},{@link #λ1}) to the end point
({@link #φ2},{@link #λ2}).
+     * The shortest 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 #GEODESIC_DISTANCE
@@ -140,6 +140,15 @@ public class GeodeticCalculator {
     private double geodesicDistance;
 
     /**
+     * Length of the rhumb line 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.
+     *
+     * @see #RHUMBLINE_LENGTH
+     * @see #getDistanceUnit()
+     */
+    private double rhumblineLength;
+
+    /**
      * A bitmask specifying which information are valid. For example if the {@link #END_POINT}
bit is not set,
      * then {@link #φ2} and {@link #λ2} need to be computed, which implies the computation
of {@link #α1} and
      * {@link #α2} as well. If the {@link #GEODESIC_DISTANCE} bit is not set, then {@link
#geodesicDistance}
@@ -150,7 +159,8 @@ public class GeodeticCalculator {
     /**
      * Bitmask specifying which information are valid.
      */
-    private static final int START_POINT = 1, END_POINT = 2, STARTING_AZIMUTH = 4, ENDING_AZIMUTH
= 8, GEODESIC_DISTANCE = 16;
+    private static final int START_POINT = 1, END_POINT = 2, STARTING_AZIMUTH = 4, ENDING_AZIMUTH
= 8,
+            GEODESIC_DISTANCE = 16, RHUMBLINE_LENGTH = 32;
 
     /**
      * Constructs a new geodetic calculator expecting coordinates in the supplied CRS.
@@ -189,6 +199,13 @@ public class GeodeticCalculator {
     }
 
     /**
+     * Returns {@code true} if at least one of the properties identified by the given mask
is invalid.
+     */
+    private boolean isInvalid(final int mask) {
+        return (validity & mask) != mask;
+    }
+
+    /**
      * 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.
@@ -215,7 +232,8 @@ public class GeodeticCalculator {
     /**
      * 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}
+     * the {@linkplain #getEndPoint() end point}, the {@linkplain #getGeodesicDistance()
geodesic distance}
+     * and the {@linkplain #getRhumblineLength() rhumb line length}
      * are discarded by this method call; some of them will need to be specified again.
      *
      * @param  latitude   the latitude in degrees between {@value Latitude#MIN_VALUE}° and
{@value Latitude#MAX_VALUE}°.
@@ -261,11 +279,10 @@ public class GeodeticCalculator {
      * @see #getEndPoint()
      */
     public DirectPosition getStartPoint() throws TransformException {
-        if ((validity & START_POINT) != 0) {
-            return geographic(φ1, λ1).inverseTransform();
-        } else {
+        if (isInvalid(START_POINT)) {
             throw new IllegalStateException(Resources.format(Resources.Keys.StartOrEndPointNotSet_1,
0));
         }
+        return geographic(φ1, λ1).inverseTransform();
     }
 
     /**
@@ -288,7 +305,8 @@ public class GeodeticCalculator {
     /**
      * Sets the destination as geographic (<var>latitude</var>, <var>longitude</var>)
coordinates.
      * The {@linkplain #getStartingAzimuth() starting azimuth}, {@linkplain #getEndingAzimuth()
ending azimuth}
-     * and {@linkplain #getGeodesicDistance() geodesic distance} will be updated as an effect
of this call.
+     * {@linkplain #getGeodesicDistance() geodesic distance} and {@linkplain #getRhumblineLength()
rhumb line length}
+     * will be updated as an effect of this call.
      *
      * @param  latitude   the latitude in degrees between {@value Latitude#MIN_VALUE}° and
{@value Latitude#MAX_VALUE}°.
      * @param  longitude  the longitude in degrees.
@@ -302,7 +320,7 @@ public class GeodeticCalculator {
         φ2 = toRadians(latitude);
         λ2 = toRadians(longitude);
         validity |= END_POINT;
-        validity &= ~(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE);
+        validity &= ~(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE | RHUMBLINE_LENGTH);
     }
 
     /**
@@ -336,7 +354,7 @@ public class GeodeticCalculator {
      * @see #getStartPoint()
      */
     public DirectPosition getEndPoint() throws TransformException {
-        if ((validity & END_POINT) == 0) {
+        if (isInvalid(END_POINT)) {
             computeEndPoint();
         }
         return geographic(φ2, λ2).inverseTransform();
@@ -345,7 +363,8 @@ public class GeodeticCalculator {
     /**
      * Sets the angular heading (relative to geographic North) at the starting point.
      * Azimuth is relative to geographic North with values increasing clockwise.
-     * The {@linkplain #getEndPoint() end point} and {@linkplain #getEndingAzimuth() ending
azimuth}
+     * The {@linkplain #getEndingAzimuth() ending azimuth}, {@linkplain #getEndPoint() end
point}
+     * and {@linkplain #getRhumblineLength() rhumb line length}
      * will be updated as an effect of this method call.
      *
      * @param  azimuth  the starting azimuth in degrees, with 0° toward north and values
increasing clockwise.
@@ -356,7 +375,7 @@ public class GeodeticCalculator {
         ArgumentChecks.ensureFinite("azimuth", azimuth);
         α1 = toRadians(IEEEremainder(azimuth, 360));
         validity |= STARTING_AZIMUTH;
-        validity &= ~(END_POINT | ENDING_AZIMUTH);
+        validity &= ~(END_POINT | ENDING_AZIMUTH | RHUMBLINE_LENGTH);
     }
 
     /**
@@ -370,7 +389,7 @@ public class GeodeticCalculator {
      * @throws IllegalStateException if the end point, azimuth or distance have not been
set.
      */
     public double getStartingAzimuth() {
-        if ((validity & STARTING_AZIMUTH) == 0) {
+        if (isInvalid(STARTING_AZIMUTH)) {
             computeDistance();
         }
         return toDegrees(α1);
@@ -386,8 +405,8 @@ public class GeodeticCalculator {
      * @throws IllegalStateException if the destination point, azimuth or distance have not
been set.
      */
     public double getEndingAzimuth() {
-        if ((validity & ENDING_AZIMUTH) == 0) {
-            if ((validity & END_POINT) == 0) {
+        if (isInvalid(ENDING_AZIMUTH)) {
+            if (isInvalid(END_POINT)) {
                 computeEndPoint();                      // Compute also ending azimuth from
start point and distance.
             } else {
                 computeDistance();                         // Compute also ending azimuth
from start point and end point.
@@ -397,8 +416,9 @@ public class GeodeticCalculator {
     }
 
     /**
-     * Sets the geodesic distance from the start point to the end point. The {@linkplain
#getEndPoint() end point}
-     * and {@linkplain #getEndingAzimuth() ending azimuth} will be updated as an effect of
this method call.
+     * Sets the geodesic distance from the start point to the end point. The {@linkplain
#getEndPoint() end point},
+     * {@linkplain #getEndingAzimuth() ending azimuth} and {@linkplain #getRhumblineLength()
rhumb line length}
+     * will be updated as an effect of this method call.
      *
      * @param  distance  the geodesic distance in unit of measurement given by {@link #getDistanceUnit()}.
      *
@@ -408,7 +428,7 @@ public class GeodeticCalculator {
         ArgumentChecks.ensurePositive("distance", distance);
         geodesicDistance = distance;
         validity |= GEODESIC_DISTANCE;
-        validity &= ~(END_POINT | ENDING_AZIMUTH);
+        validity &= ~(END_POINT | ENDING_AZIMUTH | RHUMBLINE_LENGTH);
     }
 
     /**
@@ -418,18 +438,32 @@ public class GeodeticCalculator {
      * the distance will be computed from the {@linkplain #getStartPoint() start point} and
current 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.
+     * @throws IllegalStateException if a point has not been set.
      *
      * @see #getDistanceUnit()
      */
     public double getGeodesicDistance() {
-        if ((validity & GEODESIC_DISTANCE) == 0) {
+        if (isInvalid(GEODESIC_DISTANCE)) {
             computeDistance();
         }
         return geodesicDistance;
     }
 
     /**
+     * Returns or computes the length of rhumb line (part of constant heading) from start
point to end point.
+     * This is sometime called "loxodrome". This is <strong>not</strong> the
shortest path between two points.
+     *
+     * @return length of rhumb line in the unit of measurement given by {@link #getDistanceUnit()}.
+     * @throws IllegalStateException if a point has not been set.
+     */
+    public double getRhumblineLength() {
+        if (isInvalid(RHUMBLINE_LENGTH)) {
+            computeRhumbLine();
+        }
+        return rhumblineLength;
+    }
+
+    /**
      * Returns the unit of measurement of all distance measurements.
      * This is the {@linkplain Ellipsoid#getAxisUnit() ellipsoid axis unit}.
      *
@@ -442,6 +476,41 @@ public class GeodeticCalculator {
     }
 
     /**
+     * Computes the length of rhumb line from start point to end point.
+     *
+     * @see <a href="https://en.wikipedia.org/wiki/Rhumb_line">Rhumb line on Wikipedia</a>
+     *
+     * @todo if {@literal Δλ > 180}, must split in two segments.
+     */
+    private void computeRhumbLine() {
+        if (isInvalid(START_POINT | END_POINT)) {
+            throw new IllegalStateException(Resources.format(
+                    Resources.Keys.StartOrEndPointNotSet_1, Integer.signum(validity &
START_POINT)));
+        }
+        final double Δλ = λ2 - λ1;
+        final double Δφ = φ2 - φ2;
+        double factor;
+        if (abs(Δφ) < 1E-7) {
+            factor = Δλ * cos((φ1 + φ2)/2);
+        } else {
+            /*
+             * Inverse of Gudermannian function is log(tan(π/4 + φ/2)).
+             * The loxodrome formula involves the following difference:
+             *
+             *   ΔG  =  log(tan(π/4 + φ₁/2)) - log(tan(π/4 + φ₂/2))
+             *       =  log(tan(π/4 + φ₁/2) / tan(π/4 + φ₂/2))
+             *
+             * Note that ΔG=0 if φ₁=φ₂, which implies cos(α)=0.
+             */
+            final double ΔG = log(tan(PI/4 + φ1/2) / tan(PI/4 + φ2/2));
+            final double α = atan(Δλ / ΔG);
+            factor = abs(Δφ) / cos(α);
+        }
+        rhumblineLength = radius * factor;
+        validity |= RHUMBLINE_LENGTH;
+    }
+
+    /**
      * 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 #STARTING_AZIMUTH}, {@link #ENDING_AZIMUTH} or {@link #GEODESIC_DISTANCE}
@@ -456,9 +525,7 @@ public class GeodeticCalculator {
      * @throws IllegalStateException if the distance or azimuth has not been set.
      */
     private void computeDistance() {
-        if ((validity & (START_POINT | END_POINT))
-                     != (START_POINT | END_POINT))
-        {
+        if (isInvalid(START_POINT | END_POINT)) {
             throw new IllegalStateException(Resources.format(
                     Resources.Keys.StartOrEndPointNotSet_1, Integer.signum(validity &
START_POINT)));
         }
@@ -498,10 +565,8 @@ public class GeodeticCalculator {
      * @throws IllegalStateException if the azimuth and the distance have not been set.
      */
     void computeEndPoint() {
-        if ((validity & (START_POINT | STARTING_AZIMUTH | GEODESIC_DISTANCE))
-                     != (START_POINT | STARTING_AZIMUTH | GEODESIC_DISTANCE))
-        {
-            throw new IllegalStateException((validity & START_POINT) == 0
+        if (isInvalid(START_POINT | STARTING_AZIMUTH | GEODESIC_DISTANCE)) {
+            throw new IllegalStateException(isInvalid(START_POINT)
                     ? Resources.format(Resources.Keys.StartOrEndPointNotSet_1, 0)
                     : Resources.format(Resources.Keys.AzimuthAndDistanceNotSet));
         }
@@ -534,7 +599,7 @@ public class GeodeticCalculator {
      * @see #setStartPoint(double, double)
      */
     public void moveToEndPoint() {
-        if ((validity & END_POINT) == 0) {
+        if (isInvalid(END_POINT)) {
             computeEndPoint();
         }
         φ1 = φ2;
@@ -573,10 +638,8 @@ public class GeodeticCalculator {
      * @throws IllegalStateException if some required properties have not been specified.
      */
     public Shape toGeodesicPath2D(double tolerance) throws TransformException {
-        if ((validity & (START_POINT | STARTING_AZIMUTH | END_POINT | ENDING_AZIMUTH
| GEODESIC_DISTANCE))
-                     != (START_POINT | STARTING_AZIMUTH | END_POINT | ENDING_AZIMUTH | GEODESIC_DISTANCE))
-        {
-            if ((validity & END_POINT) == 0) {
+        if (isInvalid(START_POINT | STARTING_AZIMUTH | END_POINT | ENDING_AZIMUTH | GEODESIC_DISTANCE))
{
+            if (isInvalid(END_POINT)) {
                 computeEndPoint();
             } else {
                 computeDistance();
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
index 1d2daff..295573b 100644
--- 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
@@ -131,11 +131,12 @@ public final strictfp class GeodeticCalculatorTest extends TestCase
{
             final double x = 180 * random.nextDouble();
             c.setEndPoint(0, x);
             assertEquals(x * r, c.getGeodesicDistance(), Formulas.LINEAR_TOLERANCE);
+            assertEquals(x * r, c.getRhumblineLength(),  Formulas.LINEAR_TOLERANCE);
         }
     }
 
     /**
-     * Tests spherical formulas with the example given in Wikipedia.
+     * Tests {@link GeodeticCalculator#getGeodesicDistance()} and azimuths 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.
@@ -143,7 +144,7 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
      * @see <a href="https://en.wikipedia.org/wiki/Great-circle_navigation#Example">Great-circle
navigation on Wikipedia</a>
      */
     @Test
-    public void testWikipediaExample() throws TransformException {
+    public void testGeodesic() throws TransformException {
         final GeodeticCalculator c = create(false);
         c.setStartPoint(-33.0, -71.6);          // Valparaíso
         c.setEndPoint  ( 31.4, 121.8);          // Shanghai
@@ -175,12 +176,12 @@ public final strictfp class GeodeticCalculatorTest extends TestCase
{
     /**
      * Tests geodetic calculator involving a coordinate operation.
      * This test uses a simple CRS with only the axis order interchanged.
-     * The coordinates are the same than {@link #testWikipediaExample()}.
+     * The coordinates are the same than {@link #testGeodesic()}.
      *
      * @throws TransformException if an error occurred while transforming coordinates.
      */
     @Test
-    @DependsOnMethod("testWikipediaExample")
+    @DependsOnMethod("testGeodesic")
     public void testUsingTransform() throws TransformException {
         final GeodeticCalculator c = create(true);
         final double φ = -33.0;
@@ -271,17 +272,17 @@ public final strictfp class GeodeticCalculatorTest extends TestCase
{
             180.00,  10000,  14142
         };
         final GeodeticCalculator c = create(false);
-        final double toTestValue = (60 * NAUTICAL_MILE * 180/PI) / 6371007 / 1000;
+        final double toTestValue = (60 * NAUTICAL_MILE * 180/PI) / 6371000 / 1000;
         for (int i=0; i<data.length; i+=3) {
             c.setStartPoint(45, 0);
             c.setEndPoint(45, data[i]);
             double geodesic  = c.getGeodesicDistance();
-//          double rhumbLine = c.getRhumbLineDistance();
+            double rhumbLine = c.getRhumblineLength();
             geodesic  *= toTestValue;
-//          rhumbLine *= toTestValue;
+            rhumbLine *= toTestValue;
             final double tolerance = data[i+1] / 1000;           // In kilometres.
             assertEquals("Geodesic distance",   data[i+1], geodesic,  tolerance);
-//          assertEquals("Rhumb line distance", data[i+2], rhumbLine, tolerance);
+            assertEquals("Rhumb line distance", data[i+2], rhumbLine, tolerance);
             assertEquals("Distance measured along geodesic path", geodesic, length(c) * toTestValue,
tolerance * 20);
         }
     }


Mime
View raw message