GH-89: initial support for NTv2 (#98)

* GH-89: initial support for NTv2 Signed-off-by:Bart Hanssens <bart.hanssens@bosa.fgov.be> * GH-89: formatting and use JDK8 method for skipping bytes Signed-off-by:Bart Hanssens <bart.hanssens@bosa.fgov.be> * GH-89: also support big-endian files Signed-off-by:Bart Hanssens <bart.hanssens@bosa.fgov.be> * GH-89: added test cases Signed-off-by:Bart Hanssens <bart.hanssens@bosa.fgov.be>
locationtech · Mar 28, 2023 · 96e395f · 96e395f
1 parent 069dc85
commit 96e395f
Show file tree

Hide file tree

Showing 4 changed files with 303 additions and 16 deletions.
diff --git a/core/src/main/java/org/locationtech/proj4j/datum/Grid.java b/core/src/main/java/org/locationtech/proj4j/datum/Grid.java
@@ -31,10 +31,6 @@
 import org.locationtech.proj4j.util.PolarCoordinate;
 import org.locationtech.proj4j.util.ProjectionMath;
 
-import java.net.URI;
-import java.net.URISyntaxException;
-import java.net.URL;
-
 /**
  * A Grid represents a geodetic datum defining some mapping between a
  * coordinate system referenced to the surface of the earth and spherical
@@ -151,11 +147,11 @@ public static final class ConversionTable implements Serializable {
          */
         public String id;
         /**
-         * Lower left corner coordinates
+         * Cell size
          */
         public PolarCoordinate del;
         /**
-         * Cell size
+         * Lower left corner coordinates
          */
         public PolarCoordinate ll;
         /**
@@ -317,7 +313,7 @@ private static PolarCoordinate nad_intr(PolarCoordinate t, ConversionTable table
 
     // This method corresponds to the pj_gridlist_from_nadgrids function in proj.4
     public static List<Grid> fromNadGrids(String grids) throws IOException {
-        List<Grid> gridlist = new ArrayList<Grid>();
+        List<Grid> gridlist = new ArrayList<>();
         synchronized (Grid.class) {
             for (String gridName : grids.split(",")) {
                 boolean optional = gridName.startsWith("@");
@@ -338,7 +334,11 @@ private static Grid gridinfoInit(String gridName) throws IOException {
         grid.gridName = gridName;
         grid.format = "missing";
         grid.gridOffset = 0;
-        if (gridName.equals("null")) return grid;
+
+        if (gridName.equals("null")) {
+            return grid;
+        }
+
         try(DataInputStream gridDefinition = resolveGridDefinition(gridName)) {
             if (gridDefinition == null) {
                 throw new IOException("Unknown grid: " + gridName);
@@ -353,15 +353,20 @@ private static Grid gridinfoInit(String gridName) throws IOException {
                 grid.table = CTABLEV2.init(gridDefinition);
                 gridDefinition.reset();
                 CTABLEV2.load(gridDefinition, grid);
+            } else if (NTV1.testHeader(header)) {
+                    grid.format = "ntv1";
+                    gridDefinition.mark(1024);
+                    grid.table = NTV1.init(gridDefinition);
+                    gridDefinition.reset();
+                    NTV1.load(gridDefinition, grid);
+            } else if (NTV2.testHeader(header)) {
+                    grid.format = "ntv2";
+                    gridDefinition.mark(1024);
+                    grid.table = NTV2.init(gridDefinition);
+                    gridDefinition.reset();
+                    NTV2.load(gridDefinition, grid);
             }
-            if (NTV1.testHeader(header)) {
-                grid.format = "ntv1";
-                gridDefinition.mark(1024);
-                grid.table = NTV1.init(gridDefinition);
-                gridDefinition.reset();
-                NTV1.load(gridDefinition, grid);
-            }
-		}
+        }
         return grid;
     }
 

diff --git a/core/src/main/java/org/locationtech/proj4j/datum/NTV2.java b/core/src/main/java/org/locationtech/proj4j/datum/NTV2.java
@@ -0,0 +1,202 @@
+/*******************************************************************************
+ * Copyright 2023 FPS BOSA
+ *
+ * Licensed 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.locationtech.proj4j.datum;
+
+import java.io.DataInputStream;
+import java.io.IOException;
+import java.nio.ByteBuffer;
+import java.nio.ByteOrder;
+import java.nio.charset.StandardCharsets;
+import java.util.Arrays;
+
+import org.locationtech.proj4j.util.FloatPolarCoordinate;
+import org.locationtech.proj4j.util.IntPolarCoordinate;
+import org.locationtech.proj4j.util.PolarCoordinate;
+
+/**
+ * Parser for the "National Transformation" v2 format
+ *
+ * Very basic implementation: only files with 1 subfile are supported, gridshift type is supposed to be in
+ * seconds
+ *
+ * Header structure:
+ * <pre>
+ * 0        8        16
+ * |NUM_OREC|iiii    |
+ * |NUM_SREC|iiii    |
+ * |NUM_FILE|iiii    |
+ * |GS_TYPE |ssssssss|
+ * |VERSION |ssssssss|
+ * |SYSTEM_F|ssssssss|
+ * |SYSTEM_T|ssssssss|
+ * |MAJOR_F |dddddddd|
+ * |MINOR_F |dddddddd|
+ * |MAJOR_T |dddddddd|
+ * |MINOR_T |dddddddd|
+ * </pre>
+ *
+ * Subfile header:
+ * <pre>
+ * |SUB_NAME|ssssssss|
+ * |PARENT  |ssssssss|
+ * |CREATED |ssssssss|
+ * |UPDATED |ssssssss|
+ * |S_LAT   |dddddddd|
+ * |N_LAT   |dddddddd|
+ * |E_LONG  |dddddddd|
+ * |W_LONG  |dddddddd|
+ * |LAT_INC |dddddddd|
+ * |LONG_INC|dddddddd|
+ * |GS_COUNT|iiii    |
+ * </pre>
+ *
+ * Grid shift records
+ * <pre>
+ * |dddd|dddd|dddd|dddd|
+ * </pre>
+ *
+ * End of File record
+ * <pre>
+ * |END     |dddddddd|
+ * </pre>
+ *
+ * @author Bart Hanssens
+ */
+public final class NTV2 {
+
+    private static final byte[] MAGIC = "NUM_OREC".getBytes(StandardCharsets.US_ASCII);
+
+    private static final double SEC_RAD = Math.PI / 180 / 3600;
+
+    private static final int HEADER_SIZE = 176;
+    private static final int SUB_HEADER_SIZE = 176;
+    private static final int VALUES_PER_CELL = 4;
+
+    private static final int NUM_OREC = 8;
+    private static final int S_LAT = 72;
+    private static final int N_LAT = 88;
+    private static final int E_LONG = 104;
+    private static final int W_LONG = 120;
+
+    private static final int LAT_INC = 136;
+    private static final int LONG_INC = 152;
+
+    /**
+     * Use header to check file type
+     *
+     * @param header
+     * @return true if format is NTv2
+     */
+    public static boolean testHeader(byte[] header) {
+        byte[] start = Arrays.copyOfRange(header, 0, MAGIC.length);
+        return Arrays.equals(start, MAGIC);
+    }
+
+    /**
+     * Initialize conversion table
+     *
+     * @param instream
+     * @return
+     * @throws IOException
+     */
+    public static Grid.ConversionTable init(DataInputStream instream) throws IOException {
+        byte[] buf = new byte[HEADER_SIZE];
+        instream.readFully(buf);
+
+        if (!testHeader(buf)) {
+            throw new Error("Not a NTv2 file");
+        }
+        ByteOrder endian = guessByteOrder(buf);
+
+        buf = new byte[SUB_HEADER_SIZE];
+        instream.readFully(buf);
+
+        Grid.ConversionTable table = new Grid.ConversionTable();
+        table.id = "NTv2 Grid Shift File";
+        // lower left
+        table.ll = new PolarCoordinate(-doubleFromBytes(buf, W_LONG, endian) * SEC_RAD, 
+                                        doubleFromBytes(buf, S_LAT, endian) * SEC_RAD);
+        // upper right
+        PolarCoordinate ur = new PolarCoordinate(-doubleFromBytes(buf, E_LONG, endian) * SEC_RAD, 
+                                                doubleFromBytes(buf, N_LAT, endian) * SEC_RAD);
+        // "creative" way to store a pair of values
+        table.del = new PolarCoordinate(doubleFromBytes(buf, LONG_INC, endian) * SEC_RAD, 
+                                        doubleFromBytes(buf, LAT_INC, endian) * SEC_RAD);
+        table.lim = new IntPolarCoordinate(
+            (int) (Math.abs(ur.lam - table.ll.lam) / table.del.lam + 0.5) + 1,
+            (int) (Math.abs(ur.phi - table.ll.phi) / table.del.phi + 0.5) + 1);
+
+        return table;
+    }
+
+    /**
+     * Load grid(sub)file into grid
+     *
+     * @param instream
+     * @param grid
+     * @throws IOException
+     */
+    public static void load(DataInputStream instream, Grid grid) throws IOException {
+        int cols = grid.table.lim.lam;
+        int rows = grid.table.lim.phi;
+
+        byte[] buf = new byte[HEADER_SIZE];
+        instream.readFully(buf);
+        ByteOrder endian = guessByteOrder(buf);
+
+        instream.skipBytes(SUB_HEADER_SIZE);
+
+        FloatPolarCoordinate[] tmp_cvs = new FloatPolarCoordinate[cols * rows];
+
+        float[] row_buff = new float[cols * VALUES_PER_CELL];
+        byte[] byteBuff = new byte[row_buff.length * Float.BYTES];
+
+        for (int row = 0; row < rows; row++) {
+            instream.readFully(byteBuff);
+            ByteBuffer.wrap(byteBuff).order(endian).asFloatBuffer().get(row_buff);
+            for (int col = 0; col < cols; col++) {
+                // only process the shift values, ignoring accuracy values
+                tmp_cvs[row * cols + (cols - col - 1)] = new FloatPolarCoordinate(
+                    (float) (row_buff[VALUES_PER_CELL * col + 1] * SEC_RAD),
+                    (float) (row_buff[VALUES_PER_CELL * col] * SEC_RAD));
+            }
+        }
+        grid.table.cvs = tmp_cvs;
+    }
+
+    /**
+     * Guess byte order / endianness by checking first bytes in header
+     * 
+     * @param header
+     * @return endianness
+     */
+    private static ByteOrder guessByteOrder(byte[] header) {
+        ByteBuffer buffer = ByteBuffer.wrap(header, NUM_OREC, Integer.BYTES);
+        if (buffer.order(ByteOrder.BIG_ENDIAN).getInt() == 11) {
+            return ByteOrder.BIG_ENDIAN;
+        }
+        buffer = ByteBuffer.wrap(header, NUM_OREC, Integer.BYTES);
+        if (buffer.order(ByteOrder.LITTLE_ENDIAN).getInt() == 11) {
+            return ByteOrder.LITTLE_ENDIAN;
+        }
+        throw new Error("Could not determine endianness");
+    }
+
+    private static double doubleFromBytes(byte[] b, int offset, ByteOrder order) {
+        return ByteBuffer.wrap(b, offset, Double.BYTES).order(order).getDouble();
+    }
+
+}
diff --git a/core/src/test/java/org/locationtech/proj4j/datum/NTV2Test.java b/core/src/test/java/org/locationtech/proj4j/datum/NTV2Test.java
@@ -0,0 +1,80 @@
+/*******************************************************************************
+ * Copyright 2023 FPS BOSA
+ *
+ * Licensed 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.locationtech.proj4j.datum;
+
+import org.junit.Assert;
+import org.junit.Test;
+
+import org.locationtech.proj4j.CRSFactory;
+import org.locationtech.proj4j.CoordinateReferenceSystem;
+import org.locationtech.proj4j.CoordinateTransform;
+import org.locationtech.proj4j.CoordinateTransformFactory;
+import org.locationtech.proj4j.ProjCoordinate;
+
+/**
+ * Using grid shifts for Catalonia
+ * @see https://geoinquiets.cat/
+ * 
+ * @author Bart Hanssens
+ */
+public class NTV2Test {
+
+    private static CRSFactory CRS = new CRSFactory();
+    private static CoordinateTransformFactory CT = new CoordinateTransformFactory();
+
+    private static CoordinateReferenceSystem cs1 = CRS.createFromParameters("23031",
+            "+proj=utm +zone=31 +ellps=intl +nadgrids=100800401.gsb +units=m +no_defs");
+    private static CoordinateReferenceSystem cs2 = CRS.createFromName("EPSG:25831");
+
+
+    @Test
+    public void gridShiftNTV2() {
+        CoordinateTransform ct = CT.createTransform(cs1, cs2);
+
+        ProjCoordinate expected1 = new ProjCoordinate(299905.060, 4499796.515);
+        ProjCoordinate result1 = new ProjCoordinate();
+        ct.transform(new ProjCoordinate(300000.0, 4500000.0), result1);
+
+        Assert.assertTrue(expected1.areXOrdinatesEqual(result1, 0.001) &&
+                          expected1.areYOrdinatesEqual(result1, 0.001));
+
+        ProjCoordinate expected2 = new ProjCoordinate(519906.767, 4679795.125);
+        ProjCoordinate result2 = new ProjCoordinate();
+        ct.transform(new ProjCoordinate(520000.0, 4680000.0), result2);
+
+        Assert.assertTrue(expected2.areXOrdinatesEqual(result2, 0.001) &&
+                          expected2.areYOrdinatesEqual(result2, 0.001));
+    }
+
+    @Test
+    public void gridShiftNTV2Inverse() {
+        CoordinateTransform ct = CT.createTransform(cs2, cs1);
+
+        ProjCoordinate expected1 = new ProjCoordinate(315093.094, 4740203.227);
+        ProjCoordinate result1 = new ProjCoordinate();
+        ct.transform(new ProjCoordinate(315000.0, 4740000.0), result1);
+
+        Assert.assertTrue(expected1.areXOrdinatesEqual(result1, 0.001) &&
+                          expected1.areYOrdinatesEqual(result1, 0.001));
+
+        ProjCoordinate expected2 = new ProjCoordinate(420093.993, 4600204.241);
+        ProjCoordinate result2 = new ProjCoordinate();
+        ct.transform(new ProjCoordinate(420000.0, 4600000.0), result2);
+
+        Assert.assertTrue(expected2.areXOrdinatesEqual(result2, 0.001) &&
+                          expected2.areYOrdinatesEqual(result2, 0.001));
+    }
+}
diff --git a/core/src/test/resources/proj4/nad/100800401.gsb b/core/src/test/resources/proj4/nad/100800401.gsb