001    //$HeadURL: svn+ssh://jwilden@svn.wald.intevation.org/deegree/base/branches/2.5_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    }