001 //$HeadURL: svn+ssh://rbezema@svn.wald.intevation.org/deegree/base/tags/2.1/src/org/deegree/model/csct/ct/TransverseMercatorProjection.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 This library is free software; you can redistribute it and/or
012 modify it under the terms of the GNU Lesser General Public
013 License as published by the Free Software Foundation; either
014 version 2.1 of the License, or (at your option) any later version.
015
016 This library is distributed in the hope that it will be useful,
017 but WITHOUT ANY WARRANTY; without even the implied warranty of
018 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
019 Lesser General Public License for more details.
020
021 You should have received a copy of the GNU Lesser General Public
022 License along with this library; if not, write to the Free Software
023 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
024
025 Contact:
026
027 Andreas Poth
028 lat/lon GmbH
029 Aennchenstr. 19
030 53115 Bonn
031 Germany
032 E-Mail: poth@lat-lon.de
033
034 Klaus Greve
035 Department of Geography
036 University of Bonn
037 Meckenheimer Allee 166
038 53115 Bonn
039 Germany
040 E-Mail: klaus.greve@uni-bonn.de
041
042
043 ---------------------------------------------------------------------------*/
044 package org.deegree.model.csct.ct;
045
046 import java.awt.geom.Point2D;
047 import java.util.ArrayList;
048
049 import org.deegree.model.csct.cs.Projection;
050 import org.deegree.model.csct.pt.Latitude;
051 import org.deegree.model.csct.resources.css.ResourceKeys;
052 import org.deegree.model.csct.resources.css.Resources;
053
054 /**
055 * Projections de Mercator tranverses Universelle et Modifi�e. Il s'agit de la projection
056 * Mercator cylindrique, mais dans lequel le cylindre a subit une rotation de 90�. Au lieu
057 * d'�tre tangeant � l'�quateur (ou � une autre latitude standard), le cylindre de la
058 * projection tranverse est tangeant � un m�ridien central. Les d�formation deviennent
059 * de plus en plus importantes � mesure que l'on s'�loigne du m�ridien central. Cette
060 * projection est appropri�e pour les r�gions qui s'�tendent d'avantage dans le sens
061 * nord-sud que dans le sens est-ouest.
062 *
063 * R�f�rence: John P. Snyder (Map Projections - A Working Manual,
064 * U.S. Geological Survey Professional Paper 1395, 1987)
065 *
066 * @version 1.0
067 * @author Andr� Gosselin
068 * @author Martin Desruisseaux
069 */
070 final class TransverseMercatorProjection extends CylindricalProjection {
071 /*
072 * Constants used to calculate {@link #en0}, {@link #en1},
073 * {@link #en2}, {@link #en3}, {@link #en4}.
074 */
075 private static final double C00 = 1.0, C02 = 0.25, C04 = 0.046875, C06 = 0.01953125,
076 C08 = 0.01068115234375, C22 = 0.75, C44 = 0.46875,
077 C46 = 0.01302083333333333333, C48 = 0.00712076822916666666,
078 C66 = 0.36458333333333333333, C68 = 0.00569661458333333333,
079 C88 = 0.3076171875;
080
081 /*
082 * Contants used for the forward and inverse transform for the eliptical
083 * case of the Transverse Mercator.
084 */
085 private static final double FC1 = 1.00000000000000000000000, // 1/1
086 FC2 = 0.50000000000000000000000, // 1/2
087 FC3 = 0.16666666666666666666666, // 1/6
088 FC4 = 0.08333333333333333333333, // 1/12
089 FC5 = 0.05000000000000000000000, // 1/20
090 FC6 = 0.03333333333333333333333, // 1/30
091 FC7 = 0.02380952380952380952380, // 1/42
092 FC8 = 0.01785714285714285714285; // 1/56
093
094 /**
095 * Relative precisions.
096 */
097 private static final double EPS10 = 1e-10;
098
099 /**
100 * Relative precisions.
101 */
102 private static final double EPS11 = 1e-11;
103
104 /**
105 * scale factor for semi mayor axis
106 */
107 private double scale_factor = 1.0;
108
109 /**
110 * Global scale factor. Value <code>ak0</code>
111 * is equals to <code>{@link #a}*k0</code>.
112 */
113 private final double ak0;
114
115 /**
116 * Constant needed for the <code>mlfn<code> method.
117 * Setup at construction time.
118 */
119 private final double en0, en1, en2, en3, en4;
120
121 /**
122 * Variante de l'eccentricit�, calcul�e par
123 * <code>e'� = (a�-b�)/b� = es/(1-es)</code>.
124 */
125 private final double esp;
126
127 /**
128 * indicates if the projection should be performed for the north hemisphere
129 * (1) or the south hemisphere (-1)
130 */
131 private int hemisphere = 1;
132
133 /**
134 * Construct a new map projection from the suplied parameters.
135 * Projection will default to Universal Transverse Mercator (UTM).
136 *
137 * @param parameters The parameter values in standard units.
138 * @throws MissingParameterException if a mandatory parameter is missing.
139 */
140 protected TransverseMercatorProjection( final Projection parameters )
141 throws MissingParameterException {
142 this( parameters, false );
143 } // Default to UTM.
144
145 /**
146 * Construct a new map projection from the suplied parameters.
147 *
148 * @param parameters The parameter values in standard units.
149 * @param modified <code>true</code> for MTM, <code>false</code> for UTM.
150 * @throws MissingParameterException if a mandatory parameter is missing.
151 */
152 protected TransverseMercatorProjection( final Projection parameters, final boolean modified )
153 throws MissingParameterException {
154
155 super( parameters );
156
157 String[] param = parameters.getParameters().getParameterListDescriptor().getParamNames();
158 ArrayList params = new ArrayList();
159
160 for ( int i = 0; i < param.length; i++ ) {
161 params.add( param[i] );
162 }
163
164 hemisphere = (int) parameters.getValue( "hemisphere", 1 );
165
166 scale_factor = 1.0;
167
168 if ( params.contains( "scale_factor" ) ) {
169 scale_factor = parameters.getParameters().getDoubleParameter( "scale_factor" );
170 }
171
172 ak0 = a * scale_factor;
173
174 double t;
175 esp = ( ( a * a ) / ( b * b ) ) - 1.0;
176 en0 = C00 - ( es * ( C02 + es * ( C04 + es * ( C06 + es * C08 ) ) ) );
177 en1 = es * ( C22 - es * ( C04 + es * ( C06 + es * C08 ) ) );
178 en2 = ( t = es * es ) * ( C44 - es * ( C46 + es * C48 ) );
179 en3 = ( t *= es ) * ( C66 - es * C68 );
180 en4 = t * es * C88;
181 }
182
183 /**
184 * Returns a human readable name localized for the specified locale.
185 */
186 public String getName() {
187 return "TransverseMercatorProjection";
188 } //Resources.getResources(locale).getString(key);
189
190 /**
191 * Calcule la distance meridionale sur un
192 * ellipso�de � la latitude <code>phi</code>.
193 */
194 private final double mlfn( final double phi, double sphi, double cphi ) {
195 cphi *= sphi;
196 sphi *= sphi;
197 return ( en0 * phi ) - ( cphi * ( en1 + sphi * ( en2 + sphi * ( en3 + sphi * en4 ) ) ) );
198 }
199
200 /**
201 * Transforms the specified (<var>x</var>,<var>y</var>) coordinate
202 * and stores the result in <code>ptDst</code>.
203 */
204 protected Point2D transform( double x, double y, final Point2D ptDst )
205 throws TransformException {
206
207 if ( Math.abs( y ) > ( Math.PI / 2.0 - EPS ) ) {
208 throw new TransformException( Resources.format( ResourceKeys.ERROR_POLE_PROJECTION_$1,
209 new Latitude( Math.toDegrees( y ) ) ) );
210 }
211
212 x -= centralMeridian;
213 y -= centralLatitude;
214 y *= hemisphere;
215
216 double sinphi = Math.sin( y );
217 double cosphi = Math.cos( y );
218
219 if ( isSpherical ) {
220 // Spherical model.
221 double b = cosphi * Math.sin( x );
222
223 if ( Math.abs( Math.abs( b ) - 1.0 ) <= EPS10 ) {
224 throw new TransformException(
225 Resources.format( ResourceKeys.ERROR_VALUE_TEND_TOWARD_INFINITY ) );
226 }
227
228 double yy = ( cosphi * Math.cos( x ) ) / Math.sqrt( 1.0 - ( b * b ) );
229 x = ( 0.5 * ak0 * Math.log( ( 1.0 + b ) / ( 1.0 - b ) ) ) + false_easting;/* 8-1 */
230
231 if ( ( b = Math.abs( yy ) ) >= 1.0 ) {
232 if ( ( b - 1.0 ) > EPS10 ) {
233 throw new TransformException(
234 Resources.format( ResourceKeys.ERROR_VALUE_TEND_TOWARD_INFINITY ) );
235 }
236 yy = 0.0;
237
238 } else {
239 yy = Math.acos( yy );
240 }
241
242 if ( y < 0 ) {
243 yy = -yy;
244 }
245
246 y = ak0 * yy;
247 } else {
248
249 // Ellipsoidal model.
250 double t = ( Math.abs( cosphi ) > EPS10 ) ? ( sinphi / cosphi ) : 0;
251 t *= t;
252
253 double al = cosphi * x;
254 double als = al * al;
255 al /= Math.sqrt( 1.0 - ( es * sinphi * sinphi ) );
256
257 double n = esp * cosphi * cosphi;
258
259 /* NOTE: meridinal distance at central latitude is always 0 */
260 y = ( ak0 * ( mlfn( y, sinphi, cosphi ) + sinphi
261 * al
262 * x
263 * FC2
264 * ( 1.0 + FC4
265 * als
266 * ( 5.0 - t
267 + ( n * ( 9.0 + 4.0 * n ) ) + FC6
268 * als
269 * ( 61.0
270 + ( t * ( t - 58.0 ) )
271 + ( n * ( 270.0 - 330.0 * t ) ) + FC8
272 * als
273 * ( 1385.0 + t
274 * ( t
275 * ( 543.0 - t ) - 3111.0 ) ) ) ) ) ) );
276 y += false_northing;
277
278 x = ( ak0 * al * ( FC1 + FC3
279 * als
280 * ( 1.0 - t + n + FC5
281 * als
282 * ( 5.0 + ( t * ( t - 18.0 ) )
283 + ( n * ( 14.0 - 58.0 * t ) ) + FC7
284 * als
285 * ( 61.0 + t
286 * ( t
287 * ( 179.0 - t ) - 479.0 ) ) ) ) ) );
288 x += false_easting;
289 }
290
291 if ( ptDst != null ) {
292 ptDst.setLocation( x, y );
293 return ptDst;
294 }
295 return new Point2D.Double( x, y );
296
297 }
298
299 /**
300 * Transforms the specified (<var>x</var>,<var>y</var>) coordinate
301 * and stores the result in <code>ptDst</code>.
302 */
303 protected Point2D inverseTransform( double x, double y, final Point2D ptDst )
304 throws TransformException {
305 x -= false_easting;
306 y -= false_northing;
307 y *= hemisphere;
308 //x = x - ( (int)( x / 1000000 ) * 1000000.0 );
309
310 if ( isSpherical ) {
311 // Spherical model.
312 double t = Math.exp( x / ak0 );
313 double d = 0.5 * ( t - 1 / t );
314 t = Math.cos( y / ak0 );
315
316 double phi = Math.asin( Math.sqrt( ( 1.0 - t * t ) / ( 1.0 + d * d ) ) );
317 y = ( y < 0.0 ) ? ( -phi ) : phi;
318 x = ( Math.abs( d ) > EPS10 || Math.abs( t ) > EPS10 ) ? ( Math.atan2( d, t ) + centralMeridian )
319 : centralMeridian;
320 } else {
321 // Ellipsoidal projection.
322 final double y_ak0 = y / ak0;
323 final double k = 1.0 - es;
324 double phi = y_ak0;
325
326 for ( int i = 20; true; ) // rarely goes over 5 iterations
327 {
328 if ( ( --i ) < 0 ) {
329 throw new TransformException(
330 Resources.format( ResourceKeys.ERROR_NO_CONVERGENCE ) );
331 }
332
333 final double s = Math.sin( phi );
334 double t = 1.0 - ( es * ( s * s ) );
335 t = ( mlfn( phi, s, Math.cos( phi ) ) - y_ak0 ) / ( k * t * Math.sqrt( t ) );
336 phi -= t;
337
338 if ( Math.abs( t ) < EPS11 ) {
339 break;
340 }
341 }
342
343 if ( Math.abs( phi ) >= ( Math.PI / 2.0 ) ) {
344 y = ( y < 0.0 ) ? ( -( Math.PI / 2.0 ) ) : ( Math.PI / 2.0 );
345 x = centralMeridian;
346 } else {
347 double sinphi = Math.sin( phi );
348 double cosphi = Math.cos( phi );
349 double t = ( Math.abs( cosphi ) > EPS10 ) ? ( sinphi / cosphi ) : 0.0;
350 double n = esp * cosphi * cosphi;
351 double con = 1.0 - ( es * sinphi * sinphi );
352 double d = ( x * Math.sqrt( con ) ) / ak0;
353 con *= t;
354 t *= t;
355
356 double ds = d * d;
357 y = phi
358 - ( ( con * ds / ( 1.0 - es ) ) * FC2 * ( 1.0 - ds
359 * FC4
360 * ( 5.0
361 + ( t * ( 3.0 - 9.0 * n ) )
362 + ( n * ( 1.0 - 4 * n ) ) - ds
363 * FC6
364 * ( 61.0
365 + ( t * ( 90.0 - ( 252.0 * n ) + 45.0 * t ) )
366 + ( 46.0 * n ) - ds
367 * FC8
368 * ( 1385.0 + t
369 * ( 3633.0 + t
370 * ( 4095.0 + 1574.0 * t ) ) ) ) ) ) )
371 + centralLatitude;
372
373 x = ( ( d * ( FC1 - ds
374 * FC3
375 * ( 1.0 + ( 2.0 * t ) + n - ds
376 * FC5
377 * ( 5.0
378 + ( t * ( 28.0 + ( 24 * t ) + 8.0 * n ) )
379 + ( 6.0 * n ) - ds
380 * FC7
381 * ( 61.0 + t
382 * ( 662.0 + t
383 * ( 1320.0 + 720.0 * t ) ) ) ) ) ) ) / cosphi )
384 + centralMeridian;
385 }
386 }
387
388 if ( ptDst != null ) {
389 ptDst.setLocation( x, y );
390 return ptDst;
391 }
392 return new Point2D.Double( x, y );
393
394 }
395
396 /**
397 * Returns a hash value for this projection.
398 */
399 public int hashCode() {
400 final long code = Double.doubleToLongBits( false_easting );
401 return ( (int) code ^ (int) ( code >>> 32 ) ) + ( 37 * super.hashCode() );
402 }
403
404 /**
405 * Compares the specified object with
406 * this map projection for equality.
407 */
408 public boolean equals( final Object object ) {
409 if ( object == this ) {
410 return true; // Slight optimization
411 }
412
413 if ( super.equals( object ) ) {
414 final TransverseMercatorProjection that = (TransverseMercatorProjection) object;
415 return ( Double.doubleToLongBits( this.false_easting ) == Double.doubleToLongBits( that.false_easting ) )
416 && ( Double.doubleToLongBits( this.ak0 ) == Double.doubleToLongBits( that.ak0 ) );
417 }
418
419 return false;
420 }
421
422 /**
423 * Impl�mentation de la partie entre crochets
424 * de la cha�ne retourn�e par {@link #toString()}.
425 */
426 void toString( final StringBuffer buffer ) {
427 super.toString( buffer );
428 }
429
430 /**
431 * Informations about a {@link TransverseMercatorProjection}.
432 *
433 * @version 1.0
434 * @author Martin Desruisseaux
435 */
436 static final class Provider extends MapProjection.Provider {
437 /**
438 * <code>true</code> for Modified Mercator Projection (MTM), or
439 * <code>false</code> for Universal Mercator Projection (UTM).
440 */
441 private final boolean modified;
442
443 /**
444 * Construct a new registration.
445 *
446 * @param modified <code>true</code> for Modified Mercator Projection (MTM),
447 * or <code>false</code> for Universal Mercator Projection (UTM).
448 */
449 public Provider( final boolean modified ) {
450 super( modified ? "Modified_Transverse_Mercator" : "Transverse_Mercator",
451 modified ? ResourceKeys.MTM_PROJECTION : ResourceKeys.UTM_PROJECTION );
452 this.modified = modified;
453 }
454
455 /**
456 * Create a new map projection.
457 */
458 protected Object create( final Projection parameters ) {
459 return new TransverseMercatorProjection( parameters, modified );
460 }
461 }
462 }