001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: * (C) 2006, Geotools Project Managment Committee (PMC)
005: *
006: * This library is free software; you can redistribute it and/or
007: * modify it under the terms of the GNU Lesser General Public
008: * License as published by the Free Software Foundation; either
009: * version 2.1 of the License, or (at your option) any later version.
010: *
011: * This library is distributed in the hope that it will be useful,
012: * but WITHOUT ANY WARRANTY; without even the implied warranty of
013: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
014: * Lesser General Public License for more details.
015: *
016: * This file is derived from NGA/NASA software available for unlimited distribution.
017: * See http://earth-info.nima.mil/GandG/wgs84/gravitymod/.
018: */
019: package org.geotools.referencing.operation.transform;
020:
021: // J2SE dependencies
022: import java.io.IOException;
023: import java.io.InputStream;
024: import java.io.InputStreamReader;
025: import java.io.LineNumberReader;
026: import java.io.FileNotFoundException;
027: import java.util.Collections;
028: import java.util.StringTokenizer;
029:
030: // OpenGIS dependencies
031: import org.opengis.parameter.ParameterValue;
032: import org.opengis.parameter.ParameterValueGroup;
033: import org.opengis.parameter.ParameterDescriptor;
034: import org.opengis.parameter.ParameterDescriptorGroup;
035: import org.opengis.parameter.ParameterNotFoundException;
036: import org.opengis.referencing.FactoryException;
037: import org.opengis.referencing.operation.MathTransform;
038: import org.opengis.referencing.operation.Transformation;
039: import org.opengis.referencing.operation.TransformException;
040:
041: // Geotools dependencies
042: import org.geotools.parameter.Parameter;
043: import org.geotools.parameter.ParameterGroup;
044: import org.geotools.parameter.DefaultParameterDescriptor;
045: import org.geotools.referencing.NamedIdentifier;
046: import org.geotools.referencing.operation.MathTransformProvider;
047: import org.geotools.resources.i18n.Errors;
048: import org.geotools.resources.i18n.ErrorKeys;
049: import org.geotools.resources.i18n.Vocabulary;
050: import org.geotools.resources.i18n.VocabularyKeys;
051: import org.geotools.metadata.iso.citation.Citations;
052:
053: /**
054: * Transforms vertical coordinates using coefficients from the
055: * <A HREF="http://earth-info.nima.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html">Earth
056: * Gravitational Model</A>.
057: * <p>
058: * <strong>Aknowledgement</strong><br>
059: * This class is an adaption of Fortran code
060: * <code><a href="http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/clenqt.for">clenqt.for</a></code>
061: * from the <cite>National Geospatial-Intelligence Agency</cite> and available in public domain. The
062: * <cite>normalized geopotential coefficients</cite> file bundled in this module is an adaptation of
063: * <code><a href="http://earth-info.nima.mil/GandG/wgs84/gravitymod/wgs84_180/egm180.nor">egm180.nor</a></code>
064: * file, with some spaces trimmed.
065: *
066: * @since 2.3
067: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/plugin/referencing3D/src/main/java/org/geotools/referencing/operation/transform/EarthGravitationalModel.java $
068: * @version $Id: EarthGravitationalModel.java 24384 2007-02-14 00:23:05Z desruisseaux $
069: * @author Pierre Cardinal
070: * @author Martin Desruisseaux
071: */
072: public final class EarthGravitationalModel extends VerticalTransform {
073: /**
074: * Pre-computed values of some square roots.
075: */
076: private static final double SQRT_03 = 1.7320508075688772935274463415059,
077: SQRT_05 = 2.2360679774997896964091736687313,
078: SQRT_13 = 3.6055512754639892931192212674705,
079: SQRT_17 = 4.1231056256176605498214098559741,
080: SQRT_21 = 4.5825756949558400065880471937280;
081:
082: /** The default value for {@link #nmax}. */
083: static final int DEFAULT_ORDER = 180;
084:
085: /** {@code true} for WGS84 model, or {@code false} for WGS72 model. */
086: private final boolean wgs84;
087:
088: /** Maximum degree and order attained. */
089: private final int nmax;
090:
091: /** WGS 84 semi-major axis. */
092: private final double semiMajor;
093:
094: /** The first Eccentricity Squared (e²) for WGS 84 ellipsoid. */
095: private final double esq;
096:
097: /** Even zonal coefficient. */
098: private final double c2;
099:
100: /** WGS 84 Earth's Gravitational Constant w/ atmosphere. */
101: private final double rkm;
102:
103: /** Theoretical (Normal) Gravity at the Equator (on the Ellipsoid). */
104: private final double grava;
105:
106: /** Theoretical (Normal) Gravity Formula Constant. */
107: private final double star;
108:
109: /**
110: * The geopotential coefficients read from the ASCII file.
111: * Those arrays are filled by the {@link #load} method.
112: */
113: private final double[] cnmGeopCoef, snmGeopCoef;
114:
115: /**
116: * Cleanshaw coefficients needed for the selected gravimetric quantities that are computed.
117: * Those arrays are computed by the {@link #initialize} method.
118: */
119: private final double[] aClenshaw, bClenshaw, as;
120:
121: /**
122: * Temporary buffer for use by {@link #heightOffset} only. Allocated once for ever
123: * for avoiding too many objects creation / destruction.
124: */
125: private final double[] cr, sr, s11, s12;
126:
127: /**
128: * Creates a model with the default maximum degree and order.
129: */
130: EarthGravitationalModel() {
131: this (DEFAULT_ORDER, true);
132: }
133:
134: /**
135: * Creates a model with the specified maximum degree and order.
136: */
137: EarthGravitationalModel(final int nmax, final boolean wgs84) {
138: this .nmax = nmax;
139: this .wgs84 = wgs84;
140: if (wgs84) {
141: /*
142: * WGS84 model values.
143: * NOTE: The Fortran program gives 3.9860015e+14 for 'rkm' constant. This value has been
144: * modified in later programs. From http://cddis.gsfc.nasa.gov/926/egm96/doc/S11.HTML :
145: *
146: * "We next need to consider the determination of GM, GM0, W0, U0. The value of GM0
147: * will be that adopted for the updated GM of the WGS84 ellipsoid. This value is
148: * 3.986004418e+14 m³/s², which is identical to that given in the IERS Numerical
149: * Standards [McCarthy, 1996, Table 4.1]. The best estimate of GM can be taken as
150: * the same value based on the recommendations of the IAG Special Commission SC3,
151: * Fundamental Constants [Bursa, 1995b, p. 381]."
152: */
153: semiMajor = 6378137.0;
154: esq = 0.00669437999013;
155: c2 = 108262.9989050e-8;
156: rkm = 3.986004418e+14;
157: grava = 9.7803267714;
158: star = 0.001931851386;
159: } else {
160: /*
161: * WGS72 model values.
162: */
163: semiMajor = 6378135.0;
164: esq = 0.006694317778;
165: c2 = 108263.0e-8;
166: rkm = 3.986005e+14;
167: grava = 9.7803327;
168: star = 0.005278994;
169: }
170: final int cleanshawLength = locatingArray(nmax + 3);
171: final int geopCoefLength = locatingArray(nmax + 1);
172: aClenshaw = new double[cleanshawLength];
173: bClenshaw = new double[cleanshawLength];
174: cnmGeopCoef = new double[geopCoefLength];
175: snmGeopCoef = new double[geopCoefLength];
176: as = new double[nmax + 1];
177: cr = new double[nmax + 1];
178: sr = new double[nmax + 1];
179: s11 = new double[nmax + 3];
180: s12 = new double[nmax + 3];
181: }
182:
183: /**
184: * Computes the index as it would be returned by the locating array {@code iv}
185: * (from the Fortran code).
186: * <p>
187: * Tip (used in some place in this class):
188: * {@code locatingArray(n+1)} == {@code locatingArray(n) + n + 1}.
189: */
190: private static int locatingArray(final int n) {
191: return ((n + 1) * n) >> 1;
192: }
193:
194: /**
195: * Loads the coefficients from the specified ASCII file and initialize the internal
196: * <cite>clenshaw arrays</cite>.
197: * <p>
198: * <strong>Note:</strong> ASCII may looks like an unefficient format for binary distribution.
199: * A binary file with coefficient values read by {@link java.io.DataInput#readDouble} would
200: * be more compact than an <u>uncompressed</u> ASCII file. However, binary files are hard to
201: * compress by the ZIP algorithm. Our experience show that a 675 kb uncompressed ASCII file
202: * is only 222 kb after ZIP or JAR compression. The same data as a binary file is 257 kb
203: * uncompressed and 248 kb compressed. So surprisingly, the ASCII file is more compact than
204: * the binary file after compression. Since it is the primary format provided by the
205: * Earth-Info web site, we use it directly in order to avoid a multiplication of formats.
206: *
207: * @param filename The filename (e.g. {@code "WGS84.cof"}, relative to this class directory.
208: * @throws IOException if the file can't be read or has an invalid content.
209: */
210: protected void load(final String filename) throws IOException {
211: final InputStream stream = EarthGravitationalModel.class
212: .getResourceAsStream(filename);
213: if (stream == null) {
214: throw new FileNotFoundException(filename);
215: }
216: final LineNumberReader in;
217: in = new LineNumberReader(new InputStreamReader(stream,
218: "ISO-8859-1"));
219: String line;
220: while ((line = in.readLine()) != null) {
221: final StringTokenizer tokens = new StringTokenizer(line);
222: try {
223: /*
224: * Note: we use 'parseShort' instead of 'parseInt' as an easy way to ensure that
225: * the values are in some reasonable range. The range is typically [0..180].
226: * We don't check that, but at least 'parseShort' disallows values greater
227: * than 32767. Additional note: we real all lines in all cases even if we
228: * discard some of them, in order to check the file format.
229: */
230: final int n = Short.parseShort(tokens.nextToken());
231: final int m = Short.parseShort(tokens.nextToken());
232: final double cbar = Double.parseDouble(tokens
233: .nextToken());
234: final double sbar = Double.parseDouble(tokens
235: .nextToken());
236: if (n <= nmax) {
237: final int ll = locatingArray(n) + m;
238: cnmGeopCoef[ll] = cbar;
239: snmGeopCoef[ll] = sbar;
240: }
241: } catch (RuntimeException cause) {
242: /*
243: * Catch the following exceptions:
244: * - NoSuchElementException if a line has too few numbers.
245: * - NumberFormatException if a number can't be parsed.
246: * - IndexOutOfBoundsException if 'n' or 'm' values are illegal.
247: */
248: final IOException exception = new IOException(Errors
249: .format(ErrorKeys.BAD_LINE_IN_FILE_$2,
250: filename, new Integer(in
251: .getLineNumber())));
252: exception.initCause(cause);
253: throw exception;
254: }
255: }
256: in.close();
257: initialize();
258: }
259:
260: /**
261: * Computes the <cite>clenshaw arrays</cite> after all coefficients have been read.
262: * We performs this step in a separated method than {@link #from} in case we wish
263: * to read the coefficient from an other source than an ASCII file in some future
264: * version.
265: */
266: private final void initialize() {
267: /*
268: * MODIFY CNM EVEN ZONAL COEFFICIENTS.
269: */
270: if (wgs84) {
271: final double[] c2n = new double[6];
272: c2n[1] = c2;
273: int sign = 1;
274: double esqi = esq;
275: for (int i = 2; i < c2n.length; i++) {
276: sign *= -1;
277: esqi *= esq;
278: c2n[i] = sign * (3 * esqi)
279: / ((2 * i + 1) * (2 * i + 3))
280: * (1 - i + (5 * i * c2 / esq));
281: }
282: /* all nmax */cnmGeopCoef[3] += c2n[1] / SQRT_05;
283: /* all nmax */cnmGeopCoef[10] += c2n[2] / 3;
284: /* all nmax */cnmGeopCoef[21] += c2n[3] / SQRT_13;
285: if (nmax > 6)
286: cnmGeopCoef[36] += c2n[4] / SQRT_17;
287: if (nmax > 9)
288: cnmGeopCoef[55] += c2n[5] / SQRT_21;
289: } else {
290: /* all nmax */cnmGeopCoef[3] += 4.841732e-04;
291: /* all nmax */cnmGeopCoef[10] += -7.8305e-07;
292: }
293: /*
294: * BUILD ALL CLENSHAW COEFFICIENT ARRAYS.
295: */
296: for (int i = 0; i <= nmax; i++) {
297: as[i] = -Math.sqrt(1.0 + 1.0 / (2 * (i + 1)));
298: }
299: for (int i = 0; i <= nmax; i++) {
300: for (int j = i + 1; j <= nmax; j++) {
301: final int ll = locatingArray(j) + i;
302: final int n = 2 * j + 1;
303: final int ji = (j - i) * (j + i);
304: aClenshaw[ll] = Math.sqrt(n * (2 * j - 1)
305: / (double) (ji));
306: bClenshaw[ll] = Math.sqrt(n * (j + i - 1) * (j - i - 1)
307: / (double) (ji * (2 * j - 3)));
308: }
309: }
310: }
311:
312: /**
313: * {@inheritDoc}
314: */
315: public double heightOffset(final double longitude,
316: final double latitude, final double height)
317: throws TransformException {
318: /*
319: * Note: no need to ensure that longitude is in [-180..+180°] range, because its value
320: * is used only in trigonometric functions (sin / cos), which roll it as we would expect.
321: * Latitude is used only in trigonometric functions as well.
322: */
323: final double phi = Math.toRadians(latitude);
324: final double sin_phi = Math.sin(phi);
325: final double sin2_phi = sin_phi * sin_phi;
326: final double rni = Math.sqrt(1.0 - esq * sin2_phi);
327: final double rn = semiMajor / rni;
328: final double t22 = (rn + height) * Math.cos(phi);
329: final double x2y2 = t22 * t22;
330: final double z1 = ((rn * (1 - esq)) + height) * sin_phi;
331: final double th = (Math.PI / 2.0)
332: - Math.atan(z1 / Math.sqrt(x2y2));
333: final double y = Math.sin(th);
334: final double t = Math.cos(th);
335: final double f1 = semiMajor / Math.sqrt(x2y2 + z1 * z1);
336: final double f2 = f1 * f1;
337: final double rlam = Math.toRadians(longitude);
338: final double gravn;
339: if (wgs84) {
340: gravn = grava * (1.0 + star * sin2_phi) / rni;
341: } else {
342: gravn = grava * (1.0 + star * sin2_phi) + 0.000023461
343: * (sin2_phi * sin2_phi);
344: }
345: sr[0] = 0;
346: sr[1] = Math.sin(rlam);
347: cr[0] = 1;
348: cr[1] = Math.cos(rlam);
349: for (int j = 2; j <= nmax; j++) {
350: sr[j] = (2.0 * cr[1] * sr[j - 1]) - sr[j - 2];
351: cr[j] = (2.0 * cr[1] * cr[j - 1]) - cr[j - 2];
352: }
353: double sht = 0, previousSht = 0;
354: for (int i = nmax; i >= 0; i--) {
355: for (int j = nmax; j >= i; j--) {
356: final int ll = locatingArray(j) + i;
357: final int ll2 = ll + j + 1;
358: final int ll3 = ll2 + j + 2;
359: final double ta = aClenshaw[ll2] * f1 * t;
360: final double tb = bClenshaw[ll3] * f2;
361: s11[j] = (ta * s11[j + 1]) - (tb * s11[j + 2])
362: + cnmGeopCoef[ll];
363: s12[j] = (ta * s12[j + 1]) - (tb * s12[j + 2])
364: + snmGeopCoef[ll];
365: }
366: previousSht = sht;
367: sht = (-as[i] * y * f1 * sht) + (s11[i] * cr[i])
368: + (s12[i] * sr[i]);
369: }
370: return ((s11[0] + s12[0]) * f1 + (previousSht * SQRT_03 * y * f2))
371: * rkm / (semiMajor * (gravn - (height * 0.3086e-5)));
372: }
373:
374: /**
375: * Returns the parameter descriptors for this math transform.
376: */
377: public ParameterDescriptorGroup getParameterDescriptors() {
378: return Provider.PARAMETERS;
379: }
380:
381: /**
382: * Returns the parameters for this math transform.
383: */
384: public ParameterValueGroup getParameterValues() {
385: return new ParameterGroup(getParameterDescriptors(),
386: new ParameterValue[] { new Parameter(Provider.ORDER,
387: new Integer(nmax)) });
388: }
389:
390: /**
391: * The provider for {@link EarthGravitationalModel}.
392: *
393: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/plugin/referencing3D/src/main/java/org/geotools/referencing/operation/transform/EarthGravitationalModel.java $
394: * @version $Id: EarthGravitationalModel.java 24384 2007-02-14 00:23:05Z desruisseaux $
395: * @author Martin Desruisseaux
396: */
397: public static class Provider extends MathTransformProvider {
398: /**
399: * The operation parameter descriptor for the maximum degree and order.
400: * The default value is 180.
401: */
402: public static final ParameterDescriptor ORDER = new DefaultParameterDescriptor(
403: Collections
404: .singletonMap(
405: NAME_KEY,
406: new NamedIdentifier(
407: Citations.GEOTOOLS,
408: Vocabulary
409: .formatInternational(VocabularyKeys.ORDER))),
410: DEFAULT_ORDER, 2, 180, false);
411:
412: /**
413: * The parameters group.
414: */
415: static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(
416: new NamedIdentifier[] { new NamedIdentifier(
417: Citations.GEOTOOLS,
418: Vocabulary
419: .formatInternational(VocabularyKeys.EARTH_GRAVITATIONAL_MODEL)) },
420: new ParameterDescriptor[] { ORDER });
421:
422: /**
423: * Constructs a math transform provider.
424: */
425: public Provider() {
426: super (3, 3, PARAMETERS);
427: }
428:
429: /**
430: * Returns the operation type for this transform.
431: */
432: public Class getOperationType() {
433: return Transformation.class;
434: }
435:
436: /**
437: * Creates a math transform from the specified group of parameter values.
438: *
439: * @param values The group of parameter values.
440: * @return The created math transform.
441: * @throws ParameterNotFoundException if a required parameter was not found.
442: * @throws FactoryException if this method failed to load the coefficient file.
443: */
444: protected MathTransform createMathTransform(
445: final ParameterValueGroup values)
446: throws ParameterNotFoundException, FactoryException {
447: int nmax = intValue(ORDER, values);
448: if (nmax == 0) {
449: nmax = DEFAULT_ORDER;
450: }
451: final EarthGravitationalModel mt = new EarthGravitationalModel(
452: nmax, true);
453: final String filename = "EGM180.nor";
454: try {
455: mt.load(filename);
456: } catch (IOException e) {
457: throw new FactoryException(Errors.format(
458: ErrorKeys.CANT_READ_$1, filename), e);
459: }
460: return mt;
461: }
462: }
463: }
|