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.JosmRuntimeException;
009
010/**
011 * Azimuthal Equidistant projection.
012 * <p>
013 * This implementation does not include the Guam or Micronesia variants.
014 *
015 * @author Gerald Evenden (original PROJ.4 implementation in C)
016 * @author Ben Caradoc-Davies (Transient Software Limited)
017 * @see <a href="https://pubs.er.usgs.gov/publication/pp1395"><em>Map Projections: A Working Manual</em>, Snyder (1987), pages 191-202</a>
018 * @see <a href="http://geotiff.maptools.org/proj_list/azimuthal_equidistant.html">PROJ.4 notes on parameters</a>
019 * @see <a href="https://github.com/OSGeo/proj.4/blob/master/src/PJ_aeqd.c">PROJ.4 implemention in C</a>
020 * @see <a href="https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection">Wikipedia</a>
021 * @see <a href="http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html">Wolfram Alpha</a>
022 * @since 13598
023 */
024public class AzimuthalEquidistant extends AbstractProj {
025
026    /**
027     * Less strict tolerance.
028     */
029    public static final double EPS10 = 1.e-10;
030
031    /**
032     * Stricter tolerance.
033     */
034    public static final double TOL = 1.e-14;
035
036    /**
037     * Half of π.
038     */
039    public static final double HALF_PI = Math.PI / 2;
040
041    /**
042     * The four possible modes or aspects of the projection.
043     */
044    public enum Mode {
045        /** North pole */
046        NORTH_POLAR,
047        /** South pole */
048        SOUTH_POLAR,
049        /** Equator */
050        EQUATORIAL,
051        /** Oblique */
052        OBLIQUE;
053    }
054
055    /**
056     * Length of semi-major axis, in metres. This is named '<var>a</var>' or '<var>R</var>'
057     * (Radius in spherical cases) in Snyder.
058     */
059    protected double semiMajor;
060
061    /**
062     * Length of semi-minor axis, in metres. This is named '<var>b</var>' in Snyder.
063     */
064    protected double semiMinor;
065
066    /**
067     * Central longitude in <u>radians</u>. Default value is 0, the Greenwich meridian.
068     * This is called '<var>lambda0</var>' in Snyder.
069     */
070    protected double centralMeridian;
071
072    /**
073     * Latitude of origin in <u>radians</u>. Default value is 0, the equator.
074     * This is called '<var>phi0</var>' in Snyder.
075     */
076    protected double latitudeOfOrigin;
077
078    /**
079     * The mode or aspect of the projection.
080     */
081    protected Mode mode;
082
083    /**
084     * Geodesic calculator used for this projection. Not used and set to null for polar projections.
085     */
086    //protected Geodesic geodesic; // See https://josm.openstreetmap.de/ticket/16129#comment:21
087
088    /**
089     * The sine of the central latitude of the projection.
090     */
091    protected double sinph0;
092
093    /**
094     * The cosine of the central latitude of the projection.
095     */
096    protected double cosph0;
097
098    /**
099     * Meridian distance from the equator to the pole. Not used and set to NaN for non-polar projections.
100     */
101    protected double mp;
102
103    @Override
104    public String getName() {
105        return tr("Azimuthal Equidistant");
106    }
107
108    @Override
109    public String getProj4Id() {
110        return "aeqd";
111    }
112
113    @Override
114    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
115        super.initialize(params);
116        if (params.lon0 == null)
117            throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lon_0"));
118        if (params.lat0 == null)
119            throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
120        centralMeridian = Math.toRadians(params.lon0);
121        latitudeOfOrigin = Math.toRadians(params.lat0);
122        semiMajor = params.ellps.a;
123        semiMinor = params.ellps.b;
124        if (Math.abs(latitudeOfOrigin - HALF_PI) < EPS10) {
125            mode = Mode.NORTH_POLAR;
126            mp = mlfn(HALF_PI, 1, 0);
127            sinph0 = 1;
128            cosph0 = 0;
129        } else if (Math.abs(latitudeOfOrigin + HALF_PI) < EPS10) {
130            mode = Mode.SOUTH_POLAR;
131            mp = mlfn(-HALF_PI, -1, 0);
132            sinph0 = -1;
133            cosph0 = 0;
134        } else if (Math.abs(latitudeOfOrigin) < EPS10) {
135            mode = Mode.EQUATORIAL;
136            mp = Double.NaN;
137            sinph0 = 0;
138            cosph0 = 1;
139            //geodesic = new Geodesic(semiMajor, (semiMajor - semiMinor) / semiMajor);
140            throw new ProjectionConfigurationException("Equatorial AzimuthalEquidistant not yet supported");
141        } else {
142            mode = Mode.OBLIQUE;
143            mp = Double.NaN;
144            sinph0 = Math.sin(latitudeOfOrigin);
145            cosph0 = Math.cos(latitudeOfOrigin);
146            //geodesic = new Geodesic(semiMajor, (semiMajor - semiMinor) / semiMajor);
147            throw new ProjectionConfigurationException("Oblique AzimuthalEquidistant not yet supported");
148        }
149    }
150
151    @Override
152    public Bounds getAlgorithmBounds() {
153        return new Bounds(-89, -180, 89, 180, false);
154    }
155
156    @Override
157    public double[] project(double latRad, double lonRad) {
158        return spherical ? projectSpherical(latRad, lonRad) : projectEllipsoidal(latRad, lonRad);
159    }
160
161    @Override
162    public double[] invproject(double east, double north) {
163        return spherical ? invprojectSpherical(east, north) : invprojectEllipsoidal(east, north);
164    }
165
166    double[] projectSpherical(double latRad, double lonRad) {
167        double x = 0;
168        double y = 0;
169        double sinphi = Math.sin(latRad);
170        double cosphi = Math.cos(latRad);
171        double coslam = Math.cos(lonRad);
172        switch (mode) {
173        case EQUATORIAL:
174        case OBLIQUE:
175            if (mode == Mode.EQUATORIAL) {
176                y = cosphi * coslam;
177            } else { // Oblique
178                y = sinph0 * sinphi + cosph0 * cosphi * coslam;
179            }
180            if (Math.abs(Math.abs(y) - 1) < TOL) {
181                if (y < 0) {
182                    throw new JosmRuntimeException("TOLERANCE_ERROR");
183                } else {
184                    x = 0;
185                    y = 0;
186                }
187            } else {
188                y = Math.acos(y);
189                y /= Math.sin(y);
190                x = y * cosphi * Math.sin(lonRad);
191                y *= (mode == Mode.EQUATORIAL) ? sinphi
192                        : (cosph0 * sinphi - sinph0 * cosphi * coslam);
193            }
194            break;
195        case NORTH_POLAR:
196            latRad = -latRad;
197            coslam = -coslam;
198            // fall through
199        case SOUTH_POLAR:
200            if (Math.abs(latRad - HALF_PI) < EPS10) {
201                throw new JosmRuntimeException("TOLERANCE_ERROR");
202            }
203            y = HALF_PI + latRad;
204            x = y * Math.sin(lonRad);
205            y *= coslam;
206            break;
207        }
208        return new double[] {x, y};
209    }
210
211    double[] invprojectSpherical(double east, double north) {
212        double x = east;
213        double y = north;
214        double lambda = 0;
215        double phi = 0;
216        double c_rh = Math.hypot(x, y);
217        if (c_rh > Math.PI) {
218            if (c_rh - EPS10 > Math.PI) {
219                throw new JosmRuntimeException("TOLERANCE_ERROR");
220            }
221        } else if (c_rh < EPS10) {
222            phi = latitudeOfOrigin;
223            lambda = 0.;
224        } else {
225            if (mode == Mode.OBLIQUE || mode == Mode.EQUATORIAL) {
226                double sinc = Math.sin(c_rh);
227                double cosc = Math.cos(c_rh);
228                if (mode == Mode.EQUATORIAL) {
229                    phi = aasin(y * sinc / c_rh);
230                    x *= sinc;
231                    y = cosc * c_rh;
232                } else { // Oblique
233                    phi = aasin(cosc * sinph0 + y * sinc * cosph0 / c_rh);
234                    y = (cosc - sinph0 * Math.sin(phi)) * c_rh;
235                    x *= sinc * cosph0;
236                }
237                lambda = (y == 0) ? 0 : Math.atan2(x, y);
238            } else if (mode == Mode.NORTH_POLAR) {
239                phi = HALF_PI - c_rh;
240                lambda = Math.atan2(x, -y);
241            } else { // South Polar
242                phi = c_rh - HALF_PI;
243                lambda = Math.atan2(x, y);
244            }
245        }
246        return new double[] {phi, lambda};
247    }
248
249    double[] projectEllipsoidal(double latRad, double lonRad) {
250        double x = 0;
251        double y = 0;
252        double coslam = Math.cos(lonRad);
253        double cosphi = Math.cos(latRad);
254        double sinphi = Math.sin(latRad);
255        switch (mode) {
256        case NORTH_POLAR:
257            coslam = -coslam;
258            // fall through
259        case SOUTH_POLAR:
260            double rho = Math.abs(mp - mlfn(latRad, sinphi, cosphi));
261            x = rho * Math.sin(lonRad);
262            y = rho * coslam;
263            break;
264        case EQUATORIAL:
265        case OBLIQUE:
266            if (Math.abs(lonRad) < EPS10 && Math.abs(latRad - latitudeOfOrigin) < EPS10) {
267                x = 0;
268                y = 0;
269                break;
270            }
271            /*GeodesicData g = geodesic.Inverse(Math.toDegrees(latitudeOfOrigin),
272                    Math.toDegrees(centralMeridian), Math.toDegrees(latRad),
273                    Math.toDegrees(lonRad + centralMeridian));
274            double azi1 = Math.toRadians(g.azi1);
275            x = g.s12 * Math.sin(azi1) / semiMajor;
276            y = g.s12 * Math.cos(azi1) / semiMajor;*/
277            break;
278        }
279        return new double[] {x, y};
280    }
281
282    double[] invprojectEllipsoidal(double east, double north) {
283        double x = east;
284        double y = north;
285        double lambda = 0;
286        double phi = 0;
287        double c = Math.hypot(x, y);
288        if (c < EPS10) {
289            phi = latitudeOfOrigin;
290            lambda = 0;
291        } else {
292            if (mode == Mode.OBLIQUE || mode == Mode.EQUATORIAL) {
293                /*double x2 = x * semiMajor;
294                double y2 = y * semiMajor;
295                double azi1 = Math.atan2(x2, y2);
296                double s12 = Math.sqrt(x2 * x2 + y2 * y2);
297                GeodesicData g = geodesic.Direct(Math.toDegrees(latitudeOfOrigin),
298                        Math.toDegrees(centralMeridian), Math.toDegrees(azi1), s12);
299                phi = Math.toRadians(g.lat2);
300                lambda = Math.toRadians(g.lon2);*/
301                lambda -= centralMeridian;
302            } else { // Polar
303                phi = invMlfn((mode == Mode.NORTH_POLAR) ? (mp - c) : (mp + c));
304                lambda = Math.atan2(x, (mode == Mode.NORTH_POLAR) ? -y : y);
305            }
306        }
307        return new double[] {phi, lambda};
308    }
309}