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