UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
AlbersEqualAreaConic.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: ALBERS
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates
9  * (latitude and longitude in radians) and Albers Equal Area Conic
10  * projection coordinates (easting and northing in meters) defined
11  * by two standard parallels.
12  *
13  * ERROR HANDLING
14  *
15  * This component checks parameters for valid values. If an invalid value
16  * is found the error code is combined with the current error code using
17  * the bitwise or. This combining allows multiple error codes to be
18  * returned. The possible error codes are:
19  *
20  * ALBERS_NO_ERROR : No errors occurred in function
21  * ALBERS_LAT_ERROR : Latitude outside of valid range
22  * (-90 to 90 degrees)
23  * ALBERS_LON_ERROR : Longitude outside of valid range
24  * (-180 to 360 degrees)
25  * ALBERS_EASTING_ERROR : Easting outside of valid range
26  * (depends on ellipsoid and projection
27  * parameters)
28  * ALBERS_NORTHING_ERROR : Northing outside of valid range
29  * (depends on ellipsoid and projection
30  * parameters)
31  * ALBERS_FIRST_STDP_ERROR : First standard parallel outside of valid
32  * range (-90 to 90 degrees)
33  * ALBERS_SECOND_STDP_ERROR : Second standard parallel outside of valid
34  * range (-90 to 90 degrees)
35  * ALBERS_ORIGIN_LAT_ERROR : Origin latitude outside of valid range
36  * (-90 to 90 degrees)
37  * ALBERS_CENT_MER_ERROR : Central meridian outside of valid range
38  * (-180 to 360 degrees)
39  * ALBERS_A_ERROR : Semi-major axis less than or equal to zero
40  * ALBERS_INV_F_ERROR : Inverse flattening outside of valid range
41  * (250 to 350)
42  * ALBERS_HEMISPHERE_ERROR : Standard parallels cannot be opposite
43  * latitudes
44  * ALBERS_FIRST_SECOND_ERROR : The 1st & 2nd standard parallels cannot
45  * both be 0
46  *
47  *
48  * REUSE NOTES
49  *
50  * ALBERS is intended for reuse by any application that performs an Albers
51  * Equal Area Conic projection or its inverse.
52  *
53  * REFERENCES
54  *
55  * Further information on ALBERS can be found in the Reuse Manual.
56  *
57  * ALBERS originated from: U.S. Army Topographic Engineering Center
58  * Geospatial Information Division
59  * 7701 Telegraph Road
60  * Alexandria, VA 22310-3864
61  *
62  * LICENSES
63  *
64  * None apply to this component.
65  *
66  * RESTRICTIONS
67  *
68  * ALBERS has no restrictions.
69  *
70  * ENVIRONMENT
71  *
72  * ALBERS was tested and certified in the following environments:
73  *
74  * 1. Solaris 2.5 with GCC, version 2.8.1
75  * 2. MSDOS with MS Visual C++, version 6
76  *
77  * MODIFICATIONS
78  *
79  * Date Description
80  * ---- -----------
81  * 07-09-99 Original Code
82  * 03-08-07 Original C++ Code
83  *
84  *
85  */
86 
87 
88 /***************************************************************************/
89 /*
90  * INCLUDES
91  */
92 
93 #include <math.h>
94 #include "AlbersEqualAreaConic.h"
95 #include "MapProjection6Parameters.h"
96 #include "MapProjectionCoordinates.h"
97 #include "GeodeticCoordinates.h"
98 #include "CoordinateConversionException.h"
99 #include "ErrorMessages.h"
100 
101 /*
102  * math.h - Standard C math library
103  * AlbersEqualAreaConic.h - Is for prototype error checking
104  * MapProjectionCoordinates.h - defines map projection coordinates
105  * GeodeticCoordinates.h - defines geodetic coordinates
106  * CoordinateConversionException.h - Exception handler
107  * ErrorMessages.h - Contains exception messages
108  */
109 
110 
111 using namespace MSP::CCS;
112 
113 
114 /***************************************************************************/
115 /* DEFINES
116  *
117  */
118 
119 const double PI = 3.14159265358979323e0; /* PI */
120 const double PI_OVER_2 = ( PI / 2.0);
121 const double TWO_PI = (2.0 * PI);
122 
123 
124 double oneMinusSqr( double x )
125 {
126  return 1.0 - x * x;
127 }
128 
129 
130 double albersM( double clat, double oneminussqressin )
131 {
132  return clat / sqrt(oneminussqressin);
133 }
134 
135 
136 /************************************************************************/
137 /* FUNCTIONS
138  *
139  */
140 
141 AlbersEqualAreaConic::AlbersEqualAreaConic(
142  double ellipsoidSemiMajorAxis,
143  double ellipsoidFlattening,
144  double centralMeridian,
145  double originLatitude,
146  double standardParallel1,
147  double standardParallel2,
148  double falseEasting,
149  double falseNorthing ) :
150  CoordinateSystem(),
151  es( 0.08181919084262188000 ),
152  es2( 0.0066943799901413800 ),
153  C( 1.4896626908850 ),
154  rho0( 6388749.3391064 ),
155  n( .70443998701755 ),
156  Albers_a_OVER_n( 9054194.9882824 ),
157  one_MINUS_es2( .99330562000986 ),
158  two_es( .16363838168524 ),
159  Albers_Origin_Long( 0.0 ),
160  Albers_Origin_Lat( (45 * PI / 180) ),
161  Albers_Std_Parallel_1( 40 * PI / 180 ),
162  Albers_Std_Parallel_2( 50 * PI / 180 ),
163  Albers_False_Easting( 0.0 ),
164  Albers_False_Northing( 0.0 ),
165  Albers_Delta_Northing( 40000000 ),
166  Albers_Delta_Easting( 40000000 )
167 {
168 /*
169  * The constructor receives the ellipsoid parameters and
170  * projection parameters as inputs, and sets the corresponding state
171  * variables. If any errors occur, an exception is thrown with a description
172  * of the error.
173  *
174  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
175  * ellipsoidFlattening : Flattening of ellipsoid (input)
176  * centralMeridian : Longitude in radians at the center of (input)
177  * the projection
178  * originLatitude : Latitude in radians at which the (input)
179  * point scale factor is 1.0
180  * standardParallel1 : First standard parallel (input)
181  * standardParallel2 : Second standard parallel (input)
182  * falseEasting : A coordinate value in meters assigned to the
183  * central meridian of the projection. (input)
184  * falseNorthing : A coordinate value in meters assigned to the
185  * origin latitude of the projection (input)
186  */
187 
188  double sin_lat, sin_lat_1, cos_lat;
189  double m1, m2, SQRm1;
190  double q0, q1, q2;
191  double es_sin, one_MINUS_SQRes_sin;
192  double nq0;
193  double inv_f = 1 / ellipsoidFlattening;
194 
195  if (ellipsoidSemiMajorAxis <= 0.0)
196  { /* Semi-major axis must be greater than zero */
197  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
198  }
199  if ((inv_f < 250) || (inv_f > 350))
200  { /* Inverse flattening must be between 250 and 350 */
201  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
202  }
203  if ((originLatitude < -PI_OVER_2) || (originLatitude > PI_OVER_2))
204  { /* origin latitude out of range */
205  throw CoordinateConversionException( ErrorMessages::originLatitude );
206  }
207  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
208  { /* origin longitude out of range */
209  throw CoordinateConversionException( ErrorMessages::centralMeridian );
210  }
211  if ((standardParallel1 < -PI_OVER_2) || (standardParallel1 > PI_OVER_2))
212  { /* First Standard Parallel out of range */
213  throw CoordinateConversionException( ErrorMessages::standardParallel1 );
214  }
215  if ((standardParallel2 < -PI_OVER_2) || (standardParallel2 > PI_OVER_2))
216  { /* Second Standard Parallel out of range */
217  throw CoordinateConversionException( ErrorMessages::standardParallel2 );
218  }
219  if ((standardParallel1 == 0.0) && (standardParallel2 == 0.0))
220  { /* First & Second Standard Parallels equal 0 */
221  throw CoordinateConversionException( ErrorMessages::standardParallel1_2 );
222  }
223  if (standardParallel1 == -standardParallel2)
224  { /* Parallels are opposite latitudes */
225  throw CoordinateConversionException( ErrorMessages::standardParallelHemisphere );
226  }
227 
228 
229  semiMajorAxis = ellipsoidSemiMajorAxis;
230  flattening = ellipsoidFlattening;
231 
232  Albers_Origin_Lat = originLatitude;
233  Albers_Std_Parallel_1 = standardParallel1;
234  Albers_Std_Parallel_2 = standardParallel2;
235  if (centralMeridian > PI)
236  centralMeridian -= TWO_PI;
237  Albers_Origin_Long = centralMeridian;
238  Albers_False_Easting = falseEasting;
239  Albers_False_Northing = falseNorthing;
240 
241  es2 = 2 * flattening - flattening * flattening;
242  es = sqrt(es2);
243  one_MINUS_es2 = 1 - es2;
244  two_es = 2 * es;
245 
246  sin_lat = sin(Albers_Origin_Lat);
247  es_sin = esSine(sin_lat);
248  one_MINUS_SQRes_sin = oneMinusSqr( es_sin );
249  q0 = albersQ( sin_lat, one_MINUS_SQRes_sin, es_sin );
250 
251  sin_lat_1 = sin(Albers_Std_Parallel_1);
252  cos_lat = cos(Albers_Std_Parallel_1);
253  es_sin = esSine(sin_lat_1);
254  one_MINUS_SQRes_sin = oneMinusSqr( es_sin );
255  m1 = albersM( cos_lat, one_MINUS_SQRes_sin );
256  q1 = albersQ( sin_lat_1, one_MINUS_SQRes_sin, es_sin );
257 
258  SQRm1 = m1 * m1;
259  if (fabs(Albers_Std_Parallel_1 - Albers_Std_Parallel_2) > 1.0e-10)
260  {
261  sin_lat = sin(Albers_Std_Parallel_2);
262  cos_lat = cos(Albers_Std_Parallel_2);
263  es_sin = esSine(sin_lat);
264  one_MINUS_SQRes_sin = oneMinusSqr( es_sin );
265  m2 = albersM( cos_lat, one_MINUS_SQRes_sin );
266  q2 = albersQ( sin_lat, one_MINUS_SQRes_sin, es_sin );
267  n = (SQRm1 - m2 * m2) / (q2 - q1);
268  }
269  else
270  n = sin_lat_1;
271 
272  C = SQRm1 + n * q1;
273  Albers_a_OVER_n = semiMajorAxis / n;
274  nq0 = n * q0;
275  if (C < nq0)
276  rho0 = 0;
277  else
278  rho0 = Albers_a_OVER_n * sqrt(C - nq0);
279 }
280 
281 
282 AlbersEqualAreaConic::AlbersEqualAreaConic( const AlbersEqualAreaConic &aeac )
283 {
284  semiMajorAxis = aeac.semiMajorAxis;
285  flattening = aeac.flattening;
286  es = aeac.es;
287  es2 = aeac.es2;
288  C = aeac.C;
289  rho0 = aeac.rho0;
290  n = aeac.n;
291  Albers_a_OVER_n = aeac.Albers_a_OVER_n;
292  one_MINUS_es2 = aeac.one_MINUS_es2;
293  two_es = aeac.two_es;
294  Albers_Origin_Long = aeac.Albers_Origin_Long;
295  Albers_Origin_Lat = aeac.Albers_Origin_Lat;
296  Albers_Std_Parallel_1 = aeac.Albers_Std_Parallel_1;
297  Albers_Std_Parallel_2 = aeac.Albers_Std_Parallel_2;
298  Albers_False_Easting = aeac.Albers_False_Easting;
299  Albers_False_Northing = aeac.Albers_False_Northing;
300  Albers_Delta_Northing = aeac.Albers_Delta_Northing;
301  Albers_Delta_Easting = aeac.Albers_Delta_Easting;
302 }
303 
304 
305 AlbersEqualAreaConic::~AlbersEqualAreaConic()
306 {
307 }
308 
309 
310 AlbersEqualAreaConic& AlbersEqualAreaConic::operator=(
311  const AlbersEqualAreaConic &aeac )
312 {
313  if( this != &aeac )
314  {
315  semiMajorAxis = aeac.semiMajorAxis;
316  flattening = aeac.flattening;
317  es = aeac.es;
318  es2 = aeac.es2;
319  C = aeac.C;
320  rho0 = aeac.rho0;
321  n = aeac.n;
322  Albers_a_OVER_n = aeac.Albers_a_OVER_n;
323  one_MINUS_es2 = aeac.one_MINUS_es2;
324  two_es = aeac.two_es;
325  Albers_Origin_Long = aeac.Albers_Origin_Long;
326  Albers_Origin_Lat = aeac.Albers_Origin_Lat;
327  Albers_Std_Parallel_1 = aeac.Albers_Std_Parallel_1;
328  Albers_Std_Parallel_2 = aeac.Albers_Std_Parallel_2;
329  Albers_False_Easting = aeac.Albers_False_Easting;
330  Albers_False_Northing = aeac.Albers_False_Northing;
331  Albers_Delta_Northing = aeac.Albers_Delta_Northing;
332  Albers_Delta_Easting = aeac.Albers_Delta_Easting;
333  }
334 
335  return *this;
336 }
337 
338 
339 MapProjection6Parameters* AlbersEqualAreaConic::getParameters() const
340 {
341 /*
342  * The function getParameters returns the current ellipsoid
343  * parameters, and Albers projection parameters.
344  *
345  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
346  * ellipsoidFlattening : Flattening of ellipsoid (output)
347  * centralMeridian : Longitude in radians at the center of (output)
348  * the projection
349  * originLatitude : Latitude in radians at which the (output)
350  * point scale factor is 1.0
351  * standardParallel1 : First standard parallel (output)
352  * standardParallel2 : Second standard parallel (output)
353  * falseEasting : A coordinate value in meters assigned to the
354  * central meridian of the projection. (output)
355  * falseNorthing : A coordinate value in meters assigned to the
356  * origin latitude of the projection (output)
357  */
358 
359  return new MapProjection6Parameters(
360  CoordinateType::albersEqualAreaConic,
361  Albers_Origin_Long,
362  Albers_Origin_Lat,
363  Albers_Std_Parallel_1,
364  Albers_Std_Parallel_2,
365  Albers_False_Easting,
366  Albers_False_Northing );
367 }
368 
369 
370 MSP::CCS::MapProjectionCoordinates* AlbersEqualAreaConic::convertFromGeodetic(
371  MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
372 {
373 /*
374  * The function convertFromGeodetic converts geodetic (latitude and
375  * longitude) coordinates to Albers projection (easting and northing)
376  * coordinates, according to the current ellipsoid and Albers projection
377  * parameters. If any errors occur, an exception is thrown with a description
378  * of the error.
379  *
380  * longitude : Longitude (lambda) in radians (input)
381  * latitude : Latitude (phi) in radians (input)
382  * easting : Easting (X) in meters (output)
383  * northing : Northing (Y) in meters (output)
384  */
385 
386  double dlam; /* Longitude - Central Meridan */
387  double sin_lat, cos_lat;
388  double es_sin, one_MINUS_SQRes_sin;
389  double q;
390  double rho;
391  double theta;
392  double nq;
393 
394  double longitude = geodeticCoordinates->longitude();
395  double latitude = geodeticCoordinates->latitude();
396 
397  if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
398  { /* Latitude out of range */
399  throw CoordinateConversionException( ErrorMessages::latitude );
400  }
401  if ((longitude < -PI) || (longitude > TWO_PI))
402  { /* Longitude out of range */
403  throw CoordinateConversionException( ErrorMessages::longitude );
404  }
405 
406  dlam = longitude - Albers_Origin_Long;
407  if (dlam > PI)
408  {
409  dlam -= TWO_PI;
410  }
411  if (dlam < -PI)
412  {
413  dlam += TWO_PI;
414  }
415  sin_lat = sin(latitude);
416  cos_lat = cos(latitude);
417  es_sin = esSine( sin_lat );
418  one_MINUS_SQRes_sin = oneMinusSqr( es_sin );
419  q = albersQ( sin_lat, one_MINUS_SQRes_sin, es_sin );
420  nq = n * q;
421  if (C < nq)
422  rho = 0;
423  else
424  rho = Albers_a_OVER_n * sqrt(C - nq);
425 
426 
427  theta = n * dlam;
428  double easting = rho * sin(theta) + Albers_False_Easting;
429  double northing = rho0 - rho * cos(theta) + Albers_False_Northing;
430 
431  return new MapProjectionCoordinates(
432  CoordinateType::albersEqualAreaConic, easting, northing );
433 }
434 
435 
436 MSP::CCS::GeodeticCoordinates* AlbersEqualAreaConic::convertToGeodetic(
437  MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
438 {
439 /*
440  * The function convertToGeodetic converts Albers projection
441  * (easting and northing) coordinates to geodetic (latitude and longitude)
442  * coordinates, according to the current ellipsoid and Albers projection
443  * coordinates. If any errors occur, an exception is thrown with a description
444  * of the error.
445  *
446  * easting : Easting (X) in meters (input)
447  * northing : Northing (Y) in meters (input)
448  * latitude : Latitude (phi) in radians (output)
449  * longitude : Longitude (lambda) in radians (output)
450  */
451 
452  double dy, dx;
453  double rho0_MINUS_dy;
454  double q, qconst, q_OVER_2;
455  double rho, rho_n;
456  double PHI, Delta_PHI = 1.0;
457  double sin_phi;
458  double es_sin, one_MINUS_SQRes_sin;
459  double theta = 0.0;
460  int count = 60;
461  double longitude, latitude;
462  double tolerance = 4.85e-10; /* approximately 1/1000th of
463  an arc second or 1/10th meter */
464 
465  double easting = mapProjectionCoordinates->easting();
466  double northing = mapProjectionCoordinates->northing();
467 
468  if( (easting < (Albers_False_Easting - Albers_Delta_Easting))
469  || (easting > Albers_False_Easting + Albers_Delta_Easting))
470  { /* Easting out of range */
471  throw CoordinateConversionException( ErrorMessages::easting );
472  }
473  if( (northing < (Albers_False_Northing - Albers_Delta_Northing))
474  || (northing > Albers_False_Northing + Albers_Delta_Northing))
475  { /* Northing out of range */
476  throw CoordinateConversionException( ErrorMessages::northing );
477  }
478 
479  dy = northing - Albers_False_Northing;
480  dx = easting - Albers_False_Easting;
481  rho0_MINUS_dy = rho0 - dy;
482  rho = sqrt(dx * dx + rho0_MINUS_dy * rho0_MINUS_dy);
483 
484  if (n < 0)
485  {
486  rho *= -1.0;
487  dy *= -1.0;
488  dx *= -1.0;
489  rho0_MINUS_dy *= -1.0;
490  }
491 
492  if (rho != 0.0)
493  theta = atan2(dx, rho0_MINUS_dy);
494  rho_n = rho * n;
495  q = (C - (rho_n * rho_n) / (semiMajorAxis * semiMajorAxis)) / n;
496  qconst = 1 - ((one_MINUS_es2) / (two_es)) * log((1.0 - es) / (1.0 + es));
497  if (fabs(fabs(qconst) - fabs(q)) > 1.0e-6)
498  {
499  q_OVER_2 = q / 2.0;
500  if (q_OVER_2 > 1.0)
501  latitude = PI_OVER_2;
502  else if (q_OVER_2 < -1.0)
503  latitude = -PI_OVER_2;
504  else
505  {
506  PHI = asin(q_OVER_2);
507  if (es < 1.0e-10)
508  latitude = PHI;
509  else
510  {
511  while ((fabs(Delta_PHI) > tolerance) && count)
512  {
513  sin_phi = sin(PHI);
514  es_sin = esSine( sin_phi );
515  one_MINUS_SQRes_sin = oneMinusSqr( es_sin );
516  Delta_PHI =
517  (one_MINUS_SQRes_sin * one_MINUS_SQRes_sin) / (2.0 * cos(PHI)) *
518  (q / (one_MINUS_es2) - sin_phi / one_MINUS_SQRes_sin +
519  (log((1.0 - es_sin) / (1.0 + es_sin)) / (two_es)));
520  PHI += Delta_PHI;
521  count --;
522  }
523 
524  if(!count)
525  throw CoordinateConversionException( ErrorMessages::northing );
526 
527  latitude = PHI;
528  }
529 
530  if (latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
531  latitude = PI_OVER_2;
532  else if (latitude < -PI_OVER_2)
533  latitude = -PI_OVER_2;
534 
535  }
536  }
537  else
538  {
539  if (q >= 0.0)
540  latitude = PI_OVER_2;
541  else
542  latitude = -PI_OVER_2;
543  }
544  longitude = Albers_Origin_Long + theta / n;
545 
546  if (longitude > PI)
547  longitude -= TWO_PI;
548  if (longitude < -PI)
549  longitude += TWO_PI;
550 
551  if (longitude > PI) /* force distorted values to 180, -180 degrees */
552  longitude = PI;
553  else if (longitude < -PI)
554  longitude = -PI;
555 
556  return new GeodeticCoordinates(
557  CoordinateType::geodetic, longitude, latitude );
558 }
559 
560 
561 double AlbersEqualAreaConic::esSine( double sinlat )
562 {
563  return es * sinlat;
564 }
565 
566 
567 double AlbersEqualAreaConic::albersQ(
568  double slat, double oneminussqressin, double essin )
569 {
570  return (one_MINUS_es2)*(slat / (oneminussqressin) -
571  (1 / (two_es)) *log((1 - essin) / (1 + essin)));
572 }
573 
574 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:01 for GeographicTranslator by doxygen 1.8.2