001: /*
002: * The Unified Mapping Platform (JUMP) is an extensible, interactive GUI
003: * for visualizing and manipulating spatial features with geometry and attributes.
004: *
005: * Copyright (C) 2003 Vivid Solutions
006: *
007: * This program is free software; you can redistribute it and/or
008: * modify it under the terms of the GNU General Public License
009: * as published by the Free Software Foundation; either version 2
010: * of the License, or (at your option) any later version.
011: *
012: * This program is distributed in the hope that it will be useful,
013: * but WITHOUT ANY WARRANTY; without even the implied warranty of
014: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
015: * GNU General Public License for more details.
016: *
017: * You should have received a copy of the GNU General Public License
018: * along with this program; if not, write to the Free Software
019: * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
020: *
021: * For more information, contact:
022: *
023: * Vivid Solutions
024: * Suite #1A
025: * 2328 Government Street
026: * Victoria BC V8T 5G5
027: * Canada
028: *
029: * (250)385-6040
030: * www.vividsolutions.com
031: */
032:
033: package com.vividsolutions.jump.coordsys.impl;
034:
035: import java.util.Date;
036:
037: import com.vividsolutions.jump.coordsys.*;
038: import com.vividsolutions.jump.coordsys.Geographic;
039: import com.vividsolutions.jump.coordsys.Planar;
040:
041: /**
042:
043: * This class implements the Transverse Mercator Projection.
044:
045: * @version $Revision: 4 $
046: * @author $Author: javamap $
047:
048: *<pre>
049: * $Id: TransverseMercator.java 4 2005-06-16 15:27:48Z javamap $
050: * $Date: 2005-06-16 08:27:48 -0700 (Thu, 16 Jun 2005) $
051:
052: * $Log$
053: * Revision 1.1 2005/06/16 15:25:29 javamap
054: * *** empty log message ***
055: *
056: * Revision 1.2 2005/05/03 15:23:55 javamap
057: * *** empty log message ***
058: *
059: * Revision 1.3 2003/11/05 05:16:00 dkim
060: * Added global header; cleaned up Javadoc.
061: *
062: * Revision 1.2 2003/09/15 21:43:27 jaquino
063: * Global flag for enabling/disabling
064: * CoordinateSystemSupport
065: *
066: * Revision 1.1 2003/09/15 20:26:11 jaquino
067: * Reprojection
068: *
069: * Revision 1.2 2003/07/25 17:01:03 gkostadinov
070: * Moved classses reponsible for performing the basic projection to a new
071: * package -- base.
072: *
073: * Revision 1.1 2003/07/24 23:14:43 gkostadinov
074: * adding base projection classes
075: *
076: * Revision 1.1 2003/06/20 18:34:30 gkostadinov
077: * Entering the source code into the CVS.
078: *</pre>
079: */
080:
081: public class TransverseMercator extends Projection {
082:
083: double L0;// central meridian
084: double k0;
085:
086: public TransverseMercator() {
087: }
088:
089: /**
090: *@param centralMeridian in degrees
091: */
092: public void setParameters(double centralMeridian) {
093: L0 = centralMeridian / 180.0 * Math.PI;
094: }
095:
096: /**
097: *@param q in degrees
098: */
099: public Geographic asGeographic(Planar p, Geographic q) {
100: planarToGeographicInRadians(p, q);
101: q.lat = q.lat * 180.0 / Math.PI;
102: q.lon = q.lon * 180.0 / Math.PI;
103: return q;
104: }//scale factor
105:
106: /**
107: *@param q0 in degrees
108: */
109: public Planar asPlanar(Geographic q0, Planar p) {
110: Geographic q = new Geographic();
111: q.lat = q0.lat / 180.0 * Math.PI;
112: q.lon = q0.lon / 180.0 * Math.PI;
113: geographicInRadiansToPlanar(q, p);
114: return p;
115: }
116:
117: /**
118: *@param q in radians
119: */
120: void planarToGeographicInRadians(Planar p, Geographic q) {
121: double L1;
122: L1 = footPointLatitude(p.y);
123: double a;
124: double b;
125: a = currentSpheroid.getA();
126: b = currentSpheroid.getB();
127: double ep2;
128: double N1;
129: double M1;
130: ep2 = (a * a - b * b) / (b * b);
131: // N1 = the radius of curvature of the spheroid in the prime vertical plane
132: // at the foot point latitude
133: N1 = currentSpheroid.primeVerticalRadiusOfCurvature(L1);
134: // M1 = meridian radius of curvature at the foot point latitude
135: M1 = currentSpheroid.meridianRadiusOfCurvature(L1);
136: double n1;
137: double n12;
138: double n14;
139: double n16;
140: double n18;
141: n12 = ep2 * Math.pow(Math.cos(L1), 2.0);
142: n1 = Math.sqrt(n12);
143: n14 = n12 * n12;
144: n16 = n14 * n12;
145: n18 = n14 * n14;
146: double t1;
147: double t12;
148: double t14;
149: double t16;
150: t1 = Math.tan(L1);
151: t12 = t1 * t1;
152: t14 = t12 * t12;
153: t16 = t14 * t12;
154: double u0;
155: double u1;
156: double v1;
157: double u2;
158: double v2;
159: double u3;
160: double v3;
161: u0 = t1 * Math.pow(p.x, 2.0) / (2.0 * M1 * N1);
162: u1 = t1 * Math.pow(p.x, 4.0) / (24.0 * M1 * Math.pow(N1, 3.0));
163: u2 = t1 * Math.pow(p.x, 6.0) / (720.0 * M1 * Math.pow(N1, 5.0));
164: u3 = t1 * Math.pow(p.x, 8.0)
165: / (40320.0 * M1 * Math.pow(N1, 7.0));
166: v1 = 5.0 + 3.0 * t12 + n12 - 4.0 * n14 - 9.0 * n12 * t12;
167: v2 = 61.0 - 90.0 * t12 + 46.0 * n12 + 45.0 * t14 - 252.0 * t12
168: * n12 - 3.0 * n14 + 100.0 * n16 - 66.0 * t12 * n14
169: - 90.0 * t14 * n12 + 88.0 * n18 + 225.0 * t14 * n14
170: + 84.0 * t12 * n16 - 192.0 * t12 * n18;
171: v3 = 1385.0 + 3633.0 * t12 + 4095.0 * t14 + 1575.0 * t16;
172:
173: q.lat = L1 - u0 + u1 * v1 - u2 * v2 + u3 * v3;
174: double XdN1;
175: XdN1 = p.x / N1;
176: u0 = XdN1;
177: u1 = Math.pow(XdN1, 3.0) / 6.0;
178: u2 = Math.pow(XdN1, 5.0) / 120.0;
179: u3 = Math.pow(XdN1, 7.0) / 5040.0;
180: v1 = 1.0 + 2.0 * t12 + n12;
181: v2 = 5.0 + 6.0 * n12 + 28.0 * t12 - 3.0 * n14 + 8.0 * t12 * n12
182: + 24.0 * t14 - 4.0 * n16 + 4.0 * t12 * n14 + 24.0 * t12
183: * n16;
184: v3 = 61.0 + 662.0 * t12 + 1320.0 * t14 + 720.0 * t16;
185: q.lon = 1.0 / Math.cos(L1) * (u0 - u1 * v1 + u2 * v2 - u3 * v3)
186: + L0;
187: }
188:
189: private MeridianArcLength S = new MeridianArcLength();
190:
191: /**
192: *@param q in radians
193: */
194: void geographicInRadiansToPlanar(Geographic q, Planar p) {
195: double a;
196: double b;
197: a = currentSpheroid.getA();
198: b = currentSpheroid.getB();
199: double ep2;
200: double N;
201: // ep2 = the second eccentricity squared.
202: ep2 = (a * a - b * b) / (b * b);
203: // N = the radius of curvature of the spheroid in the prime vertical plane
204: N = currentSpheroid.primeVerticalRadiusOfCurvature(q.lat);
205: double n;
206: double n2;
207: double n4;
208: double n6;
209: double n8;
210: n2 = ep2 * Math.pow(Math.cos(q.lat), 2.0);
211: n = Math.sqrt(n2);
212: n4 = n2 * n2;
213: n6 = n4 * n2;
214: n8 = n4 * n4;
215: double t;
216: double t2;
217: double t4;
218: double t6;
219: t = Math.tan(q.lat);
220: t2 = t * t;
221: t4 = t2 * t2;
222: t6 = t4 * t2;
223: S.compute(currentSpheroid, q.lat, 0);
224: double cosLat;
225: double sinLat;
226: cosLat = Math.cos(q.lat);
227: sinLat = Math.sin(q.lat);
228: double L;
229: double L2;
230: double L3;
231: double L4;
232: double L5;
233: double L6;
234: double L7;
235: double L8;
236: L = q.lon - L0;// 'L' for lambda (longitude) - must be in radians
237: L2 = L * L;
238: L3 = L2 * L;
239: L4 = L2 * L2;
240: L5 = L4 * L;
241: L6 = L4 * L2;
242: L7 = L5 * L2;
243: L8 = L4 * L4;
244: double u0;
245: double u1;
246: double v1;
247: double u2;
248: double v2;
249: double u3;
250: double v3;
251: u0 = L * cosLat;
252: u1 = L3 * Math.pow(cosLat, 3.0) / 6.0;
253: u2 = L5 * Math.pow(cosLat, 5.0) / 120.0;
254: u3 = L7 * Math.pow(cosLat, 7.0) / 5040.0;
255: v1 = 1.0 - t2 + n2;
256: v2 = 5.0 - 18.0 * t2 + t4 + 14.0 * n2 - 58.0 * t2 * n2 + 13.0
257: * n4 + 4.0 * n6 - 64.0 * n4 * t2 - 24.0 * n6 * t2;
258: v3 = 61.0 - 479.0 * t2 + 179.0 * t4 - t6;
259: p.x = u0 + u1 * v1 + u2 * v2 + u3 * v3;
260:
261: u0 = L2 / 2.0 * sinLat * cosLat;
262: u1 = L4 / 24.0 * sinLat * Math.pow(cosLat, 3.0);
263: u2 = L6 / 720.0 * sinLat * Math.pow(cosLat, 5.0);
264: u3 = L8 / 40320.0 * sinLat * Math.pow(cosLat, 7.0);
265: v1 = 5.0 - t2 + 9.0 * n2 + 4.0 * n4;
266: v2 = 61.0 - 58.0 * t2 + t4 + 270.0 * n2 - 330.0 * t2 * n2
267: + 445.0 * n4 + 324.0 * n6 - 680.0 * n4 * t2 + 88.0 * n8
268: - 600.0 * n6 * t2 - 192.0 * n8 * t2;
269: v3 = 1385.0 - 311.0 * t2 + 543.0 * t4 - t6;
270: p.y = S.s / N + u0 + u1 * v1 + u2 * v2 + u3 * v3;
271:
272: p.x = N * p.x;
273: p.y = N * p.y;
274: }
275:
276: private double footPointLatitude(double y) {
277: // returns the footpoint Latitude given the y coordinate
278: double newlat;
279: // returns the footpoint Latitude given the y coordinate
280: double Lat1;
281: // returns the footpoint Latitude given the y coordinate
282: double flat;
283: // returns the footpoint Latitude given the y coordinate
284: double dflat;
285: double a;
286: a = currentSpheroid.getA();
287: newlat = y / a;
288: int i = 0;
289: do {
290: Lat1 = newlat;
291: i++;
292: if (i == 100) {
293: //Prevent infinite loop. I observed that a typical number of iterations is 5. [Jon Aquino]
294: break;
295: }
296: S.compute(currentSpheroid, Lat1, 0);
297: flat = S.s - y;
298: dflat = a
299: * (S.a0 - 2.0 * S.a2 * Math.cos(2.0 * Lat1) + 4.0
300: * S.a4 * Math.cos(4.0 * Lat1) - 6.0 * S.a6
301: * Math.cos(6.0 * Lat1) + 8.0 * S.a8
302: * Math.cos(8.0 * Lat1));
303: newlat = Lat1 - flat / dflat;
304: //Increased tolerance from 1E-16 to 1E-15. 1E-16 was causing an infinite loop.
305: //JA 6 Nov 2001.
306: } while (Math.abs(newlat - Lat1) > 1.0e-015);
307: Lat1 = newlat;
308: return Lat1;
309: }// END - footPointLatitude
310:
311: }
|