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.projection.ProjectionConfigurationException;
008import org.openstreetmap.josm.tools.Utils;
009
010/**
011 * The polar case of the stereographic projection.
012 * <p>
013 * In the proj.4 library, the code "stere" covers several variants of the
014 * Stereographic projection, depending on the latitude of natural origin
015 * (parameter lat_0).
016 * <p>
017 *
018 * In this file, only the polar case is implemented. This corresponds to
019 * EPSG:9810 (Polar Stereographic Variant A) and EPSG:9829 (Polar Stereographic
020 * Variant B).
021 * <p>
022 *
023 * It is required, that the latitude of natural origin has the value +/-90 degrees.
024 * <p>
025 *
026 * This class has been derived from the implementation of the Geotools project;
027 * git 8cbf52d, org.geotools.referencing.operation.projection.PolarStereographic
028 * at the time of migration.
029 * <p>
030 *
031 * <b>References:</b>
032 * <ul>
033 *   <li>John P. Snyder (Map Projections - A Working Manual,<br>
034 *       U.S. Geological Survey Professional Paper 1395, 1987)</li>
035 *   <li>"Coordinate Conversions and Transformations including Formulas",<br>
036 *       EPSG Guidence Note Number 7, Version 19.</li>
037 *   <li>Gerald Evenden. <A HREF="http://members.bellatlantic.net/~vze2hc4d/proj4/sterea.pdf">
038 *       "Supplementary PROJ.4 Notes - Oblique Stereographic Alternative"</A></li>
039 *   <li>Krakiwsky, E.J., D.B. Thomson, and R.R. Steeves. 1977. A Manual
040 *       For Geodetic Coordinate Transformations in the Maritimes.
041 *       Geodesy and Geomatics Engineering, UNB. Technical Report No. 48.</li>
042 *   <li>Thomson, D.B., M.P. Mepham and R.R. Steeves. 1977.
043 *       The Stereographic Double Projection.
044 *       Geodesy and Geomatics Engineereng, UNB. Technical Report No. 46.</li>
045 * </ul>
046 *
047 * @author André Gosselin
048 * @author Martin Desruisseaux (PMO, IRD)
049 * @author Rueben Schulz
050 *
051 * @see <A HREF="http://mathworld.wolfram.com/StereographicProjection.html">Stereographic projection on MathWorld</A>
052 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/polar_stereographic.html">Polar_Stereographic</A>
053 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/oblique_stereographic.html">Oblique_Stereographic</A>
054 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/stereographic.html">Stereographic</A>
055 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/random_issues.html#stereographic">Some Random Stereographic Issues</A>
056 *
057 * @see DoubleStereographic
058 * @since 9419
059 */
060public class PolarStereographic extends AbstractProj {
061    /**
062     * Maximum number of iterations for iterative computations.
063     */
064    private static final int MAXIMUM_ITERATIONS = 15;
065
066    /**
067     * Difference allowed in iterative computations.
068     */
069    private static final double ITERATION_TOLERANCE = 1E-10;
070
071    /**
072     * Maximum difference allowed when comparing real numbers.
073     */
074    private static final double EPSILON = 1E-8;
075
076    /**
077     * A constant used in the transformations.
078     */
079    private double k0;
080
081    /**
082     * {@code true} if this projection is for the south pole, or {@code false}
083     * if it is for the north pole.
084     */
085    boolean southPole;
086
087    @Override
088    public String getName() {
089        return tr("Polar Stereographic");
090    }
091
092    @Override
093    public String getProj4Id() {
094        return "stere";
095    }
096
097    @Override
098    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
099        super.initialize(params);
100        if (params.lat0 == null)
101            throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
102        if (params.lat0 != 90.0 && params.lat0 != -90.0)
103            throw new ProjectionConfigurationException(
104                    tr("Polar Stereographic: Parameter ''{0}'' must be 90 or -90.", "lat_0"));
105        // Latitude of true scale, in radians
106        double latitudeTrueScale;
107        if (params.lat_ts == null) {
108            latitudeTrueScale = (params.lat0 < 0) ? -Math.PI/2 : Math.PI/2;
109        } else {
110            latitudeTrueScale = Utils.toRadians(params.lat_ts);
111        }
112        southPole = latitudeTrueScale < 0;
113
114        // Computes coefficients
115        double latitudeTrueScaleAbs = Math.abs(latitudeTrueScale);
116        if (Math.abs(latitudeTrueScaleAbs - Math.PI/2) >= EPSILON) {
117            final double t = Math.sin(latitudeTrueScaleAbs);
118            k0 = msfn(t, Math.cos(latitudeTrueScaleAbs)) /
119                 tsfn(latitudeTrueScaleAbs, t); // Derives from (21-32 and 21-33)
120        } else {
121            // True scale at pole (part of (21-33))
122            k0 = 2.0 / Math.sqrt(Math.pow(1+e, 1+e)*
123                                 Math.pow(1-e, 1-e));
124        }
125    }
126
127    @Override
128    public double[] project(double y, double x) {
129        final double sinlat = Math.sin(y);
130        final double coslon = Math.cos(x);
131        final double sinlon = Math.sin(x);
132        if (southPole) {
133            final double rho = k0 * tsfn(-y, -sinlat);
134            x = rho * sinlon;
135            y = rho * coslon;
136        } else {
137            final double rho = k0 * tsfn(y, sinlat);
138            x = rho * sinlon;
139            y = -rho * coslon;
140        }
141        return new double[] {x, y};
142    }
143
144    @Override
145    public double[] invproject(double x, double y) {
146        final double rho = Math.hypot(x, y);
147        if (southPole) {
148            y = -y;
149        }
150        /*
151         * Compute latitude using iterative technique.
152         */
153        final double t = rho/k0;
154        final double halfe = e/2.0;
155        double phi0 = 0;
156        for (int i = MAXIMUM_ITERATIONS;;) {
157            final double esinphi = e * Math.sin(phi0);
158            final double phi = (Math.PI/2) - 2.0*Math.atan(t*Math.pow((1-esinphi)/(1+esinphi), halfe));
159            if (Math.abs(phi-phi0) < ITERATION_TOLERANCE) {
160                x = (Math.abs(rho) < EPSILON) ? 0.0 : Math.atan2(x, -y);
161                y = southPole ? -phi : phi;
162                break;
163            }
164            phi0 = phi;
165            if (--i < 0) {
166                throw new IllegalStateException("no convergence for x="+x+", y="+y);
167            }
168        }
169        return new double[] {y, x};
170    }
171
172    @Override
173    public Bounds getAlgorithmBounds() {
174        final double cut = 60;
175        if (southPole) {
176            return new Bounds(-90, -180, cut, 180, false);
177        } else {
178            return new Bounds(-cut, -180, 90, 180, false);
179        }
180    }
181}