sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 03/03: Add GeodeticCalculator.createProjectionAroundStart() method.
Date Wed, 01 Jul 2020 14:33:23 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

commit ffb556af0af6c74bbef68fa1015136a8730b6db0
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed Jul 1 16:31:57 2020 +0200

    Add GeodeticCalculator.createProjectionAroundStart() method.
---
 .../sis/referencing/GeodesicsOnEllipsoid.java      |   9 ++
 .../apache/sis/referencing/GeodeticCalculator.java | 138 ++++++++++++++++++++-
 .../sis/referencing/GeodeticCalculatorTest.java    |  35 +++++-
 3 files changed, 180 insertions(+), 2 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
index 98302b0..753379a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
@@ -993,4 +993,13 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
         final double cosθ = (cosφ + sinφ) * (cosφ - sinφ);      // cos(2φ)  =  cos²φ
- sin²φ
         return R0*φ + sinθ*(R2 + cosθ*(R4 + cosθ*R6));
     }
+
+    /**
+     * The operation method to use for creating a map projection. For the ellipsoidal case
we use EPSG::9832.
+     * According EPSG documentation the precision is acceptable withing 800 km of projection
natural origin.
+     */
+    @Override
+    final String getProjectionMethod() {
+        return "Modified Azimuthal Equidistant";
+    }
 }
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 4dd9582..9872ee8 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
@@ -23,10 +23,14 @@ import java.io.IOException;
 import java.io.UncheckedIOException;
 import java.text.NumberFormat;
 import javax.measure.Unit;
+import javax.measure.UnitConverter;
 import javax.measure.quantity.Length;
 
+import org.opengis.util.FactoryException;
 import org.opengis.referencing.datum.Ellipsoid;
 import org.opengis.referencing.operation.Matrix;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.MathTransformFactory;
 import org.opengis.referencing.operation.TransformException;
 import org.opengis.referencing.crs.GeographicCRS;
 import org.opengis.referencing.crs.CoordinateReferenceSystem;
@@ -38,12 +42,18 @@ import org.apache.sis.measure.Latitude;
 import org.apache.sis.measure.Units;
 import org.apache.sis.measure.Quantities;
 import org.apache.sis.geometry.CoordinateFormat;
+import org.apache.sis.parameter.Parameters;
+import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.referencing.operation.transform.DefaultMathTransformFactory;
+import org.apache.sis.internal.referencing.provider.MapProjection;
 import org.apache.sis.internal.referencing.PositionTransformer;
 import org.apache.sis.internal.referencing.ReferencingUtilities;
 import org.apache.sis.internal.referencing.j2d.ShapeUtilities;
 import org.apache.sis.internal.referencing.j2d.Bezier;
 import org.apache.sis.internal.referencing.Resources;
 import org.apache.sis.internal.referencing.Formulas;
+import org.apache.sis.internal.system.DefaultFactories;
+import org.apache.sis.internal.util.Constants;
 import org.apache.sis.internal.util.Numerics;
 import org.apache.sis.util.resources.Vocabulary;
 import org.apache.sis.util.resources.Errors;
@@ -52,6 +62,8 @@ import org.apache.sis.io.TableAppender;
 
 import static java.lang.Math.*;
 import static org.apache.sis.internal.metadata.ReferencingServices.NAUTICAL_MILE;
+import static org.apache.sis.internal.referencing.provider.ModifiedAzimuthalEquidistant.LATITUDE_OF_ORIGIN;
+import static org.apache.sis.internal.referencing.provider.ModifiedAzimuthalEquidistant.LONGITUDE_OF_ORIGIN;
 
 
 /**
@@ -113,7 +125,7 @@ import static org.apache.sis.internal.metadata.ReferencingServices.NAUTICAL_MILE
  * then a distinct instance of {@code GeodeticCalculator} needs to be created for each thread.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  * @since   1.0
  * @module
  */
