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

Last change on this file since 3689 was 3473, checked in by bastiK, 14 years ago

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