UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
Stereographic.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: STEREOGRAPHIC
5  *
6  *
7  * ABSTRACT
8  *
9  * This component provides conversions between geodetic (latitude and
10  * longitude) coordinates and Stereographic (easting and northing)
11  * coordinates.
12  *
13  * ERROR HANDLING
14  *
15  * This component checks parameters for valid values. If an invalid
16  * value is found the error code is combined with the current error code
17  * using the bitwise or. This combining allows multiple error codes to
18  * be returned. The possible error codes are:
19  *
20  * STEREO_NO_ERROR : No errors occurred in function
21  * STEREO_LAT_ERROR : Latitude outside of valid range
22  * (-90 to 90 degrees)
23  * STEREO_LON_ERROR : Longitude outside of valid range
24  * (-180 to 360 degrees)
25  * STEREO_ORIGIN_LAT_ERROR : Origin latitude outside of valid range
26  * (-90 to 90 degrees)
27  * STEREO_CENT_MER_ERROR : Central meridian outside of valid range
28  * (-180 to 360 degrees)
29  * STEREO_EASTING_ERROR : Easting outside of valid range,
30  * (False_Easting +/- ~1,460,090,226 m,
31  * depending on ellipsoid and projection
32  * parameters)
33  * STEREO_NORTHING_ERROR : Northing outside of valid range,
34  * (False_Northing +/- ~1,460,090,226 m,
35  * depending on ellipsoid and projection
36  * parameters)
37  * STEREO_A_ERROR : Semi-major axis less than or equal to zero
38  * STEREO_INV_F_ERROR : Inverse flattening outside of valid range
39  * (250 to 350)
40  *
41  *
42  * REUSE NOTES
43  *
44  * STEREOGRAPHIC is intended for reuse by any application that
45  * performs a Stereographic projection.
46  *
47  *
48  * REFERENCES
49  *
50  * Further information on STEREOGRAPHIC can be found in the
51  * Reuse Manual.
52  *
53  *
54  * STEREOGRAPHIC originated from :
55  * U.S. Army Topographic Engineering Center
56  * Geospatial Information Division
57  * 7701 Telegraph Road
58  * Alexandria, VA 22310-3864
59  *
60  *
61  * LICENSES
62  *
63  * None apply to this component.
64  *
65  *
66  * RESTRICTIONS
67  *
68  * STEREOGRAPHIC has no restrictions.
69  *
70  *
71  * ENVIRONMENT
72  *
73  * STEREOGRAPHIC was tested and certified in the following
74  * environments:
75  *
76  * 1. Solaris 2.5 with GCC, version 2.8.1
77  * 2. Window 95 with MS Visual C++, version 6
78  *
79  *
80  * MODIFICATIONS
81  *
82  * Date Description
83  * ---- -----------
84  * 3-1-07 Original Code
85  *
86  */
87 
88 
89 /***************************************************************************/
90 /*
91  * INCLUDES
92  */
93 
94 #include <math.h>
95 #include "Stereographic.h"
96 #include "MapProjection4Parameters.h"
97 #include "MapProjectionCoordinates.h"
98 #include "GeodeticCoordinates.h"
99 #include "CoordinateConversionException.h"
100 #include "ErrorMessages.h"
101 
102 /*
103  * math.h - Standard C math library
104  * Stereographic.h - Is for prototype error checking
105  * MapProjectionCoordinates.h - defines map projection coordinates
106  * GeodeticCoordinates.h - defines geodetic coordinates
107  * CoordinateConversionException.h - Exception handler
108  * ErrorMessages.h - Contains exception messages
109  */
110 
111 
112 using namespace MSP::CCS;
113 
114 
115 /***************************************************************************/
116 /* DEFINES
117  *
118  */
119 
120 const double PI = 3.14159265358979323e0; /* PI */
121 const double PI_OVER_2 = (PI / 2.0);
122 const double PI_OVER_4 = (PI / 4.0);
123 const double TWO_PI = (2.0 * PI);
124 const double ONE = (1.0 * PI / 180.0); /* One degree */
125 
126 
127 /************************************************************************/
128 /* FUNCTIONS
129  *
130  */
131 
132 Stereographic::Stereographic( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double originLatitude, double falseEasting, double falseNorthing ) :
133  Stereo_Ra( 6371007.1810824 ),
134  Two_Stereo_Ra( 12742014.3621648 ),
135  Stereo_At_Pole( 0 ),
136  Stereo_Origin_Long( 0.0 ),
137  Stereo_Origin_Lat( 0.0 ),
138  Stereo_False_Easting( 0.0 ),
139  Stereo_False_Northing( 0.0 ),
140  Sin_Stereo_Origin_Lat( 0.0 ),
141  Cos_Stereo_Origin_Lat( 1.0 ),
142  Stereo_Delta_Easting( 1460090226.0 ),
143  Stereo_Delta_Northing( 1460090226.0 )
144 {
145 /*
146  * The constructor receives the ellipsoid
147  * parameters and Stereograpic projection parameters as inputs, and
148  * sets the corresponding state variables. If any errors occur, an
149  * exception is thrown with a description of the error.
150  *
151  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
152  * ellipsoidFlattening : Flattening of ellipsoid (input)
153  * centralMeridian : Longitude, in radians, at the center of (input)
154  * the projection
155  * originLatitude : Latitude, in radians, at the center of (input)
156  * the projection
157  * falseEasting : Easting (X) at center of projection, in meters (input)
158  * falseNorthing : Northing (Y) at center of projection, in meters (input)
159  */
160 
161  double es2, es4, es6;
162  double temp = 0;
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  Stereo_Ra = semiMajorAxis *
189  (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 /3024.0);
190  Two_Stereo_Ra = 2.0 * Stereo_Ra;
191  Stereo_Origin_Lat = originLatitude;
192  Sin_Stereo_Origin_Lat = sin(Stereo_Origin_Lat);
193  Cos_Stereo_Origin_Lat = cos(Stereo_Origin_Lat);
194  if (centralMeridian > PI)
195  centralMeridian -= TWO_PI;
196  Stereo_Origin_Long = centralMeridian;
197  Stereo_False_Easting = falseEasting;
198  Stereo_False_Northing = falseNorthing;
199  if(fabs(fabs(Stereo_Origin_Lat) - PI_OVER_2) < 1.0e-10)
200  Stereo_At_Pole = 1;
201  else
202  Stereo_At_Pole = 0;
203 
204  if ((Stereo_At_Pole) || (fabs(Stereo_Origin_Lat) < 1.0e-10))
205  {
206  Stereo_Delta_Easting = 1460090226.0;
207  }
208  else
209  {
210  MapProjectionCoordinates* tempCoordinates;
211  if (Stereo_Origin_Long <= 0)
212  {
213  GeodeticCoordinates gcTemp( CoordinateType::geodetic, PI + Stereo_Origin_Long - ONE, -Stereo_Origin_Lat );
214  tempCoordinates = convertFromGeodetic( &gcTemp );
215  }
216  else
217  {
218  GeodeticCoordinates gcTemp( CoordinateType::geodetic, Stereo_Origin_Long - PI - ONE, -Stereo_Origin_Lat );
219  tempCoordinates = convertFromGeodetic( &gcTemp );
220  }
221  Stereo_Delta_Easting = tempCoordinates->easting();
222  delete tempCoordinates;
223 
224  if(Stereo_False_Easting)
225  Stereo_Delta_Easting -= Stereo_False_Easting;
226  if (Stereo_Delta_Easting < 0)
227  Stereo_Delta_Easting = -Stereo_Delta_Easting;
228  }
229 }
230 
231 
232 Stereographic::Stereographic( const Stereographic &s )
233 {
234  semiMajorAxis = s.semiMajorAxis;
235  flattening = s.flattening;
236  Stereo_Ra = s.Stereo_Ra;
237  Two_Stereo_Ra = s.Two_Stereo_Ra;
238  Stereo_At_Pole = s.Stereo_At_Pole;
239  Stereo_Origin_Long = s.Stereo_Origin_Long;
240  Stereo_Origin_Lat = s.Stereo_Origin_Lat;
241  Stereo_False_Easting = s.Stereo_False_Easting;
242  Stereo_False_Northing = s.Stereo_False_Northing;
243  Sin_Stereo_Origin_Lat = s.Sin_Stereo_Origin_Lat;
244  Cos_Stereo_Origin_Lat = s.Cos_Stereo_Origin_Lat;
245  Stereo_Delta_Easting = s.Stereo_Delta_Easting;
246  Stereo_Delta_Northing = s.Stereo_Delta_Northing;
247 }
248 
249 
250 Stereographic::~Stereographic()
251 {
252 }
253 
254 
255 Stereographic& Stereographic::operator=( const Stereographic &s )
256 {
257  if( this != &s )
258  {
259  semiMajorAxis = s.semiMajorAxis;
260  flattening = s.flattening;
261  Stereo_Ra = s.Stereo_Ra;
262  Two_Stereo_Ra = s.Two_Stereo_Ra;
263  Stereo_At_Pole = s.Stereo_At_Pole;
264  Stereo_Origin_Long = s.Stereo_Origin_Long;
265  Stereo_Origin_Lat = s.Stereo_Origin_Lat;
266  Stereo_False_Easting = s.Stereo_False_Easting;
267  Stereo_False_Northing = s.Stereo_False_Northing;
268  Sin_Stereo_Origin_Lat = s.Sin_Stereo_Origin_Lat;
269  Cos_Stereo_Origin_Lat = s.Cos_Stereo_Origin_Lat;
270  Stereo_Delta_Easting = s.Stereo_Delta_Easting;
271  Stereo_Delta_Northing = s.Stereo_Delta_Northing;
272  }
273 
274  return *this;
275 }
276 
277 
278 MapProjection4Parameters* Stereographic::getParameters() const
279 {
280 /*
281  * The function getParameters returns the current ellipsoid
282  * parameters and Stereographic projection parameters.
283  *
284  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
285  * ellipsoidFlattening : Flattening of ellipsoid (output)
286  * centralMeridian : Longitude, in radians, at the center of (output)
287  * the projection
288  * originLatitude : Latitude, in radians, at the center of (output)
289  * the projection
290  * falseEasting : A coordinate value, in meters, assigned to the
291  * central meridian of the projection. (output)
292  * falseNorthing : A coordinate value, in meters, assigned to the
293  * origin latitude of the projection (output)
294  */
295 
296  return new MapProjection4Parameters( CoordinateType::stereographic, Stereo_Origin_Long, Stereo_Origin_Lat, Stereo_False_Easting, Stereo_False_Northing );
297 }
298 
299 
300 MSP::CCS::MapProjectionCoordinates* Stereographic::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
301 {
302 /*
303  * The function convertFromGeodetic converts geodetic
304  * coordinates (latitude and longitude) to Stereographic coordinates
305  * (easting and northing), according to the current ellipsoid
306  * and Stereographic projection parameters. If any errors occur,
307  * an exception is thrown with a description of the error.
308  *
309  * longitude : Longitude, in radians (input)
310  * latitude : Latitude, in radians (input)
311  * easting : Easting (X), in meters (output)
312  * northing : Northing (Y), in meters (output)
313  */
314 
315  double g, k;
316  double num = 0;
317  double Ra_k = 0;
318  double dlam; /* Longitude - Central Meridan */
319  double cos_dlam;
320  double easting, northing;
321 
322  double longitude = geodeticCoordinates->longitude();
323  double latitude = geodeticCoordinates->latitude();
324  double slat = sin(latitude);
325  double clat = cos(latitude);
326 
327  if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
328  { /* Latitude out of range */
329  throw CoordinateConversionException( ErrorMessages::latitude );
330  }
331  if ((longitude < -PI) || (longitude > TWO_PI))
332  { /* Longitude out of range */
333  throw CoordinateConversionException( ErrorMessages::longitude );
334  }
335 
336  dlam = longitude - Stereo_Origin_Long;
337  if (dlam > PI)
338  {
339  dlam -= TWO_PI;
340  }
341  if (dlam < -PI)
342  {
343  dlam += TWO_PI;
344  }
345 
346  cos_dlam = cos(dlam);
347  g = 1.0 + Sin_Stereo_Origin_Lat * slat + Cos_Stereo_Origin_Lat * clat * cos_dlam;
348  if (fabs(g) <= 1.0e-10)
349  { /* Point is out of view. Will return longitude out of range message
350  since no point out of view is implemented. */
351  throw CoordinateConversionException( ErrorMessages::longitude );
352  }
353  else
354  {
355  if (Stereo_At_Pole)
356  {
357  if (fabs(fabs(latitude) - PI_OVER_2) < 1.0e-10)
358  {
359  easting = Stereo_False_Easting;
360  northing = Stereo_False_Northing;
361  }
362  else
363  {
364  if (Stereo_Origin_Lat > 0)
365  {
366  num = Two_Stereo_Ra * tan(PI_OVER_4 - latitude / 2.0);
367  easting = Stereo_False_Easting + num * sin(dlam);
368  northing = Stereo_False_Northing + (-num * cos_dlam);
369  }
370  else
371  {
372  num = Two_Stereo_Ra * tan(PI_OVER_4 + latitude / 2.0);
373  easting = Stereo_False_Easting + num * sin(dlam);
374  northing = Stereo_False_Northing + num * cos_dlam;
375  }
376  }
377  }
378  else
379  {
380  if (fabs(Stereo_Origin_Lat) <= 1.0e-10)
381  {
382  k = 2.0 / (1.0 + clat * cos_dlam);
383  Ra_k = Stereo_Ra * k;
384  northing = Stereo_False_Northing + Ra_k * slat;
385  }
386  else
387  {
388  k = 2.0 / g;
389  Ra_k = Stereo_Ra * k;
390  northing = Stereo_False_Northing + Ra_k * (Cos_Stereo_Origin_Lat * slat - Sin_Stereo_Origin_Lat * clat * cos_dlam);
391  }
392  easting = Stereo_False_Easting + Ra_k * clat * sin(dlam);
393  }
394  }
395 
396  return new MapProjectionCoordinates( CoordinateType::stereographic, easting, northing );
397 }
398 
399 
400 MSP::CCS::GeodeticCoordinates* Stereographic::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
401 {
402 /*
403  * The function convertToGeodetic converts Stereographic projection
404  * (easting and northing) coordinates to geodetic (latitude and longitude)
405  * coordinates, according to the current ellipsoid and Stereographic projection
406  * coordinates. If any errors occur, an exception is thrown with a description
407  * of the error.
408  *
409  * easting : Easting (X), in meters (input)
410  * northing : Northing (Y), in meters (input)
411  * longitude : Longitude (lambda), in radians (output)
412  * latitude : Latitude (phi), in radians (output)
413  */
414 
415  double dx, dy;
416  double rho, c;
417  double sin_c, cos_c;
418  double dy_sin_c;
419  double longitude, latitude;
420 
421  double easting = mapProjectionCoordinates->easting();
422  double northing = mapProjectionCoordinates->northing();
423 
424  if ((easting < (Stereo_False_Easting - Stereo_Delta_Easting))
425  ||(easting > (Stereo_False_Easting + Stereo_Delta_Easting)))
426  { /* Easting out of range */
427  throw CoordinateConversionException( ErrorMessages::easting );
428  }
429  if ((northing < (Stereo_False_Northing - Stereo_Delta_Northing))
430  || (northing > (Stereo_False_Northing + Stereo_Delta_Northing)))
431  { /* Northing out of range */
432  throw CoordinateConversionException( ErrorMessages::northing );
433  }
434 
435  dy = northing - Stereo_False_Northing;
436  dx = easting - Stereo_False_Easting;
437  rho = sqrt(dx * dx + dy * dy);
438  if (fabs(rho) <= 1.0e-10)
439  {
440  latitude = Stereo_Origin_Lat;
441  longitude = Stereo_Origin_Long;
442  }
443  else
444  {
445  c = 2.0 * atan(rho / (Two_Stereo_Ra));
446  sin_c = sin(c);
447  cos_c = cos(c);
448  dy_sin_c = dy * sin_c;
449  if (Stereo_At_Pole)
450  {
451  if (Stereo_Origin_Lat > 0)
452  longitude = Stereo_Origin_Long + atan2(dx, -dy);
453  else
454  longitude = Stereo_Origin_Long + atan2(dx, dy);
455  }
456  else
457  longitude = Stereo_Origin_Long + atan2(dx * sin_c, (rho * Cos_Stereo_Origin_Lat * cos_c - dy_sin_c * Sin_Stereo_Origin_Lat));
458  latitude = asin(cos_c * Sin_Stereo_Origin_Lat + ((dy_sin_c * Cos_Stereo_Origin_Lat) / rho));
459  }
460 
461  if (fabs(latitude) < 2.2e-8) /* force lat to 0 to avoid -0 degrees */
462  latitude = 0.0;
463  if (latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
464  latitude = PI_OVER_2;
465  else if (latitude < -PI_OVER_2)
466  latitude = -PI_OVER_2;
467 
468  if (longitude > PI)
469  {
470  if (longitude - PI < 3.5e-6)
471  longitude = PI;
472  else
473  longitude -= TWO_PI;
474  }
475  if (longitude < -PI)
476  {
477  if (fabs(longitude + PI) < 3.5e-6)
478  longitude = -PI;
479  else
480  longitude += TWO_PI;
481  }
482 
483  if (fabs(longitude) < 2.0e-7) /* force lon to 0 to avoid -0 degrees */
484  longitude = 0.0;
485  if (longitude > PI) /* force distorted values to 180, -180 degrees */
486  longitude = PI;
487  else if (longitude < -PI)
488  longitude = -PI;
489 
490  return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
491 }
492 
493 
494 
495 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:02 for GeographicTranslator by doxygen 1.8.2