GeocentricConverter equality check after WGS param overrides (#77)

* GeocentricConverter equality check after WGS param overrides * comments, geocentric converter comparisons * newly fixed tests * removed no longer failing test
locationtech · Nov 3, 2021 · 45fa5a1 · 45fa5a1
1 parent 03d3123
commit 45fa5a1
Show file tree

Hide file tree

Showing 8 changed files with 61 additions and 36 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 - Adjustment to OSGB36 datum transform e.g. EPSG: 27700
+- GeocentricConverter equality check after grid shift WGS param override e.g. EPSG: 27700 [#32]
 - +nadgrids=@null support e.g. EPSG: 3857
 
 ## [1.1.3] - 2021-06-17

diff --git a/src/main/java/org/locationtech/proj4j/BasicCoordinateTransform.java b/src/main/java/org/locationtech/proj4j/BasicCoordinateTransform.java
@@ -79,28 +79,45 @@ public BasicCoordinateTransform(CoordinateReferenceSystem srcCRS,
         doDatumTransform = doInverseProjection && doForwardProjection
                 && srcCRS.getDatum() != tgtCRS.getDatum();
 
+        boolean geocentric = false;
+
         if (doDatumTransform) {
 
             boolean isEllipsoidEqual = srcCRS.getDatum().getEllipsoid().isEqual(tgtCRS.getDatum().getEllipsoid());
-            transformViaGeocentric = ! isEllipsoidEqual || srcCRS.getDatum().hasTransformToWGS84()
+            geocentric = ! isEllipsoidEqual || srcCRS.getDatum().hasTransformToWGS84()
                     || tgtCRS.getDatum().hasTransformToWGS84();
 
-            if (transformViaGeocentric) {
+            if (geocentric) {
                 srcGeoConv = new GeocentricConverter(srcCRS.getDatum().getEllipsoid());
                 tgtGeoConv = new GeocentricConverter(tgtCRS.getDatum().getEllipsoid());
 
-                if (srcCRS.getDatum().getTransformType() == Datum.TYPE_GRIDSHIFT) {
-                    srcGeoConv.overrideWithWGS84Params();
-                }
+                int srcTransformType = srcCRS.getDatum().getTransformType();
+                int tgtTransformType = tgtCRS.getDatum().getTransformType();
+
+                if (srcTransformType == Datum.TYPE_GRIDSHIFT || tgtTransformType == Datum.TYPE_GRIDSHIFT) {
+
+	                if (srcTransformType == Datum.TYPE_GRIDSHIFT) {
+	                    srcGeoConv.overrideWithWGS84Params();
+	                }
+
+	                if (tgtTransformType == Datum.TYPE_GRIDSHIFT) {
+	                    tgtGeoConv.overrideWithWGS84Params();
+	                }
+
+	                // After WGS84 params override, check if geocentric transform is still required
+	                // https://github.com/OSGeo/PROJ/blob/5.2.0/src/pj_transform.c#L892
+	                if(srcGeoConv.isEqual(tgtGeoConv)) {
+	                    geocentric = false;
+	                    srcGeoConv = null;
+	                    tgtGeoConv = null;
+	                }
 
-                if (tgtCRS.getDatum().getTransformType() == Datum.TYPE_GRIDSHIFT) {
-                    tgtGeoConv.overrideWithWGS84Params();
                 }
             }
 
-        } else {
-        	transformViaGeocentric=false;
         }
+
+        transformViaGeocentric = geocentric;
     }
 
     @Override
@@ -172,6 +189,10 @@ private void datumTransform(ProjCoordinate pt) {
                 || tgtCRS.getDatum().getTransformType() == Datum.TYPE_UNKNOWN)
             return;
 
+        /* -------------------------------------------------------------------- */
+        /*	If this datum requires grid shifts, then apply it to geodetic    */
+        /*      coordinates.                                                    */
+        /* -------------------------------------------------------------------- */
         if (srcCRS.getDatum().getTransformType() == Datum.TYPE_GRIDSHIFT) {
             srcCRS.getDatum().shift(pt);
         }
@@ -202,7 +223,9 @@ private void datumTransform(ProjCoordinate pt) {
             tgtGeoConv.convertGeocentricToGeodetic(pt);
         }
 
-
+        /* -------------------------------------------------------------------- */
+        /*      Apply grid shift to destination if required.                    */
+        /* -------------------------------------------------------------------- */
         if (tgtCRS.getDatum().getTransformType() == Datum.TYPE_GRIDSHIFT) {
             tgtCRS.getDatum().inverseShift(pt);
         }

diff --git a/src/main/java/org/locationtech/proj4j/datum/GeocentricConverter.java b/src/main/java/org/locationtech/proj4j/datum/GeocentricConverter.java
@@ -60,15 +60,16 @@ public class GeocentricConverter implements java.io.Serializable {
     double ep2;
 
     public GeocentricConverter(Ellipsoid ellipsoid) {
-        this(ellipsoid.getA(), ellipsoid.getB());
+        // Preserve the ellipsoid value precisions
+        this(ellipsoid.getA(), ellipsoid.getB(), ellipsoid.getEccentricitySquared());
     }
 
-    public GeocentricConverter(double a, double b) {
+    public GeocentricConverter(double a, double b, double e2) {
         this.a = a;
         this.b = b;
         a2 = a * a;
         b2 = b * b;
-        e2 = (a2 - b2) / a2;
+        this.e2 = e2;
         ep2 = (a2 - b2) / b2;
     }
 
@@ -77,6 +78,12 @@ public void overrideWithWGS84Params() {
         this.e2 = Ellipsoid.WGS84.getEccentricitySquared();
     }
 
+    public boolean isEqual(GeocentricConverter gc) {
+        // Check if geocentricly equal
+        // https://github.com/OSGeo/PROJ/blob/5.2.0/src/pj_transform.c#L892
+        return this.a == gc.a && this.e2 == gc.e2;
+    }
+
     /**
      * Converts geodetic coordinates
      * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),

diff --git a/src/main/java/org/locationtech/proj4j/datum/Grid.java b/src/main/java/org/locationtech/proj4j/datum/Grid.java
@@ -91,6 +91,9 @@ public static void shift(List<Grid> grids, boolean inverse, ProjCoordinate in) {
 
         for (Grid grid : grids) {
             ConversionTable table = grid.table;
+            // don't shift if the grid is invalid
+            // https://github.com/OSGeo/PROJ/blob/5.2.0/src/pj_gridlist.c#L88
+            if(table == null) continue;
             double epsilon = (Math.abs(table.del.phi) + Math.abs(table.del.lam)) / 10000d;
             // Skip tables that don't match our point at all
             if (table.ll.phi - epsilon > input.phi
@@ -314,7 +317,6 @@ public static List<Grid> fromNadGrids(String grids) throws IOException {
             for (String gridName : grids.split(",")) {
                 boolean optional = gridName.startsWith("@");
                 if (optional) gridName = gridName.substring(1);
-                if (gridName.equals("null")) return null;
                 try {
                     mergeGridFile(gridName, gridlist);
                 } catch (IOException e) {
@@ -331,6 +333,7 @@ private static Grid gridinfoInit(String gridName) throws IOException {
         grid.gridName = gridName;
         grid.format = "missing";
         grid.gridOffset = 0;
+        if (gridName.equals("null")) return grid;
         try(DataInputStream gridDefinition = resolveGridDefinition(gridName)) {
             if (gridDefinition == null) {
                 throw new IOException("Unknown grid: " + gridName);

diff --git a/src/test/java/org/locationtech/proj4j/CoordinateTransformTest.java b/src/test/java/org/locationtech/proj4j/CoordinateTransformTest.java
@@ -42,8 +42,8 @@ public void testFirst() {
   */
     }
 
-//    @Test
-    public void FAIL_testEPSG_27700() {
+    @Test
+    public void testEPSG_27700() {
         checkTransform("EPSG:4326", -2.89, 55.4, "EPSG:27700", 343733.1404, 612144.530677, 0.1);
         checkTransformAndInverse(
                 "EPSG:4326", -2.0301713578021983, 53.35168607080468,
@@ -219,8 +219,14 @@ public void testPROJ4() {
         checkTransformToWGS84("EPSG:27700", 612435.55, 1234954.16, 1.9200000236235546, 60.93999999543101, 0.0);
         checkTransformToWGS84("EPSG:27700", 327420.988668, 690284.547110, -3.1683134533969364, 56.0998025292667, 0.0);
         checkTransformFromWGS84("EPSG:3857", -3.1683134533969364, 56.0998025292667, -352695.04030562507, 7578309.225014557, 0.0);
-        // TODO https://github.com/locationtech/proj4j/issues/32
-        //checkTransform("EPSG:27700", 327420.988668, 690284.547110, "EPSG:3857", -352695.04030562507, 7578309.225014557, 0.0);
+        checkTransform("EPSG:27700", 327420.988668, 690284.547110, "EPSG:3857", -352695.04030562507, 7578309.225014557, 0.0);
+        checkTransform("EPSG:3857", -352695.04030562507, 7578309.225014557, "EPSG:27700", 327420.988668, 690284.547110, 0.001);
+        checkTransform("EPSG:31469", 5439627.33, 5661628.09, "EPSG:3857", 1573657.37, 6636624.41, 0.01);
+        checkTransform("EPSG:3857", 1573657.37, 6636624.41, "EPSG:31469", 5439627.33, 5661628.09, 0.01);
+        checkTransform("EPSG:2056", 2600670.52, 1199667.32, "EPSG:3857", 829045.23, 5933605.15, 0.01);
+        checkTransform("EPSG:3857", 829045.23, 5933605.15, "EPSG:2056", 2600670.52, 1199667.32, 0.01);
+        checkTransform("EPSG:3857", -20037508.342789244, -20037366.780895382, "EPSG:4055", -180.0, -85.01794318500549, 0.001);
+        checkTransform("EPSG:4055", -180.0, -85.01794318500549, "EPSG:3857", -20037508.342789244, -20037366.780895382, 0.0);
     }
 
     @Test

diff --git a/src/test/java/org/locationtech/proj4j/ExampleTest.java b/src/test/java/org/locationtech/proj4j/ExampleTest.java
@@ -223,15 +223,12 @@ public void latLonToStereBidirectionalTransform() {
     @Test
     public void epsgWebMercatorLegacyTest() {
         CRSFactory csFactory = new CRSFactory();
-        String parameters = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs";
         try {
-            String code = csFactory.readEpsgFromParameters(parameters);
+            String code = csFactory.readEpsgFromParameters("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs");
             Assert.assertEquals(Integer.parseInt(code), 3857);
         } catch (IOException e) {
             e.printStackTrace();
         }
-        CoordinateReferenceSystem crs = csFactory.createFromParameters("EPSG:3857", parameters);
-        assertTrue(crs.getDatum().getTransformType() != Datum.TYPE_GRIDSHIFT);
     }
 
     private boolean isInTolerance(ProjCoordinate p, double x, double y, double tolerance) {

diff --git a/src/test/java/org/locationtech/proj4j/Proj4VariousTest.java b/src/test/java/org/locationtech/proj4j/Proj4VariousTest.java
@@ -59,9 +59,8 @@ public void test3ParamToRawSameEllipsoid() {
                 "+proj=latlong +ellps=bessel", p("0dE  0dN 4.000"), 1e-5);
     }
 
-//    @Test
-    public void FAIL_test3ParamToRawSameEllipsoid2() {
-        // fails - not sure why, possibly missing towgs not handled in same way as PROJ4?
+    @Test
+    public void test3ParamToRawSameEllipsoid2() {
         checkTransform(
                 "+proj=latlong +ellps=bessel +towgs84=5,0,0", p("79d00'00.000W 45d00'00.000N 0.0"),
                 "+proj=latlong +ellps=bessel", p("79dW  45dN 0.000"), 1e-5);

diff --git a/src/test/java/org/locationtech/proj4j/TransformFailures.java b/src/test/java/org/locationtech/proj4j/TransformFailures.java
@@ -15,8 +15,6 @@
  *******************************************************************************/
 package org.locationtech.proj4j;
 
-import org.junit.Test;
-
 /**
  * A class to run tests which are known to be failures.
  * This prevents Maven from running them automatically and reporting failures.
@@ -25,13 +23,4 @@
  */
 public class TransformFailures extends BaseCoordinateTransformTest {
 
-    @Test
-    public void testEPSG_27700() {
-        checkTransform("EPSG:4326", -2.89, 55.4, "EPSG:27700", 343733.1404, 612144.530677, 0.1);
-        checkTransformAndInverse(
-                "EPSG:4326", -2.0301713578021983, 53.35168607080468,
-                "EPSG:27700", 398089, 383867,
-                0.001, 0.2 * APPROX_METRE_IN_DEGREES);
-    }
-
 }