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 futher from the central meridian. The Transverse Mercator 019 * projection is appropriate for region wich 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 ±10 degrees from 026 * the central meridian, a few mm ±15 degrees from the 027 * central meridian and a few cm ±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 * Contants used for the forward and inverse transform for the eliptical 085 * case of the Transverse Mercator. 086 */ 087 private static final double FC1 = 1.00000000000000000000000, // 1/1 088 FC2 = 0.50000000000000000000000, // 1/2 089 FC3 = 0.16666666666666666666666, // 1/6 090 FC4 = 0.08333333333333333333333, // 1/12 091 FC5 = 0.05000000000000000000000, // 1/20 092 FC6 = 0.03333333333333333333333, // 1/30 093 FC7 = 0.02380952380952380952380, // 1/42 094 FC8 = 0.01785714285714285714285; // 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 104 * length. 105 */ 106 private double eb2; 107 108 /** 109 * Latitude of origin in <u>radians</u>. Default value is 0, the equator. 110 * This is called '<var>phi0</var>' in Snyder. 111 * <p> 112 * <strong>Consider this field as final</strong>. It is not final only 113 * because some classes need to modify it at construction time. 114 */ 115 protected double latitudeOfOrigin; 116 117 /** 118 * Meridian distance at the {@code latitudeOfOrigin}. 119 * Used for calculations for the ellipsoid. 120 */ 121 private double ml0; 122 123 /** 124 * The rectified bearing of the central line, in radians. 125 */ 126 protected double rectifiedGridAngle; 127 128 /** 129 * Sine and Cosine values for the coordinate system rotation angle 130 */ 131 private double sinrot, cosrot; 132 133 @Override 134 public String getName() { 135 return tr("Transverse Mercator"); 136 } 137 138 @Override 139 public String getProj4Id() { 140 return "tmerc"; 141 } 142 143 @Override 144 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 145 super.initialize(params); 146 CheckParameterUtil.ensureParameterNotNull(params, "params"); 147 CheckParameterUtil.ensureParameterNotNull(params.ellps, "params.ellps"); 148 eb2 = params.ellps.eb2; 149 latitudeOfOrigin = params.lat0 == null ? 0 : Utils.toRadians(params.lat0); 150 ml0 = mlfn(latitudeOfOrigin, Math.sin(latitudeOfOrigin), Math.cos(latitudeOfOrigin)); 151 152 if (params.gamma != null) { 153 rectifiedGridAngle = Utils.toRadians(params.gamma); 154 } else { 155 rectifiedGridAngle = 0.0; 156 } 157 sinrot = Math.sin(rectifiedGridAngle); 158 cosrot = Math.cos(rectifiedGridAngle); 159 160 } 161 162 @Override 163 public double[] project(double y, double x) { 164 double sinphi = Math.sin(y); 165 double cosphi = Math.cos(y); 166 double u, v; 167 168 double t = (Math.abs(cosphi) > EPSILON) ? sinphi/cosphi : 0; 169 t *= t; 170 double al = cosphi*x; 171 double als = al*al; 172 al /= Math.sqrt(1.0 - e2 * sinphi*sinphi); 173 double n = eb2 * cosphi*cosphi; 174 175 /* NOTE: meridinal distance at latitudeOfOrigin is always 0 */ 176 y = mlfn(y, sinphi, cosphi) - ml0 + 177 sinphi * al * x * 178 FC2 * (1.0 + 179 FC4 * als * (5.0 - t + n*(9.0 + 4.0*n) + 180 FC6 * als * (61.0 + t * (t - 58.0) + n*(270.0 - 330.0*t) + 181 FC8 * als * (1385.0 + t * (t*(543.0 - t) - 3111.0))))); 182 183 x = al*(FC1 + FC3 * als*(1.0 - t + n + 184 FC5 * als * (5.0 + t*(t - 18.0) + n*(14.0 - 58.0*t) + 185 FC7 * als * (61.0+ t*(t*(179.0 - t) - 479.0))))); 186 187 u = y; 188 v = x; 189 x = v * cosrot + u * sinrot; 190 y = u * cosrot - v * sinrot; 191 192 return new double[] {x, y}; 193 } 194 195 @Override 196 public double[] invproject(double x, double y) { 197 double v = x * cosrot - y * sinrot; 198 double u = y * cosrot + x * sinrot; 199 x = v; 200 y = u; 201 202 double phi = invMlfn(ml0 + y); 203 204 if (Math.abs(phi) >= Math.PI/2) { 205 y = y < 0.0 ? -(Math.PI/2) : (Math.PI/2); 206 x = 0.0; 207 } else { 208 double sinphi = Math.sin(phi); 209 double cosphi = Math.cos(phi); 210 double t = (Math.abs(cosphi) > EPSILON) ? sinphi/cosphi : 0.0; 211 double n = eb2 * cosphi*cosphi; 212 double con = 1.0 - e2 * sinphi*sinphi; 213 double d = x * Math.sqrt(con); 214 con *= t; 215 t *= t; 216 double ds = d*d; 217 218 y = phi - (con*ds / (1.0 - e2)) * 219 FC2 * (1.0 - ds * 220 FC4 * (5.0 + t*(3.0 - 9.0*n) + n*(1.0 - 4*n) - ds * 221 FC6 * (61.0 + t*(90.0 - 252.0*n + 45.0*t) + 46.0*n - ds * 222 FC8 * (1385.0 + t*(3633.0 + t*(4095.0 + 1574.0*t)))))); 223 224 x = d*(FC1 - ds * FC3 * (1.0 + 2.0*t + n - 225 ds*FC5*(5.0 + t*(28.0 + 24* t + 8.0*n) + 6.0*n - 226 ds*FC7*(61.0 + t*(662.0 + t*(1320.0 + 720.0*t))))))/cosphi; 227 } 228 return new double[] {y, x}; 229 } 230 231 @Override 232 public Bounds getAlgorithmBounds() { 233 return new Bounds(-89, -7, 89, 7, false); 234 } 235 236 /** 237 * Determines the UTM zone of a given lat/lon. 238 * @param ll lat/lon to locate in the UTM grid. 239 * @return the UTM zone of {@code ll} 240 * @since 13167 241 */ 242 public static Pair<Integer, Hemisphere> locateUtmZone(LatLon ll) { 243 return new Pair<>((int) Math.floor((ll.lon() + 180d) / 6d) + 1, 244 ll.lat() > 0 ? Hemisphere.North : Hemisphere.South); 245 } 246}