source: josm/trunk/src/org/openstreetmap/josm/data/projection/proj/TransverseMercator.java @ 11549

Last change on this file since 11549 was 11549, checked in by bastiK, 7 weeks ago

applied #14346 Rotation angle for Transverse Mercator projection (based on patch by anonymous)

  • Property svn:eol-style set to native
File size: 8.5 KB
RevLine 
[3635]1// License: GPL. For details, see LICENSE file.
[4285]2package org.openstreetmap.josm.data.projection.proj;
[3635]3
[4285]4import static org.openstreetmap.josm.tools.I18n.tr;
5
[9124]6import org.openstreetmap.josm.data.Bounds;
[5066]7import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
[10250]8import org.openstreetmap.josm.tools.CheckParameterUtil;
[4285]9
[3635]10/**
[9112]11 * Transverse Mercator Projection (EPSG code 9807). This
12 * is a cylindrical projection, in which the cylinder has been rotated 90°.
13 * Instead of being tangent to the equator (or to an other standard latitude),
14 * it is tangent to a central meridian. Deformation are more important as we
15 * are going futher from the central meridian. The Transverse Mercator
16 * projection is appropriate for region wich have a greater extent north-south
17 * than east-west.
18 * <p>
[3635]19 *
[9112]20 * The elliptical equations used here are series approximations, and their accuracy
21 * decreases as points move farther from the central meridian of the projection.
22 * The forward equations here are accurate to a less than a mm &plusmn;10 degrees from
23 * the central meridian, a few mm &plusmn;15 degrees from the
24 * central meridian and a few cm &plusmn;20 degrees from the central meridian.
25 * The spherical equations are not approximations and should always give the
26 * correct values.
27 * <p>
[5066]28 *
[9112]29 * There are a number of versions of the transverse mercator projection
30 * including the Universal (UTM) and Modified (MTM) Transverses Mercator
31 * projections. In these cases the earth is divided into zones. For the UTM
32 * the zones are 6 degrees wide, numbered from 1 to 60 proceeding east from
33 * 180 degrees longitude, and between lats 84 degrees North and 80
34 * degrees South. The central meridian is taken as the center of the zone
35 * and the latitude of origin is the equator. A scale factor of 0.9996 and
36 * false easting of 500000m is used for all zones and a false northing of 10000000m
37 * is used for zones in the southern hemisphere.
38 * <p>
39 *
40 * NOTE: formulas used below are not those of Snyder, but rather those
41 *       from the {@code proj4} package of the USGS survey, which
42 *       have been reproduced verbatim. USGS work is acknowledged here.
43 * <p>
44 *
45 * This class has been derived from the implementation of the Geotools project;
46 * git 8cbf52d, org.geotools.referencing.operation.projection.TransverseMercator
47 * at the time of migration.
48 * <p>
[11549]49 * The non-standard parameter <code>gamma</code> has been added as a method
50 * to rotate the projected coordinates by a certain angle (clockwise, see
51 * {@link ObliqueMercator}).
52 * <p>
[9112]53 * <b>References:</b>
54 * <ul>
55 *   <li> Proj-4.4.6 available at <A HREF="http://www.remotesensing.org/proj">www.remotesensing.org/proj</A><br>
56 *        Relevent files are: {@code PJ_tmerc.c}, {@code pj_mlfn.c}, {@code pj_fwd.c} and {@code pj_inv.c}.</li>
57 *   <li> John P. Snyder (Map Projections - A Working Manual,
58 *        U.S. Geological Survey Professional Paper 1395, 1987).</li>
59 *   <li> "Coordinate Conversions and Transformations including Formulas",
60 *        EPSG Guidence Note Number 7, Version 19.</li>
61 * </ul>
62 *
63 * @author André Gosselin
64 * @author Martin Desruisseaux (PMO, IRD)
65 * @author Rueben Schulz
[9117]66 *
67 * @see <A HREF="http://mathworld.wolfram.com/MercatorProjection.html">Transverse Mercator projection on MathWorld</A>
68 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/transverse_mercator.html">"Transverse_Mercator" on RemoteSensing.org</A>
[3635]69 */
[9112]70public class TransverseMercator extends AbstractProj {
[3635]71
[9112]72    /**
73     * Contants used for the forward and inverse transform for the eliptical
74     * case of the Transverse Mercator.
75     */
[9117]76    private static final double FC1 = 1.00000000000000000000000,  // 1/1
77                                FC2 = 0.50000000000000000000000,  // 1/2
78                                FC3 = 0.16666666666666666666666,  // 1/6
79                                FC4 = 0.08333333333333333333333,  // 1/12
80                                FC5 = 0.05000000000000000000000,  // 1/20
81                                FC6 = 0.03333333333333333333333,  // 1/30
82                                FC7 = 0.02380952380952380952380,  // 1/42
83                                FC8 = 0.01785714285714285714285;  // 1/56
[3635]84
[9112]85    /**
86     * Maximum difference allowed when comparing real numbers.
87     */
88    private static final double EPSILON = 1E-6;
89
90    /**
91     * A derived quantity of excentricity, computed by <code>e'² = (a²-b²)/b² = es/(1-es)</code>
[9132]92     * where <var>a</var> is the semi-major axis length and <var>b</var> is the semi-minor axis
[9112]93     * length.
94     */
95    private double eb2;
96
97    /**
98     * Latitude of origin in <u>radians</u>. Default value is 0, the equator.
99     * This is called '<var>phi0</var>' in Snyder.
100     * <p>
101     * <strong>Consider this field as final</strong>. It is not final only
102     * because some classes need to modify it at construction time.
103     */
104    protected double latitudeOfOrigin;
105
106    /**
107     * Meridian distance at the {@code latitudeOfOrigin}.
108     * Used for calculations for the ellipsoid.
109     */
110    private double ml0;
111
[11549]112    /**
113     * The rectified bearing of the central line, in radians.
114     */
115    protected double rectifiedGridAngle;
116
117    /**
118     * Sine and Cosine values for the coordinate system rotation angle
119     */
120    private double sinrot, cosrot;
121
[4285]122    @Override
123    public String getName() {
124        return tr("Transverse Mercator");
[3635]125    }
126
[4285]127    @Override
128    public String getProj4Id() {
129        return "tmerc";
[3635]130    }
131
[5066]132    @Override
133    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
[9112]134        super.initialize(params);
[10250]135        CheckParameterUtil.ensureParameterNotNull(params, "params");
136        CheckParameterUtil.ensureParameterNotNull(params.ellps, "params.ellps");
[9112]137        eb2 = params.ellps.eb2;
138        latitudeOfOrigin = params.lat0 == null ? 0 : Math.toRadians(params.lat0);
139        ml0 = mlfn(latitudeOfOrigin, Math.sin(latitudeOfOrigin), Math.cos(latitudeOfOrigin));
[11549]140
141        if (params.gamma != null) {
142                rectifiedGridAngle = Math.toRadians(params.gamma);
143        } else {
144                rectifiedGridAngle = 0.0;
145        }
146        sinrot = Math.sin(rectifiedGridAngle);
147        cosrot = Math.cos(rectifiedGridAngle);
148
[5066]149    }
150
[4285]151    @Override
[9112]152    public double[] project(double y, double x) {
153        double sinphi = Math.sin(y);
154        double cosphi = Math.cos(y);
[11549]155        double u, v;
[5066]156
[9112]157        double t = (Math.abs(cosphi) > EPSILON) ? sinphi/cosphi : 0;
158        t *= t;
159        double al = cosphi*x;
160        double als = al*al;
161        al /= Math.sqrt(1.0 - e2 * sinphi*sinphi);
162        double n = eb2 * cosphi*cosphi;
[3635]163
[9112]164        /* NOTE: meridinal distance at latitudeOfOrigin is always 0 */
[9970]165        y = mlfn(y, sinphi, cosphi) - ml0 +
[9112]166            sinphi * al * x *
[9117]167            FC2 * (1.0 +
[9112]168            FC4 * als * (5.0 - t + n*(9.0 + 4.0*n) +
169            FC6 * als * (61.0 + t * (t - 58.0) + n*(270.0 - 330.0*t) +
[9970]170            FC8 * als * (1385.0 + t * (t*(543.0 - t) - 3111.0)))));
[3635]171
[9112]172        x = al*(FC1 + FC3 * als*(1.0 - t + n +
173            FC5 * als * (5.0 + t*(t - 18.0) + n*(14.0 - 58.0*t) +
[9117]174            FC7 * als * (61.0+ t*(t*(179.0 - t) - 479.0)))));
[3635]175
[11549]176        u=y; v=x;
177        x = v * cosrot + u * sinrot;
178        y = u * cosrot - v * sinrot;
179
[9117]180        return new double[] {x, y};
[3635]181    }
182
[4285]183    @Override
184    public double[] invproject(double x, double y) {
[11549]185        double v = x * cosrot - y * sinrot;
186        double u = y * cosrot + x * sinrot;
187        x=v; y=u;
188
[10748]189        double phi = invMlfn(ml0 + y);
[3635]190
[9112]191        if (Math.abs(phi) >= Math.PI/2) {
[9117]192            y = y < 0.0 ? -(Math.PI/2) : (Math.PI/2);
[9112]193            x = 0.0;
194        } else {
195            double sinphi = Math.sin(phi);
196            double cosphi = Math.cos(phi);
197            double t = (Math.abs(cosphi) > EPSILON) ? sinphi/cosphi : 0.0;
198            double n = eb2 * cosphi*cosphi;
199            double con = 1.0 - e2 * sinphi*sinphi;
200            double d = x * Math.sqrt(con);
201            con *= t;
202            t *= t;
203            double ds = d*d;
[3635]204
[9112]205            y = phi - (con*ds / (1.0 - e2)) *
206                FC2 * (1.0 - ds *
207                FC4 * (5.0 + t*(3.0 - 9.0*n) + n*(1.0 - 4*n) - ds *
208                FC6 * (61.0 + t*(90.0 - 252.0*n + 45.0*t) + 46.0*n - ds *
209                FC8 * (1385.0 + t*(3633.0 + t*(4095.0 + 1574.0*t))))));
[3635]210
[9112]211            x = d*(FC1 - ds * FC3 * (1.0 + 2.0*t + n -
212                ds*FC5*(5.0 + t*(28.0 + 24* t + 8.0*n) + 6.0*n -
213                ds*FC7*(61.0 + t*(662.0 + t*(1320.0 + 720.0*t))))))/cosphi;
214        }
[9117]215        return new double[] {y, x};
[3635]216    }
[9124]217
218    @Override
219    public Bounds getAlgorithmBounds() {
[9139]220        return new Bounds(-89, -7, 89, 7, false);
[9124]221    }
[3635]222}
Note: See TracBrowser for help on using the repository browser.