001 //$HeadURL: https://svn.wald.intevation.org/svn/deegree/base/branches/2.3_testing/src/org/deegree/crs/projections/azimuthal/StereographicAzimuthal.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.calcPhiFromConformalLatitude;
044 import static org.deegree.crs.projections.ProjectionUtils.conformalLatitudeInnerPart;
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 import org.deegree.crs.exceptions.ProjectionException;
055
056 /**
057 * The <code>StereographicAzimuthal</code> class allows for Stereographic Projections of the Poles, equator as well as
058 * oblique. This projection has following properties (Snyder p. 154):
059 * <ul>
060 * <li>Azimuthal</li>
061 * <li>Conformal</li>
062 * <li>The central meridian an a particular parallel (if shown) are straight lines.</li>
063 * <li>All meridians on the polar aspects and the equator on the equatorial aspect are straight lines.</li>
064 * <li>All other meidians and parallels are shown as arcs of circles</li>
065 * <li>A perspective projections for the sphere</li>
066 * <li>Directions from the center of the projection are true (except on ellipsoidal oblique and equatorial aspects)</li>
067 * <li>Scale increases away from the center of the projection</li>
068 * <li>Point opposite the center of the projection cannot be plotted</li>
069 * <li>Used for polar maps and miscellaneous special maps</li>
070 * </ul>
071 *
072 * <p>
073 * Like Orthographic, the stereographic projection is a true perspective in its isSpherical() form. It is the only known
074 * true perspective projection of any kind that is also conformal. Its point of projection is on the the surface of the
075 * sphere at a point jus opposite the oint of tangency of the plane or the center point of the projection. Thus, if the
076 * north pole is the center of the map, the projection is from the south-pole.
077 * </p>
078 * <p>
079 * It is known to be used by following epsg transformations:
080 * <ul>
081 * <li>EPSG:28992</li>
082 * </ul>
083 * </p>
084 *
085 * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
086 *
087 * @author last edited by: $Author: rbezema $
088 *
089 * @version $Revision: 18388 $, $Date: 2009-07-08 15:50:37 +0200 (Mi, 08. Jul 2009) $
090 *
091 */
092
093 public class StereographicAzimuthal extends AzimuthalProjection {
094
095 private static final long serialVersionUID = 5399110969291553925L;
096
097 /**
098 * This variable will hold different values, for the ellipsoidal projection:
099 * <ul>
100 * <li>for Oblique projection it is the first Part (2*a*k_0*m_1, hence it's name) of A (p.160 21-27)</li>
101 * <li>for the north-south it may hold the rho of (21-33) or (21-34)
102 * <li>for the equatorial it holds 2*scale.</li>
103 * <li>For the
104 * </ul>
105 * For the spherical projection:
106 * <ul>
107 * <li>For oblique and equatorial: 2*scale, which is compliant with 2*k_0 from Snyder (p.157 21-4)</li>
108 * <li>For north and south: cos(truelat) / tan( pi/4 - phi/2) ????</li>
109 * </ul>
110 */
111 private double akm1;
112
113 private double trueScaleLatitude;
114
115 private double[] preCalcedPhiSeries;
116
117 // use for the OBLIQUE projection instead of the projectionLatitude
118 private double conformalLatitude = 0;
119
120 // sine of the conformal latitude
121 private double sinCL = 0;
122
123 // cosine of the conformal latitude
124 private double cosCL = 0;
125
126 /**
127 * @param trueScaleLatitude
128 * the latitude (in radians) of a circle around the projection point, which contains the true scale.
129 * @param geographicCRS
130 * @param falseNorthing
131 * @param falseEasting
132 * @param naturalOrigin
133 * @param units
134 * @param scale
135 * @param id
136 * an identifiable instance containing information about this projection
137 */
138 public StereographicAzimuthal( double trueScaleLatitude, GeographicCRS geographicCRS, double falseNorthing,
139 double falseEasting, Point2d naturalOrigin, Unit units, double scale, Identifiable id ) {
140 super( geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale, true,
141 false/* not equal area */, id );
142
143 preCalcedPhiSeries = preCalcedThetaSeries( getSquaredEccentricity() );
144
145 this.trueScaleLatitude = trueScaleLatitude;
146 if ( !isSpherical() ) {
147 // double X;
148 double tmp = 0;
149 switch ( getMode() ) {
150 case NORTH_POLE:
151 case SOUTH_POLE:
152 /**
153 * The t from 21-33 and 21-34 is still to be calculated.
154 */
155 // If true scale or known scale factor k_0 is to occur at the pole, the following applies:
156 if ( Double.isNaN( trueScaleLatitude )
157 || Math.abs( Math.abs( this.trueScaleLatitude ) - HALFPI ) < EPS10 ) {
158 // Snyder (p.161 21-33) akm1 = rho .
159 akm1 = 2.
160 * getScaleFactor()
161 / Math.sqrt( Math.pow( 1 + getEccentricity(), 1 + getEccentricity() )
162 * Math.pow( 1 - getEccentricity(), 1 - getEccentricity() ) );
163
164 } else {
165 // For true scale along the circle represented by trueScaleLatitude (phi_c in snyder)..
166 // akm1 will hold m_c/t_c
167 // Calculate the rho from Snyder (p.161 21-34) and place in akm1
168 tmp = Math.sin( this.trueScaleLatitude );
169
170 // First part of m of Snyder (p.160 14-15)
171 akm1 = Math.cos( this.trueScaleLatitude )
172 / tanHalfCoLatitude( this.trueScaleLatitude, tmp, getEccentricity() );
173 tmp *= getEccentricity();
174 // Second part of m of Snyder (p.160 14-15), t = (e*sin(phi_c))
175 akm1 /= Math.sqrt( 1. - tmp * tmp );
176
177 }
178 break;
179 case EQUATOR:
180 akm1 = 2. * getScaleFactor();
181 break;
182 case OBLIQUE:
183 tmp = getSinphi0();// Math.sin( getProjectionLatitude() );
184 // Calculate the ConformalLatitude for this ellipsoid Snyder (p.160 3-1) the X
185 conformalLatitude = ( 2. * Math.atan( conformalLatitudeInnerPart( getProjectionLatitude(), tmp,
186 getEccentricity() ) ) )
187 - HALFPI;
188
189 sinCL = Math.sin( conformalLatitude );
190 cosCL = Math.cos( conformalLatitude );
191
192 tmp *= getEccentricity();
193 // the first part (2a*k_0*m) of the largeA in Snyder (p.160 21-27) is stored in akm1 it still must be
194 // devided by the cos X ... etc. to get largeA
195 akm1 = 2. * getScaleFactor() * getCosphi0() / Math.sqrt( 1. - tmp * tmp );
196
197 // Setting the sinPhi0 to the conformal Latitude X.
198 // setProjectionLatitude( X );
199 break;
200 }
201 } else {
202 akm1 = 2. * getScaleFactor();
203 }
204 }
205
206 /**
207 * Sets the id to "Snyder-StereoGraphic"
208 *
209 * @param trueScaleLatitude
210 * the latitude (in radians) of a circle around the projection point, which contains the true scale.
211 * @param geographicCRS
212 * @param falseNorthing
213 * @param falseEasting
214 * @param naturalOrigin
215 * @param units
216 * @param scale
217 */
218 public StereographicAzimuthal( double trueScaleLatitude, GeographicCRS geographicCRS, double falseNorthing,
219 double falseEasting, Point2d naturalOrigin, Unit units, double scale ) {
220 this( trueScaleLatitude, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale,
221 new Identifiable( "Snyder-StereoGraphic" ) );
222 }
223
224 /**
225 * Create a {@link StereographicAzimuthal} which has a true scale latitude at MapUtils.HALFPI.
226 *
227 * @param geographicCRS
228 * @param falseNorthing
229 * @param falseEasting
230 * @param naturalOrigin
231 * @param units
232 * @param scale
233 * @param id
234 * an identifiable instance containing information about this projection
235 */
236 public StereographicAzimuthal( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
237 Point2d naturalOrigin, Unit units, double scale, Identifiable id ) {
238 this( HALFPI, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale, id );
239 }
240
241 /**
242 * Create a {@link StereographicAzimuthal} which has a true scale latitude at MapUtils.HALFPI. Sets the id to
243 * "Snyder-StereoGraphic"
244 *
245 * @param geographicCRS
246 * @param falseNorthing
247 * @param falseEasting
248 * @param naturalOrigin
249 * @param units
250 * @param scale
251 */
252 public StereographicAzimuthal( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
253 Point2d naturalOrigin, Unit units, double scale ) {
254 this( HALFPI, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale );
255 }
256
257 /**
258 * Create a {@link StereographicAzimuthal} which has a scale of 1 and a true scale latitude,
259 *
260 * @param trueScaleLatitude
261 * the latitude (in radians) of a circle around the projection point, which contains the true scale.
262 * @param geographicCRS
263 * @param falseNorthing
264 * @param falseEasting
265 * @param naturalOrigin
266 * @param units
267 * @param id
268 * an identifiable instance containing information about this projection
269 */
270 public StereographicAzimuthal( double trueScaleLatitude, GeographicCRS geographicCRS, double falseNorthing,
271 double falseEasting, Point2d naturalOrigin, Unit units, Identifiable id ) {
272 this( trueScaleLatitude, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, 1, id );
273 }
274
275 /**
276 * Create a {@link StereographicAzimuthal} which has a scale of 1 and a true scale latitude. Sets the id to
277 * "Snyder-StereoGraphic".
278 *
279 * @param trueScaleLatitude
280 * the latitude (in radians) of a circle around the projection point, which contains the true scale.
281 * @param geographicCRS
282 * @param falseNorthing
283 * @param falseEasting
284 * @param naturalOrigin
285 * @param units
286 */
287 public StereographicAzimuthal( double trueScaleLatitude, GeographicCRS geographicCRS, double falseNorthing,
288 double falseEasting, Point2d naturalOrigin, Unit units ) {
289 this( trueScaleLatitude, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, 1 );
290 }
291
292 /**
293 * Create a {@link StereographicAzimuthal} which is conformal, has a scale of 1 and a truescale latitude at pi*0.5.
294 *
295 * @param geographicCRS
296 * @param falseNorthing
297 * @param falseEasting
298 * @param naturalOrigin
299 * @param units
300 * @param id
301 * an identifiable instance containing information about this projection
302 */
303 public StereographicAzimuthal( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
304 Point2d naturalOrigin, Unit units, Identifiable id ) {
305 this( HALFPI, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, 1, id );
306 }
307
308 /**
309 * Create a {@link StereographicAzimuthal} which is conformal, has a scale of 1 and a truescale latitude at pi*0.5.
310 * Sets the id to "Snyder-StereoGraphic".
311 *
312 * @param geographicCRS
313 * @param falseNorthing
314 * @param falseEasting
315 * @param naturalOrigin
316 * @param units
317 */
318 public StereographicAzimuthal( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
319 Point2d naturalOrigin, Unit units ) {
320 this( HALFPI, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, 1 );
321 }
322
323 @Override
324 public Point2d doInverseProjection( double x, double y ) {
325 Point2d result = new Point2d( 0, 0 );
326 x -= getFalseEasting();
327 y -= getFalseNorthing();
328 double rho = length( x, y );
329 if ( isSpherical() ) {
330 // akm1 = 2*scale.
331 double c = 2. * Math.atan( rho / akm1 );
332 double sinc = Math.sin( c );
333 double cosc = Math.cos( c );
334 switch ( getMode() ) {
335 case OBLIQUE:
336 // First find theta from Snyder (p.158 20-14).
337 if ( Math.abs( rho ) <= EPS10 ) {
338 result.y = getProjectionLatitude();
339 } else {
340 result.y = Math.asin( cosc * getSinphi0() + ( ( y * sinc * getCosphi0() ) / rho ) );
341 }
342 // And now lambda but instead of using the atan, use atan2. Is this correct????
343 c = cosc - getSinphi0() * Math.sin( result.y );
344 if ( Math.abs( c ) > EPS10 || x != 0. ) {
345 result.x = Math.atan2( x * sinc * getCosphi0(), c * rho );
346 }
347 break;
348 case EQUATOR:
349 if ( Math.abs( rho ) > EPS10 ) {
350 result.y = Math.asin( y * sinc / rho );
351 }
352 if ( cosc != 0. || x != 0. ) {
353 result.x = Math.atan2( x * sinc, cosc * rho );
354 }
355 break;
356 case NORTH_POLE:
357 y = -y;
358 case SOUTH_POLE:
359 if ( Math.abs( rho ) <= EPS10 ) {
360 result.y = getProjectionLatitude();
361 } else {
362 result.y = Math.asin( getMode() == SOUTH_POLE ? -cosc : cosc );
363 }
364 result.x = Math.atan2( x, y );
365 break;
366 }
367 } else {
368 double cosCe = 0;
369 double sinCe = 0;
370 double X = 0;
371 // for oblique and equatorial projections, akm1 = first part of Snyder (p.160 21-27). which is 2*scale*m_1
372 // double ce = 2. * Math.atan2( rho * getCosphi0(), akm1 );
373 // if ( getMode() == OBLIQUE || getMode() == EQUATOR ) {
374 // cosCe = Math.cos( ce );
375 // sinCe = Math.sin( ce );
376 // }
377 double ce = 0;
378 switch ( getMode() ) {
379 case OBLIQUE:
380 ce = 2. * Math.atan2( rho * cosCL, akm1 );
381 cosCe = Math.cos( ce );
382 sinCe = Math.sin( ce );
383 X = Math.asin( cosCe * sinCL + ( ( y * sinCe * cosCL ) / rho ) );
384 result.x = Math.atan2( x * sinCe, rho * cosCL * cosCe - y * sinCL * sinCe );
385 break;
386 case EQUATOR:
387 ce = 2. * Math.atan2( rho * getCosphi0(), akm1 );
388 cosCe = Math.cos( ce );
389 sinCe = Math.sin( ce );
390 X = Math.asin( cosCe * getSinphi0() + ( ( y * sinCe * getCosphi0() ) / rho ) );
391 result.x = Math.atan2( x * sinCe, rho * getCosphi0() * cosCe - y * getSinphi0() * sinCe );
392 break;
393 case NORTH_POLE:
394 y = -y;
395 case SOUTH_POLE:
396 /**
397 * akm1 can hold rho from 21-34 or 21-33 (if no truescale latitude was given) without the parameter t.
398 * Either way it must only be multiplied with t to fulfill the formula. Independent of the true scale
399 * latitude, the value of akm1 must therefore only be inverted to calculate Snyder (p.162 21-39 or
400 * 21-40).
401 */
402 double t = rho * 1. / this.akm1;
403 X = HALFPI - 2 * Math.atan( t );
404 if ( getMode() == SOUTH_POLE ) {
405 X *= -1;
406 }
407 result.x = Math.atan2( x, y );
408 break;
409 }
410 result.y = calcPhiFromConformalLatitude( X, preCalcedPhiSeries );
411 result.x += getProjectionLongitude();
412 }
413 return result;
414 }
415
416 @Override
417 public Point2d doProjection( double lambda, double phi )
418 throws ProjectionException {
419 Point2d result = new Point2d( 0, 0 );
420 lambda -= getProjectionLongitude();
421 double cosLamda = Math.cos( lambda );
422 double sinLamda = Math.sin( lambda );
423 double sinPhi = Math.sin( phi );
424
425 if ( isSpherical() ) {
426 double cosphi = Math.cos( phi );
427
428 switch ( getMode() ) {
429 case OBLIQUE:
430 // Calc k from Snyder (p.158 21-14) and store in xy.y
431 result.y = 1. + getSinphi0() * sinPhi + getCosphi0() * cosphi * cosLamda;
432 if ( result.y <= EPS11 ) {
433 throw new ProjectionException( "Division by zero" );
434 }
435 // akm1 was set to 2*scale.
436 result.y = akm1 / result.y;
437
438 // Calc x (p.157 21-2) and y (p.157 21-3)
439 result.x = result.y * cosphi * sinLamda;
440 result.y *= getCosphi0() * sinPhi - getSinphi0() * cosphi * cosLamda;
441 break;
442 case EQUATOR:
443 // Calc k from Snyder (p.158 21-14) and store in xy.y
444 result.y = 1. + cosphi * cosLamda;
445 if ( result.y <= EPS11 ) {
446 throw new ProjectionException( "Division by zero" );
447 }
448 // akm1 was set to 2*scale.
449 result.y = akm1 / result.y;
450
451 // Calc x (p.157 21-2) and y (p.158 21-13)
452 result.x = result.y * cosphi * sinLamda;
453 result.y *= sinPhi;
454 break;
455 case NORTH_POLE:
456 cosLamda = -cosLamda;
457 phi = -phi;
458 case SOUTH_POLE:
459 if ( Math.abs( phi - HALFPI ) < EPS10 ) {
460 throw new ProjectionException( "The requested phi: " + phi + " lies on the singularity ("
461 + ( ( getMode() == SOUTH_POLE ) ? "South-Pole" : "North-Pole" )
462 + ") of this projection's mappable area." );
463 }
464 // akm1 was set to 2*scale.
465 result.y = akm1 * Math.tan( QUARTERPI + .5 * phi );
466 result.x = sinLamda * ( result.y );
467 result.y *= cosLamda;
468 break;
469 }
470 } else {
471 double sinX = 0;
472 double cosX = 0;
473 double largeA;
474
475 if ( getMode() == OBLIQUE || getMode() == EQUATOR ) {
476 // Calculate the conformal latitue Snyder (p.160 3-1).
477 double X = 2. * Math.atan( conformalLatitudeInnerPart( phi, sinPhi, getEccentricity() ) ) - HALFPI;
478 sinX = Math.sin( X );
479 cosX = Math.cos( X );
480 }
481 switch ( getMode() ) {
482 case OBLIQUE:
483 // akm1 was set to the first part of Snyder (p.160 21-27)
484 // sinPhi0 and cosPhi0 where set to the conformal latitude of the natural origin.
485 largeA = akm1 / ( cosCL * ( 1. + sinCL * sinX + cosCL * cosX * cosLamda ) );
486 result.y = largeA * ( cosCL * sinX - sinCL * cosX * cosLamda );
487 result.x = largeA * cosX;
488 break;
489 case EQUATOR:
490 // akm1 was set to 2. * scale;
491 // Calculate the A of Snyder (p.161 21-29).
492 largeA = 2. * akm1 / ( 1. + cosX * cosLamda );
493 result.y = largeA * sinX;
494 result.x = largeA * cosX;
495 break;
496 case SOUTH_POLE:
497 phi = -phi;
498 cosLamda = -cosLamda;
499 sinPhi = -sinPhi;
500 case NORTH_POLE:
501 // akm1 holds the rho's from 21-33 or 21-34, but without the t (==halfColatitude of the conformal
502 // latitude).
503 result.x = akm1 * tanHalfCoLatitude( phi, sinPhi, getEccentricity() );
504 result.y = -result.x * cosLamda;
505 break;
506 }
507 // Sinus(Lambda) was still to be multiplied on the x.
508 result.x = result.x * sinLamda;
509 }
510 result.x += getFalseEasting();
511 result.y += getFalseNorthing();
512 return result;
513 }
514
515 @Override
516 public String getImplementationName() {
517 return "stereographicAzimuthal";
518 }
519
520 /**
521 * @return the trueScaleLatitude.
522 */
523 public final double getTrueScaleLatitude() {
524 return trueScaleLatitude;
525 }
526 }