source: josm/trunk/src/org/openstreetmap/josm/data/projection/UTM_France_DOM.java@ 4153

Last change on this file since 4153 was 3779, checked in by Upliner, 13 years ago

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