001    //$HeadURL: svn+ssh://jwilden@svn.wald.intevation.org/deegree/base/branches/2.5_testing/src/org/deegree/crs/projections/conic/LambertConformalConic.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.conic;
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.calcMFromSnyder;
044    import static org.deegree.crs.projections.ProjectionUtils.calcPhiFromConformalLatitude;
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    
055    /**
056     * The <code>LambertConformalConic</code> projection has following properties <q>(Snyder p. 104)</q>
057     * <ul>
058     * <li>Conic</li>
059     * <li>Conformal</li>
060     * <li>Parallels are unequally spaced arcs of concentric circles, more closely spaced near the center of the map</li>
061     * <li>Meridians are equally spaced radii of the same circles, thereby cutting paralles at right angles.</li>
062     * <li>Scale is true along two standard parallels, normally or along just one.</li>
063     * <li>The Pole in the same hemisphere as the standard parallels is a point; other pole is at infinity</li>
064     * <li>Used for maps of countries and regions with predominant east-west expanse.</li>
065     * <li>Presented by Lambert in 1772.</li>
066     * </ul>
067     * <p>
068     * <q>from: http://lists.maptools.org/pipermail/proj/2003-January/000592.html</q>
069     * For east-west regions, the Lambert Conformal Conic is slightly better than the Transverse Mercator because of the
070     * ability to go farther in an east-west direction and still be able to have "round-trip" transformation accuracy.
071     * Geodetically speaking, it is NOT as good as the transverse Mercator.
072     * </p>
073     * <p>
074     * It is known to be used by following epsg transformations:
075     * <ul>
076     * <li>EPSG:3034</li>
077     * </ul>
078     * </p>
079     *
080     * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
081     *
082     * @author last edited by: $Author: mschneider $
083     *
084     * @version $Revision: 18195 $, $Date: 2009-06-18 17:55:39 +0200 (Do, 18 Jun 2009) $
085     *
086     */
087    
088    public class LambertConformalConic extends ConicProjection {
089    
090        private static final long serialVersionUID = 2179059901890738599L;
091    
092        /**
093         * Will contain snyder's variable 'n' from formula (15-3) for the spherical projection or (15-8) for the ellipsoidal
094         * projection.
095         */
096        private double n;
097    
098        /**
099         * Snyder (p.108 15-7). or 0 if the projectionlatitude is on one of the poles e.g.± pi*0.5.
100         */
101        private final double rho0;
102    
103        /**
104         * Will contain snyder's variable 'F' from formula (15-2) for the spherical projection or (15-10) for the
105         * ellipsoidal projection.
106         */
107        private final double largeF;
108    
109        /**
110         * used for the calculation of phi (in the inverse projection with an ellipsoid) by applying the pre calculated
111         * values to the series of Snyder (p.15 3-5), thus avoiding iteration.
112         */
113        private double[] preCalcedPhiSeries;
114    
115        /**
116         *
117         * @param firstParallelLatitude
118         *            the latitude (in radians) of the first parallel. (Snyder phi_1).
119         * @param secondParallelLatitude
120         *            the latitude (in radians) of the second parallel. (Snyder phi_2).
121         * @param geographicCRS
122         * @param falseNorthing
123         * @param falseEasting
124         * @param naturalOrigin
125         * @param units
126         * @param scale
127         * @param id
128         *            an identifiable instance containing information about this projection
129         */
130        public LambertConformalConic( double firstParallelLatitude, double secondParallelLatitude,
131                                      GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
132                                      Point2d naturalOrigin, Unit units, double scale, Identifiable id ) {
133            super( firstParallelLatitude, secondParallelLatitude, geographicCRS, falseNorthing, falseEasting,
134                   naturalOrigin, units, scale, true/* conformal */, false /* not equalArea */, id );
135    
136            double cosphi, sinphi;
137            boolean secant;
138    
139            // If only one tangential parallel is used, the firstparallelLatitude will also have the same value as the
140            // projectionLatitude, in this case the constant 'n' from Snyder will have the value sin(phi).
141            n = sinphi = Math.sin( getFirstParallelLatitude() );
142            cosphi = Math.cos( getFirstParallelLatitude() );
143            secant = Math.abs( getFirstParallelLatitude() - getSecondParallelLatitude() ) >= EPS10;
144            if ( isSpherical() ) {
145                if ( secant ) {
146                    // two parallels are used, calc snyder (p.107 15-3), else n will contain sin(firstParallelLatitude),
147                    // according to Snyder (p.107 just before 15-4).
148                    n = Math.log( cosphi / Math.cos( getSecondParallelLatitude() ) )
149                        / Math.log( Math.tan( QUARTERPI + ( .5 * getSecondParallelLatitude() ) )
150                                    / Math.tan( QUARTERPI + ( .5 * getFirstParallelLatitude() ) ) );
151                }
152                // Snyder (p.107 15-2)
153                largeF = ( cosphi * Math.pow( Math.tan( QUARTERPI + ( .5 * getFirstParallelLatitude() ) ), n ) ) / n;
154    
155                // Snyder (p.106 15-1a) pay attention to the '-n' power term...
156                rho0 = ( Math.abs( Math.abs( getProjectionLatitude() ) - HALFPI ) < EPS10 ) ? 0.
157                                                                                           : largeF
158                                                                                             * Math.pow(
159                                                                                                         Math.tan( QUARTERPI
160                                                                                                                   + ( .5 * getProjectionLatitude() ) ),
161                                                                                                         -n );
162            } else {
163                preCalcedPhiSeries = preCalcedThetaSeries( getSquaredEccentricity() );
164                // Calc
165                double m1 = calcMFromSnyder( sinphi, cosphi, getSquaredEccentricity() );
166                double t1 = tanHalfCoLatitude( getFirstParallelLatitude(), sinphi, getEccentricity() );
167                if ( secant ) {
168                    sinphi = Math.sin( getSecondParallelLatitude() );
169                    cosphi = Math.cos( getSecondParallelLatitude() );
170                    // Basic math, the log ( x/ y ) = log(x) - log(y) if the base is the same.
171                    n = Math.log( m1 / calcMFromSnyder( sinphi, cosphi, getSquaredEccentricity() ) );
172                    n /= Math.log( t1 / tanHalfCoLatitude( getSecondParallelLatitude(), sinphi, getEccentricity() ) );
173                }
174                // Snyder (p.108 15-10), n will contain sin(getFirstLatitudePhi()) if only a tangential cone is used.
175                largeF = ( m1 * Math.pow( t1, -n ) ) / n;
176    
177                // Snyder (p.108 15-7). or 0 if the projectionlatitude is on one of the poles e.g.± pi*0.5.
178                rho0 = ( Math.abs( Math.abs( getProjectionLatitude() ) - HALFPI ) < EPS10 ) ? 0.
179                                                                                           : largeF
180                                                                                             * Math.pow(
181                                                                                                         tanHalfCoLatitude(
182                                                                                                                            getProjectionLatitude(),
183                                                                                                                            getSinphi0(),
184                                                                                                                            getEccentricity() ),
185                                                                                                         n );
186    
187            }
188        }
189    
190        /**
191         *
192         * @param firstParallelLatitude
193         *            the latitude (in radians) of the first parallel. (Snyder phi_1).
194         * @param secondParallelLatitude
195         *            the latitude (in radians) of the second parallel. (Snyder phi_2).
196         * @param geographicCRS
197         * @param falseNorthing
198         * @param falseEasting
199         * @param naturalOrigin
200         * @param units
201         * @param scale
202         */
203        public LambertConformalConic( double firstParallelLatitude, double secondParallelLatitude,
204                                      GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
205                                      Point2d naturalOrigin, Unit units, double scale ) {
206            this( firstParallelLatitude, secondParallelLatitude, geographicCRS, falseNorthing, falseEasting, naturalOrigin,
207                  units, scale, new Identifiable( "EPSG::9802" ) );
208        }
209    
210        /**
211         * Creates a Lambert Conformal projection with a tangential cone at the naturalOrigin.y's latitude.
212         *
213         * @param geographicCRS
214         * @param falseNorthing
215         * @param falseEasting
216         * @param naturalOrigin
217         * @param units
218         * @param scale
219         * @param id
220         *            an identifiable instance containing information about this projection
221         */
222        public LambertConformalConic( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
223                                      Point2d naturalOrigin, Unit units, double scale, Identifiable id ) {
224            this( Double.NaN, Double.NaN, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale, id );
225        }
226    
227        /**
228         * Creates a Lambert Conformal projection with a tangential cone at the naturalOrigin.y's latitude.
229         *
230         * @param geographicCRS
231         * @param falseNorthing
232         * @param falseEasting
233         * @param naturalOrigin
234         * @param units
235         * @param scale
236         */
237        public LambertConformalConic( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
238                                      Point2d naturalOrigin, Unit units, double scale ) {
239            this( Double.NaN, Double.NaN, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale );
240        }
241    
242        /**
243         * Creates a Lambert Conformal projection with a intersecting cone at the given parallel latitudes. and a scale of
244         * 1.
245         *
246         * @param firstParallelLatitude
247         *            the latitude (in radians) of the first parallel. (Snyder phi_1).
248         * @param secondParallelLatitude
249         *            the latitude (in radians) of the second parallel. (Snyder phi_2).
250         * @param geographicCRS
251         * @param falseNorthing
252         * @param falseEasting
253         * @param naturalOrigin
254         * @param units
255         * @param id
256         *            an identifiable instance containing information about this projection
257         */
258        public LambertConformalConic( double firstParallelLatitude, double secondParallelLatitude,
259                                      GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
260                                      Point2d naturalOrigin, Unit units, Identifiable id ) {
261            this( firstParallelLatitude, secondParallelLatitude, geographicCRS, falseNorthing, falseEasting, naturalOrigin,
262                  units, 1., id );
263        }
264    
265        /**
266         * Creates a Lambert Conformal projection with a intersecting cone at the given parallel latitudes. and a scale of
267         * 1.
268         *
269         * @param firstParallelLatitude
270         *            the latitude (in radians) of the first parallel. (Snyder phi_1).
271         * @param secondParallelLatitude
272         *            the latitude (in radians) of the second parallel. (Snyder phi_2).
273         * @param geographicCRS
274         * @param falseNorthing
275         * @param falseEasting
276         * @param naturalOrigin
277         * @param units
278         */
279        public LambertConformalConic( double firstParallelLatitude, double secondParallelLatitude,
280                                      GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
281                                      Point2d naturalOrigin, Unit units ) {
282            this( firstParallelLatitude, secondParallelLatitude, geographicCRS, falseNorthing, falseEasting, naturalOrigin,
283                  units, 1. );
284        }
285    
286        /**
287         * Creates a Lambert Conformal projection with a tangential cone at the naturalOrigin.y's latitude. And a scale of
288         * 1.
289         *
290         * @param geographicCRS
291         * @param falseNorthing
292         * @param falseEasting
293         * @param naturalOrigin
294         * @param units
295         * @param id
296         *            an identifiable instance containing information about this projection
297         */
298        public LambertConformalConic( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
299                                      Point2d naturalOrigin, Unit units, Identifiable id ) {
300            this( Double.NaN, Double.NaN, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, 1, id );
301        }
302    
303        /**
304         * Creates a Lambert Conformal projection with a tangential cone at the naturalOrigin.y's latitude. And a scale of
305         * 1.
306         *
307         * @param geographicCRS
308         * @param falseNorthing
309         * @param falseEasting
310         * @param naturalOrigin
311         * @param units
312         */
313        public LambertConformalConic( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
314                                      Point2d naturalOrigin, Unit units ) {
315            this( Double.NaN, Double.NaN, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, 1 );
316        }
317    
318        /**
319         *
320         * @see org.deegree.crs.projections.Projection#doInverseProjection(double, double)
321         */
322        @Override
323        public Point2d doInverseProjection( double x, double y ) {
324            Point2d out = new Point2d( 0, 0 );
325            x -= getFalseEasting();
326            y -= getFalseNorthing();
327            // why divide by the scale????
328            x /= getScaleFactor();
329            y = rho0 - ( y / getScaleFactor() );
330            double rho = length( x, y );
331            if ( rho > EPS11 ) {
332                if ( n < 0.0 ) {
333                    // if using the atan2 the values must be inverted.
334                    rho = -rho;
335                    x = -x;
336                    y = -y;
337                }
338                if ( isSpherical() ) {
339                    // Snyder (p.107 15-5).
340                    out.y = ( 2.0 * Math.atan( Math.pow( largeF / rho, 1.0 / n ) ) ) - HALFPI;
341                } else {
342                    // out.y = MapUtils.phi2( Math.pow( rho / largeF, 1.0 / n ), getEccentricity() );
343                    double t = Math.pow( rho / largeF, 1.0 / n );
344                    double chi = HALFPI - ( 2 * Math.atan( t ) );
345                    out.y = calcPhiFromConformalLatitude( chi, preCalcedPhiSeries );
346                }
347                // Combine Snyder (P.107/109 14-9) with (p.107/109 14-11), please pay attention to the remark of snyder on
348                // the atan2 at p.107!!!
349                out.x = Math.atan2( x, y ) / n;
350            } else {
351                out.x = 0;
352                out.y = ( n > 0.0 ) ? HALFPI : -HALFPI;
353            }
354            out.x += getProjectionLongitude();
355            return out;
356        }
357    
358        /**
359         *
360         * @see org.deegree.crs.projections.Projection#doProjection(double, double)
361         */
362        @Override
363        public Point2d doProjection( double lambda, double phi ) {
364            lambda -= getProjectionLongitude();
365            double rho = 0;
366            if ( Math.abs( Math.abs( phi ) - HALFPI ) > EPS10 ) {
367                // For spherical see Snyder (p.106 15-1) for ellipitical Snyder (p.108 15-7), pay attention to the '-n'
368                rho = largeF
369                      * ( isSpherical() ? Math.pow( Math.tan( QUARTERPI + ( .5 * phi ) ), -n )
370                                       : Math.pow( tanHalfCoLatitude( phi, Math.sin( phi ), getEccentricity() ), n ) );
371            }
372            // calc theta Snyder (p.106/108 14-4) multiply lambda with the 'n' constant.
373            double theta = lambda * n;
374    
375            Point2d out = new Point2d( 0, 0 );
376            out.x = getScaleFactor() * ( rho * Math.sin( theta ) ) + getFalseEasting();
377            out.y = getScaleFactor() * ( rho0 - ( rho * Math.cos( theta ) ) ) + getFalseNorthing();
378            return out;
379        }
380    
381        @Override
382        public String getImplementationName() {
383            return "lambertConformalConic";
384        }
385    
386        @Override
387        public boolean equals( Object other ) {
388            if ( other != null && other instanceof LambertConformalConic ) {
389                final LambertConformalConic that = (LambertConformalConic) other;
390                return super.equals( that ) && ( Math.abs( this.n - that.n ) < EPS11 )
391                       && ( Math.abs( this.largeF - that.largeF ) < EPS11 ) && ( Math.abs( this.rho0 - that.rho0 ) < EPS11 );
392            }
393            return false;
394        }
395    
396    }