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