001    //$HeadURL: svn+ssh://jwilden@svn.wald.intevation.org/deegree/base/branches/2.5_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    }