Ticket #4641: UTM_France_DOM.java

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