001    //$HeadURL: svn+ssh://jwilden@svn.wald.intevation.org/deegree/base/branches/2.5_testing/src/org/deegree/crs/projections/azimuthal/LambertAzimuthalEqualArea.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.calcPhiFromAuthalicLatitude;
044    import static org.deegree.crs.projections.ProjectionUtils.calcQForAuthalicLatitude;
045    import static org.deegree.crs.projections.ProjectionUtils.getAuthalicLatitudeSeriesValues;
046    import static org.deegree.crs.projections.ProjectionUtils.length;
047    
048    import javax.vecmath.Point2d;
049    
050    import org.deegree.crs.Identifiable;
051    import org.deegree.crs.components.Unit;
052    import org.deegree.crs.coordinatesystems.GeographicCRS;
053    import org.deegree.crs.exceptions.ProjectionException;
054    
055    /**
056     * The <code>LambertAzimuthalEqualArea</code> projection has following properties (From J.S. Snyder, Map Projections a
057     * Working Manual p. 182):
058     * <ul>
059     * <li>Azimuthal</li>
060     * <li>Equal-Area</li>
061     * <li>All meridians in the polar aspect, the central meridian in other aspects, and the Equator in the equatorial
062     * aspect are straight lines</li>
063     * <li>The outer meridian of a hemisphere in the equatorial aspect (for the sphere) and the parallels in the polar
064     * aspect (sphere or ellipsoid) are circles.</li>
065     * <li>All other meridians and the parallels are complex curves</li>
066     * <li>Not a perspective projection</li>
067     * <li>Scale decreases radially as the distance increases from the center, the only point without distortion</li>
068     * <li>Directions from the center are true for the sphere and the polar ellipsoidal forms.</li>
069     * <li>Point opposite the center is shown as a circle surrounding the map (for the sphere).</li>
070     * <li>Used for maps of continents and hemispheres</li>
071     * <li>presented by lambert in 1772</li>
072     * </ul>
073     *
074     * <p>
075     * The difference to orthographic and stereographic projection, comes from the spacing between the parallels. The space
076     * decreases with increasing distance from the pole. The opposite pole not visible on either the orthographic or
077     * stereographic may be shown on the lambert as a large circle surrounding the map, almost half again as far as the
078     * equator from the center. Normally the projectction is not shown beyond one hemisphere (or beyond the equator in the
079     * polar aspect).
080     * </p>
081     *
082     * <p>
083     * It is known to be used by following epsg transformations:
084     * <ul>
085     * <li>EPSG:3035</li>
086     * </ul>
087     * </p>
088     *
089     * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
090     *
091     * @author last edited by: $Author: mschneider $
092     *
093     * @version $Revision: 18195 $, $Date: 2009-06-18 17:55:39 +0200 (Do, 18 Jun 2009) $
094     *
095     */
096    
097    public class LambertAzimuthalEqualArea extends AzimuthalProjection {
098    
099        private static final long serialVersionUID = -5805086267437551594L;
100    
101        private double sinb1;
102    
103        private double cosb1;
104    
105        /**
106         * qp is q (needed for authalicLatitude Snyder 3-12) evaluated for a phi of 90°.
107         */
108        private double qp;
109    
110        /**
111         * Will hold the value D (A slide adjustment for the standardpoint, to achieve a correct scale in all directions at
112         * the center of the projection) calculated by Snyder (24-20).
113         */
114        private double dd;
115    
116        /**
117         * Radius for the sphere having the same surface area as the ellipsoid. Calculated with Snyder (3-13).
118         */
119        private double rq;
120    
121        /**
122         * precalculated series values to calculate the authalic latitude value from.
123         */
124        private double[] apa;
125    
126        /**
127         * A variable to hold a precalculated value for the x parameter of the oblique projection
128         */
129        private double xMultiplyForward;
130    
131        /**
132         * A variable to hold a precalculated value for the y parameter of the oblique or equatorial projection
133         */
134        private double yMultiplyForward;
135    
136        /**
137         * @param geographicCRS
138         * @param falseNorthing
139         * @param falseEasting
140         * @param naturalOrigin
141         * @param units
142         * @param scale
143         * @param id
144         *            an identifiable instance containing information about this projection
145         */
146        public LambertAzimuthalEqualArea( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
147                                          Point2d naturalOrigin, Unit units, double scale, Identifiable id ) {
148            super( geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale, false/* not conformal */,
149                   true/* equals-area */, id );
150            if ( !isSpherical() ) {
151                // sin(rad(90)) = 1;
152                qp = calcQForAuthalicLatitude( 1., getEccentricity() );
153                rq = getSemiMajorAxis() * Math.sqrt( .5 * qp );// Snyder (3-13)
154                apa = getAuthalicLatitudeSeriesValues( getSquaredEccentricity() );
155    
156                switch ( getMode() ) {
157                case NORTH_POLE:
158                case SOUTH_POLE:
159                    xMultiplyForward = getSemiMajorAxis();
160                    yMultiplyForward = ( getMode() == NORTH_POLE ) ? -getSemiMajorAxis() : getSemiMajorAxis();
161                    dd = 1.;
162                    break;
163                case EQUATOR:
164                    dd = 1. / ( rq );
165                    xMultiplyForward = getSemiMajorAxis();
166                    yMultiplyForward = getSemiMajorAxis() * .5 * qp;
167                    break;
168                case OBLIQUE:
169                    double sinphi = getSinphi0();
170                    // arcsin( q/ qp) = beta . Snyder (3-11)
171                    sinb1 = calcQForAuthalicLatitude( sinphi, getEccentricity() ) / qp;
172                    // sin*sin + cos*cos = 1
173                    cosb1 = Math.sqrt( 1. - sinb1 * sinb1 );
174                    // (24-20) D = a*m_1 / (Rq*cos(beta_1) )
175                    double m_1 = getCosphi0() / Math.sqrt( 1. - getSquaredEccentricity() * sinphi * sinphi );
176                    dd = getSemiMajorAxis() * m_1 / ( rq * cosb1 );
177                    xMultiplyForward = rq * dd;
178                    yMultiplyForward = rq / dd;
179                    break;
180                }
181            }
182        }
183    
184        /**
185         * @param geographicCRS
186         * @param falseNorthing
187         * @param falseEasting
188         * @param naturalOrigin
189         * @param units
190         * @param scale
191         */
192        public LambertAzimuthalEqualArea( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
193                                          Point2d naturalOrigin, Unit units, double scale ) {
194            this( geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale, new Identifiable( "EPSG::9820" ) );
195        }
196    
197        /**
198         * @param geographicCRS
199         * @param falseNorthing
200         * @param falseEasting
201         * @param naturalOrigin
202         * @param units
203         * @param id
204         *            an identifiable instance containing information about this projection
205         */
206        public LambertAzimuthalEqualArea( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
207                                          Point2d naturalOrigin, Unit units, Identifiable id ) {
208            this( geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, 1, id );
209        }
210    
211        /**
212         * @param geographicCRS
213         * @param falseNorthing
214         * @param falseEasting
215         * @param naturalOrigin
216         * @param units
217         */
218        public LambertAzimuthalEqualArea( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
219                                          Point2d naturalOrigin, Unit units ) {
220            this( geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, 1 );
221        }
222    
223        /*
224         * (non-Javadoc)
225         *
226         * @see org.deegree.crs.projections.Projection#doInverseProjection(double, double)
227         */
228        @Override
229        public Point2d doInverseProjection( double x, double y )
230                                throws ProjectionException {
231            Point2d lp = new Point2d( 0, 0 );
232            x -= getFalseEasting();
233            y -= getFalseNorthing();
234            // Snyder (20-18)
235            double rho = length( x, y );
236            if ( isSpherical() ) {
237                double cosC = 0;
238                double sinC = 0;
239    
240                // Snyder (20-18)
241                if ( rho * .5 > 1. ) {
242                    throw new ProjectionException( "The Y value is beyond the maximum mappable area" );
243                }
244                // lp.y = 2. * Math.asin( rho*.5 );
245                // Snyder (24-16)
246                double c = 2 * Math.asin( rho * 0.5 );
247    
248                if ( getMode() == OBLIQUE || getMode() == EQUATOR ) {
249                    sinC = Math.sin( c );
250                    cosC = Math.cos( c );
251                }
252                switch ( getMode() ) {
253                case OBLIQUE:
254                    // the theta from Snyder (20-14)
255                    lp.y = ( rho <= EPS10 ) ? getProjectionLatitude() : Math.asin( cosC * getSinphi0()
256                                                                                   + ( ( y * sinC * getCosphi0() ) / rho ) );
257    
258                    // For the calculation of the Lamda (Snyder[20-15]) proj4 obviously uses the atan2 method, I don't know
259                    // if this is correct.
260    
261                    // x = the radial coordinate (usually denoted as r) it denotes the point's distance from a central point
262                    // known as the pole (equivalent to the origin in the Cartesian system)
263                    x *= sinC * getCosphi0();
264    
265                    // y = The angular coordinate (also known as the polar angle or the azimuth angle, and usually denoted
266                    // by θ or t) denotes the positive or anticlockwise (counterclockwise) angle required to reach the point
267                    // from the 0° ray or polar axis (which is equivalent to the positive x-axis in the Cartesian coordinate
268                    // plane).
269                    y = ( cosC - Math.sin( lp.y ) * getSinphi0() ) * rho;
270                    /**
271                     * it could be something like this too.
272                     */
273                    // lp.x = Math.atan( ( x * sinC ) / ( ( rho * getCosphi0() * cosC ) - ( y * sinC * getSinphi0() ) ) );
274                    break;
275                case EQUATOR:
276                    lp.y = ( rho <= EPS10 ) ? 0. : Math.asin( y * sinC / rho );
277                    x *= sinC;
278                    y = cosC * rho;
279                    break;
280                case NORTH_POLE:
281                    y = -y;
282                    // cos(90 or -90) = 0, therefore the last term from (20-14) is null
283                    // sin( 90 or -90 = + or - 1 what is
284                    // left is asin( (-or+) cosC )
285                    // from cos(c) = sin( HALFPI - c ) follows
286                    lp.y = HALFPI - c;
287                    break;
288                case SOUTH_POLE:
289                    lp.y = c - HALFPI;
290                    break;
291                }
292                // calculation of the Lamda (Snyder[20-15]) is this correct???
293                lp.x = ( ( y == 0. && ( getMode() == EQUATOR || getMode() == OBLIQUE ) ) ? 0. : Math.atan2( x, y ) );
294            } else {
295                double q = 0;
296                double arcSinusBeta = 0;
297                switch ( getMode() ) {
298                case EQUATOR:
299                case OBLIQUE:
300                    // Snyder (p.189 24-28)
301                    x /= dd;
302                    y *= dd;
303                    rho = length( x, y );
304                    if ( rho < EPS10 ) {
305                        return new Point2d( 0, getProjectionLatitude() );
306                    }
307                    // Snyder (p.189 24-29).
308                    double ce = 2. * Math.asin( ( .5 * rho ) / rq );
309                    double cosinusCe = Math.cos( ce );
310                    double sinusCe = Math.sin( ce );
311    
312                    x *= sinusCe;
313                    if ( getMode() == OBLIQUE ) {
314                        // Snyder (p.189 24-30)
315                        arcSinusBeta = cosinusCe * sinb1 + ( ( y * sinusCe * cosb1 ) / rho );
316                        // Snyder (p.188 24-27)
317                        q = qp * ( arcSinusBeta );
318                        // calculate the angular coordinate to be used in the atan2 method.
319                        y = rho * cosb1 * cosinusCe - y * sinb1 * sinusCe;
320                    } else {
321                        arcSinusBeta = y * sinusCe / rho;
322                        q = qp * ( arcSinusBeta );
323                        y = rho * cosinusCe;
324                    }
325                    break;
326                case NORTH_POLE:
327                    y = -y;
328                case SOUTH_POLE:
329                    // will be used to calc q.
330                    q = ( x * x + y * y );
331                    if ( Math.abs( q ) < EPS10 ) {
332                        return new Point2d( 0, getProjectionLatitude() );
333                    }
334                    // Simplified Snyder (p.190 24-32), because sin(phi) = 1, the qp can be used to calc arcSinusBeta.
335                    // xMultiplyForward = getSemiMajorAxis();
336                    double aSquare = xMultiplyForward * xMultiplyForward;
337                    arcSinusBeta = 1. - ( q / ( aSquare * qp ) );
338                    if ( getMode() == SOUTH_POLE ) {
339                        arcSinusBeta = -arcSinusBeta;
340                    }
341                    break;
342                }
343                lp.x = Math.atan2( x, y );
344                lp.y = calcPhiFromAuthalicLatitude( Math.asin( arcSinusBeta ), apa );
345            }
346            lp.x += getProjectionLongitude();
347            return lp;
348        }
349    
350        /*
351         * (non-Javadoc)
352         *
353         * @see org.deegree.crs.projections.Projection#doProjection(double, double)
354         */
355        @Override
356        public Point2d doProjection( double lambda, double phi )
357                                throws ProjectionException {
358            Point2d result = new Point2d( 0, 0 );
359            lambda -= getProjectionLongitude();
360            double sinphi = Math.sin( phi );
361            double cosLamda = Math.cos( lambda );
362            double sinLambda = Math.sin( lambda );
363    
364            if ( isSpherical() ) {
365                double cosphi = Math.cos( phi );
366                double kAccent = 0;
367                switch ( getMode() ) {
368                case OBLIQUE:
369                    // Calculation of k' Snyder (24-2)
370                    kAccent = 1. + ( getSinphi0() * sinphi ) + ( getCosphi0() * cosphi * cosLamda );
371                    if ( Math.abs( kAccent ) <= EPS11 ) {
372                        throw new ProjectionException(
373                                                       "The scalefactor (k') in the perpendicular direction to the radius from the center of the map equals: "
374                                                                               + kAccent
375                                                                               + " this will cause a divide by zero." );
376                    }
377                    kAccent = Math.sqrt( 2. / kAccent );
378    
379                    result.x = kAccent * cosphi * sinLambda;
380                    result.y = kAccent * ( ( getCosphi0() * sinphi ) - ( getSinphi0() * cosphi * cosLamda ) );
381                    break;
382                case EQUATOR:
383                    kAccent = 1. + cosphi * cosLamda;
384                    if ( kAccent <= EPS11 ) {
385                        throw new ProjectionException(
386                                                       "The scalefactor (k') in the perpendicular direction to the radius from the center of the map equals: "
387                                                                               + kAccent
388                                                                               + " this will cause a divide by zero." );
389                    }
390                    kAccent = Math.sqrt( 2. / kAccent );
391                    result.x = kAccent * cosphi * sinLambda;
392                    result.y = kAccent * sinphi;
393                    break;
394                case NORTH_POLE:
395                    cosLamda = -cosLamda;
396                case SOUTH_POLE:
397                    if ( Math.abs( phi + getProjectionLatitude() ) < EPS11 ) {
398                        throw new ProjectionException( "The requested phi: " + phi + " lies on the singularity ("
399                                                       + ( ( getMode() == SOUTH_POLE ) ? "South-Pole" : "North-Pole" )
400                                                       + ") of this projection's mappable area." );
401                    }
402                    result.y = QUARTERPI - ( phi * .5 );
403                    result.y = 2. * ( getMode() == SOUTH_POLE ? Math.cos( result.y ) : Math.sin( result.y ) );
404                    result.x = result.y * sinLambda;
405                    result.y *= cosLamda;
406                    break;
407                }
408                // the radius is stil to be multiplied.
409                result.x *= getSemiMajorAxis();
410                result.y *= getSemiMajorAxis();
411            } else {
412                double sinb = 0;
413                double cosb = 0;
414    
415                // The big B of snyder (24-19), but will also be used as a place holder for the none oblique calculations.
416                double bigB = 0;
417    
418                double q = calcQForAuthalicLatitude( sinphi, getEccentricity() );
419                if ( getMode() == OBLIQUE || getMode() == EQUATOR ) {
420                    // snyder ( 3-11 )
421                    sinb = q / qp;
422                    cosb = Math.sqrt( 1. - sinb * sinb );
423                }
424                switch ( getMode() ) {
425                case OBLIQUE:
426                    // Snyder (24-19)
427                    bigB = 1. + sinb1 * sinb + cosb1 * cosb * cosLamda;
428                    break;
429                case EQUATOR:
430                    // dd, sin as well as cos(beta1) fall out Snyder (24-21).
431                    bigB = 1. + cosb * cosLamda;
432                    break;
433                case NORTH_POLE:
434                    bigB = HALFPI + phi;
435                    q = qp - q;
436                    break;
437                case SOUTH_POLE:
438                    bigB = phi - HALFPI;
439                    q = qp + q;
440                    break;
441                }
442                /**
443                 * Test to see if the projection point is 0, -> divide by zero.
444                 */
445                if ( Math.abs( bigB ) < EPS10 ) {
446                    throw new ProjectionException(
447                                                   "The projectionPoint B from the authalic latitude beta: "
448                                                                           + ( Math.toDegrees( Math.asin( sinb ) ) )
449                                                                           + "° lies on the singularity of this projection's mappable area, resulting in a divide by zero." );
450                }
451                switch ( getMode() ) {
452                case OBLIQUE:
453                    bigB = Math.sqrt( 2 / bigB );
454                    result.x = xMultiplyForward * bigB * cosb * sinLambda;
455                    result.y = yMultiplyForward * bigB * ( cosb1 * sinb - sinb1 * cosb * cosLamda );
456                    break;
457                case EQUATOR:
458                    bigB = Math.sqrt( 2 / bigB );
459                    // dd, sin as well as cosbeta1 fall out Snyder (24-21), xMulti = getSemimajorAxis(), yMulti =
460                    // getSemiMajorAxsis() * 0.5 * qp
461                    result.x = xMultiplyForward * bigB * cosb * sinLambda;
462                    result.y = yMultiplyForward * bigB * sinb;
463                    break;
464                case NORTH_POLE:
465                case SOUTH_POLE:
466                    if ( q >= 0. ) {
467                        bigB = Math.sqrt( q );
468                        // xMulti = yMulti = getSemimajorAxis()
469                        result.x = xMultiplyForward * bigB * sinLambda;
470                        // if NORTH, yMultiplyForward = -getSemiMajorAxis.
471                        result.y = yMultiplyForward * cosLamda * bigB;
472                    } else {
473                        result.x = 0;
474                        result.y = 0;
475                    }
476                    break;
477                }
478            }
479            result.x += getFalseEasting();
480            result.y += getFalseNorthing();
481            return result;
482        }
483    
484        @Override
485        public String getImplementationName() {
486            return "lambertAzimuthalEqualArea";
487        }
488    
489    }