001// License: GPL. For details, see LICENSE file.
002package org.openstreetmap.josm.data.projection.proj;
003
004import static java.lang.Math.PI;
005import static java.lang.Math.abs;
006import static java.lang.Math.atan;
007import static java.lang.Math.cos;
008import static java.lang.Math.exp;
009import static java.lang.Math.log;
010import static java.lang.Math.pow;
011import static java.lang.Math.sin;
012import static java.lang.Math.sqrt;
013import static java.lang.Math.tan;
014import static org.openstreetmap.josm.tools.I18n.tr;
015
016import org.openstreetmap.josm.data.Bounds;
017import org.openstreetmap.josm.data.projection.CustomProjection.Param;
018import org.openstreetmap.josm.data.projection.Ellipsoid;
019import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
020import org.openstreetmap.josm.tools.Utils;
021
022/**
023 * Implementation of the Lambert Conformal Conic projection.
024 *
025 * @author Pieren
026 */
027public class LambertConformalConic extends AbstractProj {
028
029    /** ellipsoid */
030    protected Ellipsoid ellps;
031
032    /**
033     * Base class of Lambert Conformal Conic parameters.
034     */
035    public static class Parameters {
036        /** latitude of origin */
037        public final double latitudeOrigin;
038
039        /**
040         * Constructs a new {@code Parameters}.
041         * @param latitudeOrigin latitude of origin
042         */
043        protected Parameters(double latitudeOrigin) {
044            this.latitudeOrigin = latitudeOrigin;
045        }
046    }
047
048    /**
049     * Parameters with a single standard parallel.
050     */
051    public static class Parameters1SP extends Parameters {
052        /**
053         * Constructs a new {@code Parameters1SP}.
054         * @param latitudeOrigin latitude of origin
055         */
056        public Parameters1SP(double latitudeOrigin) {
057            super(latitudeOrigin);
058        }
059    }
060
061    /**
062     * Parameters with two standard parallels.
063     */
064    public static class Parameters2SP extends Parameters {
065        /** first standard parallel */
066        public final double standardParallel1;
067        /** second standard parallel */
068        public final double standardParallel2;
069
070        /**
071         * Constructs a new {@code Parameters2SP}.
072         * @param latitudeOrigin latitude of origin
073         * @param standardParallel1 first standard parallel
074         * @param standardParallel2 second standard parallel
075         */
076        public Parameters2SP(double latitudeOrigin, double standardParallel1, double standardParallel2) {
077            super(latitudeOrigin);
078            this.standardParallel1 = standardParallel1;
079            this.standardParallel2 = standardParallel2;
080        }
081    }
082
083    private Parameters params;
084
085    /**
086     * projection exponent
087     */
088    protected double n;
089    /**
090     * projection factor
091     */
092    protected double f;
093    /**
094     * radius of the parallel of latitude of the false origin (2SP) or at
095     * natural origin (1SP)
096     */
097    protected double r0;
098
099    /**
100     * precision in iterative schema
101     */
102    protected static final double epsilon = 1e-12;
103
104    @Override
105    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
106        super.initialize(params);
107        ellps = params.ellps;
108        if (params.lat0 == null)
109            throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", Param.lat_0.key));
110        if (params.lat1 != null && params.lat2 != null) {
111            initialize2SP(params.lat0, params.lat1, params.lat2);
112        } else {
113            initialize1SP(params.lat0);
114        }
115    }
116
117    /**
118     * Initialize for LCC with 2 standard parallels.
119     *
120     * @param lat0 latitude of false origin (in degrees)
121     * @param lat1 latitude of first standard parallel (in degrees)
122     * @param lat2 latitude of second standard parallel (in degrees)
123     */
124    private void initialize2SP(double lat0, double lat1, double lat2) {
125        this.params = new Parameters2SP(lat0, lat1, lat2);
126
127        final double m1 = m(Utils.toRadians(lat1));
128        final double m2 = m(Utils.toRadians(lat2));
129
130        final double t1 = t(Utils.toRadians(lat1));
131        final double t2 = t(Utils.toRadians(lat2));
132        final double tf = t(Utils.toRadians(lat0));
133
134        n = (log(m1) - log(m2)) / (log(t1) - log(t2));
135        f = m1 / (n * pow(t1, n));
136        r0 = f * pow(tf, n);
137    }
138
139    /**
140     * Initialize for LCC with 1 standard parallel.
141     *
142     * @param lat0 latitude of natural origin (in degrees)
143     */
144    private void initialize1SP(double lat0) {
145        this.params = new Parameters1SP(lat0);
146        final double lat0rad = Utils.toRadians(lat0);
147
148        final double m0 = m(lat0rad);
149        final double t0 = t(lat0rad);
150
151        n = sin(lat0rad);
152        f = m0 / (n * pow(t0, n));
153        r0 = f * pow(t0, n);
154    }
155
156    /**
157     * auxiliary function t
158     * @param latRad latitude in radians
159     * @return result
160     */
161    protected double t(double latRad) {
162        return tan(PI/4 - latRad / 2.0)
163            / pow((1.0 - e * sin(latRad)) / (1.0 + e * sin(latRad)), e/2);
164    }
165
166    /**
167     * auxiliary function m
168     * @param latRad latitude in radians
169     * @return result
170     */
171    protected double m(double latRad) {
172        return cos(latRad) / (sqrt(1 - e * e * pow(sin(latRad), 2)));
173    }
174
175    @Override
176    public String getName() {
177        return tr("Lambert Conformal Conic");
178    }
179
180    @Override
181    public String getProj4Id() {
182        return "lcc";
183    }
184
185    @Override
186    public double[] project(double phi, double lambda) {
187        double sinphi = sin(phi);
188        double l = (0.5*log((1+sinphi)/(1-sinphi))) - e/2*log((1+e*sinphi)/(1-e*sinphi));
189        double r = f*exp(-n*l);
190        double gamma = n*lambda;
191        double x = r*sin(gamma);
192        double y = r0 - r*cos(gamma);
193        return new double[] {x, y};
194    }
195
196    @Override
197    public double[] invproject(double east, double north) {
198        double r = sqrt(pow(east, 2) + pow(north-r0, 2));
199        double gamma = atan(east / (r0-north));
200        double lambda = gamma/n;
201        double latIso = (-1/n) * log(abs(r/f));
202        double phi = ellps.latitude(latIso, e, epsilon);
203        return new double[] {phi, lambda};
204    }
205
206    /**
207     * Returns projection parameters.
208     * @return projection parameters
209     */
210    public final Parameters getParameters() {
211        return params;
212    }
213
214    @Override
215    public Bounds getAlgorithmBounds() {
216        double lat;
217        if (params instanceof Parameters2SP) {
218            Parameters2SP p2p = (Parameters2SP) params;
219            lat = (p2p.standardParallel1 + p2p.standardParallel2) / 2;
220        } else {
221            lat = params.latitudeOrigin;
222        }
223        double minlat = Math.max(lat - 60, -89);
224        double maxlat = Math.min(lat + 60, 89);
225        return new Bounds(minlat, -85, maxlat, 85, false);
226    }
227}