UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
PolarStereographic.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: POLAR STEREOGRAPHIC
5  *
6  *
7  * ABSTRACT
8  *
9  * This component provides conversions between geodetic (latitude and
10  * longitude) coordinates and Polar Stereographic (easting and northing)
11  * coordinates.
12  *
13  * ERROR HANDLING
14  *
15  * This component checks parameters for valid values. If an invalid
16  * value is found the error code is combined with the current error code
17  * using the bitwise or. This combining allows multiple error codes to
18  * be returned. The possible error codes are:
19  *
20  * POLAR_NO_ERROR : No errors occurred in function
21  * POLAR_LAT_ERROR : Latitude outside of valid range
22  * (-90 to 90 degrees)
23  * POLAR_LON_ERROR : Longitude outside of valid range
24  * (-180 to 360 degrees)
25  * POLAR_ORIGIN_LAT_ERROR : Latitude of true scale outside of valid
26  * range (-90 to 90 degrees)
27  * POLAR_ORIGIN_LON_ERROR : Longitude down from pole outside of valid
28  * range (-180 to 360 degrees)
29  * POLAR_EASTING_ERROR : Easting outside of valid range,
30  * depending on ellipsoid and
31  * projection parameters
32  * POLAR_NORTHING_ERROR : Northing outside of valid range,
33  * depending on ellipsoid and
34  * projection parameters
35  * POLAR_RADIUS_ERROR : Coordinates too far from pole,
36  * depending on ellipsoid and
37  * projection parameters
38  * POLAR_A_ERROR : Semi-major axis less than or equal to zero
39  * POLAR_INV_F_ERROR : Inverse flattening outside of valid range
40  * (250 to 350)
41  *
42  *
43  * REUSE NOTES
44  *
45  * POLAR STEREOGRAPHIC is intended for reuse by any application that
46  * performs a Polar Stereographic projection.
47  *
48  *
49  * REFERENCES
50  *
51  * Further information on POLAR STEREOGRAPHIC can be found in the
52  * Reuse Manual.
53  *
54  *
55  * POLAR STEREOGRAPHIC originated from :
56  * U.S. Army Topographic Engineering Center
57  * Geospatial Information Division
58  * 7701 Telegraph Road
59  * Alexandria, VA 22310-3864
60  *
61  *
62  * LICENSES
63  *
64  * None apply to this component.
65  *
66  *
67  * RESTRICTIONS
68  *
69  * POLAR STEREOGRAPHIC has no restrictions.
70  *
71  *
72  * ENVIRONMENT
73  *
74  * POLAR STEREOGRAPHIC was tested and certified in the following
75  * environments:
76  *
77  * 1. Solaris 2.5 with GCC, version 2.8.1
78  * 2. Window 95 with MS Visual C++, version 6
79  *
80  *
81  * MODIFICATIONS
82  *
83  * Date Description
84  * ---- -----------
85  * 2-27-07 Original Code
86  *
87  *
88  */
89 
90 
91 /************************************************************************/
92 /*
93  * INCLUDES
94  */
95 
96 #include <stdio.h>
97 #include <math.h>
98 #include "PolarStereographic.h"
99 #include "PolarStereographicStandardParallelParameters.h"
100 #include "PolarStereographicScaleFactorParameters.h"
101 #include "MapProjectionCoordinates.h"
102 #include "GeodeticCoordinates.h"
103 #include "CoordinateConversionException.h"
104 #include "ErrorMessages.h"
105 
106 /*
107  * math.h - Standard C math library
108  * PolarStereographic.h - Is for prototype error checking
109  * MapProjectionCoordinates.h - defines map projection coordinates
110  * GeodeticCoordinates.h - defines geodetic coordinates
111  * CoordinateConversionException.h - Exception handler
112  * ErrorMessages.h - Contains exception messages
113  */
114 
115 
116 using namespace MSP::CCS;
117 
118 
119 /************************************************************************/
120 /* DEFINES
121  *
122  */
123 
124 const double PI = 3.14159265358979323e0; /* PI */
125 const double PI_OVER_2 = (PI / 2.0);
126 const double PI_OVER_4 = (PI / 4.0);
127 const double TWO_PI = (2.0 * PI);
128 
129 #define MIN_SCALE_FACTOR 0.1
130 #define MAX_SCALE_FACTOR 3.0
131 
132 
133 /************************************************************************/
134 /* FUNCTIONS
135  *
136  */
137 
138 PolarStereographic::PolarStereographic(
139  double ellipsoidSemiMajorAxis,
140  double ellipsoidFlattening,
141  double centralMeridian,
142  double standardParallel,
143  double falseEasting,
144  double falseNorthing ) :
145  CoordinateSystem(),
146  coordinateType( CoordinateType::polarStereographicStandardParallel ),
147  es( 0.08181919084262188000 ),
148  es_OVER_2( .040909595421311 ),
149  Southern_Hemisphere( 0 ),
150  Polar_tc( 1.0 ),
151  Polar_k90( 1.0033565552493 ),
152  Polar_a_mc( 6378137.0 ),
153  two_Polar_a( 12756274.0 ),
154  Polar_Central_Meridian( 0.0 ),
155  Polar_Standard_Parallel( ((PI * 90) / 180) ),
156  Polar_False_Easting( 0.0 ),
157  Polar_False_Northing( 0.0 ),
158  Polar_Scale_Factor( 1.0 ),
159  Polar_Delta_Easting( 12713601.0 ),
160  Polar_Delta_Northing( 12713601.0 )
161 {
162 /*
163  * The constructor receives the ellipsoid
164  * parameters and Polar Stereograpic (Standard Parallel) projection parameters as inputs, and
165  * sets the corresponding state variables. If any errors occur, an
166  * exception is thrown with a description of the error.
167  *
168  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
169  * ellipsoidFlattening : Flattening of ellipsoid (input)
170  * centralMeridian : Longitude down from pole, in radians (input)
171  * standardParallel : Latitude of true scale, in radians (input)
172  * falseEasting : Easting (X) at center of projection, in meters (input)
173  * falseNorthing : Northing (Y) at center of projection, in meters (input)
174  */
175 
176  double es2;
177  double slat, sinolat, cosolat;
178  double essin;
179  double one_PLUS_es, one_MINUS_es;
180  double one_PLUS_es_sinolat, one_MINUS_es_sinolat;
181  double pow_es;
182  double inv_f = 1 / ellipsoidFlattening;
183 
184  if (ellipsoidSemiMajorAxis <= 0.0)
185  { /* Semi-major axis must be greater than zero */
186  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
187  }
188  if ((inv_f < 250) || (inv_f > 350))
189  { /* Inverse flattening must be between 250 and 350 */
190  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
191  }
192  if ((standardParallel < -PI_OVER_2) || (standardParallel > PI_OVER_2))
193  { /* Origin Latitude out of range */
194  throw CoordinateConversionException( ErrorMessages::originLatitude );
195  }
196  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
197  { /* Origin Longitude out of range */
198  throw CoordinateConversionException( ErrorMessages::centralMeridian );
199  }
200 
201  semiMajorAxis = ellipsoidSemiMajorAxis;
202  flattening = ellipsoidFlattening;
203 
204  two_Polar_a = 2.0 * semiMajorAxis;
205 
206  if (centralMeridian > PI)
207  centralMeridian -= TWO_PI;
208  if (standardParallel < 0)
209  {
210  Southern_Hemisphere = 1;
211  Polar_Standard_Parallel = -standardParallel;
212  Polar_Central_Meridian = -centralMeridian;
213  }
214  else
215  {
216  Southern_Hemisphere = 0;
217  Polar_Standard_Parallel = standardParallel;
218  Polar_Central_Meridian = centralMeridian;
219  }
220  Polar_False_Easting = falseEasting;
221  Polar_False_Northing = falseNorthing;
222 
223  es2 = 2 * flattening - flattening * flattening;
224  es = sqrt(es2);
225  es_OVER_2 = es / 2.0;
226 
227  if (fabs(fabs(Polar_Standard_Parallel) - PI_OVER_2) > 1.0e-10)
228  {
229  sinolat = sin(Polar_Standard_Parallel);
230  essin = es * sinolat;
231  pow_es = polarPow(essin);
232  cosolat = cos(Polar_Standard_Parallel);
233  double mc = cosolat / sqrt(1.0 - essin * essin);
234  Polar_a_mc = semiMajorAxis * mc;
235  Polar_tc = tan(PI_OVER_4 - Polar_Standard_Parallel / 2.0) / pow_es;
236  }
237 
238  one_PLUS_es = 1.0 + es;
239  one_MINUS_es = 1.0 - es;
240  Polar_k90 = sqrt(pow(one_PLUS_es, one_PLUS_es) * pow(one_MINUS_es, one_MINUS_es));
241 
242  slat = sin(fabs(standardParallel));
243  one_PLUS_es_sinolat = 1.0 + es * slat;
244  one_MINUS_es_sinolat = 1.0 - es * slat;
245  Polar_Scale_Factor = ((1 + slat) / 2) * (Polar_k90 / sqrt(pow(one_PLUS_es_sinolat, one_PLUS_es) * pow(one_MINUS_es_sinolat, one_MINUS_es)));
246 
247  /* Calculate Radius */
248  GeodeticCoordinates tempGeodeticCoordinates = GeodeticCoordinates( CoordinateType::geodetic, centralMeridian, 0 );
249  //MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( (GeodeticCoordinates&) GeodeticCoordinates( CoordinateType::geodetic, centralMeridian, 0 ) );
250  MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &tempGeodeticCoordinates );
251  Polar_Delta_Northing = tempCoordinates->northing();
252  delete tempCoordinates;
253 
254  if(Polar_False_Northing)
255  Polar_Delta_Northing -= Polar_False_Northing;
256  if (Polar_Delta_Northing < 0)
257  Polar_Delta_Northing = -Polar_Delta_Northing;
258  Polar_Delta_Northing *= 1.01;
259 
260  Polar_Delta_Easting = Polar_Delta_Northing;
261 }
262 
263 
264 PolarStereographic::PolarStereographic( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double scaleFactor, char hemisphere, double falseEasting, double falseNorthing ) :
265  CoordinateSystem(),
266  coordinateType( CoordinateType::polarStereographicScaleFactor ),
267  es( 0.08181919084262188000 ),
268  es_OVER_2( .040909595421311 ),
269  Southern_Hemisphere( 0 ),
270  Polar_tc( 1.0 ),
271  Polar_k90( 1.0033565552493 ),
272  Polar_a_mc( 6378137.0 ),
273  two_Polar_a( 12756274.0 ),
274  Polar_Central_Meridian( 0.0 ),
275  Polar_Standard_Parallel( ((PI * 90) / 180) ),
276  Polar_False_Easting( 0.0 ),
277  Polar_False_Northing( 0.0 ),
278  Polar_Scale_Factor( 1.0 ),
279  Polar_Delta_Easting( 12713601.0 ),
280  Polar_Delta_Northing( 12713601.0 )
281 {
282 /*
283  * The constructor receives the ellipsoid
284  * parameters and Polar Stereograpic (Scale Factor) projection parameters
285  * as inputs, and sets the corresponding state variables. If any errors occur,
286  * an exception is thrown with a description of the error.
287  *
288  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
289  * ellipsoidFlattening : Flattening of ellipsoid (input)
290  * centralMeridian : Longitude down from pole, in radians (input)
291  * scaleFactor : Scale Factor (input)
292  * falseEasting : Easting (X) at center of projection, in meters (input)
293  * falseNorthing : Northing (Y) at center of projection, in meters (input)
294  */
295 
296  double es2;
297  double sinolat, cosolat;
298  double essin;
299  double pow_es;
300  double mc;
301  double one_PLUS_es, one_MINUS_es;
302  double one_PLUS_es_sk, one_MINUS_es_sk;
303  double sk, sk_PLUS_1;
304  double tolerance = 1.0e-15;
305  int count = 30;
306  double inv_f = 1 / ellipsoidFlattening;
307 
308  if (ellipsoidSemiMajorAxis <= 0.0)
309  { /* Semi-major axis must be greater than zero */
310  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
311  }
312  if ((inv_f < 250) || (inv_f > 350))
313  { /* Inverse flattening must be between 250 and 350 */
314  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
315  }
316  if ((scaleFactor < MIN_SCALE_FACTOR) || (scaleFactor > MAX_SCALE_FACTOR))
317  {
318  throw CoordinateConversionException( ErrorMessages::scaleFactor );
319  }
320  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
321  { /* Origin Longitude out of range */
322  throw CoordinateConversionException( ErrorMessages::centralMeridian );
323  }
324  if ((hemisphere != 'N') && (hemisphere != 'S'))
325  throw CoordinateConversionException( ErrorMessages::hemisphere );
326 
327  semiMajorAxis = ellipsoidSemiMajorAxis;
328  flattening = ellipsoidFlattening;
329  Polar_Scale_Factor = scaleFactor;
330  Polar_False_Easting = falseEasting;
331  Polar_False_Northing = falseNorthing;
332 
333  two_Polar_a = 2.0 * semiMajorAxis;
334  es2 = 2 * flattening - flattening * flattening;
335  es = sqrt(es2);
336  es_OVER_2 = es / 2.0;
337 
338  one_PLUS_es = 1.0 + es;
339  one_MINUS_es = 1.0 - es;
340  Polar_k90 = sqrt(pow(one_PLUS_es, one_PLUS_es) * pow(one_MINUS_es, one_MINUS_es));
341 
342  sk = 0;
343  sk_PLUS_1 = -1 + 2 * Polar_Scale_Factor;
344  while (fabs(sk_PLUS_1 - sk) > tolerance && count)
345  {
346  sk = sk_PLUS_1;
347  one_PLUS_es_sk = 1.0 + es * sk;
348  one_MINUS_es_sk = 1.0 - es * sk;
349  sk_PLUS_1 = ((2 * Polar_Scale_Factor * sqrt(pow(one_PLUS_es_sk, one_PLUS_es) * pow(one_MINUS_es_sk, one_MINUS_es))) / Polar_k90) - 1;
350  count --;
351  }
352 
353  if(!count)
354  throw CoordinateConversionException( ErrorMessages::originLatitude );
355 
356  double standardParallel = 0.0;
357  if(sk_PLUS_1 >= -1.0 && sk_PLUS_1 <= 1.0)
358  standardParallel = asin(sk_PLUS_1);
359  else
360  throw CoordinateConversionException( ErrorMessages::originLatitude );
361 
362  if (hemisphere == 'S')
363  standardParallel *= -1.0;
364 
365  if (centralMeridian > PI)
366  centralMeridian -= TWO_PI;
367  if (standardParallel < 0)
368  {
369  Southern_Hemisphere = 1;
370  Polar_Standard_Parallel = -standardParallel;
371  Polar_Central_Meridian = -centralMeridian;
372  }
373  else
374  {
375  Southern_Hemisphere = 0;
376  Polar_Standard_Parallel = standardParallel;
377  Polar_Central_Meridian = centralMeridian;
378  }
379 
380  sinolat = sin(Polar_Standard_Parallel);
381 
382  if (fabs(fabs(Polar_Standard_Parallel) - PI_OVER_2) > 1.0e-10)
383  {
384  essin = es * sinolat;
385  pow_es = polarPow(essin);
386  cosolat = cos(Polar_Standard_Parallel);
387  mc = cosolat / sqrt(1.0 - essin * essin);
388  Polar_a_mc = semiMajorAxis * mc;
389  Polar_tc = tan(PI_OVER_4 - Polar_Standard_Parallel / 2.0) / pow_es;
390  }
391 
392  /* Calculate Radius */
393  GeodeticCoordinates tempGeodeticCoordinates = GeodeticCoordinates( CoordinateType::geodetic, centralMeridian, 0 ) ;
394  MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &tempGeodeticCoordinates );
395  Polar_Delta_Northing = tempCoordinates->northing();
396  delete tempCoordinates;
397 
398  if(Polar_False_Northing)
399  Polar_Delta_Northing -= Polar_False_Northing;
400  if (Polar_Delta_Northing < 0)
401  Polar_Delta_Northing = -Polar_Delta_Northing;
402  Polar_Delta_Northing *= 1.01;
403 
404  Polar_Delta_Easting = Polar_Delta_Northing;
405 }
406 
407 
408 PolarStereographic::PolarStereographic( const PolarStereographic &ps )
409 {
410  coordinateType = ps.coordinateType;
411  semiMajorAxis = ps.semiMajorAxis;
412  flattening = ps.flattening;
413  es = ps.es;
414  es_OVER_2 = ps.es_OVER_2;
415  Southern_Hemisphere = ps.Southern_Hemisphere;
416  Polar_tc = ps.Polar_tc;
417  Polar_k90 = ps.Polar_k90;
418  Polar_a_mc = ps.Polar_a_mc;
419  two_Polar_a = ps.two_Polar_a;
420  Polar_Central_Meridian = ps.Polar_Central_Meridian;
421  Polar_Standard_Parallel = ps.Polar_Standard_Parallel;
422  Polar_False_Easting = ps.Polar_False_Easting;
423  Polar_False_Northing = ps.Polar_False_Northing;
424  Polar_Scale_Factor = ps.Polar_Scale_Factor;
425  Polar_Delta_Easting = ps.Polar_Delta_Easting;
426  Polar_Delta_Northing = ps.Polar_Delta_Northing;
427 }
428 
429 
430 PolarStereographic::~PolarStereographic()
431 {
432 }
433 
434 
435 PolarStereographic& PolarStereographic::operator=( const PolarStereographic &ps )
436 {
437  if( this != &ps )
438  {
439  coordinateType = ps.coordinateType;
440  semiMajorAxis = ps.semiMajorAxis;
441  flattening = ps.flattening;
442  es = ps.es;
443  es_OVER_2 = ps.es_OVER_2;
444  Southern_Hemisphere = ps.Southern_Hemisphere;
445  Polar_tc = ps.Polar_tc;
446  Polar_k90 = ps.Polar_k90;
447  Polar_a_mc = ps.Polar_a_mc;
448  two_Polar_a = ps.two_Polar_a;
449  Polar_Central_Meridian = ps.Polar_Central_Meridian;
450  Polar_Standard_Parallel = ps.Polar_Standard_Parallel;
451  Polar_False_Easting = ps.Polar_False_Easting;
452  Polar_False_Northing = ps.Polar_False_Northing;
453  Polar_Scale_Factor = ps.Polar_Scale_Factor;
454  Polar_Delta_Easting = ps.Polar_Delta_Easting;
455  Polar_Delta_Northing = ps.Polar_Delta_Northing;
456  }
457 
458  return *this;
459 }
460 
461 
462 PolarStereographicStandardParallelParameters* PolarStereographic::getStandardParallelParameters() const
463 {
464 /*
465  * The function getStandardParallelParameters returns the current
466  * ellipsoid parameters and Polar (Standard Parallel) projection parameters.
467  *
468  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
469  * ellipsoidFlattening : Flattening of ellipsoid (output)
470  * centralMeridian : Longitude down from pole, in radians (output)
471  * standardParallel : Latitude of true scale, in radians (output)
472  * falseEasting : Easting (X) at center of projection, in meters (output)
473  * falseNorthing : Northing (Y) at center of projection, in meters (output)
474  */
475 
476  return new PolarStereographicStandardParallelParameters( CoordinateType::polarStereographicStandardParallel, Polar_Central_Meridian, Polar_Standard_Parallel, Polar_False_Easting, Polar_False_Northing );
477 }
478 
479 
480 PolarStereographicScaleFactorParameters* PolarStereographic::getScaleFactorParameters() const
481 {
482 /*
483  * The function getScaleFactorParameters returns the current
484  * ellipsoid parameters and Polar (Scale Factor) projection parameters.
485  *
486  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
487  * ellipsoidFlattening : Flattening of ellipsoid (output)
488  * centralMeridian : Longitude down from pole, in radians (output)
489  * scaleFactor : Scale factor (output)
490  * falseEasting : Easting (X) at center of projection, in meters (output)
491  * falseNorthing : Northing (Y) at center of projection, in meters (output)
492  */
493 
494  if(Southern_Hemisphere == 0)
495  return new PolarStereographicScaleFactorParameters( CoordinateType::polarStereographicScaleFactor, Polar_Central_Meridian, Polar_Scale_Factor, 'N', Polar_False_Easting, Polar_False_Northing );
496  else
497  return new PolarStereographicScaleFactorParameters( CoordinateType::polarStereographicScaleFactor, Polar_Central_Meridian, Polar_Scale_Factor, 'S', Polar_False_Easting, Polar_False_Northing );
498 }
499 
500 
501 MSP::CCS::MapProjectionCoordinates* PolarStereographic::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
502 {
503 /*
504  * The function convertFromGeodetic converts geodetic
505  * coordinates (latitude and longitude) to Polar Stereographic coordinates
506  * (easting and northing), according to the current ellipsoid
507  * and Polar Stereographic projection parameters. If any errors occur,
508  * an exception is thrown with a description of the error.
509  *
510  * longitude : Longitude, in radians (input)
511  * latitude : Latitude, in radians (input)
512  * easting : Easting (X), in meters (output)
513  * northing : Northing (Y), in meters (output)
514  */
515 
516  double dlam;
517  double slat;
518  double essin;
519  double t;
520  double rho;
521  double pow_es;
522  double easting, northing;
523 
524  double longitude = geodeticCoordinates->longitude();
525  double latitude = geodeticCoordinates->latitude();
526 
527  if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
528  { /* latitude out of range */
529  throw CoordinateConversionException( ErrorMessages::latitude );
530  }
531  else if ((latitude < 0) && (Southern_Hemisphere == 0))
532  { /* latitude and Origin Latitude in different hemispheres */
533  throw CoordinateConversionException( ErrorMessages::latitude );
534  }
535  else if ((latitude > 0) && (Southern_Hemisphere == 1))
536  { /* latitude and Origin Latitude in different hemispheres */
537  throw CoordinateConversionException( ErrorMessages::latitude );
538  }
539  if ((longitude < -PI) || (longitude > TWO_PI))
540  { /* longitude out of range */
541  throw CoordinateConversionException( ErrorMessages::longitude );
542  }
543 
544  if (fabs(fabs(latitude) - PI_OVER_2) < 1.0e-10)
545  {
546  easting = Polar_False_Easting;
547  northing = Polar_False_Northing;
548  }
549  else
550  {
551  if (Southern_Hemisphere != 0)
552  {
553  longitude *= -1.0;
554  latitude *= -1.0;
555  }
556  dlam = longitude - Polar_Central_Meridian;
557  if (dlam > PI)
558  {
559  dlam -= TWO_PI;
560  }
561  if (dlam < -PI)
562  {
563  dlam += TWO_PI;
564  }
565  slat = sin(latitude);
566  essin = es * slat;
567  pow_es = polarPow(essin);
568  t = tan(PI_OVER_4 - latitude / 2.0) / pow_es;
569 
570  if (fabs(fabs(Polar_Standard_Parallel) - PI_OVER_2) > 1.0e-10)
571  rho = Polar_a_mc * t / Polar_tc;
572  else
573  rho = two_Polar_a * t / Polar_k90;
574 
575 
576  if (Southern_Hemisphere != 0)
577  {
578  easting = -(rho * sin(dlam) - Polar_False_Easting);
579  northing = rho * cos(dlam) + Polar_False_Northing;
580  }
581  else
582  {
583  easting = rho * sin(dlam) + Polar_False_Easting;
584  northing = -rho * cos(dlam) + Polar_False_Northing;
585  }
586  }
587 
588  return new MapProjectionCoordinates( coordinateType, easting, northing );
589 }
590 
591 
592 MSP::CCS::GeodeticCoordinates* PolarStereographic::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
593 {
594 /*
595  * The function convertToGeodetic converts Polar
596  * Stereographic coordinates (easting and northing) to geodetic
597  * coordinates (latitude and longitude) according to the current ellipsoid
598  * and Polar Stereographic projection Parameters. If any errors occur,
599  * an exception is thrown with a description of the error.
600  *
601  * easting : Easting (X), in meters (input)
602  * northing : Northing (Y), in meters (input)
603  * longitude : Longitude, in radians (output)
604  * latitude : Latitude, in radians (output)
605  *
606  */
607 
608  double dy = 0, dx = 0;
609  double rho = 0;
610  double t;
611  double PHI, sin_PHI;
612  double tempPHI = 0.0;
613  double essin;
614  double pow_es;
615  double delta_radius;
616  double longitude, latitude;
617 
618  double easting = mapProjectionCoordinates->easting();
619  double northing = mapProjectionCoordinates->northing();
620 
621  double min_easting = Polar_False_Easting - Polar_Delta_Easting;
622  double max_easting = Polar_False_Easting + Polar_Delta_Easting;
623  double min_northing = Polar_False_Northing - Polar_Delta_Northing;
624  double max_northing = Polar_False_Northing + Polar_Delta_Northing;
625 
626  if (easting > max_easting || easting < min_easting)
627  { /* easting out of range */
628  throw CoordinateConversionException( ErrorMessages::easting );
629  }
630  if (northing > max_northing || northing < min_northing)
631  { /* northing out of range */
632  throw CoordinateConversionException( ErrorMessages::northing );
633  }
634 
635  dy = northing - Polar_False_Northing;
636  dx = easting - Polar_False_Easting;
637 
638  /* Radius of point with origin of false easting, false northing */
639  rho = sqrt(dx * dx + dy * dy);
640 
641  delta_radius = sqrt(
642  Polar_Delta_Easting * Polar_Delta_Easting +
643  Polar_Delta_Northing * Polar_Delta_Northing);
644 
645  if(rho > delta_radius)
646  { /* Point is outside of projection area */
647  throw CoordinateConversionException( ErrorMessages::radius );
648  }
649 
650  if ((dy == 0.0) && (dx == 0.0))
651  {
652  latitude = PI_OVER_2;
653  longitude = Polar_Central_Meridian;
654  }
655  else
656  {
657  if (Southern_Hemisphere != 0)
658  {
659  dy *= -1.0;
660  dx *= -1.0;
661  }
662 
663  if (fabs(fabs(Polar_Standard_Parallel) - PI_OVER_2) > 1.0e-10)
664  t = rho * Polar_tc / (Polar_a_mc);
665  else
666  t = rho * Polar_k90 / (two_Polar_a);
667  PHI = PI_OVER_2 - 2.0 * atan(t);
668  while (fabs(PHI - tempPHI) > 1.0e-10)
669  {
670  tempPHI = PHI;
671  sin_PHI = sin(PHI);
672  essin = es * sin_PHI;
673  pow_es = polarPow(essin);
674  PHI = PI_OVER_2 - 2.0 * atan(t * pow_es);
675  }
676  latitude = PHI;
677  longitude = Polar_Central_Meridian + atan2(dx, -dy);
678 
679  if (longitude > PI)
680  longitude -= TWO_PI;
681  else if (longitude < -PI)
682  longitude += TWO_PI;
683 
684 
685  if (latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
686  latitude = PI_OVER_2;
687  else if (latitude < -PI_OVER_2)
688  latitude = -PI_OVER_2;
689 
690  if (longitude > PI) /* force distorted values to 180, -180 degrees */
691  longitude = PI;
692  else if (longitude < -PI)
693  longitude = -PI;
694 
695  }
696  if (Southern_Hemisphere != 0)
697  {
698  latitude *= -1.0;
699  longitude *= -1.0;
700  }
701 
702  return new GeodeticCoordinates(
703  CoordinateType::geodetic, longitude, latitude );
704 }
705 
706 
707 double PolarStereographic::polarPow( double esSin )
708 {
709  return pow((1.0 - esSin) / (1.0 + esSin), es_OVER_2);
710 }
711 
712 
713 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:02 for GeographicTranslator by doxygen 1.8.2