[3635] | 1 | // License: GPL. For details, see LICENSE file.
|
---|
[4285] | 2 | package org.openstreetmap.josm.data.projection.proj;
|
---|
[3635] | 3 |
|
---|
[4285] | 4 | import static org.openstreetmap.josm.tools.I18n.tr;
|
---|
| 5 |
|
---|
[9124] | 6 | import org.openstreetmap.josm.data.Bounds;
|
---|
[5066] | 7 | import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
|
---|
[10250] | 8 | import 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 ±10 degrees from
|
---|
| 23 | * the central meridian, a few mm ±15 degrees from the
|
---|
| 24 | * central meridian and a few cm ±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] | 70 | public 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 | }
|
---|