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