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