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