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