NADCON changes for revisions 7333:7843

/*
 *    Geotools 2 - OpenSource mapping toolkit
 *    http://www.geotools.org
 *    (C) 2004 Geotools Project Managment Committee (PMC)
 *
 *    This library is free software; you can redistribute it and/or
 *    modify it under the terms of the GNU Lesser General Public
 *    License as published by the Free Software Foundation; either
 *    version 2.1 of the License, or (at your option) any later version.
 *
 *    This library is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *    Lesser General Public License for more details.
 *
 *    You should have received a copy of the GNU Lesser General Public
 *    License along with this library; if not, write to the Free Software
 *    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 *
 *    This package contains documentation from OpenGIS specifications.
 *    OpenGIS consortium's work is fully acknowledged here.
 */
package org.geotools.referencing.operation.transform;

// Geotools dependencies
import org.geotools.metadata.citation.Citation;
import org.geotools.referencing.Identifier;
import org.geotools.referencing.operation.MathTransformProvider;
import org.geotools.resources.Arguments;
import org.geotools.resources.cts.ResourceKeys;
import org.geotools.resources.cts.Resources;
import org.opengis.parameter.GeneralParameterValue;
import org.opengis.parameter.ParameterDescriptor;
import org.opengis.parameter.ParameterDescriptorGroup;
import org.opengis.parameter.ParameterNotFoundException;

// OpenGIS dependencies
import org.opengis.parameter.ParameterValue;
import org.opengis.parameter.ParameterValueGroup;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;

// J2SE dependencies
import java.io.BufferedReader;
import java.io.EOFException;
import java.io.File;
import java.io.FileInputStream;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.io.ObjectInputStream;
import java.io.PrintWriter;
import java.io.Serializable;
import java.net.URL;
import java.nio.ByteBuffer;
import java.nio.ByteOrder;
import java.nio.channels.Channels;
import java.nio.channels.ReadableByteChannel;
import java.util.StringTokenizer;
import java.util.prefs.Preferences;


/**
 * The North American Datum Conversion (NADCON) Transform (EPSG code 9613) is a
 * two dimentional datum shift method, created by the National Geodetic Survey
 * (NGS), that uses interpolated values from two grid shift files. This
 * method is used  to transform NAD27 (EPSG code 4267) datum coordinates (and
 * some others) to  NAD83 (EPSG code 4267) within the United States. There are
 * two sent of grid  shift files: NADCON and High Accuracy Reference Networks
 * (HARN).  NADCON shfts from NAD27 (and some others) to NAD83 while HARN
 * shifts from  the NADCON NAD83 to an improved NAD83. Both sets of grid shift
 * files may be  downloaded from  <a href="http://www.ngs.noaa.gov/PC_PROD/NADCON/">www.ngs.noaa.gov/PC_PROD/NADCON/</a>.
 * <br><br>
 *
 * Some of the NADCON grids, their areas of use, and source datums are shown
 * in the following table.
 * <br><br>
 *
 * <table>
 *   <tr><th>Shift File Name</td><th>Area</td><th>Source Datum</td><th>Accuracy at 67% confidence (m)</td></tr>
 *   <tr><td>CONUS</td><td>Conterminous U S (lower 48 states)</td><td>NAD27</td><td>0.15</td></tr>
 *   <tr><td>ALASKA</td><td>Alaska, incl. Aleutian Islands</td><td>NAD27</td><td>0.5</td></tr>
 *   <tr><td>HAWAII</td><td>Hawaiian Islands</td><td>Old Hawaiian (4135)</td><td>0.2</td></tr>
 *   <tr><td>STLRNC</td><td>St. Lawrence Is., AK</td><td>St. Lawrence Island (4136)</td><td>--</td></tr>
 *   <tr><td>STPAUL </td><td>St. Paul Is., AK</td><td>St. Paul Island (4137)</td><td>--</td></tr>
 *   <tr><td>STGEORGE</td><td>St. George Is., AK</td><td>St. George Island (4138)</td><td>--</td></tr>
 *   <tr><td>PRVI</td><td>Puerto Rico and the Virgin Islands</td><td>Puerto Rico (4139)</td><td>0.05</td></tr>
 * </table>
 * <br><br>
 *
 * Grid shift files come in two formats: binary and text. The files from  the
 * NGS are binary and have .las (latitude shift) and .los (longitude shift)
 * extentions. Text grids may be created with the NGS nadgrd program and have
 * .laa (latitude shift) and .loa (longitude shift) file extentions. Both
 * types of  files may be used here.
 * <br><br>
 *
 * The grid names to use for transforming are parameters of this
 * MathTransform.  This parameter may be the full name and path to the grids
 * or just the name  of the grids if the default location of the grids was set
 * as a preference.  This preference may be set with the main method of this
 * class. <br><br>
 *
 * Transformations here have been tested to be within 0.00001 seconds of
 * values  given by the NGS ndcon210 program for NADCON grids. American Samoa
 * and HARN  shifts have not yet been tested.  <strong>References:</strong>
 *
 * <ul>
 *   <li><a href="http://www.ngs.noaa.gov/PC_PROD/NADCON/Readme.htm">NADCONreadme</a></li>
 *   <li>American Samoa Grids for NADCON - Samoa_Readme.txt</li>
 *   <li><a href="http://www.ngs.noaa.gov/PUBS_LIB/NGS50.pdf">NADCON - The
 *       Application of Minimum-Curvature-Derived  Surfaces in the Transformation of
 *       Positional Data From the North American  Datum of 1927 to the North
 *       American Datum of 1983</a> - NOAA TM.</li>
 *   <li>ndcon210.for - NGS fortran source code for NADCON conversions. See the
 *       following subroutines: TRANSF, TO83, FGRID, INTRP, COEFF and SURF</li>
 *   <li>nadgrd.for - NGS fortran source code to export/import binary and text grid
 *       formats</li>
 *   <li>EPSG Geodesy Parameters database version 6.5</li>
 * </ul>
 *
 * @see <a href="http://www.ngs.noaa.gov/TOOLS/Nadcon/Nadcon.html"> NADCON -
 *      North American Datum Conversion Utility</a>
 *
 * @version $Id$
 * @author Rueben Schulz
 *
 * @TODO the transform code does not deal with the case where grids
 *       cross +- 180 degrees
 */
