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