001 //$HeadURL: https://svn.wald.intevation.org/svn/deegree/base/branches/2.3_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 }