public class NADCONTransform extends AbstractMathTransform implements Serializable {

    /** Serial number for interoperability with different versions. */
    private static final long serialVersionUID = -4707304160205218546L;

    /**
     * Preference node for the grid shift file location, and its default value.
     */
    private static final String GRID_LOCATION = "Grid location";

    /**
     * Preference node for the grid shift file location, and its default value.
     */
    private static final String DEFAULT_GRID_LOCATION = ".";

    /**
     * Difference allowed in iterative computations. This is half the value
     * used in the NGS fortran code (so all tests pass).
     */
    private static final double TOL = 5.0E-10;

    /** Maximum number of itterations for iterative computations. */
    private static final int MAX_ITER = 10;

    /** Conversion factor from seconds to degrees. */
    private static final double SEC_2_DEG = 3600.0;

    /** latitude and longitude grid shift file names. Output in WKT. */
    private final String latGridName;

    /** latitude and longitude grid shift file names. Output in WKT. */
    private final String longGridName;

    /*
     * the minimum longitude value covered by this grid (degrees)
     */
    private double xmin;

    /*
     * the minimum latiude value covered by this grid (degrees)
     */
    private double ymin;

    /*
     * the maximum longitude value covered by this grid (degrees)
     */
    private double xmax;

    /*
     * the maximum latitude value covered by this grid (degrees)
     */
    private double ymax;

    /*
     * the difference between longitude grid points (degrees)
     */
    private double dx;

    /*
     * the difference between latitude grid points (degrees)
     */
    private double dy;

    /**
     * longitude and latitude grid shift values. Values are organized from  low
     * to high longitude (low x index to high) and low to high latitude (low y
     * index to high).
     */
    private LocalizationGrid gridShift = null;

    /**
     * The {@link gridShift} values as a {@link LocalizationGridTransform2D}.
     * Used for interpolating shift values.
     */
    private MathTransform gridShiftTransform = null;

    /** The inverse of this transform. Will be created only when needed. */
    private transient MathTransform inverse;

