KickJava   Java API By Example, From Geeks To Geeks.

Java > Open Source Codes > org > jscience > geography > coordinates > UTM


1 /*
2  * JScience - Java(TM) Tools and Libraries for the Advancement of Sciences.
3  * Copyright (C) 2006 - JScience (http://jscience.org/)
4  * All rights reserved.
5  *
6  * Permission to use, copy, modify, and distribute this software is
7  * freely granted, provided that this notice is preserved.
8  */
9 package org.jscience.geography.coordinates;
10
11 import static org.jscience.geography.coordinates.crs.ReferenceEllipsoid.WGS84;
12
13 import javax.measure.converters.UnitConverter;
14 import javax.measure.quantities.*;
15 import javax.measure.units.*;
16
17 import org.jscience.geography.coordinates.crs.*;
18 import org.opengis.referencing.cs.CoordinateSystem;
19
20 /**
21  * This class represents the {@link ProjectedCRS projected}
22  * Universal Transverse Mercator (UTM) coordinates onto the WGS84 ellipsoid.
23  *
24  * <p>The UTM system is limited to values between -80 and +84 degrees latitude.
25  * Values beyond these limits (i.e., the polar regions) are projected
26  * using the Universal Polar Stereographic (UPS) projection. Although
27  * mathematically distinct, the two projections are represented identically.
28  * This class returns correct results for both UTM and UPS projections.
29  *
30  * The conversion routines for this class were derived from formulas described
31  * in the
32  * <a HREF="http://earth-info.nga.mil/GandG/publications/tm8358.2/TM8358_2.pdf">
33  * Defense Mapping Agency Technical Manual 8358.2.</a>
34  *
35  * @author Paul D. Anderson
36  * @version 3.0, February 25, 2006
37  * @see <a HREF="http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system">
38  * Wikipedia: Universal Transverse Mercator Coordinate System</a>
39  */
40 public class UTM extends Coordinates<ProjectedCRS> {
41
42     /**
43      * The UTM scale factor. This the exact scale factor only on a pair of
44      * lines lying either side of the central meridian, but the effect is to
45      * reduce overall distortion within the UTM zone to less than one part per
46      * thousand.
47      */
48     public static final double UTM_SCALE_FACTOR = 0.9996;
49
50     /**
51      * The UTM "false easting" value. This quantity is added to the true
52      * easting to avoid using negative numbers in the coordinates.
53      */
54     public static Quantity<Length> UTM_FALSE_EASTING = new Scalar<Length>(
55             500000.0, SI.METER);
56
57     /**
58      * The UTM "false northing" value. This quantity is added to the true
59      * northing for coordinates <b>in the southern hemisphere only</b>
60      * to avoid using negative numbers in the coordinates.
61      */
62     public static Quantity<Length> UTM_FALSE_NORTHING = new Scalar<Length>(
63             10000000.0, SI.METER);
64
65     /**
66      * The northern limit of the UTM grid. Beyond this limit the distortion
67      * introduced by the transverse Mercator projection is impractically
68      * large, and the UPS grid is used instead.
69      */
70     public static final Quantity<Angle> UTM_NORTHERN_LIMIT = new Scalar<Angle>(
71             84.0, NonSI.DEGREE_ANGLE);
72
73     /**
74      * The southern limit of the UTM grid. Beyond this limit the distortion
75      * introduced by the transverse Mercator projection is impractically
76      * large, and the UPS grid is used instead.
77      */
78     public static final Quantity<Angle> UTM_SOUTHERN_LIMIT = new Scalar<Angle>(
79             -80.0, NonSI.DEGREE_ANGLE);
80
81     /**
82      * The UPS scale factor.
83      */
84     public static final double UPS_SCALE_FACTOR = 0.994;
85
86     /**
87      * The UPS "false easting" value. This quantity is added to the true
88      * easting to avoid using negative numbers in the coordinates.
89      */
90     public static Quantity<Length> UPS_FALSE_EASTING = new Scalar<Length>(
91             2000000.0, SI.METER);
92
93     /**
94      * The UPS "false northing" value. This quantity is added to the true
95      * northing to avoid using negative numbers in the coordinates.
96      * The UPS system, unlike the UTM system, always includes the false northing.
97      */
98     public static Quantity<Length> UPS_FALSE_NORTHING = new Scalar<Length>(
99             2000000.0, SI.METER);
100
101     /*
102      * NOTE: The calculations in this class use power series expansions.
103      * The naming convention is to include the power in the name
104      * of the term, so that the square of K0 is 'K02', the cube
105      * is 'K03', etc.
106      */
107     private static final double K0 = UTM_SCALE_FACTOR;
108
109     private static final double K02 = K0 * K0;
110
111     private static final double K03 = K02 * K0;
112
113     private static final double K04 = K03 * K0;
114
115     private static final double K05 = K04 * K0;
116
117     private static final double K06 = K05 * K0;
118
119     private static final double K07 = K06 * K0;
120
121     private static final double K08 = K07 * K0;
122
123     /**
124      * Holds the coordinate reference system for all instances of this class.
125      */
126     public static final ProjectedCRS<UTM> CRS = new ProjectedCRS<UTM>() {
127
128         private final double NORTHERN_LIMIT_IN_DEGREES = UTM_NORTHERN_LIMIT
129                 .doubleValue(NonSI.DEGREE_ANGLE);
130
131         private final double SOUTHERN_LIMIT_IN_DEGREES = UTM_SOUTHERN_LIMIT
132                 .doubleValue(NonSI.DEGREE_ANGLE);
133
134         @Override JavaDoc
135         protected UTM coordinatesOf(AbsolutePosition position) {
136             LatLong latLong = LatLong.valueOf(position.latitudeWGS84
137                     .doubleValue(SI.RADIAN), position.longitudeWGS84
138                     .doubleValue(SI.RADIAN), SI.RADIAN);
139             // UTM or UPS?
140 final double latitude = position.latitudeWGS84
141                     .doubleValue(NonSI.DEGREE_ANGLE);
142             if (latitude > SOUTHERN_LIMIT_IN_DEGREES
143                     && latitude < NORTHERN_LIMIT_IN_DEGREES) {
144                 return latLongToUtm(latLong, WGS84);
145             } else {
146                 return latLongToUps(latLong, WGS84);
147             }
148         }
149
150         @Override JavaDoc
151         protected AbsolutePosition positionOf(UTM coordinates,
152                 AbsolutePosition position) {
153             final LatLong latLong;
154             if (coordinates.latitudeZone() < 'C'
155                     || coordinates.latitudeZone() > 'X') {
156                 latLong = upsToLatLong(coordinates, WGS84);
157             } else {
158                 latLong = utmToLatLong(coordinates, WGS84);
159             }
160             position.latitudeWGS84 = new Scalar<Angle>(latLong
161                     .latitudeValue(SI.RADIAN), SI.RADIAN);
162             position.longitudeWGS84 = new Scalar<Angle>(latLong
163                     .longitudeValue(SI.RADIAN), SI.RADIAN);
164             return position;
165         }
166
167         @Override JavaDoc
168         public CoordinateSystem getCoordinateSystem() {
169             return ProjectedCRS.EASTING_NORTHING_CS;
170         }
171     };
172
173     /**
174      * Holds the longitude zone identifier.
175      */
176     private int _longitudeZone;
177
178     /**
179      * Holds the latitude zone identifier.
180      */
181     private char _latitudeZone;
182
183     /**
184      * Holds the easting in meters.
185      */
186     private double _eastingInMeter;
187
188     /**
189      * Holds the northing in meters.
190      */
191     private double _northingInMeter;
192
193     /**
194      * Returns the projected UTM position corresponding to the specified
195      * coordinates.
196      *
197      * @param longitudeZone the longitude zone number.
198      * @param latitudeZone the longitude zone character.
199      * @param easting the easting value stated in the specified unit.
200      * @param northing the northing value stated in the specified unit.
201      * @param unit the easting/northing length unit.
202      * @return the corresponding surface position.
203      */
204     public static UTM valueOf(int longitudeZone, char latitudeZone,
205             double easting, double northing, Unit<Length> unit) {
206         return new UTM(longitudeZone, latitudeZone, easting, northing, unit);
207     }
208
209     /**
210      * Creates the projected UTM position corresponding to the specified
211      * coordinates.
212      *
213      * @param longitudeZone the longitude zone number.
214      * @param latitudeZone the longitude zone character.
215      * @param easting the easting value stated in the specified unit.
216      * @param northing the northing value stated in the specified unit.
217      * @param unit the easting/northing length unit.
218      */
219     public UTM(int longitudeZone, char latitudeZone, double easting,
220             double northing, Unit<Length> unit) {
221         _longitudeZone = longitudeZone;
222         _latitudeZone = latitudeZone;
223         if (unit == SI.METER) {
224             _eastingInMeter = easting;
225             _northingInMeter = northing;
226         } else {
227             UnitConverter toMeter = unit.getConverterTo(SI.METER);
228             _eastingInMeter = toMeter.convert(easting);
229             _northingInMeter = toMeter.convert(northing);
230         }
231     }
232
233     /**
234      * Returns the longitude zone identifier.
235      *
236      * @return the longitude zone number.
237      */
238     public final int longitudeZone() {
239         return _longitudeZone;
240     }
241
242     /**
243      * Returns the latitude zone identifier.
244      *
245      * @return the latitude zone character.
246      */
247     public final char latitudeZone() {
248         return _latitudeZone;
249     }
250
251     /**
252      * Returns the projected distance of the position from the central meridian.
253      *
254      * @param unit the length unit of the easting to return.
255      * @return the easting stated in the specified unit.
256      */
257     public final double eastingValue(Unit<Length> unit) {
258         return unit.equals(SI.METER) ? _eastingInMeter : SI.METER
259                 .getConverterTo(unit).convert(_eastingInMeter);
260     }
261
262     /**
263      * Returns the projected distance of the point from the equator.
264      *
265      * @param unit the length unit of the northing to return.
266      * @return the northing stated in the specified unit.
267      */
268     public final double northingValue(Unit<Length> unit) {
269         return unit.equals(SI.METER) ? _northingInMeter : SI.METER
270                 .getConverterTo(unit).convert(_northingInMeter);
271     }
272
273     @Override JavaDoc
274     public ProjectedCRS<UTM> getCoordinateReferenceSystem() {
275         return CRS;
276     }
277
278     // OpenGIS Interface.
279 public int getDimension() {
280         return 2;
281     }
282
283     // OpenGIS Interface.
284 public double getOrdinate(int dimension) throws IndexOutOfBoundsException JavaDoc {
285         if (dimension == 0) {
286             Unit u = ProjectedCRS.EASTING_NORTHING_CS.getAxis(0).getUnit();
287             return SI.METER.getConverterTo(u).convert(_eastingInMeter);
288         } else if (dimension == 1) {
289             Unit u = ProjectedCRS.EASTING_NORTHING_CS.getAxis(1).getUnit();
290             return SI.METER.getConverterTo(u).convert(_northingInMeter);
291         } else {
292             throw new IndexOutOfBoundsException JavaDoc();
293         }
294     }
295
296     /**
297      * Returns true if the position indicated by the coordinates is
298      * north of the northern limit of the UTM grid (84 degrees).
299      *
300      * @param latLong The coordinates.
301      * @return True if the latitude is greater than 84 degrees.
302      */
303     public static boolean isNorthPolar(final LatLong latLong) {
304         return latLong.latitudeValue(NonSI.DEGREE_ANGLE) > 84.0;
305     }
306
307     /**
308      * Returns true if the position indicated by the coordinates is
309      * south of the southern limit of the UTM grid (-80 degrees).
310      *
311      * @param latLong The coordinates.
312      * @return True if the latitude is less than -80 degrees.
313      */
314     public static boolean isSouthPolar(final LatLong latLong) {
315         return latLong.latitudeValue(NonSI.DEGREE_ANGLE) < -80.0;
316     }
317
318     /**
319      * Returns the UTM/UPS latitude zone identifier for the specified coordinates.
320      *
321      * @param latLong The coordinates.
322      * @return the latitude zone character.
323      */
324     public static char getLatitudeZone(final LatLong latLong) {
325         if (isNorthPolar(latLong)) {
326             if (latLong.longitudeValue(SI.RADIAN) < 0) {
327                 return 'Y';
328             } else {
329                 return 'Z';
330             }
331         }
332         if (isSouthPolar(latLong)) {
333             if (latLong.longitudeValue(SI.RADIAN) < 0) {
334                 return 'A';
335             } else {
336                 return 'B';
337             }
338         }
339         final int degreesLatitude = (int) latLong
340                 .latitudeValue(NonSI.DEGREE_ANGLE);
341         char zone = (char) ((degreesLatitude + 80) / 8 + 'C');
342         if (zone > 'H') {
343             zone++;
344         }
345         if (zone > 'N') {
346             zone++;
347         }
348         if (zone > 'X') {
349             zone = 'X';
350         }
351         return zone;
352     }
353
354     /**
355      * Returns the UTM/UPS longitude zone number for the specified
356      * coordinates.
357      *
358      * @param latLong The coordinates.
359      * @return the longitude zone number.
360      */
361     public static int getLongitudeZone(LatLong latLong) {
362
363         final double degreesLongitude = latLong
364                 .longitudeValue(NonSI.DEGREE_ANGLE);
365
366         // UPS longitude zones
367 if (isNorthPolar(latLong) || isSouthPolar(latLong)) {
368             if (degreesLongitude < 0.0) {
369                 return 30;
370             } else {
371                 return 31;
372             }
373         }
374
375         final char latitudeZone = getLatitudeZone(latLong);
376         // X latitude exceptions
377 if (latitudeZone == 'X' && degreesLongitude > 0.0
378                 && degreesLongitude < 42.0) {
379             if (degreesLongitude < 9.0) {
380                 return 31;
381             }
382             if (degreesLongitude < 21.0) {
383                 return 33;
384             }
385             if (degreesLongitude < 33.0) {
386                 return 35;
387             } else {
388                 return 37;
389             }
390         }
391         // V latitude exceptions
392 if (latitudeZone == 'V' && degreesLongitude > 0.0
393                 && degreesLongitude < 12.0) {
394             if (degreesLongitude < 3.0) {
395                 return 31;
396             } else {
397                 return 32;
398             }
399         }
400
401         return ((int) degreesLongitude + 180) / 6 + 1;
402     }
403
404     /**
405      * Returns the central meridian (in radians) for the specified
406      * UTM/UPS zone.
407      * @param longitudeZone The UTM/UPS longitude zone number.
408      * @param latitudeZone The UTM/UPS latitude zone character.
409      * @return The central meridian for the specified zone.
410      */
411     public static double getCentralMeridian(int longitudeZone, char latitudeZone) {
412         // polar zones
413 if (latitudeZone < 'C' || latitudeZone > 'X') {
414             return 0.0;
415         }
416         // X latitude zone exceptions
417 if (latitudeZone == 'X' && longitudeZone > 31 && longitudeZone <= 37) {
418             return Math.toRadians((longitudeZone - 1) * 6 - 180 + 4.5);
419         }
420         // V latitude zone exceptions
421 if (longitudeZone == 'V') {
422             if (latitudeZone == 31) {
423                 return Math.toRadians(1.5);
424             } else if (latitudeZone == 32) {
425                 return Math.toRadians(7.5);
426             }
427         }
428         return Math.toRadians((longitudeZone - 1) * 6 - 180 + 3);
429     }
430
431     /**
432      * Converts latitude/longitude coordinates to UTM coordinates based on
433      * the specified reference ellipsoid.
434      *
435      * @param latLong The latitude/longitude coordinates.
436      * @param ellipsoid The reference ellipsoid.
437      * @return The UTM coordinates.
438      */
439     public static UTM latLongToUtm(LatLong latLong, ReferenceEllipsoid ellipsoid) {
440         final char latitudeZone = getLatitudeZone(latLong);
441         final int longitudeZone = getLongitudeZone(latLong);
442
443         final double phi = latLong.latitudeValue(SI.RADIAN);
444
445         final double cosPhi = Math.cos(phi);
446         final double cos2Phi = cosPhi * cosPhi;
447         final double cos3Phi = cos2Phi * cosPhi;
448         final double cos4Phi = cos3Phi * cosPhi;
449         final double cos5Phi = cos4Phi * cosPhi;
450         final double cos6Phi = cos5Phi * cosPhi;
451         final double cos7Phi = cos6Phi * cosPhi;
452         final double cos8Phi = cos7Phi * cosPhi;
453
454         final double tanPhi = Math.tan(phi);
455         final double tan2Phi = tanPhi * tanPhi;
456         final double tan4Phi = tan2Phi * tan2Phi;
457         final double tan6Phi = tan4Phi * tan2Phi;
458
459         final double eb2 = ellipsoid.getSecondEccentricitySquared();
460         final double eb4 = eb2 * eb2;
461         final double eb6 = eb4 * eb2;
462         final double eb8 = eb6 * eb2;
463
464         final double e2c2 = eb2 * cos2Phi;
465         final double e4c4 = eb4 * cos4Phi;
466         final double e6c6 = eb6 * cos6Phi;
467         final double e8c8 = eb8 * cos8Phi;
468
469         final double t2e2c2 = tan2Phi * e2c2;
470         final double t2e4c4 = tan2Phi * e4c4;
471         final double t2e6c6 = tan2Phi * e6c6;
472         final double t2e8c8 = tan2Phi * e8c8;
473
474         final double nu = ellipsoid.verticalRadiusOfCurvature(phi);
475         final double kn1 = K0 * nu * Math.sin(phi);
476         final double t1 = K0 * ellipsoid.meridionalArc(phi);
477         final double t2 = kn1 * cosPhi / 2.0;
478         final double t3 = (kn1 * cos3Phi / 24.0)
479                 * (5.0 - tan2Phi + 9.0 * e2c2 + 4.0 * e4c4);
480         final double t4 = (kn1 * cos5Phi / 720.0)
481                 * (61.0 - 58.0 * tan2Phi + tan4Phi + 270.0 * e2c2 - 330.0
482                         * t2e2c2 + 445.0 * e4c4 - 680.0 * t2e4c4 + 324.0 * e6c6
483                         - 600.0 * t2e6c6 + 88.0 * e8c8 - 192.0 * t2e8c8);
484         final double t5 = (kn1 * cos7Phi / 40320.0)
485                 * (1385.0 - 3111.0 * tan2Phi + 543.0 * tan4Phi - tan6Phi);
486
487         final double kn2 = K0 * nu;
488         final double t6 = kn2 * cosPhi;
489         final double t7 = (kn2 * cos3Phi / 6.0) * (1.0 - tan2Phi + e2c2);
490         final double t8 = (kn2 * cos5Phi / 120.0)
491                 * (5.0 - 18.0 * tan2Phi + tan4Phi + 14.0 * e2c2 - 58.0 * t2e2c2
492                         + 13.0 * e4c4 - 64.0 * t2e4c4 + 4.0 * e6c6 - 24.0 * t2e6c6);
493         final double t9 = (kn2 * cos7Phi / 50.40)
494                 * (61.0 - 479.0 * tan2Phi + 179.0 * tan4Phi - tan6Phi);
495
496         final double lambda = latLong.longitudeValue(SI.RADIAN);
497         final double lambda0 = getCentralMeridian(longitudeZone, latitudeZone);
498         final double dL = lambda - lambda0;
499         final double dL2 = dL * dL;
500         final double dL3 = dL2 * dL;
501         final double dL4 = dL3 * dL;
502         final double dL5 = dL4 * dL;
503         final double dL6 = dL5 * dL;
504         final double dL7 = dL6 * dL;
505         final double dL8 = dL7 * dL;
506
507         final double falseNorthing;
508         if ((phi < 0.0)) {
509             // southern hemisphere -- add false northing
510 falseNorthing = UTM_FALSE_NORTHING.doubleValue(SI.METER);
511         } else {
512             // northern hemisphere -- no false northing
513 falseNorthing = 0.0;
514         }
515         final double falseEasting = UTM_FALSE_EASTING.doubleValue(SI.METER);
516         final double northing = falseNorthing + t1 + dL2 * t2 + dL4 * t3 + dL6
517                 * t4 + dL8 * t5;
518         final double easting = falseEasting + dL * t6 + dL3 * t7 + dL5 * t8
519                 + dL7 * t9;
520
521         return UTM.valueOf(longitudeZone, latitudeZone, easting, northing,
522                 SI.METER);
523     }
524
525     /**
526      * Converts latitude/longitude coordinates to UPS coordinates based on
527      * the specified reference ellipsoid.
528      *
529      * @param latLong The latitude/longitude coordinates.
530      * @param ellipsoid The reference ellipsoid.
531      * @return The UPS coordinates.
532      */
533     public static UTM latLongToUps(LatLong latLong, ReferenceEllipsoid ellipsoid) {
534
535         final char latitudeZone = getLatitudeZone(latLong);
536         final int longitudeZone = getLongitudeZone(latLong);
537
538         final double latitude = latLong.latitudeValue(SI.RADIAN);
539         final double sign = Math.signum(latitude);
540         final double phi = Math.abs(latitude);
541         final double lambda = latLong.longitudeValue(SI.RADIAN);
542
543         final double a = ellipsoid.getSemimajorAxis().doubleValue(SI.METER);
544         final double e = ellipsoid.getEccentricity();
545         final double e2 = ellipsoid.getEccentricitySquared();
546
547         final double c0 = ((2.0 * a) / Math.sqrt(1.0 - e2))
548                 * Math.pow((1.0 - e) / (1.0 + e), e / 2.0);
549         final double eSinPhi = e * Math.sin(phi);
550         final double tz = Math.pow((1 + eSinPhi) / (1 - eSinPhi), e / 2.0)
551                 * Math.tan(Math.PI / 4.0 - phi / 2.0);
552         final double radius = UPS_SCALE_FACTOR * c0 * tz;
553         final double falseNorthing = UPS_FALSE_NORTHING.doubleValue(SI.METER);
554         final double northing;
555         if (sign > 0) {
556             northing = falseNorthing - radius * Math.cos(lambda);
557         } else {
558             northing = falseNorthing + radius * Math.cos(lambda);
559         }
560         final double falseEasting = UPS_FALSE_EASTING.doubleValue(SI.METER);
561         final double easting = falseEasting + radius * Math.sin(lambda);
562
563         return UTM.valueOf(longitudeZone, latitudeZone, easting, northing,
564                 SI.METER);
565     }
566
567     /**
568      * Converts the UTM coordinates to latitude/longitude coordinates,
569      * based on the specified reference ellipsoid.
570      *
571      * @param utm The UTM coordinates.
572      * @param ellipsoid The reference ellipsoid.
573      * @return The latitude/longitude coordinates.
574      */
575     public static LatLong utmToLatLong(UTM utm, ReferenceEllipsoid ellipsoid) {
576         final double northing;
577         if ((utm.latitudeZone() < 'N')) {
578             // southern hemisphere
579 northing = utm._northingInMeter
580                     - UTM_FALSE_NORTHING.doubleValue(SI.METER);
581         } else {
582             // northern hemisphere
583 northing = utm._northingInMeter;
584         }
585
586         // footpoint latitude
587 final double arc0 = northing / K0;
588         double rho = ellipsoid.meridionalRadiusOfCurvature(0.0);
589         double phi = arc0 / rho;
590         for (int i = 0; i < 5; i++) {
591             double arc = ellipsoid.meridionalArc(phi);
592             rho = ellipsoid.meridionalRadiusOfCurvature(phi);
593             double diff = (arc0 - arc) / rho;
594             if (Math.abs(diff) < Math.ulp(phi)) {
595                 break;
596             }
597             phi += diff;
598         }
599
600         final double cosPhi = Math.cos(phi);
601         final double cos2Phi = cosPhi * cosPhi;
602         final double cos3Phi = cos2Phi * cosPhi;
603         final double cos4Phi = cos3Phi * cosPhi;
604         final double cos5Phi = cos4Phi * cosPhi;
605         final double cos6Phi = cos5Phi * cosPhi;
606         final double cos7Phi = cos6Phi * cosPhi;
607         final double cos8Phi = cos7Phi * cosPhi;
608
609         final double tanPhi = Math.tan(phi);
610         final double tan2Phi = tanPhi * tanPhi;
611         final double tan4Phi = tan2Phi * tan2Phi;
612         final double tan6Phi = tan4Phi * tan2Phi;
613
614         final double eb2 = ellipsoid.getSecondEccentricitySquared();
615         final double eb4 = eb2 * eb2;
616         final double eb6 = eb4 * eb2;
617         final double eb8 = eb6 * eb2;
618         final double e2c2 = eb2 * cos2Phi;
619         final double e4c4 = eb4 * cos4Phi;
620         final double e6c6 = eb6 * cos6Phi;
621         final double e8c8 = eb8 * cos8Phi;
622
623         final double t2e2c2 = tan2Phi * e2c2;
624         final double t2e4c4 = tan2Phi * e4c4;
625         final double t2e6c6 = tan2Phi * e6c6;
626         final double t2e8c8 = tan2Phi * e8c8;
627         final double t4e2c2 = tan4Phi * e2c2;
628         final double t4e4c4 = tan4Phi * e4c4;
629
630         final double nu = ellipsoid.verticalRadiusOfCurvature(phi);
631         final double nu2 = nu * nu;
632         final double nu3 = nu2 * nu;
633         final double nu5 = nu3 * nu2;
634         final double nu7 = nu5 * nu2;
635
636         final double lambda0 = getCentralMeridian(utm.longitudeZone(), utm
637                 .latitudeZone());
638         final double dE = utm._eastingInMeter
639                 - UTM_FALSE_EASTING.doubleValue(SI.METER);
640         final double dE2 = dE * dE;
641         final double dE3 = dE2 * dE;
642         final double dE4 = dE3 * dE;
643         final double dE5 = dE4 * dE;
644         final double dE6 = dE5 * dE;
645         final double dE7 = dE6 * dE;
646         final double dE8 = dE7 * dE;
647
648         final double t10 = tanPhi / (2.0 * rho * nu * K02);
649         final double t11 = tanPhi / (24.0 * rho * nu3 * K04)
650                 * (5.0 + 3.0 * tan2Phi + e2c2 - 9.0 * t2e2c2 - 4.0 * e4c4);
651         final double t12 = tanPhi
652                 / (720.0 * rho * nu5 * K06)
653                 * (61.0 + 90.0 * tan2Phi + 45.0 * tan4Phi + 46.0 * e2c2 - 252.0
654                         * t2e2c2 - 90.0 * t4e2c2 - 3.0 * e4c4 - 66.0 * t2e4c4
655                         + 225.0 * t4e4c4 + 100.0 * e6c6 + 84.0 * t2e6c6 + 88.0
656                         * e8c8 - 192.0 * t2e8c8);
657         final double t13 = tanPhi
658                 / (40320.0 * rho * nu7 * K08)
659                 * (1385.0 + 3633.0 * tan2Phi + 4095.0 * tan4Phi + 1575.0 * tan6Phi);
660         final double t14 = 1.0 / (cosPhi * nu * K0);
661         final double t15 = 1.0 / (6.0 * cosPhi * nu3 * K03)
662                 * (1.0 + 2.0 * tan2Phi + e2c2);
663         final double t16 = 1.0
664                 / (120.0 * cosPhi * nu5 * K05)
665                 * (5.0 + 28.0 * tan2Phi + 24.0 * tan4Phi + 6.0 * e2c2 + 8.0
666                         * t2e2c2 - 3.0 * e4c4 + 4.0 * t2e4c4 - 4.0 * e6c6 + 24.0 * t2e6c6);
667         final double t17 = 1.0 / (5040.0 * cosPhi * nu7 * K07)
668                 * (61.0 + 662.0 * tan2Phi + 1320.0 * tan4Phi + 720.0 * tan6Phi);
669
670         final double latitude = phi - dE2 * t10 + dE4 * t11 - dE6 * t12 + dE8
671                 * t13;
672         final double longitude = lambda0 + dE * t14 - dE3 * t15 + dE5 * t16
673                 - dE7 * t17;
674         return LatLong.valueOf(latitude, longitude, SI.RADIAN);
675     }
676
677     /**
678      * Converts the UPS coordinates to latitude/longitude coordinates,
679      * based on the specified reference ellipsoid.
680      *
681      * @param ups The UPS coordinates.
682      * @param ellipsoid The reference ellipsoid.
683      * @return The latitude/longitude coordinates.
684      */
685     public static LatLong upsToLatLong(UTM ups, ReferenceEllipsoid ellipsoid) {
686         final boolean northernHemisphere = ups.latitudeZone() > 'N';
687         final double dN = ups.northingValue(SI.METER)
688                 - UPS_FALSE_NORTHING.doubleValue(SI.METER);
689         final double dE = ups.eastingValue(SI.METER)
690                 - UPS_FALSE_EASTING.doubleValue(SI.METER);
691         // check for zeroes (the poles)
692 if (dE == 0.0 && dN == 0.0) {
693             if (northernHemisphere) {
694                 return LatLong.valueOf(90.0, 0.0, NonSI.DEGREE_ANGLE);
695             } else {
696                 return LatLong.valueOf(-90.0, 0.0, NonSI.DEGREE_ANGLE);
697             }
698         }
699         // compute longitude
700 final double longitude;
701         if (northernHemisphere) {
702             longitude = Math.atan2(dE, -dN);
703         } else {
704             longitude = Math.atan2(dE, dN);
705         }
706
707         // compute latitude
708 final double a = ellipsoid.getSemimajorAxis().doubleValue(SI.METER);
709         final double e = ellipsoid.getEccentricity();
710         final double e2 = ellipsoid.getEccentricitySquared();
711         final double e4 = e2 * e2;
712         final double e6 = e4 * e2;
713         final double e8 = e6 * e2;
714         final double aBar = e2 / 2.0 + 5.0 * e4 / 24.0 + e6 / 12.0 + 13 * e8
715                 / 360.0;
716         final double bBar = 7.0 * e4 / 48.0 + 29.0 * e6 / 240.0 + 811.0 * e8
717                 / 11520.0;
718         final double cBar = 7.0 * e6 / 120.0 + 81.0 * e8 / 1120.0;
719         final double dBar = 4279 * e8 / 161280.0;
720         final double c0 = ((2.0 * a) / Math.sqrt(1.0 - e2))
721                 * Math.pow((1.0 - e) / (1.0 + e), e / 2.0);
722         final double r;
723         if (dE == 0.0) {
724             r = dN;
725         } else if (dN == 0.0) {
726             r = dE;
727         } else if (dN < dE) {
728             r = dE / Math.sin(longitude);
729         } else {
730             r = dN / Math.cos(longitude);
731         }
732         final double radius = Math.abs(r);
733
734         final double chi = (Math.PI / 2.0) - 2.0
735                 * Math.atan2(radius, UPS_SCALE_FACTOR * c0);
736         final double phi = chi + aBar * Math.sin(2.0 * chi) + bBar
737                 * Math.sin(4.0 * chi) + cBar * Math.sin(6.0 * chi) + dBar
738                 * Math.sin(8.0 * chi);
739         final double latitude;
740         if (northernHemisphere) {
741             latitude = phi;
742         } else {
743             latitude = -phi;
744         }
745         return LatLong.valueOf(latitude, longitude, SI.RADIAN);
746     }
747
748 }
Popular Tags