    001    //$HeadURL: https://svn.wald.intevation.org/svn/deegree/base/branches/2.4_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    }