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 * Lambert Azimuthal Equal Area (EPSG code 9820).
012 * <p>
013 * This class has been derived from the implementation of the Geotools project;
014 * git 8cbf52d, org.geotools.referencing.operation.projection.LambertAzimuthalEqualArea
015 * at the time of migration.
016 * <p>
017 * <b>References:</b>
018 * <ul>
019 *   <li> A. Annoni, C. Luzet, E.Gubler and J. Ihde - Map Projections for Europe</li>
020 *   <li> John P. Snyder (Map Projections - A Working Manual,
021 *        U.S. Geological Survey Professional Paper 1395)</li>
022 * </ul>
023 *
024 * @author Gerald Evenden  (for original code in Proj4)
025 * @author Beate Stollberg
026 * @author Martin Desruisseaux
027 *
028 * @see <A HREF="http://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html">Lambert Azimuthal Equal-Area Projection</A>
029 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_azimuthal_equal_area.html">"Lambert_Azimuthal_Equal_Area"</A>
030 */
031public class LambertAzimuthalEqualArea extends AbstractProj {
032
033    /** Maximum difference allowed when comparing real numbers. */
034    private static final double EPSILON = 1E-7;
035
036    /** Epsilon for the comparison of small quantities. */
037    private static final double FINE_EPSILON = 1E-10;
038
039    /** Epsilon for the comparison of latitudes. */
040    private static final double EPSILON_LATITUDE = 1E-10;
041
042    /** Constants for authalic latitude. */
043    private static final double P00 = 3.33333333333333E-01;
044    private static final double P01 = 1.72222222222222E-01;
045    private static final double P02 = 1.02579365079365E-01;
046    private static final double P10 = 6.38888888888889E-02;
047    private static final double P11 = 6.64021164021164E-02;
048    private static final double P20 = 1.64150129421915E-02;
049
050    /** The projection mode. */
051    private enum Mode { OBLIQUE, EQUATORIAL, NORTH_POLE, SOUTH_POLE }
052
053    /** The projection mode for this particular instance. */
054    private Mode mode;
055
056    /** Constant parameters. */
057    private double sinb1, cosb1, xmf, ymf, qp, dd, rq;
058
059    /** Coefficients for authalic latitude. */
060    private double aPA0, aPA1, aPA2;
061
062    private double latitudeOfOrigin;
063
064    @Override
065    public String getName() {
066        return tr("Lambert Azimuthal Equal Area");
067    }
068
069    @Override
070    public String getProj4Id() {
071        return "laea";
072    }
073
074    @Override
075    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
076        super.initialize(params);
077
078        if (params.lat0 == null)
079            throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
080
081        latitudeOfOrigin = Utils.toRadians(params.lat0);
082        /*
083         * Detects the mode (oblique, etc.).
084         */
085        final double t = Math.abs(latitudeOfOrigin);
086        if (Math.abs(t - Math.PI/2) < EPSILON_LATITUDE) {
087            mode = latitudeOfOrigin < 0.0 ? Mode.SOUTH_POLE : Mode.NORTH_POLE;
088        } else if (Math.abs(t) < EPSILON_LATITUDE) {
089            mode = Mode.EQUATORIAL;
090        } else {
091            mode = Mode.OBLIQUE;
092        }
093        /*
094         * Computes the constants for authalic latitude.
095         */
096        final double es2 = e2 * e2;
097        final double es3 = e2 * es2;
098        aPA0 = P02 * es3 + P01 * es2 + P00 * e2;
099        aPA1 = P11 * es3 + P10 * es2;
100        aPA2 = P20 * es3;
101
102        final double sinphi;
103        qp = qsfn(1);
104        rq = Math.sqrt(0.5 * qp);
105        sinphi = Math.sin(latitudeOfOrigin);
106        sinb1 = qsfn(sinphi) / qp;
107        cosb1 = Math.sqrt(1.0 - sinb1 * sinb1);
108        switch (mode) {
109            case NORTH_POLE:  // Fall through
110            case SOUTH_POLE:
111                dd = 1.0;
112                xmf = ymf = rq;
113                break;
114            case EQUATORIAL:
115                dd = 1.0 / rq;
116                xmf = 1.0;
117                ymf = 0.5 * qp;
118                break;
119            case OBLIQUE:
120                dd = Math.cos(latitudeOfOrigin) / (Math.sqrt(1.0 - e2 * sinphi * sinphi) * rq * cosb1);
121                xmf = rq * dd;
122                ymf = rq / dd;
123                break;
124            default:
125                throw new AssertionError(mode);
126        }
127    }
128
129    @Override
130    public double[] project(final double phi, final double lambda) {
131        final double coslam = Math.cos(lambda);
132        final double sinlam = Math.sin(lambda);
133        final double sinphi = Math.sin(phi);
134        double q = qsfn(sinphi);
135        final double sinb, cosb, b, c, x, y;
136        switch (mode) {
137            case OBLIQUE:
138                sinb = q / qp;
139                cosb = Math.sqrt(1.0 - sinb * sinb);
140                c = 1.0 + sinb1 * sinb + cosb1 * cosb * coslam;
141                b = Math.sqrt(2.0 / c);
142                y = ymf * b * (cosb1 * sinb - sinb1 * cosb * coslam);
143                x = xmf * b * cosb * sinlam;
144                break;
145            case EQUATORIAL:
146                sinb = q / qp;
147                cosb = Math.sqrt(1.0 - sinb * sinb);
148                c = 1.0 + cosb * coslam;
149                b = Math.sqrt(2.0 / c);
150                y = ymf * b * sinb;
151                x = xmf * b * cosb * sinlam;
152                break;
153            case NORTH_POLE:
154                c = (Math.PI / 2) + phi;
155                q = qp - q;
156                if (q >= 0.0) {
157                    b = Math.sqrt(q);
158                    x = b * sinlam;
159                    y = coslam * -b;
160                } else {
161                    x = y = 0.;
162                }
163                break;
164            case SOUTH_POLE:
165                c = phi - (Math.PI / 2);
166                q = qp + q;
167                if (q >= 0.0) {
168                    b = Math.sqrt(q);
169                    x = b * sinlam;
170                    y = coslam * +b;
171                } else {
172                    x = y = 0.;
173                }
174                break;
175            default:
176                throw new AssertionError(mode);
177        }
178        if (Math.abs(c) < EPSILON_LATITUDE) {
179            return new double[] {0, 0}; // this is an error, we should handle it somehow
180        }
181        return new double[] {x, y};
182    }
183
184    @Override
185    public double[] invproject(double x, double y) {
186        switch (mode) {
187            case EQUATORIAL: // Fall through
188            case OBLIQUE:
189                return invprojectEO(x, y);
190            case NORTH_POLE:
191                return invprojectNS(x, -y);
192            case SOUTH_POLE:
193                return invprojectNS(x, y);
194            default:
195                throw new AssertionError(mode);
196        }
197    }
198
199    private double[] invprojectEO(double x, double y) {
200        final double lambda;
201        final double phi;
202        x /= dd;
203        y *= dd;
204        final double rho = Math.hypot(x, y);
205        if (rho < FINE_EPSILON) {
206            lambda = 0.0;
207            phi = latitudeOfOrigin;
208        } else {
209            final double ab;
210            double sCe = 2.0 * Math.asin(0.5 * rho / rq);
211            double cCe = Math.cos(sCe);
212            sCe = Math.sin(sCe);
213            x *= sCe;
214            if (mode == Mode.OBLIQUE) {
215                ab = cCe * sinb1 + y * sCe * cosb1 / rho;
216                y = rho * cosb1 * cCe - y * sinb1 * sCe;
217            } else {
218                ab = y * sCe / rho;
219                y = rho * cCe;
220            }
221            lambda = Math.atan2(x, y);
222            phi = authlat(Math.asin(ab));
223        }
224        return new double[] {phi, lambda};
225    }
226
227    private double[] invprojectNS(double x, double y) {
228        final double lambda;
229        final double phi;
230        final double q = x*x + y*y;
231        if (q == 0) {
232            lambda = 0.;
233            phi = latitudeOfOrigin;
234        } else {
235            double ab = 1.0 - q / qp;
236            if (mode == Mode.SOUTH_POLE) {
237                ab = -ab;
238            }
239            lambda = Math.atan2(x, y);
240            phi = authlat(Math.asin(ab));
241        }
242        return new double[] {phi, lambda};
243    }
244
245    /**
246     * Calculates <var>q</var>, Snyder equation (3-12)
247     *
248     * @param sinphi sin of the latitude <var>q</var> is calculated for.
249     * @return <var>q</var> from Snyder equation (3-12).
250     */
251    private double qsfn(final double sinphi) {
252        if (e >= EPSILON) {
253            final double con = e * sinphi;
254            return (1.0 - e2) * (sinphi / (1.0 - con*con) -
255                    (0.5 / e) * Math.log((1.0 - con) / (1.0 + con)));
256        } else {
257            return sinphi + sinphi;
258        }
259    }
260
261    /**
262     * Determines latitude from authalic latitude.
263     * @param beta authalic latitude
264     * @return corresponding latitude
265     */
266    private double authlat(final double beta) {
267        final double t = beta + beta;
268        return beta + aPA0 * Math.sin(t) + aPA1 * Math.sin(t+t) + aPA2 * Math.sin(t+t+t);
269    }
270
271    @Override
272    public Bounds getAlgorithmBounds() {
273        return new Bounds(-89, -174, 89, 174, false);
274    }
275}