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: Document the approximation done regarding CoordinateOperationContext when changing the objective CRS or the Point Of Interest (POI). Replace the use of "Authalic radius" by "Radius of conformal sphere" at the latitude where the spatial resolution is computed.
Date Tue, 11 Feb 2020 23:32:31 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 c060095  Document the approximation done regarding CoordinateOperationContext when
changing the objective CRS or the Point Of Interest (POI). Replace the use of "Authalic radius"
by "Radius of conformal sphere" at the latitude where the spatial resolution is computed.
c060095 is described below

commit c060095a58ce2a34c72172ec6a33606bfdffffa3
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed Feb 12 00:30:45 2020 +0100

    Document the approximation done regarding CoordinateOperationContext when changing the
objective CRS or the Point Of Interest (POI).
    Replace the use of "Authalic radius" by "Radius of conformal sphere" at the latitude where
the spatial resolution is computed.
---
 .../java/org/apache/sis/internal/map/Canvas.java    | 16 +++++++++++++++-
 .../org/apache/sis/internal/map/CanvasContext.java  | 19 +++++++------------
 .../apache/sis/internal/referencing/Formulas.java   | 21 ++++++++++++++++++++-
 .../sis/internal/referencing/FormulasTest.java      | 18 +++++++++++++++++-
 4 files changed, 59 insertions(+), 15 deletions(-)