    /**
     * Constructs a NADCONTransform from the specified grid shift files.
     *
     * @param latGridName path and name (or just name if {@link GRID_LOCATION}
     *        is set) to the latitude differenc file. This will have a .las or
     *        .laa file extention.
     * @param longGridName path and name (or just name if {@link GRID_LOCATION}
     *        is set) to the longitude differenc file. This will have a .los
     *        or .loa file extention.
     *
     * @throws ParameterNotFoundException if a math transform parameter cannot
     *         be found.
     * @throws FactoryException if there is a problem creating this math transform
     *         (ie file extentions are unknown or there is an error reading the
     *          grid files)
     */
    public NADCONTransform(String latGridName, String longGridName)
            throws ParameterNotFoundException, FactoryException
    {
        this.latGridName = latGridName;
        this.longGridName = longGridName;

        //decide if text or binary grid will be used
        try {
            final URL latGridURL = makeURL(latGridName);
            final URL longGridURL = makeURL(longGridName);

            if ((latGridName.endsWith(".las") && longGridName.endsWith(".los"))
                    || (latGridName.endsWith(".LAS") && longGridName.endsWith(".LOS"))) {
                loadBinaryGrid(latGridURL, longGridURL);
            } else if ((latGridName.endsWith(".laa") && longGridName.endsWith(".loa"))
                    || (latGridName.endsWith(".LAA") && longGridName.endsWith(".LOA"))) {
                loadTextGrid(latGridURL, longGridURL);
            } else {
                throw new FactoryException(Resources.format(
                        ResourceKeys.ERROR_UNSUPPORTED_FILE_TYPE_$2,
                        latGridName.substring(latGridName.lastIndexOf('.')),
                        longGridName.substring(longGridName.lastIndexOf('.'))));
            }

            gridShiftTransform = gridShift.getMathTransform();
        } catch (IOException exception) {
            final Throwable cause = exception.getCause();
            if (cause instanceof FactoryException) {
                throw (FactoryException) cause;
            }
            throw new FactoryException(exception.getLocalizedMessage(),
                exception);
        }
    }

    /**
     * Returns the parameter values for this math transform.
     *
     * @return A copy of the parameter values for this math transform.
     */
    public ParameterValueGroup getParameterValues() {
        final ParameterValue lat_diff_file = new org.geotools.parameter.Parameter(Provider.LAT_DIFF_FILE);
        lat_diff_file.setValue(latGridName);

        final ParameterValue long_diff_file = new org.geotools.parameter.Parameter(Provider.LONG_DIFF_FILE);
        long_diff_file.setValue(longGridName);

        return new org.geotools.parameter.ParameterGroup(
            Provider.PARAMETERS,
            new GeneralParameterValue[] { lat_diff_file, long_diff_file }
        );
    }

    /**
     * Gets the dimension of input points (always 2).
     *
     * @return the source dimention.
     */
    public int getDimSource() {
        return 2;
    }

    /**
     * Gets the dimension of output points (always 2).
     *
     * @return the target dimention.
     */
    public int getDimTarget() {
        return 2;
    }

    /**
     * Returns a URL from the string representation. If the string has  no
     * path, the default path preferece is added.
     *
     * @param str a string representation of a URL
     * @return a URL created from the string representation
     * @throws java.net.MalformedURLException if the URL cannot be created
     */
    private URL makeURL(String str) throws java.net.MalformedURLException {
        //has '/' or '\' or ':', so probably full path to file
        if ((str.indexOf("\\") != -1) || (str.indexOf("/") != -1)
                || (str.indexOf(":") != -1)) {
            return makeURLfromString(str);
        } else {
            //just a file name , prepend base location
            final Preferences prefs = Preferences.systemNodeForPackage(NADCONTransform.class);
            final String baseLocation = prefs.get(GRID_LOCATION,
                    DEFAULT_GRID_LOCATION);
            return makeURLfromString(baseLocation + "/" + str);
        }
    }

    /**
     * Returns a URL based on a string representation. If no protocol is given,
     * it is assumed to be a local file.
     *
     * @param str a string representation of a URL
     * @return a URL created from the string representation
     * @throws java.net.MalformedURLException if the URL cannot be created
     */
    private URL makeURLfromString(String str)
            throws java.net.MalformedURLException
    {
        try {
            return new URL(str);
        } catch (java.net.MalformedURLException e) {
            //try making this with a file protocal
            return new URL("file", "", str);
        }
    }

