UNCLASSIFIED

GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
LocalCartesian.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: LOCAL CARTESIAN
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates (latitude,
9  * longitude in radians and height in meters) and Local Cartesian coordinates
10  * (X, Y, Z).
11  *
12  * ERROR HANDLING
13  *
14  * This component checks parameters for valid values. If an invalid value
15  * is found, the error code is combined with the current error code using
16  * the bitwise or. This combining allows multiple error codes to be
17  * returned. The possible error codes are:
18  *
19  * LOCCART_NO_ERROR : No errors occurred in function
20  * LOCCART_LAT_ERROR : Latitude out of valid range
21  * (-90 to 90 degrees)
22  * LOCCART_LON_ERROR : Longitude out of valid range
23  * (-180 to 360 degrees)
24  * LOCCART_A_ERROR : Semi-major axis less than or equal to zero
25  * LOCCART_INV_F_ERROR : Inverse flattening outside of valid range
26  * (250 to 350)
27  * LOCCART_ORIGIN_LAT_ERROR : Origin Latitude out of valid range
28  * (-90 to 90 degrees)
29  * LOCCART_ORIGIN_LON_ERROR : Origin Longitude out of valid range
30  * (-180 to 360 degrees)
31  * LOCCART_ORIENTATION_ERROR : Orientation angle out of valid range
32  * (-360 to 360 degrees)
33  *
34  *
35  * REUSE NOTES
36  *
37  * LOCCART is intended for reuse by any application that performs
38  * coordinate conversions between geodetic coordinates or geocentric
39  * coordinates and local cartesian coordinates..
40  *
41  *
42  * REFERENCES
43  *
44  * Further information on LOCAL CARTESIAN can be found in the Reuse Manual.
45  *
46  * LOCCART originated from : U.S. Army Topographic Engineering Center
47  * Geospatial Inforamtion Division
48  * 7701 Telegraph Road
49  * Alexandria, VA 22310-3864
50  *
51  * LICENSES
52  *
53  * None apply to this component.
54  *
55  * RESTRICTIONS
56  *
57  * LOCCART has no restrictions.
58  *
59  * ENVIRONMENT
60  *
61  * LOCCART was tested and certified in the following environments:
62  *
63  * 1. Solaris 2.5 with GCC version 2.8.1
64  * 2. Windows 95 with MS Visual C++ version 6
65  *
66  * MODIFICATIONS
67  *
68  * Date Description
69  * ---- -----------
70  * 07-16-99 Original Code
71  * 03-2-07 Original C++ Code
72  * 3/23/11 N. Lundgren BAEts28583 Updated for memory leaks in
73  * convertFromGeodetic and convertToGeodetic
74  *
75  */
76 
77 
78 /***************************************************************************/
79 /*
80  * INCLUDES
81  */
82 
83 #include <math.h>
84 #include "Geocentric.h"
85 #include "LocalCartesian.h"
87 #include "CartesianCoordinates.h"
88 #include "GeodeticCoordinates.h"
90 #include "ErrorMessages.h"
91 
92 /*
93  * math.h - Standard C math library
94  * Geocentric.h - Is needed to call the convertToGeocentric and
95  * convertFromGeocentric functions
96  * LocalCartesian.h - Is for prototype error checking.
97  * CartesianCoordinates.h - defines cartesian coordinates
98  * GeodeticCoordinates.h - defines geodetic coordinates
99  * CoordinateConversionException.h - Exception handler
100  * ErrorMessages.h - Contains exception messages
101  */
102 
103 
104 using namespace MSP::CCS;
105 
106 
107 /***************************************************************************/
108 /* DEFINES
109  *
110  */
111 
112 const double PI = 3.14159265358979323e0; /* PI */
113 const double PI_OVER_2 = ( PI / 2.0e0);
114 const double TWO_PI = (2.0 * PI);
115 
116 
117 /************************************************************************/
118 /* FUNCTIONS
119  *
120  */
121 
122 LocalCartesian::LocalCartesian( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double originLongitude, double originLatitude, double originHeight, double orientation ) :
124  geocentric( 0 ),
125  es2( 0.0066943799901413800 ),
126  u0( 6378137.0 ),
127  v0( 0.0 ),
128  w0( 0.0 ),
129  LocalCart_Origin_Long( 0.0 ),
130  LocalCart_Origin_Lat( 0.0 ),
131  LocalCart_Origin_Height( 0.0 ),
132  LocalCart_Orientation( 0.0 ),
133  Sin_LocalCart_Origin_Lat( 0.0 ),
134  Cos_LocalCart_Origin_Lat( 1.0 ),
135  Sin_LocalCart_Origin_Lon( 0.0 ),
136  Cos_LocalCart_Origin_Lon( 1.0 ),
137  Sin_LocalCart_Orientation( 0.0 ),
138  Cos_LocalCart_Orientation( 1.0 ),
139  Sin_Lat_Sin_Orient( 0.0 ),
140  Sin_Lat_Cos_Orient( 0.0 )
141 {
142 /*
143  * The constructor receives the ellipsoid parameters
144  * and local origin parameters as inputs and sets the corresponding state variables.
145  *
146  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
147  * ellipsoidFlattening : Flattening of ellipsoid (input)
148  * originLongitude : Longitude of the local origin, in radians (input)
149  * originLatitude : Latitude of the local origin, in radians (input)
150  * originHeight : Ellipsoid height of the local origin, in meters (input)
151  * orientation : Orientation angle of the local cartesian coordinate system,
152  * in radians (input)
153  */
154 
155  double N0;
156  double inv_f = 1 / ellipsoidFlattening;
157  double val;
158 
159  if (ellipsoidSemiMajorAxis <= 0.0)
160  { /* Semi-major axis must be greater than zero */
162  }
163  if ((inv_f < 250) || (inv_f > 350))
164  { /* Inverse flattening must be between 250 and 350 */
166  }
167  if ((originLatitude < -PI_OVER_2) || (originLatitude > PI_OVER_2))
168  { /* origin latitude out of range */
170  }
171  if ((originLongitude < -PI) || (originLongitude > TWO_PI))
172  { /* origin longitude out of range */
174  }
175  if ((orientation < -PI) || (orientation > TWO_PI))
176  { /* orientation angle out of range */
178  }
179 
180  semiMajorAxis = ellipsoidSemiMajorAxis;
181  flattening = ellipsoidFlattening;
182 
183  LocalCart_Origin_Lat = originLatitude;
184  if (originLongitude > PI)
185  originLongitude -= TWO_PI;
186  LocalCart_Origin_Long = originLongitude;
187  LocalCart_Origin_Height = originHeight;
188  if (orientation > PI)
189  orientation -= TWO_PI;
190  LocalCart_Orientation = orientation;
191  es2 = 2 * flattening - flattening * flattening;
192 
193  Sin_LocalCart_Origin_Lat = sin(LocalCart_Origin_Lat);
194  Cos_LocalCart_Origin_Lat = cos(LocalCart_Origin_Lat);
195  Sin_LocalCart_Origin_Lon = sin(LocalCart_Origin_Long);
196  Cos_LocalCart_Origin_Lon = cos(LocalCart_Origin_Long);
197  Sin_LocalCart_Orientation = sin(LocalCart_Orientation);
198  Cos_LocalCart_Orientation = cos(LocalCart_Orientation);
199 
200  Sin_Lat_Sin_Orient = Sin_LocalCart_Origin_Lat * Sin_LocalCart_Orientation;
201  Sin_Lat_Cos_Orient = Sin_LocalCart_Origin_Lat * Cos_LocalCart_Orientation;
202 
203  N0 = semiMajorAxis / sqrt(1 - es2 * Sin_LocalCart_Origin_Lat * Sin_LocalCart_Origin_Lat);
204 
205  val = (N0 + LocalCart_Origin_Height) * Cos_LocalCart_Origin_Lat;
206  u0 = val * Cos_LocalCart_Origin_Lon;
207  v0 = val * Sin_LocalCart_Origin_Lon;
208  w0 = ((N0 * (1 - es2)) + LocalCart_Origin_Height) * Sin_LocalCart_Origin_Lat;
209 
210  geocentric = new Geocentric( semiMajorAxis, flattening );
211 }
212 
213 
215 {
216  geocentric = new Geocentric( *( lc.geocentric ) );
218  flattening = lc.flattening;
219  es2 = lc.es2;
220  u0 = lc.u0;
221  v0 = lc.v0;
222  w0 = lc.w0;
223  LocalCart_Origin_Long = lc.LocalCart_Origin_Long;
224  LocalCart_Origin_Lat = lc.LocalCart_Origin_Lat;
225  LocalCart_Origin_Height = lc.LocalCart_Origin_Height;
226  LocalCart_Orientation = lc.LocalCart_Orientation;
227  Sin_LocalCart_Origin_Lat = lc.Sin_LocalCart_Origin_Lat;
228  Cos_LocalCart_Origin_Lat = lc.Cos_LocalCart_Origin_Lat;
229  Sin_LocalCart_Origin_Lon = lc.Sin_LocalCart_Origin_Lon;
230  Cos_LocalCart_Origin_Lon = lc.Cos_LocalCart_Origin_Lon;
231  Sin_LocalCart_Orientation = lc.Sin_LocalCart_Orientation;
232  Cos_LocalCart_Orientation = lc.Cos_LocalCart_Orientation;
233  Sin_Lat_Sin_Orient = lc.Sin_Lat_Sin_Orient;
234  Sin_Lat_Cos_Orient = lc.Sin_Lat_Cos_Orient;
235 }
236 
237 
239 {
240  delete geocentric;
241  geocentric = 0;
242 }
243 
244 
246 {
247  if( this != &lc )
248  {
249  geocentric->operator=( *lc.geocentric );
251  flattening = lc.flattening;
252  es2 = lc.es2;
253  u0 = lc.u0;
254  v0 = lc.v0;
255  w0 = lc.w0;
256  LocalCart_Origin_Long = lc.LocalCart_Origin_Long;
257  LocalCart_Origin_Lat = lc.LocalCart_Origin_Lat;
258  LocalCart_Origin_Height = lc.LocalCart_Origin_Height;
259  LocalCart_Orientation = lc.LocalCart_Orientation;
260  Sin_LocalCart_Origin_Lat = lc.Sin_LocalCart_Origin_Lat;
261  Cos_LocalCart_Origin_Lat = lc.Cos_LocalCart_Origin_Lat;
262  Sin_LocalCart_Origin_Lon = lc.Sin_LocalCart_Origin_Lon;
263  Cos_LocalCart_Origin_Lon = lc.Cos_LocalCart_Origin_Lon;
264  Sin_LocalCart_Orientation = lc.Sin_LocalCart_Orientation;
265  Cos_LocalCart_Orientation = lc.Cos_LocalCart_Orientation;
266  Sin_Lat_Sin_Orient = lc.Sin_Lat_Sin_Orient;
267  Sin_Lat_Cos_Orient = lc.Sin_Lat_Cos_Orient;
268  }
269 
270  return *this;
271 }
272 
273 
275 {
276 /*
277  * The function getParameters returns the ellipsoid parameters
278  * and local origin parameters.
279  *
280  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
281  * ellipsoidFlattening : Flattening of ellipsoid (output)
282  * originLongitude : Longitude of the local origin, in radians (output)
283  * originLatitude : Latitude of the local origin, in radians (output)
284  * originHeight : Ellipsoid height of the local origin, in meters (output)
285  * orientation : Orientation angle of the local cartesian coordinate system,
286  * in radians (output)
287  */
288 
289  return new LocalCartesianParameters( CoordinateType::localCartesian, LocalCart_Origin_Long, LocalCart_Origin_Lat, LocalCart_Origin_Height, LocalCart_Orientation );
290 }
291 
292 
294 {
295 /*
296  * The function convertFromGeodetic converts geodetic coordinates
297  * (latitude, longitude, and height) to local cartesian coordinates (X, Y, Z),
298  * according to the current ellipsoid and local origin parameters.
299  *
300  * longitude : Geodetic longitude, in radians (input)
301  * latitude : Geodetic latitude, in radians (input)
302  * Height : Geodetic height, in meters (input)
303  * X : Calculated local cartesian X coordinate, in meters (output)
304  * Y : Calculated local cartesian Y coordinate, in meters (output)
305  * Z : Calculated local cartesian Z coordinate, in meters (output)
306  *
307  */
308 
309  double longitude = geodeticCoordinates->longitude();
310  double latitude = geodeticCoordinates->latitude();
311  double height = geodeticCoordinates->height();
312 
313  if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
314  { /* geodetic latitude out of range */
316  }
317  if ((longitude < -PI) || (longitude > TWO_PI))
318  { /* geodetic longitude out of range */
320  }
321 
322  CartesianCoordinates* geocentricCoordinates = 0;
323  CartesianCoordinates* cartesianCoordinates = 0;
324  try
325  {
326  geocentricCoordinates = geocentric->convertFromGeodetic(
327  geodeticCoordinates );
328  cartesianCoordinates = convertFromGeocentric( geocentricCoordinates );
329  delete geocentricCoordinates;
330  }
332  {
333  delete geocentricCoordinates;
334  delete cartesianCoordinates;
335  throw e;
336  }
337 
338  return cartesianCoordinates;
339 }
340 
341 
343  MSP::CCS::CartesianCoordinates* cartesianCoordinates )
344 {
345 /*
346  * The function convertToGeodetic converts local cartesian
347  * coordinates (X, Y, Z) to geodetic coordinates (latitude, longitude,
348  * and height), according to the current ellipsoid and local origin parameters.
349  *
350  * X : Local cartesian X coordinate, in meters (input)
351  * Y : Local cartesian Y coordinate, in meters (input)
352  * Z : Local cartesian Z coordinate, in meters (input)
353  * longitude : Calculated longitude value, in radians (output)
354  * latitude : Calculated latitude value, in radians (output)
355  * Height : Calculated height value, in meters (output)
356  */
357  CartesianCoordinates* geocentricCoordinates = 0;
358  GeodeticCoordinates* geodeticCoordinates = 0;
359  try {
360  geocentricCoordinates = convertToGeocentric( cartesianCoordinates );
361  geodeticCoordinates = geocentric->convertToGeodetic(
362  geocentricCoordinates );
363 
364  double longitude = geodeticCoordinates->longitude();
365 
366  if (longitude > PI)
367  geodeticCoordinates->setLongitude( longitude -= TWO_PI );
368 
369  longitude = geodeticCoordinates->longitude();
370 
371  if (longitude < -PI)
372  geodeticCoordinates->setLongitude( longitude += TWO_PI );
373 
374  delete geocentricCoordinates;
375  }
377  {
378  delete geocentricCoordinates;
379  delete geodeticCoordinates;
380  throw e;
381  }
382 
383  return geodeticCoordinates;
384 }
385 
386 
388  const MSP::CCS::CartesianCoordinates* cartesianCoordinates )
389 {
390 /*
391  * The function convertFromGeocentric converts geocentric
392  * coordinates according to the current ellipsoid and local origin parameters.
393  *
394  * U : Geocentric latitude, in meters (input)
395  * V : Geocentric longitude, in meters (input)
396  * W : Geocentric height, in meters (input)
397  * X : Calculated local cartesian X coordinate, in meters (output)
398  * Y : Calculated local cartesian Y coordinate, in meters (output)
399  * Z : Calculated local cartesian Z coordinate, in meters (output)
400  *
401  */
402 
403  double X, Y, Z;
404  double u_MINUS_u0, v_MINUS_v0, w_MINUS_w0;
405 
406  double U = cartesianCoordinates->x();
407  double V = cartesianCoordinates->y();
408  double W = cartesianCoordinates->z();
409 
410  u_MINUS_u0 = U - u0;
411  v_MINUS_v0 = V - v0;
412  w_MINUS_w0 = W - w0;
413 
414  if (LocalCart_Orientation == 0.0)
415  {
416  double cos_lon_u_MINUS_u0 = Cos_LocalCart_Origin_Lon * u_MINUS_u0;
417  double sin_lon_v_MINUS_v0 = Sin_LocalCart_Origin_Lon * v_MINUS_v0;
418 
419  X = -Sin_LocalCart_Origin_Lon * u_MINUS_u0 + Cos_LocalCart_Origin_Lon * v_MINUS_v0;
420  Y = -Sin_LocalCart_Origin_Lat * cos_lon_u_MINUS_u0 + -Sin_LocalCart_Origin_Lat * sin_lon_v_MINUS_v0 + Cos_LocalCart_Origin_Lat * w_MINUS_w0;
421  Z = Cos_LocalCart_Origin_Lat * cos_lon_u_MINUS_u0 + Cos_LocalCart_Origin_Lat * sin_lon_v_MINUS_v0 + Sin_LocalCart_Origin_Lat * w_MINUS_w0;
422  }
423  else
424  {
425  double cos_lat_w_MINUS_w0 = Cos_LocalCart_Origin_Lat * w_MINUS_w0;
426 
427  X = (-Cos_LocalCart_Orientation * Sin_LocalCart_Origin_Lon + Sin_Lat_Sin_Orient * Cos_LocalCart_Origin_Lon) * u_MINUS_u0 +
428  (Cos_LocalCart_Orientation * Cos_LocalCart_Origin_Lon + Sin_Lat_Sin_Orient * Sin_LocalCart_Origin_Lon) * v_MINUS_v0 +
429  (-Sin_LocalCart_Orientation * cos_lat_w_MINUS_w0);
430 
431  Y = (-Sin_LocalCart_Orientation * Sin_LocalCart_Origin_Lon - Sin_Lat_Cos_Orient * Cos_LocalCart_Origin_Lon) * u_MINUS_u0 +
432  (Sin_LocalCart_Orientation * Cos_LocalCart_Origin_Lon - Sin_Lat_Cos_Orient * Sin_LocalCart_Origin_Lon) * v_MINUS_v0 +
433  (Cos_LocalCart_Orientation * cos_lat_w_MINUS_w0);
434 
435  Z = (Cos_LocalCart_Origin_Lat * Cos_LocalCart_Origin_Lon) * u_MINUS_u0 +
436  (Cos_LocalCart_Origin_Lat * Sin_LocalCart_Origin_Lon) * v_MINUS_v0 +
437  Sin_LocalCart_Origin_Lat * w_MINUS_w0;
438  }
439 
441 }
442 
443 
445 {
446 /*
447  * The function Convert_Local_Cartesian_To_Geocentric converts local cartesian
448  * coordinates (x, y, z) to geocentric coordinates (X, Y, Z) according to the
449  * current ellipsoid and local origin parameters.
450  *
451  * X : Local cartesian X coordinate, in meters (input)
452  * Y : Local cartesian Y coordinate, in meters (input)
453  * Z : Local cartesian Z coordinate, in meters (input)
454  * U : Calculated U value, in meters (output)
455  * V : Calculated v value, in meters (output)
456  * W : Calculated w value, in meters (output)
457  */
458 
459  double U, V, W;
460 
461  double X = cartesianCoordinates->x();
462  double Y = cartesianCoordinates->y();
463  double Z = cartesianCoordinates->z();
464 
465  if (LocalCart_Orientation == 0.0)
466  {
467  double sin_lat_y = Sin_LocalCart_Origin_Lat * Y;
468  double cos_lat_z = Cos_LocalCart_Origin_Lat * Z;
469 
470  U = -Sin_LocalCart_Origin_Lon * X - sin_lat_y * Cos_LocalCart_Origin_Lon + cos_lat_z * Cos_LocalCart_Origin_Lon + u0;
471  V = Cos_LocalCart_Origin_Lon * X - sin_lat_y * Sin_LocalCart_Origin_Lon + cos_lat_z * Sin_LocalCart_Origin_Lon + v0;
472  W = Cos_LocalCart_Origin_Lat * Y + Sin_LocalCart_Origin_Lat * Z + w0;
473  }
474  else
475  {
476  double rotated_x, rotated_y;
477  double rotated_y_sin_lat, z_cos_lat;
478 
479  rotated_x = Cos_LocalCart_Orientation * X + Sin_LocalCart_Orientation * Y;
480  rotated_y = -Sin_LocalCart_Orientation * X + Cos_LocalCart_Orientation * Y;
481 
482  rotated_y_sin_lat = rotated_y * Sin_LocalCart_Origin_Lat;
483  z_cos_lat = Z * Cos_LocalCart_Origin_Lat;
484 
485  U = -Sin_LocalCart_Origin_Lon * rotated_x - Cos_LocalCart_Origin_Lon * rotated_y_sin_lat + Cos_LocalCart_Origin_Lon * z_cos_lat + u0;
486  V = Cos_LocalCart_Origin_Lon * rotated_x - Sin_LocalCart_Origin_Lon * rotated_y_sin_lat + Sin_LocalCart_Origin_Lon * z_cos_lat + v0;
487  W = Cos_LocalCart_Origin_Lat * rotated_y + Sin_LocalCart_Origin_Lat * Z + w0;
488  }
489 
490  return new CartesianCoordinates( CoordinateType::geocentric, U, V, W );
491 }
492 
493 // CLASSIFICATION: UNCLASSIFIED