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
052 /**
053 * The <code>GeocentricTransform</code> class is used to create a transformation between a geocentric CRS (having
054 * lat-lon coordinates) and it's geodetic CRS (having x-y-z) coordinates and vice versa.
055 *
056 * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
057 *
058 * @author last edited by: $Author:$
059 *
060 * @version $Revision:$, $Date:$
061 *
062 */
063
064 public class GeocentricTransform extends CRSTransformation {
065
066 /**
067 * Cosine of 67.5 degrees.
068 */
069 private static final double COS_67P5 = 0.38268343236508977;
070
071 /**
072 * Toms region 1 constant, which will allow for a difference between ellipsoid of 2000 Kilometers. <quote>Under this
073 * policy the maximum error is less than 42 centimeters for altitudes less then 10.000.000 Kilometers</quote>
074 */
075 private static final double AD_C = 1.0026;
076
077 /**
078 * Semi-major axis of ellipsoid in meters.
079 */
080 private final double semiMajorAxis;
081
082 /**
083 * Semi-minor axis of ellipsoid in meters.
084 */
085 private final double semiMinorAxis;
086
087 /**
088 * Square of semi-major axis {@link #semiMajorAxis}.
089 */
090 private final double squaredSemiMajorAxis;
091
092 /**
093 * Square of semi-minor axis {@link #semiMinorAxis}.
094 */
095 private final double squaredSemiMinorAxis;
096
097 /**
098 * Eccentricity squared.
099 */
100 private final double squaredEccentricity;
101
102 /**
103 * 2nd eccentricity squared.
104 */
105 private final double ep2;
106
107 /**
108 * true if the given points will use heights (e.g. have a z/height-value).
109 */
110 private boolean hasHeight;
111
112 /**
113 * @param source
114 * the geographic crs.
115 * @param target
116 * the geocentric crs.
117 * @param identifier
118 * @param name
119 * @param version
120 * @param description
121 * @param areaOfUse
122 */
123 public GeocentricTransform( GeographicCRS source, GeocentricCRS target, String identifier, String version,
124 String description, String areaOfUse ) {
125 super( source, target, identifier, "Geocentric-Transform", version, description, areaOfUse );
126 this.hasHeight = source.getDimension() == 3;
127 Ellipsoid ellipsoid = source.getGeodeticDatum().getEllipsoid();
128 semiMajorAxis = Unit.METRE.convert( ellipsoid.getSemiMajorAxis(), ellipsoid.getUnits() );
129 semiMinorAxis = Unit.METRE.convert( ellipsoid.getSemiMinorAxis(), ellipsoid.getUnits() );
130 squaredSemiMajorAxis = semiMajorAxis * semiMajorAxis;
131 squaredSemiMinorAxis = semiMinorAxis * semiMinorAxis;
132 squaredEccentricity = ellipsoid.getSquaredEccentricity();
133 // e2 = ( a2 - b2 ) / a2;
134 ep2 = ( squaredSemiMajorAxis - squaredSemiMinorAxis ) / squaredSemiMinorAxis;
135 }
136
137 /**
138 * @param source
139 * the geographic crs.
140 * @param target
141 * the geocentric crs.
142 * @param identifier
143 */
144 public GeocentricTransform( GeographicCRS source, GeocentricCRS target, String identifier ) {
145 this( source, target, identifier, null, null, null );
146 }
147
148 /**
149 * This constructor creates an identifier with following form: <code>
150 * FROM_id1_TO_id2.</code> The name will be
151 * generated the same way (with getName() instead).
152 *
153 * @param source
154 * the geographic crs.
155 * @param target
156 * the geocentric crs.
157 */
158 public GeocentricTransform( GeographicCRS source, GeocentricCRS target ) {
159 this( source, target, new StringBuilder( "FROM_" ).append( source.getIdentifier() )
160 .append( "_TO_" )
161 .append( target.getIdentifier() )
162 .toString() );
163 }
164
165 @Override
166 public List<Point3d> doTransform( List<Point3d> srcPts ) {
167 List<Point3d> result = new ArrayList<Point3d>( srcPts );
168 if ( isInverseTransform() ) {
169 // System.out.println( "An inverse geocentric transform with incoming points: " + srcPts );
170 toGeographic( result );
171 } else {
172 // System.out.println( "A geocentric transform with incoming points: x: " + Math.toDegrees( srcPts.get(0).x)
173 // + ", y: " + Math.toDegrees( srcPts.get(0).y ) );
174 toGeoCentric( result );
175 }
176 return result;
177 }
178
179 /**
180 * Converts geocentric coordinates (x, y, z) to geodetic coordinates (longitude, latitude, height), according to the
181 * current ellipsoid parameters. The method used here is derived from "An Improved Algorithm for Geocentric to
182 * Geodetic Coordinate Conversion", by Ralph Toms, Feb 1996 UCRL-JC-123138.
183 *
184 * @param srcPts
185 * the points which must be transformed.
186 */
187 protected void toGeographic( List<Point3d> srcPts ) {
188 for ( Point3d p : srcPts ) {
189 // Note: Variable names follow the notation used in Toms, Feb 1996
190 // final double W2 = x * x + y * y; // square of distance from Z axis
191
192 final double T0 = p.z * AD_C; // initial estimate of vertical component
193 final double W = ProjectionUtils.length( p.x, p.y );// Math.sqrt( W2 ); // distance from Z axis
194 final double S0 = ProjectionUtils.length( T0, W );// Math.sqrt( T0 * T0 + W*W ); // initial estimate of
195 // horizontal
196 // component
197
198 final double sin_B0 = T0 / S0; // sin(B0), B0 is estimate of Bowring variable
199 final double cos_B0 = W / S0; // cos(B0)
200 final double sin3_B0 = sin_B0 * sin_B0 * sin_B0; // cube of sin(B0)
201 final double T1 = p.z + semiMinorAxis * ep2 * sin3_B0; // corrected estimate of vertical component
202
203 // numerator of cos(phi1)
204 final double sum = W - semiMajorAxis * squaredEccentricity * ( cos_B0 * cos_B0 * cos_B0 );
205
206 // corrected estimate of horizontal component
207 final double S1 = ProjectionUtils.length( T1, sum );// Math.sqrt( T1 * T1 + sum * sum );
208
209 // sin(phi), phi is estimated latitude
210 final double sinPhi = T1 / S1;
211 final double cosPhi = sum / S1; // cos(phi)
212
213 //Lambda in tom.
214 p.x = Math.atan2( p.y, p.x );// longitude;
215 p.y = Math.atan( sinPhi / cosPhi );// latitude;
216 if ( hasHeight ) {
217 double height = 1;
218 //rn = radius of curvature of the prime vertical, of the ellipsoid at location
219 final double rn = semiMajorAxis / Math.sqrt( 1 - squaredEccentricity * ( sinPhi * sinPhi ) );
220
221 if ( cosPhi >= +COS_67P5 ) {
222 height = W / +cosPhi - rn;
223 } else if ( cosPhi <= -COS_67P5 ) {
224 height = W / -cosPhi - rn;
225 } else {
226 height = p.z / sinPhi + rn * ( squaredEccentricity - 1.0 );
227 }
228 p.z = height;
229 }
230 }
231 }
232
233 protected void toGeoCentric( List<Point3d> srcPts ) {
234 for ( Point3d p : srcPts ) {
235 final double lambda = p.x; // Longitude
236 final double phi = p.y; // Latitude
237 final double h = hasHeight ? p.z : 0; // Height above the ellipsoid (metres).
238
239 final double cosPhi = Math.cos( phi );
240 final double sinPhi = Math.sin( phi );
241 final double rn = semiMajorAxis / Math.sqrt( 1 - squaredEccentricity * ( sinPhi * sinPhi ) );
242
243 p.x = ( rn + h ) * cosPhi * Math.cos( lambda );
244 p.y = ( rn + h ) * cosPhi * Math.sin( lambda );
245 p.z = ( rn * ( 1 - squaredEccentricity ) + h ) * sinPhi;
246 }
247 }
248
249 @Override
250 public boolean isIdentity() {
251 return false;
252 }
253
254 /**
255 * @return the semiMajorAxis.
256 */
257 public final double getSemiMajorAxis() {
258 return semiMajorAxis;
259 }
260
261 /**
262 * @return the semiMinorAxis.
263 */
264 public final double getSemiMinorAxis() {
265 return semiMinorAxis;
266 }
267
268 }