    /**
     * Reads latitude and longitude binary grid shift file data into {@link
     * grid}.  The file is organized into records, with the first record
     * containing the  header information, followed by the shift data. The
     * header values are: text describing grid (64 bytes), num. columns (int),
     * num. rows (int),  num. z (int), min x (float), delta x (float), min y
     * (float), delta y (float)  and angle (float). Each record is num.
     * columns  4 bytes + 4 byte separator long and the file contains num.
     * rows + 1 (for the header) records. The data records (with the grid
     * shift values) are all floats and have a 4 byte  separator (0's) before
     * the data. Row records are organized from low  y (latitude) to high and
     * columns are orderd from low longitude to high.  Everything is written
     * in low byte order.
     *
     * @param latGridUrl URL to the binary latitude shift file (.las
     *        extention).
     * @param longGridUrl URL to the binary longitude shift file (.los
     *        extention).
     * @throws IOException if the data files cannot be read.
     * @throws FactoryException if there is an inconsistency in the data
     */
    private void loadBinaryGrid(URL latGridUrl, URL longGridUrl)
            throws IOException, FactoryException
    {
        final int HEADER_BYTES = 96;
        final int SEPARATOR_BYTES = 4;
        final int DESCRIPTION_LENGTH = 64;
        ReadableByteChannel latChannel;
        ReadableByteChannel longChannel;
        ByteBuffer latBuffer;
        ByteBuffer longBuffer;

        ////////////////////////
        //setup
        ////////////////////////
        latChannel = getReadChannel(latGridUrl);
        latBuffer = fillBuffer(latChannel, HEADER_BYTES);

        longChannel = getReadChannel(longGridUrl);
        longBuffer = fillBuffer(longChannel, HEADER_BYTES);

        ////////////////////////
        //read header info
        ////////////////////////
        //skip the header description
        latBuffer.position(latBuffer.position() + DESCRIPTION_LENGTH);

        int nc = latBuffer.getInt();
        int nr = latBuffer.getInt();
        int nz = latBuffer.getInt();

        xmin = latBuffer.getFloat();
        dx = latBuffer.getFloat();
        ymin = latBuffer.getFloat();
        dy = latBuffer.getFloat();

        float angle = latBuffer.getFloat();
        xmax = xmin + ((nc - 1) * dx);
        ymax = ymin + ((nr - 1) * dy);

        //skip the longitude header description
        longBuffer.position(longBuffer.position() + DESCRIPTION_LENGTH);

        //check that latitude grid header is the same as for latitude grid
        if ((nc != longBuffer.getInt()) || (nr != longBuffer.getInt())
                || (nz != longBuffer.getInt())
                || (xmin != longBuffer.getFloat())
                || (dx != longBuffer.getFloat())
                || (ymin != longBuffer.getFloat())
                || (dy != longBuffer.getFloat())
                || (angle != longBuffer.getFloat())) {
            throw new FactoryException(Resources.format(
                    ResourceKeys.ERROR_GRID_LOCATIONS_UNEQUAL));
        }

        ////////////////////////
        //read grid shift data into LocalizationGrid
        ////////////////////////
        final int RECORD_LENGTH = (nc * 4) + SEPARATOR_BYTES;
        final int NUM_BYTES_LEFT = ((nr + 1) * RECORD_LENGTH) - HEADER_BYTES;
        final int START_OF_DATA = RECORD_LENGTH - HEADER_BYTES;

        latBuffer = fillBuffer(latChannel, NUM_BYTES_LEFT);
        latBuffer.position(START_OF_DATA); //start of second record (data)

        longBuffer = fillBuffer(longChannel, NUM_BYTES_LEFT);
        longBuffer.position(START_OF_DATA);

        gridShift = new LocalizationGrid(nc, nr);

        int i = 0;
        int j = 0;
        for (i = 0; i < nr; i++) {
            latBuffer.position(latBuffer.position() + SEPARATOR_BYTES); //skip record separator
            longBuffer.position(longBuffer.position() + SEPARATOR_BYTES);

            for (j = 0; j < nc; j++) {
                gridShift.setLocalizationPoint(j, i,
                    (new Float(longBuffer.getFloat())).doubleValue(),
                    (new Float(latBuffer.getFloat())).doubleValue());
            }
        }

        assert i == nr : i;
        assert j == nc : j;
    }

    /**
     * Returns a new bytebuffer, of numBytes length and little endian byte
     * order, filled from the channel.
     *
     * @param channel the channel to fill the buffer from
     * @param numBytes number of bytes to read
     * @return a new bytebuffer filled from the channel
     * @throws IOException if there is a problem reading the channel
     * @throws EOFException if the end of the channel is reached
     */
    private ByteBuffer fillBuffer(ReadableByteChannel channel, int numBytes)
        throws IOException {
        ByteBuffer buf = ByteBuffer.allocateDirect(numBytes);

        if (fill(buf, channel) == -1) {
            throw new EOFException(Resources.format(
                    ResourceKeys.ERROR_END_OF_DATA_FILE));
        }

        buf.flip();
        buf.order(ByteOrder.LITTLE_ENDIAN);

        return buf;
    }

