    001    //$HeadURL: https://svn.wald.intevation.org/svn/deegree/base/branches/2.4_testing/src/org/deegree/crs/projections/azimuthal/StereographicAlternative.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.projections.azimuthal;
    038
    039    import static org.deegree.crs.projections.ProjectionUtils.EPS11;
    040    import static org.deegree.crs.projections.ProjectionUtils.HALFPI;
    041    import static org.deegree.crs.projections.ProjectionUtils.QUARTERPI;
    042
    043    import javax.vecmath.Point2d;
    044
    045    import org.deegree.crs.Identifiable;
    046    import org.deegree.crs.components.Unit;
    047    import org.deegree.crs.coordinatesystems.GeographicCRS;
    048    import org.deegree.crs.exceptions.ProjectionException;
    049    import org.deegree.framework.log.ILogger;
    050    import org.deegree.framework.log.LoggerFactory;
    051
    052    /**
    053     * <code>StereographicAlternative</code> projection may be imagined to be a projection of the earth's surface onto a
    054     * plane in contact with the earth at a single tangent point from the opposite end of the diameter through that tangent
    055     * point.
    056     * <p>
    057     * An alternative approach is given by Snyder {@link StereographicAzimuthal}, where, instead of defining a single
    058     * conformal sphere at the origin point, the conformal latitude at each point on the ellipsoid is computed. The
    059     * conformal longitude is then always equivalent to the geodetic longitude. This approach is a valid alternative to the
    060     * above, but gives slightly different results away from the origin point. It is therefore considered by EPSG to be a
    061     * different projection method. Hence this implementation.
    062     * </p>
    063     * <p>
    064     * This projection is best known in its polar form and is frequently used for mapping polar areas where it complements
    065     * the Universal Transverse Mercator used for lower latitudes. Its spherical form has also been widely used by the US
    066     * Geological Survey for planetary mapping and the mapping at small scale of continental hydrocarbon provinces. In its
    067     * transverse or oblique ellipsoidal forms it is useful for mapping limited areas centered on the point where the plane
    068     * of the projection is regarded as tangential to the ellipsoid., e.g. the Netherlands. The tangent point is the origin
    069     * of the projected coordinate system and the meridian through it is regarded as the central meridian. In order to
    070     * reduce the scale error at the extremities of the projection area it is usual to introduce a scale factor of less than
    071     * unity at the origin such that a unit scale factor applies on a near circle centered at the origin and some distance
    072     * from it.
    073     * </p>
    074     *
    075     * <p>
    076     * The coordinate transformation from geographical to projected coordinates is executed via the distance and azimuth of
    077     * the point from the center point or origin. For a sphere the formulas are relatively simple. For the ellipsoid the
    078     * same formulas are used but with auxiliary latitudes, known as conformal latitudes, substituted for the geodetic
    079     * latitudes of the spherical formulas for the origin and the point .
    080     * </p>
    081     * <quote>from <a
    082     * href="http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34j.html">http://www.posc.org/</a></quote>
    083     *
    084     * <p>
    085     * Determinations of oblique projections on an ellipsoid can be difficult to solve and result in long, complex
    086     * computations. Because conformal transformations can be made multiple time without loss of the conformal property a
    087     * method of determining oblique projections involves conformal transformation of the elliptical coordinates to
    088     * coordinates on a conformal sphere. The transformed coordinates can now be translated/rotated on the sphere and then
    089     * converted to planar coordinates with a conformal spherical projection. is Ce/2 1 − e sin φ 2 arctan K tanC (π/4 +
    090     * φ/2) − π/2 (3.6) χ = 1 + e sin φ λc = Cλ (3.7) √ 1 − e2 Rc = (3.8) 1 − e2 sin2 φ0
    091     * </p>
    092     * From the <a href="http://members.verizon.net/~gerald.evenden/proj4/manual.pdf">libproj4-manual</a> by Gerald I.
    093     * Evenden
    094     *
    095     *
    096     * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
    097     *
    098     * @author last edited by: $Author: rbezema $
    099     *
    100     * @version $Revision: 19653 $, $Date: 2009-09-15 14:56:30 +0200 (Di, 15. Sep 2009) $
    101     *
    102     */
    103    public class StereographicAlternative extends AzimuthalProjection {
    104
    105        private final static ILogger LOG = LoggerFactory.getLogger( StereographicAlternative.class );
    106
    107        private static final long serialVersionUID = -444957611911311171L;
    108
    109        private double sinc0 = 0;
    110
    111        private double cosc0 = 0;
    112
    113        /**
    114         * Two times the radius of the conformal sphere, which is denoted by Gerald I. Evenden as R_c
    115         */
    116        private final double R2;
    117
    118        /**
    119         * THe latitude of on the conformal sphere which is denoted by Gerald I. Evenden as Chi_0
    120         */
    121        private final double latitudeOnCS;
    122
    123        /**
    124         * The exponent of the calculation of conformal latitude Chi in Gerald I. Evenden (3.6)
    125         */
    126        private final double clExponent;
    127
    128        /**
    129         * The value K in Gerald I. Evenden (3.11)
    130         */
    131        private final double K;
    132
    133        /**
    134         * The central geographic latitude of the projection on the conformal sphere, denoted by Gerald I. Evenden (3.9)
    135         * with C
    136         */
    137        private final double centralGeographicLatitude;
    138
    139        // private double chi = 0;
    140
    141        // private double radiusOfCS = 0;
    142
    143        /**
    144         * @param geographicCRS
    145         * @param falseNorthing
    146         * @param falseEasting
    147         * @param naturalOrigin
    148         * @param units
    149         * @param scale
    150         * @param id
    151         *            an identifiable instance containing information about this projection
    152         */
    153        public StereographicAlternative( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
    154                                         Point2d naturalOrigin, Unit units, double scale, Identifiable id ) {
    155            super( geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale, true/* conformal */, false, id );
    156
    157            // if (!(P->en = pj_gauss_ini(P->e, P->phi0, &(P->phic0), &R))) E_ERROR_0;
    158            double es = getSquaredEccentricity();
    159            double sinPhi0 = getSinphi0();
    160            double cosPhiSquare = getCosphi0() * getCosphi0();// Math.cos( getProjectionLatitude() );
    161
    162            // R_c
    163            double radiusOfCS = Math.sqrt( 1. - es ) / ( 1. - es * sinPhi0 * sinPhi0 );
    164            centralGeographicLatitude = Math.sqrt( 1. + es * cosPhiSquare * cosPhiSquare / ( 1. - es ) );
    165            // xsi_0
    166            latitudeOnCS = Math.asin( sinPhi0 / centralGeographicLatitude );
    167            clExponent = 0.5 * centralGeographicLatitude * getEccentricity();
    168            K = Math.tan( .5 * latitudeOnCS + QUARTERPI )
    169                / ( Math.pow( Math.tan( .5 * getProjectionLatitude() + QUARTERPI ), centralGeographicLatitude ) * srat(
    170                                                                                                                        getEccentricity()
    171                                                                                                                                                * sinPhi0,
    172                                                                                                                        clExponent ) );
    173            sinc0 = Math.sin( latitudeOnCS );
    174            cosc0 = Math.cos( latitudeOnCS );
    175            R2 = 2 * radiusOfCS;
    176        }
    177
    178        /**
    179         * Sets the id of this projection to epsg::9809 (Oblique Stereographic)
    180         *
    181         * @param geographicCRS
    182         * @param falseNorthing
    183         * @param falseEasting
    184         * @param naturalOrigin
    185         * @param units
    186         * @param scale
    187         */
    188        public StereographicAlternative( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
    189                                         Point2d naturalOrigin, Unit units, double scale ) {
    190            this( geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale, new Identifiable( "EPSG::9809" ) );
    191        }
    192
    193        /*
    194         * (non-Javadoc)
    195         *
    196         * @see org.deegree.crs.projections.Projection#doInverseProjection(double, double)
    197         */
    198        @Override
    199        public Point2d doInverseProjection( double x, double y )
    200                                throws ProjectionException {
    201            if ( LOG.isDebug() ) {
    202                LOG.logDebug( "(sterea incoming inv proj) x:" + x + ", y:" + y );
    203            }
    204            Point2d result = new Point2d();
    205            x -= getFalseEasting();
    206            y -= getFalseNorthing();
    207
    208            x /= getScaleFactor();
    209            y /= getScaleFactor();
    210            double rho = Math.hypot( x, y );
    211
    212            if ( rho > EPS11 ) {
    213                double c = 2 * Math.atan2( rho, R2 );
    214                double sinc = Math.sin( c );
    215                double cosc = Math.cos( c );
    216                result.y = Math.asin( cosc * sinc0 + y * sinc * cosc0 / rho );
    217                result.x = Math.atan2( x * sinc, rho * cosc0 * cosc - y * sinc0 * sinc );
    218            } else {
    219                result.y = latitudeOnCS;
    220                result.x = 0;
    221            }
    222            result = pj_inv_gauss( result );
    223            result.x += getProjectionLongitude();
    224            if ( LOG.isDebug() ) {
    225                LOG.logDebug( "(outgoing) lam:" + Math.toDegrees( result.x ) + ", phi:" + Math.toDegrees( result.y ) );
    226            }
    227            return result;
    228        }
    229
    230        /*
    231         * (non-Javadoc)
    232         *
    233         * @see org.deegree.crs.projections.Projection#doProjection(double, double)
    234         */
    235        @Override
    236        public Point2d doProjection( double lambda, double phi )
    237                                throws ProjectionException {
    238            if ( LOG.isDebug() ) {
    239                LOG.logDebug( "(sterea incoming proj) lam:" + Math.toDegrees( lambda ) + ", phi:" + Math.toDegrees( phi ) );
    240            }
    241            Point2d result = new Point2d();
    242            lambda -= getProjectionLongitude();
    243            Point2d lp = pj_gauss( lambda, phi );
    244            double sinc = Math.sin( lp.y );
    245            double cosc = Math.cos( lp.y );
    246            double cosl = Math.cos( lp.x );
    247
    248            double k = getScaleFactor() * ( R2 / ( 1. + sinc0 * sinc + cosc0 * cosc * cosl ) );
    249            result.x = k * cosc * Math.sin( lp.x );
    250            result.y = k * ( cosc0 * sinc - sinc0 * cosc * cosl );
    251
    252            result.x += getFalseEasting();
    253            result.y += getFalseNorthing();
    254            if ( LOG.isDebug() ) {
    255                LOG.logDebug( "(outgoing) x:" + result.x + ", y:" + result.y );
    256            }
    257            return result;
    258        }
    259
    260        /*
    261         * (non-Javadoc)
    262         *
    263         * @see org.deegree.crs.projections.Projection#getDeegreeSpecificName()
    264         */
    265        @Override
    266        public String getImplementationName() {
    267            return "stereographicAlternative";
    268        }
    269
    270        private static double srat( double esinp, double exp ) {
    271            return Math.pow( ( 1. - esinp ) / ( 1. + esinp ), exp );
    272        }
    273
    274        /**
    275         * To determine the inverse solution, geographic coordinates from Gaussian sphere coordinates, execute with the
    276         * initial value of φi−1 = χ and φi−1 iteratively replaced by φ until |φ − φi−1 | is less than an acceptable error
    277         * value. (taken from the proj-lib manual.
    278         */
    279        private Point2d pj_inv_gauss( Point2d slp )
    280                                throws ProjectionException {
    281            Point2d elp = new Point2d();
    282
    283            elp.x = slp.x / centralGeographicLatitude;
    284            double num = Math.pow( Math.tan( .5 * slp.y + QUARTERPI ) / K, 1. / centralGeographicLatitude );
    285            int MAX_ITER = 20;
    286            int i = MAX_ITER;
    287            for ( ; i > 0; --i ) {
    288                elp.y = 2. * Math.atan( num * srat( getEccentricity() * Math.sin( slp.y ), -.5 * getEccentricity() ) )
    289                        - HALFPI;
    290                if ( Math.abs( ( elp.y - slp.y ) ) < EPS11 ) {
    291                    break;
    292                }
    293                slp.y = elp.y;
    294            }
    295            /* convergence failed */
    296            if ( i == 0 ) {
    297                throw new ProjectionException( "No convertgence while calculation the inverse gaus approximation" );
    298            }
    299            return elp;
    300        }
    301
    302        /**
    303         * The conformal transformation of ellipsoid coordinates (φ, λ) to conformal sphere coordinates (χ, λ_c ), where λ
    304         * is relative to the longitude of projection origin, R_c is radius of the conformal sphere. χ_0 is the latitude on
    305         * the conformal sphere at the central geographic latitude of the projection.
    306         */
    307        private Point2d pj_gauss( double lambda, double phi ) {
    308            Point2d slp = new Point2d();
    309            slp.y = 2.
    310                    * Math.atan( K * Math.pow( Math.tan( .5 * phi + QUARTERPI ), centralGeographicLatitude )
    311                                 * srat( getEccentricity() * Math.sin( phi ), clExponent ) ) - HALFPI;
    312            slp.x = centralGeographicLatitude * ( lambda );
    313            return slp;
    314        }
    315
    316    }