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