    /**
     * Fills the bytebuffer from the channel. Code was lifted from
     * ShapefileDataStore.
     *
     * @param buffer bytebuffer to fill from the channel
     * @param channel channel to fill the buffer from
     * @return number of bytes read
     * @throws IOException if there is a problem reading the channel
     */
    private int fill(ByteBuffer buffer, ReadableByteChannel channel)
        throws IOException {
        int r = buffer.remaining();

        // channel reads return -1 when EOF or other error
        // because they a non-blocking reads, 0 is a valid return value!!
        while ((buffer.remaining() > 0) && (r != -1)) {
            r = channel.read(buffer);
        }

        if (r == -1) {
            buffer.limit(buffer.position());
        }

        return r;
    }

    /**
     * Obtain a ReadableByteChannel from the given URL. If the url protocol is
     * file, a FileChannel will be returned. Otherwise a generic channel will
     * be obtained from the urls input stream. Code swiped from
     * ShapefileDataStore.
     *
     * @param url URL to create the channel from
     * @return a new PeadableByteChannel from the input url
     * @throws IOException if there is a problem creating the channel
     */
    private ReadableByteChannel getReadChannel(URL url)
        throws IOException {
        ReadableByteChannel channel = null;

        if (url.getProtocol().equals("file")) {
            File file = new File(url.getFile());

            if (!file.exists() || !file.canRead()) {
                throw new IOException(Resources.format(
                        ResourceKeys.ERROR_FILE_DOES_NOT_EXIST_$1, file));
            }

            FileInputStream in = new FileInputStream(file);
            channel = in.getChannel();
        } else {
            InputStream in = url.openConnection().getInputStream();
            channel = Channels.newChannel(in);
        }

        return channel;
    }

    /**
     * Reads latitude and longitude text grid shift file data into {@link
     * grid}.    The first two lines of the shift data file contain the
     * header, with the first being a description of the grid. The second line
     * contains 8 values separated by spaces: num. columns, num. rows, num. z,
     * min x, delta x, min y, delta y and angle. Shift data values follow this
     * and are also  separated by spaces. Row records are organized from low y
     * (latitude) to  high and columns are orderd from low longitude to high.
     *
     * @param latGridUrl URL to the text latitude shift file (.laa extention).
     * @param longGridUrl URL to the text longitude shift file (.loa
     *        extention).
     * @throws IOException if the data files cannot be read.
     * @throws FactoryException if there is an inconsistency in the data
     */
    private void loadTextGrid(URL latGridUrl, URL longGridUrl)
        throws IOException, FactoryException {
        String latLine;
        String longLine;
        StringTokenizer latSt;
        StringTokenizer longSt;

        ////////////////////////
        //setup
        ////////////////////////
        InputStreamReader latIsr = new InputStreamReader(latGridUrl.openStream());
        BufferedReader latBr = new BufferedReader(latIsr);

        InputStreamReader longIsr = new InputStreamReader(longGridUrl.openStream());
        BufferedReader longBr = new BufferedReader(longIsr);

        ////////////////////////
        //read header info
        ////////////////////////
        latLine = latBr.readLine(); //skip header description
        latLine = latBr.readLine();
        latSt = new StringTokenizer(latLine, " ");

        if (latSt.countTokens() != 8) {
            throw new FactoryException(Resources.format(
                    ResourceKeys.ERROR_HEADER_UNEXPECTED_LENGTH_$1,
                    String.valueOf(latSt.countTokens())));
        }

        int nc = Integer.parseInt(latSt.nextToken());
        int nr = Integer.parseInt(latSt.nextToken());
        int nz = Integer.parseInt(latSt.nextToken());

        xmin = Float.parseFloat(latSt.nextToken());
        dx = Float.parseFloat(latSt.nextToken());
        ymin = Float.parseFloat(latSt.nextToken());
        dy = Float.parseFloat(latSt.nextToken());

        float angle = Float.parseFloat(latSt.nextToken());
        xmax = xmin + ((nc - 1) * dx);
        ymax = ymin + ((nr - 1) * dy);

        //now read long shift grid
        longLine = longBr.readLine(); //skip header description
        longLine = longBr.readLine();
        longSt = new StringTokenizer(longLine, " ");

        if (longSt.countTokens() != 8) {
            throw new FactoryException(Resources.format(
                    ResourceKeys.ERROR_HEADER_UNEXPECTED_LENGTH_$1,
                    String.valueOf(longSt.countTokens())));
        }

        //check that latitude grid header is the same as for latitude grid
        if ((nc != Integer.parseInt(longSt.nextToken()))
                || (nr != Integer.parseInt(longSt.nextToken()))
                || (nz != Integer.parseInt(longSt.nextToken()))
                || (xmin != Float.parseFloat(longSt.nextToken()))
                || (dx != Float.parseFloat(longSt.nextToken()))
                || (ymin != Float.parseFloat(longSt.nextToken()))
                || (dy != Float.parseFloat(longSt.nextToken()))
                || (angle != Float.parseFloat(longSt.nextToken()))) {
            throw new FactoryException(Resources.format(
                    ResourceKeys.ERROR_GRID_LOCATIONS_UNEQUAL));
        }

        ////////////////////////
        //read grid shift data into LocalizationGrid
        ////////////////////////
        gridShift = new LocalizationGrid(nc, nr);

        int i = 0;
        int j = 0;
        for (i = 0; i < nr; i++) {
            for (j = 0; j < nc;) {
                latLine = latBr.readLine();
                latSt = new StringTokenizer(latLine, " ");
                longLine = longBr.readLine();
                longSt = new StringTokenizer(longLine, " ");

                while (latSt.hasMoreTokens() && longSt.hasMoreTokens()) {
                    gridShift.setLocalizationPoint(j, i,
                        (double) Float.parseFloat(longSt.nextToken()),
                        (double) Float.parseFloat(latSt.nextToken()));
                    ++j;
                }
            }
        }

        assert i == nr : i;
        assert j == nc : j;
    }

