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.
|
---|
5 | package org.openstreetmap.josm.data.projection;
|
---|
6 |
|
---|
7 | import static org.openstreetmap.josm.tools.I18n.tr;
|
---|
8 |
|
---|
9 | import javax.swing.JOptionPane;
|
---|
10 |
|
---|
11 | import org.openstreetmap.josm.Main;
|
---|
12 | import org.openstreetmap.josm.data.coor.EastNorth;
|
---|
13 | import org.openstreetmap.josm.data.coor.LatLon;
|
---|
14 |
|
---|
15 | public 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° and 57° 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 positionned 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 getCacheDirectoryName() {
|
---|
162 | return "lambert";
|
---|
163 | }
|
---|
164 |
|
---|
165 | public double scaleFactor() {
|
---|
166 | return 1.0 / 360;
|
---|
167 | }
|
---|
168 |
|
---|
169 | @Override
|
---|
170 | public boolean equals(Object o) {
|
---|
171 | return o instanceof Lambert;
|
---|
172 | }
|
---|
173 |
|
---|
174 | @Override
|
---|
175 | public int hashCode() {
|
---|
176 | return Lambert.class.hashCode();
|
---|
177 | }
|
---|
178 |
|
---|
179 | /**
|
---|
180 | * Initializes from geographic coordinates. Note that reference ellipsoid
|
---|
181 | * used by Lambert is the Clark ellipsoid.
|
---|
182 | *
|
---|
183 | * @param lat latitude in grad
|
---|
184 | * @param lon longitude in grad
|
---|
185 | * @param Xs false east (coordinate system origin) in meters
|
---|
186 | * @param Ys false north (coordinate system origin) in meters
|
---|
187 | * @param c projection constant
|
---|
188 | * @param n projection exponent
|
---|
189 | * @return EastNorth projected coordinates in meter
|
---|
190 | */
|
---|
191 | private EastNorth ConicProjection(double lat, double lon, double Xs, double Ys, double c, double n) {
|
---|
192 | double eslt = Ellipsoid.clarke.e * Math.sin(lat);
|
---|
193 | double l = Math.log(Math.tan(Math.PI / 4.0 + (lat / 2.0))
|
---|
194 | * Math.pow((1.0 - eslt) / (1.0 + eslt), Ellipsoid.clarke.e / 2.0));
|
---|
195 | double east = Xs + c * Math.exp(-n * l) * Math.sin(n * (lon - lg0));
|
---|
196 | double north = Ys - c * Math.exp(-n * l) * Math.cos(n * (lon - lg0));
|
---|
197 | return new EastNorth(east, north);
|
---|
198 | }
|
---|
199 |
|
---|
200 | /**
|
---|
201 | * Initializes from projected coordinates (conic projection). Note that
|
---|
202 | * reference ellipsoid used by Lambert is Clark
|
---|
203 | *
|
---|
204 | * @param coord projected coordinates pair in meters
|
---|
205 | * @param Xs false east (coordinate system origin) in meters
|
---|
206 | * @param Ys false north (coordinate system origin) in meters
|
---|
207 | * @param c projection constant
|
---|
208 | * @param n projection exponent
|
---|
209 | * @return LatLon in radian
|
---|
210 | */
|
---|
211 | private LatLon Geographic(EastNorth eastNorth, double Xs, double Ys, double c, double n) {
|
---|
212 | double dx = eastNorth.east() - Xs;
|
---|
213 | double dy = Ys - eastNorth.north();
|
---|
214 | double R = Math.sqrt(dx * dx + dy * dy);
|
---|
215 | double gamma = Math.atan(dx / dy);
|
---|
216 | double l = -1.0 / n * Math.log(Math.abs(R / c));
|
---|
217 | l = Math.exp(l);
|
---|
218 | double lon = lg0 + gamma / n;
|
---|
219 | double lat = 2.0 * Math.atan(l) - Math.PI / 2.0;
|
---|
220 | double delta = 1.0;
|
---|
221 | while (delta > epsilon) {
|
---|
222 | double eslt = Ellipsoid.clarke.e * Math.sin(lat);
|
---|
223 | double nlt = 2.0 * Math.atan(Math.pow((1.0 + eslt) / (1.0 - eslt), Ellipsoid.clarke.e / 2.0) * l) - Math.PI
|
---|
224 | / 2.0;
|
---|
225 | delta = Math.abs(nlt - lat);
|
---|
226 | lat = nlt;
|
---|
227 | }
|
---|
228 | return new LatLon(lat, lon); // in radian
|
---|
229 | }
|
---|
230 |
|
---|
231 | /**
|
---|
232 | * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to Lambert
|
---|
233 | * geographic, (ellipsoid Clark)
|
---|
234 | *
|
---|
235 | * @param wgs
|
---|
236 | * @return
|
---|
237 | */
|
---|
238 | private LatLon GRS802Clark(LatLon wgs) {
|
---|
239 | double lat = Math.toRadians(wgs.lat()); // degree to radian
|
---|
240 | double lon = Math.toRadians(wgs.lon());
|
---|
241 | // WGS84 geographic => WGS84 cartesian
|
---|
242 | double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat)));
|
---|
243 | double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
|
---|
244 | double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
|
---|
245 | double Z = (N * (1.0 - Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat);
|
---|
246 | // WGS84 => Lambert ellipsoide similarity
|
---|
247 | X += 168.0;
|
---|
248 | Y += 60.0;
|
---|
249 | Z += -320.0;
|
---|
250 | // Lambert cartesian => Lambert geographic
|
---|
251 | return Geographic(X, Y, Z, Ellipsoid.clarke);
|
---|
252 | }
|
---|
253 |
|
---|
254 | /**
|
---|
255 | * @param lambert
|
---|
256 | * @return
|
---|
257 | */
|
---|
258 | private LatLon Clark2GRS80(LatLon lambert) {
|
---|
259 | double lat = lambert.lat(); // in radian
|
---|
260 | double lon = lambert.lon();
|
---|
261 | // Lambert geographic => Lambert cartesian
|
---|
262 | double N = Ellipsoid.clarke.a / (Math.sqrt(1.0 - Ellipsoid.clarke.e2 * Math.sin(lat) * Math.sin(lat)));
|
---|
263 | double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
|
---|
264 | double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
|
---|
265 | double Z = (N * (1.0 - Ellipsoid.clarke.e2)/* + height */) * Math.sin(lat);
|
---|
266 | // Lambert => WGS84 ellipsoide similarity
|
---|
267 | X += -168.0;
|
---|
268 | Y += -60.0;
|
---|
269 | Z += 320.0;
|
---|
270 | // WGS84 cartesian => WGS84 geographic
|
---|
271 | return Geographic(X, Y, Z, Ellipsoid.GRS80);
|
---|
272 | }
|
---|
273 |
|
---|
274 | /**
|
---|
275 | * initializes from cartesian coordinates
|
---|
276 | *
|
---|
277 | * @param X
|
---|
278 | * 1st coordinate in meters
|
---|
279 | * @param Y
|
---|
280 | * 2nd coordinate in meters
|
---|
281 | * @param Z
|
---|
282 | * 3rd coordinate in meters
|
---|
283 | * @param ell
|
---|
284 | * reference ellipsoid
|
---|
285 | */
|
---|
286 | private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) {
|
---|
287 | double norm = Math.sqrt(X * X + Y * Y);
|
---|
288 | double lg = 2.0 * Math.atan(Y / (X + norm));
|
---|
289 | double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
|
---|
290 | double delta = 1.0;
|
---|
291 | while (delta > epsilon) {
|
---|
292 | double s2 = Math.sin(lt);
|
---|
293 | s2 *= s2;
|
---|
294 | double l = Math.atan((Z / norm)
|
---|
295 | / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
|
---|
296 | delta = Math.abs(l - lt);
|
---|
297 | lt = l;
|
---|
298 | }
|
---|
299 | double s2 = Math.sin(lt);
|
---|
300 | s2 *= s2;
|
---|
301 | // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
|
---|
302 | return new LatLon(lt, lg);
|
---|
303 | }
|
---|
304 |
|
---|
305 | }
|
---|