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.cylindric;
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.CRSException;
046 import org.deegree.crs.projections.ProjectionUtils;
047 import org.deegree.framework.log.ILogger;
048 import org.deegree.framework.log.LoggerFactory;
049
050 /**
051 * The <code>TransverseMercator</code> projection has following properties:
052 * <ul>
053 * <li>Cylindrical (transverse)</li>
054 * <li>Conformal</li>
055 * <li>The central meridian, each meridian 90° from central meridian and the equator are straight lines</li>
056 * <li>All other meridians and parallels are complex curves</li>
057 * <li>Scale is true along central meridian or along two straight lines equidistant from and parallel to central
058 * merdian. (These lines are only approximately straight for the ellipsoid)</li>
059 * <li>Scale becomes infinite on sphere 90° from central meridian</li>
060 * <li>Used extensively for quadrangle maps at scales from 1:24.000 to 1:250.000</li>
061 * <li>Often used to show regions with greater north-south extent</li>
062 * <li>presented by lambert in 1772</li>
063 * </ul>
064 *
065 * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
066 *
067 * @author last edited by: $Author:$
068 *
069 * @version $Revision:$, $Date:$
070 *
071 */
072
073 public class TransverseMercator extends CylindricalProjection {
074
075 private static ILogger LOG = LoggerFactory.getLogger( TransverseMercator.class );
076 /**
077 * Contants used for the forward and inverse transform for the eliptical case of the Transverse Mercator.
078 */
079 private static final double FC1 = 1., // 1/1
080 FC2 = 0.5, // 1/2
081 FC3 = 0.16666666666666666666666, // 1/6
082 FC4 = 0.08333333333333333333333, // 1/12
083 FC5 = 0.05, // 1/20
084 FC6 = 0.03333333333333333333333, // 1/30
085 FC7 = 0.02380952380952380952380, // 1/42
086 FC8 = 0.01785714285714285714285; // 1/56
087
088 // 1 for northern hemisphere -1 for southern.
089 private int hemisphere;
090
091 // esp will can hold two values, for the sphere it will hold the scale, for the ellipsoid Snyder (p.61 8-12).
092 private double esp;
093
094 private double ml0;
095
096 private double[] en;
097
098 /**
099 * @param northernHemisphere
100 * true if on the northernhemisphere false otherwhise.
101 * @param geographicCRS
102 * @param falseNorthing
103 * @param falseEasting
104 * @param naturalOrigin
105 * @param units
106 * @param scale
107 * @param identifiers
108 * @param names
109 * @param versions
110 * @param descriptions
111 * @param areasOfUse
112 */
113 public TransverseMercator( boolean northernHemisphere, GeographicCRS geographicCRS, double falseNorthing,
114 double falseEasting, Point2d naturalOrigin, Unit units, double scale,
115 String[] identifiers, String[] names, String[] versions, String[] descriptions,
116 String[] areasOfUse ) {
117 super( geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, scale, true,// always conformal
118 false, // not equalArea
119 identifiers,
120 names,
121 versions,
122 descriptions,
123 areasOfUse );
124 this.hemisphere = ( northernHemisphere ) ? 1 : -1;
125 if ( isSpherical() ) {
126 esp = getScale();
127 ml0 = .5 * esp;
128 } else {
129 en = ProjectionUtils.getRectifiyingLatitudeValues( getSquaredEccentricity() );
130 ml0 = ProjectionUtils.getDistanceAlongMeridian( getProjectionLatitude(), getSinphi0(), getCosphi0(), en );
131 esp = getSquaredEccentricity() / ( 1. - getSquaredEccentricity() );
132 // esp = ( ( getSemiMajorAxis() * getSemiMajorAxis() ) / ( getSemiMinorAxis() * getSemiMinorAxis() ) ) -
133 // 1.0;
134 }
135
136 }
137
138 /**
139 * @param northernHemisphere
140 * true if on the northernhemisphere false otherwhise.
141 * @param geographicCRS
142 * @param falseNorthing
143 * @param falseEasting
144 * @param naturalOrigin
145 * @param units
146 * @param scale
147 * @param identifier
148 * @param name
149 * @param version
150 * @param description
151 * @param areaOfUse
152 */
153 public TransverseMercator( boolean northernHemisphere, GeographicCRS geographicCRS, double falseNorthing,
154 double falseEasting, Point2d naturalOrigin, Unit units, double scale, String identifier,
155 String name, String version, String description, String areaOfUse ) {
156 this( northernHemisphere,
157 geographicCRS,
158 falseNorthing,
159 falseEasting,
160 naturalOrigin,
161 units,
162 scale,
163 new String[] { identifier },
164 new String[] { name },
165 new String[] { version },
166 new String[] { description },
167 new String[] { areaOfUse } );
168 }
169
170 /**
171 * @param northernHemisphere
172 * true if on the northernhemisphere false otherwhise.
173 * @param geographicCRS
174 * @param falseNorthing
175 * @param falseEasting
176 * @param naturalOrigin
177 * @param units
178 * @param scale
179 * @param identifiers
180 */
181 public TransverseMercator( boolean northernHemisphere, GeographicCRS geographicCRS, double falseNorthing,
182 double falseEasting, Point2d naturalOrigin, Unit units, double scale,
183 String[] identifiers ) {
184 this( northernHemisphere,
185 geographicCRS,
186 falseNorthing,
187 falseEasting,
188 naturalOrigin,
189 units,
190 scale,
191 identifiers,
192 null,
193 null,
194 null,
195 null );
196 }
197
198 /**
199 * @param northernHemisphere
200 * true if on the northernhemisphere false otherwhise.
201 * @param geographicCRS
202 * @param falseNorthing
203 * @param falseEasting
204 * @param naturalOrigin
205 * @param units
206 * @param scale
207 * @param identifier
208 */
209 public TransverseMercator( boolean northernHemisphere, GeographicCRS geographicCRS, double falseNorthing,
210 double falseEasting, Point2d naturalOrigin, Unit units, double scale, String identifier ) {
211 this( northernHemisphere,
212 geographicCRS,
213 falseNorthing,
214 falseEasting,
215 naturalOrigin,
216 units,
217 scale,
218 new String[] { identifier } );
219 }
220
221 /**
222 * A northern hemisphere conformal transverse mercator projection with a scale of one. Using the given datum.
223 *
224 * @param geographicCRS
225 * @param falseNorthing
226 * @param falseEasting
227 * @param naturalOrigin
228 * @param units
229 * @param identifiers
230 * @param name
231 */
232 public TransverseMercator( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
233 Point2d naturalOrigin, Unit units, String[] identifiers ) {
234 this( true, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, 1., identifiers );
235 }
236
237 /**
238 * A northern hemisphere conformal transverse mercator projection with a scale of one. Using the given datum.
239 *
240 * @param geographicCRS
241 * @param falseNorthing
242 * @param falseEasting
243 * @param naturalOrigin
244 * @param units
245 * @param identifier
246 * @param name
247 */
248 public TransverseMercator( GeographicCRS geographicCRS, double falseNorthing, double falseEasting,
249 Point2d naturalOrigin, Unit units, String identifier ) {
250 this( true, geographicCRS, falseNorthing, falseEasting, naturalOrigin, units, 1., identifier );
251 }
252
253 @Override
254 public Point2d doInverseProjection( double x, double y )
255 throws CRSException {
256 Point2d result = new Point2d( 0, 0 );
257 LOG.logDebug( "InverseProjection, incoming points x: " + x + " y: " + y );
258 x -= getFalseEasting();
259 y -= getFalseNorthing();
260 y *= hemisphere;
261
262 if ( isSpherical() ) {
263 // h holds e^x, the sinh = 0.5*(e^x - e^-x), cosh = 0.5(e^x + e^-x)
264 double h = Math.exp( x / getScaleFactor() );
265 // sinh holds the sinh from Snyder (p.60 8-7)
266 double sinh = .5 * ( h - 1. / h );
267
268 // Snyder (p.60 8-8)
269 // reuse variable
270 double cosD = Math.cos( getProjectionLatitude() + ( y / getScaleFactor() ) );
271 /**
272 * To calc phi from Snyder (p.60 8-6), use following trick! sin^2(D) + cos^2(D) = 1 => sin(D) = sqrt( 1-
273 * cos^2(D) ) and cosh^2(x) - sin^2(x) = 1 => cosh(x) = sqrt( 1+sin^2(x) )
274 */
275 result.y = ProjectionUtils.asinScaled( Math.sqrt( ( 1. - cosD * cosD ) / ( 1. + sinh * sinh ) ) );
276 // if ( y < 0 ) {// southern hemisphere
277 // out.y = -out.y;
278 // }
279 result.x = Math.atan2( sinh, cosD );
280 } else {
281 // out.y will hold the phi_1 from Snyder (p.63 8-18).
282 result.y = ProjectionUtils.calcPhiFromMeridianDistance( ml0 + ( y / getScaleFactor() ),
283 getSquaredEccentricity(),
284 en );
285 if ( Math.abs( result.y ) >= ProjectionUtils.HALFPI ) {
286 result.y = y < 0. ? -ProjectionUtils.HALFPI : ProjectionUtils.HALFPI;
287 result.x = 0;
288 } else {
289
290 double sinphi = Math.sin( result.y );
291 double cosphi = Math.cos( result.y );
292 // largeT Will hold the tan^2(phi) Snyder (p.64 8-22).
293 double largeT = ( Math.abs( cosphi ) > ProjectionUtils.EPS10 ) ? sinphi / cosphi : 0;
294
295 // will hold the C_1 from Synder (p.64 8-21)
296 double largeC = esp * cosphi * cosphi;
297
298 // Holds a modified N from Synder (p.64 8-23), multiplied with the largeT, it is the first term fo the
299 // calculation of phi e.g. N*T/R
300 double con = 1. - ( getSquaredEccentricity() * sinphi * sinphi );
301 // largeD holds the D from Snyder (p.64 8-25). (x/(1/N) = x*N)
302 double largeD = x * Math.sqrt( con ) / getScaleFactor();
303 con *= largeT;
304 largeT *= largeT;
305 double ds = largeD * largeD;
306
307 /**
308 * As for the forward projection, I'm not sure if this is correct, this should be checked!
309 */
310 result.y -= ( con * ds / ( 1. - getSquaredEccentricity() ) ) * FC2
311 * ( 1. - ds * FC4
312 * ( 5. + largeT * ( 3. - 9. * largeC ) + largeC * ( 1. - 4 * largeC ) - ds * FC6
313 * ( 61. + largeT
314 * ( 90. - 252. * largeC + 45. * largeT )
315 + 46.
316 * largeC - ds * FC8
317 * ( 1385. + largeT * ( 3633. + largeT * ( 4095. + 1574. * largeT ) ) ) ) ) );
318 result.x = largeD * ( FC1 - ds * FC3
319 * ( 1. + 2. * largeT + largeC - ds * FC5
320 * ( 5. + largeT
321 * ( 28. + 24. * largeT + 8. * largeC )
322 + 6.
323 * largeC - ds * FC7
324 * ( 61. + largeT * ( 662. + largeT * ( 1320. + 720. * largeT ) ) ) ) ) )
325 / cosphi;
326 }
327 }
328 result.y += getProjectionLatitude();
329 result.x += getProjectionLongitude();
330
331 return result;
332 }
333
334 /*
335 * (non-Javadoc)
336 *
337 * @see org.deegree.crs.projections.Projection#doProjection(double, double)
338 */
339 @Override
340 public Point2d doProjection( double lambda, double phi )
341 throws CRSException {
342 LOG.logDebug( "Projection, incoming points lambda: " + lambda + " phi: " + phi );
343 Point2d result = new Point2d( 0, 0 );
344 lambda -= getProjectionLongitude();
345 phi -= getProjectionLatitude();
346 phi *= hemisphere;
347 double cosphi = Math.cos( phi );
348 if ( isSpherical() ) {
349 double b = cosphi * Math.sin( lambda );
350
351 // Snyder (p.58 8-1)
352 result.x = ml0 * getScaleFactor() * Math.log( ( 1. + b ) / ( 1. - b ) );
353
354 // reformed and inserted the k from (p.58 8-4), so no tangens has to be calculated.
355 double ty = cosphi * Math.cos( lambda ) / Math.sqrt( 1. - b * b );
356 ty = ProjectionUtils.acosScaled( ty );
357 if ( phi < 0.0 ) {
358 ty = -ty;
359 }
360 // esp just holds the scale
361 result.y = esp * ( ty - getProjectionLatitude() );
362 } else {
363 double sinphi = Math.sin( phi );
364 double largeT = ( Math.abs( cosphi ) > ProjectionUtils.EPS10 ) ? sinphi / cosphi : 0.0;
365 // largeT holds Snyder (p.61 8-13).
366 largeT *= largeT;
367 double largeA = cosphi * lambda;
368 double squaredLargeA = largeA * largeA;
369 // largeA now holds A/N Snyder (p.61 4-20 and 8-15)
370 largeA /= Math.sqrt( 1. - getSquaredEccentricity() * sinphi * sinphi );
371
372 // largeC will hold Snyder (p.61 8-14), esp holds Snyder (p.61 8-12).
373 double largeC = esp * cosphi * cosphi;
374 double largeM = ProjectionUtils.getDistanceAlongMeridian( phi, sinphi, cosphi, en );
375
376 result.x = getScaleFactor() * largeA
377 * ( FC1 + FC3 * squaredLargeA
378 * ( 1. - largeT + largeC + FC5 * squaredLargeA
379 * ( 5. + largeT
380 * ( largeT - 18. )
381 + largeC
382 * ( 14. - 58. * largeT ) + FC7 * squaredLargeA
383 * ( 61. + largeT * ( largeT * ( 179. - largeT ) - 479. ) ) ) ) );
384
385 result.y = getScaleFactor() * ( largeM - ml0 + sinphi * largeA
386 * lambda
387 * FC2
388 * ( 1. + FC4 * squaredLargeA
389 * ( 5. - largeT + largeC * ( 9. + 4. * largeC ) + FC6 * squaredLargeA
390 * ( 61. + largeT
391 * ( largeT - 58. )
392 + largeC
393 * ( 270. - 330 * largeT ) + FC8 * squaredLargeA
394 * ( 1385. + largeT * ( largeT * ( 543. - largeT ) - 3111. ) ) ) ) ) );
395
396 }
397 result.x += getFalseEasting();
398 result.y += getFalseNorthing();
399 return result;
400 }
401
402 /**
403 * @param latitude
404 * to get the nearest paralles to.
405 * @return the nearest parallel in radians of given latitude
406 */
407 public int getRowFromNearestParallel( double latitude ) {
408 int degrees = (int) Math.round( Math.toDegrees( ProjectionUtils.normalizeLatitude( latitude ) ) );
409 if ( degrees < -80 || degrees > 84 ) {
410 return 0;
411 }
412 if ( degrees > 80 ) {
413 return 24; // last parallel
414 }
415 return ( ( degrees + 80 ) / 8 ) + 3;
416 }
417
418 /**
419 * the utm zone from a given meridian
420 *
421 * @param longitude
422 * in radians
423 * @return the utm zone.
424 */
425 public int getZoneFromNearestMeridian( double longitude ) {
426 int zone = (int) Math.floor( ( ProjectionUtils.normalizeLongitude( longitude ) + Math.PI ) * 30.0 / Math.PI ) + 1;
427 if ( zone < 1 ) {
428 zone = 1;
429 } else if ( zone > 60 ) {
430 zone = 60;
431 }
432 return zone;
433 }
434
435 /**
436 * Set the zone for this transverse mercator, not yet used.
437 *
438 * @param zone
439 * the new utm zone.
440 */
441 public void setUTMZone( int zone ) {
442 zone--;
443 setProjectionLongitude( ( zone + .5 ) * Math.PI / 30. - Math.PI );
444
445 setProjectionLatitude( 0.0 );
446 if ( getScale() != 0.9996 ) {
447 setScale( 0.9996 );
448 }
449 if ( Math.abs( getFalseEasting() - 500000 ) > 0.000001 ) {
450 setFalseEasting( 500000 );
451 }
452 if ( !isSpherical() ) {
453 // recalculate the rectifying latitudes and the distance along the meridian.
454 en = ProjectionUtils.getRectifiyingLatitudeValues( getSquaredEccentricity() );
455 ml0 = ProjectionUtils.getDistanceAlongMeridian( getProjectionLatitude(), getSinphi0(), getCosphi0(), en );
456 }
457 }
458
459 @Override
460 public String getDeegreeSpecificName() {
461 return "transverseMercator";
462 }
463
464 /**
465 * @return the true if defined on the northern hemisphere.
466 */
467 public final boolean getHemisphere() {
468 return ( hemisphere == 1 );
469 }
470 }