    /**
     * Transforms a list of coordinate point ordinal values. This method is
     * provided for efficiently transforming many points. The supplied array
     * of ordinal values will contain packed ordinal values.  For example, if
     * the source dimension is 3, then the ordinals will be packed in this
     * order:
     * (<var>x<sub>0</sub></var>,<var>y<sub>0</sub></var>,<var>z<sub>0</sub></var>,
     *
     * <var>x<sub>1</sub></var>,<var>y<sub>1</sub></var>,<var>z<sub>1</sub></var>
     * ...).  All input and output values are in degrees.
     *
     * @param srcPts the array containing the source point coordinates.
     * @param srcOff the offset to the first point to be transformed in the
     *        source array.
     * @param dstPts the array into which the transformed point coordinates are
     *        returned. May be the same than <code>srcPts</code>.
     * @param dstOff the offset to the location of the first transformed point
     *        that is stored in the destination array.
     * @param numPts the number of point objects to be transformed.
     *
     * @throws TransformException if the input point is outside the area
     *         covered by this grid.
     */
    public void transform(final double[] srcPts, int srcOff,
        final double[] dstPts, int dstOff, int numPts)
        throws TransformException {
        int step = 0;

        if ((srcPts == dstPts) && (srcOff < dstOff)
                && ((srcOff + (numPts * getDimSource())) > dstOff)) {
            step = -getDimSource();
            srcOff -= ((numPts - 1) * step);
            dstOff -= ((numPts - 1) * step);
        }

        while (--numPts >= 0) {
            double x = srcPts[srcOff++];
            double y = srcPts[srcOff++];

            //check bounding box
//issue of bbox crossing +- 180 degrees (ie input point of -188 longitude);
//abs(x - xmin) > 0 , rollLongitude() ???
            if (((x < xmin) || (x > xmax)) || ((y < ymin) || (y > ymax))) {
                throw new TransformException(Resources.format(
                        ResourceKeys.ERROR_POINT_OUTSIDE_GRID));
            }

            //find the grid the point is in (index is 0 based)
            final double xgrid = (x - xmin) / dx;
            final double ygrid = (y - ymin) / dy;
            double[] array = new double[] { xgrid, ygrid };

            //use the LocalizationGridTransform2D transform method (bilineal interpolation)
            //returned shift values are in seconds, longitude shift values are + west
            gridShiftTransform.transform(array, 0, array, 0, 1);

            dstPts[dstOff++] = x - (array[0] / SEC_2_DEG);
            dstPts[dstOff++] = y + (array[1] / SEC_2_DEG);
            srcOff += step;
            dstOff += step;
        }
    }

