1 | /*
|
---|
2 | * GeoTools - The Open Source Java GIS Toolkit
|
---|
3 | * http://geotools.org
|
---|
4 | *
|
---|
5 | * (C) 2006-2008, Open Source Geospatial Foundation (OSGeo)
|
---|
6 | *
|
---|
7 | * This library is free software; you can redistribute it and/or
|
---|
8 | * modify it under the terms of the GNU Lesser General Public
|
---|
9 | * License as published by the Free Software Foundation;
|
---|
10 | * version 2.1 of the License.
|
---|
11 | *
|
---|
12 | * This library is distributed in the hope that it will be useful,
|
---|
13 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
14 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
---|
15 | * Lesser General Public License for more details.
|
---|
16 | *
|
---|
17 | * This package contains formulas from the PROJ package of USGS.
|
---|
18 | * USGS's work is fully acknowledged here. This derived work has
|
---|
19 | * been relicensed under LGPL with Frank Warmerdam's permission.
|
---|
20 | */
|
---|
21 | package org.geotools.referencing.operation.projection;
|
---|
22 |
|
---|
23 | import java.awt.geom.Point2D;
|
---|
24 | import java.util.Collection;
|
---|
25 | import javax.measure.unit.NonSI;
|
---|
26 | import org.opengis.parameter.GeneralParameterDescriptor;
|
---|
27 | import org.opengis.parameter.ParameterDescriptor;
|
---|
28 | import org.opengis.parameter.ParameterDescriptorGroup;
|
---|
29 | import org.opengis.parameter.ParameterNotFoundException;
|
---|
30 | import org.opengis.parameter.ParameterValueGroup;
|
---|
31 | import org.opengis.referencing.operation.MathTransform;
|
---|
32 | import org.geotools.metadata.iso.citation.Citations;
|
---|
33 | import org.geotools.referencing.NamedIdentifier;
|
---|
34 | import org.geotools.resources.i18n.ErrorKeys;
|
---|
35 |
|
---|
36 | import static java.lang.Math.*;
|
---|
37 |
|
---|
38 |
|
---|
39 | /**
|
---|
40 | * Lambert Azimuthal Equal Area (EPSG code 9820).
|
---|
41 | * <p>
|
---|
42 | * <b>References:</b>
|
---|
43 | * <ul>
|
---|
44 | * <li> A. Annoni, C. Luzet, E.Gubler and J. Ihde - Map Projections for Europe</li>
|
---|
45 | * <li> John P. Snyder (Map Projections - A Working Manual,
|
---|
46 | * U.S. Geological Survey Professional Paper 1395)</li>
|
---|
47 | * </ul>
|
---|
48 | *
|
---|
49 | * @see <A HREF="http://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html">Lambert Azimuthal Equal-Area Projection on MathWorld</A>
|
---|
50 | * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_azimuthal_equal_area.html">"Lambert_Azimuthal_Equal_Area" on RemoteSensing.org</A>
|
---|
51 | *
|
---|
52 | * @since 2.4
|
---|
53 | * @version $Id: LambertAzimuthalEqualArea.java 37299 2011-05-25 05:21:24Z mbedward $
|
---|
54 | *
|
---|
55 | * @source $URL: http://svn.osgeo.org/geotools/branches/2.7.x/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/LambertAzimuthalEqualArea.java $
|
---|
56 | * @author Gerald Evenden (for original code in Proj4)
|
---|
57 | * @author Beate Stollberg
|
---|
58 | * @author Martin Desruisseaux
|
---|
59 | */
|
---|
60 | public class LambertAzimuthalEqualArea extends MapProjection {
|
---|
61 | /** For cross-version compatibility. */
|
---|
62 | private static final long serialVersionUID = 1639914708790574760L;
|
---|
63 |
|
---|
64 | /** Maximum difference allowed when comparing real numbers. */
|
---|
65 | private static final double EPSILON = 1E-7;
|
---|
66 |
|
---|
67 | /** Epsilon for the comparaison of small quantities. */
|
---|
68 | private static final double FINE_EPSILON = 1E-10;
|
---|
69 |
|
---|
70 | /** Epsilon for the comparaison of latitudes. */
|
---|
71 | private static final double EPSILON_LATITUDE = 1E-10;
|
---|
72 |
|
---|
73 | /** Constants for authalic latitude. */
|
---|
74 | private static final double P00 = 0.33333333333333333333,
|
---|
75 | P01 = 0.17222222222222222222,
|
---|
76 | P02 = 0.10257936507936507936,
|
---|
77 | P10 = 0.06388888888888888888,
|
---|
78 | P11 = 0.06640211640211640211,
|
---|
79 | P20 = 0.01641501294219154443;
|
---|
80 |
|
---|
81 | /** The projection mode. */
|
---|
82 | static final int OBLIQUE=0, EQUATORIAL=1, NORTH_POLE=2, SOUTH_POLE=3;
|
---|
83 |
|
---|
84 | /** The projection mode for this particular instance. */
|
---|
85 | final int mode;
|
---|
86 |
|
---|
87 | /** Constant parameters. */
|
---|
88 | final double sinb1, cosb1, xmf, ymf, mmf, qp, dd, rq;
|
---|
89 |
|
---|
90 | /** Coefficients for authalic latitude. */
|
---|
91 | private final double APA0, APA1, APA2;
|
---|
92 |
|
---|
93 | /**
|
---|
94 | * Constructs a new map projection from the supplied parameters.
|
---|
95 | *
|
---|
96 | * @param parameters The parameter values in standard units.
|
---|
97 | * @throws ParameterNotFoundException if a mandatory parameter is missing.
|
---|
98 | */
|
---|
99 | protected LambertAzimuthalEqualArea(final ParameterValueGroup parameters)
|
---|
100 | throws ParameterNotFoundException
|
---|
101 | {
|
---|
102 | // Fetch parameters
|
---|
103 | super(parameters);
|
---|
104 | final Collection<GeneralParameterDescriptor> expected = getParameterDescriptors().descriptors();
|
---|
105 | latitudeOfOrigin = doubleValue(expected, Provider.LATITUDE_OF_CENTRE, parameters);
|
---|
106 | centralMeridian = doubleValue(expected, Provider.LONGITUDE_OF_CENTRE, parameters);
|
---|
107 | ensureLatitudeInRange (Provider.LATITUDE_OF_CENTRE, latitudeOfOrigin, true);
|
---|
108 | ensureLongitudeInRange(Provider.LONGITUDE_OF_CENTRE, centralMeridian, true);
|
---|
109 | /*
|
---|
110 | * Detects the mode (oblique, etc.).
|
---|
111 | */
|
---|
112 | final double t = abs(latitudeOfOrigin);
|
---|
113 | if (abs(t - PI/2) < EPSILON_LATITUDE) {
|
---|
114 | mode = latitudeOfOrigin < 0.0 ? SOUTH_POLE : NORTH_POLE;
|
---|
115 | } else if (abs(t) < EPSILON_LATITUDE) {
|
---|
116 | mode = EQUATORIAL;
|
---|
117 | } else {
|
---|
118 | mode = OBLIQUE;
|
---|
119 | }
|
---|
120 | /*
|
---|
121 | * Computes the constants for authalic latitude.
|
---|
122 | */
|
---|
123 | final double es2 = excentricitySquared * excentricitySquared;
|
---|
124 | final double es3 = excentricitySquared * es2;
|
---|
125 | APA0 = P02 * es3 + P01 * es2 + P00 * excentricitySquared;
|
---|
126 | APA1 = P11 * es3 + P10 * es2;
|
---|
127 | APA2 = P20 * es3;
|
---|
128 |
|
---|
129 | final double sinphi;
|
---|
130 | qp = qsfn(1);
|
---|
131 | rq = sqrt(0.5 * qp);
|
---|
132 | mmf = 0.5 / (1 - excentricitySquared);
|
---|
133 | sinphi = sin(latitudeOfOrigin);
|
---|
134 | if (isSpherical) {
|
---|
135 | sinb1 = sin(latitudeOfOrigin);
|
---|
136 | cosb1 = cos(latitudeOfOrigin);
|
---|
137 | } else {
|
---|
138 | sinb1 = qsfn(sinphi) / qp;
|
---|
139 | cosb1 = sqrt(1.0 - sinb1 * sinb1);
|
---|
140 | }
|
---|
141 | switch (mode) {
|
---|
142 | case NORTH_POLE: // Fall through
|
---|
143 | case SOUTH_POLE: {
|
---|
144 | dd = 1.0;
|
---|
145 | xmf = ymf = rq;
|
---|
146 | break;
|
---|
147 | }
|
---|
148 | case EQUATORIAL: {
|
---|
149 | dd = 1.0 / rq;
|
---|
150 | xmf = 1.0;
|
---|
151 | ymf = 0.5 * qp;
|
---|
152 | break;
|
---|
153 | }
|
---|
154 | case OBLIQUE: {
|
---|
155 | dd = cos(latitudeOfOrigin) /
|
---|
156 | (sqrt(1.0 - excentricitySquared * sinphi * sinphi) * rq * cosb1);
|
---|
157 | xmf = rq * dd;
|
---|
158 | ymf = rq / dd;
|
---|
159 | break;
|
---|
160 | }
|
---|
161 | default: {
|
---|
162 | throw new AssertionError(mode);
|
---|
163 | }
|
---|
164 | }
|
---|
165 | }
|
---|
166 |
|
---|
167 | /**
|
---|
168 | * {@inheritDoc}
|
---|
169 | */
|
---|
170 | public ParameterDescriptorGroup getParameterDescriptors() {
|
---|
171 | return Provider.PARAMETERS;
|
---|
172 | }
|
---|
173 |
|
---|
174 | /**
|
---|
175 | * {@inheritDoc}
|
---|
176 | */
|
---|
177 | @Override
|
---|
178 | public ParameterValueGroup getParameterValues() {
|
---|
179 | final ParameterValueGroup values = super.getParameterValues();
|
---|
180 | final Collection<GeneralParameterDescriptor> expected = getParameterDescriptors().descriptors();
|
---|
181 | set(expected, Provider.LATITUDE_OF_CENTRE, values, latitudeOfOrigin);
|
---|
182 | set(expected, Provider.LONGITUDE_OF_CENTRE, values, centralMeridian);
|
---|
183 | return values;
|
---|
184 | }
|
---|
185 |
|
---|
186 | /**
|
---|
187 | * Transforms the specified (<var>λ</var>,<var>φ</var>) coordinates
|
---|
188 | * (units in radians) and stores the result in {@code ptDst} (linear distance
|
---|
189 | * on a unit sphere).
|
---|
190 | */
|
---|
191 | protected Point2D transformNormalized(final double lambda, final double phi, Point2D ptDst)
|
---|
192 | throws ProjectionException
|
---|
193 | {
|
---|
194 | final double coslam = cos(lambda);
|
---|
195 | final double sinlam = sin(lambda);
|
---|
196 | final double sinphi = sin(phi);
|
---|
197 | double q = qsfn(sinphi);
|
---|
198 | final double sinb, cosb, b, c, x, y;
|
---|
199 | switch (mode) {
|
---|
200 | case OBLIQUE: {
|
---|
201 | sinb = q / qp;
|
---|
202 | cosb = sqrt(1.0 - sinb * sinb);
|
---|
203 | c = 1.0 + sinb1 * sinb + cosb1 * cosb * coslam;
|
---|
204 | b = sqrt(2.0 / c);
|
---|
205 | y = ymf * b * (cosb1 * sinb - sinb1 * cosb * coslam);
|
---|
206 | x = xmf * b * cosb * sinlam;
|
---|
207 | break;
|
---|
208 | }
|
---|
209 | case EQUATORIAL: {
|
---|
210 | sinb = q / qp;
|
---|
211 | cosb = sqrt(1.0 - sinb * sinb);
|
---|
212 | c = 1.0 + cosb * coslam;
|
---|
213 | b = sqrt(2.0 / c);
|
---|
214 | y = ymf * b * sinb;
|
---|
215 | x = xmf * b * cosb * sinlam;
|
---|
216 | break;
|
---|
217 | }
|
---|
218 | case NORTH_POLE: {
|
---|
219 | c = (PI / 2) + phi;
|
---|
220 | q = qp - q;
|
---|
221 | if (q >= 0.0) {
|
---|
222 | b = sqrt(q);
|
---|
223 | x = b * sinlam;
|
---|
224 | y = coslam * -b;
|
---|
225 | } else {
|
---|
226 | x = y = 0.;
|
---|
227 | }
|
---|
228 | break;
|
---|
229 | }
|
---|
230 | case SOUTH_POLE: {
|
---|
231 | c = phi - (PI / 2);
|
---|
232 | q = qp + q;
|
---|
233 | if (q >= 0.0) {
|
---|
234 | b = sqrt(q);
|
---|
235 | x = b * sinlam;
|
---|
236 | y = coslam * +b;
|
---|
237 | } else {
|
---|
238 | x = y = 0.;
|
---|
239 | }
|
---|
240 | break;
|
---|
241 | }
|
---|
242 | default: {
|
---|
243 | throw new AssertionError(mode);
|
---|
244 | }
|
---|
245 | }
|
---|
246 | if (abs(c) < EPSILON_LATITUDE) {
|
---|
247 | throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
|
---|
248 | }
|
---|
249 | if (ptDst != null) {
|
---|
250 | ptDst.setLocation(x,y);
|
---|
251 | return ptDst;
|
---|
252 | }
|
---|
253 | return new Point2D.Double(x,y);
|
---|
254 | }
|
---|
255 |
|
---|
256 | /**
|
---|
257 | * Transforms the specified (<var>x</var>,<var>y</var>) coordinate
|
---|
258 | * and stores the result in {@code ptDst}.
|
---|
259 | */
|
---|
260 | @Override
|
---|
261 | @SuppressWarnings("fallthrough")
|
---|
262 | protected Point2D inverseTransformNormalized(double x, double y, Point2D ptDst)
|
---|
263 | throws ProjectionException
|
---|
264 | {
|
---|
265 | final double lambda, phi;
|
---|
266 | switch (mode) {
|
---|
267 | case EQUATORIAL: // Fall through
|
---|
268 | case OBLIQUE: {
|
---|
269 | x /= dd;
|
---|
270 | y *= dd;
|
---|
271 | final double rho = hypot(x, y);
|
---|
272 | if (rho < FINE_EPSILON) {
|
---|
273 | lambda = 0.0;
|
---|
274 | phi = latitudeOfOrigin;
|
---|
275 | } else {
|
---|
276 | double sCe, cCe, /*q,*/ ab;
|
---|
277 | sCe = 2.0 * asin(0.5 * rho / rq);
|
---|
278 | cCe = cos(sCe);
|
---|
279 | sCe = sin(sCe);
|
---|
280 | x *= sCe;
|
---|
281 | if (mode == OBLIQUE) {
|
---|
282 | ab = cCe * sinb1 + y * sCe * cosb1 / rho;
|
---|
283 | //q = qp * ab;
|
---|
284 | y = rho * cosb1 * cCe - y * sinb1 * sCe;
|
---|
285 | } else {
|
---|
286 | ab = y * sCe / rho;
|
---|
287 | //q = qp * ab;
|
---|
288 | y = rho * cCe;
|
---|
289 | }
|
---|
290 | lambda = atan2(x, y);
|
---|
291 | phi = authlat(asin(ab));
|
---|
292 | }
|
---|
293 | break;
|
---|
294 | }
|
---|
295 | case NORTH_POLE: {
|
---|
296 | y = -y;
|
---|
297 | // Fall through
|
---|
298 | }
|
---|
299 | case SOUTH_POLE: {
|
---|
300 | final double q = x*x + y*y;
|
---|
301 | if (q == 0) {
|
---|
302 | lambda = 0.;
|
---|
303 | phi = latitudeOfOrigin;
|
---|
304 | } else {
|
---|
305 | double ab = 1.0 - q / qp;
|
---|
306 | if (mode == SOUTH_POLE) {
|
---|
307 | ab = -ab;
|
---|
308 | }
|
---|
309 | lambda = atan2(x, y);
|
---|
310 | phi = authlat(asin(ab));
|
---|
311 | }
|
---|
312 | break;
|
---|
313 | }
|
---|
314 | default: {
|
---|
315 | throw new AssertionError(mode);
|
---|
316 | }
|
---|
317 | }
|
---|
318 | if (ptDst != null) {
|
---|
319 | ptDst.setLocation(lambda, phi);
|
---|
320 | return ptDst;
|
---|
321 | }
|
---|
322 | return new Point2D.Double(lambda, phi);
|
---|
323 | }
|
---|
324 |
|
---|
325 |
|
---|
326 | /**
|
---|
327 | * Provides the transform equations for the spherical case.
|
---|
328 | *
|
---|
329 | * @version $Id: LambertAzimuthalEqualArea.java 37299 2011-05-25 05:21:24Z mbedward $
|
---|
330 | * @author Martin Desruisseaux
|
---|
331 | */
|
---|
332 | private static final class Spherical extends LambertAzimuthalEqualArea {
|
---|
333 | /**
|
---|
334 | * For cross-version compatibility.
|
---|
335 | */
|
---|
336 | private static final long serialVersionUID = 2091431369806844342L;
|
---|
337 |
|
---|
338 | /**
|
---|
339 | * Constructs a new map projection from the suplied parameters.
|
---|
340 | *
|
---|
341 | * @param parameters The parameter values in standard units.
|
---|
342 | * @throws ParameterNotFoundException if a mandatory parameter is missing.
|
---|
343 | */
|
---|
344 | protected Spherical(final ParameterValueGroup parameters)
|
---|
345 | throws ParameterNotFoundException
|
---|
346 | {
|
---|
347 | super(parameters);
|
---|
348 | ensureSpherical();
|
---|
349 | }
|
---|
350 |
|
---|
351 | /**
|
---|
352 | * Transforms the specified (<var>λ</var>,<var>φ</var>) coordinates
|
---|
353 | * (units in radians) and stores the result in {@code ptDst} (linear distance
|
---|
354 | * on a unit sphere).
|
---|
355 | */
|
---|
356 | @Override
|
---|
357 | protected Point2D transformNormalized(final double lambda, final double phi, Point2D ptDst)
|
---|
358 | throws ProjectionException
|
---|
359 | {
|
---|
360 | // Compute using ellipsoidal formulas, for comparaison later.
|
---|
361 | assert (ptDst = super.transformNormalized(lambda, phi, ptDst)) != null;
|
---|
362 |
|
---|
363 | final double sinphi = sin(phi);
|
---|
364 | final double cosphi = cos(phi);
|
---|
365 | final double coslam = cos(lambda);
|
---|
366 | double x,y;
|
---|
367 | switch (mode) {
|
---|
368 | case EQUATORIAL: {
|
---|
369 | y = 1.0 + cosphi * coslam;
|
---|
370 | if (y <= FINE_EPSILON) {
|
---|
371 | throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
|
---|
372 | }
|
---|
373 | y = sqrt(2.0 / y);
|
---|
374 | x = y * cosphi * sin(lambda);
|
---|
375 | y *= sinphi;
|
---|
376 | break;
|
---|
377 | }
|
---|
378 | case OBLIQUE: {
|
---|
379 | y = 1.0 + sinb1 * sinphi + cosb1 * cosphi * coslam;
|
---|
380 | if (y <= FINE_EPSILON) {
|
---|
381 | throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
|
---|
382 | }
|
---|
383 | y = sqrt(2.0 / y);
|
---|
384 | x = y * cosphi * sin(lambda);
|
---|
385 | y *= cosb1 * sinphi - sinb1 * cosphi * coslam;
|
---|
386 | break;
|
---|
387 | }
|
---|
388 | case NORTH_POLE: {
|
---|
389 | if (abs(phi + latitudeOfOrigin) < EPSILON_LATITUDE) {
|
---|
390 | throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
|
---|
391 | }
|
---|
392 | y = (PI/4) - phi * 0.5;
|
---|
393 | y = 2.0 * sin(y);
|
---|
394 | x = y * sin(lambda);
|
---|
395 | y *= -coslam;
|
---|
396 | break;
|
---|
397 | }
|
---|
398 | case SOUTH_POLE: {
|
---|
399 | if (abs(phi + latitudeOfOrigin) < EPSILON_LATITUDE) {
|
---|
400 | throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
|
---|
401 | }
|
---|
402 | y = (PI/4) - phi * 0.5;
|
---|
403 | y = 2.0 * cos(y);
|
---|
404 | x = y * sin(lambda);
|
---|
405 | y *= +coslam;
|
---|
406 | break;
|
---|
407 | }
|
---|
408 | default: {
|
---|
409 | throw new AssertionError(mode);
|
---|
410 | }
|
---|
411 | }
|
---|
412 | assert checkTransform(x, y, ptDst);
|
---|
413 | if (ptDst != null) {
|
---|
414 | ptDst.setLocation(x,y);
|
---|
415 | return ptDst;
|
---|
416 | }
|
---|
417 | return new Point2D.Double(x,y);
|
---|
418 | }
|
---|
419 |
|
---|
420 | /**
|
---|
421 | * Transforms the specified (<var>x</var>,<var>y</var>) coordinate
|
---|
422 | * and stores the result in {@code ptDst} using equations for a sphere.
|
---|
423 | */
|
---|
424 | @Override
|
---|
425 | protected Point2D inverseTransformNormalized(double x, double y, Point2D ptDst)
|
---|
426 | throws ProjectionException
|
---|
427 | {
|
---|
428 | // Compute using ellipsoidal formulas, for comparaison later.
|
---|
429 | assert (ptDst = super.inverseTransformNormalized(x, y, ptDst)) != null;
|
---|
430 |
|
---|
431 | double lambda, phi;
|
---|
432 | final double rh = hypot(x, y);
|
---|
433 | phi = rh * 0.5;
|
---|
434 | if (phi > 1.0) {
|
---|
435 | throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
|
---|
436 | }
|
---|
437 | phi = 2.0 * asin(phi);
|
---|
438 | switch (mode) {
|
---|
439 | case EQUATORIAL: {
|
---|
440 | final double sinz = sin(phi);
|
---|
441 | final double cosz = cos(phi);
|
---|
442 | phi = abs(rh) <= FINE_EPSILON ? 0.0 : asin(y * sinz / rh);
|
---|
443 | x *= sinz;
|
---|
444 | y = cosz * rh;
|
---|
445 | lambda = (y == 0) ? 0.0 : atan2(x, y);
|
---|
446 | break;
|
---|
447 | }
|
---|
448 | case OBLIQUE: {
|
---|
449 | final double sinz = sin(phi);
|
---|
450 | final double cosz = cos(phi);
|
---|
451 | phi = abs(rh) <= FINE_EPSILON ? latitudeOfOrigin :
|
---|
452 | asin(cosz * sinb1 + y * sinz * cosb1 / rh);
|
---|
453 | x *= sinz * cosb1;
|
---|
454 | y = (cosz - sin(phi) * sinb1) * rh;
|
---|
455 | lambda = (y == 0) ? 0.0 : atan2(x, y);
|
---|
456 | break;
|
---|
457 | }
|
---|
458 | case NORTH_POLE: {
|
---|
459 | phi = (PI / 2) - phi;
|
---|
460 | lambda = atan2(x, -y);
|
---|
461 | break;
|
---|
462 | }
|
---|
463 | case SOUTH_POLE: {
|
---|
464 | phi -= (PI / 2);
|
---|
465 | lambda = atan2(x, y);
|
---|
466 | break;
|
---|
467 | }
|
---|
468 | default: {
|
---|
469 | throw new AssertionError(mode);
|
---|
470 | }
|
---|
471 | }
|
---|
472 | assert checkInverseTransform(lambda, phi, ptDst);
|
---|
473 | if (ptDst != null) {
|
---|
474 | ptDst.setLocation(lambda, phi);
|
---|
475 | return ptDst;
|
---|
476 | }
|
---|
477 | return new Point2D.Double(lambda, phi);
|
---|
478 | }
|
---|
479 | }
|
---|
480 |
|
---|
481 | /**
|
---|
482 | * Calculates <var>q</var>, Snyder equation (3-12)
|
---|
483 | *
|
---|
484 | * @param sinphi sin of the latitude <var>q</var> is calculated for.
|
---|
485 | * @return <var>q</var> from Snyder equation (3-12).
|
---|
486 | */
|
---|
487 | private double qsfn(final double sinphi) {
|
---|
488 | if (excentricity >= EPSILON) {
|
---|
489 | final double con = excentricity * sinphi;
|
---|
490 | return ((1.0 - excentricitySquared) * (sinphi / (1.0 - con*con) -
|
---|
491 | (0.5 / excentricity) * log((1.0 - con) / (1.0 + con))));
|
---|
492 | } else {
|
---|
493 | return sinphi + sinphi;
|
---|
494 | }
|
---|
495 | }
|
---|
496 |
|
---|
497 | /**
|
---|
498 | * Determines latitude from authalic latitude.
|
---|
499 | */
|
---|
500 | private double authlat(final double beta) {
|
---|
501 | final double t = beta + beta;
|
---|
502 | return beta + APA0 * sin(t) + APA1 * sin(t+t) + APA2 * sin(t+t+t);
|
---|
503 | }
|
---|
504 |
|
---|
505 |
|
---|
506 |
|
---|
507 |
|
---|
508 | //////////////////////////////////////////////////////////////////////////////////////////
|
---|
509 | //////////////////////////////////////////////////////////////////////////////////////////
|
---|
510 | //////// ////////
|
---|
511 | //////// PROVIDERS ////////
|
---|
512 | //////// ////////
|
---|
513 | //////////////////////////////////////////////////////////////////////////////////////////
|
---|
514 | //////////////////////////////////////////////////////////////////////////////////////////
|
---|
515 |
|
---|
516 | /**
|
---|
517 | * The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
|
---|
518 | * provider} for an {@linkplain LambertAzimuthalEqualArea Lambert Equal Area} projection
|
---|
519 | * (EPSG code 9820).
|
---|
520 | *
|
---|
521 | * @since 2.4
|
---|
522 | * @version $Id: LambertAzimuthalEqualArea.java 37299 2011-05-25 05:21:24Z mbedward $
|
---|
523 | * @author Beate Stollberg
|
---|
524 | *
|
---|
525 | * @see org.geotools.referencing.operation.DefaultMathTransformFactory
|
---|
526 | */
|
---|
527 | public static class Provider extends AbstractProvider {
|
---|
528 | /**
|
---|
529 | * For cross-version compatibility.
|
---|
530 | */
|
---|
531 | private static final long serialVersionUID = 3877793025552244132L;
|
---|
532 |
|
---|
533 | /**
|
---|
534 | * The operation parameter descriptor for the {@link #latitudeOfOrigin}
|
---|
535 | * parameter value. Valid values range is from -90 to 90°. Default value is 0.
|
---|
536 | */
|
---|
537 | public static final ParameterDescriptor LATITUDE_OF_CENTRE = createDescriptor(
|
---|
538 | new NamedIdentifier[] {
|
---|
539 | new NamedIdentifier(Citations.OGC, "latitude_of_center"),
|
---|
540 | new NamedIdentifier(Citations.EPSG, "Latitude of natural origin"),
|
---|
541 | new NamedIdentifier(Citations.EPSG, "Spherical latitude of origin"),
|
---|
542 | new NamedIdentifier(Citations.ESRI, "Latitude_Of_Origin"),
|
---|
543 | new NamedIdentifier(Citations.GEOTIFF, "ProjCenterLat")
|
---|
544 | },
|
---|
545 | 0, -90, 90, NonSI.DEGREE_ANGLE);
|
---|
546 |
|
---|
547 | /**
|
---|
548 | * The operation parameter descriptor for the {@link #centralMeridian}
|
---|
549 | * parameter value. Valid values range is from -180 to 180°. Default value is 0.
|
---|
550 | */
|
---|
551 | public static final ParameterDescriptor LONGITUDE_OF_CENTRE = createDescriptor(
|
---|
552 | new NamedIdentifier[] {
|
---|
553 | new NamedIdentifier(Citations.OGC, "longitude_of_center"),
|
---|
554 | new NamedIdentifier(Citations.EPSG, "Longitude of natural origin"),
|
---|
555 | new NamedIdentifier(Citations.EPSG, "Spherical longitude of origin"),
|
---|
556 | new NamedIdentifier(Citations.ESRI, "Central_Meridian"),
|
---|
557 | new NamedIdentifier(Citations.GEOTIFF, "ProjCenterLong")
|
---|
558 | },
|
---|
559 | 0, -180, 180, NonSI.DEGREE_ANGLE);
|
---|
560 |
|
---|
561 | /**
|
---|
562 | * The parameters group.
|
---|
563 | */
|
---|
564 | static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(new NamedIdentifier[] {
|
---|
565 | new NamedIdentifier(Citations.OGC, "Lambert_Azimuthal_Equal_Area"),
|
---|
566 | new NamedIdentifier(Citations.EPSG, "Lambert Azimuthal Equal Area"),
|
---|
567 | new NamedIdentifier(Citations.EPSG, "Lambert Azimuthal Equal Area (Spherical)"),
|
---|
568 | new NamedIdentifier(Citations.GEOTIFF, "CT_LambertAzimEqualArea"),
|
---|
569 | new NamedIdentifier(Citations.EPSG, "9820"),
|
---|
570 | }, new ParameterDescriptor[] {
|
---|
571 | SEMI_MAJOR, SEMI_MINOR,
|
---|
572 | LATITUDE_OF_CENTRE, LONGITUDE_OF_CENTRE,
|
---|
573 | FALSE_EASTING, FALSE_NORTHING
|
---|
574 | });
|
---|
575 |
|
---|
576 | /**
|
---|
577 | * Constructs a new provider.
|
---|
578 | */
|
---|
579 | public Provider() {
|
---|
580 | super(PARAMETERS);
|
---|
581 | }
|
---|
582 |
|
---|
583 | /**
|
---|
584 | * Creates a transform from the specified group of parameter values.
|
---|
585 | *
|
---|
586 | * @param parameters The group of parameter values.
|
---|
587 | * @return The created math transform.
|
---|
588 | * @throws ParameterNotFoundException if a required parameter was not found.
|
---|
589 | */
|
---|
590 | public MathTransform createMathTransform(final ParameterValueGroup parameters)
|
---|
591 | throws ParameterNotFoundException
|
---|
592 | {
|
---|
593 | return isSpherical(parameters) ? new Spherical(parameters) :
|
---|
594 | new LambertAzimuthalEqualArea(parameters);
|
---|
595 | }
|
---|
596 | }
|
---|
597 | }
|
---|