UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
Polyconic.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: POLYCONIC
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates
9  * (latitude and longitude in radians) and Polyconic projection coordinates
10  * (easting and northing in meters).
11  *
12  * ERROR HANDLING
13  *
14  * This component checks parameters for valid values. If an invalid value
15  * is found, the error code is combined with the current error code using
16  * the bitwise or. This combining allows multiple error codes to be
17  * returned. The possible error codes are:
18  *
19  * POLY_NO_ERROR : No errors occurred in function
20  * POLY_LAT_ERROR : Latitude outside of valid range
21  * (-90 to 90 degrees)
22  * POLY_LON_ERROR : Longitude outside of valid range
23  * (-180 to 360 degrees)
24  * POLY_EASTING_ERROR : Easting outside of valid range
25  * (False_Easting +/- ~20,500,000 m,
26  * depending on ellipsoid parameters
27  * and Origin_Latitude)
28  * POLY_NORTHING_ERROR : Northing outside of valid range
29  * (False_Northing +/- ~15,500,000 m,
30  * depending on ellipsoid parameters
31  * and Origin_Latitude)
32  * POLY_ORIGIN_LAT_ERROR : Origin latitude outside of valid range
33  * (-90 to 90 degrees)
34  * POLY_CENT_MER_ERROR : Central meridian outside of valid range
35  * (-180 to 360 degrees)
36  * POLY_A_ERROR : Semi-major axis less than or equal to zero
37  * POLY_INV_F_ERROR : Inverse flattening outside of valid range
38  * (250 to 350)
39  * POLY_LON_WARNING : Distortion will result if longitude is more
40  * than 90 degrees from the Central Meridian
41  *
42  * REUSE NOTES
43  *
44  * POLYCONIC is intended for reuse by any application that performs a
45  * Polyconic projection or its inverse.
46  *
47  * REFERENCES
48  *
49  * Further information on POLYCONIC can be found in the Reuse Manual.
50  *
51  * POLYCONIC originated from : U.S. Army Topographic Engineering Center
52  * Geospatial Information Division
53  * 7701 Telegraph Road
54  * Alexandria, VA 22310-3864
55  *
56  * LICENSES
57  *
58  * None apply to this component.
59  *
60  * RESTRICTIONS
61  *
62  * POLYCONIC has no restrictions.
63  *
64  * ENVIRONMENT
65  *
66  * POLYCONIC was tested and certified in the following environments:
67  *
68  * 1. Solaris 2.5 with GCC, version 2.8.1
69  * 2. Windows 95 with MS Visual C++, version 6
70  *
71  * MODIFICATIONS
72  *
73  * Date Description
74  * ---- -----------
75  * 10-06-99 Original Code
76  * 03-05-07 Original C++ Code
77  *
78  */
79 
80 
81 /***************************************************************************/
82 /*
83  * INCLUDES
84  */
85 
86 #include <math.h>
87 #include "Polyconic.h"
88 #include "MapProjection4Parameters.h"
89 #include "MapProjectionCoordinates.h"
90 #include "GeodeticCoordinates.h"
91 #include "CoordinateConversionException.h"
92 #include "ErrorMessages.h"
93 #include "WarningMessages.h"
94 
95 /*
96  * math.h - Standard C math library
97  * Polyconic.h - Is for prototype error checking
98  * MapProjectionCoordinates.h - defines map projection coordinates
99  * GeodeticCoordinates.h - defines geodetic coordinates
100  * CoordinateConversionException.h - Exception handler
101  * ErrorMessages.h - Contains exception messages
102  * WarningMessages.h - Contains warning messages
103  */
104 
105 
106 using namespace MSP::CCS;
107 
108 
109 /***************************************************************************/
110 /* DEFINES
111  *
112  */
113 
114 const double PI = 3.14159265358979323e0; /* PI */
115 const double PI_OVER_2 = ( PI / 2.0);
116 const double TWO_PI = (2.0 * PI);
117 const double FOURTY_ONE (41.0 * PI / 180); /* 41 degrees in radians */
118 
119 
120 double polyCoeffTimesSine( double coeff, double x, double latit )
121 {
122  return coeff * (sin (x * latit));
123 }
124 
125 
126 /************************************************************************/
127 /* FUNCTIONS
128  *
129  */
130 
131 Polyconic::Polyconic( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double originLatitude, double falseEasting, double falseNorthing ) :
132  CoordinateSystem(),
133  es2( 0.0066943799901413800 ),
134  es4( 4.4814723452405e-005 ),
135  es6( 3.0000678794350e-007 ),
136  M0( 0.0 ),
137  c0( .99832429845280 ),
138  c1( .0025146070605187 ),
139  c2( 2.6390465943377e-006 ),
140  c3( 3.4180460865959e-009 ),
141  Poly_Origin_Long( 0.0 ),
142  Poly_Origin_Lat( 0.0 ),
143  Poly_False_Easting( 0.0 ),
144  Poly_False_Northing( 0.0 ),
145  Poly_Max_Easting( 20037509.0 ),
146  Poly_Max_Northing( 15348215.0 ),
147  Poly_Min_Easting( -20037509.0 ),
148  Poly_Min_Northing( -15348215.0 )
149 {
150 /*
151  * The constructor receives the ellipsoid parameters and
152  * Polyconic projection parameters as inputs, and sets the corresponding state
153  * variables. If any errors occur, an exception is thrown with a description
154  * of the error.
155  *
156  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
157  * ellipsoidFlattening : Flattening of ellipsoid (input)
158  * centralMeridian : Longitude in radians at the center of (input)
159  * the projection
160  * originLatitude : Latitude in radians at which the (input)
161  * point scale factor is 1.0
162  * falseEasting : A coordinate value in meters assigned to the
163  * central meridian of the projection. (input)
164  * falseNorthing : A coordinate value in meters assigned to the
165  * origin latitude of the projection (input)
166  */
167 
168  double j, three_es4;
169  double lat, sin2lat, sin4lat, sin6lat;
170  double inv_f = 1 / ellipsoidFlattening;
171 
172  if (ellipsoidSemiMajorAxis <= 0.0)
173  { /* Semi-major axis must be greater than zero */
174  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
175  }
176  if ((inv_f < 250) || (inv_f > 350))
177  { /* Inverse flattening must be between 250 and 350 */
178  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
179  }
180  if ((originLatitude < -PI_OVER_2) || (originLatitude > PI_OVER_2))
181  { /* origin latitude out of range */
182  throw CoordinateConversionException( ErrorMessages::originLatitude );
183  }
184  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
185  { /* origin longitude out of range */
186  throw CoordinateConversionException( ErrorMessages::centralMeridian );
187  }
188 
189  semiMajorAxis = ellipsoidSemiMajorAxis;
190  flattening = ellipsoidFlattening;
191 
192  Poly_Origin_Lat = originLatitude;
193  if (centralMeridian > PI)
194  centralMeridian -= TWO_PI;
195  Poly_Origin_Long = centralMeridian;
196  Poly_False_Northing = falseNorthing;
197  Poly_False_Easting = falseEasting;
198  es2 = 2 * flattening - flattening * flattening;
199  es4 = es2 * es2;
200  es6 = es4 * es2;
201 
202  j = 45.0 * es6 / 1024.0;
203  three_es4 = 3.0 * es4;
204  c0 = 1.0 - es2 / 4.0 - three_es4 / 64.0 - 5.0 * es6 / 256.0;
205  c1 = 3.0 * es2 / 8.0 + three_es4 / 32.0 + j;
206  c2 = 15.0 * es4 / 256.0 + j;
207  c3 = 35.0 * es6 / 3072.0;
208 
209  lat = c0 * Poly_Origin_Lat;
210  sin2lat = polyCoeffTimesSine(c1, 2.0, Poly_Origin_Lat);
211  sin4lat = polyCoeffTimesSine(c2, 4.0, Poly_Origin_Lat);
212  sin6lat = polyCoeffTimesSine(c3, 6.0, Poly_Origin_Lat);
213  M0 = polyM(lat, sin2lat, sin4lat, sin6lat);
214 
215  if (Poly_Origin_Long > 0)
216  {
217  GeodeticCoordinates gcMax( CoordinateType::geodetic, Poly_Origin_Long - PI, FOURTY_ONE );
218  MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcMax );
219  Poly_Max_Northing = tempCoordinates->northing();
220  delete tempCoordinates;
221 
222  GeodeticCoordinates gcMin( CoordinateType::geodetic, Poly_Origin_Long - PI, -FOURTY_ONE );
223  tempCoordinates = convertFromGeodetic( &gcMin );
224  Poly_Min_Northing = tempCoordinates->northing();
225  delete tempCoordinates;
226 
227  Poly_Max_Easting = 19926189.0;
228  Poly_Min_Easting = -20037509.0;
229  }
230  else if (Poly_Origin_Long < 0)
231  {
232  GeodeticCoordinates gcMax( CoordinateType::geodetic, Poly_Origin_Long + PI, FOURTY_ONE );
233  MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcMax );
234  Poly_Max_Northing = tempCoordinates->northing();
235  delete tempCoordinates;
236 
237  GeodeticCoordinates gcMin( CoordinateType::geodetic, Poly_Origin_Long + PI, -FOURTY_ONE );
238  tempCoordinates = convertFromGeodetic( &gcMin );
239  Poly_Min_Northing = tempCoordinates->northing();
240  delete tempCoordinates;
241 
242  Poly_Max_Easting = 20037509.0;
243  Poly_Min_Easting = -19926189.0;
244  }
245  else
246  {
247  GeodeticCoordinates gcMax( CoordinateType::geodetic, PI, FOURTY_ONE );
248  MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcMax );
249  Poly_Max_Northing = tempCoordinates->northing();
250  delete tempCoordinates;
251 
252  GeodeticCoordinates gcMin( CoordinateType::geodetic, PI, -FOURTY_ONE );
253  tempCoordinates = convertFromGeodetic( &gcMin ); Poly_Min_Northing = tempCoordinates->northing();
254  delete tempCoordinates;
255 
256  Poly_Max_Easting = 20037509.0;
257  Poly_Min_Easting = -20037509.0;
258  }
259 
260  if(Poly_False_Northing)
261  {
262  Poly_Max_Northing -= Poly_False_Northing;
263  Poly_Min_Northing -= Poly_False_Northing;
264  }
265 }
266 
267 
268 Polyconic::Polyconic( const Polyconic &p )
269 {
270  semiMajorAxis = p.semiMajorAxis;
271  flattening = p.flattening;
272  es2 = p.es2;
273  es4 = p.es4;
274  es6 = p.es6;
275  M0 = p.M0;
276  c0 = p.c0;
277  c1 = p.c1;
278  c2 = p.c2;
279  c3 = p.c3;
280  Poly_Origin_Long = p.Poly_Origin_Long;
281  Poly_Origin_Lat = p.Poly_Origin_Lat;
282  Poly_False_Easting = p.Poly_False_Easting;
283  Poly_False_Northing = p.Poly_False_Northing;
284  Poly_Max_Easting = p.Poly_Max_Easting;
285  Poly_Max_Northing = p.Poly_Max_Northing;
286  Poly_Min_Easting = p.Poly_Min_Easting;
287  Poly_Min_Northing = p.Poly_Min_Northing;
288 }
289 
290 
291 Polyconic::~Polyconic()
292 {
293 }
294 
295 
296 Polyconic& Polyconic::operator=( const Polyconic &p )
297 {
298  if( this != &p )
299  {
300  semiMajorAxis = p.semiMajorAxis;
301  flattening = p.flattening;
302  es2 = p.es2;
303  es4 = p.es4;
304  es6 = p.es6;
305  M0 = p.M0;
306  c0 = p.c0;
307  c1 = p.c1;
308  c2 = p.c2;
309  c3 = p.c3;
310  Poly_Origin_Long = p.Poly_Origin_Long;
311  Poly_Origin_Lat = p.Poly_Origin_Lat;
312  Poly_False_Easting = p.Poly_False_Easting;
313  Poly_False_Northing = p.Poly_False_Northing;
314  Poly_Max_Easting = p.Poly_Max_Easting;
315  Poly_Max_Northing = p.Poly_Max_Northing;
316  Poly_Min_Easting = p.Poly_Min_Easting;
317  Poly_Min_Northing = p.Poly_Min_Northing;
318  }
319 
320  return *this;
321 }
322 
323 
324 MapProjection4Parameters* Polyconic::getParameters() const
325 {
326 /*
327  * The function getParameters returns the current ellipsoid
328  * parameters, and Polyconic projection parameters.
329  *
330  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
331  * ellipsoidFlattening : Flattening of ellipsoid (output)
332  * centralMeridian : Longitude in radians at the center of (output)
333  * the projection
334  * originLatitude : Latitude in radians at which the (output)
335  * point scale factor is 1.0
336  * falseEasting : A coordinate value in meters assigned to the
337  * central meridian of the projection. (output)
338  * falseNorthing : A coordinate value in meters assigned to the
339  * origin latitude of the projection (output)
340  */
341 
342  return new MapProjection4Parameters( CoordinateType::polyconic, Poly_Origin_Long, Poly_Origin_Lat, Poly_False_Easting, Poly_False_Northing );
343 }
344 
345 
346 MSP::CCS::MapProjectionCoordinates* Polyconic::convertFromGeodetic(
347  MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
348 {
349 /*
350  * The function convertFromGeodetic converts geodetic (latitude and
351  * longitude) coordinates to Polyconic projection (easting and northing)
352  * coordinates, according to the current ellipsoid and Polyconic projection
353  * parameters. If any errors occur, an exception is thrown with a description
354  * of the error.
355  *
356  * longitude : Longitude (lambda) in radians (input)
357  * latitude : Latitude (phi) in radians (input)
358  * easting : Easting (X) in meters (output)
359  * northing : Northing (Y) in meters (output)
360  */
361 
362  double lat, sin2lat, sin4lat, sin6lat;
363  double dlam; /* Longitude - Central Meridan */
364  double NN;
365  double NN_OVER_tlat;
366  double MM;
367  double EE;
368  double easting, northing;
369 
370  double longitude = geodeticCoordinates->longitude();
371  double latitude = geodeticCoordinates->latitude();
372  double slat = sin(latitude);
373 
374  if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
375  { /* Latitude out of range */
376  throw CoordinateConversionException( ErrorMessages::latitude );
377  }
378  if ((longitude < -PI) || (longitude > TWO_PI))
379  { /* Longitude out of range */
380  throw CoordinateConversionException( ErrorMessages::longitude );
381  }
382 
383  char warning[256];
384  warning[0] = '\0';
385  dlam = longitude - Poly_Origin_Long;
386  if (fabs(dlam) > (PI / 2))
387  { /* Distortion results if Longitude is > 90 deg from the Central Meridian */
388  strcat( warning, MSP::CCS::WarningMessages::longitude );
389  }
390  if (dlam > PI)
391  {
392  dlam -= TWO_PI;
393  }
394  if (dlam < -PI)
395  {
396  dlam += TWO_PI;
397  }
398  if (latitude == 0.0)
399  {
400  easting = semiMajorAxis * dlam + Poly_False_Easting;
401  northing = -M0 + Poly_False_Northing;
402  }
403  else
404  {
405  NN = semiMajorAxis / sqrt(1.0 - es2 * (slat * slat));
406  NN_OVER_tlat = NN / tan(latitude);
407  lat = c0 * latitude;
408  sin2lat = polyCoeffTimesSine(c1, 2.0, latitude);
409  sin4lat = polyCoeffTimesSine(c2, 4.0, latitude);
410  sin6lat = polyCoeffTimesSine(c3, 6.0, latitude);
411  MM = polyM(lat, sin2lat, sin4lat, sin6lat);
412  EE = dlam * slat;
413  easting = NN_OVER_tlat * sin(EE) + Poly_False_Easting;
414  northing = MM - M0 + NN_OVER_tlat * (1.0 - cos(EE)) +
415  Poly_False_Northing;
416  }
417 
418  return new MapProjectionCoordinates(
419  CoordinateType::polyconic, warning, easting, northing );
420 }
421 
422 
423 MSP::CCS::GeodeticCoordinates* Polyconic::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
424 {
425 /*
426  * The function convertToGeodetic converts Polyconic projection
427  * (easting and northing) coordinates to geodetic (latitude and longitude)
428  * coordinates, according to the current ellipsoid and Polyconic projection
429  * coordinates. If any errors occur, an exception is thrown with a description
430  * of the error.
431  *
432  * easting : Easting (X) in meters (input)
433  * northing : Northing (Y) in meters (input)
434  * longitude : Longitude (lambda) in radians (output)
435  * latitude : Latitude (phi) in radians (output)
436  */
437 
438  double dx; /* Delta easting - Difference in easting (easting-FE) */
439  double dy; /* Delta northing - Difference in northing (northing-FN) */
440  double dx_OVER_Poly_a;
441  double AA;
442  double BB;
443  double CC = 0.0;
444  double PHIn, Delta_PHI = 1.0;
445  double sin_PHIn;
446  double PHI, sin2PHI,sin4PHI, sin6PHI;
447  double Mn, Mn_prime, Ma;
448  double AA_Ma;
449  double Ma2_PLUS_BB;
450  double AA_MINUS_Ma;
451  double tolerance = 1.0e-12; /* approximately 1/1000th of
452  an arc second or 1/10th meter */
453  double longitude, latitude;
454  int count = 45000;
455 
456  double easting = mapProjectionCoordinates->easting();
457  double northing = mapProjectionCoordinates->northing();
458 
459  if ((easting < (Poly_False_Easting + Poly_Min_Easting))
460  || (easting > (Poly_False_Easting + Poly_Max_Easting)))
461  { /* Easting out of range */
462  throw CoordinateConversionException( ErrorMessages::easting );
463  }
464  if ((northing < (Poly_False_Northing + Poly_Min_Northing))
465  || (northing > (Poly_False_Northing + Poly_Max_Northing)))
466  { /* Northing out of range */
467  throw CoordinateConversionException( ErrorMessages::northing );
468  }
469 
470  dy = northing - Poly_False_Northing;
471  dx = easting - Poly_False_Easting;
472  dx_OVER_Poly_a = dx / semiMajorAxis;
473  if (floatEq(dy,-M0,1))
474  {
475  latitude = 0.0;
476  longitude = dx_OVER_Poly_a + Poly_Origin_Long;
477  }
478  else
479  {
480  AA = (M0 + dy) / semiMajorAxis;
481  BB = dx_OVER_Poly_a * dx_OVER_Poly_a + (AA * AA);
482  PHIn = AA;
483 
484  while (fabs(Delta_PHI) > tolerance && count)
485  {
486  sin_PHIn = sin(PHIn);
487  CC = sqrt(1.0 - es2 * sin_PHIn * sin_PHIn) * tan(PHIn);
488  PHI = c0 * PHIn;
489  sin2PHI = polyCoeffTimesSine(c1, 2.0, PHIn);
490  sin4PHI = polyCoeffTimesSine(c2, 4.0, PHIn);
491  sin6PHI = polyCoeffTimesSine(c3, 6.0, PHIn);
492  Mn = polyM(PHI, sin2PHI, sin4PHI, sin6PHI);
493  Mn_prime = c0 - 2.0 * c1 * cos(2.0 * PHIn) + 4.0 * c2 * cos(4.0 * PHIn) -
494  6.0 * c3 * cos(6.0 * PHIn);
495  Ma = Mn / semiMajorAxis;
496  AA_Ma = AA * Ma;
497  Ma2_PLUS_BB = Ma * Ma + BB;
498  AA_MINUS_Ma = AA - Ma;
499  Delta_PHI = (AA_Ma * CC + AA_MINUS_Ma - 0.5 * (Ma2_PLUS_BB) * CC) /
500  (es2 * sin2PHI * (Ma2_PLUS_BB - 2.0 * AA_Ma) /
501  4.0 * CC + (AA_MINUS_Ma) * (CC * Mn_prime - 2.0 / sin2PHI) - Mn_prime);
502  PHIn -= Delta_PHI;
503  count --;
504  }
505 
506  if(!count)
507  throw CoordinateConversionException( ErrorMessages::northing );
508 
509  latitude = PHIn;
510 
511  if (latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
512  latitude = PI_OVER_2;
513  else if (latitude < -PI_OVER_2)
514  latitude = -PI_OVER_2;
515 
516  if (floatEq(fabs(latitude),PI_OVER_2,.00001) || (latitude == 0))
517  longitude = Poly_Origin_Long;
518 
519  else
520  {
521  longitude = (asin(dx_OVER_Poly_a * CC)) / sin(latitude) +
522  Poly_Origin_Long;
523  }
524  }
525  if (longitude > PI)
526  longitude -= TWO_PI;
527  if (longitude < -PI)
528  longitude += TWO_PI;
529 
530  if (longitude > PI) /* force distorted values to 180, -180 degrees */
531  longitude = PI;
532  else if (longitude < -PI)
533  longitude = -PI;
534 
535  return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
536 }
537 
538 
539 double Polyconic::polyM( double c0lat, double c1s2lat, double c2s4lat, double c3s6lat )
540 {
541  return semiMajorAxis * (c0lat - c1s2lat + c2s4lat - c3s6lat);
542 }
543 
544 
545 double Polyconic::floatEq( double x, double v, double epsilon )
546 {
547  return ((v - epsilon) < x) && (x < (v + epsilon));
548 }
549 
550 
551 
552 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:02 for GeographicTranslator by doxygen 1.8.2