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