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