@@ -1072,6 +1084,130 @@ public class GeodeticCalculator {
     }
 
     /**
+     * The factory for map projections created by {@link #createProjectionAroundStart()},
fetched when first needed.
+     * {@link DefaultMathTransformFactory#caching(boolean) Caching is disabled} on this factory
because profiling
+     * shows that {@link DefaultMathTransformFactory#unique(MathTransform)} consumes a lot
of time when projections
+     * are created frequently. Since each projection is specific to current {@linkplain #getStartPoint()
start point},
+     * they are unlikely to be shared anyway.
+     *
+     * @see #createProjectionAroundStart()
+     */
+    private DefaultMathTransformFactory projectionFactory;
+
+    /**
+     * The provider of "Azimuthal Equidistant (Spherical)" or "Modified Azimuthal Equidistant"
projection.
+     * Usually it is not necessary to keep a reference to the provider because {@link #projectionFactory}
+     * finds them automatically. However by keeping a reference to it, we save the search
phase.
+     *
+     * @see #createProjectionAroundStart()
+     */
+    private MapProjection projectionProvider;
+
+    /**
+     * Parameters of the projection created by {@link #createProjectionAroundStart()}, saved
for reuse when
+     * new projection is requested. Only the "Latitude of natural origin" and "Longitude
of natural origin"
+     * parameter values will change for different projections.
+     *
+     * @see #createProjectionAroundStart()
+     */
+    private Parameters projectionParameters;
+
+    /**
+     * Conversion from {@linkplain #getPositionCRS() position CRS} to projection base CRS.
+     * Computed when first needed. This transform does not change after creation.
+     */
+    private MathTransform toProjectionBase;
+
+    /**
+     * The operation method to use for creating a map projection. In the spherical case this
is
+     * "Azimuthal Equidistant (Spherical)". In the ellipsoidal case it become "Modified Azimuthal
Equidistant".
+     */
+    String getProjectionMethod() {
+        return "Azimuthal Equidistant (Spherical)";
+    }
+
+    /**
+     * Creates an <cite>Azimuthal Equidistant</cite> projection centered on current
starting point. On input,
+     * the {@code MathTransform} expects coordinates expressed in the {@linkplain #getPositionCRS()
position CRS}.
+     * On output, the {@code MathTransform} produces coordinates in a {@link org.opengis.referencing.crs.ProjectedCRS}
+     * having the following characteristics:
+     *
+     * <ul>
+     *   <li>Coordinate system is a two-dimensional {@link org.opengis.referencing.cs.CartesianCS}
+     *       with (Easting, Northing) axis order and directions.</li>
+     *   <li>Unit of measurement is the same as {@linkplain #getPositionCRS() position
CRS}
+     *       if those units are linear, or {@link Units#METRE} otherwise.
+     *   <li>Projection of the {@linkplain #getStartPoint() start point} results in
(0,0).</li>
+     *   <li>Distances relative to (0,0) are approximately exact for distances less
than 800 km.</li>
+     *   <li>Azimuths from (0,0) to other points are approximately exact for points
located at less than 800 km.</li>
+     * </ul>
+     *
+     * Given above characteristics, the following calculations are satisfying approximations
when using
+     * (<var>x</var>, <var>y</var>) coordinates in the output space
for <var>D</var> &lt; 800 km:
+     *
+     * <blockquote>
+     * <var>D</var> = √(<var>x</var>² + <var>y</var>²)
 — distance from projection center.<br>
+     * <var>θ</var> = atan2(<var>y</var>, <var>x</var>)
— arithmetic angle from projection center to (<var>x</var>, <var>y</var>).<br>
+     * <var>x</var> = <var>D</var>⋅cos <var>θ</var><br>
+     * <var>y</var> = <var>D</var>⋅sin <var>θ</var>
    — end point for a distance and angle from start point.<br>
+     * </blockquote>
+     *
+     * The following calculations are <strong>not</strong> exacts, because distances
and azimuths are approximately
+     * exacts only when measured from (0,0) coordinates:
+     *
+     * <blockquote>
+     * <var>D</var> = √[(<var>x₂</var> − <var>x₁</var>)²
+ (<var>y₂</var> − <var>y₁</var>)²]
+     *          — distances between points other then projection center are not valid.<br>
+     * <var>θ</var> = atan2(<var>y₂</var> − <var>y₁</var>,
<var>x₂</var> − <var>x₁</var>)
+     *          — azimuths between points other then projection center are not valid.<br>
+     * <i>etc.</i>
+     * </blockquote>
+     *
+     * This method can be invoked repetitively for doing calculations around different points.
+     * All returned {@link MathTransform} instances are immutable;
+     * changing {@code GeodeticCalculator} state does not affect those transforms.
+     *
+     * @return transform from {@linkplain #getPositionCRS() position CRS} to <cite>Azimuthal
Equidistant</cite>
+     *         projected CRS centered on current {@linkplain #getStartPoint() start point}.
+     * @throws IllegalStateException if the start point has not been set.
+     * @throws GeodeticException if the projection can not be computed.
+     *
+     * @see #moveToEndPoint()
+     * @see org.apache.sis.referencing.operation.projection.AzimuthalEquidistant
+     *
+     * @since 1.1
+     */
+    public MathTransform createProjectionAroundStart() {
+        if (isInvalid(START_POINT)) {
+            throw new IllegalStateException(Resources.format(Resources.Keys.StartOrEndPointNotSet_1,
0));
+        }
+        try {
+            if (projectionParameters == null) {
+                final CoordinateReferenceSystem positionCRS, baseCRS;
+                final Unit<?>       crsUnit;
+                final UnitConverter toLinearUnit;
+
+                positionCRS           = getPositionCRS();
+                baseCRS               = ReferencingUtilities.toNormalizedGeographicCRS(positionCRS,
false, false);
+                crsUnit               = ReferencingUtilities.getUnit(positionCRS.getCoordinateSystem());
+                toLinearUnit          = ellipsoid.getAxisUnit().getConverterTo(Units.isLinear(crsUnit)
? crsUnit.asType(Length.class) : Units.METRE);
+                toProjectionBase      = CRS.findOperation(positionCRS, baseCRS, null).getMathTransform();
+                projectionFactory     = DefaultFactories.forBuildin(MathTransformFactory.class,
DefaultMathTransformFactory.class).caching(false);
+                projectionProvider    = (MapProjection) projectionFactory.getOperationMethod(getProjectionMethod());
+                projectionParameters  = Parameters.castOrWrap(projectionProvider.getParameters().createValue());
+                projectionParameters.parameter(Constants.SEMI_MAJOR).setValue(toLinearUnit.convert(ellipsoid.getSemiMajorAxis()));
+                projectionParameters.parameter(Constants.SEMI_MINOR).setValue(toLinearUnit.convert(ellipsoid.getSemiMinorAxis()));
+            }
+            projectionParameters.getOrCreate(LATITUDE_OF_ORIGIN) .setValue(φ1, Units.RADIAN);
+            projectionParameters.getOrCreate(LONGITUDE_OF_ORIGIN).setValue(λ1, Units.RADIAN);
+            return MathTransforms.concatenate(toProjectionBase,
+                    projectionProvider.createMathTransform(projectionFactory, projectionParameters));
+        } catch (FactoryException e) {
+            throw new GeodeticException(e.getMessage(), e);
+        }
+    }
+
+    /**
      * Returns a string representation of start point, end point, azimuths and distance.
      * The text representation is implementation-specific and may change in any future version.
      * Current implementation is like below:
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 a6be0e5..d2850c3 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
@@ -27,6 +27,8 @@ import java.io.IOException;
 import java.io.LineNumberReader;
 import org.opengis.geometry.DirectPosition;
 import org.opengis.referencing.cs.AxisDirection;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.internal.referencing.j2d.ShapeUtilitiesExt;
 import org.apache.sis.internal.referencing.Formulas;
 import org.apache.sis.referencing.crs.HardCodedCRS;
@@ -61,7 +63,7 @@ import static org.apache.sis.internal.metadata.ReferencingServices.AUTHALIC_RADI
  * This base class tests calculator using spherical formulas.
  * Subclass executes the same test but using ellipsoidal formulas.
  *
- * @version 1.0
+ * @version 1.1
  * @since   1.0
  * @module
  */
