diff --git a/core/src/main/java/org/locationtech/proj4j/datum/Grid.java b/core/src/main/java/org/locationtech/proj4j/datum/Grid.java index b000050..fb84969 100755 --- 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 fromNadGrids(String grids) throws IOException { - List gridlist = new ArrayList(); + List 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 new file mode 100755 index 0000000..303152d --- /dev/null +++ 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: + *
+ * 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|
+ * 
+ * + * Subfile header: + *
+ * |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    |
+ * 
+ * + * Grid shift records + *
+ * |dddd|dddd|dddd|dddd|
+ * 
+ * + * End of File record + *
+ * |END     |dddddddd|
+ * 
+ * + * @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 new file mode 100644 index 0000000..de08d52 --- /dev/null +++ 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 new file mode 100644 index 0000000..cc1e29e Binary files /dev/null and b/core/src/test/resources/proj4/nad/100800401.gsb differ