001 //$HeadURL: $
002 /*---------------- FILE HEADER ------------------------------------------
003 This file is part of deegree.
004 Copyright (C) 2001-2008 by:
005 Department of Geography, University of Bonn
006 http://www.giub.uni-bonn.de/deegree/
007 lat/lon GmbH
008 http://www.lat-lon.de
009
010 This library is free software; you can redistribute it and/or
011 modify it under the terms of the GNU Lesser General Public
012 License as published by the Free Software Foundation; either
013 version 2.1 of the License, or (at your option) any later version.
014 This library is distributed in the hope that it will be useful,
015 but WITHOUT ANY WARRANTY; without even the implied warranty of
016 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
017 Lesser General Public License for more details.
018 You should have received a copy of the GNU Lesser General Public
019 License along with this library; if not, write to the Free Software
020 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
021 Contact:
022
023 Andreas Poth
024 lat/lon GmbH
025 Aennchenstr. 19
026 53177 Bonn
027 Germany
028 E-Mail: poth@lat-lon.de
029
030 Prof. Dr. Klaus Greve
031 Department of Geography
032 University of Bonn
033 Meckenheimer Allee 166
034 53115 Bonn
035 Germany
036 E-Mail: greve@giub.uni-bonn.de
037 ---------------------------------------------------------------------------*/
038
039 package org.deegree.crs.projections.conic;
040
041
042 import javax.vecmath.Point2d;
043
044 import org.deegree.crs.components.Unit;
045 import org.deegree.crs.coordinatesystems.GeographicCRS;
046 import org.deegree.crs.projections.ProjectionUtils;
047
048 /**
049 * The <code>LambertConformalConic</code> projection has following properties
050 * <q>(Snyder p. 104)</q>
051 * <ul>
052 * <li>Conic</li>
053 * <li>Conformal</li>
054 * <li>Parallels are unequally spaced arcs of concentric circles, more closely spaced near the center of the map</li>
055 * <li>Meridians are equally spaced radii of the same circles, thereby cutting paralles at right angles.</li>
056 * <li>Scale is true along two standard parallels, normally or along just one.</li>
057 * <li>The Pole in the same hemisphere as the standard parallels is a point; other pole is at infinity</li>
058 * <li>Used for maps of countries and regions with predominant east-west expanse.</li>
059 * <li>Presented by Lambert in 1772.</li>
060 * </ul>
061 * <p>
062 * <q>from: http://lists.maptools.org/pipermail/proj/2003-January/000592.html </q>
063 * For east-west regions, the Lambert Conformal Conic is slightly better than the Transverse Mercator because of the
064 * ability to go farther in an east-west direction and still be able to have "round-trip" transformation accuracy.
065 * Geodetically speaking, it is NOT as good as the transverse Mercator.
066 * </p>
067 * <p>
068 * It is known to be used by following epsg transformations:
069 * <ul>
070 * <li>EPSG:3034</li>
071 * </ul>
072 * </p>
073 *
074 * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
075 *
076 * @author last edited by: $Author:$
077 *
078 * @version $Revision:$, $Date:$
079 *
080 */
081
082 public class LambertConformalConic extends ConicProjection {
083
084 // Will contain snyder's variable 'n' from fromula (15-3) for the spherical projection or (15-8) for the ellipsoidal
085 // projection.
086 private double n;
087
088 private double rho0;
089
090 // Will contain snyder's variable 'F' from fromula (15-2) for the spherical projection or (15-10) for the
091 // ellipsoidal projection.
092 private double largeF;
093
094 // used for the calcualtion of phi (in the inverse projection with an ellipsoid) by applying the pre calculated
095 // values to the series of Snyder (p.15 3-5), thus avoiding iteration.
096 private double[] preCalcedPhiSeries;
097
098 /**
099 *
100 * @param firstParallelLatitude
101 * the latitiude (in radians) of the first parallel. (Snyder phi_1).
102 * @param secondParallelLatitude
103 * the latitiude (in radians) of the second parallel. (Snyder phi_2).
104 * @param geographicCRS
105 * @param falseNorthing
106 * @param falseEasting
107 * @param naturalOrigin
108 * @param units
109 * @param scale
110 * @param identifiers
111 * @param names
112 * @param versions
113 * @param descriptions
114 * @param areasOfUse
115 */
116 public LambertConformalConic( double firstParallelLatitude, double secondParallelLatitude,
117 GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
118 Point2d naturalOrigin, Unit units, double scale, String[] identifiers,
119 String[] names, String[] versions, String[] descriptions, String[] areasOfUse ) {
120 super( firstParallelLatitude,
121 secondParallelLatitude,
122 geographicCRS,
123 falseNorthing,
124 falseEasting,
125 naturalOrigin,
126 units,
127 scale,
128 true, //conformal
129 false, //not equalArea
130 identifiers,
131 names,
132 versions,
133 descriptions,
134 areasOfUse );
135
136 double cosphi, sinphi;
137 boolean secant;
138
139 // if ( Math.abs( getFirstParallelLatitude() + getSecondParallelLatitude() ) < )
140 // throw new ProjectionException();
141
142 // If only one tangential parallel is used, the firstparallelLatitude will also have the same value as the
143 // projectionLatitude, in this case the constant 'n' from Snyder will have the value sin(phi).
144 n = sinphi = Math.sin( getFirstParallelLatitude() );
145 cosphi = Math.cos( getFirstParallelLatitude() );
146 secant = Math.abs( getFirstParallelLatitude() - getSecondParallelLatitude() ) >= ProjectionUtils.EPS10;
147 if ( isSpherical() ) {
148 if ( secant ) {
149 // two parallels are used, calc snyder (p.107 15-3), else n will contain sin(firstParallelLatitude),
150 // according to Snyder (p.107 just before 15-4).
151 n = Math.log( cosphi / Math.cos( getSecondParallelLatitude() ) ) / Math.log( Math.tan( ProjectionUtils.QUARTERPI + ( .5 * getSecondParallelLatitude() ) ) / Math.tan( ProjectionUtils.QUARTERPI + ( .5 * getFirstParallelLatitude() ) ) );
152 }
153 // Snyder (p.107 15-2)
154 largeF = ( cosphi * Math.pow( Math.tan( ProjectionUtils.QUARTERPI + ( .5 * getFirstParallelLatitude() ) ), n ) ) / n;
155
156 // Snyder (p.106 15-1a) pay attention to the '-n' power term...
157 rho0 = ( Math.abs( Math.abs( getProjectionLatitude() ) - ProjectionUtils.HALFPI ) < ProjectionUtils.EPS10 ) ? 0.
158 : largeF * Math.pow( Math.tan( ProjectionUtils.QUARTERPI + ( .5 * getProjectionLatitude() ) ),
159 -n );
160 } else {
161 preCalcedPhiSeries = ProjectionUtils.preCalcedThetaSeries( getSquaredEccentricity() );
162 // Calc
163 double m1 = ProjectionUtils.calcMFromSnyder( sinphi, cosphi, getSquaredEccentricity() );
164 double t1 = ProjectionUtils.tanHalfCoLatitude( getFirstParallelLatitude(), sinphi, getEccentricity() );
165 if ( secant ) {
166 sinphi = Math.sin( getSecondParallelLatitude() );
167 cosphi = Math.cos( getSecondParallelLatitude() );
168 // Basic math, the log ( x/ y ) = log(x) - log(y) if the base is the same.
169 n = Math.log( m1 / ProjectionUtils.calcMFromSnyder( sinphi, cosphi, getSquaredEccentricity() ) );
170 n /= Math.log( t1 / ProjectionUtils.tanHalfCoLatitude( getSecondParallelLatitude(), sinphi, getEccentricity() ) );
171 }
172 // Snyder (p.108 15-10), n will contain sin(getFirstLatitudePhi()) if only a tangential cone is used.
173 largeF = ( m1 * Math.pow( t1, -n ) ) / n;
174
175 // Snyder (p.108 15-7). or 0 if the projectionlatitude is on one of the poles e.g.± pi*0.5.
176 rho0 = ( Math.abs( Math.abs( getProjectionLatitude() ) - ProjectionUtils.HALFPI ) < ProjectionUtils.EPS10 ) ? 0.
177 : largeF * Math.pow( ProjectionUtils.tanHalfCoLatitude( getProjectionLatitude(),
178 getSinphi0(),
179 getEccentricity() ),
180 n );
181
182 }
183 }
184
185 /**
186 *
187 * @param firstParallelLatitude
188 * the latitiude (in radians) of the first parallel. (Snyder phi_1).
189 * @param secondParallelLatitude
190 * the latitiude (in radians) of the second parallel. (Snyder phi_2).
191 * @param geographicCRS
192 * @param falseNorthing
193 * @param falseEasting
194 * @param naturalOrigin
195 * @param units
196 * @param scale
197 * @param identifier
198 * @param name
199 * @param version
200 * @param description
201 * @param areaOfUse
202 */
203 public LambertConformalConic( double firstParallelLatitude, double secondParallelLatitude,
204 GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
205 Point2d naturalOrigin, Unit units, double scale, String identifier,
206 String name, String version, String description, String areaOfUse ) {
207 this( firstParallelLatitude,
208 secondParallelLatitude,
209 geographicCRS,
210 falseNorthing,
211 falseEasting,
212 naturalOrigin,
213 units,
214 scale,
215 new String[]{identifier},
216 new String[]{name},
217 new String[]{version},
218 new String[]{description},
219 new String[]{areaOfUse} );
220 }
221
222 /**
223 * Creates a Lambert Conformal projection with a intersecting cone at the given parallel latitudes.
224 *
225 * @param firstParallelLatitude
226 * the latitiude (in radians) of the first parallel. (Snyder phi_1).
227 * @param secondParallelLatitude
228 * the latitiude (in radians) of the second parallel. (Snyder phi_2).
229 * @param geographicCRS
230 * @param falseNorthing
231 * @param falseEasting
232 * @param naturalOrigin
233 * @param units
234 * @param scale
235 * @param identifiers
236 */
237 public LambertConformalConic( double firstParallelLatitude, double secondParallelLatitude,
238 GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
239 Point2d naturalOrigin, Unit units, double scale, String[] identifiers ) {
240 this( firstParallelLatitude,
241 secondParallelLatitude,
242 geographicCRS,
243 falseNorthing,
244 falseEasting,
245 naturalOrigin,
246 units,
247 scale,
248 identifiers, null,null,null,null );
249 }
250
251 /**
252 * Creates a Lambert Conformal projection with a intersecting cone at the given parallel latitudes.
253 *
254 * @param firstParallelLatitude
255 * the latitiude (in radians) of the first parallel. (Snyder phi_1).
256 * @param secondParallelLatitude
257 * the latitiude (in radians) of the second parallel. (Snyder phi_2).
258 * @param geographicCRS
259 * @param falseNorthing
260 * @param falseEasting
261 * @param naturalOrigin
262 * @param units
263 * @param scale
264 * @param identifier
265 */
266 public LambertConformalConic( double firstParallelLatitude, double secondParallelLatitude,
267 GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
268 Point2d naturalOrigin, Unit units, double scale, String identifier ) {
269 this( firstParallelLatitude,
270 secondParallelLatitude,
271 geographicCRS,
272 falseNorthing,
273 falseEasting,
274 naturalOrigin,
275 units,
276 scale,
277 new String[]{identifier} );
278 }
279
280 /**
281 * Creates a Lambert Conformal projection with a tangential cone at the naturalOrigin.y's latitude.
282 *
283 * @param geographicCRS
284 * @param falseNorthing
285 * @param falseEasting
286 * @param naturalOrigin
287 * @param units
288 * @param scale
289 * @param identifiers
290 * @param name
291 */
292 public LambertConformalConic( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
293 Point2d naturalOrigin, Unit units, double scale, String[] identifiers ) {
294 this( Double.NaN,
295 Double.NaN,
296 geographicCRS,
297 falseNorthing,
298 falseEasting,
299 naturalOrigin,
300 units,
301 scale,
302 identifiers );
303 }
304
305 /**
306 * Creates a Lambert Conformal projection with a tangential cone at the naturalOrigin.y's latitude.
307 *
308 * @param geographicCRS
309 * @param falseNorthing
310 * @param falseEasting
311 * @param naturalOrigin
312 * @param units
313 * @param scale
314 * @param identifier
315 * @param name
316 */
317 public LambertConformalConic( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
318 Point2d naturalOrigin, Unit units, double scale, String identifier ) {
319 this( Double.NaN,
320 Double.NaN,
321 geographicCRS,
322 falseNorthing,
323 falseEasting,
324 naturalOrigin,
325 units,
326 scale,
327 identifier );
328 }
329
330 /**
331 * Creates a Lambert Conformal projection with a intersecting cone at the given parallel latitudes. and a scale of
332 * 1.
333 *
334 * @param firstParallelLatitude
335 * the latitiude (in radians) of the first parallel. (Snyder phi_1).
336 * @param secondParallelLatitude
337 * the latitiude (in radians) of the second parallel. (Snyder phi_2).
338 * @param geographicCRS
339 * @param falseNorthing
340 * @param falseEasting
341 * @param naturalOrigin
342 * @param units
343 * @param identifiers
344 * @param name
345 */
346 public LambertConformalConic( double firstParallelLatitude, double secondParallelLatitude,
347 GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
348 Point2d naturalOrigin, Unit units, String[] identifiers ) {
349 this( firstParallelLatitude,
350 secondParallelLatitude,
351 geographicCRS,
352 falseNorthing,
353 falseEasting,
354 naturalOrigin,
355 units,
356 1,
357 identifiers );
358
359 }
360
361 /**
362 * Creates a Lambert Conformal projection with a intersecting cone at the given parallel latitudes. and a scale of
363 * 1.
364 *
365 * @param firstParallelLatitude
366 * the latitiude (in radians) of the first parallel. (Snyder phi_1).
367 * @param secondParallelLatitude
368 * the latitiude (in radians) of the second parallel. (Snyder phi_2).
369 * @param geographicCRS
370 * @param falseNorthing
371 * @param falseEasting
372 * @param naturalOrigin
373 * @param units
374 * @param identifier
375 * @param name
376 */
377 public LambertConformalConic( double firstParallelLatitude, double secondParallelLatitude,
378 GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
379 Point2d naturalOrigin, Unit units, String identifier ) {
380 this( firstParallelLatitude,
381 secondParallelLatitude,
382 geographicCRS,
383 falseNorthing,
384 falseEasting,
385 naturalOrigin,
386 units,
387 1,
388 identifier);
389
390 }
391
392 /**
393 * Creates a Lambert Conformal projection with a tangential cone at the naturalOrigin.y's latitude. And a scale of
394 * 1.
395 *
396 * @param geographicCRS
397 * @param falseNorthing
398 * @param falseEasting
399 * @param naturalOrigin
400 * @param units
401 * @param identifiers
402 * @param name
403 */
404 public LambertConformalConic( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
405 Point2d naturalOrigin, Unit units, String[] identifiers) {
406 this( Double.NaN,
407 Double.NaN,
408 geographicCRS,
409 falseNorthing,
410 falseEasting,
411 naturalOrigin,
412 units,
413 1,
414 identifiers );
415
416 }
417
418 /**
419 * Creates a Lambert Conformal projection with a tangential cone at the naturalOrigin.y's latitude. And a scale of
420 * 1.
421 *
422 * @param geographicCRS
423 * @param falseNorthing
424 * @param falseEasting
425 * @param naturalOrigin
426 * @param units
427 * @param identifier
428 * @param name
429 */
430 public LambertConformalConic( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
431 Point2d naturalOrigin, Unit units, String identifier ) {
432 this( Double.NaN,
433 Double.NaN,
434 geographicCRS,
435 falseNorthing,
436 falseEasting,
437 naturalOrigin,
438 units,
439 1,
440 identifier );
441
442 }
443
444 /**
445 *
446 * @see org.deegree.crs.projections.Projection#doInverseProjection(double, double)
447 */
448 @Override
449 public Point2d doInverseProjection( double x, double y ) {
450 Point2d out = new Point2d( 0, 0 );
451 x -= getFalseEasting();
452 y -= getFalseNorthing();
453 // why divide by the scale????
454 x /= getScaleFactor();
455 y = rho0 - (y / getScaleFactor());
456 double rho = ProjectionUtils.length( x, y );
457 if ( rho > ProjectionUtils.EPS11 ) {
458 if ( n < 0.0 ) {
459 // if using the atan2 the values must be inverted.
460 rho = -rho;
461 x = -x;
462 y = -y;
463 }
464 if ( isSpherical() ) {
465 // Snyder (p.107 15-5).
466 out.y = ( 2.0 * Math.atan( Math.pow( largeF / rho, 1.0 / n ) ) ) - ProjectionUtils.HALFPI;
467 } else {
468 // out.y = MapUtils.phi2( Math.pow( rho / largeF, 1.0 / n ), getEccentricity() );
469 double t = Math.pow( rho / largeF, 1.0 / n );
470 double chi = ProjectionUtils.HALFPI - ( 2 * Math.atan( t ) );
471 out.y = ProjectionUtils.calcPhiFromConformalLatitude( chi, preCalcedPhiSeries );
472 }
473 // Combine Snyder (P.107/109 14-9) with (p.107/109 14-11), please pay attention to the remark of snyder on
474 // the atan2 at p.107!!!
475 out.x = Math.atan2( x, y ) / n;
476 } else {
477 out.x = 0;
478 out.y = ( n > 0.0 ) ? ProjectionUtils.HALFPI : -ProjectionUtils.HALFPI;
479 }
480 out.x += getProjectionLongitude();
481 return out;
482 }
483
484 /**
485 *
486 * @see org.deegree.crs.projections.Projection#doProjection(double, double)
487 */
488 @Override
489 public Point2d doProjection( double lambda, double phi ) {
490 lambda -= getProjectionLongitude();
491 double rho = 0;
492 if ( Math.abs( Math.abs( phi ) - ProjectionUtils.HALFPI ) > ProjectionUtils.EPS10 ) {
493 // For spherical see Snyder (p.106 15-1) for ellipitical Snyder (p.108 15-7), pay attention to the '-n'
494 rho = largeF * ( isSpherical() ? Math.pow( Math.tan( ProjectionUtils.QUARTERPI + ( .5 * phi ) ), -n )
495 : Math.pow( ProjectionUtils.tanHalfCoLatitude( phi,
496 Math.sin( phi ),
497 getEccentricity() ), n ) );
498 }
499 // calc theta Snyder (p.106/108 14-4) multiply lambda with the 'n' constant.
500 double theta = lambda * n;
501
502 Point2d out = new Point2d( 0, 0 );
503 out.x = getScaleFactor() * ( rho * Math.sin( theta ) ) + getFalseEasting();
504 out.y = getScaleFactor() * ( rho0 - ( rho * Math.cos( theta ) ) ) + getFalseNorthing();
505 return out;
506 }
507
508 @Override
509 public String getDeegreeSpecificName() {
510 return "lambertConformalConic";
511 }
512
513 }