UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
Gnomonic.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: GNOMONIC
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates
9  * (latitude and longitude in radians) and Gnomonic
10  * projection coordinates (easting and northing in meters). This projection
11  * employs a spherical Earth model. The spherical radius used is the radius
12  * of 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  * GNOM_NO_ERROR : No errors occurred in function
22  * GNOM_LAT_ERROR : Latitude outside of valid range
23  * (-90 to 90 degrees)
24  * GNOM_LON_ERROR : Longitude outside of valid range
25  * (-180 to 360 degrees)
26  * GNOM_EASTING_ERROR : Easting outside of valid range
27  * (depends on ellipsoid and projection
28  * parameters)
29  * GNOM_NORTHING_ERROR : Northing outside of valid range
30  * (depends on ellipsoid and projection
31  * parameters)
32  * GNOM_ORIGIN_LAT_ERROR : Origin latitude outside of valid range
33  * (-90 to 90 degrees)
34  * GNOM_CENT_MER_ERROR : Central meridian outside of valid range
35  * (-180 to 360 degrees)
36  * GNOM_A_ERROR : Semi-major axis less than or equal to zero
37  * GNOM_INV_F_ERROR : Inverse flattening outside of valid range
38  * (250 to 350)
39  *
40  *
41  * REUSE NOTES
42  *
43  * GNOMONIC is intended for reuse by any application that
44  * performs a Gnomonic projection or its inverse.
45  *
46  * REFERENCES
47  *
48  * Further information on GNOMONIC can be found in the Reuse Manual.
49  *
50  * GNOMONIC originated from: 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  * GNOMONIC has no restrictions.
62  *
63  * ENVIRONMENT
64  *
65  * GNOMONIC was tested and certified in the following environments:
66  *
67  * 1. Solaris 2.5 with GCC, version 2.8.1
68  * 2. MSDOS with MS Visual C++, version 6
69  *
70  * MODIFICATIONS
71  *
72  * Date Description
73  * ---- -----------
74  * 05-22-00 Original Code
75  * 03-05-07 Original C++ Code
76  * 05-31-11 Fixed dicsontinuitty at south pole. DR 27958
77  *
78  *
79  */
80 
81 
82 /***************************************************************************/
83 /*
84  * INCLUDES
85  */
86 
87 #include <math.h>
88 #include "Gnomonic.h"
89 #include "MapProjection4Parameters.h"
90 #include "MapProjectionCoordinates.h"
91 #include "GeodeticCoordinates.h"
92 #include "CoordinateConversionException.h"
93 #include "ErrorMessages.h"
94 
95 /*
96  * math.h - Standard C math library
97  * Gnomonic.h - Is for prototype error checking
98  * MapProjectionCoordinates.h - defines map projection coordinates
99  * GeodeticCoordinates.h - defines geodetic coordinates
100  * CoordinateConversionException.h - Exception handler
101  * ErrorMessages.h - Contains exception messages
102  */
103 
104 using namespace MSP::CCS;
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 
115 
116 /************************************************************************/
117 /* FUNCTIONS
118  *
119  */
120 
121 Gnomonic::Gnomonic(
122  double ellipsoidSemiMajorAxis,
123  double ellipsoidFlattening,
124  double centralMeridian,
125  double originLatitude,
126  double falseEasting,
127  double falseNorthing ) :
128  CoordinateSystem(),
129  Ra( 6371007.1810824 ),
130  Sin_Gnom_Origin_Lat( 0.0 ),
131  Cos_Gnom_Origin_Lat( 1.0 ),
132  Gnom_Origin_Long( 0.0 ),
133  Gnom_Origin_Lat( 0.0 ),
134  Gnom_False_Easting( 0.0 ),
135  Gnom_False_Northing( 0.0 ),
136  abs_Gnom_Origin_Lat( 0.0 ),
137  Gnom_Delta_Northing( 40000000 ),
138  Gnom_Delta_Easting( 40000000 )
139 {
140 /*
141  * The constructor receives the ellipsoid parameters and
142  * projection parameters as inputs, and sets the corresponding state
143  * variables. If any errors occur, an exception is thrown with a description
144  * of the error.
145  *
146  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
147  * ellipsoidFlattening : Flattening of ellipsoid (input)
148  * centralMeridian : Longitude in radians at the center of (input)
149  * the projection
150  * originLatitude : Latitude in radians at which the (input)
151  * point scale factor is 1.0
152  * falseEasting : A coordinate value in meters assigned to the
153  * central meridian of the projection. (input)
154  * falseNorthing : A coordinate value in meters assigned to the
155  * origin latitude of the projection (input)
156  */
157 
158  double es2, es4, es6;
159  double inv_f = 1 / ellipsoidFlattening;
160 
161  if (ellipsoidSemiMajorAxis <= 0.0)
162  { /* Semi-major axis must be greater than zero */
163  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
164  }
165  if ((inv_f < 250) || (inv_f > 350))
166  { /* Inverse flattening must be between 250 and 350 */
167  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
168  }
169  if ((originLatitude < -PI_OVER_2) || (originLatitude > PI_OVER_2))
170  { /* origin latitude out of range */
171  throw CoordinateConversionException( ErrorMessages::originLatitude );
172  }
173  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
174  { /* origin longitude out of range */
175  throw CoordinateConversionException( ErrorMessages::centralMeridian );
176  }
177 
178  semiMajorAxis = ellipsoidSemiMajorAxis;
179  flattening = ellipsoidFlattening;
180 
181  es2 = 2 * flattening - flattening * flattening;
182  es4 = es2 * es2;
183  es6 = es4 * es2;
184  /* spherical radius */
185  Ra = semiMajorAxis *
186  (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 / 3024.0);
187  Gnom_Origin_Lat = originLatitude;
188  Sin_Gnom_Origin_Lat = sin(Gnom_Origin_Lat);
189  Cos_Gnom_Origin_Lat = cos(Gnom_Origin_Lat);
190  abs_Gnom_Origin_Lat = fabs(Gnom_Origin_Lat);
191  if (centralMeridian > PI)
192  centralMeridian -= TWO_PI;
193  Gnom_Origin_Long = centralMeridian;
194  Gnom_False_Northing = falseNorthing;
195  Gnom_False_Easting = falseEasting;
196 }
197 
198 
199 Gnomonic::Gnomonic( const Gnomonic &g )
200 {
201  semiMajorAxis = g.semiMajorAxis;
202  flattening = g.flattening;
203  Ra = g.Ra;
204  Sin_Gnom_Origin_Lat = g.Sin_Gnom_Origin_Lat;
205  Cos_Gnom_Origin_Lat = g.Cos_Gnom_Origin_Lat;
206  Gnom_Origin_Long = g.Gnom_Origin_Long;
207  Gnom_Origin_Lat = g.Gnom_Origin_Lat;
208  Gnom_False_Easting = g.Gnom_False_Easting;
209  Gnom_False_Northing = g.Gnom_False_Northing;
210  abs_Gnom_Origin_Lat = g.abs_Gnom_Origin_Lat;
211  Gnom_Delta_Northing = g.Gnom_Delta_Northing;
212  Gnom_Delta_Easting = g.Gnom_Delta_Easting;
213 }
214 
215 
216 Gnomonic::~Gnomonic()
217 {
218 }
219 
220 
221 Gnomonic& Gnomonic::operator=( const Gnomonic &g )
222 {
223  if( this != &g )
224  {
225  semiMajorAxis = g.semiMajorAxis;
226  flattening = g.flattening;
227  Ra = g.Ra;
228  Sin_Gnom_Origin_Lat = g.Sin_Gnom_Origin_Lat;
229  Cos_Gnom_Origin_Lat = g.Cos_Gnom_Origin_Lat;
230  Gnom_Origin_Long = g.Gnom_Origin_Long;
231  Gnom_Origin_Lat = g.Gnom_Origin_Lat;
232  Gnom_False_Easting = g.Gnom_False_Easting;
233  Gnom_False_Northing = g.Gnom_False_Northing;
234  abs_Gnom_Origin_Lat = g.abs_Gnom_Origin_Lat;
235  Gnom_Delta_Northing = g.Gnom_Delta_Northing;
236  Gnom_Delta_Easting = g.Gnom_Delta_Easting;
237  }
238 
239  return *this;
240 }
241 
242 
243 MapProjection4Parameters* Gnomonic::getParameters() const
244 {
245 /*
246  * The function getParameters returns the current ellipsoid
247  * parameters and Gnomonic projection parameters.
248  *
249  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
250  * ellipsoidFlattening : Flattening of ellipsoid (output)
251  * centralMeridian : Longitude in radians at the center of (output)
252  * the projection
253  * originLatitude : Latitude in radians at which the (output)
254  * point scale factor is 1.0
255  * falseEasting : A coordinate value in meters assigned to the
256  * central meridian of the projection. (output)
257  * falseNorthing : A coordinate value in meters assigned to the
258  * origin latitude of the projection (output)
259  */
260 
261  return new MapProjection4Parameters(
262  CoordinateType::gnomonic, Gnom_Origin_Long, Gnom_Origin_Lat,
263  Gnom_False_Easting, Gnom_False_Northing );
264 }
265 
266 
267 MSP::CCS::MapProjectionCoordinates* Gnomonic::convertFromGeodetic(
268  MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
269 {
270 /*
271  * The function convertFromGeodetic converts geodetic (latitude and
272  * longitude) coordinates to Gnomonic projection (easting and northing)
273  * coordinates, according to the current ellipsoid and Gnomonic projection
274  * parameters. If any errors occur, an exception is thrown with a description
275  * of the error.
276  *
277  * longitude : Longitude (lambda) in radians (input)
278  * latitude : Latitude (phi) in radians (input)
279  * easting : Easting (X) in meters (output)
280  * northing : Northing (Y) in meters (output)
281  */
282 
283  double dlam; /* Longitude - Central Meridan */
284  double cos_c;
285  double k_prime; /* scale factor */
286  double Ra_kprime;
287  double Ra_cotlat;
288  double sin_dlam, cos_dlam;
289  double temp_Easting, temp_Northing;
290  double easting, northing;
291 
292  double longitude = geodeticCoordinates->longitude();
293  double latitude = geodeticCoordinates->latitude();
294  double slat = sin(latitude);
295  double clat = cos(latitude);
296 
297  if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
298  { /* Latitude out of range */
299  throw CoordinateConversionException( ErrorMessages::latitude );
300  }
301  if ((longitude < -PI) || (longitude > TWO_PI))
302  { /* Longitude out of range */
303  throw CoordinateConversionException( ErrorMessages::longitude );
304  }
305  dlam = longitude - Gnom_Origin_Long;
306  sin_dlam = sin(dlam);
307  cos_dlam = cos(dlam);
308  cos_c = Sin_Gnom_Origin_Lat * slat + Cos_Gnom_Origin_Lat * clat * cos_dlam;
309  if (cos_c <= 1.0e-10)
310  { /* Point is out of view. Will return longitude out of range message
311  since no point out of view is implemented. */
312  throw CoordinateConversionException( ErrorMessages::longitude );
313  }
314 
315  if (dlam > PI)
316  {
317  dlam -= TWO_PI;
318  }
319  if (dlam < -PI)
320  {
321  dlam += TWO_PI;
322  }
323 
324  if (fabs(abs_Gnom_Origin_Lat - PI_OVER_2) < 1.0e-10)
325  {
326  Ra_cotlat = Ra * (clat / slat);
327  temp_Easting = Ra_cotlat * sin_dlam;
328  temp_Northing = Ra_cotlat * cos_dlam;
329  if (Gnom_Origin_Lat >= 0.0)
330  {
331  easting = temp_Easting + Gnom_False_Easting;
332  northing = -1.0 * temp_Northing + Gnom_False_Northing;
333  }
334  else
335  {
336  easting = -1.0 * temp_Easting + Gnom_False_Easting;
337  northing = -1.0 * temp_Northing + Gnom_False_Northing;
338  }
339  }
340  else if (abs_Gnom_Origin_Lat <= 1.0e-10)
341  {
342  easting = Ra * tan(dlam) + Gnom_False_Easting;
343  northing = Ra * tan(latitude) / cos_dlam + Gnom_False_Northing;
344  }
345  else
346  {
347  k_prime = 1 / cos_c;
348  Ra_kprime = Ra * k_prime;
349  easting = Ra_kprime * clat * sin_dlam + Gnom_False_Easting;
350  northing = Ra_kprime *
351  (Cos_Gnom_Origin_Lat * slat - Sin_Gnom_Origin_Lat * clat * cos_dlam)
352  + Gnom_False_Northing;
353  }
354 
355  return new MapProjectionCoordinates(
356  CoordinateType::gnomonic, easting, northing );
357 }
358 
359 
360 MSP::CCS::GeodeticCoordinates* Gnomonic::convertToGeodetic(
361  MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
362 {
363 /*
364  * The function convertToGeodetic converts Gnomonic projection
365  * (easting and northing) coordinates to geodetic (latitude and longitude)
366  * coordinates, according to the current ellipsoid and Gnomonic projection
367  * coordinates. If any errors occur, an exception is thrown with a description
368  * of the error.
369  *
370  * easting : Easting (X) in meters (input)
371  * northing : Northing (Y) in meters (input)
372  * longitude : Longitude (lambda) in radians (output)
373  * latitude : Latitude (phi) in radians (output)
374  */
375 
376  double dx, dy;
377  double rho;
378  double c;
379  double sin_c, cos_c;
380  double dy_sinc;
381  double longitude, latitude;
382 
383  double easting = mapProjectionCoordinates->easting();
384  double northing = mapProjectionCoordinates->northing();
385 
386  if ((easting < (Gnom_False_Easting - Gnom_Delta_Easting))
387  || (easting > (Gnom_False_Easting + Gnom_Delta_Easting)))
388  { /* Easting out of range */
389  throw CoordinateConversionException( ErrorMessages::easting );
390  }
391  if ((northing < (Gnom_False_Northing - Gnom_Delta_Northing))
392  || (northing > (Gnom_False_Northing + Gnom_Delta_Northing)))
393  { /* Northing out of range */
394  throw CoordinateConversionException( ErrorMessages::northing );
395  }
396 
397  dy = northing - Gnom_False_Northing;
398  dx = easting - Gnom_False_Easting;
399  rho = sqrt(dx * dx + dy * dy);
400  if (fabs(rho) <= 1.0e-10)
401  {
402  latitude = Gnom_Origin_Lat;
403  longitude = Gnom_Origin_Long;
404  }
405  else
406  {
407  c = atan(rho / Ra);
408  sin_c = sin(c);
409  cos_c = cos(c);
410  dy_sinc = dy * sin_c;
411  latitude = asin((cos_c * Sin_Gnom_Origin_Lat)
412  + ((dy_sinc * Cos_Gnom_Origin_Lat) / rho));
413 
414  if (fabs(abs_Gnom_Origin_Lat - PI_OVER_2) < 1.0e-10)
415  {
416  if (Gnom_Origin_Lat >= 0.0)
417  longitude = Gnom_Origin_Long + atan2(dx, -dy);
418  else
419  longitude = Gnom_Origin_Long + atan2(dx, dy);
420  }
421  else
422  longitude = Gnom_Origin_Long
423  + atan2((dx * sin_c),
424  (rho * Cos_Gnom_Origin_Lat * cos_c - dy_sinc *Sin_Gnom_Origin_Lat));
425  }
426  if (latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
427  latitude = PI_OVER_2;
428  else if (latitude < -PI_OVER_2)
429  latitude = -PI_OVER_2;
430  if (longitude > PI)
431  longitude -= TWO_PI;
432  if (longitude < -PI)
433  longitude += TWO_PI;
434  if (longitude > PI) /* force distorted values to 180, -180 degrees */
435  longitude = PI;
436  else if (longitude < -PI)
437  longitude = -PI;
438 
439  return new GeodeticCoordinates(
440  CoordinateType::geodetic, longitude, latitude );
441 }
442 
443 
444 
445 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:02 for GeographicTranslator by doxygen 1.8.2