001 //$HeadURL: svn+ssh://rbezema@svn.wald.intevation.org/deegree/base/tags/2.1/src/org/deegree/model/csct/ct/StereographicProjection.java $
002 /*---------------- FILE HEADER ------------------------------------------
003
004 This file is part of deegree.
005 Copyright (C) 2001 by:
006 EXSE, Department of Geography, University of Bonn
007 http://www.giub.uni-bonn.de/exse/
008 lat/lon GmbH
009 http://www.lat-lon.de
010
011 It has been implemented within SEAGIS - An OpenSource implementation of OpenGIS specification
012 (C) 2001, Institut de Recherche pour le D�veloppement (http://sourceforge.net/projects/seagis/)
013 SEAGIS Contacts: Surveillance de l'Environnement Assist�e par Satellite
014 Institut de Recherche pour le D�veloppement / US-Espace
015 mailto:seasnet@teledetection.fr
016
017
018 This library is free software; you can redistribute it and/or
019 modify it under the terms of the GNU Lesser General Public
020 License as published by the Free Software Foundation; either
021 version 2.1 of the License, or (at your option) any later version.
022
023 This library is distributed in the hope that it will be useful,
024 but WITHOUT ANY WARRANTY; without even the implied warranty of
025 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
026 Lesser General Public License for more details.
027
028 You should have received a copy of the GNU Lesser General Public
029 License along with this library; if not, write to the Free Software
030 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
031
032 Contact:
033
034 Andreas Poth
035 lat/lon GmbH
036 Aennchenstr. 19
037 53115 Bonn
038 Germany
039 E-Mail: poth@lat-lon.de
040
041 Klaus Greve
042 Department of Geography
043 University of Bonn
044 Meckenheimer Allee 166
045 53115 Bonn
046 Germany
047 E-Mail: klaus.greve@uni-bonn.de
048
049
050 ---------------------------------------------------------------------------*/
051 package org.deegree.model.csct.ct;
052
053 // OpenGIS (SEAS) dependencies
054 import java.awt.geom.Point2D;
055
056 import org.deegree.model.csct.cs.Projection;
057 import org.deegree.model.csct.resources.css.ResourceKeys;
058 import org.deegree.model.csct.resources.css.Resources;
059
060 /**
061 * Projection st�r�ographique. Les directions � partir du point central sont vrais, mais
062 * les aires et les longueurs deviennent de plus en plus d�form�es � mesure que l'on
063 * s'�loigne du centre. Cette projection est utilis�e pour repr�senter des r�gions
064 * polaires. Elle peut �tre appropri�e pour d'autres r�gions ayant une forme circulaire.
065 * <br>
066 * <br>
067 *
068 * R�f�rence: John P. Snyder (Map Projections - A Working Manual, U.S. Geological Survey
069 * Professional Paper 1395, 1987)
070 *
071 * @version 1.0
072 * @author Andr� Gosselin
073 * @author Martin Desruisseaux
074 */
075 final class StereographicProjection extends PlanarProjection {
076 /**
077 * Nombre maximal d'it�rations permises lors du calcul de la projection inverse.
078 */
079 private static final int MAX_ITER = 10;
080
081 /** Projection mode for switch statement. */
082 private static final int SPHERICAL_NORTH = 0;
083
084 /** Projection mode for switch statement. */
085 private static final int SPHERICAL_SOUTH = 1;
086
087 /** Projection mode for switch statement. */
088 private static final int ELLIPSOIDAL_SOUTH = 2;
089
090 /** Projection mode for switch statement. */
091 private static final int ELLIPSOIDAL_NORTH = 3;
092
093 /** Projection mode for switch statement. */
094 private static final int SPHERICAL_OBLIQUE = 4;
095
096 /** Projection mode for switch statement. */
097 private static final int SPHERICAL_EQUATORIAL = 5;
098
099 /** Projection mode for switch statement. */
100 private static final int ELLIPSOIDAL_EQUATORIAL = 6;
101
102 /** Projection mode for switch statement. */
103 private static final int ELLIPSOIDAL_OBLIQUE = 7;
104
105 /**
106 * Projection mode. It must be one of the following constants:
107 * {@link #SPHERICAL_NORTH},{@link #SPHERICAL_SOUTH},{@link #ELLIPSOIDAL_NORTH},
108 * {@link #ELLIPSOIDAL_SOUTH}.{@link #SPHERICAL_OBLIQUE},
109 * {@link #SPHERICAL_EQUATORIAL},{@link #ELLIPSOIDAL_OBLIQUE}or
110 * {@link #ELLIPSOIDAL_EQUATORIAL}.
111 */
112 private final int mode;
113
114 /**
115 * Global scale factor. Value <code>ak0</code> is equals to
116 * <code>{@link #a}*k0</code>.
117 */
118 private final double k0, ak0;
119
120 /**
121 * Facteurs utilis�s lors des projections obliques et equatorialles.
122 */
123 private final double sinphi0, cosphi0, chi1, sinChi1, cosChi1;
124
125 /**
126 * Latitude of true scale, in radians. Used for {@link #toString}implementation.
127 */
128 private final double latitudeTrueScale;
129
130 /**
131 * Construct a new map projection from the suplied parameters.
132 *
133 * @param parameters The parameter values in standard units.
134 * @throws MissingParameterException if a mandatory parameter is missing.
135 */
136 protected StereographicProjection( final Projection parameters )
137 throws MissingParameterException {
138 this( parameters, true, true );
139 }
140
141 /**
142 * Construct a new map projection from the suplied parameters.
143 *
144 * @param parameters The parameter values in standard units.
145 * @param polar <code>true</code> for polar projection.
146 * @param auto <code>true</code> if projection (polar vs oblique) can be selected
147 * automatically.
148 * @throws MissingParameterException if a mandatory parameter is missing.
149 */
150 private StereographicProjection( final Projection parameters, final boolean polar,
151 final boolean auto ) throws MissingParameterException {
152 //////////////////////////
153 // Fetch parameters //
154 //////////////////////////
155 super( parameters );
156 final double defaultLatitude = parameters.getValue( "latitude_of_origin", polar ? 90 : 0 );
157 latitudeTrueScale = Math.abs( latitudeToRadians(
158 parameters.getValue(
159 "latitude_true_scale",
160 defaultLatitude ),
161 true ) );
162
163 //////////////////////////
164 // Compute constants //
165 //////////////////////////
166 if ( auto ? ( Math.abs( Math.abs( centralLatitude ) - ( Math.PI / 2 ) ) < EPS ) : polar ) {
167 if ( centralLatitude < 0 ) {
168 centralLatitude = -( Math.PI / 2 );
169 mode = ( isSpherical ) ? SPHERICAL_SOUTH : ELLIPSOIDAL_SOUTH;
170 } else {
171 centralLatitude = +( Math.PI / 2 );
172 mode = ( isSpherical ) ? SPHERICAL_NORTH : ELLIPSOIDAL_NORTH;
173 }
174 } else if ( Math.abs( centralLatitude ) < EPS ) {
175 centralLatitude = 0;
176 mode = ( isSpherical ) ? SPHERICAL_EQUATORIAL : ELLIPSOIDAL_EQUATORIAL;
177 } else {
178 mode = ( isSpherical ) ? SPHERICAL_OBLIQUE : ELLIPSOIDAL_OBLIQUE;
179 }
180 switch ( mode ) {
181 default: {
182 cosphi0 = Math.cos( centralLatitude );
183 sinphi0 = Math.sin( centralLatitude );
184 chi1 = 2.0 * Math.atan( ssfn( centralLatitude, sinphi0 ) ) - ( Math.PI / 2 );
185 cosChi1 = Math.cos( chi1 );
186 sinChi1 = Math.sin( chi1 );
187 break;
188 }
189 case SPHERICAL_EQUATORIAL:
190 case ELLIPSOIDAL_EQUATORIAL: {
191 cosphi0 = 1.0;
192 sinphi0 = 0.0;
193 chi1 = 0.0;
194 cosChi1 = 1.0;
195 sinChi1 = 0.0;
196 break;
197 }
198 }
199
200 //////////////////////////
201 // Compute k0 and ak0 //
202 //////////////////////////
203 switch ( mode ) {
204 default: {
205 System.out.println( "StereographicProjection: line 201: should not happen" );
206 }
207 case ELLIPSOIDAL_NORTH:
208 case ELLIPSOIDAL_SOUTH: {
209 if ( Math.abs( latitudeTrueScale - ( Math.PI / 2 ) ) >= EPS ) {
210 final double t = Math.sin( latitudeTrueScale );
211 k0 = ( Math.cos( latitudeTrueScale ) / ( Math.sqrt( 1 - es * t * t ) ) )
212 / tsfn( latitudeTrueScale, t );
213 } else {
214 // True scale at pole
215 k0 = 2.0 / Math.sqrt( Math.pow( 1 + e, 1 + e ) * Math.pow( 1 - e, 1 - e ) );
216 }
217 break;
218 }
219
220 case SPHERICAL_NORTH:
221 case SPHERICAL_SOUTH: {
222 if ( Math.abs( latitudeTrueScale - ( Math.PI / 2 ) ) >= EPS )
223 k0 = 1 + Math.sin( latitudeTrueScale );
224 else
225 k0 = 2;
226 break;
227 }
228
229 case ELLIPSOIDAL_OBLIQUE:
230 case ELLIPSOIDAL_EQUATORIAL: {
231 k0 = 2.0 * Math.cos( centralLatitude ) / Math.sqrt( 1 - es * sinphi0 * sinphi0 );
232 break;
233 }
234
235 case SPHERICAL_OBLIQUE:
236 case SPHERICAL_EQUATORIAL: {
237 k0 = 2;
238 break;
239 }
240 }
241 ak0 = a * k0;
242 }
243
244 /**
245 * Returns a human readable name localized for the specified locale.
246 */
247 public String getName( ) {
248 return Resources.getResources( null ).getString( ResourceKeys.STEREOGRAPHIC_PROJECTION );
249 }
250
251 /**
252 * Transforms the specified ( <var>x </var>, <var>y </var>) coordinate and stores the
253 * result in <code>ptDst</code>.
254 */
255 protected Point2D transform( double x, double y, final Point2D ptDst )
256 throws TransformException {
257 x -= centralMeridian;
258 final double coslat = Math.cos( y );
259 final double sinlat = Math.sin( y );
260 final double coslon = Math.cos( x );
261 final double sinlon = Math.sin( x );
262
263 switch ( mode ) {
264 default: {
265 System.out.println( "StereographicProjection: line 262: should not happen" );
266 }
267 case ELLIPSOIDAL_NORTH: {
268 final double rho = ak0 * tsfn( y, sinlat );
269 x = rho * sinlon;
270 y = -rho * coslon;
271 break;
272 }
273 case ELLIPSOIDAL_SOUTH: {
274 final double rho = ak0 * tsfn( -y, -sinlat );
275 x = rho * sinlon;
276 y = rho * coslon;
277 break;
278 }
279 case SPHERICAL_NORTH: {
280 if ( !( Math.abs( 1 + sinlat ) >= TOL ) ) {
281 throw new TransformException(
282 Resources.format( ResourceKeys.ERROR_VALUE_TEND_TOWARD_INFINITY ) );
283 }
284 // (21-8)
285 final double f = ak0 * coslat / ( 1 + sinlat );// == tan (pi/4 - phi/2)
286 x = f * sinlon; // (21-5)
287 y = -f * coslon; // (21-6)
288 break;
289 }
290 case SPHERICAL_SOUTH: {
291 if ( !( Math.abs( 1 - sinlat ) >= TOL ) ) {
292 throw new TransformException(
293 Resources.format( ResourceKeys.ERROR_VALUE_TEND_TOWARD_INFINITY ) );
294 }
295 // (21-12)
296 final double f = ak0 * coslat / ( 1 - sinlat );// == tan (pi/4 + phi/2)
297 x = f * sinlon; // (21-9)
298 y = f * coslon; // (21-10)
299 break;
300 }
301 case SPHERICAL_EQUATORIAL: {
302 double f = 1.0 + coslat * coslon;
303 if ( !( f >= TOL ) ) {
304 throw new TransformException(
305 Resources.format( ResourceKeys.ERROR_VALUE_TEND_TOWARD_INFINITY ) );
306 }
307 f = ak0 / f;
308 x = f * coslat * sinlon;
309 y = f * sinlat;
310 break;
311 }
312 case SPHERICAL_OBLIQUE: {
313 double f = 1.0 + sinphi0 * sinlat + cosphi0 * coslat * coslon; // (21-4)
314 if ( !( f >= TOL ) ) {
315 throw new TransformException(
316 Resources.format( ResourceKeys.ERROR_VALUE_TEND_TOWARD_INFINITY ) );
317 }
318 f = ak0 / f;
319 x = f * coslat * sinlon; // (21-2)
320 y = f * ( cosphi0 * sinlat - sinphi0 * coslat * coslon );// (21-3)
321 break;
322 }
323 case ELLIPSOIDAL_EQUATORIAL: {
324 final double chi = 2.0 * Math.atan( ssfn( y, sinlat ) ) - ( Math.PI / 2 );
325 final double sinChi = Math.sin( chi );
326 final double cosChi = Math.cos( chi );
327 final double A = ak0 / ( 1.0 + cosChi * coslon );
328 x = A * cosChi * sinlon;
329 y = A * sinChi;
330 break;
331 }
332 case ELLIPSOIDAL_OBLIQUE: {
333 final double chi = 2.0 * Math.atan( ssfn( y, sinlat ) ) - ( Math.PI / 2 );
334 final double sinChi = Math.sin( chi );
335 final double cosChi = Math.cos( chi );
336 final double cosChi_coslon = cosChi * coslon;
337 final double A = ak0 / cosChi1 / ( 1 + sinChi1 * sinChi + cosChi1 * cosChi_coslon );
338 x = A * cosChi * sinlon;
339 y = A * ( cosChi1 * sinChi - sinChi1 * cosChi_coslon );
340 break;
341 }
342 }
343 y += false_northing;
344 x += false_easting;
345 if ( ptDst != null ) {
346 ptDst.setLocation( x, y );
347 return ptDst;
348 }
349 return new Point2D.Double( x, y );
350 }
351
352 /**
353 * Transforms the specified ( <var>x </var>, <var>y </var>) coordinate and stores the
354 * result in <code>ptDst</code>.
355 */
356 protected Point2D inverseTransform( double x, double y, final Point2D ptDst )
357 throws TransformException {
358 x -= false_easting;
359 y -= false_northing;
360 x /= a;
361 y /= a;
362 final double rho = Math.sqrt( x * x + y * y );
363 choice: switch ( mode ) {
364 default: {
365
366 }
367 case SPHERICAL_NORTH: {
368 y = -y;
369 // fallthrough
370 }
371 case SPHERICAL_SOUTH: {
372 // (20-17) call atan2(x,y) to properly deal with y==0
373 x = ( Math.abs( x ) < TOL && Math.abs( y ) < TOL ) ? centralMeridian
374 : Math.atan2( x, y )
375 + centralMeridian;
376 if ( Math.abs( rho ) < TOL )
377 y = centralLatitude;
378 else {
379 final double c = 2.0 * Math.atan( rho / k0 );
380 final double cosc = Math.cos( c );
381 y = ( mode == SPHERICAL_NORTH ) ? Math.asin( cosc ) : Math.asin( -cosc ); // (20-14)
382 // with
383 // phi1=90
384 }
385 break;
386 }
387 case SPHERICAL_EQUATORIAL: {
388 if ( Math.abs( rho ) < TOL ) {
389 y = 0.0;
390 x = centralMeridian;
391 } else {
392 final double c = 2.0 * Math.atan( rho / k0 );
393 final double cosc = Math.cos( c );
394 final double sinc = Math.sin( c );
395 y = Math.asin( y * sinc / rho ); // (20-14) with phi1=0
396 final double t = x * sinc;
397 final double ct = rho * cosc;
398 x = ( Math.abs( t ) < TOL && Math.abs( ct ) < TOL ) ? centralMeridian
399 : Math.atan2( t, ct )
400 + centralMeridian;
401 }
402 break;
403 }
404 case SPHERICAL_OBLIQUE: {
405 if ( Math.abs( rho ) < TOL ) {
406 y = centralLatitude;
407 x = centralMeridian;
408 } else {
409 final double c = 2.0 * Math.atan( rho / k0 );
410 final double cosc = Math.cos( c );
411 final double sinc = Math.sin( c );
412 final double ct = rho * cosphi0 * cosc - y * sinphi0 * sinc; // (20-15)
413 final double t = x * sinc;
414 y = Math.asin( cosc * sinphi0 + y * sinc * cosphi0 / rho );
415 x = ( Math.abs( ct ) < TOL && Math.abs( t ) < TOL ) ? centralMeridian
416 : Math.atan2( t, ct )
417 + centralMeridian;
418 }
419 break;
420 }
421 case ELLIPSOIDAL_SOUTH: {
422 y = -y;
423 // fallthrough
424 }
425 case ELLIPSOIDAL_NORTH: {
426 final double t = rho / k0;
427 /*
428 * Compute lat using iterative technique.
429 */
430 final double halfe = e / 2.0;
431 double phi0 = 0;
432 for ( int i = MAX_ITER; --i >= 0; ) {
433 final double esinphi = e * Math.sin( phi0 );
434 final double phi = ( Math.PI / 2 )
435 - 2.0
436 * Math.atan( t
437 * Math.pow( ( 1 - esinphi ) / ( 1 + esinphi ),
438 halfe ) );
439 if ( Math.abs( phi - phi0 ) < TOL ) {
440 x = ( Math.abs( rho ) < TOL ) ? centralMeridian : Math.atan2( x, -y )
441 + centralMeridian;
442 y = ( mode == ELLIPSOIDAL_NORTH ) ? phi : -phi;
443 break choice;
444 }
445 phi0 = phi;
446 }
447 throw new TransformException( Resources.format( ResourceKeys.ERROR_NO_CONVERGENCE ) );
448 }
449 case ELLIPSOIDAL_OBLIQUE: {
450 // fallthrough
451 }
452 case ELLIPSOIDAL_EQUATORIAL: {
453 final double ce = 2.0 * Math.atan2( rho * cosChi1, k0 );
454 final double cosce = Math.cos( ce );
455 final double since = Math.sin( ce );
456 final double chi = ( Math.abs( rho ) >= TOL ) ? Math.asin( cosce
457 * sinChi1
458 + ( y * since * cosChi1 / rho ) )
459 : chi1;
460 final double t = Math.tan( Math.PI / 4.0 + chi / 2.0 );
461 /*
462 * Compute lat using iterative technique.
463 */
464 final double halfe = e / 2.0;
465 double phi0 = chi;
466 for ( int i = MAX_ITER; --i >= 0; ) {
467 final double esinphi = e * Math.sin( phi0 );
468 final double phi = 2.0
469 * Math.atan( t
470 * Math.pow( ( 1 + esinphi ) / ( 1 - esinphi ),
471 halfe ) ) - ( Math.PI / 2 );
472 if ( Math.abs( phi - phi0 ) < TOL ) {
473 x = ( Math.abs( rho ) < TOL ) ? centralMeridian
474 : Math.atan2( x * since, rho * cosChi1 * cosce - y
475 * sinChi1 * since )
476 + centralMeridian;
477 y = phi;
478 break choice;
479 }
480 phi0 = phi;
481 }
482 throw new TransformException( Resources.format( ResourceKeys.ERROR_NO_CONVERGENCE ) );
483 }
484 }
485 if ( ptDst != null ) {
486 ptDst.setLocation( x, y );
487 return ptDst;
488 }
489 return new Point2D.Double( x, y );
490 }
491
492 /**
493 * Compute part of function (3-1) from Snyder
494 */
495 private double ssfn( double phi, double sinphi ) {
496 sinphi *= e;
497 return Math.tan( ( Math.PI / 4.0 ) + phi / 2.0 )
498 * Math.pow( ( 1 - sinphi ) / ( 1 + sinphi ), e / 2.0 );
499 }
500
501 /**
502 * Returns a hash value for this map projection.
503 */
504 public int hashCode() {
505 final long code = Double.doubleToLongBits( k0 );
506 return ( (int) code ^ (int) ( code >>> 32 ) ) + 37 * super.hashCode();
507 }
508
509 /**
510 * Compares the specified object with this map projection for equality.
511 */
512 public boolean equals( final Object object ) {
513 if ( object == this )
514 return true; // Slight optimization
515 if ( super.equals( object ) ) {
516 final StereographicProjection that = (StereographicProjection) object;
517 return Double.doubleToLongBits( this.k0 ) == Double.doubleToLongBits( that.k0 )
518 && Double.doubleToLongBits( this.ak0 ) == Double.doubleToLongBits( that.ak0 )
519 && Double.doubleToLongBits( this.sinphi0 ) == Double.doubleToLongBits( that.sinphi0 )
520 && Double.doubleToLongBits( this.cosphi0 ) == Double.doubleToLongBits( that.cosphi0 )
521 && Double.doubleToLongBits( this.chi1 ) == Double.doubleToLongBits( that.chi1 )
522 && Double.doubleToLongBits( this.sinChi1 ) == Double.doubleToLongBits( that.sinChi1 )
523 && Double.doubleToLongBits( this.cosChi1 ) == Double.doubleToLongBits( that.cosChi1 );
524 }
525 return false;
526 }
527
528 /**
529 * Impl�mentation de la partie entre crochets de la cha�ne retourn�e par
530 * {@link #toString()}.
531 */
532 void toString( final StringBuffer buffer ) {
533 super.toString( buffer );
534 addParameter( buffer, "latitude_true_scale", Math.toDegrees( latitudeTrueScale ) );
535 }
536
537 /**
538 * Informations about a {@link StereographicProjection}.
539 *
540 * @version 1.0
541 * @author Martin Desruisseaux
542 */
543 static final class Provider extends MapProjection.Provider {
544 /**
545 * <code>true</code> for polar stereographic, or <code>false</code> for
546 * equatorial and oblique stereographic.
547 */
548 private final boolean polar;
549
550 /**
551 * <code>true</code> if polar/oblique/equatorial stereographic can be
552 * automatically choosen.
553 */
554 private final boolean auto;
555
556 /**
557 * Construct a new provider. The type (polar, oblique or equatorial) will be
558 * choosen automatically according the latitude or origin.
559 */
560 public Provider() {
561 super( "Stereographic", ResourceKeys.STEREOGRAPHIC_PROJECTION );
562 put( "latitude_true_scale", 90.0, LATITUDE_RANGE );
563 polar = true;
564 auto = true;
565 }
566
567 /**
568 * Construct an object for polar or oblique stereographic.
569 *
570 * @param polar <code>true</code> for polar stereographic, or <code>false</code>
571 * for equatorial and oblique stereographic.
572 */
573 public Provider( final boolean polar ) {
574 super( polar ? "Polar_Stereographic" : "Oblique_Stereographic",
575 ResourceKeys.STEREOGRAPHIC_PROJECTION );
576 put( "latitude_true_scale", polar ? 90.0 : 0.0, LATITUDE_RANGE );
577 this.polar = polar;
578 this.auto = false;
579 }
580
581 /**
582 * Create a new map projection.
583 */
584 protected Object create( final Projection parameters ) {
585 return new StereographicProjection( parameters, polar, auto );
586 }
587 }
588 }