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