001: //$HeadURL: $
002: /*---------------- FILE HEADER ------------------------------------------
003: This file is part of deegree.
004: Copyright (C) 2001-2008 by:
005: Department of Geography, University of Bonn
006: http://www.giub.uni-bonn.de/deegree/
007: lat/lon GmbH
008: http://www.lat-lon.de
009:
010: This library is free software; you can redistribute it and/or
011: modify it under the terms of the GNU Lesser General Public
012: License as published by the Free Software Foundation; either
013: version 2.1 of the License, or (at your option) any later version.
014: This library is distributed in the hope that it will be useful,
015: but WITHOUT ANY WARRANTY; without even the implied warranty of
016: MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
017: Lesser General Public License for more details.
018: You should have received a copy of the GNU Lesser General Public
019: License along with this library; if not, write to the Free Software
020: Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
021: Contact:
022:
023: Andreas Poth
024: lat/lon GmbH
025: Aennchenstr. 19
026: 53177 Bonn
027: Germany
028: E-Mail: poth@lat-lon.de
029:
030: Prof. Dr. Klaus Greve
031: Department of Geography
032: University of Bonn
033: Meckenheimer Allee 166
034: 53115 Bonn
035: Germany
036: E-Mail: greve@giub.uni-bonn.de
037: ---------------------------------------------------------------------------*/
038:
039: package org.deegree.crs.transformations;
040:
041: import java.util.ArrayList;
042: import java.util.List;
043:
044: import javax.vecmath.Point3d;
045:
046: import org.deegree.crs.components.Ellipsoid;
047: import org.deegree.crs.components.Unit;
048: import org.deegree.crs.coordinatesystems.GeocentricCRS;
049: import org.deegree.crs.coordinatesystems.GeographicCRS;
050: import org.deegree.crs.projections.ProjectionUtils;
051: import org.deegree.framework.log.ILogger;
052: import org.deegree.framework.log.LoggerFactory;
053:
054: /**
055: * The <code>GeocentricTransform</code> class is used to create a transformation between a geocentric CRS (having
056: * lat-lon coordinates) and it's geodetic CRS (having x-y-z) coordinates and vice versa.
057: *
058: * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
059: *
060: * @author last edited by: $Author:$
061: *
062: * @version $Revision:$, $Date:$
063: *
064: */
065:
066: public class GeocentricTransform extends CRSTransformation {
067:
068: private static ILogger LOG = LoggerFactory
069: .getLogger(GeocentricTransform.class);
070:
071: /**
072: * Cosine of 67.5 degrees.
073: */
074: private static final double COS_67P5 = 0.38268343236508977;
075:
076: /**
077: * Toms region 1 constant, which will allow for a difference between ellipsoid of 2000 Kilometers. <quote>Under this
078: * policy the maximum error is less than 42 centimeters for altitudes less then 10.000.000 Kilometers</quote>
079: */
080: private static final double AD_C = 1.0026;
081:
082: /**
083: * Semi-major axis of ellipsoid in meters.
084: */
085: private final double semiMajorAxis;
086:
087: /**
088: * Semi-minor axis of ellipsoid in meters.
089: */
090: private final double semiMinorAxis;
091:
092: /**
093: * Square of semi-major axis {@link #semiMajorAxis}.
094: */
095: private final double squaredSemiMajorAxis;
096:
097: /**
098: * Square of semi-minor axis {@link #semiMinorAxis}.
099: */
100: private final double squaredSemiMinorAxis;
101:
102: /**
103: * Eccentricity squared.
104: */
105: private final double squaredEccentricity;
106:
107: /**
108: * 2nd eccentricity squared.
109: */
110: private final double ep2;
111:
112: /**
113: * true if the given points will use heights (e.g. have a z/height-value).
114: */
115: private boolean hasHeight;
116:
117: /**
118: * @param source
119: * the geographic crs.
120: * @param target
121: * the geocentric crs.
122: * @param hasHeight a flag signaling the inverse transformation to calculate a height value too.
123: */
124: public GeocentricTransform(GeographicCRS source,
125: GeocentricCRS target, boolean hasHeight) {
126: super (source, target);
127: this .hasHeight = hasHeight;
128: Ellipsoid ellipsoid = source.getGeodeticDatum().getEllipsoid();
129: semiMajorAxis = Unit.METRE.convert(
130: ellipsoid.getSemiMajorAxis(), ellipsoid.getUnits());
131: semiMinorAxis = Unit.METRE.convert(
132: ellipsoid.getSemiMinorAxis(), ellipsoid.getUnits());
133: squaredSemiMajorAxis = semiMajorAxis * semiMajorAxis;
134: squaredSemiMinorAxis = semiMinorAxis * semiMinorAxis;
135: squaredEccentricity = ellipsoid.getSquaredEccentricity();
136: // e2 = ( a2 - b2 ) / a2;
137: ep2 = (squaredSemiMajorAxis - squaredSemiMinorAxis)
138: / squaredSemiMinorAxis;
139: }
140:
141: @Override
142: public List<Point3d> doTransform(List<Point3d> srcPts) {
143: List<Point3d> result = new ArrayList<Point3d>(srcPts);
144: if (ILogger.LOG_DEBUG == LOG.getLevel()) {
145: StringBuilder sb = new StringBuilder(
146: isInverseTransform() ? "An inverse" : "A");
147: sb.append(getName());
148: sb.append(" with incoming points: ");
149: sb.append(srcPts);
150: LOG.logDebug(sb.toString());
151: }
152: if (isInverseTransform()) {
153: toGeographic(result);
154: } else {
155: toGeoCentric(result);
156: }
157: return result;
158: }
159:
160: /**
161: * Converts geocentric coordinates (x, y, z) to geodetic coordinates (longitude, latitude, height), according to the
162: * current ellipsoid parameters. The method used here is derived from "An Improved Algorithm for Geocentric to
163: * Geodetic Coordinate Conversion", by Ralph Toms, Feb 1996 UCRL-JC-123138.
164: *
165: * @param srcPts
166: * the points which must be transformed.
167: */
168: protected void toGeographic(List<Point3d> srcPts) {
169: for (Point3d p : srcPts) {
170: // Note: Variable names follow the notation used in Toms, Feb 1996
171: // final double W2 = x * x + y * y; // square of distance from Z axis
172:
173: final double T0 = p.z * AD_C; // initial estimate of vertical component
174: final double W = ProjectionUtils.length(p.x, p.y);// Math.sqrt( W2 ); // distance from Z axis
175: final double S0 = ProjectionUtils.length(T0, W);// Math.sqrt( T0 * T0 + W*W ); // initial estimate of
176: // horizontal
177: // component
178:
179: final double sin_B0 = T0 / S0; // sin(B0), B0 is estimate of Bowring variable
180: final double cos_B0 = W / S0; // cos(B0)
181: final double sin3_B0 = sin_B0 * sin_B0 * sin_B0; // cube of sin(B0)
182: final double T1 = p.z + semiMinorAxis * ep2 * sin3_B0; // corrected estimate of vertical component
183:
184: // numerator of cos(phi1)
185: final double sum = W - semiMajorAxis * squaredEccentricity
186: * (cos_B0 * cos_B0 * cos_B0);
187:
188: // corrected estimate of horizontal component
189: final double S1 = ProjectionUtils.length(T1, sum);// Math.sqrt( T1 * T1 + sum * sum );
190:
191: // sin(phi), phi is estimated latitude
192: final double sinPhi = T1 / S1;
193: final double cosPhi = sum / S1; // cos(phi)
194:
195: //Lambda in tom.
196: p.x = Math.atan2(p.y, p.x);// longitude;
197: p.y = Math.atan(sinPhi / cosPhi);// latitude;
198: if (hasHeight) {
199: double height = 1;
200: //rn = radius of curvature of the prime vertical, of the ellipsoid at location
201: final double rn = semiMajorAxis
202: / Math.sqrt(1 - squaredEccentricity
203: * (sinPhi * sinPhi));
204:
205: if (cosPhi >= +COS_67P5) {
206: height = W / +cosPhi - rn;
207: } else if (cosPhi <= -COS_67P5) {
208: height = W / -cosPhi - rn;
209: } else {
210: height = p.z / sinPhi + rn
211: * (squaredEccentricity - 1.0);
212: }
213: p.z = height;
214: }
215: }
216: }
217:
218: protected void toGeoCentric(List<Point3d> srcPts) {
219: for (Point3d p : srcPts) {
220: final double lambda = p.x; // Longitude
221: final double phi = p.y; // Latitude
222: final double h = hasHeight ? p.z : 0; // Height above the ellipsoid (metres).
223:
224: final double cosPhi = Math.cos(phi);
225: final double sinPhi = Math.sin(phi);
226: final double rn = semiMajorAxis
227: / Math.sqrt(1 - squaredEccentricity
228: * (sinPhi * sinPhi));
229:
230: p.x = (rn + h) * cosPhi * Math.cos(lambda);
231: p.y = (rn + h) * cosPhi * Math.sin(lambda);
232: p.z = (rn * (1 - squaredEccentricity) + h) * sinPhi;
233: }
234: }
235:
236: @Override
237: public boolean isIdentity() {
238: return false;
239: }
240:
241: /**
242: * @return the semiMajorAxis.
243: */
244: public final double getSemiMajorAxis() {
245: return semiMajorAxis;
246: }
247:
248: /**
249: * @return the semiMinorAxis.
250: */
251: public final double getSemiMinorAxis() {
252: return semiMinorAxis;
253: }
254:
255: @Override
256: public String getName() {
257: return "Geocentric-Transform";
258: }
259:
260: }
|