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 * Albers Equal Area Projection (EPSG code 9822). This is a conic projection with parallels being
012 * unequally spaced arcs of concentric circles, more closely spaced at north and south edges of the
013 * map. Merideans are equally spaced radii of the same circles and intersect parallels at right
014 * angles. As the name implies, this projection minimizes distortion in areas.
015 * <p>
016 * The {@code "standard_parallel_2"} parameter is optional and will be given the same value as
017 * {@code "standard_parallel_1"} if not set (creating a 1 standard parallel projection).
018 * <p>
019 * <b>NOTE:</b>
020 * formulae used below are from a port, to Java, of the {@code proj4}
021 * package of the USGS survey. USGS work is acknowledged here.
022 * <p>
023 * This class has been derived from the implementation of the Geotools project;
024 * git 8cbf52d, org.geotools.referencing.operation.projection.AlbersEqualArea
025 * at the time of migration.
026 * <p>
027 * <b>References:</b>
028 * <ul>
029 *   <li> Proj-4.4.7 available at <A HREF="http://www.remotesensing.org/proj">www.remotesensing.org/proj</A><br>
030 *        Relevent files are: PJ_aea.c, pj_fwd.c and pj_inv.c </li>
031 *   <li> John P. Snyder (Map Projections - A Working Manual,
032 *        U.S. Geological Survey Professional Paper 1395, 1987)</li>
033 *   <li> "Coordinate Conversions and Transformations including Formulas",
034 *        EPSG Guidence Note Number 7, Version 19.</li>
035 * </ul>
036 *
037 * @author Gerald I. Evenden (for original code in Proj4)
038 * @author Rueben Schulz
039 *
040 * @see <A HREF="http://mathworld.wolfram.com/AlbersEqual-AreaConicProjection.html">Albers Equal-Area Conic Projection on MathWorld</A>
041 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/albers_equal_area_conic.html">"Albers_Conic_Equal_Area" on RemoteSensing.org</A>
042 * @see <A HREF="http://srmwww.gov.bc.ca/gis/bceprojection.html">British Columbia Albers Standard Projection</A>
043 *
044 * @since 9419
045 */
046public class AlbersEqualArea extends AbstractProj {
047
048    /**
049     * Maximum number of iterations for iterative computations.
050     */
051    private static final int MAXIMUM_ITERATIONS = 15;
052
053    /**
054     * Difference allowed in iterative computations.
055     */
056    private static final double ITERATION_TOLERANCE = 1E-10;
057
058    /**
059     * Maximum difference allowed when comparing real numbers.
060     */
061    private static final double EPSILON = 1E-6;
062
063    /**
064     * Constants used by the spherical and elliptical Albers projection.
065     */
066    private double n, n2, c, rho0;
067
068    /**
069     * An error condition indicating iteration will not converge for the
070     * inverse ellipse. See Snyder (14-20)
071     */
072    private double ec;
073
074    @Override
075    public String getName() {
076        return tr("Albers Equal Area");
077    }
078
079    @Override
080    public String getProj4Id() {
081        return "aea";
082    }
083
084    @Override
085    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
086        super.initialize(params);
087        if (params.lat0 == null)
088            throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
089        if (params.lat1 == null)
090            throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_1"));
091
092        double lat0 = Utils.toRadians(params.lat0);
093        // Standards parallels in radians.
094        double phi1 = Utils.toRadians(params.lat1);
095        double phi2 = params.lat2 == null ? phi1 : Utils.toRadians(params.lat2);
096
097        // Compute Constants
098        if (Math.abs(phi1 + phi2) < EPSILON) {
099            throw new ProjectionConfigurationException(tr("standard parallels are opposite"));
100        }
101        double sinphi = Math.sin(phi1);
102        double cosphi = Math.cos(phi1);
103        double n = sinphi;
104        boolean secant = Math.abs(phi1 - phi2) >= EPSILON;
105        double m1 = msfn(sinphi, cosphi);
106        double q1 = qsfn(sinphi);
107        if (secant) { // secant cone
108            sinphi = Math.sin(phi2);
109            cosphi = Math.cos(phi2);
110            double m2 = msfn(sinphi, cosphi);
111            double q2 = qsfn(sinphi);
112            n = (m1 * m1 - m2 * m2) / (q2 - q1);
113        }
114        c = m1 * m1 + n * q1;
115        rho0 = Math.sqrt(c - n * qsfn(Math.sin(lat0))) /n;
116        ec = 1.0 - .5 * (1.0-e2) *
117             Math.log((1.0 - e) / (1.0 + e)) / e;
118        n2 = n * n;
119        this.n = n;
120    }
121
122    @Override
123    public double[] project(double y, double x) {
124        double theta = n * x;
125        double rho = c - (spherical ? n2 * Math.sin(y) : n * qsfn(Math.sin(y)));
126        if (rho < 0.0) {
127            if (rho > -EPSILON) {
128                rho = 0.0;
129            } else {
130                throw new AssertionError();
131            }
132        }
133        rho = Math.sqrt(rho) / n;
134        // CHECKSTYLE.OFF: SingleSpaceSeparator
135        y = rho0 - rho * Math.cos(theta);
136        x =        rho * Math.sin(theta);
137        // CHECKSTYLE.ON: SingleSpaceSeparator
138        return new double[] {x, y};
139    }
140
141    @Override
142    public double[] invproject(double x, double y) {
143        y = rho0 - y;
144        double rho = Math.hypot(x, y);
145        if (rho > EPSILON) {
146            if (n < 0.0) {
147                rho = -rho;
148                x = -x;
149                y = -y;
150            }
151            x = Math.atan2(x, y) / n;
152            y = rho * n;
153            if (spherical) {
154                y = aasin((c - y*y) / n2);
155            } else {
156                y = (c - y*y) / n;
157                if (Math.abs(y) <= ec) {
158                    y = phi1(y);
159                } else {
160                    y = (y < 0.0) ? -Math.PI/2.0 : Math.PI/2.0;
161                }
162            }
163        } else {
164            x = 0.0;
165            y = n > 0.0 ? Math.PI/2.0 : -Math.PI/2.0;
166        }
167        return new double[] {y, x};
168    }
169
170    /**
171     * Iteratively solves equation (3-16) from Snyder.
172     *
173     * @param qs arcsin(q/2), used in the first step of iteration
174     * @return the latitude
175     */
176    public double phi1(final double qs) {
177        final double toneEs = 1 - e2;
178        double phi = Math.asin(0.5 * qs);
179        if (e < EPSILON) {
180            return phi;
181        }
182        for (int i = 0; i < MAXIMUM_ITERATIONS; i++) {
183            final double sinpi = Math.sin(phi);
184            final double cospi = Math.cos(phi);
185            final double con = e * sinpi;
186            final double com = 1.0 - con*con;
187            final double dphi = 0.5 * com*com / cospi *
188                    (qs/toneEs - sinpi / com + 0.5/e * Math.log((1. - con) / (1. + con)));
189            phi += dphi;
190            if (Math.abs(dphi) <= ITERATION_TOLERANCE) {
191                return phi;
192            }
193        }
194        throw new IllegalStateException("no convergence for qs="+qs);
195    }
196
197    /**
198     * Calculates q, Snyder equation (3-12)
199     *
200     * @param sinphi sin of the latitude q is calculated for
201     * @return q from Snyder equation (3-12)
202     */
203    private double qsfn(final double sinphi) {
204        final double oneEs = 1 - e2;
205        if (e >= EPSILON) {
206            final double con = e * sinphi;
207            return oneEs * (sinphi / (1. - con*con) -
208                    (0.5/e) * Math.log((1.-con) / (1.+con)));
209        } else {
210            return sinphi + sinphi;
211        }
212    }
213
214    @Override
215    public Bounds getAlgorithmBounds() {
216        return new Bounds(-89, -180, 89, 180, false);
217    }
218}