source: josm/trunk/src/org/openstreetmap/josm/data/projection/UTM_20N_France_DOM.java@ 2516

Last change on this file since 2516 was 2516, checked in by stoecker, 14 years ago

close #4015 - Zoomlevel changes whenever the preference dialog is closed
File size: 16.1 KB
Line 
1// License: GPL. For details, see LICENSE file.
2package org.openstreetmap.josm.data.projection;
3
4/**
5 * This class implements all projections for French departements in the Caribbean Sea using
6 * UTM zone 20N transvers Mercator and specific geodesic settings (7 parameters transformation algorithm).
7 */
8import static org.openstreetmap.josm.tools.I18n.tr;
9
10import java.awt.GridBagLayout;
11import java.util.Collection;
12import java.util.Collections;
13
14import javax.swing.JComboBox;
15import javax.swing.JLabel;
16import javax.swing.JPanel;
17
18import org.openstreetmap.josm.data.Bounds;
19import org.openstreetmap.josm.data.coor.EastNorth;
20import org.openstreetmap.josm.data.coor.LatLon;
21import org.openstreetmap.josm.tools.GBC;
22
23public class UTM_20N_France_DOM implements Projection, ProjectionSubPrefs {
24
25 private static String FortMarigotName = tr("Guadeloupe Fort-Marigot 1949");
26 private static String SainteAnneName = tr("Guadeloupe Ste-Anne 1948");
27 private static String MartiniqueName = tr("Martinique Fort Desaix 1952");
28 public static String[] utmGeodesicsNames = { FortMarigotName, SainteAnneName, MartiniqueName};
29
30 private Bounds FortMarigotBounds = new Bounds( new LatLon(17.6,-63.25), new LatLon(18.5,-62.5));
31 private Bounds SainteAnneBounds = new Bounds( new LatLon(15.8,-61.9), new LatLon(16.6,-60.9));
32 private Bounds MartiniqueBounds = new Bounds( new LatLon(14.25,-61.25), new LatLon(15.025,-60.725));
33 private Bounds[] utmBounds = { FortMarigotBounds, SainteAnneBounds, MartiniqueBounds};
34
35 private String FortMarigotEPSG = "EPSG::2969";
36 private String SainteAnneEPSG = "EPSG::2970";
37 private String MartiniqueEPSG = "EPSG::2973";
38 private String[] utmEPSGs = { FortMarigotEPSG, SainteAnneEPSG, MartiniqueEPSG};
39
40 /**
41 * false east in meters (constant)
42 */
43 private static final double Xs = 500000.0;
44 /**
45 * false north in meters (0 in northern hemisphere, 10000000 in southern
46 * hemisphere)
47 */
48 private static double Ys = 0;
49 /**
50 * origin meridian longitude
51 */
52 protected double lg0;
53 /**
54 * UTM zone (from 1 to 60)
55 */
56 private static int ZONE = 20;
57 /**
58 * whether north or south hemisphere
59 */
60 private boolean isNorth = true;
61
62 public static final int DEFAULT_GEODESIC = 0;
63
64 private static int currentGeodesic = DEFAULT_GEODESIC;
65
66 /**
67 * 7 parameters transformation
68 */
69 private static double tx = 0.0;
70 private static double ty = 0.0;
71 private static double tz = 0.0;
72 private static double rx = 0;
73 private static double ry = 0;
74 private static double rz = 0;
75 private static double scaleDiff = 0;
76 /**
77 * precision in iterative schema
78 */
79 public static final double epsilon = 1e-11;
80 public final static double DEG_TO_RAD = Math.PI / 180;
81 public final static double RAD_TO_DEG = 180 / Math.PI;
82
83 private void refresh7ParametersTranslation() {
84 //System.out.println("Current UTM geodesic system: " + utmGeodesicsNames[currentGeodesic]);
85 if (currentGeodesic == 0) { // UTM_20N_Guadeloupe_Fort_Marigot
86 set7ParametersTranslation(new double[]{136.596, 248.148, -429.789},
87 new double[]{0, 0, 0},
88 0);
89 } else if (currentGeodesic == 1) { // UTM_20N_Guadeloupe_Ste_Anne
90 set7ParametersTranslation(new double[]{-472.29, -5.63, -304.12},
91 new double[]{0.4362, -0.8374, 0.2563},
92 1.8984E-6);
93 } else { // UTM_20N_Martinique_Fort_Desaix
94 set7ParametersTranslation(new double[]{126.926, 547.939, 130.409},
95 new double[]{-2.78670, 5.16124, -0.85844},
96 13.82265E-6);
97 }
98 }
99
100 private void set7ParametersTranslation(double[] translation, double[] rotation, double scalediff) {
101 tx = translation[0];
102 ty = translation[1];
103 tz = translation[2];
104 rx = rotation[0]/206264.806247096355; // seconds to radian
105 ry = rotation[1]/206264.806247096355;
106 rz = rotation[2]/206264.806247096355;
107 scaleDiff = scalediff;
108 }
109
110 public EastNorth latlon2eastNorth(LatLon p) {
111 // translate ellipsoid GRS80 (WGS83) => reference ellipsoid geographic
112 LatLon geo = GRS802Hayford(p);
113
114 // reference ellipsoid geographic => UTM projection
115 return MTProjection(geo, Ellipsoid.hayford.a, Ellipsoid.hayford.e);
116 }
117
118 /**
119 * Translate latitude/longitude in WGS84, (ellipsoid GRS80) to UTM
120 * geographic, (ellipsoid Hayford)
121 */
122 private LatLon GRS802Hayford(LatLon wgs) {
123 double lat = Math.toRadians(wgs.lat()); // degree to radian
124 double lon = Math.toRadians(wgs.lon());
125 // WGS84 geographic => WGS84 cartesian
126 double N = Ellipsoid.GRS80.a / (Math.sqrt(1.0 - Ellipsoid.GRS80.e2 * Math.sin(lat) * Math.sin(lat)));
127 double X = (N/* +height */) * Math.cos(lat) * Math.cos(lon);
128 double Y = (N/* +height */) * Math.cos(lat) * Math.sin(lon);
129 double Z = (N * (1.0 - Ellipsoid.GRS80.e2)/* + height */) * Math.sin(lat);
130 // translation
131 double coord[] = invSevenParametersTransformation(X, Y, Z);
132 // UTM cartesian => UTM geographic
133 return Geographic(coord[0], coord[1], coord[2], Ellipsoid.hayford);
134 }
135
136 /**
137 * initializes from cartesian coordinates
138 *
139 * @param X
140 * 1st coordinate in meters
141 * @param Y
142 * 2nd coordinate in meters
143 * @param Z
144 * 3rd coordinate in meters
145 * @param ell
146 * reference ellipsoid
147 */
148 private LatLon Geographic(double X, double Y, double Z, Ellipsoid ell) {
149 double norm = Math.sqrt(X * X + Y * Y);
150 double lg = 2.0 * Math.atan(Y / (X + norm));
151 double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
152 double delta = 1.0;
153 while (delta > epsilon) {
154 double s2 = Math.sin(lt);
155 s2 *= s2;
156 double l = Math.atan((Z / norm)
157 / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
158 delta = Math.abs(l - lt);
159 lt = l;
160 }
161 double s2 = Math.sin(lt);
162 s2 *= s2;
163 // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
164 return new LatLon(lt, lg);
165 }
166
167 /**
168 * initalizes from geographic coordinates
169 *
170 * @param coord geographic coordinates triplet
171 * @param a reference ellipsoid long axis
172 * @param e reference ellipsoid excentricity
173 */
174 private EastNorth MTProjection(LatLon coord, double a, double e) {
175 double n = 0.9996 * a;
176 Ys = (coord.lat() >= 0.0) ? 0.0 : 10000000.0;
177 double r6d = Math.PI / 30.0;
178 //zone = (int) Math.floor((coord.lon() + Math.PI) / r6d) + 1;
179 lg0 = r6d * (ZONE - 0.5) - Math.PI;
180 double e2 = e * e;
181 double e4 = e2 * e2;
182 double e6 = e4 * e2;
183 double e8 = e4 * e4;
184 double C[] = {
185 1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0 - 175.0*e8/16384.0,
186 e2/8.0 - e4/96.0 - 9.0*e6/1024.0 - 901.0*e8/184320.0,
187 13.0*e4/768.0 + 17.0*e6/5120.0 - 311.0*e8/737280.0,
188 61.0*e6/15360.0 + 899.0*e8/430080.0,
189 49561.0*e8/41287680.0
190 };
191 double s = e * Math.sin(coord.lat());
192 double l = Math.log(Math.tan(Math.PI/4.0 + coord.lat()/2.0) *
193 Math.pow((1.0 - s) / (1.0 + s), e/2.0));
194 double phi = Math.asin(Math.sin(coord.lon() - lg0) /
195 ((Math.exp(l) + Math.exp(-l)) / 2.0));
196 double ls = Math.log(Math.tan(Math.PI/4.0 + phi/2.0));
197 double lambda = Math.atan(((Math.exp(l) - Math.exp(-l)) / 2.0) /
198 Math.cos(coord.lon() - lg0));
199
200 double north = C[0] * lambda;
201 double east = C[0] * ls;
202 for(int k = 1; k < 5; k++) {
203 double r = 2.0 * k * lambda;
204 double m = 2.0 * k * ls;
205 double em = Math.exp(m);
206 double en = Math.exp(-m);
207 double sr = Math.sin(r)/2.0 * (em + en);
208 double sm = Math.cos(r)/2.0 * (em - en);
209 north += C[k] * sr;
210 east += C[k] * sm;
211 }
212 east *= n;
213 east += Xs;
214 north *= n;
215 north += Ys;
216 return new EastNorth(east, north);
217 }
218
219 public LatLon eastNorth2latlon(EastNorth p) {
220 MTProjection(p.east(), p.north(), ZONE, isNorth);
221 LatLon geo = Geographic(p, Ellipsoid.hayford.a, Ellipsoid.hayford.e, 0.0 /* z */);
222
223 // reference ellipsoid geographic => reference ellipsoid cartesian
224 double N = Ellipsoid.hayford.a / (Math.sqrt(1.0 - Ellipsoid.hayford.e2 * Math.sin(geo.lat()) * Math.sin(geo.lat())));
225 double X = (N /*+ h*/) * Math.cos(geo.lat()) * Math.cos(geo.lon());
226 double Y = (N /*+ h*/) * Math.cos(geo.lat()) * Math.sin(geo.lon());
227 double Z = (N * (1.0-Ellipsoid.hayford.e2) /*+ h*/) * Math.sin(geo.lat());
228
229 // translation
230 double coord[] = sevenParametersTransformation(X, Y, Z);
231
232 // WGS84 cartesian => WGS84 geographic
233 LatLon wgs = cart2LatLon(coord[0], coord[1], coord[2], Ellipsoid.GRS80);
234 return new LatLon(Math.toDegrees(wgs.lat()), Math.toDegrees(wgs.lon()));
235 }
236
237 /**
238 * initializes new projection coordinates (in north hemisphere)
239 *
240 * @param east east from origin in meters
241 * @param north north from origin in meters
242 * @param zone zone number (from 1 to 60)
243 * @param isNorth true in north hemisphere, false in south hemisphere
244 */
245 private void MTProjection(double east, double north, int zone, boolean isNorth) {
246 Ys = isNorth ? 0.0 : 10000000.0;
247 double r6d = Math.PI / 30.0;
248 lg0 = r6d * (zone - 0.5) - Math.PI;
249 }
250
251 public double scaleFactor() {
252 return 1/Math.PI/2;
253 }
254
255 /**
256 * initalizes from projected coordinates (Mercator transverse projection)
257 *
258 * @param coord projected coordinates pair
259 * @param e reference ellipsoid excentricity
260 * @param a reference ellipsoid long axis
261 * @param z altitude in meters
262 */
263 private LatLon Geographic(EastNorth coord, double a, double e, double z) {
264 double n = 0.9996 * a;
265 double e2 = e * e;
266 double e4 = e2 * e2;
267 double e6 = e4 * e2;
268 double e8 = e4 * e4;
269 double C[] = {
270 1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0 - 175.0*e8/16384.0,
271 e2/8.0 + e4/48.0 + 7.0*e6/2048.0 + e8/61440.0,
272 e4/768.0 + 3.0*e6/1280.0 + 559.0*e8/368640.0,
273 17.0*e6/30720.0 + 283.0*e8/430080.0,
274 4397.0*e8/41287680.0
275 };
276 double l = (coord.north() - Ys) / (n * C[0]);
277 double ls = (coord.east() - Xs) / (n * C[0]);
278 double l0 = l;
279 double ls0 = ls;
280 for(int k = 1; k < 5; k++) {
281 double r = 2.0 * k * l0;
282 double m = 2.0 * k * ls0;
283 double em = Math.exp(m);
284 double en = Math.exp(-m);
285 double sr = Math.sin(r)/2.0 * (em + en);
286 double sm = Math.cos(r)/2.0 * (em - en);
287 l -= C[k] * sr;
288 ls -= C[k] * sm;
289 }
290 double lon = lg0 + Math.atan(((Math.exp(ls) - Math.exp(-ls)) / 2.0) /
291 Math.cos(l));
292 double phi = Math.asin(Math.sin(l) /
293 ((Math.exp(ls) + Math.exp(-ls)) / 2.0));
294 l = Math.log(Math.tan(Math.PI/4.0 + phi/2.0));
295 double lat = 2.0 * Math.atan(Math.exp(l)) - Math.PI / 2.0;
296 double lt0;
297 do {
298 lt0 = lat;
299 double s = e * Math.sin(lat);
300 lat = 2.0 * Math.atan(Math.pow((1.0 + s) / (1.0 - s), e/2.0) *
301 Math.exp(l)) - Math.PI / 2.0;
302 }
303 while(Math.abs(lat-lt0) >= epsilon);
304 //h = z;
305
306 return new LatLon(lat, lon);
307 }
308
309 /**
310 * initializes from cartesian coordinates
311 *
312 * @param X 1st coordinate in meters
313 * @param Y 2nd coordinate in meters
314 * @param Z 3rd coordinate in meters
315 * @param ell reference ellipsoid
316 */
317 private LatLon cart2LatLon(double X, double Y, double Z, Ellipsoid ell) {
318 double norm = Math.sqrt(X * X + Y * Y);
319 double lg = 2.0 * Math.atan(Y / (X + norm));
320 double lt = Math.atan(Z / (norm * (1.0 - (ell.a * ell.e2 / Math.sqrt(X * X + Y * Y + Z * Z)))));
321 double delta = 1.0;
322 while (delta > epsilon) {
323 double s2 = Math.sin(lt);
324 s2 *= s2;
325 double l = Math.atan((Z / norm)
326 / (1.0 - (ell.a * ell.e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - ell.e2 * s2)))));
327 delta = Math.abs(l - lt);
328 lt = l;
329 }
330 double s2 = Math.sin(lt);
331 s2 *= s2;
332 // h = norm / Math.cos(lt) - ell.a / Math.sqrt(1.0 - ell.e2 * s2);
333 return new LatLon(lt, lg);
334 }
335
336 /**
337 * 7 parameters transformation
338 * @param coord X, Y, Z in array
339 * @return transformed X, Y, Z in array
340 */
341 private double[] sevenParametersTransformation(double Xa, double Ya, double Za){
342 double Xb = tx + Xa*(1+scaleDiff) + Za*ry - Ya*rz;
343 double Yb = ty + Ya*(1+scaleDiff) + Xa*rz - Za*rx;
344 double Zb = tz + Za*(1+scaleDiff) + Ya*rx - Xa*ry;
345 return new double[]{Xb, Yb, Zb};
346 }
347
348 /**
349 * 7 parameters inverse transformation
350 * @param coord X, Y, Z in array
351 * @return transformed X, Y, Z in array
352 */
353 private double[] invSevenParametersTransformation(double Xa, double Ya, double Za){
354 double Xb = (1-scaleDiff)*(-tx + Xa + ((-tz+Za)*(-ry) - (-ty+Ya)*(-rz)));
355 double Yb = (1-scaleDiff)*(-ty + Ya + ((-tx+Xa)*(-rz) - (-tz+Za)*(-rx)));
356 double Zb = (1-scaleDiff)*(-tz + Za + ((-ty+Ya)*(-rx) - (-tx+Xa)*(-ry)));
357 return new double[]{Xb, Yb, Zb};
358 }
359
360 public String getCacheDirectoryName() {
361 return this.toString();
362 }
363
364 /**
365 * Returns the default zoom scale in pixel per degree ({@see #NavigatableComponent#scale}))
366 */
367 public double getDefaultZoomInPPD() {
368 // this will set the map scaler to about 1000 m (in default scale, 1 pixel will be 10 meters)
369 return 10.0;
370 }
371
372 public Bounds getWorldBoundsLatLon() {
373 return utmBounds[currentGeodesic];
374 }
375
376 public String toCode() {
377 return utmEPSGs[currentGeodesic];
378 }
379
380 @Override
381 public int hashCode() {
382 return getClass().getName().hashCode()+currentGeodesic; // our only real variable
383 }
384
385 @Override public String toString() {
386 return (tr("UTM 20N (France)"));
387 }
388
389 public int getCurrentGeodesic() {
390 return currentGeodesic;
391 }
392
393 public void setupPreferencePanel(JPanel p) {
394 JComboBox prefcb = new JComboBox(utmGeodesicsNames);
395
396 prefcb.setSelectedIndex(currentGeodesic);
397 p.setLayout(new GridBagLayout());
398 p.add(new JLabel(tr("UTM20 North Geodesic system")), GBC.std().insets(5,5,0,5));
399 p.add(GBC.glue(1, 0), GBC.std().fill(GBC.HORIZONTAL));
400 p.add(prefcb, GBC.eop().fill(GBC.HORIZONTAL));
401 p.add(GBC.glue(1, 1), GBC.eol().fill(GBC.BOTH));
402 }
403
404 public Collection<String> getPreferences(JPanel p) {
405 Object prefcb = p.getComponent(2);
406 if(!(prefcb instanceof JComboBox))
407 return null;
408 currentGeodesic = ((JComboBox)prefcb).getSelectedIndex();
409 refresh7ParametersTranslation();
410 return Collections.singleton(Integer.toString(currentGeodesic+1));
411 }
412
413 public Collection<String> getPreferencesFromCode(String code) {
414 for (int i=0; i < utmEPSGs.length; i++ )
415 if (utmEPSGs[i].endsWith(code))
416 return Collections.singleton(Integer.toString(i));
417 return null;
418 }
419
420 public void setPreferences(Collection<String> args) {
421 currentGeodesic = DEFAULT_GEODESIC;
422 if (args != null) {
423 try {
424 for(String s : args)
425 {
426 currentGeodesic = Integer.parseInt(s)-1;
427 if(currentGeodesic < 0 || currentGeodesic > 2) {
428 currentGeodesic = DEFAULT_GEODESIC;
429 }
430 break;
431 }
432 } catch(NumberFormatException e) {}
433 }
434 refresh7ParametersTranslation();
435 }
436
437}
Note: See TracBrowser for help on using the repository browser.