diff --git a/core/sis-portrayal/src/main/java/org/apache/sis/internal/map/Canvas.java b/core/sis-portrayal/src/main/java/org/apache/sis/internal/map/Canvas.java
index 59b501c..d70de30 100644
--- a/core/sis-portrayal/src/main/java/org/apache/sis/internal/map/Canvas.java
+++ b/core/sis-portrayal/src/main/java/org/apache/sis/internal/map/Canvas.java
@@ -513,6 +513,12 @@ public class Canvas extends Observable implements Localized {
                  * Compute the change unconditionally as a way to verify that the new CRS
is compatible with
                  * data currently shown. Another reason is that checking identity transform
is more reliable
                  * than the `compareIgnoreMetadata(oldValue, newValue)` check.
+                 *
+                 * Note: we are invoking `findTransform(…)` with a CoordinateOperationContext
computed from
+                 * the old CRS. But it is okay because the context information are geographic
area (degrees)
+                 * and approximate resolution (metres), which should not change a lot since
we will continue
+                 * to view the same area after the CRS change. Those information only need
to be approximate
+                 * anyway, and in many cases will be totally ignored.
                  */
                 final MathTransform newToOld = findTransform(newValue, oldValue);
                 if (pointOfInterest != null && !newToOld.isIdentity()) {
@@ -745,7 +751,15 @@ public class Canvas extends Observable implements Localized {
             /*
              * Transform the Point Of Interest to the objective CRS as a way to test its
validity.
              * All canvas fields will be updated only if this operation succeeds.
-             * Note: `oldValue` can not be null if `mt` is non-null.
+             *
+             * Note 1: in the CoordinateOperationContext used for selecting a MathTransform,
the geographic area is
+             * still the same but the spatial resolution could be slightly different because
computed at a new point
+             * of interest. But we can not use the new point of interest now, because we
need the MathTransform for
+             * computing it. However in practice the resolution is often ignored, or does
not vary a lot in regions
+             * where it matter. So we assume it is okay to keep the CoordinateOperationContext
with old resolution
+             * in the following call to `findTransform(…)` or usage of `multidimToObjective`.
+             *
+             * Note 2: `oldValue` can not be null if `multidimToObjective` is non-null.
              */
             MathTransform mt = multidimToObjective;
             if (mt == null || !Utilities.equalsIgnoreMetadata(crs, oldValue.getCoordinateReferenceSystem()))
{
diff --git a/core/sis-portrayal/src/main/java/org/apache/sis/internal/map/CanvasContext.java
b/core/sis-portrayal/src/main/java/org/apache/sis/internal/map/CanvasContext.java
index 5e673bd..59ae628 100644
--- a/core/sis-portrayal/src/main/java/org/apache/sis/internal/map/CanvasContext.java
+++ b/core/sis-portrayal/src/main/java/org/apache/sis/internal/map/CanvasContext.java
@@ -69,12 +69,6 @@ final class CanvasContext extends CoordinateOperationContext {
     private CoordinateOperation objectiveToGeographic;
 
     /**
-     * Authalic radius of the geographic CRS which is the target of {@link #objectiveToGeographic},
in metres.
-     * This is computed in same time than {@link #objectiveToGeographic}.
-     */
-    private double radius;
-
-    /**
      * The geographic area, computed when first requested and saved for reuse. This is reset
to {@code null} every time
      * that {@link Canvas#objectiveCRS}, {@link Canvas#objectiveToDisplay} or {@link Canvas#displayBounds}
is modified.
      *
@@ -103,9 +97,7 @@ final class CanvasContext extends CoordinateOperationContext {
      * Sets the operation from {@link Canvas#objectiveCRS} to geographic CRS.
      */
     final void setObjectiveToGeographic(final CoordinateOperation op) {
-        final Ellipsoid ellipsoid = ((GeographicCRS) op.getTargetCRS()).getDatum().getEllipsoid();
-        radius = ellipsoid.getAxisUnit().getConverterTo(Units.METRE).convert(Formulas.getAuthalicRadius(ellipsoid));
-        objectiveToGeographic = op;             // Set only if above line succeeded.
+        objectiveToGeographic = op;
         clear();
     }
 
@@ -192,18 +184,21 @@ final class CanvasContext extends CoordinateOperationContext {
             final Matrix derivative = MathTransforms.derivativeAndTransform(displayToGeographic,
poi, 0, poi, 1);
             final double[] vector   = new double[derivative.getNumCol()];
             final double[] combined = new double[derivative.getNumRow()];
-            for (int j=0; j<combined.length; j++) {
+            for (int j=combined.length; --j>=0;) {                      // Process latitude
(1) before longitude (0).
                 for (int i=0; i<vector.length; i++) {
                     vector[i] = derivative.getElement(j,i);
                 }
                 double m = MathFunctions.magnitude(vector);
                 switch (j) {
-                    case 0: m *= Math.cos(Math.toRadians(poi[1]));      // Adjust longitude,
then fall through.
-                    case 1: m  = Math.toRadians(m) * radius; break;     // Latitude (this
case) or Longitude (case 0).
+                    case 0: m *= Math.cos(combined[1]);                 // Adjust longitude,
then fall through.
+                    case 1: m  = Math.toRadians(m); break;              // Latitude (this
case) or longitude (case 0).
                     // Other cases: assume value already in metres.
                 }
                 combined[j] = m;
             }
+            final Ellipsoid ellipsoid = ((GeographicCRS) objectiveToGeographic.getTargetCRS()).getDatum().getEllipsoid();
+            double radius = Formulas.radiusOfConformalSphere(ellipsoid, combined[1]);
+            radius = ellipsoid.getAxisUnit().getConverterTo(Units.METRE).convert(radius);
             resolution = MathFunctions.magnitude(combined) * radius;
         }
     }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
index 6cce1df..c1d38d2 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
@@ -33,7 +33,7 @@ import static org.apache.sis.internal.metadata.ReferencingServices.NAUTICAL_MILE
  * do not want to expose publicly those arbitrary values (or at least not in a too direct
way).
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  * @since   0.4
  * @module
  */
@@ -161,6 +161,25 @@ public final class Formulas extends Static {
     }
 
     /**
+     * Returns the radius of the conformal sphere (assuming a semi-major axis length of 1)
at a given latitude.
+     * The radius of conformal sphere is computed as below:
+     *
+     * <blockquote>Rc = √(1 – ℯ²) / (1 – ℯ²sin²φ)  where  ℯ² = 1 -
(b/a)²</blockquote>
+     *
+     * This is a function of latitude and therefore not constant.
+     *
+     * @param  ellipsoid  the ellipsoid for which to compute the radius of conformal sphere.
+     * @param  φ          the latitude in radians where to compute the radius of conformal
sphere.
+     * @return radius of the conformal sphere at latitude φ.
+     */
+    public static double radiusOfConformalSphere(final Ellipsoid ellipsoid, final double
φ) {
+        final double sinφ = Math.sin(φ);
+        final double a = ellipsoid.getSemiMajorAxis();
+        final double r = ellipsoid.getSemiMinorAxis() / a;
+        return a * (r / (1 - (1 - r*r) * (sinφ*sinφ)));
+    }
+
+    /**
      * Computes the semi-minor axis length from the given semi-major axis and inverse flattening
factor.
      *
      * @param  semiMajorAxis      the semi-major axis length.
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/FormulasTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/FormulasTest.java
index 46b81ff..e237feb 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/FormulasTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/FormulasTest.java
@@ -16,8 +16,11 @@
  */
 package org.apache.sis.internal.referencing;
 
+import java.util.Collections;
 import org.apache.sis.internal.metadata.ReferencingServices;
 import org.apache.sis.measure.Longitude;
+import org.apache.sis.measure.Units;
+import org.apache.sis.referencing.datum.DefaultEllipsoid;
 import org.apache.sis.referencing.datum.HardCodedDatum;
 import org.apache.sis.test.TestCase;
 import org.junit.Test;
@@ -29,7 +32,7 @@ import static org.junit.Assert.*;
  * Tests {@link Formulas}.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  * @since   0.4
  * @module
  */
@@ -106,4 +109,17 @@ public final strictfp class FormulasTest extends TestCase {
         assertEquals("International 1924", 297, Formulas.getInverseFlattening(6378388, 6356911.9461279465),
1E-11);
         assertEquals("Clarke 1858", 294.26067636926103, Formulas.getInverseFlattening(20926348,
20855233), 1E-11);
     }
+
+    /**
+     * Tests {@link Formulas#radiusOfConformalSphere(Ellipsoid, double)}.
+     * This test computes the Radius of Conformal Sphere using the values given
+     * by the <a href="http://www.iogp.org/pubs/373-07-2.pdf">EPSG guide</a>
for
+     * the <cite>Amersfoort / RD New</cite> projection (a Stereographic one).
+     */
+    @Test
+    public void testRadiusOfConformalSphere() {
+        final DefaultEllipsoid ellipsoid = DefaultEllipsoid.createFlattenedSphere(
+                Collections.singletonMap(DefaultEllipsoid.NAME_KEY, "Bessel 1841"), 6377397.155,
299.1528128, Units.METRE);
+        assertEquals(6382644.571, Formulas.radiusOfConformalSphere(ellipsoid, StrictMath.toRadians(52.156160556)),
0.001);
+    }
 }


Mime
View raw message