1 | // License: GPL. For details, see LICENSE file.
|
---|
2 | package org.openstreetmap.josm.data.projection.proj;
|
---|
3 |
|
---|
4 | import static java.lang.Math.abs;
|
---|
5 | import static java.lang.Math.atan;
|
---|
6 | import static java.lang.Math.cos;
|
---|
7 | import static java.lang.Math.exp;
|
---|
8 | import static java.lang.Math.log;
|
---|
9 | import static java.lang.Math.pow;
|
---|
10 | import static java.lang.Math.sin;
|
---|
11 | import static java.lang.Math.sqrt;
|
---|
12 | import static org.openstreetmap.josm.tools.I18n.tr;
|
---|
13 | import static org.openstreetmap.josm.tools.Utils.toRadians;
|
---|
14 |
|
---|
15 | import org.openstreetmap.josm.data.Bounds;
|
---|
16 | import org.openstreetmap.josm.data.projection.CustomProjection.Param;
|
---|
17 | import org.openstreetmap.josm.data.projection.Ellipsoid;
|
---|
18 | import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
|
---|
19 |
|
---|
20 | /**
|
---|
21 | * Lambert Conical Conformal Projection. Areas and shapes are deformed as one moves away from standard parallels.
|
---|
22 | * The angles are true in a limited area. This projection is used for the charts of North America, France and Belgium.
|
---|
23 | * <p>
|
---|
24 | * This implementation provides transforms for two cases of the lambert conic conformal projection:
|
---|
25 | * </p>
|
---|
26 | * <ul>
|
---|
27 | * <li>{@code Lambert_Conformal_Conic_1SP} (EPSG code 9801)</li>
|
---|
28 | * <li>{@code Lambert_Conformal_Conic_2SP} (EPSG code 9802)</li>
|
---|
29 | * </ul>
|
---|
30 | * <p>
|
---|
31 | * For the 1SP case the latitude of origin is used as the standard parallel (SP).
|
---|
32 | * To use 1SP with a latitude of origin different from the SP, use the 2SP and set the SP1 to the single SP.
|
---|
33 | * The {@code "standard_parallel_2"} parameter is optional and will be given the same value
|
---|
34 | * as {@code "standard_parallel_1"} if not set (creating a 1 standard parallel projection).
|
---|
35 | * </p>
|
---|
36 | * <b>References:</b>
|
---|
37 | * <ul>
|
---|
38 | * <li>John P. Snyder (Map Projections - A Working Manual,<br>U.S. Geological Survey Professional Paper 1395, 1987)</li>
|
---|
39 | * <li>"Coordinate Conversions and Transformations including Formulas",<br>EPSG Guidence Note Number 7, Version 19.</li>
|
---|
40 | * </ul>
|
---|
41 | *
|
---|
42 | * @author Pieren
|
---|
43 | * @author André Gosselin
|
---|
44 | * @author Martin Desruisseaux (PMO, IRD)
|
---|
45 | * @author Rueben Schulz
|
---|
46 | *
|
---|
47 | * @see <A HREF="http://mathworld.wolfram.com/LambertConformalConicProjection.html">Lambert conformal conic projection on MathWorld</A>
|
---|
48 | * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html">lambert_conic_conformal_1sp</A>
|
---|
49 | * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_2sp.html">lambert_conic_conformal_2sp</A>
|
---|
50 | *
|
---|
51 | * @since 13639 (align implementation with proj.4 / GeoTools)
|
---|
52 | * @since 4285 (reworked from Lambert / LambertCC9Zones)
|
---|
53 | * @since 2304 (initial implementation by Pieren)
|
---|
54 | */
|
---|
55 | public class LambertConformalConic extends AbstractProj {
|
---|
56 |
|
---|
57 | /** ellipsoid */
|
---|
58 | protected Ellipsoid ellps;
|
---|
59 |
|
---|
60 | /**
|
---|
61 | * Base class of Lambert Conformal Conic parameters.
|
---|
62 | */
|
---|
63 | public static class Parameters {
|
---|
64 | /** latitude of origin */
|
---|
65 | public final double latitudeOrigin;
|
---|
66 |
|
---|
67 | /**
|
---|
68 | * Constructs a new {@code Parameters}.
|
---|
69 | * @param latitudeOrigin latitude of origin
|
---|
70 | */
|
---|
71 | protected Parameters(double latitudeOrigin) {
|
---|
72 | this.latitudeOrigin = latitudeOrigin;
|
---|
73 | }
|
---|
74 | }
|
---|
75 |
|
---|
76 | /**
|
---|
77 | * Parameters with a single standard parallel.
|
---|
78 | */
|
---|
79 | public static class Parameters1SP extends Parameters {
|
---|
80 | /**
|
---|
81 | * Constructs a new {@code Parameters1SP}.
|
---|
82 | * @param latitudeOrigin latitude of origin
|
---|
83 | */
|
---|
84 | public Parameters1SP(double latitudeOrigin) {
|
---|
85 | super(latitudeOrigin);
|
---|
86 | }
|
---|
87 | }
|
---|
88 |
|
---|
89 | /**
|
---|
90 | * Parameters with two standard parallels.
|
---|
91 | */
|
---|
92 | public static class Parameters2SP extends Parameters {
|
---|
93 | /** first standard parallel */
|
---|
94 | public final double standardParallel1;
|
---|
95 | /** second standard parallel */
|
---|
96 | public final double standardParallel2;
|
---|
97 |
|
---|
98 | /**
|
---|
99 | * 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 | }
|
---|