root / branches / F2 / libraries / libJCRS / src / org / geotools / referencing / operation / projection / TransverseMercator.java @ 12186
History | View | Annotate | Download (29.8 KB)
1 | 12186 | jlgomez | /*
|
---|---|---|---|
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$
|
||
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$
|
||
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$
|
||
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$
|
||
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 | } |