001/*
002 * Copyright (c) 2003 Objectix Pty Ltd  All rights reserved.
003 *
004 * This library is free software; you can redistribute it and/or
005 * modify it under the terms of the GNU Lesser General Public
006 * License as published by the Free Software Foundation.
007 *
008 * THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESSED OR IMPLIED
009 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
010 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
011 * DISCLAIMED.  IN NO EVENT SHALL OBJECTIX PTY LTD BE LIABLE FOR ANY
012 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
013 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
014 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
015 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
016 * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
017 * OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
018 * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
019 */
020package org.openstreetmap.josm.data.projection.datum;
021
022import java.io.IOException;
023import java.io.InputStream;
024import java.io.Serializable;
025import java.nio.charset.StandardCharsets;
026import java.util.Arrays;
027import java.util.Objects;
028
029import org.openstreetmap.josm.tools.Logging;
030import org.openstreetmap.josm.tools.Utils;
031
032/**
033 * Models the NTv2 Sub Grid within a Grid Shift File.
034 *
035 * @author Peter Yuill
036 * Modified for JOSM :
037 * - removed the RandomAccessFile mode (Pieren)
038 * - read grid file by single bytes. Workaround for a bug in some VM not supporting
039 *   file reading by group of 4 bytes from a jar file.
040 * - removed the Cloneable interface
041 * @since 2507
042 */
043public class NTV2SubGrid implements Serializable {
044
045    private static final long serialVersionUID = 1L;
046
047    private final String subGridName;
048    private final String parentSubGridName;
049    private final String created;
050    private final String updated;
051    private final double minLat;
052    private final double maxLat;
053    private final double minLon;
054    private final double maxLon;
055    private final double latInterval;
056    private final double lonInterval;
057    private final int nodeCount;
058
059    private final int lonColumnCount;
060    private final int latRowCount;
061    private final float[] latShift;
062    private final float[] lonShift;
063    private float[] latAccuracy;
064    private float[] lonAccuracy;
065
066    private NTV2SubGrid[] subGrid;
067
068    /**
069     * Construct a Sub Grid from an InputStream, loading the node data into
070     * arrays in this object.
071     *
072     * @param in GridShiftFile InputStream
073     * @param bigEndian is the file bigEndian?
074     * @param loadAccuracy is the node Accuracy data to be loaded?
075     * @throws IOException if any I/O error occurs
076     */
077    public NTV2SubGrid(InputStream in, boolean bigEndian, boolean loadAccuracy) throws IOException {
078        byte[] b8 = new byte[8];
079        byte[] b4 = new byte[4];
080        byte[] b1 = new byte[1];
081        readBytes(in, b8);
082        readBytes(in, b8);
083        subGridName = new String(b8, StandardCharsets.UTF_8).trim();
084        readBytes(in, b8);
085        readBytes(in, b8);
086        parentSubGridName = new String(b8, StandardCharsets.UTF_8).trim();
087        readBytes(in, b8);
088        readBytes(in, b8);
089        created = new String(b8, StandardCharsets.UTF_8);
090        readBytes(in, b8);
091        readBytes(in, b8);
092        updated = new String(b8, StandardCharsets.UTF_8);
093        readBytes(in, b8);
094        readBytes(in, b8);
095        minLat = NTV2Util.getDouble(b8, bigEndian);
096        readBytes(in, b8);
097        readBytes(in, b8);
098        maxLat = NTV2Util.getDouble(b8, bigEndian);
099        readBytes(in, b8);
100        readBytes(in, b8);
101        minLon = NTV2Util.getDouble(b8, bigEndian);
102        readBytes(in, b8);
103        readBytes(in, b8);
104        maxLon = NTV2Util.getDouble(b8, bigEndian);
105        readBytes(in, b8);
106        readBytes(in, b8);
107        latInterval = NTV2Util.getDouble(b8, bigEndian);
108        readBytes(in, b8);
109        readBytes(in, b8);
110        lonInterval = NTV2Util.getDouble(b8, bigEndian);
111        lonColumnCount = 1 + (int) ((maxLon - minLon) / lonInterval);
112        latRowCount = 1 + (int) ((maxLat - minLat) / latInterval);
113        readBytes(in, b8);
114        readBytes(in, b8);
115        nodeCount = NTV2Util.getInt(b8, bigEndian);
116        if (nodeCount != lonColumnCount * latRowCount)
117            throw new IllegalStateException("SubGrid " + subGridName + " has inconsistent grid dimesions");
118        latShift = new float[nodeCount];
119        lonShift = new float[nodeCount];
120        if (loadAccuracy) {
121            latAccuracy = new float[nodeCount];
122            lonAccuracy = new float[nodeCount];
123        }
124
125        for (int i = 0; i < nodeCount; i++) {
126            // Read the grid file byte after byte. This is a workaround about a bug in
127            // certain VM which are not able to read byte blocks when the resource file is in a .jar file (Pieren)
128            readBytes(in, b1); b4[0] = b1[0];
129            readBytes(in, b1); b4[1] = b1[0];
130            readBytes(in, b1); b4[2] = b1[0];
131            readBytes(in, b1); b4[3] = b1[0];
132            latShift[i] = NTV2Util.getFloat(b4, bigEndian);
133            readBytes(in, b1); b4[0] = b1[0];
134            readBytes(in, b1); b4[1] = b1[0];
135            readBytes(in, b1); b4[2] = b1[0];
136            readBytes(in, b1); b4[3] = b1[0];
137            lonShift[i] = NTV2Util.getFloat(b4, bigEndian);
138            readBytes(in, b1); b4[0] = b1[0];
139            readBytes(in, b1); b4[1] = b1[0];
140            readBytes(in, b1); b4[2] = b1[0];
141            readBytes(in, b1); b4[3] = b1[0];
142            if (loadAccuracy) {
143                latAccuracy[i] = NTV2Util.getFloat(b4, bigEndian);
144            }
145            readBytes(in, b1); b4[0] = b1[0];
146            readBytes(in, b1); b4[1] = b1[0];
147            readBytes(in, b1); b4[2] = b1[0];
148            readBytes(in, b1); b4[3] = b1[0];
149            if (loadAccuracy) {
150                lonAccuracy[i] = NTV2Util.getFloat(b4, bigEndian);
151            }
152        }
153    }
154
155    private static void readBytes(InputStream in, byte[] b) throws IOException {
156        if (in.read(b) < b.length) {
157            Logging.error("Failed to read expected amount of bytes ("+ b.length +") from stream");
158        }
159    }
160
161    /**
162     * Tests if a specified coordinate is within this Sub Grid
163     * or one of its Sub Grids. If the coordinate is outside
164     * this Sub Grid, null is returned. If the coordinate is
165     * within this Sub Grid, but not within any of its Sub Grids,
166     * this Sub Grid is returned. If the coordinate is within
167     * one of this Sub Grid's Sub Grids, the method is called
168     * recursively on the child Sub Grid.
169     *
170     * @param lon Longitude in Positive West Seconds
171     * @param lat Latitude in Seconds
172     * @return the Sub Grid containing the Coordinate or null
173     */
174    public NTV2SubGrid getSubGridForCoord(double lon, double lat) {
175        return !isCoordWithin(lon, lat)
176                ? null
177                : subGrid == null
178                ? this
179                : Arrays.stream(subGrid)
180                .filter(aSubGrid -> aSubGrid.isCoordWithin(lon, lat))
181                .map(aSubGrid -> aSubGrid.getSubGridForCoord(lon, lat))
182                .filter(Objects::nonNull)
183                .findFirst().orElse(this);
184    }
185
186    /**
187     * Tests if a specified coordinate is within this Sub Grid.
188     * A coordinate on either outer edge (maximum Latitude or
189     * maximum Longitude) is deemed to be outside the grid.
190     *
191     * @param lon Longitude in Positive West Seconds
192     * @param lat Latitude in Seconds
193     * @return true or false
194     */
195    private boolean isCoordWithin(double lon, double lat) {
196        return (lon >= minLon) && (lon < maxLon) && (lat >= minLat) && (lat < maxLat);
197    }
198
199    /**
200     * Bi-Linear interpolation of four nearest node values as described in
201     * 'GDAit Software Architecture Manual' produced by the <a
202     * href='http://www.dtpli.vic.gov.au/property-and-land-titles/geodesy/geocentric-datum-of-australia-1994-gda94/gda94-useful-tools'>
203     * Geomatics Department of the University of Melbourne</a>
204     * @param a value at the A node
205     * @param b value at the B node
206     * @param c value at the C node
207     * @param d value at the D node
208     * @param x Longitude factor
209     * @param y Latitude factor
210     * @return interpolated value
211     */
212    private static double interpolate(float a, float b, float c, float d, double x, double y) {
213        return a + (((double) b - (double) a) * x) + (((double) c - (double) a) * y) +
214        (((double) a + (double) d - b - c) * x * y);
215    }
216
217    /**
218     * Interpolate shift and accuracy values for a coordinate in the 'from' datum
219     * of the GridShiftFile. The algorithm is described in
220     * 'GDAit Software Architecture Manual' produced by the <a
221     * href='http://www.dtpli.vic.gov.au/property-and-land-titles/geodesy/geocentric-datum-of-australia-1994-gda94/gda94-useful-tools'>
222     * Geomatics Department of the University of Melbourne</a>
223     * <p>This method is thread safe for both memory based and file based node data.
224     * @param gs GridShift object containing the coordinate to shift and the shift values
225     */
226    public void interpolateGridShift(NTV2GridShift gs) {
227        int lonIndex = (int) ((gs.getLonPositiveWestSeconds() - minLon) / lonInterval);
228        int latIndex = (int) ((gs.getLatSeconds() - minLat) / latInterval);
229
230        double x = (gs.getLonPositiveWestSeconds() - (minLon + (lonInterval * lonIndex))) / lonInterval;
231        double y = (gs.getLatSeconds() - (minLat + (latInterval * latIndex))) / latInterval;
232
233        // Find the nodes at the four corners of the cell
234
235        int indexA = lonIndex + (latIndex * lonColumnCount);
236        int indexB = indexA + 1;
237        int indexC = indexA + lonColumnCount;
238        int indexD = indexC + 1;
239
240        gs.setLonShiftPositiveWestSeconds(interpolate(
241                lonShift[indexA], lonShift[indexB], lonShift[indexC], lonShift[indexD], x, y));
242
243        gs.setLatShiftSeconds(interpolate(
244                latShift[indexA], latShift[indexB], latShift[indexC], latShift[indexD], x, y));
245
246        if (lonAccuracy == null) {
247            gs.setLonAccuracyAvailable(false);
248        } else {
249            gs.setLonAccuracyAvailable(true);
250            gs.setLonAccuracySeconds(interpolate(
251                    lonAccuracy[indexA], lonAccuracy[indexB], lonAccuracy[indexC], lonAccuracy[indexD], x, y));
252        }
253
254        if (latAccuracy == null) {
255            gs.setLatAccuracyAvailable(false);
256        } else {
257            gs.setLatAccuracyAvailable(true);
258            gs.setLatAccuracySeconds(interpolate(
259                    latAccuracy[indexA], latAccuracy[indexB], latAccuracy[indexC], latAccuracy[indexD], x, y));
260        }
261    }
262
263    /**
264     * Returns the parent sub grid name.
265     * @return the parent sub grid name
266     */
267    public String getParentSubGridName() {
268        return parentSubGridName;
269    }
270
271    /**
272     * Returns the sub grid name.
273     * @return the sub grid name
274     */
275    public String getSubGridName() {
276        return subGridName;
277    }
278
279    /**
280     * Returns the node count.
281     * @return the node count
282     */
283    public int getNodeCount() {
284        return nodeCount;
285    }
286
287    /**
288     * Returns the sub grid count.
289     * @return the sub grid count
290     */
291    public int getSubGridCount() {
292        return subGrid == null ? 0 : subGrid.length;
293    }
294
295    /**
296     * Set an array of Sub Grids of this sub grid
297     * @param subGrid subgrids
298     */
299    public void setSubGridArray(NTV2SubGrid... subGrid) {
300        this.subGrid = Utils.copyArray(subGrid);
301    }
302
303    @Override
304    public String toString() {
305        return subGridName;
306    }
307
308    /**
309     * Returns textual details about the sub grid.
310     * @return textual details about the sub grid
311     */
312    public String getDetails() {
313        return new StringBuilder(256)
314            .append("Sub Grid : ")
315            .append(subGridName)
316            .append("\nParent   : ")
317            .append(parentSubGridName)
318            .append("\nCreated  : ")
319            .append(created)
320            .append("\nUpdated  : ")
321            .append(updated)
322            .append("\nMin Lat  : ")
323            .append(minLat)
324            .append("\nMax Lat  : ")
325            .append(maxLat)
326            .append("\nMin Lon  : ")
327            .append(minLon)
328            .append("\nMax Lon  : ")
329            .append(maxLon)
330            .append("\nLat Intvl: ")
331            .append(latInterval)
332            .append("\nLon Intvl: ")
333            .append(lonInterval)
334            .append("\nNode Cnt : ")
335            .append(nodeCount)
336            .toString();
337    }
338
339    /**
340     * Get maximum latitude value
341     * @return maximum latitude
342     */
343    public double getMaxLat() {
344        return maxLat;
345    }
346
347    /**
348     * Get maximum longitude value
349     * @return maximum longitude
350     */
351    public double getMaxLon() {
352        return maxLon;
353    }
354
355    /**
356     * Get minimum latitude value
357     * @return minimum latitude
358     */
359    public double getMinLat() {
360        return minLat;
361    }
362
363    /**
364     * Get minimum longitude value
365     * @return minimum longitude
366     */
367    public double getMinLon() {
368        return minLon;
369    }
370}