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