UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
Eckert4.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: ECKERT4
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates
9  * (latitude and longitude in radians) and Eckert IV projection coordinates
10  * (easting and northing in meters). This projection employs a spherical
11  * Earth model. The spherical radius used is the radius of the sphere
12  * having the same area as the ellipsoid.
13  *
14  * ERROR HANDLING
15  *
16  * This component checks parameters for valid values. If an invalid value
17  * is found, the error code is combined with the current error code using
18  * the bitwise or. This combining allows multiple error codes to be
19  * returned. The possible error codes are:
20  *
21  * ECK4_NO_ERROR : No errors occurred in function
22  * ECK4_LAT_ERROR : Latitude outside of valid range
23  * (-90 to 90 degrees)
24  * ECK4_LON_ERROR : Longitude outside of valid range
25  * (-180 to 360 degrees)
26  * ECK4_EASTING_ERROR : Easting outside of valid range
27  * (False_Easting +/- ~17,000,000 m,
28  * depending on ellipsoid parameters)
29  * ECK4_NORTHING_ERROR : Northing outside of valid range
30  * (False_Northing +/- 0 to 8,000,000 m,
31  * depending on ellipsoid parameters)
32  * ECK4_CENT_MER_ERROR : Central_Meridian outside of valid range
33  * (-180 to 360 degrees)
34  * ECK4_A_ERROR : Semi-major axis less than or equal to zero
35  * ECK4_INV_F_ERROR : Inverse flattening outside of valid range
36  * (250 to 350)
37  *
38  * REUSE NOTES
39  *
40  * ECKERT4 is intended for reuse by any application that performs a
41  * Eckert IV projection or its inverse.
42  *
43  * REFERENCES
44  *
45  * Further information on ECKERT4 can be found in the Reuse Manual.
46  *
47  * ECKERT4 originated from : U.S. Army Topographic Engineering Center
48  * Geospatial Information Division
49  * 7701 Telegraph Road
50  * Alexandria, VA 22310-3864
51  *
52  * LICENSES
53  *
54  * None apply to this component.
55  *
56  * RESTRICTIONS
57  *
58  * ECKERT4 has no restrictions.
59  *
60  * ENVIRONMENT
61  *
62  * ECKERT4 was tested and certified in the following environments:
63  *
64  * 1. Solaris 2.5 with GCC 2.8.1
65  * 2. MS Windows 95 with MS Visual C++ 6
66  *
67  * MODIFICATIONS
68  *
69  * Date Description
70  * ---- -----------
71  * 04-16-99 Original Code
72  * 03-06-07 Original C++ Code
73  *
74  */
75 
76 
77 /***************************************************************************/
78 /*
79  * INCLUDES
80  */
81 
82 #include <math.h>
83 #include "Eckert4.h"
84 #include "MapProjection3Parameters.h"
85 #include "MapProjectionCoordinates.h"
86 #include "GeodeticCoordinates.h"
87 #include "CoordinateConversionException.h"
88 #include "ErrorMessages.h"
89 
90 /*
91  * math.h - Standard C math library
92  * Eckert4.h - Is for prototype error checking.
93  * MapProjectionCoordinates.h - defines map projection coordinates
94  * GeodeticCoordinates.h - defines geodetic coordinates
95  * CoordinateConversionException.h - Exception handler
96  * ErrorMessages.h - Contains exception messages
97  */
98 
99 
100 using namespace MSP::CCS;
101 
102 
103 /***************************************************************************/
104 /* DEFINES
105  *
106  */
107 
108 const double PI = 3.14159265358979323e0; /* PI */
109 const double PI_OVER_2 = ( PI / 2.0);
110 const double TWO_PI = (2.0 * PI) ;
111 const double two_PLUS_PI_OVER_2 = (2.0 + PI / 2.0);
112 
113 
114 double calculateNum( double theta, double sinTheta, double cosTheta )
115 {
116  return theta + sinTheta * cosTheta + 2.0 * sinTheta;
117 }
118 
119 
120 /************************************************************************/
121 /* FUNCTIONS
122  *
123  */
124 
125 Eckert4::Eckert4( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double falseEasting, double falseNorthing ) :
126  CoordinateSystem(),
127  es2( 0.0066943799901413800 ),
128  es4( 4.4814723452405e-005 ),
129  es6( 3.0000678794350e-007 ),
130  Ra0( 2690082.6043273 ),
131  Ra1( 8451143.5741087 ),
132  Eck4_Origin_Long( 0.0 ),
133  Eck4_False_Easting( 0.0 ),
134  Eck4_False_Northing( 0.0 ),
135  Eck4_Delta_Northing( 8451144.0 ),
136  Eck4_Max_Easting( 16902288.0 ),
137  Eck4_Min_Easting( -16902288.0 )
138 {
139 /*
140  * The constructor receives the ellipsoid parameters and
141  * projection parameters as inputs, and sets the corresponding state
142  * variables. If any errors occur, an exception is thrown with a description
143  * of the error.
144  *
145  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
146  * ellipsoidFlattening : Flattening of ellipsoid (input)
147  * centralMeridian : Longitude in radians at the center of (input)
148  * the projection
149  * falseEasting : A coordinate value in meters assigned to the
150  * central meridian of the projection. (input)
151  * falseNorthing : A coordinate value in meters assigned to the
152  * origin latitude of the projection (input)
153  */
154 
155  double Ra; /* Spherical radius */
156  double inv_f = 1 / ellipsoidFlattening;
157 
158  if (ellipsoidSemiMajorAxis <= 0.0)
159  { /* Semi-major axis must be greater than zero */
160  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
161  }
162  if ((inv_f < 250) || (inv_f > 350))
163  { /* Inverse flattening must be between 250 and 350 */
164  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
165  }
166  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
167  { /* origin longitude out of range */
168  throw CoordinateConversionException( ErrorMessages::centralMeridian );
169  }
170 
171  semiMajorAxis = ellipsoidSemiMajorAxis;
172  flattening = ellipsoidFlattening;
173 
174  es2 = 2 * flattening - flattening * flattening;
175  es4 = es2 * es2;
176  es6 = es4 * es2;
177  /* spherical radius */
178  Ra = semiMajorAxis *
179  (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 / 3024.0);
180  Ra0 = 0.4222382 * Ra;
181  Ra1 = 1.3265004 * Ra;
182  if (centralMeridian > PI)
183  centralMeridian -= TWO_PI;
184  Eck4_Origin_Long = centralMeridian;
185  Eck4_False_Easting = falseEasting;
186  Eck4_False_Northing = falseNorthing;
187  if (Eck4_Origin_Long > 0)
188  {
189  Eck4_Max_Easting = 16808386.0;
190  Eck4_Min_Easting = -16902288.0;
191  }
192  else if (Eck4_Origin_Long < 0)
193  {
194  Eck4_Max_Easting = 16902288.0;
195  Eck4_Min_Easting = -16808386.0;
196  }
197  else
198  {
199  Eck4_Max_Easting = 16902288.0;
200  Eck4_Min_Easting = -16902288.0;
201  }
202 }
203 
204 
205 Eckert4::Eckert4( const Eckert4 &e )
206 {
207  semiMajorAxis = e.semiMajorAxis;
208  flattening = e.flattening;
209  es2 = e.es2;
210  es4 = e.es4;
211  es6 = e.es6;
212  Ra0 = e.Ra0;
213  Ra1 = e.Ra1;
214  Eck4_Origin_Long = e.Eck4_Origin_Long;
215  Eck4_False_Easting = e.Eck4_False_Easting;
216  Eck4_False_Northing = e.Eck4_False_Northing;
217  Eck4_Delta_Northing = e.Eck4_Delta_Northing;
218  Eck4_Max_Easting = e.Eck4_Max_Easting;
219  Eck4_Min_Easting = e.Eck4_Min_Easting;
220 }
221 
222 
223 Eckert4::~Eckert4()
224 {
225 }
226 
227 
228 Eckert4& Eckert4::operator=( const Eckert4 &e )
229 {
230  if( this != &e )
231  {
232  semiMajorAxis = e.semiMajorAxis;
233  flattening = e.flattening;
234  es2 = e.es2;
235  es4 = e.es4;
236  es6 = e.es6;
237  Ra0 = e.Ra0;
238  Ra1 = e.Ra1;
239  Eck4_Origin_Long = e.Eck4_Origin_Long;
240  Eck4_False_Easting = e.Eck4_False_Easting;
241  Eck4_False_Northing = e.Eck4_False_Northing;
242  Eck4_Delta_Northing = e.Eck4_Delta_Northing;
243  Eck4_Max_Easting = e.Eck4_Max_Easting;
244  Eck4_Min_Easting = e.Eck4_Min_Easting;
245  }
246 
247  return *this;
248 }
249 
250 
251 MapProjection3Parameters* Eckert4::getParameters() const
252 {
253 /*
254  * The function getParameters returns the current ellipsoid
255  * parameters and Eckert IV projection parameters.
256  *
257  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
258  * ellipsoidFlattening : Flattening of ellipsoid (output)
259  * centralMeridian : Longitude in radians at the center of (output)
260  * the projection
261  * falseEasting : A coordinate value in meters assigned to the
262  * central meridian of the projection. (output)
263  * falseNorthing : A coordinate value in meters assigned to the
264  * origin latitude of the projection (output)
265  */
266 
267  return new MapProjection3Parameters( CoordinateType::eckert4, Eck4_Origin_Long, Eck4_False_Easting, Eck4_False_Northing );
268 }
269 
270 
271 MSP::CCS::MapProjectionCoordinates* Eckert4::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
272 {
273 /*
274  * The function convertFromGeodetic converts geodetic (latitude and
275  * longitude) coordinates to Eckert IV projection (easting and northing)
276  * coordinates, according to the current ellipsoid, spherical radius and
277  * Eckert IV projection parameters.
278  * If any errors occur, an exception is thrown with a description
279  * of the error.
280  *
281  * longitude : Longitude (lambda) in radians (input)
282  * latitude : Latitude (phi) in radians (input)
283  * easting : Easting (X) in meters (output)
284  * northing : Northing (Y) in meters (output)
285  */
286 
287  double sin_theta, cos_theta;
288  double num;
289  double dlam; /* Longitude - Central Meridan */
290  double delta_theta = 1.0;
291  double dt_tolerance = 4.85e-10; /* approximately 1/1000th of
292  an arc second or 1/10th meter */
293  int count = 200;
294 
295  double longitude = geodeticCoordinates->longitude();
296  double latitude = geodeticCoordinates->latitude();
297  double slat = sin(latitude);
298  double theta = latitude / 2.0;
299 
300  if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
301  { /* Latitude out of range */
302  throw CoordinateConversionException( ErrorMessages::latitude );
303  }
304  if ((longitude < -PI) || (longitude > TWO_PI))
305  { /* Longitude out of range */
306  throw CoordinateConversionException( ErrorMessages::longitude );
307  }
308 
309  dlam = longitude - Eck4_Origin_Long;
310  if (dlam > PI)
311  {
312  dlam -= TWO_PI;
313  }
314  if (dlam < -PI)
315  {
316  dlam += TWO_PI;
317  }
318  while (fabs(delta_theta) > dt_tolerance && count)
319  {
320  sin_theta = sin(theta);
321  cos_theta = cos(theta);
322  num = calculateNum(theta, sin_theta, cos_theta);
323  delta_theta = -(num - two_PLUS_PI_OVER_2 * slat) /
324  (2.0 * cos_theta * (1.0 + cos_theta));
325  theta += delta_theta;
326  count --;
327  }
328 
329  if(!count)
330  throw CoordinateConversionException( ErrorMessages::northing );
331 
332  double easting = Ra0 * dlam * (1.0 + cos(theta)) + Eck4_False_Easting;
333  double northing = Ra1 * sin(theta) + Eck4_False_Northing;
334 
335  return new MapProjectionCoordinates( CoordinateType::eckert4, easting, northing );
336 }
337 
338 
339 MSP::CCS::GeodeticCoordinates* Eckert4::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
340 {
341 /*
342  * The function convertToGeodetic converts Eckert IV projection
343  * (easting and northing) coordinates to geodetic (latitude and longitude)
344  * coordinates, according to the current ellipsoid, spherical radius and
345  * Eckert IV projection coordinates.
346  * If any errors occur, an exception is thrown with a description
347  * of the error.
348  *
349  * easting : Easting (X) in meters (input)
350  * northing : Northing (Y) in meters (input)
351  * longitude : Longitude (lambda) in radians (output)
352  * latitude : Latitude (phi) in radians (output)
353  */
354 
355  double theta;
356  double sin_theta, cos_theta;
357  double num;
358  double dx, dy;
359  double i;
360 
361  double easting = mapProjectionCoordinates->easting();
362  double northing = mapProjectionCoordinates->northing();
363 
364  if ((easting < (Eck4_False_Easting + Eck4_Min_Easting))
365  || (easting > (Eck4_False_Easting + Eck4_Max_Easting)))
366  { /* Easting out of range */
367  throw CoordinateConversionException( ErrorMessages::easting );
368  }
369  if ((northing < (Eck4_False_Northing - Eck4_Delta_Northing))
370  || (northing > (Eck4_False_Northing + Eck4_Delta_Northing)))
371  { /* Northing out of range */
372  throw CoordinateConversionException( ErrorMessages::northing );
373  }
374 
375  dy = northing - Eck4_False_Northing;
376  dx = easting - Eck4_False_Easting;
377  i = dy / Ra1;
378  if (i > 1.0)
379  i = 1.0;
380  else if (i < -1.0)
381  i = -1.0;
382 
383  theta = asin(i);
384  sin_theta = sin(theta);
385  cos_theta = cos(theta);
386  num = calculateNum(theta, sin_theta, cos_theta);
387 
388  double latitude = asin(num / two_PLUS_PI_OVER_2);
389  double longitude = Eck4_Origin_Long + dx / (Ra0 * (1 + cos_theta));
390 
391  if (latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
392  latitude = PI_OVER_2;
393  else if (latitude < -PI_OVER_2)
394  latitude = -PI_OVER_2;
395 
396  if (longitude > PI)
397  longitude -= TWO_PI;
398  if (longitude < -PI)
399  longitude += TWO_PI;
400 
401  if (longitude > PI) /* force distorted values to 180, -180 degrees */
402  longitude = PI;
403  else if (longitude < -PI)
404  longitude = -PI;
405 
406  return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
407 }
408 
409 
410 
411 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:01 for GeographicTranslator by doxygen 1.8.2