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