001// License: GPL. For details, see LICENSE file.
002package org.openstreetmap.josm.data.projection.proj;
003
004import static org.openstreetmap.josm.tools.I18n.tr;
005
006import org.openstreetmap.josm.data.Bounds;
007import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
008
009/**
010 * Transverse Mercator Projection (EPSG code 9807). This
011 * is a cylindrical projection, in which the cylinder has been rotated 90°.
012 * Instead of being tangent to the equator (or to an other standard latitude),
013 * it is tangent to a central meridian. Deformation are more important as we
014 * are going futher from the central meridian. The Transverse Mercator
015 * projection is appropriate for region wich have a greater extent north-south
016 * than east-west.
017 * <p>
018 *
019 * The elliptical equations used here are series approximations, and their accuracy
020 * decreases as points move farther from the central meridian of the projection.
021 * The forward equations here are accurate to a less than a mm &plusmn;10 degrees from
022 * the central meridian, a few mm &plusmn;15 degrees from the
023 * central meridian and a few cm &plusmn;20 degrees from the central meridian.
024 * The spherical equations are not approximations and should always give the
025 * correct values.
026 * <p>
027 *
028 * There are a number of versions of the transverse mercator projection
029 * including the Universal (UTM) and Modified (MTM) Transverses Mercator
030 * projections. In these cases the earth is divided into zones. For the UTM
031 * the zones are 6 degrees wide, numbered from 1 to 60 proceeding east from
032 * 180 degrees longitude, and between lats 84 degrees North and 80
033 * degrees South. The central meridian is taken as the center of the zone
034 * and the latitude of origin is the equator. A scale factor of 0.9996 and
035 * false easting of 500000m is used for all zones and a false northing of 10000000m
036 * is used for zones in the southern hemisphere.
037 * <p>
038 *
039 * NOTE: formulas used below are not those of Snyder, but rather those
040 *       from the {@code proj4} package of the USGS survey, which
041 *       have been reproduced verbatim. USGS work is acknowledged here.
042 * <p>
043 *
044 * This class has been derived from the implementation of the Geotools project;
045 * git 8cbf52d, org.geotools.referencing.operation.projection.TransverseMercator
046 * at the time of migration.
047 * <p>
048 *
049 * <b>References:</b>
050 * <ul>
051 *   <li> Proj-4.4.6 available at <A HREF="http://www.remotesensing.org/proj">www.remotesensing.org/proj</A><br>
052 *        Relevent files are: {@code PJ_tmerc.c}, {@code pj_mlfn.c}, {@code pj_fwd.c} and {@code pj_inv.c}.</li>
053 *   <li> John P. Snyder (Map Projections - A Working Manual,
054 *        U.S. Geological Survey Professional Paper 1395, 1987).</li>
055 *   <li> "Coordinate Conversions and Transformations including Formulas",
056 *        EPSG Guidence Note Number 7, Version 19.</li>
057 * </ul>
058 *
059 * @author André Gosselin
060 * @author Martin Desruisseaux (PMO, IRD)
061 * @author Rueben Schulz
062 *
063 * @see <A HREF="http://mathworld.wolfram.com/MercatorProjection.html">Transverse Mercator projection on MathWorld</A>
064 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/transverse_mercator.html">"Transverse_Mercator" on RemoteSensing.org</A>
065 */
066public class TransverseMercator extends AbstractProj {
067
068    /**
069     * Contants used for the forward and inverse transform for the eliptical
070     * case of the Transverse Mercator.
071     */
072    private static final double FC1 = 1.00000000000000000000000,  // 1/1
073                                FC2 = 0.50000000000000000000000,  // 1/2
074                                FC3 = 0.16666666666666666666666,  // 1/6
075                                FC4 = 0.08333333333333333333333,  // 1/12
076                                FC5 = 0.05000000000000000000000,  // 1/20
077                                FC6 = 0.03333333333333333333333,  // 1/30
078                                FC7 = 0.02380952380952380952380,  // 1/42
079                                FC8 = 0.01785714285714285714285;  // 1/56
080
081    /**
082     * Maximum difference allowed when comparing real numbers.
083     */
084    private static final double EPSILON = 1E-6;
085
086    /**
087     * A derived quantity of excentricity, computed by <code>e'² = (a²-b²)/b² = es/(1-es)</code>
088     * where <var>a</var> is the semi-major axis length and <var>b</var> is the semi-minor axis
089     * length.
090     */
091    private double eb2;
092
093    /**
094     * Latitude of origin in <u>radians</u>. Default value is 0, the equator.
095     * This is called '<var>phi0</var>' in Snyder.
096     * <p>
097     * <strong>Consider this field as final</strong>. It is not final only
098     * because some classes need to modify it at construction time.
099     */
100    protected double latitudeOfOrigin;
101
102    /**
103     * Meridian distance at the {@code latitudeOfOrigin}.
104     * Used for calculations for the ellipsoid.
105     */
106    private double ml0;
107
108    @Override
109    public String getName() {
110        return tr("Transverse Mercator");
111    }
112
113    @Override
114    public String getProj4Id() {
115        return "tmerc";
116    }
117
118    @Override
119    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
120        super.initialize(params);
121        eb2 = params.ellps.eb2;
122        latitudeOfOrigin = params.lat0 == null ? 0 : Math.toRadians(params.lat0);
123        ml0 = mlfn(latitudeOfOrigin, Math.sin(latitudeOfOrigin), Math.cos(latitudeOfOrigin));
124    }
125
126    @Override
127    public double[] project(double y, double x) {
128        x = normalizeLon(x);
129        double sinphi = Math.sin(y);
130        double cosphi = Math.cos(y);
131
132        double t = (Math.abs(cosphi) > EPSILON) ? sinphi/cosphi : 0;
133        t *= t;
134        double al = cosphi*x;
135        double als = al*al;
136        al /= Math.sqrt(1.0 - e2 * sinphi*sinphi);
137        double n = eb2 * cosphi*cosphi;
138
139        /* NOTE: meridinal distance at latitudeOfOrigin is always 0 */
140        y = (mlfn(y, sinphi, cosphi) - ml0 +
141            sinphi * al * x *
142            FC2 * (1.0 +
143            FC4 * als * (5.0 - t + n*(9.0 + 4.0*n) +
144            FC6 * als * (61.0 + t * (t - 58.0) + n*(270.0 - 330.0*t) +
145            FC8 * als * (1385.0 + t * (t*(543.0 - t) - 3111.0))))));
146
147        x = al*(FC1 + FC3 * als*(1.0 - t + n +
148            FC5 * als * (5.0 + t*(t - 18.0) + n*(14.0 - 58.0*t) +
149            FC7 * als * (61.0+ t*(t*(179.0 - t) - 479.0)))));
150
151        return new double[] {x, y};
152    }
153
154    @Override
155    public double[] invproject(double x, double y) {
156        double phi = inv_mlfn(ml0 + y);
157
158        if (Math.abs(phi) >= Math.PI/2) {
159            y = y < 0.0 ? -(Math.PI/2) : (Math.PI/2);
160            x = 0.0;
161        } else {
162            double sinphi = Math.sin(phi);
163            double cosphi = Math.cos(phi);
164            double t = (Math.abs(cosphi) > EPSILON) ? sinphi/cosphi : 0.0;
165            double n = eb2 * cosphi*cosphi;
166            double con = 1.0 - e2 * sinphi*sinphi;
167            double d = x * Math.sqrt(con);
168            con *= t;
169            t *= t;
170            double ds = d*d;
171
172            y = phi - (con*ds / (1.0 - e2)) *
173                FC2 * (1.0 - ds *
174                FC4 * (5.0 + t*(3.0 - 9.0*n) + n*(1.0 - 4*n) - ds *
175                FC6 * (61.0 + t*(90.0 - 252.0*n + 45.0*t) + 46.0*n - ds *
176                FC8 * (1385.0 + t*(3633.0 + t*(4095.0 + 1574.0*t))))));
177
178            x = d*(FC1 - ds * FC3 * (1.0 + 2.0*t + n -
179                ds*FC5*(5.0 + t*(28.0 + 24* t + 8.0*n) + 6.0*n -
180                ds*FC7*(61.0 + t*(662.0 + t*(1320.0 + 720.0*t))))))/cosphi;
181        }
182        return new double[] {y, x};
183    }
184
185    @Override
186    public Bounds getAlgorithmBounds() {
187        return new Bounds(-89, -7, 89, 7, false);
188    }
189}