UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
CylindricalEqualArea.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: CYLINDRICAL EQUAL AREA
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates
9  * (latitude and longitude in radians) and Cylindrical Equal Area projection
10  * coordinates (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  * CYEQ_NO_ERROR : No errors occurred in function
20  * CYEQ_LAT_ERROR : Latitude outside of valid range
21  * (-90 to 90 degrees)
22  * CYEQ_LON_ERROR : Longitude outside of valid range
23  * (-180 to 360 degrees)
24  * CYEQ_EASTING_ERROR : Easting outside of valid range
25  * (False_Easting +/- ~20,000,000 m,
26  * depending on ellipsoid parameters
27  * and Origin_Latitude)
28  * CYEQ_NORTHING_ERROR : Northing outside of valid range
29  * (False_Northing +/- ~6,000,000 m,
30  * depending on ellipsoid parameters
31  * and Origin_Latitude)
32  * CYEQ_ORIGIN_LAT_ERROR : Origin latitude outside of valid range
33  * (-90 to 90 degrees)
34  * CYEQ_CENT_MER_ERROR : Central meridian outside of valid range
35  * (-180 to 360 degrees)
36  * CYEQ_A_ERROR : Semi-major axis less than or equal to zero
37  * CYEQ_INV_F_ERROR : Inverse flattening outside of valid range
38  * (250 to 350)
39  *
40  * REUSE NOTES
41  *
42  * CYLINDRICAL EQUAL AREA is intended for reuse by any application that performs a
43  * Cylindrical Equal Area projection or its inverse.
44  *
45  * REFERENCES
46  *
47  * Further information on CYLINDRICAL EQUAL AREA can be found in the Reuse Manual.
48  *
49  * CYLINDRICAL EQUAL AREA originated from :
50  * U.S. Army Topographic Engineering Center
51  * Geospatial Information Division
52  * 7701 Telegraph Road
53  * Alexandria, VA 22310-3864
54  *
55  * LICENSES
56  *
57  * None apply to this component.
58  *
59  * RESTRICTIONS
60  *
61  * CYLINDRICAL EQUAL AREA has no restrictions.
62  *
63  * ENVIRONMENT
64  *
65  * CYLINDRICAL EQUAL AREA was tested and certified in the following environments:
66  *
67  * 1. Solaris 2.5 with GCC 2.8.1
68  * 2. MS Windows 95 with MS Visual C++ 6
69  *
70  * MODIFICATIONS
71  *
72  * Date Description
73  * ---- -----------
74  * 04-16-99 Original Code
75  * 03-07-07 Original C++ Code
76  *
77  */
78 
79 
80 /***************************************************************************/
81 /*
82  * INCLUDES
83  */
84 
85 #include <math.h>
86 #include "CylindricalEqualArea.h"
87 #include "MapProjection4Parameters.h"
88 #include "MapProjectionCoordinates.h"
89 #include "GeodeticCoordinates.h"
90 #include "CoordinateConversionException.h"
91 #include "ErrorMessages.h"
92 
93 /*
94  * math.h - Standard C math library
95  * CylindricalEqualArea.h - Is for prototype error checking
96  * MapProjectionCoordinates.h - defines map projection coordinates
97  * GeodeticCoordinates.h - defines geodetic coordinates
98  * CoordinateConversionException.h - Exception handler
99  * ErrorMessages.h - Contains exception messages
100  */
101 
102 
103 using namespace MSP::CCS;
104 
105 
106 /***************************************************************************/
107 /* DEFINES
108  *
109  */
110 
111 const double PI = 3.14159265358979323e0; /* PI */
112 const double PI_OVER_2 = ( PI / 2.0);
113 const double TWO_PI = (2.0 * PI);
114 const double ONE = (1.0 * PI / 180); /* 1 degree in radians */
115 
116 
117 double cyeqCoeffTimesSine( double coeff, double c, double Beta )
118 {
119  return coeff * sin(c * Beta);
120 }
121 
122 
123 /************************************************************************/
124 /* FUNCTIONS
125  *
126  */
127 
128 CylindricalEqualArea::CylindricalEqualArea( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double originLatitude, double falseEasting, double falseNorthing ) :
129  CoordinateSystem(),
130  es( .081819190842622 ),
131  es2( 0.0066943799901413800 ),
132  es4( 4.4814723452405e-005 ),
133  es6( 3.0000678794350e-007 ),
134  k0( 1.0 ),
135  Cyeq_a_k0( 6378137.0 ),
136  two_k0( 2.0 ),
137  c0( .0022392088624809 ),
138  c1( 2.8830839728915e-006 ),
139  c2( 5.0331826636906e-009 ),
140  Cyeq_Origin_Long( 0.0 ),
141  Cyeq_Origin_Lat( 0.0 ),
142  Cyeq_False_Easting( 0.0 ),
143  Cyeq_False_Northing( 0.0 ),
144  Cyeq_Max_Easting( 20037509.0 ),
145  Cyeq_Min_Easting( -20037509.0 ),
146  Cyeq_Delta_Northing( 6363886.0 )
147 {
148 /*
149  * The constructor receives the ellipsoid parameters and
150  * Cylindrical Equal Area projcetion parameters as inputs, and sets the corresponding
151  * state variables. If any errors occur, an exception is thrown with a description
152  * of the error.
153  *
154  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
155  * ellipsoidFlattening : Flattening of ellipsoid (input)
156  * centralMeridian : Longitude in radians at the center of (input)
157  * the projection
158  * originLatitude : Latitude in radians at which the (input)
159  * point scale factor is 1.0
160  * falseEasting : A coordinate value in meters assigned to the
161  * central meridian of the projection. (input)
162  * falseNorthing : A coordinate value in meters assigned to the
163  * origin latitude of the projection (input)
164  */
165 
166  double Sin_Cyeq_Origin_Lat;
167  double inv_f = 1 / ellipsoidFlattening;
168 
169  if (ellipsoidSemiMajorAxis <= 0.0)
170  { /* Semi-major axis must be greater than zero */
171  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
172  }
173  if ((inv_f < 250) || (inv_f > 350))
174  { /* Inverse flattening must be between 250 and 350 */
175  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
176  }
177  if ((originLatitude < -PI_OVER_2) || (originLatitude > PI_OVER_2))
178  { /* origin latitude out of range */
179  throw CoordinateConversionException( ErrorMessages::originLatitude );
180  }
181  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
182  { /* origin longitude out of range */
183  throw CoordinateConversionException( ErrorMessages::centralMeridian );
184  }
185 
186  semiMajorAxis = ellipsoidSemiMajorAxis;
187  flattening = ellipsoidFlattening;
188 
189  Cyeq_Origin_Lat = originLatitude;
190  if (centralMeridian > PI)
191  centralMeridian -= TWO_PI;
192  Cyeq_Origin_Long = centralMeridian;
193  Cyeq_False_Northing = falseNorthing;
194  Cyeq_False_Easting = falseEasting;
195  es2 = 2 * flattening - flattening * flattening;
196  es4 = es2 * es2;
197  es6 = es4 * es2;
198  es = sqrt(es2);
199  c0 = es2 / 3.0 + 31.0 * es4 / 180.0 + 517.0 * es6 / 5040.0;
200  c1 = 23.0 * es4 / 360.0 + 251.0 * es6 / 3780.0;
201  c2 = 761.0 * es6 / 45360.0;
202  Sin_Cyeq_Origin_Lat = sin(Cyeq_Origin_Lat);
203  k0 = cos(Cyeq_Origin_Lat) / sqrt(1.0 - es2 * Sin_Cyeq_Origin_Lat * Sin_Cyeq_Origin_Lat);
204  Cyeq_a_k0 = semiMajorAxis * k0;
205  two_k0 = 2.0 * k0;
206 
207  if (Cyeq_Origin_Long > 0)
208  {
209  GeodeticCoordinates gcMax( CoordinateType::geodetic, Cyeq_Origin_Long - PI - ONE, PI_OVER_2 );
210  MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcMax );
211  Cyeq_Max_Easting = tempCoordinates->easting();
212  delete tempCoordinates;
213 
214  GeodeticCoordinates gcMin( CoordinateType::geodetic, Cyeq_Origin_Long - PI, PI_OVER_2 );
215  tempCoordinates = convertFromGeodetic( &gcMin );
216  Cyeq_Min_Easting = tempCoordinates->easting();
217  delete tempCoordinates;
218 
219  GeodeticCoordinates gcDelta( CoordinateType::geodetic, PI, PI_OVER_2 );
220  tempCoordinates = convertFromGeodetic( &gcDelta );
221  Cyeq_Delta_Northing = tempCoordinates->northing();
222  delete tempCoordinates;
223  }
224  else if (Cyeq_Origin_Long < 0)
225  {
226  GeodeticCoordinates gcMax( CoordinateType::geodetic, Cyeq_Origin_Long + PI, PI_OVER_2 );
227  MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcMax );
228  Cyeq_Max_Easting = tempCoordinates->easting();
229  delete tempCoordinates;
230 
231  GeodeticCoordinates gcMin( CoordinateType::geodetic, Cyeq_Origin_Long + PI + ONE, PI_OVER_2 );
232  tempCoordinates = convertFromGeodetic( &gcMin );
233  Cyeq_Min_Easting = tempCoordinates->easting();
234  delete tempCoordinates;
235 
236  GeodeticCoordinates gcDelta( CoordinateType::geodetic, PI, PI_OVER_2 );
237  tempCoordinates = convertFromGeodetic( &gcDelta );
238  Cyeq_Delta_Northing = tempCoordinates->northing();
239  delete tempCoordinates;
240  }
241  else
242  {
243  GeodeticCoordinates gcTemp( CoordinateType::geodetic, PI, PI_OVER_2 );
244  MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcTemp );
245  Cyeq_Max_Easting = tempCoordinates->easting();
246  Cyeq_Delta_Northing = tempCoordinates->northing();
247  delete tempCoordinates;
248  Cyeq_Min_Easting = -Cyeq_Max_Easting;
249  }
250 
251  if(Cyeq_False_Northing)
252  Cyeq_Delta_Northing -= Cyeq_False_Northing;
253  if (Cyeq_Delta_Northing < 0)
254  Cyeq_Delta_Northing = -Cyeq_Delta_Northing;
255 
256  if(Cyeq_False_Easting)
257  {
258  Cyeq_Min_Easting -= Cyeq_False_Easting;
259  Cyeq_Max_Easting -= Cyeq_False_Easting;
260  }
261 }
262 
263 
264 CylindricalEqualArea::CylindricalEqualArea( const CylindricalEqualArea &cea )
265 {
266  semiMajorAxis = cea.semiMajorAxis;
267  flattening = cea.flattening;
268  es = cea.es;
269  es2 = cea.es2;
270  es4 = cea.es4;
271  es6 = cea.es6;
272  k0 = cea.k0;
273  Cyeq_a_k0 = cea.Cyeq_a_k0;
274  two_k0 = cea.two_k0;
275  c0 = cea.c0;
276  c1 = cea.c1;
277  c2 = cea.c2;
278  Cyeq_Origin_Long = cea.Cyeq_Origin_Long;
279  Cyeq_Origin_Lat = cea.Cyeq_Origin_Lat;
280  Cyeq_False_Easting = cea.Cyeq_False_Easting;
281  Cyeq_False_Northing = cea.Cyeq_False_Northing;
282  Cyeq_Max_Easting = cea.Cyeq_Max_Easting;
283  Cyeq_Min_Easting = cea.Cyeq_Min_Easting;
284  Cyeq_Delta_Northing = cea.Cyeq_Delta_Northing;
285 }
286 
287 
288 CylindricalEqualArea::~CylindricalEqualArea()
289 {
290 }
291 
292 
293 CylindricalEqualArea& CylindricalEqualArea::operator=( const CylindricalEqualArea &cea )
294 {
295  if( this != &cea )
296  {
297  semiMajorAxis = cea.semiMajorAxis;
298  flattening = cea.flattening;
299  es = cea.es;
300  es2 = cea.es2;
301  es4 = cea.es4;
302  es6 = cea.es6;
303  k0 = cea.k0;
304  Cyeq_a_k0 = cea.Cyeq_a_k0;
305  two_k0 = cea.two_k0;
306  c0 = cea.c0;
307  c1 = cea.c1;
308  c2 = cea.c2;
309  Cyeq_Origin_Long = cea.Cyeq_Origin_Long;
310  Cyeq_Origin_Lat = cea.Cyeq_Origin_Lat;
311  Cyeq_False_Easting = cea.Cyeq_False_Easting;
312  Cyeq_False_Northing = cea.Cyeq_False_Northing;
313  Cyeq_Max_Easting = cea.Cyeq_Max_Easting;
314  Cyeq_Min_Easting = cea.Cyeq_Min_Easting;
315  Cyeq_Delta_Northing = cea.Cyeq_Delta_Northing;
316  }
317 
318  return *this;
319 }
320 
321 
322 MapProjection4Parameters* CylindricalEqualArea::getParameters() const
323 {
324 /*
325  * The function getParameters returns the current ellipsoid
326  * parameters, and Cylindrical Equal Area projection parameters.
327  *
328  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
329  * ellipsoidFlattening : Flattening of ellipsoid (output)
330  * centralMeridian : Longitude in radians at the center of (output)
331  * the projection
332  * originLatitude : Latitude in radians at which the (output)
333  * point scale factor is 1.0
334  * falseEasting : A coordinate value in meters assigned to the
335  * central meridian of the projection. (output)
336  * falseNorthing : A coordinate value in meters assigned to the
337  * origin latitude of the projection (output)
338  */
339 
340  return new MapProjection4Parameters( CoordinateType::cylindricalEqualArea, Cyeq_Origin_Long, Cyeq_Origin_Lat, Cyeq_False_Easting, Cyeq_False_Northing );
341 }
342 
343 
344 MSP::CCS::MapProjectionCoordinates* CylindricalEqualArea::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
345 {
346 /*
347  * The function convertFromGeodetic converts geodetic (latitude and
348  * longitude) coordinates to Cylindrical Equal Area projection (easting and northing)
349  * coordinates, according to the current ellipsoid and Cylindrical Equal Area projection
350  * parameters. If any errors occur, an exception is thrown with a description
351  * of the error.
352  *
353  * longitude : Longitude (lambda) in radians (input)
354  * latitude : Latitude (phi) in radians (input)
355  * easting : Easting (X) in meters (output)
356  * northing : Northing (Y) in meters (output)
357  */
358 
359  double dlam; /* Longitude - Central Meridan */
360  double qq;
361  double x;
362 
363  double longitude = geodeticCoordinates->longitude();
364  double latitude = geodeticCoordinates->latitude();
365  double sin_lat = sin(latitude);
366 
367  if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
368  { /* Latitude out of range */
369  throw CoordinateConversionException( ErrorMessages::latitude );
370  }
371  if ((longitude < -PI) || (longitude > TWO_PI))
372  { /* Longitude out of range */
373  throw CoordinateConversionException( ErrorMessages::longitude );
374  }
375 
376  dlam = longitude - Cyeq_Origin_Long;
377  if (dlam > PI)
378  {
379  dlam -= TWO_PI;
380  }
381  if (dlam < -PI)
382  {
383  dlam += TWO_PI;
384  }
385  x = es * sin_lat;
386  qq = cyleqarQ( sin_lat, x );
387 
388  double easting = Cyeq_a_k0 * dlam + Cyeq_False_Easting;
389  double northing = semiMajorAxis * qq / two_k0 + Cyeq_False_Northing;
390 
391  return new MapProjectionCoordinates( CoordinateType::cylindricalEqualArea, easting, northing );
392 }
393 
394 
395 MSP::CCS::GeodeticCoordinates* CylindricalEqualArea::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
396 {
397 /*
398  * The function convertToGeodetic converts
399  * Cylindrical Equal Area projection (easting and northing) coordinates
400  * to geodetic (latitude and longitude) coordinates, according to the
401  * current ellipsoid and Cylindrical Equal Area projection
402  * coordinates. If any errors occur, an exception is thrown with a description
403  * of the error.
404  *
405  * easting : Easting (X) in meters (input)
406  * northing : Northing (Y) in meters (input)
407  * longitude : Longitude (lambda) in radians (output)
408  * latitude : Latitude (phi) in radians (output)
409  */
410 
411  double sin2beta, sin4beta, sin6beta;
412  double dx; /* Delta easting - Difference in easting (easting-FE) */
413  double dy; /* Delta northing - Difference in northing (northing-FN) */
414  double qp;
415  double beta;
416  double sin_lat = sin(PI_OVER_2);
417  double i;
418  double x;
419 
420  double easting = mapProjectionCoordinates->easting();
421  double northing = mapProjectionCoordinates->northing();
422 
423  if ((easting < (Cyeq_False_Easting + Cyeq_Min_Easting))
424  || (easting > (Cyeq_False_Easting + Cyeq_Max_Easting)))
425  { /* Easting out of range */
426  throw CoordinateConversionException( ErrorMessages::easting );
427  }
428  if ((northing < (Cyeq_False_Northing - fabs(Cyeq_Delta_Northing)))
429  || (northing > (Cyeq_False_Northing + fabs(Cyeq_Delta_Northing))))
430  { /* Northing out of range */
431  throw CoordinateConversionException( ErrorMessages::northing );
432  }
433 
434  dy = northing - Cyeq_False_Northing;
435  dx = easting - Cyeq_False_Easting;
436  x = es * sin_lat;
437  qp = cyleqarQ( sin_lat, x );
438  i = two_k0 * dy / (semiMajorAxis * qp);
439  if (i > 1.0)
440  i = 1.0;
441  else if (i < -1.0)
442  i = -1.0;
443  beta = asin(i);
444  sin2beta = cyeqCoeffTimesSine(c0, 2.0, beta);
445  sin4beta = cyeqCoeffTimesSine(c1, 4.0, beta);
446  sin6beta = cyeqCoeffTimesSine(c2, 6.0, beta);
447 
448  double latitude = beta + sin2beta + sin4beta + sin6beta;
449  double longitude = Cyeq_Origin_Long + dx / Cyeq_a_k0;
450 
451  if (latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
452  latitude = PI_OVER_2;
453  else if (latitude < -PI_OVER_2)
454  latitude = -PI_OVER_2;
455 
456  if (longitude > PI)
457  longitude -= TWO_PI;
458  if (longitude < -PI)
459  longitude += TWO_PI;
460 
461  if (longitude > PI) /* force distorted values to 180, -180 degrees */
462  longitude = PI;
463  else if (longitude < -PI)
464  longitude = -PI;
465 
466  return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
467 }
468 
469 
470 double CylindricalEqualArea::cyleqarQ( double slat, double x )
471 {
472  return (1.0-es2)*(slat/(1.0-x*x)-(1.0/(2.0*es)) * log((1.0-x)/(1.0+x)));
473 }
474 
475 
476 
477 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:01 for GeographicTranslator by doxygen 1.8.2