[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 | } |
---|