@@ -584,4 +586,35 @@ public strictfp class GeodeticCalculatorTest extends TestCase {
         assertEquals(4892987, c.getRhumblineLength(), 1);
         assertEquals(90, c.getConstantAzimuth(), STRICT);
     }
+
+    /**
+     * Tests {@link GeodeticCalculator#createProjectionAroundStart()}.
+     * This method tests self-consistency.
+     *
+     * @throws TransformException if an error occurred while projection the test point.
+     */
+    @Test
+    public void test() throws TransformException {
+        final GeodeticCalculator c = create(false);
+        final double distance = 600000;                         // In metres.
+        final double azimuth  = 37;                             // Geographic angle (degrees
relative to North).
+        final double theta    = Math.toRadians(90 - azimuth);   // Azimuth converted to arithmetic
angle.
+        c.setStartGeographicPoint(60, 20);
+        /*
+         * Compute the end point using map projection instead than GeodeticCalculator.
+         * This is valid only for distance within 800 km when using ellipsoidal formulas.
+         */
+        final MathTransform projection = c.createProjectionAroundStart();
+        DirectPosition endPoint = new DirectPosition2D(distance*Math.cos(theta), distance*Math.sin(theta));
+        endPoint = projection.inverse().transform(endPoint, endPoint);
+        /*
+         * Compute the geodesic distance between the point compute by GeodeticCalculator
+         * and the point computed by the Azimuthal Equidistant projection.
+         */
+        c.setStartingAzimuth(azimuth);
+        c.setGeodesicDistance(distance);
+        c.moveToEndPoint();
+        c.setEndPoint(endPoint);
+        assertEquals(0, c.getGeodesicDistance(), 1);
+    }
 }


Mime
View raw message