source: josm/trunk/src/org/openstreetmap/josm/data/projection/Lambert.java@ 1583

Last change on this file since 1583 was 1583, checked in by stoecker, 15 years ago

close #2535 - patch by podolsir - NPE for missing icon, update translations, added et, nb, gl
File size: 11.9 KB
Line 
1//License: GPL. For details, see LICENSE file.
2//Thanks to Johan Montagnat and its geoconv java converter application
3//(http://www.i3s.unice.fr/~johan/gps/ , published under GPL license)
4//from which some code and constants have been reused here.
5package org.openstreetmap.josm.data.projection;
6
7import static org.openstreetmap.josm.tools.I18n.tr;
8
9import javax.swing.JOptionPane;
10
11import org.openstreetmap.josm.Main;
12import org.openstreetmap.josm.data.coor.EastNorth;
13import org.openstreetmap.josm.data.coor.LatLon;
14
15public class Lambert implements Projection {
16 /**
17 * Lambert I, II, III, and IV projection exponents
18 */
19 public static final double n[] = { 0.7604059656, 0.7289686274, 0.6959127966, 0.6712679322 };
20
21 /**
22 * Lambert I, II, III, and IV projection constants
23 */
24 public static final double c[] = { 11603796.98, 11745793.39, 11947992.52, 12136281.99 };
25
26 /**
27 * Lambert I, II, III, and IV false east
28 */
29 public static final double Xs[] = { 600000.0, 600000.0, 600000.0, 234.358 };
30
31 /**
32 * Lambert I, II, III, and IV false north
33 */
34 public static final double Ys[] = { 5657616.674, 6199695.768, 6791905.085, 7239161.542 };
35
36 /**
37 * Lambert I, II, III, and IV longitudinal offset to Greenwich meridian
38 */
39 public static final double lg0 = 0.04079234433198; // 2deg20'14.025"
40
41 /**
42 * precision in iterative schema
43 */
44
45 public static final double epsilon = 1e-11;
46
47 /**
48 * France is divided in 4 Lambert projection zones (1,2,3 + 4th for Corsica)
49 */
50 public static final double cMaxLatZone1 = Math.toRadians(57 * 0.9);
51
52 public static final double zoneLimits[] = { Math.toRadians(53.5 * 0.9), // between Zone 1 and Zone 2 (in grad *0.9)
53 Math.toRadians(50.5 * 0.9), // between Zone 2 and Zone 3
54 Math.toRadians(47.51963 * 0.9), // between Zone 3 and Zone 4
55 Math.toRadians(46.17821 * 0.9) };// lowest latitude of Zone 4
56
57 public static final double cMinLonZones = Math.toRadians(-4.9074074074074059 * 0.9);
58
59 public static final double cMaxLonZones = Math.toRadians(10.2 * 0.9);
60
61 /**
62 * Because josm cannot work correctly if two zones are displayed, we allow some overlapping
63 */
64 public static final double cMaxOverlappingZones = Math.toRadians(1.5 * 0.9);
65
66 public static int layoutZone = -1;
67
68 private static int currentZone = 0;
69
70 private static boolean dontDisplayErrors = false;
71
72 /**
73 * @param p WGS84 lat/lon (ellipsoid GRS80) (in degree)
74 * @return eastnorth projection in Lambert Zone (ellipsoid Clark)
75 */
76 public EastNorth latlon2eastNorth(LatLon p) {
77 // translate ellipsoid GRS80 (WGS83) => Clark
78 LatLon geo = GRS802Clark(p);
79 double lt = geo.lat(); // in radian
80 double lg = geo.lon();
81
82 // check if longitude and latitude are inside the French Lambert zones
83 currentZone = 0;
84 boolean outOfLambertZones = false;
85 if (lt >= zoneLimits[3] && lt <= cMaxLatZone1 && lg >= cMinLonZones && lg <= cMaxLonZones) {
86 // zone I
87 if (lt > zoneLimits[0])
88 currentZone = 0;
89 // zone II
90 else if (lt > zoneLimits[1])
91 currentZone = 1;
92 // zone III
93 else if (lt > zoneLimits[2])
94 currentZone = 2;
95 // zone III or IV
96 else if (lt > zoneLimits[3])
97 // Note: zone IV is dedicated to Corsica island and extends from 47.8 to
98 // 45.9 degrees of latitude. There is an overlap with zone III that can be
99 // solved only with longitude (covers Corsica if lon > 7.2 degree)
100 if (lg < Math.toRadians(8 * 0.9))
101 currentZone = 2;
102 else
103 currentZone = 3;
104 } else {
105 outOfLambertZones = true; // possible when MAX_LAT is used
106 if (p.lat() != 0 && Math.abs(p.lat()) != Projection.MAX_LAT
107 && p.lon() != 0 && Math.abs(p.lon()) != Projection.MAX_LON
108 && dontDisplayErrors == false) {
109 JOptionPane.showMessageDialog(Main.parent,
110 tr("The projection \"{0}\" is designed for\n"
111 + "latitudes between 46.1\u00b0 and 57\u00b0 only.\n"
112 + "Use another projection system if you are not using\n"
113 + "a French WMS server.\n"
114 + "Do not upload any data after this message.", this.toString()));
115 dontDisplayErrors = true;
116 }
117 }
118 if (!outOfLambertZones) {
119 if (layoutZone == -1) {
120 layoutZone = currentZone;
121 dontDisplayErrors = false;
122 } else if (layoutZone != currentZone) {
123 if ((currentZone < layoutZone && Math.abs(zoneLimits[currentZone] - lt) > cMaxOverlappingZones)
124 || (currentZone > layoutZone && Math.abs(zoneLimits[layoutZone] - lt) > cMaxOverlappingZones)) {
125 JOptionPane.showMessageDialog(Main.parent,
126 tr("IMPORTANT : data positioned far away from\n"
127 + "the current Lambert zone limits.\n"
128 + "Do not upload any data after this message.\n"
129 + "Undo your last action, save your work\n"
130 + "and start a new layer on the new zone."));
131 layoutZone = -1;
132 dontDisplayErrors = true;
133 } else {
134 System.out.println("temporarily extend Lambert zone " + layoutZone + " projection at lat,lon:"
135 + lt + "," + lg);
136 }
137 }
138 }
139 if (layoutZone == -1) {
140 return ConicProjection(lt, lg, Xs[currentZone], Ys[currentZone], c[currentZone], n[currentZone]);
141 } // else
142 return ConicProjection(lt, lg, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]);
143 }
144
145 public LatLon eastNorth2latlon(EastNorth p) {
146 LatLon geo;
147 if (layoutZone == -1)
148 // possible until the Lambert zone is determined by latlon2eastNorth() with a valid LatLon
149 geo = Geographic(p, Xs[currentZone], Ys[currentZone], c[currentZone], n[currentZone]);
150 else
151 geo = Geographic(p, Xs[layoutZone], Ys[layoutZone], c[layoutZone], n[layoutZone]);
152 // translate ellipsoid Clark => GRS80 (WGS83)
153 LatLon wgs = Clark2GRS80(geo);
154 return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon()));
155 }
156
157 @Override public String toString() {
158 return tr("Lambert Zone (France)");
159 }
160
161 public String toCode() {
162 return "EPSG::"+(27571+currentZone);
163 }
164
165 public String getCacheDirectoryName() {
166 return "lambert";
167 }
168
169 public double scaleFactor() {
170 return 1.0 / 360;
171 }
172
173 @Override
174 public boolean equals(Object o) {
175 return o instanceof Lambert;
176 }
177
178 @Override
179 public int hashCode() {
180 return Lambert.class.hashCode();
181 }
182
183 /**
184 * Initializes from geographic coordinates. Note that reference ellipsoid
185 * used by Lambert is the Clark ellipsoid.
186 *
187 * @param lat latitude in grad
188 * @param lon longitude in grad
189 * @param Xs false east (coordinate system origin) in meters
190 * @param Ys false north (coordinate system origin) in meters
191 * @param c projection constant
192 * @param n projection exponent
193 * @return EastNorth projected coordinates in meter
194 */
195 private EastNorth ConicProjection(double lat, double lon, double Xs, double Ys, double c, double n) {
196 double eslt = Ellipsoid.clarke.e * Math.sin(lat);
197 double l = Math.log(Math.tan(Math.PI / 4.0 + (lat / 2.0))
198 * Math.pow((1.0 - eslt) / (1.0 + eslt), Ellipsoid.clarke.e / 2.0));
199 double east = Xs + c * Math.exp(-n * l) * Math.sin(n * (lon - lg0));
200 double north = Ys - c * Math.exp(-n * l) * Math.cos(n * (lon - lg0));
201 return new EastNorth(east, north);
202 }
203
204 /**
205 * Initializes from projected coordinates (conic projection). Note that
206 * reference ellipsoid used by Lambert is Clark
207 *
208 * @param eastNorth projected coordinates pair in meters
209 * @param Xs false east (coordinate system origin) in meters
210 * @param Ys false north (coordinate system origin) in meters
211 * @param c projection constant
212 * @param n projection exponent
213 * @return LatLon in radian
214 */
215 private LatLon Geographic(EastNorth eastNorth, double Xs, double Ys, double c, double n) {
216 double dx = eastNorth.east() - Xs;
217 double dy = Ys - eastNorth.north();
218 double R = Math.sqrt(dx * dx + dy * dy);
219 double gamma = Math.atan(dx / dy);
220 double l = -1.0 / n * Math.log(Math.abs(R / c));
221 l = Math.exp(l);
222 double lon = lg0 + gamma / n;
223 double lat = 2.0 * Math.atan(l) - Math.PI / 2.0;
224 double delta = 1.0;
225 while (delta > epsilon) {
226 double eslt = Ellipsoid.clarke.e * Math.sin(lat);
227 double nlt = 2.0 * Math.atan(Math.pow((1.0 + eslt) / (1.0 - eslt), Ellipsoid.clarke.e / 2.0) * l) - Math.PI
228 / 2.0;
229 delta = Math.abs(nlt - lat);
230 lat = nlt;
231 }
232 return new LatLon(lat, lon); // in radian
233 }
234
235 /**
236 * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to Lambert
237 * geographic, (ellipsoid Clark)
238 */
239 private LatLon GRS802Clark(LatLon wgs) {
240 double lat = Math.toRadians(wgs.lat()); // degree to radian
241 double lon = Math.toRadians(wgs.lon());
242 // WGS84 geographic => WGS84 cartesian
243 double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat)));
244 double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
245 double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
246 double Z = (N * (1.0 - Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat);
247 // WGS84 => Lambert ellipsoide similarity
248 X += 168.0;
249 Y += 60.0;
250 Z += -320.0;
251 // Lambert cartesian => Lambert geographic
252 return Geographic(X, Y, Z, Ellipsoid.clarke);
253 }
254
255 private LatLon Clark2GRS80(LatLon lambert) {
256 double lat = lambert.lat(); // in radian
257 double lon = lambert.lon();
258 // Lambert geographic => Lambert cartesian
259 double N = Ellipsoid.clarke.a / (Math.sqrt(1.0 - Ellipsoid.clarke.e2 * Math.sin(lat) * Math.sin(lat)));
260 double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
261 double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
262 double Z = (N * (1.0 - Ellipsoid.clarke.e2)/* + height */) * Math.sin(lat);
263 // Lambert => WGS84 ellipsoide similarity
264 X += -168.0;
265 Y += -60.0;
266 Z += 320.0;
267 // WGS84 cartesian => WGS84 geographic
268 return Geographic(X, Y, Z, Ellipsoid.GRS80);
269 }
270
271 /**
272 * initializes from cartesian coordinates
273 *
274 * @param X
275 * 1st coordinate in meters
276 * @param Y
277 * 2nd coordinate in meters
278 * @param Z
279 * 3rd coordinate in meters
280 * @param ell
281 * reference ellipsoid
282 */
283 private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) {
284 double norm = Math.sqrt(X * X + Y * Y);
285 double lg = 2.0 * Math.atan(Y / (X + norm));
286 double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
287 double delta = 1.0;
288 while (delta > epsilon) {
289 double s2 = Math.sin(lt);
290 s2 *= s2;
291 double l = Math.atan((Z / norm)
292 / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
293 delta = Math.abs(l - lt);
294 lt = l;
295 }
296 double s2 = Math.sin(lt);
297 s2 *= s2;
298 // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
299 return new LatLon(lt, lg);
300 }
301
302}
Note: See TracBrowser for help on using the repository browser.