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.coor.LatLon;
008import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
009import org.openstreetmap.josm.tools.CheckParameterUtil;
010import org.openstreetmap.josm.tools.Pair;
011import org.openstreetmap.josm.tools.Utils;
012
013/**
014 * Transverse Mercator Projection (EPSG code 9807). This
015 * is a cylindrical projection, in which the cylinder has been rotated 90°.
016 * Instead of being tangent to the equator (or to an other standard latitude),
017 * it is tangent to a central meridian. Deformation are more important as we
018 * are going further from the central meridian. The Transverse Mercator
019 * projection is appropriate for region which have a greater extent north-south
020 * than east-west.
021 * <p>
022 *
023 * The elliptical equations used here are series approximations, and their accuracy
024 * decreases as points move farther from the central meridian of the projection.
025 * The forward equations here are accurate to a less than a mm &plusmn;10 degrees from
026 * the central meridian, a few mm &plusmn;15 degrees from the
027 * central meridian and a few cm &plusmn;20 degrees from the central meridian.
028 * The spherical equations are not approximations and should always give the
029 * correct values.
030 * <p>
031 *
032 * There are a number of versions of the transverse mercator projection
033 * including the Universal (UTM) and Modified (MTM) Transverses Mercator
034 * projections. In these cases the earth is divided into zones. For the UTM
035 * the zones are 6 degrees wide, numbered from 1 to 60 proceeding east from
036 * 180 degrees longitude, and between lats 84 degrees North and 80
037 * degrees South. The central meridian is taken as the center of the zone
038 * and the latitude of origin is the equator. A scale factor of 0.9996 and
039 * false easting of 500000m is used for all zones and a false northing of 10000000m
040 * is used for zones in the southern hemisphere.
041 * <p>
042 *
043 * NOTE: formulas used below are not those of Snyder, but rather those
044 *       from the {@code proj4} package of the USGS survey, which
045 *       have been reproduced verbatim. USGS work is acknowledged here.
046 * <p>
047 *
048 * This class has been derived from the implementation of the Geotools project;
049 * git 8cbf52d, org.geotools.referencing.operation.projection.TransverseMercator
050 * at the time of migration.
051 * <p>
052 * The non-standard parameter <code>gamma</code> has been added as a method
053 * to rotate the projected coordinates by a certain angle (clockwise, see
054 * {@link ObliqueMercator}).
055 * <p>
056 * <b>References:</b>
057 * <ul>
058 *   <li> Proj-4.4.6 available at <A HREF="http://www.remotesensing.org/proj">www.remotesensing.org/proj</A><br>
059 *        Relevent files are: {@code PJ_tmerc.c}, {@code pj_mlfn.c}, {@code pj_fwd.c} and {@code pj_inv.c}.</li>
060 *   <li> John P. Snyder (Map Projections - A Working Manual,
061 *        U.S. Geological Survey Professional Paper 1395, 1987).</li>
062 *   <li> "Coordinate Conversions and Transformations including Formulas",
063 *        EPSG Guidence Note Number 7, Version 19.</li>
064 * </ul>
065 *
066 * @author André Gosselin
067 * @author Martin Desruisseaux (PMO, IRD)
068 * @author Rueben Schulz
069 *
070 * @see <A HREF="http://mathworld.wolfram.com/MercatorProjection.html">Transverse Mercator projection on MathWorld</A>
071 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/transverse_mercator.html">"Transverse_Mercator" on RemoteSensing.org</A>
072 */
073public class TransverseMercator extends AbstractProj {
074
075    /** Earth emispheres **/
076    public enum Hemisphere {
077        /** North emisphere */
078        North,
079        /** South emisphere */
080        South
081    }
082
083    /**
084     * Constants used for the forward and inverse transform for the elliptical
085     * case of the Transverse Mercator.
086     */
087    private static final double FC1 = 1. / 1;
088    private static final double FC2 = 1. / 2;
089    private static final double FC3 = 1. / 6;
090    private static final double FC4 = 1. / 12;
091    private static final double FC5 = 1. / 20;
092    private static final double FC6 = 1. / 30;
093    private static final double FC7 = 1. / 42;
094    private static final double FC8 = 1. / 56;
095
096    /**
097     * Maximum difference allowed when comparing real numbers.
098     */
099    private static final double EPSILON = 1E-6;
100
101    /**
102     * A derived quantity of excentricity, computed by <code>e'² = (a²-b²)/b² = es/(1-es)</code>
103     * where <var>a</var> is the semi-major axis length and <var>b</var> is the semi-minor axis length.
104     */
105    private double eb2;
106
107    /**
108     * Latitude of origin in <u>radians</u>. Default value is 0, the equator.
109     * This is called '<var>phi0</var>' in Snyder.
110     */
111    protected double latitudeOfOrigin;
112
113    /**
114     * Meridian distance at the {@code latitudeOfOrigin}.
115     * Used for calculations for the ellipsoid.
116     */
117    private double ml0;
118
119    /**
120     * The rectified bearing of the central line, in radians.
121     */
122    protected double rectifiedGridAngle;
123
124    /**
125     * Sine and Cosine values for the coordinate system rotation angle
126     */
127    private double sinrot, cosrot;
128
129    @Override
130    public String getName() {
131        return tr("Transverse Mercator");
132    }
133
134    @Override
135    public String getProj4Id() {
136        return "tmerc";
137    }
138
139    @Override
140    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
141        super.initialize(params);
142        CheckParameterUtil.ensureParameterNotNull(params, "params");
143        CheckParameterUtil.ensureParameterNotNull(params.ellps, "params.ellps");
144        eb2 = params.ellps.eb2;
145        latitudeOfOrigin = params.lat0 == null ? 0 : Utils.toRadians(params.lat0);
146        ml0 = mlfn(latitudeOfOrigin, Math.sin(latitudeOfOrigin), Math.cos(latitudeOfOrigin));
147
148        if (params.gamma != null) {
149            rectifiedGridAngle = Utils.toRadians(params.gamma);
150        } else {
151            rectifiedGridAngle = 0.0;
152        }
153        sinrot = Math.sin(rectifiedGridAngle);
154        cosrot = Math.cos(rectifiedGridAngle);
155    }
156
157    @Override
158    public double[] project(double y, double x) {
159        double sinphi = Math.sin(y);
160        double cosphi = Math.cos(y);
161        double u, v;
162
163        double t = (Math.abs(cosphi) > EPSILON) ? sinphi/cosphi : 0;
164        t *= t;
165        double al = cosphi*x;
166        double als = al*al;
167        al /= Math.sqrt(1.0 - e2 * sinphi*sinphi);
168        double n = eb2 * cosphi*cosphi;
169
170        /* NOTE: meridinal distance at latitudeOfOrigin is always 0 */
171        y = mlfn(y, sinphi, cosphi) - ml0 +
172            sinphi * al * x *
173            FC2 * (1.0 +
174            FC4 * als * (5.0 - t + n*(9.0 + 4.0*n) +
175            FC6 * als * (61.0 + t * (t - 58.0) + n*(270.0 - 330.0*t) +
176            FC8 * als * (1385.0 + t * (t*(543.0 - t) - 3111.0)))));
177
178        x = al*(FC1 + FC3 * als*(1.0 - t + n +
179            FC5 * als * (5.0 + t*(t - 18.0) + n*(14.0 - 58.0*t) +
180            FC7 * als * (61.0+ t*(t*(179.0 - t) - 479.0)))));
181
182        u = y;
183        v = x;
184        x = v * cosrot + u * sinrot;
185        y = u * cosrot - v * sinrot;
186
187        return new double[] {x, y};
188    }
189
190    @Override
191    public double[] invproject(double x, double y) {
192        double v = x * cosrot - y * sinrot;
193        double u = y * cosrot + x * sinrot;
194        x = v;
195        y = u;
196
197        double phi = invMlfn(ml0 + y);
198
199        if (Math.abs(phi) >= Math.PI/2) {
200            y = y < 0.0 ? -(Math.PI/2) : (Math.PI/2);
201            x = 0.0;
202        } else {
203            double sinphi = Math.sin(phi);
204            double cosphi = Math.cos(phi);
205            double t = (Math.abs(cosphi) > EPSILON) ? sinphi/cosphi : 0.0;
206            double n = eb2 * cosphi*cosphi;
207            double con = 1.0 - e2 * sinphi*sinphi;
208            double d = x * Math.sqrt(con);
209            con *= t;
210            t *= t;
211            double ds = d*d;
212
213            y = phi - (con*ds / (1.0 - e2)) *
214                FC2 * (1.0 - ds *
215                FC4 * (5.0 + t*(3.0 - 9.0*n) + n*(1.0 - 4*n) - ds *
216                FC6 * (61.0 + t*(90.0 - 252.0*n + 45.0*t) + 46.0*n - ds *
217                FC8 * (1385.0 + t*(3633.0 + t*(4095.0 + 1575.0*t))))));
218
219            x = d*(FC1 - ds * FC3 * (1.0 + 2.0*t + n -
220                ds*FC5*(5.0 + t*(28.0 + 24* t + 8.0*n) + 6.0*n -
221                ds*FC7*(61.0 + t*(662.0 + t*(1320.0 + 720.0*t))))))/cosphi;
222        }
223        return new double[] {y, x};
224    }
225
226    @Override
227    public Bounds getAlgorithmBounds() {
228        return new Bounds(-89, -7, 89, 7, false);
229    }
230
231    /**
232     * Determines the UTM zone of a given lat/lon.
233     * @param ll lat/lon to locate in the UTM grid.
234     * @return the UTM zone of {@code ll}
235     * @since 13167
236     */
237    public static Pair<Integer, Hemisphere> locateUtmZone(LatLon ll) {
238        return new Pair<>((int) Math.floor((ll.lon() + 180d) / 6d) + 1,
239                ll.lat() > 0 ? Hemisphere.North : Hemisphere.South);
240    }
241}