    /**
     * Transforms nad83 values to nad27. Input and output values are in
     * degrees.  This is done by itteratively finding a nad27 value that
     * shifts to the  input nad83 value. The input nad83 value is used as the
     * first  approximation.
     *
     * @param srcPts the array containing the source point coordinates.
     * @param srcOff the offset to the first point to be transformed in the
     *        source array.
     * @param dstPts the array into which the transformed point coordinates are
     *        returned. May be the same than <code>srcPts</code>.
     * @param dstOff the offset to the location of the first transformed point
     *        that is stored in the destination array.
     * @param numPts the number of point objects to be transformed.
     *
     * @throws TransformException if the input point is outside the area
     *         covered by this grid.
     */
    public void inverseTransform(final double[] srcPts, int srcOff,
        final double[] dstPts, int dstOff, int numPts)
        throws TransformException {
        int step = 0;

        if ((srcPts == dstPts) && (srcOff < dstOff)
                && ((srcOff + (numPts * getDimSource())) > dstOff)) {
            step = -getDimSource();
            srcOff -= ((numPts - 1) * step);
            dstOff -= ((numPts - 1) * step);
        }

        while (--numPts >= 0) {
            final double x = srcPts[srcOff++];
            final double y = srcPts[srcOff++];
            double xtemp = x;
            double ytemp = y;

            for (int i = MAX_ITER;;) {
                double[] array = { xtemp, ytemp };
                transform(array, 0, array, 0, 1);
                double xdif = array[0] - x;
                double ydif = array[1] - y;

                if (Math.abs(xdif) > TOL) {
                    xtemp = xtemp - xdif;
                }
                if (Math.abs(ydif) > TOL) {
                    ytemp = ytemp - ydif;
                }

                if ((Math.abs(xdif) <= TOL) && (Math.abs(ydif) <= TOL)) {
                    dstPts[dstOff++] = xtemp;
                    dstPts[dstOff++] = ytemp;
                    break;
                }
                if (--i < 0) {
                    throw new TransformException(Resources.format(
                            ResourceKeys.ERROR_NO_CONVERGENCE));
                }
            }

            srcOff += step;
            dstOff += step;
        }
    }

    /**
     * Returns the inverse of this transform.
     *
     * @return the inverse of this transform
     */
    public MathTransform inverse() {
        if (inverse == null) {
            // No need to synchronize; this is not a big deal if this object is created twice.
            inverse = new Inverse();
        }

        return inverse;
    }

    /**
     * Returns a hash value for this transform. To make this faster it does not
     * check the grid values.
     *
     * @return a hash value for this transform.
     */
    public final int hashCode() {
        final long code = Double.doubleToLongBits(xmin)
            + (37 * (Double.doubleToLongBits(ymin)
            + (37 * (Double.doubleToLongBits(xmax)
            + (37 * (Double.doubleToLongBits(ymax)
            + (37 * (Double.doubleToLongBits(dx)
            + (37 * (Double.doubleToLongBits(dy)))))))))));

        return (int) code ^ (int) (code >>> 32);
    }

    /**
     * Compares the specified object with this math transform for equality.
     *
     * @param object the object to compare to
     * @return <code>true</code> if the objects are equal.
     */
    public final boolean equals(final Object object) {
        if (object == this) {
            // Slight optimization
            return true;
        }

        if (super.equals(object)) {
            final NADCONTransform that = (NADCONTransform) object;

            return (Double.doubleToLongBits(this.xmin) == Double.doubleToLongBits(that.xmin))
                && (Double.doubleToLongBits(this.ymin) == Double.doubleToLongBits(that.ymin))
                && (Double.doubleToLongBits(this.xmax) == Double.doubleToLongBits(that.xmax))
                && (Double.doubleToLongBits(this.ymax) == Double.doubleToLongBits(that.ymax))
                && (Double.doubleToLongBits(this.dx) == Double.doubleToLongBits(that.dx))
                && (Double.doubleToLongBits(this.dy) == Double.doubleToLongBits(that.dy))
                && (this.gridShiftTransform).equals(that.gridShiftTransform);
        }

        return false;
    }

    /**
     * Used to set the preference for the default grid shift file location.
     * This allows grids parameters to be specified by name only, without the
     * full path. This needs to be done only once, by the superuser (root or
     * admin). Path values may be simple file system paths or more complex
     * text  representations of a url. A value of "default" resets this
     * preference to its default value.
     * <br><br>
     *
     * Example:
     * <blockquote>
     * <pre>
     * java org.geotools.referencing.operation.transform.NADCONTransform file:///home/rschulz/GIS/NADCON/data
     * </pre>
     * </blockquote>
     *
     * @param args a single argument for the defualt location of grid shift
     *        files
     */
    public static void main(String[] args) {
        final Arguments arguments = new Arguments(args);
        final PrintWriter out = arguments.out;
        final Preferences prefs = Preferences.systemNodeForPackage(NADCONTransform.class);

        if (args.length == 1) {
            if (args[0].equalsIgnoreCase("default")) {
                prefs.remove(GRID_LOCATION);
            } else {
                prefs.put(GRID_LOCATION, args[0]);
            }

            return;
        } else {
            final String location = prefs.get(GRID_LOCATION,
                    DEFAULT_GRID_LOCATION);
            out.println(
                "Usage: java org.geotools.referencing.operation.transform.NADCONTransform "
                + "<defalult grid file location (path)>");
            out.print("Grid location: \"");
            out.print(location);
            out.println('"');

            return;
        }
    }

