    001    //$HeadURL: https://svn.wald.intevation.org/svn/deegree/base/branches/2.4_testing/src/org/deegree/crs/projections/azimuthal/StereographicAzimuthal.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.EPS10;
    040    import static org.deegree.crs.projections.ProjectionUtils.EPS11;
    041    import static org.deegree.crs.projections.ProjectionUtils.HALFPI;
    042    import static org.deegree.crs.projections.ProjectionUtils.QUARTERPI;
    043    import static org.deegree.crs.projections.ProjectionUtils.calcPhiFromConformalLatitude;
    044    import static org.deegree.crs.projections.ProjectionUtils.conformalLatitudeInnerPart;
    045    import static org.deegree.crs.projections.ProjectionUtils.length;
    046    import static org.deegree.crs.projections.ProjectionUtils.preCalcedThetaSeries;
    047    import static org.deegree.crs.projections.ProjectionUtils.tanHalfCoLatitude;
    048
    049    import javax.vecmath.Point2d;
    050
    051    import org.deegree.crs.Identifiable;
    052    import org.deegree.crs.components.Unit;
    053    import org.deegree.crs.coordinatesystems.GeographicCRS;
    054    import org.deegree.crs.exceptions.ProjectionException;
    055
    056    /**
    057     * The <code>StereographicAzimuthal</code> class allows for Stereographic Projections of the Poles, equator as well as
    058     * oblique. This projection has following properties (Snyder p. 154):
    059     * <ul>
    060     * <li>Azimuthal</li>
    061     * <li>Conformal</li>
    062     * <li>The central meridian an a particular parallel (if shown) are straight lines.</li>
    063     * <li>All meridians on the polar aspects and the equator on the equatorial aspect are straight lines.</li>
    064     * <li>All other meidians and parallels are shown as arcs of circles</li>
    065     * <li>A perspective projections for the sphere</li>
    066     * <li>Directions from the center of the projection are true (except on ellipsoidal oblique and equatorial aspects)</li>
    067     * <li>Scale increases away from the center of the projection</li>
    068     * <li>Point opposite the center of the projection cannot be plotted</li>
    069     * <li>Used for polar maps and miscellaneous special maps</li>
    070     * </ul>
    071     *
    072     * <p>
    073     * Like Orthographic, the stereographic projection is a true perspective in its isSpherical() form. It is the only known
    074     * true perspective projection of any kind that is also conformal. Its point of projection is on the the surface of the
    075     * sphere at a point jus opposite the oint of tangency of the plane or the center point of the projection. Thus, if the
    076     * north pole is the center of the map, the projection is from the south-pole.
    077     * </p>
    078     * <p>
    079     * It is known to be used by following epsg transformations:
    080     * <ul>
    081     * <li>EPSG:28992</li>
    082     * </ul>
    083     * </p>
    084     *
    085     * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
    086     *
    087     * @author last edited by: $Author: rbezema $
    088     *
    089     * @version $Revision: 18388 $, $Date: 2009-07-08 15:50:37 +0200 (Mi, 08. Jul 2009) $
    090     *
    091     */
    092
    093    public class StereographicAzimuthal extends AzimuthalProjection {
    094
    095        private static final long serialVersionUID = 5399110969291553925L;
    096
    097        /**
    098         * This variable will hold different values, for the ellipsoidal projection:
    099         * <ul>
    100         * <li>for Oblique projection it is the first Part (2*a*k_0*m_1, hence it's name) of A (p.160 21-27)</li>
    101         * <li>for the north-south it may hold the rho of (21-33) or (21-34)
    102         * <li>for the equatorial it holds 2*scale.</li>
    103         * <li>For the
    104         * </ul>
    105         * For the spherical projection:
    106         * <ul>
    107         * <li>For oblique and equatorial: 2*scale, which is compliant with 2*k_0 from Snyder (p.157 21-4)</li>
    108         * <li>For north and south: cos(truelat) / tan( pi/4 - phi/2) ????</li>
    109         * </ul>
    110         */
    111        private double akm1;
    112
    113        private double trueScaleLatitude;
    114
    115        private double[] preCalcedPhiSeries;
    116
    117        // use for the OBLIQUE projection instead of the projectionLatitude
    118        private double conformalLatitude = 0;
    119
    120        // sine of the conformal latitude
    121        private double sinCL = 0;
    122
    123        // cosine of the conformal latitude
    124        private double cosCL = 0;
    125
    126        /**
    127         * @param trueScaleLatitude
    128         *            the latitude (in radians) of a circle around the projection point, which contains the true scale.
    129         * @param geographicCRS
    130         * @param falseNorthing
    131         * @param falseEasting
    132         * @param naturalOrigin
    133         * @param units
    134         * @param scale
    135         * @param id
    136         *            an identifiable instance containing information about this projection
    137         */
    138        public StereographicAzimuthal( double trueScaleLatitude, GeographicCRS geographicCRS, double falseNorthing,
    139                                       double falseEasting, Point2d naturalOrigin, Unit units, double scale, Identifiable id ) {
    140            super( geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale, true,
    141                   false/* not equal area */, id );
    142
    143            preCalcedPhiSeries = preCalcedThetaSeries( getSquaredEccentricity() );
    144
    145            this.trueScaleLatitude = trueScaleLatitude;
    146            if ( !isSpherical() ) {
    147                // double X;
    148                double tmp = 0;
    149                switch ( getMode() ) {
    150                case NORTH_POLE:
    151                case SOUTH_POLE:
    152                    /**
    153                     * The t from 21-33 and 21-34 is still to be calculated.
    154                     */
    155                    // If true scale or known scale factor k_0 is to occur at the pole, the following applies:
    156                    if ( Double.isNaN( trueScaleLatitude )
    157                         || Math.abs( Math.abs( this.trueScaleLatitude ) - HALFPI ) < EPS10 ) {
    158                        // Snyder (p.161 21-33) akm1 = rho .
    159                        akm1 = 2.
    160                               * getScaleFactor()
    161                               / Math.sqrt( Math.pow( 1 + getEccentricity(), 1 + getEccentricity() )
    162                                            * Math.pow( 1 - getEccentricity(), 1 - getEccentricity() ) );
    163
    164                    } else {
    165                        // For true scale along the circle represented by trueScaleLatitude (phi_c in snyder)..
    166                        // akm1 will hold m_c/t_c
    167                        // Calculate the rho from Snyder (p.161 21-34) and place in akm1
    168                        tmp = Math.sin( this.trueScaleLatitude );
    169
    170                        // First part of m of Snyder (p.160 14-15)
    171                        akm1 = Math.cos( this.trueScaleLatitude )
    172                               / tanHalfCoLatitude( this.trueScaleLatitude, tmp, getEccentricity() );
    173                        tmp *= getEccentricity();
    174                        // Second part of m of Snyder (p.160 14-15), t = (e*sin(phi_c))
    175                        akm1 /= Math.sqrt( 1. - tmp * tmp );
    176
    177                    }
    178                    break;
    179                case EQUATOR:
    180                    akm1 = 2. * getScaleFactor();
    181                    break;
    182                case OBLIQUE:
    183                    tmp = getSinphi0();// Math.sin( getProjectionLatitude() );
    184                    // Calculate the ConformalLatitude for this ellipsoid Snyder (p.160 3-1) the X
    185                    conformalLatitude = ( 2. * Math.atan( conformalLatitudeInnerPart( getProjectionLatitude(), tmp,
    186                                                                                      getEccentricity() ) ) )
    187                                        - HALFPI;
    188
    189                    sinCL = Math.sin( conformalLatitude );
    190                    cosCL = Math.cos( conformalLatitude );
    191
    192                    tmp *= getEccentricity();
    193                    // the first part (2a*k_0*m) of the largeA in Snyder (p.160 21-27) is stored in akm1 it still must be
    194                    // devided by the cos X ... etc. to get largeA
    195                    akm1 = 2. * getScaleFactor() * getCosphi0() / Math.sqrt( 1. - tmp * tmp );
    196
    197                    // Setting the sinPhi0 to the conformal Latitude X.
    198                    // setProjectionLatitude( X );
    199                    break;
    200                }
    201            } else {
    202                akm1 = 2. * getScaleFactor();
    203            }
    204        }
    205
    206        /**
    207         * Sets the id to "Snyder-StereoGraphic"
    208         *
    209         * @param trueScaleLatitude
    210         *            the latitude (in radians) of a circle around the projection point, which contains the true scale.
    211         * @param geographicCRS
    212         * @param falseNorthing
    213         * @param falseEasting
    214         * @param naturalOrigin
    215         * @param units
    216         * @param scale
    217         */
    218        public StereographicAzimuthal( double trueScaleLatitude, GeographicCRS geographicCRS, double falseNorthing,
    219                                       double falseEasting, Point2d naturalOrigin, Unit units, double scale ) {
    220            this( trueScaleLatitude, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale,
    221                  new Identifiable( "Snyder-StereoGraphic" ) );
    222        }
    223
    224        /**
    225         * Create a {@link StereographicAzimuthal} which has a true scale latitude at MapUtils.HALFPI.
    226         *
    227         * @param geographicCRS
    228         * @param falseNorthing
    229         * @param falseEasting
    230         * @param naturalOrigin
    231         * @param units
    232         * @param scale
    233         * @param id
    234         *            an identifiable instance containing information about this projection
    235         */
    236        public StereographicAzimuthal( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
    237                                       Point2d naturalOrigin, Unit units, double scale, Identifiable id ) {
    238            this( HALFPI, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale, id );
    239        }
    240
    241        /**
    242         * Create a {@link StereographicAzimuthal} which has a true scale latitude at MapUtils.HALFPI. Sets the id to
    243         * "Snyder-StereoGraphic"
    244         *
    245         * @param geographicCRS
    246         * @param falseNorthing
    247         * @param falseEasting
    248         * @param naturalOrigin
    249         * @param units
    250         * @param scale
    251         */
    252        public StereographicAzimuthal( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
    253                                       Point2d naturalOrigin, Unit units, double scale ) {
    254            this( HALFPI, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale );
    255        }
    256
    257        /**
    258         * Create a {@link StereographicAzimuthal} which has a scale of 1 and a true scale latitude,
    259         *
    260         * @param trueScaleLatitude
    261         *            the latitude (in radians) of a circle around the projection point, which contains the true scale.
    262         * @param geographicCRS
    263         * @param falseNorthing
    264         * @param falseEasting
    265         * @param naturalOrigin
    266         * @param units
    267         * @param id
    268         *            an identifiable instance containing information about this projection
    269         */
    270        public StereographicAzimuthal( double trueScaleLatitude, GeographicCRS geographicCRS, double falseNorthing,
    271                                       double falseEasting, Point2d naturalOrigin, Unit units, Identifiable id ) {
    272            this( trueScaleLatitude, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, 1, id );
    273        }
    274
    275        /**
    276         * Create a {@link StereographicAzimuthal} which has a scale of 1 and a true scale latitude. Sets the id to
    277         * "Snyder-StereoGraphic".
    278         *
    279         * @param trueScaleLatitude
    280         *            the latitude (in radians) of a circle around the projection point, which contains the true scale.
    281         * @param geographicCRS
    282         * @param falseNorthing
    283         * @param falseEasting
    284         * @param naturalOrigin
    285         * @param units
    286         */
    287        public StereographicAzimuthal( double trueScaleLatitude, GeographicCRS geographicCRS, double falseNorthing,
    288                                       double falseEasting, Point2d naturalOrigin, Unit units ) {
    289            this( trueScaleLatitude, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, 1 );
    290        }
    291
    292        /**
    293         * Create a {@link StereographicAzimuthal} which is conformal, has a scale of 1 and a truescale latitude at pi*0.5.
    294         *
    295         * @param geographicCRS
    296         * @param falseNorthing
    297         * @param falseEasting
    298         * @param naturalOrigin
    299         * @param units
    300         * @param id
    301         *            an identifiable instance containing information about this projection
    302         */
    303        public StereographicAzimuthal( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
    304                                       Point2d naturalOrigin, Unit units, Identifiable id ) {
    305            this( HALFPI, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, 1, id );
    306        }
    307
    308        /**
    309         * Create a {@link StereographicAzimuthal} which is conformal, has a scale of 1 and a truescale latitude at pi*0.5.
    310         * Sets the id to "Snyder-StereoGraphic".
    311         *
    312         * @param geographicCRS
    313         * @param falseNorthing
    314         * @param falseEasting
    315         * @param naturalOrigin
    316         * @param units
    317         */
    318        public StereographicAzimuthal( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
    319                                       Point2d naturalOrigin, Unit units ) {
    320            this( HALFPI, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, 1 );
    321        }
    322
    323        @Override
    324        public Point2d doInverseProjection( double x, double y ) {
    325            Point2d result = new Point2d( 0, 0 );
    326            x -= getFalseEasting();
    327            y -= getFalseNorthing();
    328            double rho = length( x, y );
    329            if ( isSpherical() ) {
    330                // akm1 = 2*scale.
    331                double c = 2. * Math.atan( rho / akm1 );
    332                double sinc = Math.sin( c );
    333                double cosc = Math.cos( c );
    334                switch ( getMode() ) {
    335                case OBLIQUE:
    336                    // First find theta from Snyder (p.158 20-14).
    337                    if ( Math.abs( rho ) <= EPS10 ) {
    338                        result.y = getProjectionLatitude();
    339                    } else {
    340                        result.y = Math.asin( cosc * getSinphi0() + ( ( y * sinc * getCosphi0() ) / rho ) );
    341                    }
    342                    // And now lambda but instead of using the atan, use atan2. Is this correct????
    343                    c = cosc - getSinphi0() * Math.sin( result.y );
    344                    if ( Math.abs( c ) > EPS10 || x != 0. ) {
    345                        result.x = Math.atan2( x * sinc * getCosphi0(), c * rho );
    346                    }
    347                    break;
    348                case EQUATOR:
    349                    if ( Math.abs( rho ) > EPS10 ) {
    350                        result.y = Math.asin( y * sinc / rho );
    351                    }
    352                    if ( cosc != 0. || x != 0. ) {
    353                        result.x = Math.atan2( x * sinc, cosc * rho );
    354                    }
    355                    break;
    356                case NORTH_POLE:
    357                    y = -y;
    358                case SOUTH_POLE:
    359                    if ( Math.abs( rho ) <= EPS10 ) {
    360                        result.y = getProjectionLatitude();
    361                    } else {
    362                        result.y = Math.asin( getMode() == SOUTH_POLE ? -cosc : cosc );
    363                    }
    364                    result.x = Math.atan2( x, y );
    365                    break;
    366                }
    367            } else {
    368                double cosCe = 0;
    369                double sinCe = 0;
    370                double X = 0;
    371                // for oblique and equatorial projections, akm1 = first part of Snyder (p.160 21-27). which is 2*scale*m_1
    372                // double ce = 2. * Math.atan2( rho * getCosphi0(), akm1 );
    373                // if ( getMode() == OBLIQUE || getMode() == EQUATOR ) {
    374                // cosCe = Math.cos( ce );
    375                // sinCe = Math.sin( ce );
    376                // }
    377                double ce = 0;
    378                switch ( getMode() ) {
    379                case OBLIQUE:
    380                    ce = 2. * Math.atan2( rho * cosCL, akm1 );
    381                    cosCe = Math.cos( ce );
    382                    sinCe = Math.sin( ce );
    383                    X = Math.asin( cosCe * sinCL + ( ( y * sinCe * cosCL ) / rho ) );
    384                    result.x = Math.atan2( x * sinCe, rho * cosCL * cosCe - y * sinCL * sinCe );
    385                    break;
    386                case EQUATOR:
    387                    ce = 2. * Math.atan2( rho * getCosphi0(), akm1 );
    388                    cosCe = Math.cos( ce );
    389                    sinCe = Math.sin( ce );
    390                    X = Math.asin( cosCe * getSinphi0() + ( ( y * sinCe * getCosphi0() ) / rho ) );
    391                    result.x = Math.atan2( x * sinCe, rho * getCosphi0() * cosCe - y * getSinphi0() * sinCe );
    392                    break;
    393                case NORTH_POLE:
    394                    y = -y;
    395                case SOUTH_POLE:
    396                    /**
    397                     * akm1 can hold rho from 21-34 or 21-33 (if no truescale latitude was given) without the parameter t.
    398                     * Either way it must only be multiplied with t to fulfill the formula. Independent of the true scale
    399                     * latitude, the value of akm1 must therefore only be inverted to calculate Snyder (p.162 21-39 or
    400                     * 21-40).
    401                     */
    402                    double t = rho * 1. / this.akm1;
    403                    X = HALFPI - 2 * Math.atan( t );
    404                    if ( getMode() == SOUTH_POLE ) {
    405                        X *= -1;
    406                    }
    407                    result.x = Math.atan2( x, y );
    408                    break;
    409                }
    410                result.y = calcPhiFromConformalLatitude( X, preCalcedPhiSeries );
    411                result.x += getProjectionLongitude();
    412            }
    413            return result;
    414        }
    415
    416        @Override
    417        public Point2d doProjection( double lambda, double phi )
    418                                throws ProjectionException {
    419            Point2d result = new Point2d( 0, 0 );
    420            lambda -= getProjectionLongitude();
    421            double cosLamda = Math.cos( lambda );
    422            double sinLamda = Math.sin( lambda );
    423            double sinPhi = Math.sin( phi );
    424
    425            if ( isSpherical() ) {
    426                double cosphi = Math.cos( phi );
    427
    428                switch ( getMode() ) {
    429                case OBLIQUE:
    430                    // Calc k from Snyder (p.158 21-14) and store in xy.y
    431                    result.y = 1. + getSinphi0() * sinPhi + getCosphi0() * cosphi * cosLamda;
    432                    if ( result.y <= EPS11 ) {
    433                        throw new ProjectionException( "Division by zero" );
    434                    }
    435                    // akm1 was set to 2*scale.
    436                    result.y = akm1 / result.y;
    437
    438                    // Calc x (p.157 21-2) and y (p.157 21-3)
    439                    result.x = result.y * cosphi * sinLamda;
    440                    result.y *= getCosphi0() * sinPhi - getSinphi0() * cosphi * cosLamda;
    441                    break;
    442                case EQUATOR:
    443                    // Calc k from Snyder (p.158 21-14) and store in xy.y
    444                    result.y = 1. + cosphi * cosLamda;
    445                    if ( result.y <= EPS11 ) {
    446                        throw new ProjectionException( "Division by zero" );
    447                    }
    448                    // akm1 was set to 2*scale.
    449                    result.y = akm1 / result.y;
    450
    451                    // Calc x (p.157 21-2) and y (p.158 21-13)
    452                    result.x = result.y * cosphi * sinLamda;
    453                    result.y *= sinPhi;
    454                    break;
    455                case NORTH_POLE:
    456                    cosLamda = -cosLamda;
    457                    phi = -phi;
    458                case SOUTH_POLE:
    459                    if ( Math.abs( phi - HALFPI ) < EPS10 ) {
    460                        throw new ProjectionException( "The requested phi: " + phi + " lies on the singularity ("
    461                                                       + ( ( getMode() == SOUTH_POLE ) ? "South-Pole" : "North-Pole" )
    462                                                       + ") of this projection's mappable area." );
    463                    }
    464                    // akm1 was set to 2*scale.
    465                    result.y = akm1 * Math.tan( QUARTERPI + .5 * phi );
    466                    result.x = sinLamda * ( result.y );
    467                    result.y *= cosLamda;
    468                    break;
    469                }
    470            } else {
    471                double sinX = 0;
    472                double cosX = 0;
    473                double largeA;
    474
    475                if ( getMode() == OBLIQUE || getMode() == EQUATOR ) {
    476                    // Calculate the conformal latitue Snyder (p.160 3-1).
    477                    double X = 2. * Math.atan( conformalLatitudeInnerPart( phi, sinPhi, getEccentricity() ) ) - HALFPI;
    478                    sinX = Math.sin( X );
    479                    cosX = Math.cos( X );
    480                }
    481                switch ( getMode() ) {
    482                case OBLIQUE:
    483                    // akm1 was set to the first part of Snyder (p.160 21-27)
    484                    // sinPhi0 and cosPhi0 where set to the conformal latitude of the natural origin.
    485                    largeA = akm1 / ( cosCL * ( 1. + sinCL * sinX + cosCL * cosX * cosLamda ) );
    486                    result.y = largeA * ( cosCL * sinX - sinCL * cosX * cosLamda );
    487                    result.x = largeA * cosX;
    488                    break;
    489                case EQUATOR:
    490                    // akm1 was set to 2. * scale;
    491                    // Calculate the A of Snyder (p.161 21-29).
    492                    largeA = 2. * akm1 / ( 1. + cosX * cosLamda );
    493                    result.y = largeA * sinX;
    494                    result.x = largeA * cosX;
    495                    break;
    496                case SOUTH_POLE:
    497                    phi = -phi;
    498                    cosLamda = -cosLamda;
    499                    sinPhi = -sinPhi;
    500                case NORTH_POLE:
    501                    // akm1 holds the rho's from 21-33 or 21-34, but without the t (==halfColatitude of the conformal
    502                    // latitude).
    503                    result.x = akm1 * tanHalfCoLatitude( phi, sinPhi, getEccentricity() );
    504                    result.y = -result.x * cosLamda;
    505                    break;
    506                }
    507                // Sinus(Lambda) was still to be multiplied on the x.
    508                result.x = result.x * sinLamda;
    509            }
    510            result.x += getFalseEasting();
    511            result.y += getFalseNorthing();
    512            return result;
    513        }
    514
    515        @Override
    516        public String getImplementationName() {
    517            return "stereographicAzimuthal";
    518        }
    519
    520        /**
    521         * @return the trueScaleLatitude.
    522         */
    523        public final double getTrueScaleLatitude() {
    524            return trueScaleLatitude;
    525        }
    526    }