001// License: GPL. For details, see LICENSE file.
002package org.openstreetmap.josm.data.projection.proj;
003
004import static java.lang.Math.abs;
005import static java.lang.Math.atan;
006import static java.lang.Math.cos;
007import static java.lang.Math.exp;
008import static java.lang.Math.log;
009import static java.lang.Math.pow;
010import static java.lang.Math.sin;
011import static java.lang.Math.sqrt;
012import static org.openstreetmap.josm.tools.I18n.tr;
013import static org.openstreetmap.josm.tools.Utils.toRadians;
014
015import org.openstreetmap.josm.data.Bounds;
016import org.openstreetmap.josm.data.projection.CustomProjection.Param;
017import org.openstreetmap.josm.data.projection.Ellipsoid;
018import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
019
020/**
021 * Lambert Conical Conformal Projection. Areas and shapes are deformed as one moves away from standard parallels.
022 * The angles are true in a limited area. This projection is used for the charts of North America, France and Belgium.
023 * <p>
024 * This implementation provides transforms for two cases of the lambert conic conformal projection:
025 * </p>
026 * <ul>
027 *   <li>{@code Lambert_Conformal_Conic_1SP} (EPSG code 9801)</li>
028 *   <li>{@code Lambert_Conformal_Conic_2SP} (EPSG code 9802)</li>
029 * </ul>
030 * <p>
031 * For the 1SP case the latitude of origin is used as the standard parallel (SP).
032 * To use 1SP with a latitude of origin different from the SP, use the 2SP and set the SP1 to the single SP.
033 * The {@code "standard_parallel_2"} parameter is optional and will be given the same value
034 * as {@code "standard_parallel_1"} if not set (creating a 1 standard parallel projection).
035 * </p>
036 * <b>References:</b>
037 * <ul>
038 *   <li>John P. Snyder (Map Projections - A Working Manual,<br>U.S. Geological Survey Professional Paper 1395, 1987)</li>
039 *   <li>"Coordinate Conversions and Transformations including Formulas",<br>EPSG Guidence Note Number 7, Version 19.</li>
040 * </ul>
041 *
042 * @author Pieren
043 * @author André Gosselin
044 * @author Martin Desruisseaux (PMO, IRD)
045 * @author Rueben Schulz
046 *
047 * @see <A HREF="http://mathworld.wolfram.com/LambertConformalConicProjection.html">Lambert conformal conic projection on MathWorld</A>
048 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html">lambert_conic_conformal_1sp</A>
049 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_2sp.html">lambert_conic_conformal_2sp</A>
050 *
051 * @since 13639 (align implementation with proj.4 / GeoTools)
052 * @since 4285 (reworked from Lambert / LambertCC9Zones)
053 * @since 2304 (initial implementation by Pieren)
054 */
055public class LambertConformalConic extends AbstractProj {
056
057    /** ellipsoid */
058    protected Ellipsoid ellps;
059
060    /**
061     * Base class of Lambert Conformal Conic parameters.
062     */
063    public static class Parameters {
064        /** latitude of origin */
065        public final double latitudeOrigin;
066
067        /**
068         * Constructs a new {@code Parameters}.
069         * @param latitudeOrigin latitude of origin
070         */
071        protected Parameters(double latitudeOrigin) {
072            this.latitudeOrigin = latitudeOrigin;
073        }
074    }
075
076    /**
077     * Parameters with a single standard parallel.
078     */
079    public static class Parameters1SP extends Parameters {
080        /**
081         * Constructs a new {@code Parameters1SP}.
082         * @param latitudeOrigin latitude of origin
083         */
084        public Parameters1SP(double latitudeOrigin) {
085            super(latitudeOrigin);
086        }
087    }
088
089    /**
090     * Parameters with two standard parallels.
091     */
092    public static class Parameters2SP extends Parameters {
093        /** first standard parallel */
094        public final double standardParallel1;
095        /** second standard parallel */
096        public final double standardParallel2;
097
098        /**
099         * Constructs a new {@code Parameters2SP}.
100         * @param latitudeOrigin latitude of origin
101         * @param standardParallel1 first standard parallel
102         * @param standardParallel2 second standard parallel
103         */
104        public Parameters2SP(double latitudeOrigin, double standardParallel1, double standardParallel2) {
105            super(latitudeOrigin);
106            this.standardParallel1 = standardParallel1;
107            this.standardParallel2 = standardParallel2;
108        }
109    }
110
111    private Parameters params;
112
113    /**
114     * projection exponent
115     */
116    protected double n;
117
118    /**
119     * projection factor
120     */
121    protected double f;
122
123    /**
124     * radius of the parallel of latitude of the false origin (2SP) or at natural origin (1SP)
125     */
126    protected double r0;
127
128    /**
129     * precision in iterative schema
130     */
131    protected static final double epsilon = 1e-12;
132
133    @Override
134    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
135        super.initialize(params);
136        ellps = params.ellps;
137        if (params.lat0 == null)
138            throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", Param.lat_0.key));
139        if (params.lat1 != null && params.lat2 != null) {
140            this.params = new Parameters2SP(params.lat0, params.lat1, params.lat2);
141            initialize2SP(toRadians(params.lat0), toRadians(params.lat1), toRadians(params.lat2));
142        } else {
143            this.params = new Parameters1SP(params.lat0);
144            initialize1SP(toRadians(params.lat0));
145        }
146    }
147
148    /**
149     * Initialize for LCC with 2 standard parallels.
150     *
151     * @param lat0 latitude of false origin (in radians)
152     * @param lat1 latitude of first standard parallel (in radians)
153     * @param lat2 latitude of second standard parallel (in radians)
154     */
155    private void initialize2SP(double lat0, double lat1, double lat2) {
156
157        final double cosphi1 = cos(lat1);
158        final double sinphi1 = sin(lat1);
159
160        final double cosphi2 = cos(lat2);
161        final double sinphi2 = sin(lat2);
162
163        final double m1 = msfn(sinphi1, cosphi1);
164        final double m2 = msfn(sinphi2, cosphi2);
165
166        final double t0 = tsfn(lat0, sin(lat0));
167        final double t1 = tsfn(lat1, sinphi1);
168        final double t2 = tsfn(lat2, sinphi2);
169
170        n = log(m1/m2) / log(t1/t2);
171        f = m1 * pow(t1, -n) / n;
172        r0 = f * pow(t0, n);
173    }
174
175    /**
176     * Initialize for LCC with 1 standard parallel.
177     *
178     * @param lat0 latitude of natural origin (in radians)
179     */
180    private void initialize1SP(double lat0) {
181        final double m0 = msfn(sin(lat0), cos(lat0));
182        final double t0 = tsfn(lat0, sin(lat0));
183
184        n = sin(lat0);
185        f = m0 * pow(t0, -n) / n;
186        r0 = f * pow(t0, n);
187    }
188
189    @Override
190    public String getName() {
191        return tr("Lambert Conformal Conic");
192    }
193
194    @Override
195    public String getProj4Id() {
196        return "lcc";
197    }
198
199    @Override
200    public double[] project(double phi, double lambda) {
201        double sinphi = sin(phi);
202        double l = (0.5*log((1+sinphi)/(1-sinphi))) - e/2*log((1+e*sinphi)/(1-e*sinphi));
203        double r = f*exp(-n*l);
204        double gamma = n*lambda;
205        double x = r*sin(gamma);
206        double y = r0 - r*cos(gamma);
207        return new double[] {x, y};
208    }
209
210    @Override
211    public double[] invproject(double east, double north) {
212        double r = sqrt(pow(east, 2) + pow(north-r0, 2));
213        double gamma = atan(east / (r0-north));
214        double lambda = gamma/n;
215        double latIso = (-1/n) * log(abs(r/f));
216        double phi = ellps.latitude(latIso, e, epsilon);
217        return new double[] {phi, lambda};
218    }
219
220    /**
221     * Returns projection parameters.
222     * @return projection parameters
223     */
224    public final Parameters getParameters() {
225        return params;
226    }
227
228    @Override
229    public Bounds getAlgorithmBounds() {
230        double lat;
231        if (params instanceof Parameters2SP) {
232            Parameters2SP p2p = (Parameters2SP) params;
233            lat = (p2p.standardParallel1 + p2p.standardParallel2) / 2;
234        } else {
235            lat = params.latitudeOrigin;
236        }
237        double minlat = Math.max(lat - 60, -89);
238        double maxlat = Math.min(lat + 60, 89);
239        return new Bounds(minlat, -85, maxlat, 85, false);
240    }
241}