1 | // License: GPL. For details, see LICENSE file.
|
---|
2 | package org.openstreetmap.josm.data.projection.proj;
|
---|
3 |
|
---|
4 | import static org.openstreetmap.josm.tools.I18n.tr;
|
---|
5 |
|
---|
6 | import org.openstreetmap.josm.data.Bounds;
|
---|
7 | import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
|
---|
8 |
|
---|
9 | /**
|
---|
10 | * Transverse Mercator Projection (EPSG code 9807). This
|
---|
11 | * is a cylindrical projection, in which the cylinder has been rotated 90°.
|
---|
12 | * Instead of being tangent to the equator (or to an other standard latitude),
|
---|
13 | * it is tangent to a central meridian. Deformation are more important as we
|
---|
14 | * are going futher from the central meridian. The Transverse Mercator
|
---|
15 | * projection is appropriate for region wich have a greater extent north-south
|
---|
16 | * than east-west.
|
---|
17 | * <p>
|
---|
18 | *
|
---|
19 | * The elliptical equations used here are series approximations, and their accuracy
|
---|
20 | * decreases as points move farther from the central meridian of the projection.
|
---|
21 | * The forward equations here are accurate to a less than a mm ±10 degrees from
|
---|
22 | * the central meridian, a few mm ±15 degrees from the
|
---|
23 | * central meridian and a few cm ±20 degrees from the central meridian.
|
---|
24 | * The spherical equations are not approximations and should always give the
|
---|
25 | * correct values.
|
---|
26 | * <p>
|
---|
27 | *
|
---|
28 | * There are a number of versions of the transverse mercator projection
|
---|
29 | * including the Universal (UTM) and Modified (MTM) Transverses Mercator
|
---|
30 | * projections. In these cases the earth is divided into zones. For the UTM
|
---|
31 | * the zones are 6 degrees wide, numbered from 1 to 60 proceeding east from
|
---|
32 | * 180 degrees longitude, and between lats 84 degrees North and 80
|
---|
33 | * degrees South. The central meridian is taken as the center of the zone
|
---|
34 | * and the latitude of origin is the equator. A scale factor of 0.9996 and
|
---|
35 | * false easting of 500000m is used for all zones and a false northing of 10000000m
|
---|
36 | * is used for zones in the southern hemisphere.
|
---|
37 | * <p>
|
---|
38 | *
|
---|
39 | * NOTE: formulas used below are not those of Snyder, but rather those
|
---|
40 | * from the {@code proj4} package of the USGS survey, which
|
---|
41 | * have been reproduced verbatim. USGS work is acknowledged here.
|
---|
42 | * <p>
|
---|
43 | *
|
---|
44 | * This class has been derived from the implementation of the Geotools project;
|
---|
45 | * git 8cbf52d, org.geotools.referencing.operation.projection.TransverseMercator
|
---|
46 | * at the time of migration.
|
---|
47 | * <p>
|
---|
48 | *
|
---|
49 | * <b>References:</b>
|
---|
50 | * <ul>
|
---|
51 | * <li> Proj-4.4.6 available at <A HREF="http://www.remotesensing.org/proj">www.remotesensing.org/proj</A><br>
|
---|
52 | * Relevent files are: {@code PJ_tmerc.c}, {@code pj_mlfn.c}, {@code pj_fwd.c} and {@code pj_inv.c}.</li>
|
---|
53 | * <li> John P. Snyder (Map Projections - A Working Manual,
|
---|
54 | * U.S. Geological Survey Professional Paper 1395, 1987).</li>
|
---|
55 | * <li> "Coordinate Conversions and Transformations including Formulas",
|
---|
56 | * EPSG Guidence Note Number 7, Version 19.</li>
|
---|
57 | * </ul>
|
---|
58 | *
|
---|
59 | * @author André Gosselin
|
---|
60 | * @author Martin Desruisseaux (PMO, IRD)
|
---|
61 | * @author Rueben Schulz
|
---|
62 | *
|
---|
63 | * @see <A HREF="http://mathworld.wolfram.com/MercatorProjection.html">Transverse Mercator projection on MathWorld</A>
|
---|
64 | * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/transverse_mercator.html">"Transverse_Mercator" on RemoteSensing.org</A>
|
---|
65 | */
|
---|
66 | public class TransverseMercator extends AbstractProj {
|
---|
67 |
|
---|
68 | /**
|
---|
69 | * Contants used for the forward and inverse transform for the eliptical
|
---|
70 | * case of the Transverse Mercator.
|
---|
71 | */
|
---|
72 | private static final double FC1 = 1.00000000000000000000000, // 1/1
|
---|
73 | FC2 = 0.50000000000000000000000, // 1/2
|
---|
74 | FC3 = 0.16666666666666666666666, // 1/6
|
---|
75 | FC4 = 0.08333333333333333333333, // 1/12
|
---|
76 | FC5 = 0.05000000000000000000000, // 1/20
|
---|
77 | FC6 = 0.03333333333333333333333, // 1/30
|
---|
78 | FC7 = 0.02380952380952380952380, // 1/42
|
---|
79 | FC8 = 0.01785714285714285714285; // 1/56
|
---|
80 |
|
---|
81 | /**
|
---|
82 | * Maximum difference allowed when comparing real numbers.
|
---|
83 | */
|
---|
84 | private static final double EPSILON = 1E-6;
|
---|
85 |
|
---|
86 | /**
|
---|
87 | * A derived quantity of excentricity, computed by <code>e'² = (a²-b²)/b² = es/(1-es)</code>
|
---|
88 | * where <var>a</var> is the semi-major axis length and <var>b</var> is the semi-minor axis
|
---|
89 | * length.
|
---|
90 | */
|
---|
91 | private double eb2;
|
---|
92 |
|
---|
93 | /**
|
---|
94 | * Latitude of origin in <u>radians</u>. Default value is 0, the equator.
|
---|
95 | * This is called '<var>phi0</var>' in Snyder.
|
---|
96 | * <p>
|
---|
97 | * <strong>Consider this field as final</strong>. It is not final only
|
---|
98 | * because some classes need to modify it at construction time.
|
---|
99 | */
|
---|
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 | }
|
---|