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