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