source: josm/trunk/src/org/openstreetmap/josm/data/projection/proj/ObliqueMercator.java @ 12445

Last change on this file since 12445 was 12445, checked in by Don-vip, 21 months ago

update to error-prone 2.0.21, groovy 2.4.12

  • Property svn:eol-style set to native
File size: 19.7 KB
Line 
1// License: GPL. For details, see LICENSE file.
2package org.openstreetmap.josm.data.projection.proj;
3
4import static org.openstreetmap.josm.tools.I18n.tr;
5
6import org.openstreetmap.josm.data.Bounds;
7import org.openstreetmap.josm.data.coor.LatLon;
8import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
9import org.openstreetmap.josm.tools.Utils;
10
11/**
12 * Oblique Mercator Projection. A conformal, oblique, cylindrical projection with the cylinder
13 * touching the ellipsoid (or sphere) along a great circle path (the central line). The
14 * {@linkplain Mercator} and {@linkplain TransverseMercator Transverse Mercator} projections can
15 * be thought of as special cases of the oblique mercator, where the central line is along the
16 * equator or a meridian, respectively. The Oblique Mercator projection has been used in
17 * Switzerland, Hungary, Madagascar, Malaysia, Borneo and the panhandle of Alaska.
18 * <p>
19 * The Oblique Mercator projection uses a (<var>U</var>,<var>V</var>) coordinate system, with the
20 * <var>U</var> axis along the central line. During the forward projection, coordinates from the
21 * ellipsoid are projected conformally to a sphere of constant total curvature, called the
22 * "aposphere", before being projected onto the plane. The projection coordinates are further
23 * convented to a (<var>X</var>,<var>Y</var>) coordinate system by rotating the calculated
24 * (<var>u</var>,<var>v</var>) coordinates to give output (<var>x</var>,<var>y</var>) coordinates.
25 * The rotation value is usually the same as the projection azimuth (the angle, east of north, of
26 * the central line), but some cases allow a separate rotation parameter.
27 * <p>
28 * There are two forms of the oblique mercator, differing in the origin of their grid coordinates.
29 * The Hotine Oblique Mercator (EPSG code 9812) has grid coordinates start at the intersection of
30 * the central line and the equator of the aposphere.
31 * The Oblique Mercator (EPSG code 9815) is the same, except the grid coordinates begin at the
32 * central point (where the latitude of center and central line intersect). ESRI separates these
33 * two case by appending {@code "Natural_Origin"} (for the {@code "Hotine_Oblique_Mercator"}) and
34 * {@code "Center"} (for the {@code "Oblique_Mercator"}) to the projection names.
35 * <p>
36 * Two different methods are used to specify the central line for the oblique mercator:
37 * 1) a central point and an azimuth, east of north, describing the central line and
38 * 2) two points on the central line. The EPSG does not use the two point method,
39 * while ESRI separates the two cases by putting {@code "Azimuth"} and {@code "Two_Point"}
40 * in their projection names. Both cases use the point where the {@code "latitude_of_center"}
41 * parameter crosses the central line as the projection's central point.
42 * The {@linkplain #centralMeridian central meridian} is not a projection parameter,
43 * and is instead calculated as the intersection between the central line and the
44 * equator of the aposphere.
45 * <p>
46 * For the azimuth method, the central latitude cannot be &plusmn;90.0 degrees
47 * and the central line cannot be at a maximum or minimum latitude at the central point.
48 * In the two point method, the latitude of the first and second points cannot be
49 * equal. Also, the latitude of the first point and central point cannot be
50 * &plusmn;90.0 degrees. Furthermore, the latitude of the first point cannot be 0.0 and
51 * the latitude of the second point cannot be -90.0 degrees. A change of
52 * 10<sup>-7</sup> radians can allow calculation at these special cases. Snyder's restriction
53 * of the central latitude being 0.0 has been removed, since the equations appear
54 * to work correctly in this case.
55 * <p>
56 * Azimuth values of 0.0 and &plusmn;90.0 degrees are allowed (and used in Hungary
57 * and Switzerland), though these cases would usually use a Mercator or
58 * Transverse Mercator projection instead. Azimuth values &gt; 90 degrees cause
59 * errors in the equations.
60 * <p>
61 * The oblique mercator is also called the "Rectified Skew Orthomorphic" (RSO). It appears
62 * is that the only difference from the oblique mercator is that the RSO allows the rotation
63 * from the (<var>U</var>,<var>V</var>) to (<var>X</var>,<var>Y</var>) coordinate system to
64 * be different from the azimuth. This separate parameter is called
65 * {@code "rectified_grid_angle"} (or {@code "XY_Plane_Rotation"} by ESRI) and is also
66 * included in the EPSG's parameters for the Oblique Mercator and Hotine Oblique Mercator.
67 * The rotation parameter is optional in all the non-two point projections and will be
68 * set to the azimuth if not specified.
69 * <p>
70 * Projection cases and aliases implemented by the {@link ObliqueMercator} are:
71 * <ul>
72 *   <li>{@code Oblique_Mercator} (EPSG code 9815)<br>
73 *       grid coordinates begin at the central point,
74 *       has {@code "rectified_grid_angle"} parameter.</li>
75 *   <li>{@code Hotine_Oblique_Mercator_Azimuth_Center} (ESRI)<br>
76 *       grid coordinates begin at the central point.</li>
77 *   <li>{@code Rectified_Skew_Orthomorphic_Center} (ESRI)<br>
78 *       grid coordinates begin at the central point,
79 *       has {@code "rectified_grid_angle"} parameter.</li>
80 *
81 *   <li>{@code Hotine_Oblique_Mercator} (EPSG code 9812)<br>
82 *       grid coordinates begin at the interseciton of the central line and aposphere equator,
83 *       has {@code "rectified_grid_angle"} parameter.</li>
84 *   <li>{@code Hotine_Oblique_Mercator_Azimuth_Natural_Origin} (ESRI)<br>
85 *       grid coordinates begin at the interseciton of the central line and aposphere equator.</li>
86 *   <li>{@code Rectified_Skew_Orthomorphic_Natural_Origin} (ESRI)<br>
87 *       grid coordinates begin at the interseciton of the central line and aposphere equator,
88 *       has {@code "rectified_grid_angle"} parameter.</li>
89 *
90 *   <li>{@code Hotine_Oblique_Mercator_Two_Point_Center} (ESRI)<br>
91 *       grid coordinates begin at the central point.</li>
92 *   <li>{@code Hotine_Oblique_Mercator_Two_Point_Natural_Origin} (ESRI)<br>
93 *       grid coordinates begin at the interseciton of the central line and aposphere equator.</li>
94 * </ul>
95 * <p>
96 * This class has been derived from the implementation of the Geotools project;
97 * git 8cbf52d, org.geotools.referencing.operation.projection.ObliqueMercator
98 * at the time of migration.
99 * <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 throught 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}&times;{@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}
Note: See TracBrowser for help on using the repository browser.