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.Utils; 010 011/** 012 * Oblique Mercator Projection. A conformal, oblique, cylindrical projection with the cylinder 013 * touching the ellipsoid (or sphere) along a great circle path (the central line). The 014 * {@linkplain Mercator} and {@linkplain TransverseMercator Transverse Mercator} projections can 015 * be thought of as special cases of the oblique mercator, where the central line is along the 016 * equator or a meridian, respectively. The Oblique Mercator projection has been used in 017 * Switzerland, Hungary, Madagascar, Malaysia, Borneo and the panhandle of Alaska. 018 * <p> 019 * The Oblique Mercator projection uses a (<var>U</var>,<var>V</var>) coordinate system, with the 020 * <var>U</var> axis along the central line. During the forward projection, coordinates from the 021 * ellipsoid are projected conformally to a sphere of constant total curvature, called the 022 * "aposphere", before being projected onto the plane. The projection coordinates are further 023 * convented to a (<var>X</var>,<var>Y</var>) coordinate system by rotating the calculated 024 * (<var>u</var>,<var>v</var>) coordinates to give output (<var>x</var>,<var>y</var>) coordinates. 025 * The rotation value is usually the same as the projection azimuth (the angle, east of north, of 026 * the central line), but some cases allow a separate rotation parameter. 027 * <p> 028 * There are two forms of the oblique mercator, differing in the origin of their grid coordinates. 029 * The Hotine Oblique Mercator (EPSG code 9812) has grid coordinates start at the intersection of 030 * the central line and the equator of the aposphere. 031 * The Oblique Mercator (EPSG code 9815) is the same, except the grid coordinates begin at the 032 * central point (where the latitude of center and central line intersect). ESRI separates these 033 * two case by appending {@code "Natural_Origin"} (for the {@code "Hotine_Oblique_Mercator"}) and 034 * {@code "Center"} (for the {@code "Oblique_Mercator"}) to the projection names. 035 * <p> 036 * Two different methods are used to specify the central line for the oblique mercator: 037 * 1) a central point and an azimuth, east of north, describing the central line and 038 * 2) two points on the central line. The EPSG does not use the two point method, 039 * while ESRI separates the two cases by putting {@code "Azimuth"} and {@code "Two_Point"} 040 * in their projection names. Both cases use the point where the {@code "latitude_of_center"} 041 * parameter crosses the central line as the projection's central point. 042 * The {@linkplain #centralMeridian central meridian} is not a projection parameter, 043 * and is instead calculated as the intersection between the central line and the 044 * equator of the aposphere. 045 * <p> 046 * For the azimuth method, the central latitude cannot be ±90.0 degrees 047 * and the central line cannot be at a maximum or minimum latitude at the central point. 048 * In the two point method, the latitude of the first and second points cannot be 049 * equal. Also, the latitude of the first point and central point cannot be 050 * ±90.0 degrees. Furthermore, the latitude of the first point cannot be 0.0 and 051 * the latitude of the second point cannot be -90.0 degrees. A change of 052 * 10<sup>-7</sup> radians can allow calculation at these special cases. Snyder's restriction 053 * of the central latitude being 0.0 has been removed, since the equations appear 054 * to work correctly in this case. 055 * <p> 056 * Azimuth values of 0.0 and ±90.0 degrees are allowed (and used in Hungary 057 * and Switzerland), though these cases would usually use a Mercator or 058 * Transverse Mercator projection instead. Azimuth values > 90 degrees cause 059 * errors in the equations. 060 * <p> 061 * The oblique mercator is also called the "Rectified Skew Orthomorphic" (RSO). It appears 062 * is that the only difference from the oblique mercator is that the RSO allows the rotation 063 * from the (<var>U</var>,<var>V</var>) to (<var>X</var>,<var>Y</var>) coordinate system to 064 * be different from the azimuth. This separate parameter is called 065 * {@code "rectified_grid_angle"} (or {@code "XY_Plane_Rotation"} by ESRI) and is also 066 * included in the EPSG's parameters for the Oblique Mercator and Hotine Oblique Mercator. 067 * The rotation parameter is optional in all the non-two point projections and will be 068 * set to the azimuth if not specified. 069 * <p> 070 * Projection cases and aliases implemented by the {@link ObliqueMercator} are: 071 * <ul> 072 * <li>{@code Oblique_Mercator} (EPSG code 9815)<br> 073 * grid coordinates begin at the central point, 074 * has {@code "rectified_grid_angle"} parameter.</li> 075 * <li>{@code Hotine_Oblique_Mercator_Azimuth_Center} (ESRI)<br> 076 * grid coordinates begin at the central point.</li> 077 * <li>{@code Rectified_Skew_Orthomorphic_Center} (ESRI)<br> 078 * grid coordinates begin at the central point, 079 * has {@code "rectified_grid_angle"} parameter.</li> 080 * 081 * <li>{@code Hotine_Oblique_Mercator} (EPSG code 9812)<br> 082 * grid coordinates begin at the intersection of the central line and aposphere equator, 083 * has {@code "rectified_grid_angle"} parameter.</li> 084 * <li>{@code Hotine_Oblique_Mercator_Azimuth_Natural_Origin} (ESRI)<br> 085 * grid coordinates begin at the interseciton of the central line and aposphere equator.</li> 086 * <li>{@code Rectified_Skew_Orthomorphic_Natural_Origin} (ESRI)<br> 087 * grid coordinates begin at the interseciton of the central line and aposphere equator, 088 * has {@code "rectified_grid_angle"} parameter.</li> 089 * 090 * <li>{@code Hotine_Oblique_Mercator_Two_Point_Center} (ESRI)<br> 091 * grid coordinates begin at the central point.</li> 092 * <li>{@code Hotine_Oblique_Mercator_Two_Point_Natural_Origin} (ESRI)<br> 093 * grid coordinates begin at the interseciton of the central line and aposphere equator.</li> 094 * </ul> 095 * <p> 096 * This class has been derived from the implementation of the Geotools project; 097 * git 8cbf52d, org.geotools.referencing.operation.projection.ObliqueMercator 098 * at the time of migration. 099 * <p> 100 * Note that automatic calculation of bounds is very limited for this projection, 101 * since the central line can have any orientation. 102 * <p> 103 * <b>References:</b> 104 * <ul> 105 * <li>{@code libproj4} is available at 106 * <A HREF="http://members.bellatlantic.net/~vze2hc4d/proj4/">libproj4 Miscellanea</A><br> 107 * Relevent files are: {@code PJ_omerc.c}, {@code pj_tsfn.c}, 108 * {@code pj_fwd.c}, {@code pj_inv.c} and {@code lib_proj.h}</li> 109 * <li>John P. Snyder (Map Projections - A Working Manual, 110 * U.S. Geological Survey Professional Paper 1395, 1987)</li> 111 * <li>"Coordinate Conversions and Transformations including Formulas", 112 * EPSG Guidence Note Number 7 part 2, Version 24.</li> 113 * <li>Gerald Evenden, 2004, <a href="http://members.verizon.net/~vze2hc4d/proj4/omerc.pdf"> 114 * Documentation of revised Oblique Mercator</a></li> 115 * </ul> 116 * 117 * @author Gerald I. Evenden (for original code in Proj4) 118 * @author Rueben Schulz 119 * 120 * @see <A HREF="http://mathworld.wolfram.com/MercatorProjection.html">Oblique Mercator projection on MathWorld</A> 121 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/hotine_oblique_mercator.html">"hotine_oblique_mercator" on RemoteSensing.org</A> 122 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/oblique_mercator.html">"oblique_mercator" on RemoteSensing.org</A> 123 */ 124public class ObliqueMercator extends AbstractProj implements ICentralMeridianProvider { 125 126 /** 127 * Maximum difference allowed when comparing real numbers. 128 */ 129 private static final double EPSILON = 1E-6; 130 131 /** 132 * Maximum difference allowed when comparing latitudes. 133 */ 134 private static final double EPSILON_LATITUDE = 1E-10; 135 136 ////// 137 ////// Map projection parameters. The following are NOT used by the transformation 138 ////// methods, but are stored in order to restitute them in WKT formatting. They 139 ////// are made visible ('protected' access) for documentation purpose and because 140 ////// they are user-supplied parameters, not derived coefficients. 141 ////// 142 143 /** 144 * The azimuth of the central line passing through the centre of the projection, in radians. 145 */ 146 protected double azimuth; 147 148 /** 149 * The rectified bearing of the central line, in radians. This is equals to the 150 * {@linkplain #azimuth} if the parameter value is not set. 151 */ 152 protected double rectifiedGridAngle; 153 154 ////// 155 ////// Map projection coefficients computed from the above parameters. 156 ////// They are the fields used for coordinate transformations. 157 ////// 158 159 /** 160 * Constants used in the transformation. 161 */ 162 private double b, g; 163 164 /** 165 * Convenience value equal to {@code a} / {@link #b}. 166 */ 167 private double arb; 168 169 /** 170 * Convenience value equal to {@code a}×{@link #b}. 171 */ 172 private double ab; 173 174 /** 175 * Convenience value equal to {@link #b} / {@code a}. 176 */ 177 private double bra; 178 179 /** 180 * <var>v</var> values when the input latitude is a pole. 181 */ 182 private double vPoleN, vPoleS; 183 184 /** 185 * Sine and Cosine values for gamma0 (the angle between the meridian 186 * and central line at the intersection between the central line and 187 * the Earth equator on aposphere). 188 */ 189 private double singamma0, cosgamma0; 190 191 /** 192 * Sine and Cosine values for the rotation between (U,V) and 193 * (X,Y) coordinate systems 194 */ 195 private double sinrot, cosrot; 196 197 /** 198 * <var>u</var> value (in (U,V) coordinate system) of the central point. Used in 199 * the oblique mercator case. The <var>v</var> value of the central point is 0.0. 200 */ 201 private double uc; 202 203 /** 204 * Central longitude in <u>radians</u>. Default value is 0, the Greenwich meridian. 205 * This is called '<var>lambda0</var>' in Snyder. 206 */ 207 protected double centralMeridian; 208 209 /** 210 * A reference point, which is known to be on the central line. 211 */ 212 private LatLon referencePoint; 213 214 @Override 215 public String getName() { 216 return tr("Oblique Mercator"); 217 } 218 219 @Override 220 public String getProj4Id() { 221 return "omerc"; 222 } 223 224 @Override 225 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 226 super.initialize(params); 227 boolean twoPoint = params.alpha == null; 228 229 double latCenter = 0; 230 if (params.lat0 != null) { 231 latCenter = Utils.toRadians(params.lat0); 232 } 233 234 final double com = Math.sqrt(1.0 - e2); 235 double sinph0 = Math.sin(latCenter); 236 double cosph0 = Math.cos(latCenter); 237 final double con = 1. - e2 * sinph0 * sinph0; 238 double temp = cosph0 * cosph0; 239 b = Math.sqrt(1.0 + e2 * (temp * temp) / (1.0 - e2)); 240 double a = b * com / con; 241 final double d = b * com / (cosph0 * Math.sqrt(con)); 242 double f = d * d - 1.0; 243 if (f < 0.0) { 244 f = 0.0; 245 } else { 246 f = Math.sqrt(f); 247 if (latCenter < 0.0) { 248 f = -f; 249 } 250 } 251 g = f += d; 252 g = f * Math.pow(tsfn(latCenter, sinph0), b); 253 254 /* 255 * Computes the constants that depend on the "twoPoint" vs "azimuth" case. In the 256 * two points case, we compute them from (LAT_OF_1ST_POINT, LONG_OF_1ST_POINT) and 257 * (LAT_OF_2ND_POINT, LONG_OF_2ND_POINT). For the "azimuth" case, we compute them 258 * from LONGITUDE_OF_CENTRE and AZIMUTH. 259 */ 260 final double gamma0; 261 Double lonCenter = null; 262 if (twoPoint) { 263 if (params.lon1 == null) 264 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lon_1")); 265 if (params.lat1 == null) 266 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_1")); 267 if (params.lon2 == null) 268 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lon_2")); 269 if (params.lat2 == null) 270 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_2")); 271 referencePoint = new LatLon(params.lat1, params.lat2); 272 double lon1 = Utils.toRadians(params.lon1); 273 double lat1 = Utils.toRadians(params.lat1); 274 double lon2 = Utils.toRadians(params.lon2); 275 double lat2 = Utils.toRadians(params.lat2); 276 277 if (Math.abs(lat1 - lat2) <= EPSILON || 278 Math.abs(lat1) <= EPSILON || 279 Math.abs(Math.abs(lat1) - Math.PI / 2) <= EPSILON || 280 Math.abs(Math.abs(latCenter) - Math.PI / 2) <= EPSILON || 281 Math.abs(Math.abs(lat2) - Math.PI / 2) <= EPSILON) { 282 throw new ProjectionConfigurationException( 283 tr("Unsuitable parameters ''{0}'' and ''{1}'' for two point method.", "lat_1", "lat_2")); 284 } 285 /* 286 * The coefficients for the "two points" case. 287 */ 288 final double h = Math.pow(tsfn(lat1, Math.sin(lat1)), b); 289 final double l = Math.pow(tsfn(lat2, Math.sin(lat2)), b); 290 final double fp = g / h; 291 final double p = (l - h) / (l + h); 292 double j = g * g; 293 j = (j - l * h) / (j + l * h); 294 double diff = lon1 - lon2; 295 if (diff < -Math.PI) { 296 lon2 -= 2.0 * Math.PI; 297 } else if (diff > Math.PI) { 298 lon2 += 2.0 * Math.PI; 299 } 300 centralMeridian = normalizeLonRad(0.5 * (lon1 + lon2) - 301 Math.atan(j * Math.tan(0.5 * b * (lon1 - lon2)) / p) / b); 302 gamma0 = Math.atan(2.0 * Math.sin(b * normalizeLonRad(lon1 - centralMeridian)) / 303 (fp - 1.0 / fp)); 304 azimuth = Math.asin(d * Math.sin(gamma0)); 305 rectifiedGridAngle = azimuth; 306 } else { 307 if (params.lonc == null) 308 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lonc")); 309 if (params.lat0 == null) 310 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); 311 if (params.alpha == null) 312 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "alpha")); 313 referencePoint = new LatLon(params.lat0, params.lonc); 314 315 lonCenter = Utils.toRadians(params.lonc); 316 azimuth = Utils.toRadians(params.alpha); 317 // CHECKSTYLE.OFF: SingleSpaceSeparator 318 if ((azimuth > -1.5*Math.PI && azimuth < -0.5*Math.PI) || 319 (azimuth > 0.5*Math.PI && azimuth < 1.5*Math.PI)) { 320 throw new ProjectionConfigurationException( 321 tr("Illegal value for parameter ''{0}'': {1}", "alpha", Double.toString(params.alpha))); 322 } 323 // CHECKSTYLE.ON: SingleSpaceSeparator 324 if (params.gamma != null) { 325 rectifiedGridAngle = Utils.toRadians(params.gamma); 326 } else { 327 rectifiedGridAngle = azimuth; 328 } 329 gamma0 = Math.asin(Math.sin(azimuth) / d); 330 // Check for asin(+-1.00000001) 331 temp = 0.5 * (f - 1.0 / f) * Math.tan(gamma0); 332 if (Math.abs(temp) > 1.0) { 333 if (Math.abs(Math.abs(temp) - 1.0) > EPSILON) { 334 throw new ProjectionConfigurationException(tr("error in initialization")); 335 } 336 temp = (temp > 0) ? 1.0 : -1.0; 337 } 338 centralMeridian = lonCenter - Math.asin(temp) / b; 339 } 340 341 /* 342 * More coefficients common to all kind of oblique mercator. 343 */ 344 singamma0 = Math.sin(gamma0); 345 cosgamma0 = Math.cos(gamma0); 346 sinrot = Math.sin(rectifiedGridAngle); 347 cosrot = Math.cos(rectifiedGridAngle); 348 arb = a / b; 349 ab = a * b; 350 bra = b / a; 351 vPoleN = arb * Math.log(Math.tan(0.5 * (Math.PI/2.0 - gamma0))); 352 vPoleS = arb * Math.log(Math.tan(0.5 * (Math.PI/2.0 + gamma0))); 353 boolean hotine = params.no_off != null && params.no_off; 354 if (hotine) { 355 uc = 0.0; 356 } else { 357 if (Math.abs(Math.abs(azimuth) - Math.PI/2.0) < EPSILON_LATITUDE) { 358 // lonCenter == null in twoPoint, but azimuth cannot be 90 here (lat1 != lat2) 359 if (lonCenter == null) { 360 throw new ProjectionConfigurationException("assertion error"); 361 } 362 uc = a * (lonCenter - centralMeridian); 363 } else { 364 double uC = Math.abs(arb * Math.atan2(Math.sqrt(d * d - 1.0), Math.cos(azimuth))); 365 if (latCenter < 0.0) { 366 uC = -uC; 367 } 368 this.uc = uC; 369 } 370 } 371 } 372 373 private static double normalizeLonRad(double a) { 374 return Utils.toRadians(LatLon.normalizeLon(Utils.toDegrees(a))); 375 } 376 377 @Override 378 public double[] project(double y, double x) { 379 double u, v; 380 if (Math.abs(Math.abs(y) - Math.PI/2.0) > EPSILON) { 381 double q = g / Math.pow(tsfn(y, Math.sin(y)), b); 382 double temp = 1.0 / q; 383 double s = 0.5 * (q - temp); 384 double v2 = Math.sin(b * x); 385 double u2 = (s * singamma0 - v2 * cosgamma0) / (0.5 * (q + temp)); 386 if (Math.abs(Math.abs(u2) - 1.0) < EPSILON) { 387 v = 0; // this is actually an error and should be reported to the caller somehow 388 } else { 389 v = 0.5 * arb * Math.log((1.0 - u2) / (1.0 + u2)); 390 } 391 temp = Math.cos(b * x); 392 if (Math.abs(temp) < EPSILON_LATITUDE) { 393 u = ab * x; 394 } else { 395 u = arb * Math.atan2(s * cosgamma0 + v2 * singamma0, temp); 396 } 397 } else { 398 v = y > 0 ? vPoleN : vPoleS; 399 u = arb * y; 400 } 401 u -= uc; 402 x = v * cosrot + u * sinrot; 403 y = u * cosrot - v * sinrot; 404 return new double[] {x, y}; 405 } 406 407 @Override 408 public double[] invproject(double x, double y) { 409 double v = x * cosrot - y * sinrot; 410 double u = y * cosrot + x * sinrot + uc; 411 double qp = Math.exp(-bra * v); 412 double temp = 1.0 / qp; 413 double sp = 0.5 * (qp - temp); 414 double vp = Math.sin(bra * u); 415 double up = (vp * cosgamma0 + sp * singamma0) / (0.5 * (qp + temp)); 416 if (Math.abs(Math.abs(up) - 1.0) < EPSILON) { 417 return new double[] { 418 up < 0.0 ? -(Math.PI / 2.0) : (Math.PI / 2.0), 419 0.0}; 420 } else { 421 return new double[] { 422 cphi2(Math.pow(g / Math.sqrt((1. + up) / (1. - up)), 1.0 / b)), //calculate t 423 -Math.atan2(sp * cosgamma0 - vp * singamma0, Math.cos(bra * u)) / b}; 424 } 425 } 426 427 @Override 428 public Bounds getAlgorithmBounds() { 429 // The central line of this projection can be oriented in any direction. 430 // Moreover, the projection doesn't work too well very far off the central line. 431 // This makes it hard to choose proper bounds automatically. 432 // 433 // We return a small box around a reference point. This default box is 434 // probably too small for most applications. If this is the case, the 435 // bounds should be configured explicitly. 436 double lat = referencePoint.lat(); 437 double dLat = 3.0; 438 double lon = referencePoint.lon() - Utils.toDegrees(centralMeridian); 439 double dLon = 3.0; 440 return new Bounds(lat - dLat, lon - dLon, lat + dLat, lon + dLon, false); 441 } 442 443 @Override 444 public double getCentralMeridian() { 445 return Utils.toDegrees(centralMeridian); 446 } 447}