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

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

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