    /**
     * Inverse of a {@link NADCONTransform}.
     *
     * @version $Id$
     * @author Rueben Schulz
     */
    private final class Inverse extends AbstractMathTransform.Inverse implements Serializable {
        /** Serial number for interoperability with different versions. */
        private static final long serialVersionUID = -4707304160205218546L;

        /**
         * Default constructor.
         */
        public Inverse() {
            NADCONTransform.this.super();
        }

        /**
         * Returns the parameter values for this math transform.
         *
         * @return A copy of the parameter values for this math transform.
         */
        public ParameterValueGroup getParameterValues() {
            return null;
        }

        /**
         * Inverse transform an array of points.
         *
         * @param source
         * @param srcOffset
         * @param dest
         * @param dstOffset
         * @param length
         *
         * @throws TransformException if the input point is outside the area
         *         covered by this grid.
         */
        public void transform(final double[] source, final int srcOffset,
            final double[] dest, final int dstOffset, final int length)
            throws TransformException {
            NADCONTransform.this.inverseTransform(source, srcOffset, dest,
                dstOffset, length);
        }

        /**
         * Restore reference to this object after deserialization.
         *
         * @param in DOCUMENT ME!
         * @throws IOException DOCUMENT ME!
         * @throws ClassNotFoundException DOCUMENT ME!
         */
        private void readObject(ObjectInputStream in)
            throws IOException, ClassNotFoundException {
            in.defaultReadObject();
            NADCONTransform.this.inverse = this;
        }
    }

    /**
     * The provider for {@link NADCONTransform}. This provider will construct
     * transforms from {@linkPlain org.geotools.referencing.crs.GeographicCRS
     * geographic} to {@linkPlain org.geotools.referencing.crs.GeographicCRS
     * geographic} coordinate reference systems.
     *
     * @version $Id$
     * @author Rueben Schulz
     */
    public static class Provider extends MathTransformProvider {
        /** Serial number for interoperability with different versions. */
        private static final long serialVersionUID = -4707304160205218546L;

        /**
         * The operation parameter descriptor for the "Latitude_difference_file"
         * parameter value. The default value is "conus.las".
         */
        public static final ParameterDescriptor LAT_DIFF_FILE
            = new org.geotools.parameter.ParameterDescriptor(
                "Latitude_difference_file", "".getClass(), null, "conus.las");

        /**
         * The operation parameter descriptor for the "Longitude_difference_file"
         * parameter value. The default value is "conus.los".
         */
        public static final ParameterDescriptor LONG_DIFF_FILE
            = new org.geotools.parameter.ParameterDescriptor(
                "Longitude_difference_file", new String().getClass(), null, "conus.los");

        /** The parameters group. */
        static final ParameterDescriptorGroup PARAMETERS = group(new Identifier[] {
                    new Identifier(Citation.OPEN_GIS, null, "NADCON"),
                    new Identifier(Citation.EPSG, "EPSG", "9613")
                }, new ParameterDescriptor[] { LAT_DIFF_FILE, LONG_DIFF_FILE });

        /**
         * Constructs a provider.
         */
        public Provider() {
            super(2, 2, PARAMETERS);
        }

        /**
         * Returns the resources key for {@linkPlain #getName localized name}.
         * This method is for internal purpose by Geotools implementation
         * only.
         *
         * @return the resource key for this math transform
         */
        protected int getLocalizationKey() {
            return ResourceKeys.NADCON_TRANSFORM;
        }

        /**
         * Creates a math transform from the specified group of parameter
         * values.
         *
         * @param values The group of parameter values.
         * @return The created math transform.
         * @throws ParameterNotFoundException if a required parameter was not
         *         found.
         * @throws FactoryException if there is a problem creating this
         *         math transform.
         */
        protected MathTransform createMathTransform(ParameterValueGroup values)
            throws ParameterNotFoundException, FactoryException {
            return new NADCONTransform(
                stringValue(values, LAT_DIFF_FILE),
                stringValue(values, LONG_DIFF_FILE));
        }
    }
}