001: //$HeadURL: $
002: /*---------------- FILE HEADER ------------------------------------------
003: This file is part of deegree.
004: Copyright (C) 2001-2008 by:
005: Department of Geography, University of Bonn
006: http://www.giub.uni-bonn.de/deegree/
007: lat/lon GmbH
008: http://www.lat-lon.de
009:
010: This library is free software; you can redistribute it and/or
011: modify it under the terms of the GNU Lesser General Public
012: License as published by the Free Software Foundation; either
013: version 2.1 of the License, or (at your option) any later version.
014: This library is distributed in the hope that it will be useful,
015: but WITHOUT ANY WARRANTY; without even the implied warranty of
016: MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
017: Lesser General Public License for more details.
018: You should have received a copy of the GNU Lesser General Public
019: License along with this library; if not, write to the Free Software
020: Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
021: Contact:
022:
023: Andreas Poth
024: lat/lon GmbH
025: Aennchenstr. 19
026: 53177 Bonn
027: Germany
028: E-Mail: poth@lat-lon.de
029:
030: Prof. Dr. Klaus Greve
031: Department of Geography
032: University of Bonn
033: Meckenheimer Allee 166
034: 53115 Bonn
035: Germany
036: E-Mail: greve@giub.uni-bonn.de
037: ---------------------------------------------------------------------------*/
038:
039: package org.deegree.crs.projections.azimuthal;
040:
041: import javax.vecmath.Point2d;
042:
043: import org.deegree.crs.components.Unit;
044: import org.deegree.crs.coordinatesystems.GeographicCRS;
045: import org.deegree.crs.exceptions.ProjectionException;
046: import org.deegree.crs.projections.ProjectionUtils;
047:
048: /**
049: * The <code>StereographicAzimuthal</code> class allows for Stereographic Projections of the Poles, equator as well as
050: * oblique. This projection has following properties (Snyder p. 154):
051: * <ul>
052: * <li>Azimuthal</li>
053: * <li>Conformal</li>
054: * <li>The central meridian an a particular parallel (if shown) are straight lines.</li>
055: * <li>All meridians on the polar aspects and the equator on the equatorial aspect are straight lines.</li>
056: * <li>All other meidians and parallels are shown as arcs of circles</li>
057: * <li>A perspective projections for the sphere</li>
058: * <li>Directions from the center of the projection are true (except on ellipsoidal oblique and equatorial aspects)</li>
059: * <li>Scale increases away from the center of the projection</li>
060: * <li>Point opposite the center of the projection cannot be plotted</li>
061: * <li>Used for polar maps and miscellaneous special maps</li>
062: * </ul>
063: *
064: * <p>
065: * Like Orthographic, the stereographic projection is a true perspective in its isSpherical() form. It is the only known
066: * true perspective projection of any kind that is also conformal. Its point of projection is on the the surface of the
067: * sphere at a point jus opposite the oint of tangency of the plane or the center point of the projection. Thus, if the
068: * north pole is the center of the map, the projection is from the south-pole.
069: * </p>
070: * <p>
071: * It is known to be used by following epsg transformations:
072: * <ul>
073: * <li>EPSG:28992</li>
074: * </ul>
075: * </p>
076: *
077: * @author <a href="mailto:bezema@lat-lon.de">Rutger Bezema</a>
078: *
079: * @author last edited by: $Author:$
080: *
081: * @version $Revision:$, $Date:$
082: *
083: */
084:
085: public class StereographicAzimuthal extends AzimuthalProjection {
086:
087: /**
088: * This variable will hold different values, for the ellipsoidal projection:
089: * <ul>
090: * <li>for Oblique projection it is the first Part (2*a*k_0*m_1, hence it's name) of A (p.160 21-27)</li>
091: * <li>for the north-south it may hold the rho of (21-33) or (21-34)
092: * <li>for the equatorial it holds 2*scale.</li>
093: * <li>For the
094: * </ul>
095: * For the spherical projection:
096: * <ul>
097: * <li>For oblique and equatorial: 2*scale, which is compliant with 2*k_0 from Snyder (p.157 21-4) </li>
098: * <li>For north and south: cos(truelat) / tan( pi/4 - phi/2) ???? </li>
099: * </ul>
100: */
101: private double akm1;
102:
103: private double trueScaleLatitude;
104:
105: private double[] preCalcedPhiSeries;
106:
107: // use for the OBLIQUE projection instead of the projectionLatitude
108: private double conformalLatitude = 0;
109:
110: // sine of the conformal latitude
111: private double sinCL = 0;
112:
113: // cosine of the conformal latitude
114: private double cosCL = 0;
115:
116: /**
117: * @param trueScaleLatitude
118: * the latitude (in radians) of a circle around the projection point, which contains the true scale.
119: * @param geographicCRS
120: * @param falseNorthing
121: * @param falseEasting
122: * @param naturalOrigin
123: * @param units
124: * @param scale
125: * @param identifiers
126: * @param names
127: * @param versions
128: * @param descriptions
129: * @param areasOfUse
130: */
131: public StereographicAzimuthal(double trueScaleLatitude,
132: GeographicCRS geographicCRS, double falseNorthing,
133: double falseEasting, Point2d naturalOrigin, Unit units,
134: double scale, String[] identifiers, String[] names,
135: String[] versions, String[] descriptions,
136: String[] areasOfUse) {
137: super (geographicCRS, falseNorthing, falseEasting,
138: naturalOrigin, units, scale, true, false,// not equal area.
139: identifiers, names, versions, descriptions, areasOfUse);
140:
141: preCalcedPhiSeries = ProjectionUtils
142: .preCalcedThetaSeries(getSquaredEccentricity());
143:
144: this .trueScaleLatitude = Math.abs(trueScaleLatitude);
145: if (!isSpherical()) {
146: // double X;
147: double tmp = 0;
148: switch (getMode()) {
149: case NORTH_POLE:
150: case SOUTH_POLE:
151: /**
152: * The t from 21-33 and 21-34 is still to be calculated.
153: */
154: // If true scale or known scale factor k_0 is to occur at the pole, the following applies:
155: if (Math.abs(this .trueScaleLatitude
156: - ProjectionUtils.HALFPI) < ProjectionUtils.EPS10) {
157: // Snyder (p.161 21-33) akm1 = rho.
158: akm1 = 2.
159: * getScaleFactor()
160: / Math.sqrt(Math.pow(1 + getEccentricity(),
161: 1 + getEccentricity())
162: * Math.pow(1 - getEccentricity(),
163: 1 - getEccentricity()));
164:
165: } else {
166: // For true scale along the circle represented by trueScaleLatitude (phi_c in snyder)..
167: // akm1 will hold m_c/t_c
168: // Calculate the rho from Snyder (p.161 21-34) and place in akm1
169: tmp = Math.sin(this .trueScaleLatitude);
170:
171: // First part of m of Snyder (p.160 14-15)
172: akm1 = Math.cos(this .trueScaleLatitude)
173: / ProjectionUtils.tanHalfCoLatitude(
174: this .trueScaleLatitude, tmp,
175: getEccentricity());
176: tmp *= getEccentricity();
177: // Second part of m of Snyder (p.160 14-15), t = (e*sin(phi_c))
178: akm1 /= Math.sqrt(1. - tmp * tmp);
179:
180: }
181: break;
182: case EQUATOR:
183: akm1 = 2. * getScaleFactor();
184: break;
185: case OBLIQUE:
186: tmp = getSinphi0();// Math.sin( getProjectionLatitude() );
187: // Calculate the ConformalLatitude for this ellipsoid Snyder (p.160 3-1) the X
188: conformalLatitude = (2. * Math.atan(ProjectionUtils
189: .conformalLatitudeInnerPart(
190: getProjectionLatitude(), tmp,
191: getEccentricity())))
192: - ProjectionUtils.HALFPI;
193:
194: sinCL = Math.sin(conformalLatitude);
195: cosCL = Math.cos(conformalLatitude);
196:
197: tmp *= getEccentricity();
198: // the first part (2a*k_0*m) of the largeA in Snyder (p.160 21-27) is stored in akm1 it still must be
199: // devided by the cos X ... etc. to get largeA
200: akm1 = 2. * getScaleFactor() * getCosphi0()
201: / Math.sqrt(1. - tmp * tmp);
202:
203: // Setting the sinPhi0 to the conformal Latitude X.
204: // setProjectionLatitude( X );
205: break;
206: }
207: } else {
208: akm1 = 2. * getScaleFactor();
209: }
210: }
211:
212: /**
213: * @param trueScaleLatitude
214: * the latitude (in radians) of a circle around the projection point, which contains the true scale.
215: * @param geographicCRS
216: * @param falseNorthing
217: * @param falseEasting
218: * @param naturalOrigin
219: * @param units
220: * @param scale
221: * @param identifier
222: * @param name
223: * @param version
224: * @param description
225: * @param areaOfUse
226: */
227: public StereographicAzimuthal(double trueScaleLatitude,
228: GeographicCRS geographicCRS, double falseNorthing,
229: double falseEasting, Point2d naturalOrigin, Unit units,
230: double scale, String identifier, String name,
231: String version, String description, String areaOfUse) {
232: this (trueScaleLatitude, geographicCRS, falseNorthing,
233: falseEasting, naturalOrigin, units, scale,
234: new String[] { identifier }, new String[] { name },
235: new String[] { version }, new String[] { description },
236: new String[] { areaOfUse });
237: }
238:
239: /**
240: * @param trueScaleLatitude
241: * the latitude (in radians) of a circle around the projection point, which contains the true scale.
242: * @param geographicCRS
243: * @param falseNorthing
244: * @param falseEasting
245: * @param naturalOrigin
246: * @param units
247: * @param scale
248: * @param identifiers
249: */
250: public StereographicAzimuthal(double trueScaleLatitude,
251: GeographicCRS geographicCRS, double falseNorthing,
252: double falseEasting, Point2d naturalOrigin, Unit units,
253: double scale, String[] identifiers) {
254: this (trueScaleLatitude, geographicCRS, falseNorthing,
255: falseEasting, naturalOrigin, units, scale, identifiers,
256: null, null, null, null);
257: }
258:
259: /**
260: * @param trueScaleLatitude
261: * the latitude (in radians) of a circle around the projection point, which contains the true scale.
262: * @param geographicCRS
263: * @param falseNorthing
264: * @param falseEasting
265: * @param naturalOrigin
266: * @param units
267: * @param scale
268: * @param identifier
269: */
270: public StereographicAzimuthal(double trueScaleLatitude,
271: GeographicCRS geographicCRS, double falseNorthing,
272: double falseEasting, Point2d naturalOrigin, Unit units,
273: double scale, String identifier) {
274: this (trueScaleLatitude, geographicCRS, falseNorthing,
275: falseEasting, naturalOrigin, units, scale,
276: new String[] { identifier });
277: }
278:
279: /**
280: * Create a {@link StereographicAzimuthal} which has a truescale latitude at MapUtils.HALFPI.
281: *
282: * @param geographicCRS
283: * @param falseNorthing
284: * @param falseEasting
285: * @param naturalOrigin
286: * @param units
287: * @param scale
288: * @param identifiers
289: */
290: public StereographicAzimuthal(GeographicCRS geographicCRS,
291: double falseNorthing, double falseEasting,
292: Point2d naturalOrigin, Unit units, double scale,
293: String[] identifiers) {
294: this (ProjectionUtils.HALFPI, geographicCRS, falseNorthing,
295: falseEasting, naturalOrigin, units, scale, identifiers);
296: }
297:
298: /**
299: * Create a {@link StereographicAzimuthal} which has a truescale latitude at MapUtils.HALFPI.
300: *
301: * @param geographicCRS
302: * @param falseNorthing
303: * @param falseEasting
304: * @param naturalOrigin
305: * @param units
306: * @param scale
307: * @param identifier
308: */
309: public StereographicAzimuthal(GeographicCRS geographicCRS,
310: double falseNorthing, double falseEasting,
311: Point2d naturalOrigin, Unit units, double scale,
312: String identifier) {
313: this (ProjectionUtils.HALFPI, geographicCRS, falseNorthing,
314: falseEasting, naturalOrigin, units, scale, identifier);
315: }
316:
317: /**
318: * Create a {@link StereographicAzimuthal} which has a scale of 1 and a truescale latitude,
319: *
320: * @param trueScaleLatitude
321: * the latitude (in radians) of a circle around the projection point, which contains the true scale.
322: * @param geographicCRS
323: * @param falseNorthing
324: * @param falseEasting
325: * @param naturalOrigin
326: * @param units
327: * @param identifiers
328: */
329: public StereographicAzimuthal(double trueScaleLatitude,
330: GeographicCRS geographicCRS, double falseNorthing,
331: double falseEasting, Point2d naturalOrigin, Unit units,
332: String[] identifiers) {
333: this (trueScaleLatitude, geographicCRS, falseNorthing,
334: falseEasting, naturalOrigin, units, 1, identifiers);
335: }
336:
337: /**
338: * Create a {@link StereographicAzimuthal} which has a scale of 1 and a truescale latitude,
339: *
340: * @param trueScaleLatitude
341: * the latitude (in radians) of a circle around the projection point, which contains the true scale.
342: * @param geographicCRS
343: * @param falseNorthing
344: * @param falseEasting
345: * @param naturalOrigin
346: * @param units
347: * @param identifier
348: */
349: public StereographicAzimuthal(double trueScaleLatitude,
350: GeographicCRS geographicCRS, double falseNorthing,
351: double falseEasting, Point2d naturalOrigin, Unit units,
352: String identifier) {
353: this (trueScaleLatitude, geographicCRS, falseNorthing,
354: falseEasting, naturalOrigin, units, 1, identifier);
355: }
356:
357: /**
358: * Create a {@link StereographicAzimuthal} which is conformal, has a scale of 1 and a truescale latitude at pi*0.5.
359: *
360: * @param geographicCRS
361: * @param falseNorthing
362: * @param falseEasting
363: * @param naturalOrigin
364: * @param units
365: * @param identifiers
366: */
367: public StereographicAzimuthal(GeographicCRS geographicCRS,
368: double falseNorthing, double falseEasting,
369: Point2d naturalOrigin, Unit units, String[] identifiers) {
370: this (ProjectionUtils.HALFPI, geographicCRS, falseNorthing,
371: falseEasting, naturalOrigin, units, 1, identifiers);
372: }
373:
374: /**
375: * Create a {@link StereographicAzimuthal} which is conformal, has a scale of 1 and a truescale latitude at pi*0.5.
376: *
377: * @param geographicCRS
378: * @param falseNorthing
379: * @param falseEasting
380: * @param naturalOrigin
381: * @param units
382: * @param identifier
383: */
384: public StereographicAzimuthal(GeographicCRS geographicCRS,
385: double falseNorthing, double falseEasting,
386: Point2d naturalOrigin, Unit units, String identifier) {
387: this (ProjectionUtils.HALFPI, geographicCRS, falseNorthing,
388: falseEasting, naturalOrigin, units, 1, identifier);
389: }
390:
391: /*
392: * (non-Javadoc)
393: *
394: * @see org.deegree.crs.projections.Projection#doInverseProjection(double, double)
395: */
396: @Override
397: public Point2d doInverseProjection(double x, double y) {
398: Point2d lp = new Point2d(0, 0);
399: x -= getFalseEasting();
400: y -= getFalseNorthing();
401: double rho = ProjectionUtils.length(x, y);
402: if (isSpherical()) {
403: // akm1 = 2*scale.
404: double c = 2. * Math.atan(rho / akm1);
405: double sinc = Math.sin(c);
406: double cosc = Math.cos(c);
407: switch (getMode()) {
408: case OBLIQUE:
409: // First find theta from Snyder (p.158 20-14).
410: if (Math.abs(rho) <= ProjectionUtils.EPS10) {
411: lp.y = getProjectionLatitude();
412: } else {
413: lp.y = Math.asin(cosc * getSinphi0()
414: + ((y * sinc * getCosphi0()) / rho));
415: }
416: // And now lambda but instead of using the atan, use atan2. Is this correct????
417: c = cosc - getSinphi0() * Math.sin(lp.y);
418: if (Math.abs(c) > ProjectionUtils.EPS10 || x != 0.) {
419: lp.x = Math.atan2(x * sinc * getCosphi0(), c * rho);
420: }
421: break;
422: case EQUATOR:
423: if (Math.abs(rho) > ProjectionUtils.EPS10) {
424: lp.y = Math.asin(y * sinc / rho);
425: }
426: if (cosc != 0. || x != 0.) {
427: lp.x = Math.atan2(x * sinc, cosc * rho);
428: }
429: break;
430: case NORTH_POLE:
431: y = -y;
432: case SOUTH_POLE:
433: if (Math.abs(rho) <= ProjectionUtils.EPS10) {
434: lp.y = getProjectionLatitude();
435: } else {
436: lp.y = Math.asin(getMode() == SOUTH_POLE ? -cosc
437: : cosc);
438: }
439: lp.x = Math.atan2(x, y);
440: break;
441: }
442: } else {
443: double cosCe = 0;
444: double sinCe = 0;
445: double X = 0;
446: // for oblique and equatorial projections, akm1 = first part of Snyder (p.160 21-27). which is 2*scale*m_1
447: // double ce = 2. * Math.atan2( rho * getCosphi0(), akm1 );
448: // if ( getMode() == OBLIQUE || getMode() == EQUATOR ) {
449: // cosCe = Math.cos( ce );
450: // sinCe = Math.sin( ce );
451: // }
452: double ce = 0;
453: switch (getMode()) {
454: case OBLIQUE:
455: ce = 2. * Math.atan2(rho * cosCL, akm1);
456: cosCe = Math.cos(ce);
457: sinCe = Math.sin(ce);
458: X = Math.asin(cosCe * sinCL
459: + ((y * sinCe * cosCL) / rho));
460: lp.x = Math.atan2(x * sinCe, rho * cosCL * cosCe - y
461: * sinCL * sinCe);
462: break;
463: case EQUATOR:
464: ce = 2. * Math.atan2(rho * getCosphi0(), akm1);
465: cosCe = Math.cos(ce);
466: sinCe = Math.sin(ce);
467: X = Math.asin(cosCe * getSinphi0()
468: + ((y * sinCe * getCosphi0()) / rho));
469: lp.x = Math.atan2(x * sinCe, rho * getCosphi0() * cosCe
470: - y * getSinphi0() * sinCe);
471: break;
472: case NORTH_POLE:
473: y = -y;
474: case SOUTH_POLE:
475: /**
476: * akm1 can hold rho from 21-34 or 21-33 (if no truescale latitude was given) without the parameter t.
477: * Either way it must only be multiplied with t to forfill the formula. Independent of the true scale
478: * latitude, the value of akm1 must therefore only be inverted to calculate Snyder (p.162 21-39 or
479: * 21-40).
480: */
481: double t = rho * 1. / this .akm1;
482: X = ProjectionUtils.HALFPI - 2 * Math.atan(t);
483: lp.x = Math.atan2(x, y);
484: break;
485: }
486: lp.y = ProjectionUtils.calcPhiFromConformalLatitude(X,
487: preCalcedPhiSeries);
488: lp.x += getProjectionLongitude();
489: }
490: return lp;
491: }
492:
493: /*
494: * (non-Javadoc)
495: *
496: * @see org.deegree.crs.projections.Projection#doProjection(double, double)
497: */
498: @Override
499: public Point2d doProjection(double lambda, double phi)
500: throws ProjectionException {
501: Point2d xy = new Point2d(0, 0);
502: lambda -= getProjectionLongitude();
503: double cosLamda = Math.cos(lambda);
504: double sinLamda = Math.sin(lambda);
505: double sinPhi = Math.sin(phi);
506:
507: if (isSpherical()) {
508: double cosphi = Math.cos(phi);
509:
510: switch (getMode()) {
511: case OBLIQUE:
512: // Calc k from Snyder (p.158 21-14) and store in xy.y
513: xy.y = 1. + getSinphi0() * sinPhi + getCosphi0()
514: * cosphi * cosLamda;
515: if (xy.y <= ProjectionUtils.EPS11) {
516: throw new ProjectionException("Division by zero");
517: }
518: // akm1 was set to 2*scale.
519: xy.y = akm1 / xy.y;
520:
521: // Calc x (p.157 21-2) and y (p.157 21-3)
522: xy.x = xy.y * cosphi * sinLamda;
523: xy.y *= getCosphi0() * sinPhi - getSinphi0() * cosphi
524: * cosLamda;
525: break;
526: case EQUATOR:
527: // Calc k from Snyder (p.158 21-14) and store in xy.y
528: xy.y = 1. + cosphi * cosLamda;
529: if (xy.y <= ProjectionUtils.EPS11) {
530: throw new ProjectionException("Division by zero");
531: }
532: // akm1 was set to 2*scale.
533: xy.y = akm1 / xy.y;
534:
535: // Calc x (p.157 21-2) and y (p.158 21-13)
536: xy.x = xy.y * cosphi * sinLamda;
537: xy.y *= sinPhi;
538: break;
539: case NORTH_POLE:
540: cosLamda = -cosLamda;
541: phi = -phi;
542: case SOUTH_POLE:
543: if (Math.abs(phi - ProjectionUtils.HALFPI) < ProjectionUtils.EPS10) {
544: throw new ProjectionException("The requested phi: "
545: + phi
546: + " lies on the singularity ("
547: + ((getMode() == SOUTH_POLE) ? "South-Pole"
548: : "North-Pole")
549: + ") of this projection's mappable area.");
550: }
551: // akm1 was set to 2*scale.
552: xy.y = akm1
553: * Math
554: .tan(ProjectionUtils.QUARTERPI + .5
555: * phi);
556: xy.x = sinLamda * (xy.y );
557: xy.y *= cosLamda;
558: break;
559: }
560: } else {
561: double sinX = 0;
562: double cosX = 0;
563: double largeA;
564:
565: if (getMode() == OBLIQUE || getMode() == EQUATOR) {
566: // Calculate the conformal latitue Snyder (p.160 3-1).
567: double X = 2.
568: * Math.atan(ProjectionUtils
569: .conformalLatitudeInnerPart(phi,
570: sinPhi, getEccentricity()))
571: - ProjectionUtils.HALFPI;
572: sinX = Math.sin(X);
573: cosX = Math.cos(X);
574: }
575: switch (getMode()) {
576: case OBLIQUE:
577: // System.out.println( "Oblique stereographic projection" );
578: // akm1 was set to the first part of Snyder (p.160 21-27)
579: // sinPhi0 and cosPhi0 where set to the conformal latitude of the natural origin.
580: largeA = akm1
581: / (cosCL * (1. + sinCL * sinX + cosCL * cosX
582: * cosLamda));
583: xy.y = largeA
584: * (cosCL * sinX - sinCL * cosX * cosLamda);
585: xy.x = largeA * cosX;
586: break;
587: case EQUATOR:
588: // System.out.println( "equatorial stereographic projection" );
589: // akm1 was set to 2. * scale;
590: // Calculate the A of Snyder (p.161 21-29).
591: largeA = 2. * akm1 / (1. + cosX * cosLamda);
592: xy.y = largeA * sinX;
593: xy.x = largeA * cosX;
594: break;
595: case SOUTH_POLE:
596: // System.out.println( "South-pole stereographic projection" );
597: phi = -phi;
598: cosLamda = -cosLamda;
599: sinPhi = -sinPhi;
600: case NORTH_POLE:
601: // System.out.println( "North-pole stereographic projection" );
602: // akm1 holds the rho's from 21-33 or 21-34, but without the t (==halfColatitude of the conformal
603: // latitude).
604: xy.x = akm1
605: * ProjectionUtils.tanHalfCoLatitude(phi,
606: sinPhi, getEccentricity());
607: xy.y = -xy.x * cosLamda;
608: break;
609: }
610: // Sinus(Lambda) was still to be multiplied on the x.
611: xy.x = xy.x * sinLamda;
612: }
613: xy.x += getFalseEasting();
614: xy.y += getFalseNorthing();
615: return xy;
616: }
617:
618: @Override
619: public String getDeegreeSpecificName() {
620: return "stereographicAzimuthal";
621: }
622:
623: /**
624: * @return the trueScaleLatitude.
625: */
626: public final double getTrueScaleLatitude() {
627: return trueScaleLatitude;
628: }
629: }
|