0001: /*
0002: * GeoTools - OpenSource mapping toolkit
0003: * http://geotools.org
0004: *
0005: * (C) 2004-2006, Geotools Project Managment Committee (PMC)
0006: *
0007: * This library is free software; you can redistribute it and/or
0008: * modify it under the terms of the GNU Lesser General Public
0009: * License as published by the Free Software Foundation; either
0010: * version 2.1 of the License, or (at your option) any later version.
0011: *
0012: * This library is distributed in the hope that it will be useful,
0013: * but WITHOUT ANY WARRANTY; without even the implied warranty of
0014: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
0015: * Lesser General Public License for more details.
0016: */
0017: package org.geotools.referencing.operation.transform;
0018:
0019: // J2SE dependencies
0020: import java.io.BufferedReader;
0021: import java.io.EOFException;
0022: import java.io.File;
0023: import java.io.FileInputStream;
0024: import java.io.IOException;
0025: import java.io.InputStream;
0026: import java.io.InputStreamReader;
0027: import java.io.ObjectInputStream;
0028: import java.io.PrintWriter;
0029: import java.io.Serializable;
0030: import java.net.URL;
0031: import java.net.MalformedURLException;
0032: import java.nio.ByteBuffer;
0033: import java.nio.ByteOrder;
0034: import java.nio.channels.Channels;
0035: import java.nio.channels.ReadableByteChannel;
0036: import java.util.StringTokenizer;
0037: import java.util.prefs.Preferences;
0038:
0039: // OpenGIS dependencies
0040: import org.opengis.parameter.GeneralParameterValue;
0041: import org.opengis.parameter.ParameterDescriptor;
0042: import org.opengis.parameter.ParameterDescriptorGroup;
0043: import org.opengis.parameter.ParameterNotFoundException;
0044: import org.opengis.parameter.ParameterValue;
0045: import org.opengis.parameter.ParameterValueGroup;
0046: import org.opengis.referencing.FactoryException;
0047: import org.opengis.referencing.operation.MathTransform;
0048: import org.opengis.referencing.operation.MathTransform2D;
0049: import org.opengis.referencing.operation.Transformation;
0050: import org.opengis.referencing.operation.TransformException;
0051:
0052: // Geotools dependencies
0053: import org.geotools.metadata.iso.citation.Citations;
0054: import org.geotools.parameter.DefaultParameterDescriptor;
0055: import org.geotools.parameter.Parameter;
0056: import org.geotools.parameter.ParameterGroup;
0057: import org.geotools.referencing.NamedIdentifier;
0058: import org.geotools.referencing.operation.MathTransformProvider;
0059: import org.geotools.referencing.operation.builder.LocalizationGrid;
0060: import org.geotools.resources.Arguments;
0061: import org.geotools.resources.i18n.Errors;
0062: import org.geotools.resources.i18n.ErrorKeys;
0063: import org.geotools.resources.i18n.Vocabulary;
0064: import org.geotools.resources.i18n.VocabularyKeys;
0065:
0066: /**
0067: * Transform backed by the North American Datum Conversion grid.
0068: * The North American Datum Conversion (NADCON) Transform (EPSG code 9613) is a
0069: * two dimentional datum shift method, created by the National Geodetic Survey
0070: * (NGS), that uses interpolated values from two grid shift files. This
0071: * method is used to transform NAD27 (EPSG code 4267) datum coordinates (and
0072: * some others) to NAD83 (EPSG code 4267) within the United States. There are
0073: * two set of grid shift files: NADCON and High Accuracy Reference Networks
0074: * (HARN). NADCON shfts from NAD27 (and some others) to NAD83 while HARN
0075: * shifts from the NADCON NAD83 to an improved NAD83. Both sets of grid shift
0076: * files may be downloaded from
0077: * <a href="http://www.ngs.noaa.gov/PC_PROD/NADCON/">www.ngs.noaa.gov/PC_PROD/NADCON/</a>.
0078: * <p>
0079: *
0080: * Some of the NADCON grids, their areas of use, and source datums are shown
0081: * in the following table.
0082: * <p>
0083: *
0084: * <table>
0085: * <tr><th>Shift File Name</td><th>Area</td><th>Source Datum</td><th>Accuracy at 67% confidence (m)</td></tr>
0086: * <tr><td>CONUS</td><td>Conterminous U S (lower 48 states)</td><td>NAD27</td><td>0.15</td></tr>
0087: * <tr><td>ALASKA</td><td>Alaska, incl. Aleutian Islands</td><td>NAD27</td><td>0.5</td></tr>
0088: * <tr><td>HAWAII</td><td>Hawaiian Islands</td><td>Old Hawaiian (4135)</td><td>0.2</td></tr>
0089: * <tr><td>STLRNC</td><td>St. Lawrence Is., AK</td><td>St. Lawrence Island (4136)</td><td>--</td></tr>
0090: * <tr><td>STPAUL </td><td>St. Paul Is., AK</td><td>St. Paul Island (4137)</td><td>--</td></tr>
0091: * <tr><td>STGEORGE</td><td>St. George Is., AK</td><td>St. George Island (4138)</td><td>--</td></tr>
0092: * <tr><td>PRVI</td><td>Puerto Rico and the Virgin Islands</td><td>Puerto Rico (4139)</td><td>0.05</td></tr>
0093: * </table>
0094: * <p>
0095: *
0096: * Grid shift files come in two formats: binary and text. The files from the NGS are
0097: * binary and have {@code .las} (latitude shift) and {@code .los} (longitude shift)
0098: * extentions. Text grids may be created with the <cite>NGS nadgrd</cite> program and have
0099: * {@code .laa} (latitude shift) and {@code .loa} (longitude shift) file extentions.
0100: * Both types of files may be used here.
0101: * <p>
0102: *
0103: * The grid names to use for transforming are parameters of this
0104: * {@link MathTransform}. This parameter may be the full name and path to the grids
0105: * or just the name of the grids if the default location of the grids was set
0106: * as a preference. This preference may be set with the main method of this
0107: * class.
0108: * <p>
0109: *
0110: * Transformations here have been tested to be within 0.00001 seconds of
0111: * values given by the <cite>NGS ndcon210</cite> program for NADCON grids. American Samoa
0112: * and HARN shifts have not yet been tested. <strong>References:</strong>
0113: *
0114: * <ul>
0115: * <li><a href="http://www.ngs.noaa.gov/PC_PROD/NADCON/Readme.htm">NADCONreadme</a></li>
0116: * <li>American Samoa Grids for NADCON - Samoa_Readme.txt</li>
0117: * <li><a href="http://www.ngs.noaa.gov/PUBS_LIB/NGS50.pdf">NADCON - The
0118: * Application of Minimum-Curvature-Derived Surfaces in the Transformation of
0119: * Positional Data From the North American Datum of 1927 to the North
0120: * American Datum of 1983</a> - NOAA TM.</li>
0121: * <li>{@code ndcon210.for} - NGS fortran source code for NADCON conversions. See the
0122: * following subroutines: TRANSF, TO83, FGRID, INTRP, COEFF and SURF</li>
0123: * <li>{@code nadgrd.for} - NGS fortran source code to export/import binary and text grid
0124: * formats</li>
0125: * <li>EPSG Geodesy Parameters database version 6.5</li>
0126: * </ul>
0127: *
0128: * @see <a href="http://www.ngs.noaa.gov/TOOLS/Nadcon/Nadcon.html"> NADCON -
0129: * North American Datum Conversion Utility</a>
0130: *
0131: * @since 2.1
0132: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/transform/NADCONTransform.java $
0133: * @version $Id: NADCONTransform.java 29101 2008-02-06 12:11:19Z desruisseaux $
0134: * @author Rueben Schulz
0135: *
0136: * @todo the transform code does not deal with the case where grids cross +- 180 degrees.
0137: */
0138: public class NADCONTransform extends AbstractMathTransform implements
0139: MathTransform2D, Serializable {
0140: /**
0141: * Serial number for interoperability with different versions.
0142: */
0143: private static final long serialVersionUID = -4707304160205218546L;
0144:
0145: /**
0146: * Preference node for the grid shift file location.
0147: */
0148: private static final String GRID_LOCATION = "Grid location";
0149:
0150: /**
0151: * The default value for the grid shift file location.
0152: */
0153: private static final String DEFAULT_GRID_LOCATION = ".";
0154:
0155: /**
0156: * Difference allowed in iterative computations. This is half the value
0157: * used in the NGS fortran code (so all tests pass).
0158: */
0159: private static final double TOL = 5.0E-10;
0160:
0161: /**
0162: * Maximum number of iterations for iterative computations.
0163: */
0164: private static final int MAX_ITER = 10;
0165:
0166: /**
0167: * Conversion factor from seconds to decimal degrees.
0168: */
0169: private static final double SEC_2_DEG = 3600.0;
0170:
0171: /**
0172: * Latitude grid shift file names. Output in WKT.
0173: */
0174: private final String latGridName;
0175:
0176: /**
0177: * Longitude grid shift file names. Output in WKT.
0178: */
0179: private final String longGridName;
0180:
0181: /**
0182: * The minimum longitude value covered by this grid (decimal degrees)
0183: */
0184: private double xmin;
0185:
0186: /**
0187: * The minimum latitude value covered by this grid (decimal degrees)
0188: */
0189: private double ymin;
0190:
0191: /**
0192: * The maximum longitude value covered by this grid (decimal degrees)
0193: */
0194: private double xmax;
0195:
0196: /**
0197: * The maximum latitude value covered by this grid (decimal degrees)
0198: */
0199: private double ymax;
0200:
0201: /**
0202: * The difference between longitude grid points (decimal degrees)
0203: */
0204: private double dx;
0205:
0206: /**
0207: * The difference between latitude grid points (decimal degrees)
0208: */
0209: private double dy;
0210:
0211: /**
0212: * Longitude and latitude grid shift values. Values are organized from low
0213: * to high longitude (low x index to high) and low to high latitude (low y
0214: * index to high).
0215: */
0216: private LocalizationGrid gridShift;
0217:
0218: /**
0219: * The {@link #gridShift} values as a {@code LocalizationGridTransform2D}.
0220: * Used for interpolating shift values.
0221: */
0222: private MathTransform gridShiftTransform;
0223:
0224: /**
0225: * The inverse of this transform. Will be created only when needed.
0226: */
0227: private transient MathTransform inverse;
0228:
0229: /**
0230: * Constructs a {@code NADCONTransform} from the specified grid shift files.
0231: *
0232: * @param latGridName path and name (or just name if {@link #GRID_LOCATION}
0233: * is set) to the latitude difference file. This will have a {@code .las} or
0234: * {@code .laa} file extention.
0235: * @param longGridName path and name (or just name if {@link #GRID_LOCATION}
0236: * is set) to the longitude difference file. This will have a {@code .los}
0237: * or {@code .loa} file extention.
0238: *
0239: * @throws ParameterNotFoundException if a math transform parameter cannot be found.
0240: * @throws FactoryException if there is a problem creating this math transform
0241: * (ie file extentions are unknown or there is an error reading the
0242: * grid files)
0243: */
0244: public NADCONTransform(final String latGridName,
0245: final String longGridName)
0246: throws ParameterNotFoundException, FactoryException {
0247: this .latGridName = latGridName;
0248: this .longGridName = longGridName;
0249:
0250: //decide if text or binary grid will be used
0251: try {
0252: final URL latGridURL = makeURL(latGridName);
0253: final URL longGridURL = makeURL(longGridName);
0254:
0255: if ((latGridName.endsWith(".las") && longGridName
0256: .endsWith(".los"))
0257: || (latGridName.endsWith(".LAS") && longGridName
0258: .endsWith(".LOS"))) {
0259: loadBinaryGrid(latGridURL, longGridURL);
0260: } else if ((latGridName.endsWith(".laa") && longGridName
0261: .endsWith(".loa"))
0262: || (latGridName.endsWith(".LAA") && longGridName
0263: .endsWith(".LOA"))) {
0264: loadTextGrid(latGridURL, longGridURL);
0265: } else {
0266: throw new FactoryException(Errors.format(
0267: ErrorKeys.UNSUPPORTED_FILE_TYPE_$2,
0268: latGridName.substring(latGridName
0269: .lastIndexOf('.') + 1), longGridName
0270: .substring(longGridName
0271: .lastIndexOf('.') + 1)));
0272: // Note: the +1 above hide the dot, but also make sure that the code is
0273: // valid even if the path do not contains '.' at all (-1 + 1 == 0).
0274: }
0275:
0276: gridShiftTransform = gridShift.getMathTransform();
0277: } catch (IOException exception) {
0278: final Throwable cause = exception.getCause();
0279: if (cause instanceof FactoryException) {
0280: throw (FactoryException) cause;
0281: }
0282: throw new FactoryException(exception.getLocalizedMessage(),
0283: exception);
0284: }
0285: }
0286:
0287: /**
0288: * Returns the parameter descriptors for this math transform.
0289: */
0290: public ParameterDescriptorGroup getParameterDescriptors() {
0291: return Provider.PARAMETERS;
0292: }
0293:
0294: /**
0295: * Returns the parameter values for this math transform.
0296: *
0297: * @return A copy of the parameter values for this math transform.
0298: */
0299: public ParameterValueGroup getParameterValues() {
0300: final ParameterValue lat_diff_file = new Parameter(
0301: Provider.LAT_DIFF_FILE);
0302: lat_diff_file.setValue(latGridName);
0303:
0304: final ParameterValue long_diff_file = new Parameter(
0305: Provider.LONG_DIFF_FILE);
0306: long_diff_file.setValue(longGridName);
0307:
0308: return new ParameterGroup(getParameterDescriptors(),
0309: new GeneralParameterValue[] { lat_diff_file,
0310: long_diff_file });
0311: }
0312:
0313: /**
0314: * Gets the dimension of input points (always 2).
0315: *
0316: * @return the source dimensions.
0317: */
0318: public int getSourceDimensions() {
0319: return 2;
0320: }
0321:
0322: /**
0323: * Gets the dimension of output points (always 2).
0324: *
0325: * @return the target dimensions.
0326: */
0327: public int getTargetDimensions() {
0328: return 2;
0329: }
0330:
0331: /**
0332: * Returns a URL from the string representation. If the string has no
0333: * path, the default path preferece is added.
0334: *
0335: * @param str a string representation of a URL
0336: * @return a URL created from the string representation
0337: * @throws MalformedURLException if the URL cannot be created
0338: */
0339: private URL makeURL(final String str) throws MalformedURLException {
0340: //has '/' or '\' or ':', so probably full path to file
0341: if ((str.indexOf('\\') >= 0) || (str.indexOf('/') >= 0)
0342: || (str.indexOf(':') >= 0)) {
0343: return makeURLfromString(str);
0344: } else {
0345: // just a file name , prepend base location
0346: final Preferences prefs = Preferences
0347: .userNodeForPackage(NADCONTransform.class);
0348: final String baseLocation = prefs.get(GRID_LOCATION,
0349: DEFAULT_GRID_LOCATION);
0350: return makeURLfromString(baseLocation + "/" + str);
0351: }
0352: }
0353:
0354: /**
0355: * Returns a URL based on a string representation. If no protocol is given,
0356: * it is assumed to be a local file.
0357: *
0358: * @param str a string representation of a URL
0359: * @return a URL created from the string representation
0360: * @throws MalformedURLException if the URL cannot be created
0361: */
0362: private URL makeURLfromString(final String str)
0363: throws MalformedURLException {
0364: try {
0365: return new URL(str);
0366: } catch (MalformedURLException e) {
0367: //try making this with a file protocal
0368: return new URL("file", "", str);
0369: }
0370: }
0371:
0372: /**
0373: * Reads latitude and longitude binary grid shift file data into {@link
0374: * grid}. The file is organized into records, with the first record
0375: * containing the header information, followed by the shift data. The
0376: * header values are: text describing grid (64 bytes), num. columns (int),
0377: * num. rows (int), num. z (int), min x (float), delta x (float), min y
0378: * (float), delta y (float) and angle (float). Each record is num.
0379: * columns 4 bytes + 4 byte separator long and the file contains num.
0380: * rows + 1 (for the header) records. The data records (with the grid
0381: * shift values) are all floats and have a 4 byte separator (0's) before
0382: * the data. Row records are organized from low y (latitude) to high and
0383: * columns are orderd from low longitude to high. Everything is written
0384: * in low byte order.
0385: *
0386: * @param latGridUrl URL to the binary latitude shift file (.las extention).
0387: * @param longGridUrl URL to the binary longitude shift file (.los extention).
0388: * @throws IOException if the data files cannot be read.
0389: * @throws FactoryException if there is an inconsistency in the data
0390: */
0391: private void loadBinaryGrid(final URL latGridUrl,
0392: final URL longGridUrl) throws IOException, FactoryException {
0393: final int HEADER_BYTES = 96;
0394: final int SEPARATOR_BYTES = 4;
0395: final int DESCRIPTION_LENGTH = 64;
0396: ReadableByteChannel latChannel;
0397: ReadableByteChannel longChannel;
0398: ByteBuffer latBuffer;
0399: ByteBuffer longBuffer;
0400:
0401: ////////////////////////
0402: //setup
0403: ////////////////////////
0404: latChannel = getReadChannel(latGridUrl);
0405: latBuffer = fillBuffer(latChannel, HEADER_BYTES);
0406:
0407: longChannel = getReadChannel(longGridUrl);
0408: longBuffer = fillBuffer(longChannel, HEADER_BYTES);
0409:
0410: ////////////////////////
0411: //read header info
0412: ////////////////////////
0413: //skip the header description
0414: latBuffer.position(latBuffer.position() + DESCRIPTION_LENGTH);
0415:
0416: int nc = latBuffer.getInt();
0417: int nr = latBuffer.getInt();
0418: int nz = latBuffer.getInt();
0419:
0420: xmin = latBuffer.getFloat();
0421: dx = latBuffer.getFloat();
0422: ymin = latBuffer.getFloat();
0423: dy = latBuffer.getFloat();
0424:
0425: float angle = latBuffer.getFloat();
0426: xmax = xmin + ((nc - 1) * dx);
0427: ymax = ymin + ((nr - 1) * dy);
0428:
0429: //skip the longitude header description
0430: longBuffer.position(longBuffer.position() + DESCRIPTION_LENGTH);
0431:
0432: //check that latitude grid header is the same as for latitude grid
0433: if ((nc != longBuffer.getInt()) || (nr != longBuffer.getInt())
0434: || (nz != longBuffer.getInt())
0435: || (xmin != longBuffer.getFloat())
0436: || (dx != longBuffer.getFloat())
0437: || (ymin != longBuffer.getFloat())
0438: || (dy != longBuffer.getFloat())
0439: || (angle != longBuffer.getFloat())) {
0440: throw new FactoryException(Errors
0441: .format(ErrorKeys.GRID_LOCATIONS_UNEQUAL));
0442: }
0443:
0444: ////////////////////////
0445: //read grid shift data into LocalizationGrid
0446: ////////////////////////
0447: final int RECORD_LENGTH = (nc * 4) + SEPARATOR_BYTES;
0448: final int NUM_BYTES_LEFT = ((nr + 1) * RECORD_LENGTH)
0449: - HEADER_BYTES;
0450: final int START_OF_DATA = RECORD_LENGTH - HEADER_BYTES;
0451:
0452: latBuffer = fillBuffer(latChannel, NUM_BYTES_LEFT);
0453: latBuffer.position(START_OF_DATA); //start of second record (data)
0454:
0455: longBuffer = fillBuffer(longChannel, NUM_BYTES_LEFT);
0456: longBuffer.position(START_OF_DATA);
0457:
0458: gridShift = new LocalizationGrid(nc, nr);
0459:
0460: int i = 0;
0461: int j = 0;
0462: for (i = 0; i < nr; i++) {
0463: latBuffer.position(latBuffer.position() + SEPARATOR_BYTES); //skip record separator
0464: longBuffer
0465: .position(longBuffer.position() + SEPARATOR_BYTES);
0466:
0467: for (j = 0; j < nc; j++) {
0468: gridShift
0469: .setLocalizationPoint(j, i, (new Float(
0470: longBuffer.getFloat())).doubleValue(),
0471: (new Float(latBuffer.getFloat()))
0472: .doubleValue());
0473: }
0474: }
0475:
0476: assert i == nr : i;
0477: assert j == nc : j;
0478: }
0479:
0480: /**
0481: * Returns a new bytebuffer, of numBytes length and little endian byte
0482: * order, filled from the channel.
0483: *
0484: * @param channel the channel to fill the buffer from
0485: * @param numBytes number of bytes to read
0486: * @return a new bytebuffer filled from the channel
0487: * @throws IOException if there is a problem reading the channel
0488: * @throws EOFException if the end of the channel is reached
0489: */
0490: private ByteBuffer fillBuffer(ReadableByteChannel channel,
0491: int numBytes) throws IOException {
0492: ByteBuffer buf = ByteBuffer.allocateDirect(numBytes);
0493:
0494: if (fill(buf, channel) == -1) {
0495: throw new EOFException(Errors
0496: .format(ErrorKeys.END_OF_DATA_FILE));
0497: }
0498:
0499: buf.flip();
0500: buf.order(ByteOrder.LITTLE_ENDIAN);
0501:
0502: return buf;
0503: }
0504:
0505: /**
0506: * Fills the bytebuffer from the channel. Code was lifted from
0507: * ShapefileDataStore.
0508: *
0509: * @param buffer bytebuffer to fill from the channel
0510: * @param channel channel to fill the buffer from
0511: * @return number of bytes read
0512: * @throws IOException if there is a problem reading the channel
0513: */
0514: private int fill(ByteBuffer buffer, ReadableByteChannel channel)
0515: throws IOException {
0516: int r = buffer.remaining();
0517:
0518: // channel reads return -1 when EOF or other error
0519: // because they a non-blocking reads, 0 is a valid return value!!
0520: while ((buffer.remaining() > 0) && (r != -1)) {
0521: r = channel.read(buffer);
0522: }
0523:
0524: if (r == -1) {
0525: buffer.limit(buffer.position());
0526: }
0527:
0528: return r;
0529: }
0530:
0531: /**
0532: * Obtain a ReadableByteChannel from the given URL. If the url protocol is
0533: * file, a FileChannel will be returned. Otherwise a generic channel will
0534: * be obtained from the urls input stream. Code swiped from
0535: * ShapefileDataStore.
0536: *
0537: * @param url URL to create the channel from
0538: * @return a new PeadableByteChannel from the input url
0539: * @throws IOException if there is a problem creating the channel
0540: */
0541: private ReadableByteChannel getReadChannel(URL url)
0542: throws IOException {
0543: ReadableByteChannel channel = null;
0544:
0545: if (url.getProtocol().equals("file")) {
0546: File file = new File(url.getFile());
0547:
0548: if (!file.exists() || !file.canRead()) {
0549: throw new IOException(Errors.format(
0550: ErrorKeys.FILE_DOES_NOT_EXIST_$1, file));
0551: }
0552:
0553: FileInputStream in = new FileInputStream(file);
0554: channel = in.getChannel();
0555: } else {
0556: InputStream in = url.openConnection().getInputStream();
0557: channel = Channels.newChannel(in);
0558: }
0559:
0560: return channel;
0561: }
0562:
0563: /**
0564: * Reads latitude and longitude text grid shift file data into {@link
0565: * grid}. The first two lines of the shift data file contain the
0566: * header, with the first being a description of the grid. The second line
0567: * contains 8 values separated by spaces: num. columns, num. rows, num. z,
0568: * min x, delta x, min y, delta y and angle. Shift data values follow this
0569: * and are also separated by spaces. Row records are organized from low y
0570: * (latitude) to high and columns are orderd from low longitude to high.
0571: *
0572: * @param latGridUrl URL to the text latitude shift file (.laa extention).
0573: * @param longGridUrl URL to the text longitude shift file (.loa
0574: * extention).
0575: * @throws IOException if the data files cannot be read.
0576: * @throws FactoryException if there is an inconsistency in the data
0577: */
0578: private void loadTextGrid(URL latGridUrl, URL longGridUrl)
0579: throws IOException, FactoryException {
0580: String latLine;
0581: String longLine;
0582: StringTokenizer latSt;
0583: StringTokenizer longSt;
0584:
0585: ////////////////////////
0586: //setup
0587: ////////////////////////
0588: InputStreamReader latIsr = new InputStreamReader(latGridUrl
0589: .openStream());
0590: BufferedReader latBr = new BufferedReader(latIsr);
0591:
0592: InputStreamReader longIsr = new InputStreamReader(longGridUrl
0593: .openStream());
0594: BufferedReader longBr = new BufferedReader(longIsr);
0595:
0596: ////////////////////////
0597: //read header info
0598: ////////////////////////
0599: latLine = latBr.readLine(); //skip header description
0600: latLine = latBr.readLine();
0601: latSt = new StringTokenizer(latLine, " ");
0602:
0603: if (latSt.countTokens() != 8) {
0604: throw new FactoryException(Errors.format(
0605: ErrorKeys.HEADER_UNEXPECTED_LENGTH_$1, String
0606: .valueOf(latSt.countTokens())));
0607: }
0608:
0609: int nc = Integer.parseInt(latSt.nextToken());
0610: int nr = Integer.parseInt(latSt.nextToken());
0611: int nz = Integer.parseInt(latSt.nextToken());
0612:
0613: xmin = Float.parseFloat(latSt.nextToken());
0614: dx = Float.parseFloat(latSt.nextToken());
0615: ymin = Float.parseFloat(latSt.nextToken());
0616: dy = Float.parseFloat(latSt.nextToken());
0617:
0618: float angle = Float.parseFloat(latSt.nextToken());
0619: xmax = xmin + ((nc - 1) * dx);
0620: ymax = ymin + ((nr - 1) * dy);
0621:
0622: //now read long shift grid
0623: longLine = longBr.readLine(); //skip header description
0624: longLine = longBr.readLine();
0625: longSt = new StringTokenizer(longLine, " ");
0626:
0627: if (longSt.countTokens() != 8) {
0628: throw new FactoryException(Errors.format(
0629: ErrorKeys.HEADER_UNEXPECTED_LENGTH_$1, String
0630: .valueOf(longSt.countTokens())));
0631: }
0632:
0633: //check that latitude grid header is the same as for latitude grid
0634: if ((nc != Integer.parseInt(longSt.nextToken()))
0635: || (nr != Integer.parseInt(longSt.nextToken()))
0636: || (nz != Integer.parseInt(longSt.nextToken()))
0637: || (xmin != Float.parseFloat(longSt.nextToken()))
0638: || (dx != Float.parseFloat(longSt.nextToken()))
0639: || (ymin != Float.parseFloat(longSt.nextToken()))
0640: || (dy != Float.parseFloat(longSt.nextToken()))
0641: || (angle != Float.parseFloat(longSt.nextToken()))) {
0642: throw new FactoryException(Errors
0643: .format(ErrorKeys.GRID_LOCATIONS_UNEQUAL));
0644: }
0645:
0646: ////////////////////////
0647: //read grid shift data into LocalizationGrid
0648: ////////////////////////
0649: gridShift = new LocalizationGrid(nc, nr);
0650:
0651: int i = 0;
0652: int j = 0;
0653: for (i = 0; i < nr; i++) {
0654: for (j = 0; j < nc;) {
0655: latLine = latBr.readLine();
0656: latSt = new StringTokenizer(latLine, " ");
0657: longLine = longBr.readLine();
0658: longSt = new StringTokenizer(longLine, " ");
0659:
0660: while (latSt.hasMoreTokens() && longSt.hasMoreTokens()) {
0661: gridShift.setLocalizationPoint(j, i, (double) Float
0662: .parseFloat(longSt.nextToken()),
0663: (double) Float
0664: .parseFloat(latSt.nextToken()));
0665: ++j;
0666: }
0667: }
0668: }
0669:
0670: assert i == nr : i;
0671: assert j == nc : j;
0672: }
0673:
0674: /**
0675: * Transforms a list of coordinate point ordinal values. This method is
0676: * provided for efficiently transforming many points. The supplied array
0677: * of ordinal values will contain packed ordinal values. For example, if
0678: * the source dimension is 3, then the ordinals will be packed in this
0679: * order:
0680: * (<var>x<sub>0</sub></var>,<var>y<sub>0</sub></var>,<var>z<sub>0</sub></var>,
0681: *
0682: * <var>x<sub>1</sub></var>,<var>y<sub>1</sub></var>,<var>z<sub>1</sub></var>
0683: * ...). All input and output values are in decimal degrees.
0684: *
0685: * @param srcPts the array containing the source point coordinates.
0686: * @param srcOff the offset to the first point to be transformed in the
0687: * source array.
0688: * @param dstPts the array into which the transformed point coordinates are
0689: * returned. May be the same than {@code srcPts}.
0690: * @param dstOff the offset to the location of the first transformed point
0691: * that is stored in the destination array.
0692: * @param numPts the number of point objects to be transformed.
0693: *
0694: * @throws TransformException if the input point is outside the area
0695: * covered by this grid.
0696: */
0697: public void transform(final double[] srcPts, int srcOff,
0698: final double[] dstPts, int dstOff, int numPts)
0699: throws TransformException {
0700: int step = 0;
0701:
0702: if ((srcPts == dstPts)
0703: && (srcOff < dstOff)
0704: && ((srcOff + (numPts * getSourceDimensions())) > dstOff)) {
0705: step = -getSourceDimensions();
0706: srcOff -= ((numPts - 1) * step);
0707: dstOff -= ((numPts - 1) * step);
0708: }
0709:
0710: while (--numPts >= 0) {
0711: double x = srcPts[srcOff++];
0712: double y = srcPts[srcOff++];
0713:
0714: //check bounding box
0715: //issue of bbox crossing +- 180 degrees (ie input point of -188 longitude);
0716: //abs(x - xmin) > 0 , rollLongitude() ???
0717: if (((x < xmin) || (x > xmax))
0718: || ((y < ymin) || (y > ymax))) {
0719: throw new TransformException(Errors
0720: .format(ErrorKeys.POINT_OUTSIDE_GRID));
0721: }
0722:
0723: //find the grid the point is in (index is 0 based)
0724: final double xgrid = (x - xmin) / dx;
0725: final double ygrid = (y - ymin) / dy;
0726: double[] array = new double[] { xgrid, ygrid };
0727:
0728: //use the LocalizationGridTransform2D transform method (bilineal interpolation)
0729: //returned shift values are in seconds, longitude shift values are + west
0730: gridShiftTransform.transform(array, 0, array, 0, 1);
0731:
0732: dstPts[dstOff++] = x - (array[0] / SEC_2_DEG);
0733: dstPts[dstOff++] = y + (array[1] / SEC_2_DEG);
0734: srcOff += step;
0735: dstOff += step;
0736: }
0737: }
0738:
0739: /**
0740: * Transforms nad83 values to nad27. Input and output values are in
0741: * decimal degrees. This is done by itteratively finding a nad27 value that
0742: * shifts to the input nad83 value. The input nad83 value is used as the
0743: * first approximation.
0744: *
0745: * @param srcPts the array containing the source point coordinates.
0746: * @param srcOff the offset to the first point to be transformed in the
0747: * source array.
0748: * @param dstPts the array into which the transformed point coordinates are
0749: * returned. May be the same than {@code srcPts}.
0750: * @param dstOff the offset to the location of the first transformed point
0751: * that is stored in the destination array.
0752: * @param numPts the number of point objects to be transformed.
0753: *
0754: * @throws TransformException if the input point is outside the area
0755: * covered by this grid.
0756: */
0757: public void inverseTransform(final double[] srcPts, int srcOff,
0758: final double[] dstPts, int dstOff, int numPts)
0759: throws TransformException {
0760: int step = 0;
0761:
0762: if ((srcPts == dstPts)
0763: && (srcOff < dstOff)
0764: && ((srcOff + (numPts * getSourceDimensions())) > dstOff)) {
0765: step = -getSourceDimensions();
0766: srcOff -= ((numPts - 1) * step);
0767: dstOff -= ((numPts - 1) * step);
0768: }
0769:
0770: while (--numPts >= 0) {
0771: final double x = srcPts[srcOff++];
0772: final double y = srcPts[srcOff++];
0773: double xtemp = x;
0774: double ytemp = y;
0775:
0776: for (int i = MAX_ITER;;) {
0777: double[] array = { xtemp, ytemp };
0778: transform(array, 0, array, 0, 1);
0779: double xdif = array[0] - x;
0780: double ydif = array[1] - y;
0781:
0782: if (Math.abs(xdif) > TOL) {
0783: xtemp = xtemp - xdif;
0784: }
0785: if (Math.abs(ydif) > TOL) {
0786: ytemp = ytemp - ydif;
0787: }
0788:
0789: if ((Math.abs(xdif) <= TOL) && (Math.abs(ydif) <= TOL)) {
0790: dstPts[dstOff++] = xtemp;
0791: dstPts[dstOff++] = ytemp;
0792: break;
0793: }
0794: if (--i < 0) {
0795: throw new TransformException(Errors
0796: .format(ErrorKeys.NO_CONVERGENCE));
0797: }
0798: }
0799:
0800: srcOff += step;
0801: dstOff += step;
0802: }
0803: }
0804:
0805: /**
0806: * Returns the inverse of this transform.
0807: *
0808: * @return the inverse of this transform
0809: */
0810: public MathTransform inverse() {
0811: if (inverse == null) {
0812: // No need to synchronize; this is not a big deal if this object is created twice.
0813: inverse = new Inverse();
0814: }
0815:
0816: return inverse;
0817: }
0818:
0819: /**
0820: * Returns a hash value for this transform. To make this faster it does not
0821: * check the grid values.
0822: *
0823: * @return a hash value for this transform.
0824: */
0825: public final int hashCode() {
0826: final long code = Double.doubleToLongBits(xmin)
0827: + (37 * (Double.doubleToLongBits(ymin) + (37 * (Double
0828: .doubleToLongBits(xmax) + (37 * (Double
0829: .doubleToLongBits(ymax) + (37 * (Double
0830: .doubleToLongBits(dx) + (37 * (Double
0831: .doubleToLongBits(dy)))))))))));
0832:
0833: return (int) code ^ (int) (code >>> 32);
0834: }
0835:
0836: /**
0837: * Compares the specified object with this math transform for equality.
0838: *
0839: * @param object the object to compare to
0840: * @return {@code true} if the objects are equal.
0841: */
0842: public final boolean equals(final Object object) {
0843: if (object == this ) {
0844: // Slight optimization
0845: return true;
0846: }
0847:
0848: if (super .equals(object)) {
0849: final NADCONTransform that = (NADCONTransform) object;
0850:
0851: return (Double.doubleToLongBits(this .xmin) == Double
0852: .doubleToLongBits(that.xmin))
0853: && (Double.doubleToLongBits(this .ymin) == Double
0854: .doubleToLongBits(that.ymin))
0855: && (Double.doubleToLongBits(this .xmax) == Double
0856: .doubleToLongBits(that.xmax))
0857: && (Double.doubleToLongBits(this .ymax) == Double
0858: .doubleToLongBits(that.ymax))
0859: && (Double.doubleToLongBits(this .dx) == Double
0860: .doubleToLongBits(that.dx))
0861: && (Double.doubleToLongBits(this .dy) == Double
0862: .doubleToLongBits(that.dy))
0863: && (this .gridShiftTransform)
0864: .equals(that.gridShiftTransform);
0865: }
0866:
0867: return false;
0868: }
0869:
0870: /**
0871: * Used to set the preference for the default grid shift file location.
0872: * This allows grids parameters to be specified by name only, without the
0873: * full path. This needs to be done only once, by the user.
0874: * Path values may be simple file system paths or more complex
0875: * text representations of a url. A value of "default" resets this
0876: * preference to its default value.
0877: * <p>
0878: *
0879: * Example:
0880: * <blockquote>
0881: * <pre>
0882: * java org.geotools.referencing.operation.transform.NADCONTransform file:///home/rschulz/GIS/NADCON/data
0883: * </pre>
0884: * </blockquote>
0885: *
0886: * @param args a single argument for the defualt location of grid shift
0887: * files
0888: */
0889: public static void main(String[] args) {
0890: final Arguments arguments = new Arguments(args);
0891: final PrintWriter out = arguments.out;
0892: final Preferences prefs = Preferences
0893: .userNodeForPackage(NADCONTransform.class);
0894:
0895: if (args.length == 1) {
0896: if (args[0].equalsIgnoreCase("default")) {
0897: prefs.remove(GRID_LOCATION);
0898: } else {
0899: prefs.put(GRID_LOCATION, args[0]);
0900: }
0901:
0902: return;
0903: } else {
0904: final String location = prefs.get(GRID_LOCATION,
0905: DEFAULT_GRID_LOCATION);
0906: out
0907: .println("Usage: java org.geotools.referencing.operation.transform.NADCONTransform "
0908: + "<defalult grid file location (path)>");
0909: out.print("Grid location: \"");
0910: out.print(location);
0911: out.println('"');
0912:
0913: return;
0914: }
0915: }
0916:
0917: /**
0918: * Inverse of a {@link NADCONTransform}.
0919: *
0920: * @version $Id: NADCONTransform.java 29101 2008-02-06 12:11:19Z desruisseaux $
0921: * @author Rueben Schulz
0922: */
0923: private final class Inverse extends AbstractMathTransform.Inverse
0924: implements Serializable {
0925: /** Serial number for interoperability with different versions. */
0926: private static final long serialVersionUID = -4707304160205218546L;
0927:
0928: /**
0929: * Default constructor.
0930: */
0931: public Inverse() {
0932: NADCONTransform.this .super ();
0933: }
0934:
0935: /**
0936: * Returns the parameter values for this math transform.
0937: *
0938: * @return A copy of the parameter values for this math transform.
0939: */
0940: public ParameterValueGroup getParameterValues() {
0941: return null;
0942: }
0943:
0944: /**
0945: * Inverse transform an array of points.
0946: *
0947: * @param source
0948: * @param srcOffset
0949: * @param dest
0950: * @param dstOffset
0951: * @param length
0952: *
0953: * @throws TransformException if the input point is outside the area
0954: * covered by this grid.
0955: */
0956: public void transform(final double[] source,
0957: final int srcOffset, final double[] dest,
0958: final int dstOffset, final int length)
0959: throws TransformException {
0960: NADCONTransform.this .inverseTransform(source, srcOffset,
0961: dest, dstOffset, length);
0962: }
0963:
0964: /**
0965: * Restore reference to this object after deserialization.
0966: *
0967: * @param in DOCUMENT ME!
0968: * @throws IOException DOCUMENT ME!
0969: * @throws ClassNotFoundException DOCUMENT ME!
0970: */
0971: private void readObject(ObjectInputStream in)
0972: throws IOException, ClassNotFoundException {
0973: in.defaultReadObject();
0974: NADCONTransform.this .inverse = this ;
0975: }
0976: }
0977:
0978: /**
0979: * The provider for {@link NADCONTransform}. This provider will construct
0980: * transforms from {@linkplain org.geotools.referencing.crs.DefaultGeographicCRS
0981: * geographic} to {@linkplain org.geotools.referencing.crs.DefaultGeographicCRS
0982: * geographic} coordinate reference systems.
0983: *
0984: * @version $Id: NADCONTransform.java 29101 2008-02-06 12:11:19Z desruisseaux $
0985: * @author Rueben Schulz
0986: */
0987: public static class Provider extends MathTransformProvider {
0988: /** Serial number for interoperability with different versions. */
0989: private static final long serialVersionUID = -4707304160205218546L;
0990:
0991: /**
0992: * The operation parameter descriptor for the "Latitude_difference_file"
0993: * parameter value. The default value is "conus.las".
0994: */
0995: public static final ParameterDescriptor LAT_DIFF_FILE = new DefaultParameterDescriptor(
0996: "Latitude_difference_file", String.class, null,
0997: "conus.las");
0998:
0999: /**
1000: * The operation parameter descriptor for the "Longitude_difference_file"
1001: * parameter value. The default value is "conus.los".
1002: */
1003: public static final ParameterDescriptor LONG_DIFF_FILE = new DefaultParameterDescriptor(
1004: "Longitude_difference_file", String.class, null,
1005: "conus.los");
1006:
1007: /**
1008: * The parameters group.
1009: */
1010: static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(
1011: new NamedIdentifier[] {
1012: new NamedIdentifier(Citations.OGC, "NADCON"),
1013: new NamedIdentifier(Citations.EPSG, "NADCON"),
1014: new NamedIdentifier(Citations.EPSG, "9613"),
1015: new NamedIdentifier(
1016: Citations.GEOTOOLS,
1017: Vocabulary
1018: .formatInternational(VocabularyKeys.NADCON_TRANSFORM)) },
1019: new ParameterDescriptor[] { LAT_DIFF_FILE,
1020: LONG_DIFF_FILE });
1021:
1022: /**
1023: * Constructs a provider.
1024: */
1025: public Provider() {
1026: super (2, 2, PARAMETERS);
1027: }
1028:
1029: /**
1030: * Returns the operation type.
1031: */
1032: public Class getOperationType() {
1033: return Transformation.class;
1034: }
1035:
1036: /**
1037: * Creates a math transform from the specified group of parameter
1038: * values.
1039: *
1040: * @param values The group of parameter values.
1041: * @return The created math transform.
1042: * @throws ParameterNotFoundException if a required parameter was not
1043: * found.
1044: * @throws FactoryException if there is a problem creating this
1045: * math transform.
1046: */
1047: protected MathTransform createMathTransform(
1048: final ParameterValueGroup values)
1049: throws ParameterNotFoundException, FactoryException {
1050: return new NADCONTransform(stringValue(LAT_DIFF_FILE,
1051: values), stringValue(LONG_DIFF_FILE, values));
1052: }
1053: }
1054: }
|