UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
AzimuthalEquidistant.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: AZIMUTHAL EQUIDISTANT
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates
9  * (latitude and longitude in radians) and Azimuthal Equidistant
10  * projection coordinates (easting and northing in meters). This projection
11  * employs a spherical Earth model. The spherical radius used is the radius of
12  * the sphere 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  * AZEQ_NO_ERROR : No errors occurred in function
22  * AZEQ_LAT_ERROR : Latitude outside of valid range
23  * (-90 to 90 degrees)
24  * AZEQ_LON_ERROR : Longitude outside of valid range
25  * (-180 to 360 degrees)
26  * AZEQ_EASTING_ERROR : Easting outside of valid range
27  * (depends on ellipsoid and projection
28  * parameters)
29  * AZEQ_NORTHING_ERROR : Northing outside of valid range
30  * (depends on ellipsoid and projection
31  * parameters)
32  * AZEQ_ORIGIN_LAT_ERROR : Origin latitude outside of valid range
33  * (-90 to 90 degrees)
34  * AZEQ_CENT_MER_ERROR : Central meridian outside of valid range
35  * (-180 to 360 degrees)
36  * AZEQ_A_ERROR : Semi-major axis less than or equal to zero
37  * AZEQ_INV_F_ERROR : Inverse flattening outside of valid range
38  * (250 to 350)
39  * AZEQ_PROJECTION_ERROR : Point is plotted as a circle of radius PI * Ra
40  *
41  *
42  * REUSE NOTES
43  *
44  * AZIMUTHAL EQUIDISTANT is intended for reuse by any application that
45  * performs an Azimuthal Equidistant projection or its inverse.
46  *
47  * REFERENCES
48  *
49  * Further information on AZIMUTHAL EQUIDISTANT can be found in the Reuse Manual.
50  *
51  * AZIMUTHAL EQUIDISTANT originated from:
52  * 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  * AZIMUTHAL EQUIDISTANT has no restrictions.
64  *
65  * ENVIRONMENT
66  *
67  * AZIMUTHAL EQUIDISTANT was tested and certified in the following environments:
68  *
69  * 1. Solaris 2.5 with GCC, version 2.8.1
70  * 2. MSDOS with MS Visual C++, version 6
71  *
72  * MODIFICATIONS
73  *
74  * Date Description
75  * ---- -----------
76  * 05-19-00 Original Code
77  * 03-08-07 Original C++ Code
78  *
79  *
80  */
81 
82 
83 /***************************************************************************/
84 /*
85  * INCLUDES
86  */
87 
88 #include <math.h>
89 #include "AzimuthalEquidistant.h"
90 #include "MapProjection4Parameters.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  * AzimuthalEquidistant.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 AzimuthalEquidistant::AzimuthalEquidistant(
126  double ellipsoidSemiMajorAxis,
127  double ellipsoidFlattening,
128  double centralMeridian,
129  double originLatitude,
130  double falseEasting,
131  double falseNorthing ) :
132  CoordinateSystem(),
133  Ra( 6371007.1810824 ),
134  Sin_Azeq_Origin_Lat( 0.0 ),
135  Cos_Azeq_Origin_Lat( 1.0 ),
136  Azeq_Origin_Long( 0.0 ),
137  Azeq_Origin_Lat( 0.0 ),
138  Azeq_False_Easting( 0.0 ),
139  Azeq_False_Northing( 0.0 ),
140  abs_Azeq_Origin_Lat( 0.0 ),
141  Azeq_Delta_Northing( 19903915.0 ),
142  Azeq_Delta_Easting( 19903915.0 )
143 {
144 /*
145  * The function setParameters receives the ellipsoid parameters and
146  * projection parameters as inputs, and sets the corresponding state
147  * variables. If any errors occur, an exception is thrown with a description
148  * of the error.
149  *
150  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
151  * ellipsoidFlattening : Flattening of ellipsoid (input)
152  * centralMeridian : Longitude in radians at the center of (input)
153  * the projection
154  * originLatitude : Latitude in radians at which the (input)
155  * point scale factor is 1.0
156  * falseEasting : A coordinate value in meters assigned to the
157  * central meridian of the projection. (input)
158  * falseNorthing : A coordinate value in meters assigned to the
159  * origin latitude of the projection (input)
160  */
161 
162  double es2, es4, es6;
163  double inv_f = 1 / ellipsoidFlattening;
164 
165  if (ellipsoidSemiMajorAxis <= 0.0)
166  { /* Semi-major axis must be greater than zero */
167  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
168  }
169  if ((inv_f < 250) || (inv_f > 350))
170  { /* Inverse flattening must be between 250 and 350 */
171  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
172  }
173  if ((originLatitude < -PI_OVER_2) || (originLatitude > PI_OVER_2))
174  { /* origin latitude out of range */
175  throw CoordinateConversionException( ErrorMessages::originLatitude );
176  }
177  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
178  { /* origin longitude out of range */
179  throw CoordinateConversionException( ErrorMessages::centralMeridian );
180  }
181 
182  semiMajorAxis = ellipsoidSemiMajorAxis;
183  flattening = ellipsoidFlattening;
184 
185  es2 = 2 * flattening - flattening * flattening;
186  es4 = es2 * es2;
187  es6 = es4 * es2;
188  /* spherical radius */
189  Ra = semiMajorAxis *
190  (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 / 3024.0);
191  Azeq_Origin_Lat = originLatitude;
192  Sin_Azeq_Origin_Lat = sin(Azeq_Origin_Lat);
193  Cos_Azeq_Origin_Lat = cos(Azeq_Origin_Lat);
194  abs_Azeq_Origin_Lat = fabs(Azeq_Origin_Lat);
195  if (centralMeridian > PI)
196  centralMeridian -= TWO_PI;
197  Azeq_Origin_Long = centralMeridian;
198  Azeq_False_Northing = falseNorthing;
199  Azeq_False_Easting = falseEasting;
200 
201  if (fabs(abs_Azeq_Origin_Lat - PI_OVER_2) < 1.0e-10)
202  {
203  Azeq_Delta_Northing = 20015110.0;
204  Azeq_Delta_Easting = 20015110.0;
205  }
206  else if (abs_Azeq_Origin_Lat >= 1.0e-10)
207  {
208  if (Azeq_Origin_Long > 0.0)
209  {
210  GeodeticCoordinates gcDelta( CoordinateType::geodetic, Azeq_Origin_Long - PI + ONE, -Azeq_Origin_Lat );
211  MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcDelta );
212  Azeq_Delta_Easting = tempCoordinates->easting();
213  delete tempCoordinates;
214  }
215  else
216  {
217  GeodeticCoordinates gcDelta(CoordinateType::geodetic, Azeq_Origin_Long + PI - ONE, -Azeq_Origin_Lat );
218  MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcDelta );
219  Azeq_Delta_Easting = tempCoordinates->easting();
220  delete tempCoordinates;
221  }
222 
223  if(Azeq_False_Easting)
224  Azeq_Delta_Easting -= Azeq_False_Easting;
225  if (Azeq_Delta_Easting < 0)
226  Azeq_Delta_Easting = -Azeq_Delta_Easting;
227 
228  Azeq_Delta_Northing = 19903915.0;
229  }
230  else
231  {
232  Azeq_Delta_Northing = 19903915.0;
233  Azeq_Delta_Easting = 19903915.0;
234  }
235 }
236 
237 
238 AzimuthalEquidistant::AzimuthalEquidistant( const AzimuthalEquidistant &ae )
239 {
240  semiMajorAxis = ae.semiMajorAxis;
241  flattening = ae.flattening;
242  Ra = ae.Ra;
243  Sin_Azeq_Origin_Lat = ae.Sin_Azeq_Origin_Lat;
244  Cos_Azeq_Origin_Lat = ae.Cos_Azeq_Origin_Lat;
245  Azeq_Origin_Long = ae.Azeq_Origin_Long;
246  Azeq_Origin_Lat = ae.Azeq_Origin_Lat;
247  Azeq_False_Easting = ae.Azeq_False_Easting;
248  Azeq_False_Northing = ae.Azeq_False_Northing;
249  abs_Azeq_Origin_Lat = ae.abs_Azeq_Origin_Lat;
250  Azeq_Delta_Northing = ae.Azeq_Delta_Northing;
251  Azeq_Delta_Easting = ae.Azeq_Delta_Easting;
252 }
253 
254 
255 AzimuthalEquidistant::~AzimuthalEquidistant()
256 {
257 }
258 
259 
260 AzimuthalEquidistant& AzimuthalEquidistant::operator=( const AzimuthalEquidistant &ae )
261 {
262  if( this != &ae )
263  {
264  semiMajorAxis = ae.semiMajorAxis;
265  flattening = ae.flattening;
266  Ra = ae.Ra;
267  Sin_Azeq_Origin_Lat = ae.Sin_Azeq_Origin_Lat;
268  Cos_Azeq_Origin_Lat = ae.Cos_Azeq_Origin_Lat;
269  Azeq_Origin_Long = ae.Azeq_Origin_Long;
270  Azeq_Origin_Lat = ae.Azeq_Origin_Lat;
271  Azeq_False_Easting = ae.Azeq_False_Easting;
272  Azeq_False_Northing = ae.Azeq_False_Northing;
273  abs_Azeq_Origin_Lat = ae.abs_Azeq_Origin_Lat;
274  Azeq_Delta_Northing = ae.Azeq_Delta_Northing;
275  Azeq_Delta_Easting = ae.Azeq_Delta_Easting;
276  }
277 
278  return *this;
279 }
280 
281 
282 MapProjection4Parameters* AzimuthalEquidistant::getParameters() const
283 {
284 /*
285  * The function getParameters returns the current ellipsoid
286  * parameters and Azimuthal Equidistant projection parameters.
287  *
288  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
289  * ellipsoidFlattening : Flattening of ellipsoid (output)
290  * centralMeridian : Longitude in radians at the center of (output)
291  * the projection
292  * originLatitude : Latitude in radians at which the (output)
293  * point scale factor is 1.0
294  * falseEasting : A coordinate value in meters assigned to the
295  * central meridian of the projection. (output)
296  * falseNorthing : A coordinate value in meters assigned to the
297  * origin latitude of the projection (output)
298  */
299 
300  return new MapProjection4Parameters( CoordinateType::azimuthalEquidistant, Azeq_Origin_Long, Azeq_Origin_Lat, Azeq_False_Easting, Azeq_False_Northing );
301 }
302 
303 
304 MSP::CCS::MapProjectionCoordinates* AzimuthalEquidistant::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
305 {
306 /*
307  * The function convertFromGeodetic converts geodetic (latitude and
308  * longitude) coordinates to Azimuthal Equidistant projection (easting and northing)
309  * coordinates, according to the current ellipsoid and Azimuthal Equidistant projection
310  * parameters. If any errors occur, an exception is thrown with a description
311  * of the error.
312  *
313  * longitude : Longitude (lambda) in radians (input)
314  * latitude : Latitude (phi) in radians (input)
315  * easting : Easting (X) in meters (output)
316  * northing : Northing (Y) in meters (output)
317  */
318 
319  double dlam; /* Longitude - Central Meridan */
320  double k_prime; /* scale factor */
321  double c; /* angular distance from center */
322  double cos_c;
323  double sin_dlam, cos_dlam;
324  double Ra_kprime;
325  double Ra_PI_OVER_2_Lat;
326  double easting, northing;
327 
328  double longitude = geodeticCoordinates->longitude();
329  double latitude = geodeticCoordinates->latitude();
330  double slat = sin(latitude);
331  double clat = cos(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 - Azeq_Origin_Long;
343  if (dlam > PI)
344  {
345  dlam -= TWO_PI;
346  }
347  if (dlam < -PI)
348  {
349  dlam += TWO_PI;
350  }
351 
352  sin_dlam = sin(dlam);
353  cos_dlam = cos(dlam);
354  if (fabs(abs_Azeq_Origin_Lat - PI_OVER_2) < 1.0e-10)
355  {
356  if (Azeq_Origin_Lat >= 0.0)
357  {
358  Ra_PI_OVER_2_Lat = Ra * (PI_OVER_2 - latitude);
359  easting = Ra_PI_OVER_2_Lat * sin_dlam + Azeq_False_Easting;
360  northing = -1.0 * (Ra_PI_OVER_2_Lat * cos_dlam) + Azeq_False_Northing;
361  }
362  else
363  {
364  Ra_PI_OVER_2_Lat = Ra * (PI_OVER_2 + latitude);
365  easting = Ra_PI_OVER_2_Lat * sin_dlam + Azeq_False_Easting;
366  northing = Ra_PI_OVER_2_Lat * cos_dlam + Azeq_False_Northing;
367  }
368  }
369  else if (abs_Azeq_Origin_Lat <= 1.0e-10)
370  {
371  cos_c = clat * cos_dlam;
372  if (fabs(fabs(cos_c) - 1.0) < 1.0e-14)
373  {
374  if (cos_c >= 0.0)
375  {
376  easting = Azeq_False_Easting;
377  northing = Azeq_False_Northing;
378  }
379  else
380  {
381  /* if cos_c == -1 */
382  throw CoordinateConversionException( MSP::CCS::ErrorMessages::projection );
383  }
384  }
385  else
386  {
387  c = acos(cos_c);
388  k_prime = c / sin(c);
389  Ra_kprime = Ra * k_prime;
390  easting = Ra_kprime * clat * sin_dlam + Azeq_False_Easting;
391  northing = Ra_kprime * slat + Azeq_False_Northing;
392  }
393  }
394  else
395  {
396  cos_c = (Sin_Azeq_Origin_Lat * slat) + (Cos_Azeq_Origin_Lat * clat * cos_dlam);
397  if (fabs(fabs(cos_c) - 1.0) < 1.0e-14)
398  {
399  if (cos_c >= 0.0)
400  {
401  easting = Azeq_False_Easting;
402  northing = Azeq_False_Northing;
403  }
404  else
405  {
406  /* if cos_c == -1 */
407  throw CoordinateConversionException( MSP::CCS::ErrorMessages::projection );
408  }
409  }
410  else
411  {
412  c = acos(cos_c);
413  k_prime = c / sin(c);
414  Ra_kprime = Ra * k_prime;
415  easting = Ra_kprime * clat * sin_dlam + Azeq_False_Easting;
416  northing = Ra_kprime * (Cos_Azeq_Origin_Lat * slat - Sin_Azeq_Origin_Lat * clat * cos_dlam) + Azeq_False_Northing;
417  }
418  }
419 
420  return new MapProjectionCoordinates( CoordinateType::azimuthalEquidistant, easting, northing );
421 }
422 
423 
424 MSP::CCS::GeodeticCoordinates* AzimuthalEquidistant::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
425 {
426 /*
427  * The function convertToGeodetic converts Azimuthal_Equidistant projection
428  * (easting and northing) coordinates to geodetic (latitude and longitude)
429  * coordinates, according to the current ellipsoid and Azimuthal_Equidistant projection
430  * coordinates. If any errors occur, an exception is thrown with a description
431  * of the error.
432  *
433  * easting : Easting (X) in meters (input)
434  * northing : Northing (Y) in meters (input)
435  * longitude : Longitude (lambda) in radians (output)
436  * latitude : Latitude (phi) in radians (output)
437  */
438 
439  double dx, dy;
440  double rho; /* height above ellipsoid */
441  double c; /* angular distance from center */
442  double sin_c, cos_c, dy_sinc;
443  double longitude, latitude;
444 
445  double easting = mapProjectionCoordinates->easting();
446  double northing = mapProjectionCoordinates->northing();
447 
448  if ((easting < (Azeq_False_Easting - Azeq_Delta_Easting))
449  || (easting > (Azeq_False_Easting + Azeq_Delta_Easting)))
450  { /* Easting out of range */
451  throw CoordinateConversionException( ErrorMessages::easting );
452  }
453  if ((northing < (Azeq_False_Northing - Azeq_Delta_Northing))
454  || (northing > (Azeq_False_Northing + Azeq_Delta_Northing)))
455  { /* Northing out of range */
456  throw CoordinateConversionException( ErrorMessages::northing );
457  }
458 
459  dy = northing - Azeq_False_Northing;
460  dx = easting - Azeq_False_Easting;
461  rho = sqrt(dx * dx + dy * dy);
462  if (fabs(rho) <= 1.0e-10)
463  {
464  latitude = Azeq_Origin_Lat;
465  longitude = Azeq_Origin_Long;
466  }
467  else
468  {
469  c = rho / Ra;
470  sin_c = sin(c);
471  cos_c = cos(c);
472  dy_sinc = dy * sin_c;
473  latitude = asin((cos_c * Sin_Azeq_Origin_Lat) + ((dy_sinc * Cos_Azeq_Origin_Lat) / rho));
474  if (fabs(abs_Azeq_Origin_Lat - PI_OVER_2) < 1.0e-10)
475  {
476  if (Azeq_Origin_Lat >= 0.0)
477  longitude = Azeq_Origin_Long + atan2(dx, -dy);
478  else
479  longitude = Azeq_Origin_Long + atan2(dx, dy);
480  }
481  else
482  longitude = Azeq_Origin_Long + atan2((dx * sin_c), ((rho * Cos_Azeq_Origin_Lat * cos_c) - (dy_sinc * Sin_Azeq_Origin_Lat)));
483  }
484 
485  if (latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
486  latitude = PI_OVER_2;
487  else if (latitude < -PI_OVER_2)
488  latitude = -PI_OVER_2;
489 
490  if (longitude > PI)
491  longitude -= TWO_PI;
492  if (longitude < -PI)
493  longitude += TWO_PI;
494 
495  if (longitude > PI) /* force distorted values to 180, -180 degrees */
496  longitude = PI;
497  else if (longitude < -PI)
498  longitude = -PI;
499 
500  return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
501 }
502 
503 
504 
505 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:01 for GeographicTranslator by doxygen 1.8.2