root / branches / F2 / libraries / libJCRS / src / org / geotools / referencing / operation / projection / TransverseMercator.java @ 12186
History | View | Annotate | Download (29.8 KB)
1 |
/*
|
---|---|
2 |
* GeoTools - OpenSource mapping toolkit
|
3 |
* http://geotools.org
|
4 |
*
|
5 |
* (C) 2003-2006, Geotools Project Managment Committee (PMC)
|
6 |
* (C) 2002, Centre for Computational Geography
|
7 |
* (C) 2001, Institut de Recherche pour le Dveloppement
|
8 |
* (C) 2000, Frank Warmerdam
|
9 |
* (C) 1999, Fisheries and Oceans Canada
|
10 |
*
|
11 |
* This library is free software; you can redistribute it and/or
|
12 |
* modify it under the terms of the GNU Lesser General Public
|
13 |
* License as published by the Free Software Foundation; either
|
14 |
* version 2.1 of the License, or (at your option) any later version.
|
15 |
*
|
16 |
* This library is distributed in the hope that it will be useful,
|
17 |
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
18 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
19 |
* Lesser General Public License for more details.
|
20 |
*
|
21 |
* This package contains formulas from the PROJ package of USGS.
|
22 |
* USGS's work is fully acknowledged here. This derived work has
|
23 |
* been relicensed under LGPL with Frank Warmerdam's permission.
|
24 |
*/
|
25 |
package org.geotools.referencing.operation.projection; |
26 |
|
27 |
// J2SE dependencies and extensions
|
28 |
import java.awt.geom.AffineTransform; |
29 |
import java.awt.geom.Point2D; |
30 |
import java.util.Collection; |
31 |
|
32 |
// OpenGIS dependencies
|
33 |
import org.opengis.metadata.Identifier; |
34 |
import org.opengis.parameter.ParameterDescriptor; |
35 |
import org.opengis.parameter.ParameterDescriptorGroup; |
36 |
import org.opengis.parameter.ParameterNotFoundException; |
37 |
import org.opengis.parameter.ParameterValueGroup; |
38 |
import org.opengis.referencing.operation.CylindricalProjection; |
39 |
import org.opengis.referencing.operation.MathTransform; |
40 |
|
41 |
// Geotools dependencies
|
42 |
import org.geotools.metadata.iso.citation.CitationImpl; |
43 |
import org.geotools.referencing.NamedIdentifier; |
44 |
import org.geotools.referencing.operation.transform.ConcatenatedTransform; |
45 |
import org.geotools.referencing.operation.transform.ProjectiveTransform; |
46 |
//import org.geotools.resources.i18n.VocabularyKeys;
|
47 |
//import org.geotools.resources.i18n.Vocabulary;
|
48 |
import org.geotools.resources.cts.ResourceKeys; |
49 |
import org.geotools.resources.cts.Resources; |
50 |
|
51 |
|
52 |
/**
|
53 |
* Transverse Mercator Projection (EPSG code 9807). This
|
54 |
* is a cylindrical projection, in which the cylinder has been rotated 90.
|
55 |
* Instead of being tangent to the equator (or to an other standard latitude),
|
56 |
* it is tangent to a central meridian. Deformation are more important as we
|
57 |
* are going futher from the central meridian. The Transverse Mercator
|
58 |
* projection is appropriate for region wich have a greater extent north-south
|
59 |
* than east-west.
|
60 |
* <p>
|
61 |
*
|
62 |
* The elliptical equations used here are series approximations, and their accuracy
|
63 |
* decreases as points move farther from the central meridian of the projection.
|
64 |
* The forward equations here are accurate to a less than a mm ±10 degrees from
|
65 |
* the central meridian, a few mm ±15 degrees from the
|
66 |
* central meridian and a few cm ±20 degrees from the central meridian.
|
67 |
* The spherical equations are not approximations and should always give the
|
68 |
* correct values.
|
69 |
* <p>
|
70 |
*
|
71 |
* There are a number of versions of the transverse mercator projection
|
72 |
* including the Universal (UTM) and Modified (MTM) Transverses Mercator
|
73 |
* projections. In these cases the earth is divided into zones. For the UTM
|
74 |
* the zones are 6 degrees wide, numbered from 1 to 60 proceeding east from
|
75 |
* 180 degrees longitude, and between lats 84 degrees North and 80
|
76 |
* degrees South. The central meridian is taken as the center of the zone
|
77 |
* and the latitude of origin is the equator. A scale factor of 0.9996 and
|
78 |
* false easting of 500000m is used for all zones and a false northing of 10000000m
|
79 |
* is used for zones in the southern hemisphere.
|
80 |
* <p>
|
81 |
*
|
82 |
* NOTE: formulas used below are not those of Snyder, but rather those
|
83 |
* from the {@code proj4} package of the USGS survey, which
|
84 |
* have been reproduced verbatim. USGS work is acknowledged here.
|
85 |
* <p>
|
86 |
*
|
87 |
* <strong>References:</strong><ul>
|
88 |
* <li> Proj-4.4.6 available at <A HREF="http://www.remotesensing.org/proj">www.remotesensing.org/proj</A><br>
|
89 |
* Relevent files are: PJ_tmerc.c, pj_mlfn.c, pj_fwd.c and pj_inv.c </li>
|
90 |
* <li> John P. Snyder (Map Projections - A Working Manual,
|
91 |
* U.S. Geological Survey Professional Paper 1395, 1987)</li>
|
92 |
* <li> "Coordinate Conversions and Transformations including Formulas",
|
93 |
* EPSG Guidence Note Number 7, Version 19.</li>
|
94 |
* </ul>
|
95 |
*
|
96 |
* @see <A HREF="http://mathworld.wolfram.com/MercatorProjection.html">Transverse Mercator projection on MathWorld</A>
|
97 |
* @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/transverse_mercator.html">"Transverse_Mercator" on Remote Sensing</A>
|
98 |
*
|
99 |
* @since 2.1
|
100 |
* @source $URL: http://svn.geotools.org/geotools/tags/2.3.0/module/referencing/src/org/geotools/referencing/operation/projection/TransverseMercator.java $
|
101 |
* @version $Id: TransverseMercator.java 12186 2007-06-18 08:46:09Z jlgomez $
|
102 |
* @author Andr Gosselin
|
103 |
* @author Martin Desruisseaux
|
104 |
* @author Rueben Schulz
|
105 |
*/
|
106 |
public class TransverseMercator extends MapProjection { |
107 |
/**
|
108 |
* A derived quantity of excentricity, computed by <code>e' = (a鲲-b)/b = es/(1-es)</code>
|
109 |
* where <var>a</var> is the semi-major axis length and <var>b</bar> is the semi-minor axis
|
110 |
* length.
|
111 |
*/
|
112 |
private final double esp; |
113 |
|
114 |
/**
|
115 |
* Meridian distance at the {@code latitudeOfOrigin}.
|
116 |
* Used for calculations for the ellipsoid.
|
117 |
*/
|
118 |
private final double ml0; |
119 |
|
120 |
/**
|
121 |
* Constant needed for the <code>mlfn<code> method.
|
122 |
* Setup at construction time.
|
123 |
*/
|
124 |
private final double en0,en1,en2,en3,en4; |
125 |
|
126 |
/**
|
127 |
* Constants used to calculate {@link #en0}, {@link #en1},
|
128 |
* {@link #en2}, {@link #en3}, {@link #en4}.
|
129 |
*/
|
130 |
private static final double C00= 1.0, |
131 |
C02= 0.25,
|
132 |
C04= 0.046875,
|
133 |
C06= 0.01953125,
|
134 |
C08= 0.01068115234375,
|
135 |
C22= 0.75,
|
136 |
C44= 0.46875,
|
137 |
C46= 0.01302083333333333333,
|
138 |
C48= 0.00712076822916666666,
|
139 |
C66= 0.36458333333333333333,
|
140 |
C68= 0.00569661458333333333,
|
141 |
C88= 0.3076171875;
|
142 |
|
143 |
/**
|
144 |
* Contants used for the forward and inverse transform for the eliptical
|
145 |
* case of the Transverse Mercator.
|
146 |
*/
|
147 |
private static final double FC1= 1.00000000000000000000000, // 1/1 |
148 |
FC2= 0.50000000000000000000000, // 1/2 |
149 |
FC3= 0.16666666666666666666666, // 1/6 |
150 |
FC4= 0.08333333333333333333333, // 1/12 |
151 |
FC5= 0.05000000000000000000000, // 1/20 |
152 |
FC6= 0.03333333333333333333333, // 1/30 |
153 |
FC7= 0.02380952380952380952380, // 1/42 |
154 |
FC8= 0.01785714285714285714285; // 1/56 |
155 |
|
156 |
/**
|
157 |
* Relative iteration precision used in the <code>mlfn<code> method. This
|
158 |
* overrides the value in the MapProjection class.
|
159 |
*/
|
160 |
private static final double TOL = 1E-11; |
161 |
|
162 |
/**
|
163 |
* Constructs a new map projection from the supplied parameters.
|
164 |
*
|
165 |
* @param parameters The parameter values in standard units.
|
166 |
* @throws ParameterNotFoundException if a mandatory parameter is missing.
|
167 |
*/
|
168 |
protected TransverseMercator(final ParameterValueGroup parameters) |
169 |
throws ParameterNotFoundException
|
170 |
{ |
171 |
//Fetch parameters
|
172 |
super(parameters);
|
173 |
|
174 |
// Compute constants
|
175 |
esp = excentricitySquared / (1.0 - excentricitySquared);
|
176 |
|
177 |
double t;
|
178 |
en0 = C00 - excentricitySquared * (C02 + excentricitySquared * |
179 |
(C04 + excentricitySquared * (C06 + excentricitySquared * C08))); |
180 |
en1 = excentricitySquared * (C22 - excentricitySquared * |
181 |
(C04 + excentricitySquared * (C06 + excentricitySquared * C08))); |
182 |
en2 = (t = excentricitySquared * excentricitySquared) * |
183 |
(C44 - excentricitySquared * (C46 + excentricitySquared * C48)); |
184 |
en3 = (t *= excentricitySquared) * (C66 - excentricitySquared * C68); |
185 |
en4 = t * excentricitySquared * C88; |
186 |
ml0 = mlfn(latitudeOfOrigin, Math.sin(latitudeOfOrigin), Math.cos(latitudeOfOrigin)); |
187 |
} |
188 |
|
189 |
/**
|
190 |
* {@inheritDoc}
|
191 |
*/
|
192 |
public ParameterDescriptorGroup getParameterDescriptors() {
|
193 |
return Provider.PARAMETERS; |
194 |
} |
195 |
|
196 |
/**
|
197 |
* Transforms the specified (<var>x</var>,<var>y</var>) coordinate (units in radians)
|
198 |
* and stores the result in {@code ptDst} (linear distance on a unit sphere).
|
199 |
*/
|
200 |
protected Point2D transformNormalized(double x, double y, Point2D ptDst) |
201 |
throws ProjectionException
|
202 |
{ |
203 |
double sinphi = Math.sin(y); |
204 |
double cosphi = Math.cos(y); |
205 |
|
206 |
double t = (Math.abs(cosphi)>EPS) ? sinphi/cosphi : 0; |
207 |
t *= t; |
208 |
double al = cosphi*x;
|
209 |
double als = al*al;
|
210 |
al /= Math.sqrt(1.0 - excentricitySquared * sinphi*sinphi); |
211 |
double n = esp * cosphi*cosphi;
|
212 |
|
213 |
/* NOTE: meridinal distance at latitudeOfOrigin is always 0 */
|
214 |
y = (mlfn(y, sinphi, cosphi) - ml0 + |
215 |
sinphi*al*x* |
216 |
FC2 * ( 1.0 +
|
217 |
FC4 * als * (5.0 - t + n*(9.0 + 4.0*n) + |
218 |
FC6 * als * (61.0 + t * (t - 58.0) + n*(270.0 - 330.0*t) + |
219 |
FC8 * als * (1385.0 + t * ( t*(543.0 - t) - 3111.0)))))); |
220 |
|
221 |
x = al*(FC1 + FC3 * als*(1.0 - t + n +
|
222 |
FC5 * als * (5.0 + t*(t - 18.0) + n*(14.0 - 58.0*t) + |
223 |
FC7 * als * (61.0+ t*(t*(179.0 - t) - 479.0 ))))); |
224 |
|
225 |
if (ptDst != null) { |
226 |
ptDst.setLocation(x,y); |
227 |
return ptDst;
|
228 |
} |
229 |
return new Point2D.Double(x,y); |
230 |
} |
231 |
|
232 |
/**
|
233 |
* Transforms the specified (<var>x</var>,<var>y</var>) coordinate
|
234 |
* and stores the result in {@code ptDst}.
|
235 |
*/
|
236 |
protected Point2D inverseTransformNormalized(double x, double y, Point2D ptDst) |
237 |
throws ProjectionException
|
238 |
{ |
239 |
double phi = inv_mlfn(ml0 + y);
|
240 |
|
241 |
if (Math.abs(phi) >= (Math.PI/2)) { |
242 |
y = y<0.0 ? -(Math.PI/2) : (Math.PI/2); |
243 |
x = 0.0;
|
244 |
} else {
|
245 |
double sinphi = Math.sin(phi); |
246 |
double cosphi = Math.cos(phi); |
247 |
double t = (Math.abs(cosphi) > EPS) ? sinphi/cosphi : 0.0; |
248 |
double n = esp * cosphi*cosphi;
|
249 |
double con = 1.0 - excentricitySquared * sinphi*sinphi; |
250 |
double d = x*Math.sqrt(con); |
251 |
con *= t; |
252 |
t *= t; |
253 |
double ds = d*d;
|
254 |
|
255 |
y = phi - (con*ds / (1.0 - excentricitySquared)) *
|
256 |
FC2 * (1.0 - ds *
|
257 |
FC4 * (5.0 + t*(3.0 - 9.0*n) + n*(1.0 - 4*n) - ds * |
258 |
FC6 * (61.0 + t*(90.0 - 252.0*n + 45.0*t) + 46.0*n - ds * |
259 |
FC8 * (1385.0 + t*(3633.0 + t*(4095.0 + 1574.0*t)))))); |
260 |
|
261 |
x = d*(FC1 - ds * FC3 * (1.0 + 2.0*t + n - |
262 |
ds*FC5*(5.0 + t*(28.0 + 24* t + 8.0*n) + 6.0*n - |
263 |
ds*FC7*(61.0 + t*(662.0 + t*(1320.0 + 720.0*t))))))/cosphi; |
264 |
} |
265 |
|
266 |
if (ptDst != null) { |
267 |
ptDst.setLocation(x,y); |
268 |
return ptDst;
|
269 |
} |
270 |
return new Point2D.Double(x,y); |
271 |
} |
272 |
|
273 |
/**
|
274 |
* {@inheritDoc}
|
275 |
*/
|
276 |
protected double getToleranceForAssertions(final double longitude, final double latitude) { |
277 |
if (Math.abs(longitude - centralMeridian) > 0.26) { //15 degrees |
278 |
// When far from the valid area, use a larger tolerance.
|
279 |
return 2.5; |
280 |
} else if (Math.abs(longitude - centralMeridian) > 0.22) { //12.5 degrees |
281 |
return 1.0; |
282 |
} else if (Math.abs(longitude - centralMeridian) > 0.17) { //10 degrees |
283 |
return 0.5; |
284 |
} |
285 |
// a normal tolerance
|
286 |
return 1E-6; |
287 |
} |
288 |
|
289 |
|
290 |
/**
|
291 |
* Provides the transform equations for the spherical case of the
|
292 |
* TransverseMercator projection.
|
293 |
*
|
294 |
* @version $Id: TransverseMercator.java 12186 2007-06-18 08:46:09Z jlgomez $
|
295 |
* @author Andr Gosselin
|
296 |
* @author Martin Desruisseaux
|
297 |
* @author Rueben Schulz
|
298 |
*/
|
299 |
private static final class Spherical extends TransverseMercator { |
300 |
/**
|
301 |
* Constructs a new map projection from the suplied parameters.
|
302 |
*
|
303 |
* @param parameters The parameter values in standard units.
|
304 |
* @throws ParameterNotFoundException if a mandatory parameter is missing.
|
305 |
*/
|
306 |
protected Spherical(final ParameterValueGroup parameters) |
307 |
throws ParameterNotFoundException
|
308 |
{ |
309 |
super(parameters);
|
310 |
//assert isSpherical;
|
311 |
} |
312 |
|
313 |
/**
|
314 |
* {@inheritDoc}
|
315 |
*/
|
316 |
protected Point2D transformNormalized(double x, double y, Point2D ptDst) |
317 |
throws ProjectionException
|
318 |
{ |
319 |
// Compute using ellipsoidal formulas, for comparaison later.
|
320 |
double normalizedLongitude = x;
|
321 |
//assert (ptDst = super.transformNormalized(x, y, ptDst)) != null;
|
322 |
|
323 |
double cosphi = Math.cos(y); |
324 |
double b = cosphi * Math.sin(x); |
325 |
if (Math.abs(Math.abs(b) - 1.0) <= EPS) { |
326 |
throw new ProjectionException(Resources.format(ResourceKeys.ERROR_VALUE_TEND_TOWARD_INFINITY)); |
327 |
} |
328 |
|
329 |
//Using Snyder's equation for calculating y, instead of the one used in Proj4
|
330 |
//poential problems when y and x = 90 degrees, but behaves ok in tests
|
331 |
y = Math.atan2(Math.tan(y),Math.cos(x)) - latitudeOfOrigin; /* Snyder 8-3 */ |
332 |
x = 0.5 * Math.log((1.0+b)/(1.0-b)); /* Snyder 8-1 */ |
333 |
|
334 |
//assert Math.abs(ptDst.getX()-x) <= getToleranceForSphereAssertions(normalizedLongitude,0) : ptDst.getX()-x;
|
335 |
//assert Math.abs(ptDst.getY()-y) <= getToleranceForSphereAssertions(normalizedLongitude,0) : ptDst.getY()-y;
|
336 |
if (ptDst != null) { |
337 |
ptDst.setLocation(x,y); |
338 |
return ptDst;
|
339 |
} |
340 |
return new Point2D.Double(x,y); |
341 |
} |
342 |
|
343 |
/**
|
344 |
* {@inheritDoc}
|
345 |
*/
|
346 |
protected Point2D inverseTransformNormalized(double x, double y, Point2D ptDst) |
347 |
throws ProjectionException
|
348 |
{ |
349 |
// Compute using ellipsoidal formulas, for comparaison later.
|
350 |
//assert (ptDst = super.inverseTransformNormalized(x, y, ptDst)) != null;
|
351 |
|
352 |
double t = Math.exp(x); |
353 |
double sinhX = 0.5 * (t-1.0/t); //sinh(x) |
354 |
double cosD = Math.cos(latitudeOfOrigin + y); |
355 |
double phi = Math.asin(Math.sqrt((1.0-cosD*cosD)/(1.0+sinhX*sinhX))); |
356 |
// correct for the fact that we made everything positive using sqrt(x*x)
|
357 |
y = ((y + latitudeOfOrigin)<0.0) ? -phi : phi;
|
358 |
x = (Math.abs(sinhX)<=EPS && Math.abs(cosD)<=EPS) ? |
359 |
0.0 : Math.atan2(sinhX,cosD); |
360 |
|
361 |
//assert Math.abs(ptDst.getX()-x) <= getToleranceForSphereAssertions(x,0) : ptDst.getX()-x;
|
362 |
//assert Math.abs(ptDst.getY()-y) <= getToleranceForSphereAssertions(x,0) : ptDst.getY()-y;
|
363 |
if (ptDst != null) { |
364 |
ptDst.setLocation(x,y); |
365 |
return ptDst;
|
366 |
} |
367 |
return new Point2D.Double(x,y); |
368 |
} |
369 |
|
370 |
/*
|
371 |
* Allow a larger tolerance when comparing spherical to elliptical calculations
|
372 |
* when we are far from the central meridian (elliptical calculations are
|
373 |
* not valid here).
|
374 |
*
|
375 |
* longitude here is standardized (in radians) and centralMeridian has
|
376 |
* already been removed from it.
|
377 |
*/
|
378 |
protected double getToleranceForSphereAssertions(final double longitude, final double latitude) { |
379 |
if (Math.abs(Math.abs(longitude)- Math.PI/2.0) < TOL) { //90 degrees |
380 |
// elliptical equations are at their worst here
|
381 |
return 1E18; |
382 |
} |
383 |
if (Math.abs(longitude) > 0.26) { //15 degrees |
384 |
// When far from the valid area, use a very larger tolerance.
|
385 |
return 1000000; |
386 |
} |
387 |
// a normal tolerance
|
388 |
return 1E-6; |
389 |
} |
390 |
} |
391 |
|
392 |
|
393 |
/**
|
394 |
* Calculates the meridian distance. This is the distance along the central
|
395 |
* meridian from the equator to {@code phi}. Accurate to < 1e-5 meters
|
396 |
* when used in conjuction with typical major axis values.
|
397 |
*
|
398 |
* @param phi latitude to calculate meridian distance for.
|
399 |
* @param sphi sin(phi).
|
400 |
* @param cphi cos(phi).
|
401 |
* @return meridian distance for the given latitude.
|
402 |
*/
|
403 |
private final double mlfn(final double phi, double sphi, double cphi) { |
404 |
cphi *= sphi; |
405 |
sphi *= sphi; |
406 |
return en0 * phi - cphi *
|
407 |
(en1 + sphi * |
408 |
(en2 + sphi * |
409 |
(en3 + sphi * |
410 |
(en4)))); |
411 |
} |
412 |
|
413 |
/**
|
414 |
* Calculates the latitude ({@code phi}) from a meridian distance.
|
415 |
* Determines phi to TOL (1e-11) radians, about 1e-6 seconds.
|
416 |
*
|
417 |
* @param arg meridian distance to calulate latitude for.
|
418 |
* @return the latitude of the meridian distance.
|
419 |
* @throws ProjectionException if the itteration does not converge.
|
420 |
*/
|
421 |
private final double inv_mlfn(double arg) throws ProjectionException { |
422 |
double s, t, phi, k = 1.0/(1.0 - excentricitySquared); |
423 |
int i;
|
424 |
phi = arg; |
425 |
for (i=MAX_ITER; true;) { // rarely goes over 5 iterations |
426 |
if (--i < 0) { |
427 |
throw new ProjectionException(Resources.format(ResourceKeys.ERROR_NO_CONVERGENCE)); |
428 |
} |
429 |
s = Math.sin(phi);
|
430 |
t = 1.0 - excentricitySquared * s * s;
|
431 |
t = (mlfn(phi, s, Math.cos(phi)) - arg) * (t * Math.sqrt(t)) * k; |
432 |
phi -= t; |
433 |
if (Math.abs(t) < TOL) { |
434 |
return phi;
|
435 |
} |
436 |
} |
437 |
} |
438 |
|
439 |
/**
|
440 |
* Convenience method computing the zone code from the central meridian.
|
441 |
* Information about zones convention must be specified in argument. Two
|
442 |
* widely set of arguments are of Universal Transverse Mercator (UTM) and
|
443 |
* Modified Transverse Mercator (MTM) projections:<br>
|
444 |
* <p>
|
445 |
*
|
446 |
* UTM projection (zones numbered from 1 to 60):<br>
|
447 |
* <p>
|
448 |
* {@code getZone(-177, 6);}<br>
|
449 |
* <p>
|
450 |
* MTM projection (zones numbered from 1 to 120):<br>
|
451 |
* <p>
|
452 |
* {@code getZone(-52.5, -3);}<br>
|
453 |
*
|
454 |
* @param centralLongitudeZone1 Longitude in the middle of zone 1, in decimal degrees
|
455 |
* relative to Greenwich. Positive longitudes are toward east, and negative
|
456 |
* longitudes toward west.
|
457 |
* @param zoneWidth Number of degrees of longitudes in one zone. A positive value
|
458 |
* means that zones are numbered from west to east (i.e. in the direction of
|
459 |
* positive longitudes). A negative value means that zones are numbered from
|
460 |
* east to west.
|
461 |
* @return The zone number. First zone is numbered 1.
|
462 |
*/
|
463 |
private int getZone(final double centralLongitudeZone1, final double zoneWidth) { |
464 |
final double zoneCount = Math.abs(360/zoneWidth); |
465 |
double t;
|
466 |
t = centralLongitudeZone1 - 0.5*zoneWidth; // Longitude at the beginning of the first zone. |
467 |
t = Math.toDegrees(centralMeridian) - t; // Degrees of longitude between the central longitude and longitude 1. |
468 |
t = Math.floor(t/zoneWidth + EPS); // Number of zones between the central longitude and longitude 1. |
469 |
t -= zoneCount*Math.floor(t/zoneCount); // If negative, bring back to the interval 0 to (zoneCount-1). |
470 |
return ((int) t)+1; |
471 |
} |
472 |
|
473 |
/**
|
474 |
* Convenience method returning the meridian in the middle of current zone. This meridian is
|
475 |
* typically the central meridian. This method may be invoked to make sure that the central
|
476 |
* meridian is correctly set.
|
477 |
*
|
478 |
* @param centralLongitudeZone1 Longitude in the middle of zone 1, in decimal degrees
|
479 |
* relative to Greenwich. Positive longitudes are toward east, and negative
|
480 |
* longitudes toward west.
|
481 |
* @param zoneWidth Number of degrees of longitudes in one zone. A positive value
|
482 |
* means that zones are numbered from west to east (i.e. in the direction of
|
483 |
* positive longitudes). A negative value means that zones are numbered from
|
484 |
* east to west.
|
485 |
* @return The central meridian.
|
486 |
*/
|
487 |
private double getCentralMedirian(final double centralLongitudeZone1, final double zoneWidth) { |
488 |
double t;
|
489 |
t = centralLongitudeZone1 + (getZone(centralLongitudeZone1, zoneWidth)-1)*zoneWidth;
|
490 |
t -= 360*Math.floor((t+180)/360); // Bring back into [-180..+180] range. |
491 |
return t;
|
492 |
} |
493 |
|
494 |
/**
|
495 |
* Convenience method computing the zone code from the central meridian.
|
496 |
*
|
497 |
* @return The zone number, using the scalefactor and false easting
|
498 |
* to decide if this is a UTM or MTM case. Returns 0 if the
|
499 |
* case of the projection cannot be determined.
|
500 |
*/
|
501 |
public int getZone() { |
502 |
// UTM
|
503 |
if (scaleFactor == 0.9996 && falseEasting == 500000.0) { |
504 |
return getZone(-177, 6); |
505 |
} |
506 |
// MTM
|
507 |
if (scaleFactor == 0.9999 && falseEasting == 304800.0){ |
508 |
return getZone(-52.5, -3); |
509 |
} |
510 |
// unknown
|
511 |
throw new IllegalStateException(Resources.format(ResourceKeys.ERROR_UNKNOW_TYPE_$1));//.UNKNOW_PROJECTION_TYPE)); |
512 |
} |
513 |
|
514 |
/**
|
515 |
* Convenience method returning the meridian in the middle of current zone. This meridian is
|
516 |
* typically the central meridian. This method may be invoked to make sure that the central
|
517 |
* meridian is correctly set.
|
518 |
*
|
519 |
* @return The central meridian, using the scalefactor and false easting
|
520 |
* to decide if this is a UTM or MTM case. Returns {@link Double#NaN}
|
521 |
* if the case of the projection cannot be determined.
|
522 |
*/
|
523 |
public double getCentralMeridian() { |
524 |
// UTM
|
525 |
if (scaleFactor == 0.9996 && falseEasting == 500000.0) { |
526 |
return getCentralMedirian(-177, 6); |
527 |
} |
528 |
// MTM
|
529 |
if (scaleFactor == 0.9999 && falseEasting == 304800.0){ |
530 |
return getCentralMedirian(-52.5, -3); |
531 |
} |
532 |
// unknown
|
533 |
throw new IllegalStateException(Resources.format(ResourceKeys.ERROR_UNKNOW_TYPE_$1));//.UNKNOW_PROJECTION_TYPE)); |
534 |
} |
535 |
|
536 |
/**
|
537 |
* Returns a hash value for this projection.
|
538 |
*/
|
539 |
public int hashCode() { |
540 |
final long code = Double.doubleToLongBits(ml0); |
541 |
return ((int)code ^ (int)(code >>> 32)) + 37*super.hashCode(); |
542 |
} |
543 |
|
544 |
/**
|
545 |
* Compares the specified object with this map projection for equality.
|
546 |
*/
|
547 |
public boolean equals(final Object object) { |
548 |
if (object == this) { |
549 |
// Slight optimization
|
550 |
return true; |
551 |
} |
552 |
// Relevant parameters are already compared in MapProjection
|
553 |
return super.equals(object); |
554 |
} |
555 |
|
556 |
|
557 |
|
558 |
|
559 |
//////////////////////////////////////////////////////////////////////////////////////////
|
560 |
//////////////////////////////////////////////////////////////////////////////////////////
|
561 |
//////// ////////
|
562 |
//////// PROVIDER ////////
|
563 |
//////// ////////
|
564 |
//////////////////////////////////////////////////////////////////////////////////////////
|
565 |
//////////////////////////////////////////////////////////////////////////////////////////
|
566 |
|
567 |
/**
|
568 |
* The {@link org.geotools.referencing.operation.MathTransformProvider}
|
569 |
* for a {@link TransverseMercator} projection.
|
570 |
*
|
571 |
* @see org.geotools.referencing.operation.DefaultMathTransformFactory
|
572 |
*
|
573 |
* @since 2.1
|
574 |
* @version $Id: TransverseMercator.java 12186 2007-06-18 08:46:09Z jlgomez $
|
575 |
* @author Martin Desruisseaux
|
576 |
* @author Rueben Schulz
|
577 |
*/
|
578 |
public static class Provider extends AbstractProvider { |
579 |
/**
|
580 |
* Returns a descriptor group for the specified parameters.
|
581 |
*/
|
582 |
static ParameterDescriptorGroup createDescriptorGroup(final Identifier[] identifiers) { |
583 |
return createDescriptorGroup(identifiers, new ParameterDescriptor[] { |
584 |
SEMI_MAJOR, SEMI_MINOR, |
585 |
CENTRAL_MERIDIAN, LATITUDE_OF_ORIGIN, |
586 |
SCALE_FACTOR, FALSE_EASTING, |
587 |
FALSE_NORTHING |
588 |
}); |
589 |
} |
590 |
|
591 |
/**
|
592 |
* The parameters group.
|
593 |
*/
|
594 |
static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(new NamedIdentifier[] { |
595 |
new NamedIdentifier(CitationImpl.OGC, "Transverse_Mercator"), |
596 |
new NamedIdentifier(CitationImpl.EPSG, "Transverse Mercator"), |
597 |
new NamedIdentifier(CitationImpl.EPSG, "Gauss-Kruger"), |
598 |
new NamedIdentifier(CitationImpl.EPSG, "9807"), |
599 |
new NamedIdentifier(CitationImpl.GEOTIFF, "CT_TransverseMercator"), |
600 |
new NamedIdentifier(CitationImpl.ESRI, "Transverse_Mercator"), |
601 |
new NamedIdentifier(CitationImpl.ESRI, "Gauss_Kruger")//, |
602 |
//new NamedIdentifier(CitationImpl.GEOTOOLS, Vocabulary.formatInternational(
|
603 |
// VocabularyKeys.TRANSVERSE_MERCATOR_PROJECTION))
|
604 |
}); |
605 |
|
606 |
/**
|
607 |
* Constructs a new provider.
|
608 |
*/
|
609 |
public Provider() { |
610 |
super(PARAMETERS);
|
611 |
} |
612 |
|
613 |
/**
|
614 |
* Constructs a new provider with the specified parameters.
|
615 |
*/
|
616 |
Provider(final ParameterDescriptorGroup descriptor) { |
617 |
super(descriptor);
|
618 |
} |
619 |
|
620 |
/**
|
621 |
* Returns the operation type for this map projection.
|
622 |
*/
|
623 |
protected Class getOperationType() { |
624 |
return CylindricalProjection.class;
|
625 |
} |
626 |
|
627 |
/**
|
628 |
* Creates a transform from the specified group of parameter values.
|
629 |
*
|
630 |
* @param parameters The group of parameter values.
|
631 |
* @return The created math transform.
|
632 |
* @throws ParameterNotFoundException if a required parameter was not found.
|
633 |
*/
|
634 |
public MathTransform createMathTransform(final ParameterValueGroup parameters) |
635 |
throws ParameterNotFoundException
|
636 |
{ |
637 |
if (isSpherical(parameters)) {
|
638 |
return new Spherical(parameters); |
639 |
} else {
|
640 |
return new TransverseMercator(parameters); |
641 |
} |
642 |
} |
643 |
} |
644 |
|
645 |
/**
|
646 |
* The {@link org.geotools.referencing.operation.MathTransformProvider} for a South Orientated
|
647 |
* {@link TransverseMercator} projection (EPSG code 9808). Note that at the contrary of what
|
648 |
* this class name suggest, the projected coordinates are still increasing toward North. This
|
649 |
* is because all {@link MapProjection}s must complies with standard axis orientations. The
|
650 |
* real axis flip is performed by the CRS framework outside this package.
|
651 |
* See "<cite>Axis units and orientation</cite>" in
|
652 |
* {@linkplain org.geotools.referencing.operation.projection package description} for details.
|
653 |
* <p>
|
654 |
* The usual Transverse Mercator formulas are:
|
655 |
* <p>
|
656 |
* <ul>
|
657 |
* <li><var>easting</var> = <var>{@linkplain #falseEasting false easting}</var> + <var>px</var></li>
|
658 |
* <li><var>northing</var> = <var>{@linkplain #falseNorthing false northing}</var> + <var>py</var></li>
|
659 |
* </ul>
|
660 |
* <p>
|
661 |
* The Transverse Mercator South Orientated Projection formulas are:
|
662 |
* <p>
|
663 |
* <ul>
|
664 |
* <li><var>westing</var> = <var>{@linkplain #falseEasting false easting}</var> - <var>px</var></li>
|
665 |
* <li><var>southing</var> = <var>{@linkplain #falseNorthing false northing}</var> - <var>py</var></li>
|
666 |
* </ul>
|
667 |
* <p>
|
668 |
* Where the <var>px</var> and <var>py</var> terms are the same in both cases.
|
669 |
* Transforms created by this provider actually computes
|
670 |
* (<var>easting</var>,<var>northing</var>) = (-<var>westing</var>,-<var>southing</var>).
|
671 |
* This is equivalent to a {@link TransverseMercator} projection with
|
672 |
* {@link #falseEasting falseEasting} and {@link #falseNorthing falseNorthing} sign reverted.
|
673 |
* This operation is implemented as a concatenation of a North-orientated transverse mercator
|
674 |
* projection with an affine transform for (<var>false easting</var>,<var>false northing</var>)
|
675 |
* correction.
|
676 |
*
|
677 |
* @since 2.2
|
678 |
* @version $Id: TransverseMercator.java 12186 2007-06-18 08:46:09Z jlgomez $
|
679 |
* @author Martin Desruisseaux
|
680 |
*/
|
681 |
public static class Provider_SouthOrientated extends Provider { |
682 |
/**
|
683 |
* Constructs a new provider.
|
684 |
*/
|
685 |
public Provider_SouthOrientated() {
|
686 |
super(createDescriptorGroup(new NamedIdentifier[] { |
687 |
new NamedIdentifier(CitationImpl.EPSG, "Transverse Mercator (South Orientated)"), |
688 |
new NamedIdentifier(CitationImpl.EPSG, "9808") |
689 |
})); |
690 |
} |
691 |
|
692 |
/**
|
693 |
* Creates a transform from the specified group of parameter values.
|
694 |
*
|
695 |
* @param parameters The group of parameter values.
|
696 |
* @return The created math transform.
|
697 |
* @throws ParameterNotFoundException if a required parameter was not found.
|
698 |
*/
|
699 |
public MathTransform createMathTransform(final ParameterValueGroup parameters) |
700 |
throws ParameterNotFoundException
|
701 |
{ |
702 |
final MapProjection projection = (MapProjection) super.createMathTransform(parameters); |
703 |
if (projection.falseEasting==0 && projection.falseNorthing==0) { |
704 |
return projection;
|
705 |
} |
706 |
final AffineTransform step = AffineTransform.getTranslateInstance( |
707 |
-2*projection.falseEasting, -2*projection.falseNorthing); |
708 |
return ConcatenatedTransform.create(ProjectiveTransform.create(step), projection);
|
709 |
} |
710 |
} |
711 |
} |