UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
DatumLibraryImplementation.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: Datum Library Implementation
5  *
6  * ABSTRACT
7  *
8  * This component provides datum shifts for a large collection of local
9  * datums, WGS72, and WGS84. A particular datum can be accessed by using its
10  * standard 5-letter code to find its index in the datum table. The index
11  * can then be used to retrieve the name, type, ellipsoid code, and datum
12  * shift parameters, and to perform shifts to or from that datum.
13  *
14  * By sequentially retrieving all of the datum codes and/or names, a menu
15  * of the available datums can be constructed. The index values resulting
16  * from selections from this menu can then be used to access the parameters
17  * of the selected datum, or to perform datum shifts involving that datum.
18  *
19  * This component supports both 3-parameter local datums, for which only X,
20  * Y, and Z translations relative to WGS 84 have been defined, and
21  * 7-parameter local datums, for which X, Y, and Z rotations, and a scale
22  * factor, are also defined. It also includes entries for WGS 84 (with an
23  * index of 0), and WGS 72 (with an index of 1), but no shift parameter
24  * values are defined for these.
25  *
26  * This component provides datum shift functions for both geocentric and
27  * geodetic coordinates. WGS84 is used as an intermediate state when
28  * shifting from one local datum to another. When geodetic coordinates are
29  * given Molodensky's method is used, except near the poles where the 3-step
30  * step method is used instead. Specific algorithms are used for shifting
31  * between WGS72 and WGS84.
32  *
33  * This component depends on two data files, named 3_param.dat and
34  * 7_param.dat, which contain the datum parameter values. Copies of these
35  * files must be located in the directory specified by the value of the
36  * environment variable "MSPCCS_DATA", if defined, or else in the current
37  * directory whenever a program containing this component is executed.
38  *
39  * Additional datums can be added to these files, either manually or using
40  * the Create_Datum function. However, if a large number of datums are
41  * added, the datum table array sizes in this component will have to be
42  * increased.
43  *
44  * This component depends on two other components: the Ellipsoid component
45  * for access to ellipsoid parameters; and the Geocentric component for
46  * conversions between geodetic and geocentric coordinates.
47  *
48  * ERROR HANDLING
49  *
50  * This component checks for input file errors and input parameter errors.
51  * If an invalid value is found, the error code is combined with the current
52  * error code using the bitwise or. This combining allows multiple error
53  * codes to be returned. The possible error codes are:
54  *
55  * DATUM_NO_ERROR : No errors occurred in function
56  * DATUM_NOT_INITIALIZED_ERROR : Datum module has not been initialized
57  * DATUM_7PARAM_FILE_OPEN_ERROR : 7 parameter file opening error
58  * DATUM_7PARAM_FILE_PARSING_ERROR : 7 parameter file structure error
59  * DATUM_3PARAM_FILE_OPEN_ERROR : 3 parameter file opening error
60  * DATUM_3PARAM_FILE_PARSING_ERROR : 3 parameter file structure error
61  * DATUM_INVALID_INDEX_ERROR : Index out of valid range (less than one
62  * or more than Number_of_Datums)
63  * DATUM_INVALID_SRC_INDEX_ERROR : Source datum index invalid
64  * DATUM_INVALID_DEST_INDEX_ERROR : Destination datum index invalid
65  * DATUM_INVALID_CODE_ERROR : Datum code not found in table
66  * DATUM_LAT_ERROR : Latitude out of valid range (-90 to 90)
67  * DATUM_LON_ERROR : Longitude out of valid range (-180 to
68  * 360)
69  * DATUM_SIGMA_ERROR : Standard error values must be positive
70  * (or -1 if unknown)
71  * DATUM_DOMAIN_ERROR : Domain of validity not well defined
72  * DATUM_ELLIPSE_ERROR : Error in ellipsoid module
73  * DATUM_NOT_USERDEF_ERROR : Datum code is not user defined - cannot
74  * be deleted
75  *
76  *
77  * REUSE NOTES
78  *
79  * Datum is intended for reuse by any application that needs access to
80  * datum shift parameters relative to WGS 84.
81  *
82  *
83  * REFERENCES
84  *
85  * Further information on Datum can be found in the Reuse Manual.
86  *
87  * Datum originated from : U.S. Army Topographic Engineering Center (USATEC)
88  * Geospatial Information Division (GID)
89  * 7701 Telegraph Road
90  * Alexandria, VA 22310-3864
91  *
92  * LICENSES
93  *
94  * None apply to this component.
95  *
96  * RESTRICTIONS
97  *
98  * Datum has no restrictions.
99  *
100  * ENVIRONMENT
101  *
102  * Datum was tested and certified in the following environments:
103  *
104  * 1. Solaris 2.5 with GCC 2.8.1
105  * 2. MS Windows 95 with MS Visual C++ 6
106  *
107  * MODIFICATIONS
108  *
109  * Date Description
110  * ---- -----------
111  * 03/30/97 Original Code
112  * 05/28/99 Added user-definable datums (for JMTK)
113  * Added datum domain of validity checking (for JMTK)
114  * Added datum shift accuracy calculation (for JMTK)
115  * 06/27/06 Moved data files to data directory
116  * 03-14-07 Original C++ Code
117  * 03-21-08 Added check for west, east longitude values > 180
118  * 06-11-10 S. Gillis, BAEts26724, Fixed memory error problem
119  * when MSPCCS_DATA is not set
120  * 07-07-10 K.Lam, BAEts27269, Replace C functions in threads.h
121  * with C++ methods in classes CCSThreadMutex
122  * 05/17/11 T. Thompson, BAEts27393, let user know when problem
123  * is due to undefined MSPCCS_DATA
124  * 07/13/12 K.Lam, BAEts29544, fixed problem with create datum
125  * 07/17/12 S.Gillis,MSP_00029561,Fixed problem with deleting datum
126  * 08/13/12 S. Gillis, MSP_00029654, Added lat/lon to define7ParamDatum
127  */
128 
129 
130 /***************************************************************************/
131 /*
132  * INCLUDES
133  */
134 
135 #include <math.h>
136 #include <stdio.h>
137 #include <stdlib.h>
138 #include <ctype.h>
139 #include "DatumLibraryImplementation.h"
140 #include "EllipsoidLibraryImplementation.h"
141 #include "EllipsoidParameters.h"
142 #include "SevenParameterDatum.h"
143 #include "ThreeParameterDatum.h"
144 #include "Geocentric.h"
145 #include "Datum.h"
146 #include "CartesianCoordinates.h"
147 #include "GeodeticCoordinates.h"
148 #include "CoordinateConversionException.h"
149 #include "ErrorMessages.h"
150 #include "Accuracy.h"
151 #include "CCSThreadMutex.h"
152 #include "CCSThreadLock.h"
153 
154 /*
155  * math.h - standard C mathematics library
156  * EllipsoidLibrary.h - used to get ellipsoid parameters
157  * SevenParameterDatum.h - creates a 7 parameter datum
158  * ThreeParameterDatum.h - creates a 3 parameter datum
159  * Geocentric.h - converts between geodetic and geocentric coordinates
160  * Datum.h - used to store individual datum information
161  * DatumLibraryImplementation.h - for error ehecking and error codes
162  * CCSThreadMutex.h - used for thread safety
163  * CCSThreadLock.h - used for thread safety
164  * CartesianCoordinates.h - defines cartesian coordinates
165  * GeodeticCoordinates.h - defines geodetic coordinates
166  * CoordinateConversionException.h - Exception handler
167  * ErrorMessages.h - Contains exception messages
168  */
169 
170 
171 using namespace MSP::CCS;
172 using MSP::CCSThreadMutex;
173 using MSP::CCSThreadLock;
174 
175 /***************************************************************************/
176 /*
177  * DEFINES
178  */
179 
180 const double SECONDS_PER_RADIAN = 206264.8062471; /* Seconds in a radian */
181 const double PI = 3.14159265358979323e0;
182 const double PI_OVER_2 = (PI / 2.0);
183 const double PI_OVER_180 = (PI / 180.0);
184 const double _180_OVER_PI = (180.0 / PI);
185 const double TWO_PI = (2.0 * PI);
186 const double MIN_LAT = (-PI/2.0);
187 const double MAX_LAT = (+PI/2.0);
188 const double MIN_LON = -PI;
189 const double MAX_LON = (2.0 * PI);
190 const int DATUM_CODE_LENGTH = 7;
191 const int DATUM_NAME_LENGTH = 33;
192 const int ELLIPSOID_CODE_LENGTH = 3;
193 const int MAX_WGS = 2;
194 const double MOLODENSKY_MAX = (89.75 * PI_OVER_180); /* Polar limit */
195 const int FILENAME_LENGTH = 128;
196 const char *WGS84_Datum_Code = "WGE";
197 const char *WGS72_Datum_Code = "WGC";
198 
199 
200 /************************************************************************/
201 /* LOCAL FUNCTIONS
202  *
203  */
204 
205 GeodeticCoordinates* molodenskyShift(
206  const double a,
207  const double da,
208  const double f,
209  const double df,
210  const double dx,
211  const double dy,
212  const double dz,
213  const double sourceLongitude,
214  const double sourceLatitude,
215  const double sourceHeight )
216 {
217 /*
218  * The function molodenskyShift shifts geodetic coordinates
219  * using the Molodensky method.
220  *
221  * a : Semi-major axis of source ellipsoid in meters (input)
222  * da : Destination a minus source a (input)
223  * f : Flattening of source ellipsoid (input)
224  * df : Destination f minus source f (input)
225  * dx : X coordinate shift in meters (input)
226  * dy : Y coordinate shift in meters (input)
227  * dz : Z coordinate shift in meters (input)
228  * sourceLongitude : Longitude in radians (input)
229  * sourceLatitude : Latitude in radians. (input)
230  * sourceHeight : Height in meters. (input)
231  * targetLongitude : Calculated longitude in radians. (output)
232  * targetLatitude : Calculated latitude in radians. (output)
233  * targetHeight : Calculated height in meters. (output)
234  */
235 
236  double tLon_in; /* temp longitude */
237  double e2; /* Intermediate calculations for dp, dl */
238  double ep2; /* Intermediate calculations for dp, dl */
239  double sin_Lat; /* sin(Latitude_1) */
240  double sin2_Lat; /* (sin(Latitude_1))^2 */
241  double sin_Lon; /* sin(Longitude_1) */
242  double cos_Lat; /* cos(Latitude_1) */
243  double cos_Lon; /* cos(Longitude_1) */
244  double w2; /* Intermediate calculations for dp, dl */
245  double w; /* Intermediate calculations for dp, dl */
246  double w3; /* Intermediate calculations for dp, dl */
247  double m; /* Intermediate calculations for dp, dl */
248  double n; /* Intermediate calculations for dp, dl */
249  double dp; /* Delta phi */
250  double dp1; /* Delta phi calculations */
251  double dp2; /* Delta phi calculations */
252  double dp3; /* Delta phi calculations */
253  double dl; /* Delta lambda */
254  double dh; /* Delta height */
255  double dh1; /* Delta height calculations */
256  double dh2; /* Delta height calculations */
257 
258  if (sourceLongitude > PI)
259  tLon_in = sourceLongitude - TWO_PI;
260  else
261  tLon_in = sourceLongitude;
262 
263  e2 = 2 * f - f * f;
264  ep2 = e2 / (1 - e2);
265  sin_Lat = sin(sourceLatitude);
266  cos_Lat = cos(sourceLatitude);
267  sin_Lon = sin(tLon_in);
268  cos_Lon = cos(tLon_in);
269  sin2_Lat = sin_Lat * sin_Lat;
270  w2 = 1.0 - e2 * sin2_Lat;
271  w = sqrt(w2);
272  w3 = w * w2;
273  m = (a * (1.0 - e2)) / w3;
274  n = a / w;
275  dp1 = cos_Lat * dz - sin_Lat * cos_Lon * dx - sin_Lat * sin_Lon * dy;
276  dp2 = ((e2 * sin_Lat * cos_Lat) / w) * da;
277  dp3 = sin_Lat * cos_Lat * (2.0 * n + ep2 * m * sin2_Lat) * (1.0 - f) * df;
278  dp = (dp1 + dp2 + dp3) / (m + sourceHeight);
279  dl = (-sin_Lon * dx + cos_Lon * dy) / ((n + sourceHeight) * cos_Lat);
280  dh1 = (cos_Lat * cos_Lon * dx) + (cos_Lat * sin_Lon * dy) + (sin_Lat * dz);
281  dh2 = -(w * da) + ((a * (1 - f)) / w) * sin2_Lat * df;
282  dh = dh1 + dh2;
283 
284  double targetLatitude = sourceLatitude + dp;
285  double targetLongitude = sourceLongitude + dl;
286  double targetHeight = sourceHeight + dh;
287 
288  if (targetLongitude > TWO_PI)
289  targetLongitude -= TWO_PI;
290  if (targetLongitude < (- PI))
291  targetLongitude += TWO_PI;
292 
293  return new GeodeticCoordinates(
294  CoordinateType::geodetic, targetLongitude, targetLatitude, targetHeight );
295 }
296 
297 
298 /************************************************************************/
299 /* FUNCTIONS
300  *
301  */
302 
303 /* This class is a safeguard to make sure the singleton gets deleted
304  * when the application exits
305  */
306 namespace MSP
307 {
308  namespace CCS
309  {
310  class DatumLibraryImplementationCleaner
311  {
312  public:
313 
314  ~DatumLibraryImplementationCleaner()
315  {
316  CCSThreadLock lock(&DatumLibraryImplementation::mutex);
317  DatumLibraryImplementation::deleteInstance();
318  }
319 
320  } datumLibraryImplementationCleanerInstance;
321  }
322 }
323 
324 // Make this class a singleton, so the data files are only initialized once
325 CCSThreadMutex DatumLibraryImplementation::mutex;
326 DatumLibraryImplementation* DatumLibraryImplementation::instance = 0;
327 int DatumLibraryImplementation::instanceCount = 0;
328 
329 
330 DatumLibraryImplementation* DatumLibraryImplementation::getInstance()
331 {
332  CCSThreadLock lock(&mutex);
333  if( instance == 0 )
334  instance = new DatumLibraryImplementation;
335 
336  instanceCount++;
337 
338  return instance;
339 }
340 
341 
342 void DatumLibraryImplementation::removeInstance()
343 {
344 /*
345  * The function removeInstance removes this DatumLibraryImplementation
346  * instance from the total number of instances.
347  */
348  CCSThreadLock lock(&mutex);
349  if( --instanceCount < 1 )
350  {
351  deleteInstance();
352  }
353 }
354 
355 
356 void DatumLibraryImplementation::deleteInstance()
357 {
358 /*
359  * Delete the singleton.
360  */
361 
362  if( instance != 0 )
363  {
364  delete instance;
365  instance = 0;
366  }
367 }
368 
369 
370 DatumLibraryImplementation::DatumLibraryImplementation():
371  _ellipsoidLibraryImplementation( 0 ),
372  datum3ParamCount( 0 ),
373  datum7ParamCount( 0 )
374 {
375  loadDatums();
376 }
377 
378 
379 DatumLibraryImplementation::DatumLibraryImplementation(
380  const DatumLibraryImplementation &dl )
381 {
382  int size = dl.datumList.size();
383  for( int i = 0; i < size; i++ )
384  {
385  switch( dl.datumList[i]->datumType() )
386  {
387  case DatumType::threeParamDatum:
388  datumList.push_back( new ThreeParameterDatum(
389  *( ( ThreeParameterDatum* )( dl.datumList[i] ) ) ) );
390  break;
391  case DatumType::sevenParamDatum:
392  datumList.push_back( new SevenParameterDatum(
393  *( ( SevenParameterDatum* )( dl.datumList[i] ) ) ) );
394  break;
395  case DatumType::wgs84Datum:
396  case DatumType::wgs72Datum:
397  datumList.push_back( new Datum( *( dl.datumList[i] ) ) );
398  break;
399  }
400  }
401 
402  _ellipsoidLibraryImplementation = dl._ellipsoidLibraryImplementation;
403  datum3ParamCount = dl.datum3ParamCount;
404  datum7ParamCount = dl.datum7ParamCount;
405 }
406 
407 
408 DatumLibraryImplementation::~DatumLibraryImplementation()
409 {
410  std::vector<Datum*>::iterator iter = datumList.begin();
411  while( iter != datumList.end() )
412  {
413  delete( *iter );
414  iter++;
415  }
416  datumList.clear();
417 
418  _ellipsoidLibraryImplementation = 0;
419 }
420 
421 
422 DatumLibraryImplementation& DatumLibraryImplementation::operator=(
423  const DatumLibraryImplementation &dl )
424 {
425  if ( &dl == this )
426  return *this;
427 
428  int size = dl.datumList.size();
429  for( int i = 0; i < size; i++ )
430  {
431  switch( dl.datumList[i]->datumType() )
432  {
433  case DatumType::threeParamDatum:
434  datumList.push_back( new ThreeParameterDatum(
435  *( ( ThreeParameterDatum* )( dl.datumList[i] ) ) ) );
436  break;
437  case DatumType::sevenParamDatum:
438  datumList.push_back( new SevenParameterDatum(
439  *( ( SevenParameterDatum* )( dl.datumList[i] ) ) ) );
440  break;
441  case DatumType::wgs84Datum:
442  case DatumType::wgs72Datum:
443  datumList.push_back( new Datum( *( dl.datumList[i] ) ) );
444  break;
445  }
446  }
447 
448  _ellipsoidLibraryImplementation = dl._ellipsoidLibraryImplementation;
449  datum3ParamCount = dl.datum3ParamCount;
450  datum7ParamCount = dl.datum7ParamCount;
451 
452  return *this;
453 }
454 
455 
456 void DatumLibraryImplementation::define3ParamDatum(
457  const char *code,
458  const char *name,
459  const char *ellipsoidCode,
460  double deltaX,
461  double deltaY,
462  double deltaZ,
463  double sigmaX,
464  double sigmaY,
465  double sigmaZ,
466  double westLongitude,
467  double eastLongitude,
468  double southLatitude,
469  double northLatitude )
470 {
471 /*
472  * The function define3ParamDatum creates a new local 3-parameter datum with the
473  * specified code, name, and axes. If the datum table has not been initialized,
474  * the specified code is already in use, or a new version of the 3-param.dat
475  * file cannot be created, an exception is thrown.
476  * Note that the indexes of all datums in the datum table may be
477  * changed by this function.
478  *
479  * code : 5-letter new datum code. (input)
480  * name : Name of the new datum (input)
481  * ellipsoidCode : 2-letter code for the associated ellipsoid (input)
482  * deltaX : X translation to WGS84 in meters (input)
483  * deltaY : Y translation to WGS84 in meters (input)
484  * deltaZ : Z translation to WGS84 in meters (input)
485  * sigmaX : Standard error in X in meters (input)
486  * sigmaY : Standard error in Y in meters (input)
487  * sigmaZ : Standard error in Z in meters (input)
488  * westLongitude : Western edge of validity rectangle in radians (input)
489  * eastLongitude : Eastern edge of validity rectangle in radians (input)
490  * southLatitude : Southern edge of validity rectangle in radians(input)
491  * northLatitude : Northern edge of validity rectangle in radians(input)
492  */
493 
494  char datum_Code[DATUM_CODE_LENGTH];
495  long index = 0;
496  long ellipsoid_index = 0;
497  long code_length = 0;
498  char *PathName = NULL;
499  FILE *fp_3param = NULL;
500 
501  if (!(((sigmaX > 0.0) || (sigmaX == -1.0)) &&
502  ((sigmaY > 0.0) || (sigmaY == -1.0)) &&
503  ((sigmaZ > 0.0) || (sigmaZ == -1.0))))
504  throw CoordinateConversionException( ErrorMessages::datumSigma );
505 
506  if ((southLatitude < MIN_LAT) || (southLatitude > MAX_LAT))
507  throw CoordinateConversionException( ErrorMessages::latitude );
508  if ((westLongitude < MIN_LON) || (westLongitude > MAX_LON))
509  throw CoordinateConversionException( ErrorMessages::longitude );
510  if ((northLatitude < MIN_LAT) || (northLatitude > MAX_LAT))
511  throw CoordinateConversionException( ErrorMessages::latitude );
512  if ((eastLongitude < MIN_LON) || (eastLongitude > MAX_LON))
513  throw CoordinateConversionException( ErrorMessages::longitude );
514  if (southLatitude >= northLatitude)
515  throw CoordinateConversionException( ErrorMessages::datumDomain );
516  if((westLongitude >= eastLongitude) &&
517  (westLongitude >= 0 && westLongitude < 180) &&
518  (eastLongitude >= 0 && eastLongitude < 180))
519  throw CoordinateConversionException( ErrorMessages::datumDomain );
520 
521  // assume the datum code is new
522  bool isNewDatumCode = true;
523  try
524  {
525  // check if datum code exists
526  datumIndex( code, &index );
527  // get here if datum code is found in current datum table
528  isNewDatumCode = false;
529  }
530  catch(CoordinateConversionException e)
531  {
532  // the datum code is new, keep going
533  }
534 
535  // the datum code exists in current datum table, throw an error
536  if ( !isNewDatumCode )
537  throw CoordinateConversionException( ErrorMessages::invalidDatumCode );
538 
539  code_length = strlen( code );
540 
541  if( code_length > ( DATUM_CODE_LENGTH-1 ) )
542  throw CoordinateConversionException( ErrorMessages::invalidDatumCode );
543 
544  if( _ellipsoidLibraryImplementation )
545  {
546  _ellipsoidLibraryImplementation->ellipsoidIndex(
547  ellipsoidCode, &ellipsoid_index );
548  }
549  else
550  throw CoordinateConversionException( ErrorMessages::ellipse );
551 
552 
553  strcpy( datum_Code, code );
554  /* Convert code to upper case */
555  for( long i = 0; i < code_length; i++ )
556  datum_Code[i] = ( char )toupper( datum_Code[i] );
557 
558  int numDatums = datumList.size();
559  datumList.push_back( new ThreeParameterDatum(
560  numDatums, ( char* )datum_Code, ( char* )ellipsoidCode,
561  ( char* )name, DatumType::threeParamDatum, deltaX, deltaY, deltaZ,
562  westLongitude, eastLongitude, southLatitude, northLatitude, sigmaX,
563  sigmaY, sigmaZ, true ) );
564  datum3ParamCount++;
565 
566  write3ParamFile();
567 }
568 
569 
570 void DatumLibraryImplementation::define7ParamDatum(
571  const char *code,
572  const char *name,
573  const char *ellipsoidCode,
574  double deltaX,
575  double deltaY,
576  double deltaZ,
577  double rotationX,
578  double rotationY,
579  double rotationZ,
580  double scale,
581  double westLongitude,
582  double eastLongitude,
583  double southLatitude,
584  double northLatitude )
585 {
586 /*
587  * The function define7ParamDatum creates a new local 7-parameter datum with the
588  * specified code, name, and axes. If the datum table has not been initialized,
589  * the specified code is already in use, or a new version of the 7-param.dat
590  * file cannot be created, an exception is thrown.
591  * Note that the indexes of all datums in the datum table may be
592  * changed by this function.
593  *
594  * code : 5-letter new datum code. (input)
595  * name : Name of the new datum (input)
596  * ellipsoidCode : 2-letter code for the associated ellipsoid (input)
597  * deltaX : X translation to WGS84 in meters (input)
598  * deltaY : Y translation to WGS84 in meters (input)
599  * deltaZ : Z translation to WGS84 in meters (input)
600  * rotationX : X rotation to WGS84 in arc seconds (input)
601  * rotationY : Y rotation to WGS84 in arc seconds (input)
602  * rotationZ : Z rotation to WGS84 in arc seconds (input)
603  * scale : Scale factor (input)
604  * westLongitude : Western edge of validity rectangle in radians (input)
605  * eastLongitude : Eastern edge of validity rectangle in radians (input)
606  * southLatitude : Southern edge of validity rectangle in radians(input)
607  * northLatitude : Northern edge of validity rectangle in radians(input)
608  */
609 
610  char datum_Code[DATUM_CODE_LENGTH];
611  long index = 0;
612  long ellipsoid_index = 0;
613  long code_length = 0;
614  char *PathName = NULL;
615  FILE *fp_7param = NULL;
616 
617  if ((rotationX < -60.0) || (rotationX > 60.0))
618  throw CoordinateConversionException( ErrorMessages::datumRotation );
619  if ((rotationY < -60.0) || (rotationY > 60.0))
620  throw CoordinateConversionException( ErrorMessages::datumRotation );
621  if ((rotationZ < -60.0) || (rotationZ > 60.0))
622  throw CoordinateConversionException( ErrorMessages::datumRotation );
623 
624  if ((scale < -0.001) || (scale > 0.001))
625  throw CoordinateConversionException( ErrorMessages::scaleFactor );
626 
627  // assume the datum code is new
628  bool isNewDatumCode = true;
629  try
630  {
631  // check if datum code exists
632  datumIndex( code, &index );
633  // get here if datum code is found in current datum table
634  isNewDatumCode = false;
635  }
636  catch(CoordinateConversionException e)
637  {
638  // the datum code is new, keep going
639  }
640 
641  // the datum code exists in current datum table, throw an error
642  if ( !isNewDatumCode )
643  throw CoordinateConversionException( ErrorMessages::invalidDatumCode );
644 
645  code_length = strlen( code );
646  if( code_length > ( DATUM_CODE_LENGTH-1 ) )
647  throw CoordinateConversionException( ErrorMessages::invalidDatumCode );
648 
649  if( _ellipsoidLibraryImplementation )
650  {
651  _ellipsoidLibraryImplementation->ellipsoidIndex(
652  ellipsoidCode, &ellipsoid_index );
653  }
654  else
655  throw CoordinateConversionException( ErrorMessages::ellipse );
656 
657  long i;
658  strcpy( datum_Code, code );
659  /* Convert code to upper case */
660  for( i = 0; i < code_length; i++ )
661  datum_Code[i] = ( char )toupper( datum_Code[i] );
662 
663  datumList.insert( datumList.begin() + MAX_WGS + datum7ParamCount,
664  new SevenParameterDatum( datum7ParamCount, ( char* )datum_Code,
665  ( char* )ellipsoidCode, ( char* )name, DatumType::sevenParamDatum,
666  deltaX, deltaY, deltaZ,
667  westLongitude, eastLongitude, southLatitude, northLatitude,
668  rotationX / SECONDS_PER_RADIAN,
669  rotationY / SECONDS_PER_RADIAN, rotationZ / SECONDS_PER_RADIAN,
670  scale, true ) );
671  datum7ParamCount++;
672 
673  write7ParamFile();
674 }
675 
676 
677 void DatumLibraryImplementation::removeDatum( const char* code )
678 {
679 /*
680  * The function removeDatum deletes a local (3-parameter) datum with the
681  * specified code. If the datum table has not been initialized or a new
682  * version of the 3-param.dat file cannot be created, an exception is thrown,
683  * Note that the indexes of all datums
684  * in the datum table may be changed by this function.
685  *
686  * code : 5-letter datum code. (input)
687  *
688  */
689 
690  char *PathName = NULL;
691  FILE *fp_3param = NULL;
692  FILE *fp_7param = NULL;
693  long index = 0;
694  bool delete_3param_datum = true;
695 
696  datumIndex( code, &index );
697 
698  if( datumList[index]->datumType() == DatumType::threeParamDatum )
699  {
700  if( !( ( ThreeParameterDatum* )datumList[index] )->userDefined() )
701  throw CoordinateConversionException( ErrorMessages::notUserDefined );
702  }
703  else if( datumList[index]->datumType() == DatumType::sevenParamDatum )
704  {
705  delete_3param_datum = false;
706  if( !( ( SevenParameterDatum* )datumList[index] )->userDefined() )
707  throw CoordinateConversionException( ErrorMessages::notUserDefined );
708  }
709  else
710  throw CoordinateConversionException( ErrorMessages::notUserDefined );
711 
712  datumList.erase( datumList.begin() + index );
713 
714  if( !delete_3param_datum )
715  {
716  datum7ParamCount--;
717 
718  write7ParamFile();
719  }
720  else if( delete_3param_datum )
721  {
722  datum3ParamCount--;
723 
724  write3ParamFile();
725  }
726 }
727 
728 
729 void DatumLibraryImplementation::datumCount( long *count )
730 {
731 /*
732  * The function datumCount returns the number of Datums in the table
733  * if the table was initialized without error.
734  *
735  * count : number of datums in the datum table (output)
736  */
737 
738  *count = datumList.size();
739 }
740 
741 
742 void DatumLibraryImplementation::datumIndex( const char *code, long *index )
743 {
744 /*
745  * The function datumIndex returns the index of the datum with the
746  * specified code.
747  *
748  * code : The datum code being searched for. (input)
749  * index : The index of the datum in the table with the (output)
750  * specified code.
751  */
752 
753  char temp_code[DATUM_CODE_LENGTH];
754  long length;
755  long pos = 0;
756  long i = 0;
757 
758  *index = 0;
759 
760  if( !code )
761  throw CoordinateConversionException( ErrorMessages::invalidDatumCode );
762 
763  length = strlen( code );
764  if ( length > ( DATUM_CODE_LENGTH-1 ) )
765  throw CoordinateConversionException( ErrorMessages::invalidDatumCode );
766  else
767  {
768  strcpy( temp_code, code );
769 
770  /* Convert to upper case */
771  for( i=0; i < length; i++ )
772  temp_code[i] = ( char )toupper( temp_code[i] );
773 
774  /* Strip blank spaces */
775  while( pos < length )
776  {
777  if( isspace( temp_code[pos] ) )
778  {
779  for( i=pos; i <= length; i++ )
780  temp_code[i] = temp_code[i+1];
781  length -= 1;
782  }
783  else
784  pos += 1;
785  }
786 
787  int numDatums = datumList.size();
788  /* Search for code */
789  i = 0;
790  while( i < numDatums && strcmp( temp_code, datumList[i]->code() ) )
791  {
792  i++;
793  }
794  if( i == numDatums || strcmp( temp_code, datumList[i]->code() ) )
795  throw CoordinateConversionException( ErrorMessages::invalidDatumCode );
796  else
797  *index = i;
798  }
799 }
800 
801 
802 void DatumLibraryImplementation::datumCode( const long index, char *code )
803 {
804 /*
805  * The function datumCode returns the 5-letter code of the datum
806  * referenced by index.
807  *
808  * index : The index of a given datum in the datum table. (input)
809  * code : The datum Code of the datum referenced by Index. (output)
810  */
811 
812  if( index < 0 || index >= datumList.size() )
813  throw CoordinateConversionException( ErrorMessages::invalidIndex );
814  else
815  strcpy( code, datumList[index]->code() );
816 }
817 
818 
819 void DatumLibraryImplementation::datumName( const long index, char *name )
820 {
821 /*
822  * The function datumName returns the name of the datum referenced by
823  * index.
824  *
825  * index : The index of a given datum in the datum table. (input)
826  * name : The datum Name of the datum referenced by Index. (output)
827  */
828 
829  if( index < 0 || index >= datumList.size() )
830  throw CoordinateConversionException( ErrorMessages::invalidIndex );
831  else
832  strcpy( name, datumList[index]->name() );
833 }
834 
835 
836 void DatumLibraryImplementation::datumEllipsoidCode(
837  const long index, char *code )
838 {
839 /*
840  * The function datumEllipsoidCode returns the 2-letter ellipsoid code
841  * for the ellipsoid associated with the datum referenced by index.
842  *
843  * index : The index of a given datum in the datum table. (input)
844  * code : The ellipsoid code for the ellipsoid associated with (output)
845  * the datum referenced by index.
846  */
847 
848  if( index < 0 || index >= datumList.size() )
849  throw CoordinateConversionException( ErrorMessages::invalidIndex );
850  else
851  strcpy( code, datumList[index]->ellipsoidCode() );
852 }
853 
854 
855 void DatumLibraryImplementation::datumStandardErrors(
856  const long index,
857  double *sigmaX,
858  double *sigmaY,
859  double *sigmaZ )
860 {
861 /*
862  * The function datumStandardErrors returns the standard errors in X,Y, & Z
863  * for the datum referenced by index.
864  *
865  * index : The index of a given datum in the datum table (input)
866  * sigma_X : Standard error in X in meters (output)
867  * sigma_Y : Standard error in Y in meters (output)
868  * sigma_Z : Standard error in Z in meters (output)
869  */
870 
871  if( index < 0 || index >= datumList.size() )
872  throw CoordinateConversionException( ErrorMessages::invalidIndex );
873  else
874  {
875  Datum* datum = datumList[index];
876 
877  if( datum->datumType() == DatumType::threeParamDatum )
878  {
879  ThreeParameterDatum* threeParameterDatum = ( ThreeParameterDatum* )datum;
880  *sigmaX = threeParameterDatum->sigmaX();
881  *sigmaY = threeParameterDatum->sigmaY();
882  *sigmaZ = threeParameterDatum->sigmaZ();
883  }
884  else
885  throw CoordinateConversionException( ErrorMessages::invalidIndex );
886  }
887 }
888 
889 
890 void DatumLibraryImplementation::datumSevenParameters(
891  const long index,
892  double *rotationX,
893  double *rotationY,
894  double *rotationZ,
895  double *scaleFactor )
896 {
897 /*
898  * The function datumSevenParameters returns parameter values,
899  * used only by a seven parameter datum,
900  * for the datum referenced by index.
901  *
902  * index : The index of a given datum in the datum table. (input)
903  * rotationX : X rotation in radians (output)
904  * rotationY : Y rotation in radians (output)
905  * rotationZ : Z rotation in radians (output)
906  * scaleFactor : Scale factor (output)
907  */
908 
909  if( index < 0 || index >= datumList.size() )
910  throw CoordinateConversionException( ErrorMessages::invalidIndex );
911  else
912  {
913  Datum* datum = datumList[index];
914 
915  if( datum->datumType() == DatumType::sevenParamDatum )
916  {
917  SevenParameterDatum* sevenParameterDatum = ( SevenParameterDatum* )datum;
918 
919  *rotationX = sevenParameterDatum->rotationX();
920  *rotationY = sevenParameterDatum->rotationY();
921  *rotationZ = sevenParameterDatum->rotationZ();
922  *scaleFactor = sevenParameterDatum->scaleFactor();
923  }
924  else
925  throw CoordinateConversionException( ErrorMessages::invalidIndex );
926  }
927 }
928 
929 
930 void DatumLibraryImplementation::datumTranslationValues(
931  const long index,
932  double *deltaX,
933  double *deltaY,
934  double *deltaZ )
935 {
936 /*
937  * The function datumTranslationValues returns the translation values
938  * for the datum referenced by index.
939  *
940  * index : The index of a given datum in the datum table. (input)
941  * deltaX : X translation in meters (output)
942  * deltaY : Y translation in meters (output)
943  * deltaZ : Z translation in meters (output)
944  */
945 
946  if( index < 0 || index >= datumList.size() )
947  throw CoordinateConversionException( ErrorMessages::invalidIndex );
948  else
949  {
950  Datum* datum = datumList[index];
951 
952  *deltaX = datum->deltaX();
953  *deltaY = datum->deltaY();
954  *deltaZ = datum->deltaZ();
955  }
956 }
957 
958 
959 Accuracy* DatumLibraryImplementation::datumShiftError(
960  const long sourceIndex,
961  const long targetIndex,
962  double longitude,
963  double latitude,
964  Accuracy* sourceAccuracy,
965  Precision::Enum precision )
966 {
967  double sinlat = sin( latitude );
968  double coslat = cos( latitude );
969  double sinlon = sin( longitude );
970  double coslon = cos( longitude );
971  double sigma_delta_lat;
972  double sigma_delta_lon;
973  double sigma_delta_height;
974  double sx, sy, sz;
975  double ce90_in = -1.0;
976  double le90_in = -1.0;
977  double se90_in = -1.0;
978  double ce90_out = -1.0;
979  double le90_out = -1.0;
980  double se90_out = -1.0;
981  double circularError90 = sourceAccuracy->circularError90();
982  double linearError90 = sourceAccuracy->linearError90();
983  double sphericalError90 = sourceAccuracy->sphericalError90();
984 
985  int numDatums = datumList.size();
986 
987  if( ( sourceIndex < 0 ) || ( sourceIndex >= numDatums ) )
988  throw CoordinateConversionException( ErrorMessages::invalidIndex );
989  if( ( targetIndex < 0 ) || ( targetIndex >= numDatums ) )
990  throw CoordinateConversionException( ErrorMessages::invalidIndex );
991  if( ( latitude < ( -90 * PI_OVER_180 ) ) ||
992  ( latitude > ( 90 * PI_OVER_180 ) ) )
993  throw CoordinateConversionException( ErrorMessages::latitude );
994  if( ( longitude < ( -PI ) ) || ( longitude > TWO_PI ) )
995  throw CoordinateConversionException( ErrorMessages::longitude );
996 
997  if( sourceIndex == targetIndex )
998  { /* Just copy */
999  circularError90 = circularError90;
1000  linearError90 = linearError90;
1001  sphericalError90 = sphericalError90;
1002  }
1003  else
1004  {
1005  Datum* sourceDatum = datumList[sourceIndex];
1006  Datum* targetDatum = datumList[targetIndex];
1007 
1008  /* calculate input datum errors */
1009  switch( sourceDatum->datumType() )
1010  {
1011  case DatumType::wgs84Datum:
1012  case DatumType::wgs72Datum:
1013  case DatumType::sevenParamDatum:
1014  {
1015  ce90_in = 0.0;
1016  le90_in = 0.0;
1017  se90_in = 0.0;
1018  break;
1019  }
1020  case DatumType::threeParamDatum:
1021  {
1022  ThreeParameterDatum* sourceThreeParameterDatum =
1023  ( ThreeParameterDatum* )sourceDatum;
1024 
1025  if( ( sourceThreeParameterDatum->sigmaX() < 0 )
1026  || ( sourceThreeParameterDatum->sigmaY() < 0 )
1027  || ( sourceThreeParameterDatum->sigmaZ() < 0 ) )
1028  {
1029  ce90_in = -1.0;
1030  le90_in = -1.0;
1031  se90_in = -1.0;
1032  }
1033  else
1034  {
1035  sx = sourceThreeParameterDatum->sigmaX() * sinlat * coslon;
1036  sy = sourceThreeParameterDatum->sigmaY() * sinlat * sinlon;
1037  sz = sourceThreeParameterDatum->sigmaZ() * coslat;
1038  sigma_delta_lat = sqrt( ( sx * sx ) + ( sy * sy ) + ( sz * sz ) );
1039 
1040  sx = sourceThreeParameterDatum->sigmaX() * sinlon;
1041  sy = sourceThreeParameterDatum->sigmaY() * coslon;
1042  sigma_delta_lon = sqrt( ( sx * sx ) + ( sy * sy ) );
1043 
1044  sx = sourceThreeParameterDatum->sigmaX() * coslat * coslon;
1045  sy = sourceThreeParameterDatum->sigmaY() * coslat * sinlon;
1046  sz = sourceThreeParameterDatum->sigmaZ() * sinlat;
1047  sigma_delta_height = sqrt( ( sx * sx ) + ( sy * sy ) + ( sz * sz ) );
1048 
1049  ce90_in = 2.146 * ( sigma_delta_lat + sigma_delta_lon ) / 2.0;
1050  le90_in = 1.6449 * sigma_delta_height;
1051  se90_in = 2.5003 *
1052  ( sourceThreeParameterDatum->sigmaX() +
1053  sourceThreeParameterDatum->sigmaY() +
1054  sourceThreeParameterDatum->sigmaZ() ) / 3.0;
1055  }
1056  break;
1057  }
1058  }
1059 
1060  /* calculate output datum errors */
1061  switch( targetDatum->datumType() )
1062  {
1063  case DatumType::wgs84Datum:
1064  case DatumType::wgs72Datum:
1065  case DatumType::sevenParamDatum:
1066  {
1067  ce90_out = 0.0;
1068  le90_out = 0.0;
1069  se90_out = 0.0;
1070  break;
1071  }
1072  case DatumType::threeParamDatum:
1073  {
1074  ThreeParameterDatum* targetThreeParameterDatum =
1075  ( ThreeParameterDatum* )targetDatum;
1076  if ((targetThreeParameterDatum->sigmaX() < 0)
1077  ||(targetThreeParameterDatum->sigmaY() < 0)
1078  ||(targetThreeParameterDatum->sigmaZ() < 0))
1079  {
1080  ce90_out = -1.0;
1081  le90_out = -1.0;
1082  se90_out = -1.0;
1083  }
1084  else
1085  {
1086  sx = targetThreeParameterDatum->sigmaX() * sinlat * coslon;
1087  sy = targetThreeParameterDatum->sigmaY() * sinlat * sinlon;
1088  sz = targetThreeParameterDatum->sigmaZ() * coslat;
1089  sigma_delta_lat = sqrt((sx * sx) + (sy * sy) + (sz * sz));
1090 
1091  sx = targetThreeParameterDatum->sigmaX() * sinlon;
1092  sy = targetThreeParameterDatum->sigmaY() * coslon;
1093  sigma_delta_lon = sqrt( ( sx * sx ) + ( sy * sy ) );
1094 
1095  sx = targetThreeParameterDatum->sigmaX() * coslat * coslon;
1096  sy = targetThreeParameterDatum->sigmaY() * coslat * sinlon;
1097  sz = targetThreeParameterDatum->sigmaZ() * sinlat;
1098  sigma_delta_height = sqrt( ( sx * sx ) + ( sy * sy ) + ( sz * sz ) );
1099 
1100  ce90_out = 2.146 * ( sigma_delta_lat + sigma_delta_lon ) / 2.0;
1101  le90_out = 1.6449 * sigma_delta_height;
1102  se90_out = 2.5003 *
1103  ( targetThreeParameterDatum->sigmaX() +
1104  targetThreeParameterDatum->sigmaY() +
1105  targetThreeParameterDatum->sigmaZ() ) / 3.0;
1106  }
1107  break;
1108  }
1109  }
1110 
1111  /* combine errors */
1112  if( ( circularError90 < 0.0 ) || ( ce90_in < 0.0 ) || ( ce90_out < 0.0 ) )
1113  {
1114  circularError90 = -1.0;
1115  linearError90 = -1.0;
1116  sphericalError90 = -1.0;
1117  }
1118  else
1119  {
1120  circularError90 = sqrt( ( circularError90 * circularError90 ) +
1121  ( ce90_in * ce90_in ) + ( ce90_out * ce90_out ) );
1122  if( circularError90 < 1.0 )
1123  {
1124  circularError90 = 1.0;
1125  }
1126 
1127  if( ( linearError90 < 0.0 ) || ( le90_in < 0.0 ) || ( le90_out < 0.0 ) )
1128  {
1129  linearError90 = -1.0;
1130  sphericalError90 = -1.0;
1131  }
1132  else
1133  {
1134  linearError90 = sqrt( ( linearError90 * linearError90 ) +
1135  ( le90_in * le90_in ) + ( le90_out * le90_out ) );
1136  if( linearError90 < 1.0 )
1137  {
1138  linearError90 = 1.0;
1139  }
1140 
1141  if( ( sphericalError90 < 0.0 )
1142  || ( se90_in < 0.0 )
1143  || ( se90_out < 0.0 ) )
1144  sphericalError90 = -1.0;
1145  else
1146  {
1147  sphericalError90 = sqrt(
1148  ( sphericalError90 * sphericalError90 ) +
1149  ( se90_in * se90_in ) + ( se90_out * se90_out ) );
1150  if( sphericalError90 < 1.0 )
1151  {
1152  sphericalError90 = 1.0;
1153  }
1154  }
1155  }
1156  }
1157  }
1158 
1159  // Correct for limited precision of input/output coordinate.
1160 
1161  // sigma for uniform distribution due to rounding.
1162  double sigma = Precision::toMeters( precision ) / sqrt( 12.0 );
1163 
1164  // For linear error
1165  if( linearError90 > 0.0 )
1166  {
1167  double lePrec = 1.6449 * sigma;
1168  linearError90 = sqrt(
1169  linearError90 * linearError90 + lePrec * lePrec);
1170  }
1171 
1172  // For circular error
1173  if( circularError90 > 0.0 )
1174  {
1175  double cePrec = 2.146 * sigma;
1176  circularError90 = sqrt(
1177  circularError90 * circularError90 + cePrec * cePrec);
1178  }
1179 
1180  // For spherical error
1181  if( sphericalError90 > 0.0 )
1182  {
1183  // Assume sigma in all directions
1184  double sePrec = 2.5003 * sigma;
1185  sphericalError90 = sqrt(
1186  sphericalError90 * sphericalError90 + sePrec * sePrec );
1187  }
1188 
1189  return new Accuracy( circularError90, linearError90, sphericalError90 );
1190 }
1191 
1192 
1193 void DatumLibraryImplementation::datumUserDefined(
1194  const long index, long *result )
1195 {
1196 /*
1197  * The function datumUserDefined checks whether or not the specified datum is
1198  * user defined. It returns 1 if the datum is user defined, and returns
1199  * 0 otherwise. If index is valid DATUM_NO_ERROR is returned, otherwise
1200  * DATUM_INVALID_INDEX_ERROR is returned.
1201  *
1202  * index : Index of a given datum in the datum table (input)
1203  * result : Indicates whether specified datum is user defined (1)
1204  * or not (0) (output)
1205  */
1206 
1207  *result = false;
1208 
1209  if( index < 0 || index >= datumList.size() )
1210  throw CoordinateConversionException( ErrorMessages::invalidIndex );
1211  else
1212  {
1213  Datum* datum = datumList[index];
1214 
1215  if( datum->datumType() == DatumType::threeParamDatum )
1216  {
1217  ThreeParameterDatum* threeParameterDatum = ( ThreeParameterDatum* )datum;
1218 
1219  if( threeParameterDatum->userDefined() )
1220  *result = true;
1221  }
1222  else if( datum->datumType() == DatumType::sevenParamDatum )
1223  {
1224  SevenParameterDatum* sevenParamDatum = ( SevenParameterDatum* )datum;
1225 
1226  if( sevenParamDatum->userDefined() )
1227  *result = true;
1228  }
1229  else
1230  throw CoordinateConversionException( ErrorMessages::invalidIndex );
1231  }
1232 }
1233 
1234 
1235 bool DatumLibraryImplementation::datumUsesEllipsoid( const char *ellipsoidCode )
1236 {
1237 /*
1238  * The function datumUsesEllipsoid returns 1 if the ellipsoid is in use by a
1239  * user defined datum. Otherwise, 0 is returned.
1240  *
1241  * ellipsoidCode : The ellipsoid code being searched for. (input)
1242  */
1243 
1244  char temp_code[DATUM_CODE_LENGTH];
1245  long length;
1246  long pos = 0;
1247  long i = 0;
1248  bool ellipsoid_in_use = false;
1249 
1250  length = strlen( ellipsoidCode );
1251  if( length <= ( ELLIPSOID_CODE_LENGTH-1 ) )
1252  {
1253  strcpy( temp_code, ellipsoidCode );
1254 
1255  /* Convert to upper case */
1256  for( i=0;i<length;i++ )
1257  temp_code[i] = ( char )toupper( temp_code[i] );
1258 
1259  /* Strip blank spaces */
1260  while( pos < length )
1261  {
1262  if( isspace( temp_code[pos] ) )
1263  {
1264  for( i=pos; i<=length; i++ )
1265  temp_code[i] = temp_code[i+1];
1266  length -= 1;
1267  }
1268  else
1269  pos += 1;
1270  }
1271 
1272  /* Search for code */
1273  i = 0;
1274  int numDatums = datumList.size();
1275  while( ( i < numDatums ) && ( !ellipsoid_in_use ) )
1276  {
1277  if( strcmp( temp_code, datumList[i]->ellipsoidCode() ) == 0 )
1278  ellipsoid_in_use = true;
1279  i++;
1280  }
1281  }
1282 
1283  return ellipsoid_in_use;
1284 }
1285 
1286 
1287 void DatumLibraryImplementation::datumValidRectangle(
1288  const long index,
1289  double *westLongitude,
1290  double *eastLongitude,
1291  double *southLatitude,
1292  double *northLatitude )
1293 {
1294 /*
1295  * The function datumValidRectangle returns the edges of the validity
1296  * rectangle for the datum referenced by index.
1297  *
1298  * index : The index of a given datum in the datum table (input)
1299  * westLongitude : Western edge of validity rectangle in radians (output)
1300  * eastLongitude : Eastern edge of validity rectangle in radians (output)
1301  * southLatitude : Southern edge of validity rectangle in radians (output)
1302  * northLatitude : Northern edge of validity rectangle in radians (output)
1303  *
1304  */
1305 
1306  if( index < 0 && index >= datumList.size() )
1307  throw CoordinateConversionException( ErrorMessages::invalidIndex );
1308  else
1309  {
1310  *westLongitude = datumList[index]->westLongitude();
1311  *eastLongitude = datumList[index]->eastLongitude();
1312  *southLatitude = datumList[index]->southLatitude();
1313  *northLatitude = datumList[index]->northLatitude();
1314  }
1315 }
1316 
1317 
1318 CartesianCoordinates* DatumLibraryImplementation::geocentricDatumShift(
1319  const long sourceIndex,
1320  const double sourceX,
1321  const double sourceY,
1322  const double sourceZ,
1323  const long targetIndex )
1324 {
1325 /*
1326  * The function geocentricDatumShift shifts a geocentric coordinate (X, Y, Z in meters) relative
1327  * to the source datum to geocentric coordinate (X, Y, Z in meters) relative
1328  * to the destination datum.
1329  *
1330  * sourceIndex : Index of source datum (input)
1331  * sourceX : X coordinate relative to source datum (input)
1332  * sourceY : Y coordinate relative to source datum (input)
1333  * sourceZ : Z coordinate relative to source datum (input)
1334  * targetIndex : Index of destination datum (input)
1335  * targetX : X coordinate relative to destination datum (output)
1336  * targetY : Y coordinate relative to destination datum (output)
1337  * targetZ : Z coordinate relative to destination datum (output)
1338  *
1339  */
1340 
1341  int numDatums = datumList.size();
1342 
1343  if( ( sourceIndex < 0 ) || ( sourceIndex >= numDatums ) )
1344  throw CoordinateConversionException( ErrorMessages::invalidIndex );
1345  if( ( targetIndex < 0 ) || ( targetIndex >= numDatums ) )
1346  throw CoordinateConversionException( ErrorMessages::invalidIndex );
1347 
1348  if( sourceIndex == targetIndex )
1349  {
1350  return new CartesianCoordinates(
1351  CoordinateType::geocentric, sourceX, sourceY, sourceZ );
1352  }
1353  else
1354  {
1355  CartesianCoordinates* wgs84CartesianCoordinates = geocentricShiftToWGS84(
1356  sourceIndex, sourceX, sourceY, sourceZ );
1357 
1358  CartesianCoordinates* targetCartesianCoordinates = geocentricShiftFromWGS84(
1359  wgs84CartesianCoordinates->x(), wgs84CartesianCoordinates->y(),
1360  wgs84CartesianCoordinates->z(), targetIndex );
1361 
1362  delete wgs84CartesianCoordinates;
1363 
1364  return targetCartesianCoordinates;
1365  }
1366 }
1367 
1368 
1369 CartesianCoordinates* DatumLibraryImplementation::geocentricShiftFromWGS84(
1370  const double WGS84X,
1371  const double WGS84Y,
1372  const double WGS84Z,
1373  const long targetIndex )
1374 {
1375 /*
1376  * The function geocentricShiftFromWGS84 shifts a geocentric coordinate (X, Y, Z in meters) relative
1377  * to WGS84 to a geocentric coordinate (X, Y, Z in meters) relative to the
1378  * local datum referenced by index.
1379  *
1380  * WGS84X : X coordinate relative to WGS84 (input)
1381  * WGS84Y : Y coordinate relative to WGS84 (input)
1382  * WGS84Z : Z coordinate relative to WGS84 (input)
1383  * targetIndex : Index of destination datum (input)
1384  * targetX : X coordinate relative to the destination datum (output)
1385  * targetY : Y coordinate relative to the destination datum (output)
1386  * targetZ : Z coordinate relative to the destination datum (output)
1387  */
1388 
1389  int numDatums = datumList.size();
1390 
1391  if( ( targetIndex < 0 ) || ( targetIndex >= numDatums ) )
1392  throw CoordinateConversionException( ErrorMessages::invalidIndex );
1393 
1394  Datum* localDatum = datumList[targetIndex];
1395  switch( localDatum->datumType() )
1396  {
1397  case DatumType::wgs72Datum:
1398  {
1399  CartesianCoordinates* wgs72CartesianCoordinates =
1400  geocentricShiftWGS84ToWGS72( WGS84X, WGS84Y, WGS84Z );
1401 
1402  return wgs72CartesianCoordinates;
1403  }
1404  case DatumType::wgs84Datum:
1405  {
1406  return new CartesianCoordinates(
1407  CoordinateType::geocentric, WGS84X, WGS84Y, WGS84Z );
1408  }
1409  case DatumType::sevenParamDatum:
1410  {
1411  SevenParameterDatum* sevenParameterDatum =
1412  ( SevenParameterDatum* )localDatum;
1413 
1414  double targetX = WGS84X - sevenParameterDatum->deltaX() - sevenParameterDatum->rotationZ() * WGS84Y
1415  + sevenParameterDatum->rotationY() * WGS84Z - sevenParameterDatum->scaleFactor() * WGS84X;
1416 
1417  double targetY = WGS84Y - sevenParameterDatum->deltaY() + sevenParameterDatum->rotationZ() * WGS84X
1418  - sevenParameterDatum->rotationX() * WGS84Z - sevenParameterDatum->scaleFactor() * WGS84Y;
1419 
1420  double targetZ = WGS84Z - sevenParameterDatum->deltaZ() - sevenParameterDatum->rotationY() * WGS84X
1421  + sevenParameterDatum->rotationX() * WGS84Y - sevenParameterDatum->scaleFactor() * WGS84Z;
1422 
1423  return new CartesianCoordinates( CoordinateType::geocentric, targetX, targetY, targetZ );
1424  }
1425  case DatumType::threeParamDatum:
1426  {
1427  ThreeParameterDatum* threeParameterDatum = ( ThreeParameterDatum* )localDatum;
1428 
1429  double targetX = WGS84X - threeParameterDatum->deltaX();
1430  double targetY = WGS84Y - threeParameterDatum->deltaY();
1431  double targetZ = WGS84Z - threeParameterDatum->deltaZ();
1432 
1433  return new CartesianCoordinates( CoordinateType::geocentric, targetX, targetY, targetZ );
1434  }
1435  default:
1436  throw CoordinateConversionException( ErrorMessages::datumType );
1437  }
1438 }
1439 
1440 
1441 CartesianCoordinates* DatumLibraryImplementation::geocentricShiftToWGS84(
1442  const long sourceIndex,
1443  const double sourceX,
1444  const double sourceY,
1445  const double sourceZ )
1446 {
1447 /*
1448  * The function geocentricShiftToWGS84 shifts a geocentric coordinate (X, Y, Z in meters) relative
1449  * to the datum referenced by index to a geocentric coordinate (X, Y, Z in
1450  * meters) relative to WGS84.
1451  *
1452  * sourceIndex : Index of source datum (input)
1453  * sourceX : X coordinate relative to the source datum (input)
1454  * sourceY : Y coordinate relative to the source datum (input)
1455  * sourceZ : Z coordinate relative to the source datum (input)
1456  * WGS84X : X coordinate relative to WGS84 (output)
1457  * WGS84Y : Y coordinate relative to WGS84 (output)
1458  * WGS84Z : Z coordinate relative to WGS84 (output)
1459  */
1460 
1461  int numDatums = datumList.size();
1462 
1463  if( ( sourceIndex < 0 ) || (sourceIndex > numDatums ) )
1464  throw CoordinateConversionException( ErrorMessages::invalidIndex );
1465 
1466  Datum* localDatum = datumList[sourceIndex];
1467  switch( localDatum->datumType() )
1468  {
1469  case DatumType::wgs72Datum:
1470  {
1471  CartesianCoordinates* wgs84CartesianCoordinates =
1472  geocentricShiftWGS72ToWGS84( sourceX, sourceY, sourceZ );
1473 
1474  return wgs84CartesianCoordinates;
1475  }
1476  case DatumType::wgs84Datum:
1477  {
1478  return new CartesianCoordinates( CoordinateType::geocentric, sourceX, sourceY, sourceZ );
1479  }
1480  case DatumType::sevenParamDatum:
1481  {
1482  SevenParameterDatum* sevenParameterDatum = ( SevenParameterDatum* )localDatum;
1483 
1484  double wgs84X = sourceX + sevenParameterDatum->deltaX() + sevenParameterDatum->rotationZ() * sourceY
1485  - sevenParameterDatum->rotationY() * sourceZ + sevenParameterDatum->scaleFactor() * sourceX;
1486 
1487  double wgs84Y = sourceY + sevenParameterDatum->deltaY() - sevenParameterDatum->rotationZ() * sourceX
1488  + sevenParameterDatum->rotationX() * sourceZ + sevenParameterDatum->scaleFactor() * sourceY;
1489 
1490  double wgs84Z = sourceZ + sevenParameterDatum->deltaZ() + sevenParameterDatum->rotationY() * sourceX
1491  - sevenParameterDatum->rotationX() * sourceY + sevenParameterDatum->scaleFactor() * sourceZ;
1492 
1493  return new CartesianCoordinates( CoordinateType::geocentric, wgs84X, wgs84Y, wgs84Z );
1494  }
1495  case DatumType::threeParamDatum:
1496  {
1497  ThreeParameterDatum* threeParameterDatum = ( ThreeParameterDatum* )localDatum;
1498 
1499  double wgs84X = sourceX + threeParameterDatum->deltaX();
1500  double wgs84Y = sourceY + threeParameterDatum->deltaY();
1501  double wgs84Z = sourceZ + threeParameterDatum->deltaZ();
1502 
1503  return new CartesianCoordinates( CoordinateType::geocentric, wgs84X, wgs84Y, wgs84Z );
1504  }
1505  default:
1506  throw CoordinateConversionException( ErrorMessages::datumType );
1507  }
1508 }
1509 
1510 
1511 GeodeticCoordinates* DatumLibraryImplementation::geodeticDatumShift(
1512  const long sourceIndex,
1513  const GeodeticCoordinates* sourceCoordinates,
1514  const long targetIndex )
1515 {
1516 /*
1517  * The function geodeticDatumShift shifts geodetic coordinates (latitude, longitude in radians
1518  * and height in meters) relative to the source datum to geodetic coordinates
1519  * (latitude, longitude in radians and height in meters) relative to the
1520  * destination datum.
1521  *
1522  * sourceIndex : Index of source datum (input)
1523  * sourceLongitude : Longitude in radians relative to source datum (input)
1524  * sourceLatitude : Latitude in radians relative to source datum (input)
1525  * sourceHeight : Height in meters relative to source datum (input)
1526  * targetIndex : Index of destination datum (input)
1527  * targetLongitude : Longitude in radians relative to destination datum(output)
1528  * targetLatitude : Latitude in radians relative to destination datum (output)
1529  * targetHeight : Height in meters relative to destination datum (output)
1530  */
1531 
1532  long E_Index;
1533  double a;
1534  double f;
1535 
1536  double sourceLongitude = sourceCoordinates->longitude();
1537  double sourceLatitude = sourceCoordinates->latitude();
1538  double sourceHeight = sourceCoordinates->height();
1539 
1540  int numDatums = datumList.size();
1541 
1542  if( sourceIndex < 0 || sourceIndex >= numDatums )
1543  throw CoordinateConversionException( ErrorMessages::invalidIndex );
1544  if( targetIndex < 0 || targetIndex >= numDatums )
1545  throw CoordinateConversionException( ErrorMessages::invalidIndex );
1546 
1547  if(( sourceLatitude < ( -90 * PI_OVER_180 ) ) ||
1548  ( sourceLatitude > ( 90 * PI_OVER_180 ) ) )
1549  throw CoordinateConversionException( ErrorMessages::latitude );
1550  if( ( sourceLongitude < ( -PI )) || ( sourceLongitude > TWO_PI ) )
1551  throw CoordinateConversionException( ErrorMessages::longitude );
1552 
1553  Datum* sourceDatum = datumList[sourceIndex];
1554  Datum* targetDatum = datumList[targetIndex];
1555 
1556  if ( sourceIndex == targetIndex )
1557  { /* Just copy */
1558  return new GeodeticCoordinates(
1559  CoordinateType::geodetic, sourceLatitude, sourceLongitude, sourceHeight);
1560  }
1561  else
1562  {
1563  if( _ellipsoidLibraryImplementation )
1564  {
1565  if( sourceDatum->datumType() == DatumType::sevenParamDatum )
1566  {
1567  _ellipsoidLibraryImplementation->ellipsoidIndex(
1568  sourceDatum->ellipsoidCode(), &E_Index );
1569  _ellipsoidLibraryImplementation->ellipsoidParameters( E_Index, &a, &f );
1570 
1571  Geocentric geocentricFromGeodetic( a, f );
1572 
1573  CartesianCoordinates* sourceCartesianCoordinates =
1574  geocentricFromGeodetic.convertFromGeodetic( sourceCoordinates );
1575 
1576  if( targetDatum->datumType() == DatumType::sevenParamDatum )
1577  { /* Use 3-step method for both stages */
1578  CartesianCoordinates* targetCartesianCoordinates = geocentricDatumShift(
1579  sourceIndex,
1580  sourceCartesianCoordinates->x(),
1581  sourceCartesianCoordinates->y(),
1582  sourceCartesianCoordinates->z(), targetIndex );
1583 
1584  _ellipsoidLibraryImplementation->ellipsoidIndex(
1585  targetDatum->ellipsoidCode(), &E_Index );
1586  _ellipsoidLibraryImplementation->ellipsoidParameters( E_Index, &a, &f );
1587 
1588  Geocentric geocentricToGeodetic( a, f );
1589  GeodeticCoordinates* targetGeodeticCoordinates =
1590  geocentricToGeodetic.convertToGeodetic( targetCartesianCoordinates );
1591 
1592  delete targetCartesianCoordinates;
1593  delete sourceCartesianCoordinates;
1594 
1595  return targetGeodeticCoordinates;
1596  }
1597  else
1598  { /* Use 3-step method for 1st stage, Molodensky if possible for 2nd stage */
1599  CartesianCoordinates* wgs84CartesianCoordinates = geocentricShiftToWGS84(
1600  sourceIndex, sourceCartesianCoordinates->x(),
1601  sourceCartesianCoordinates->y(), sourceCartesianCoordinates->z() );
1602 
1603  long wgs84EllipsoidIndex;
1604  _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
1605  _ellipsoidLibraryImplementation->ellipsoidParameters(
1606  wgs84EllipsoidIndex, &a, &f );
1607 
1608  Geocentric geocentricToGeodetic( a, f );
1609  GeodeticCoordinates* wgs84GeodeticCoordinates =
1610  geocentricToGeodetic.convertToGeodetic( wgs84CartesianCoordinates );
1611 
1612  GeodeticCoordinates* targetGeodeticCoordinates =
1613  geodeticShiftFromWGS84( wgs84GeodeticCoordinates, targetIndex );
1614 
1615  delete wgs84GeodeticCoordinates;
1616  delete wgs84CartesianCoordinates;
1617 
1618  return targetGeodeticCoordinates;
1619  }
1620  }
1621  else if( targetDatum->datumType() == DatumType::sevenParamDatum )
1622  { /* Use Molodensky if possible for 1st stage, 3-step method for 2nd stage */
1623  GeodeticCoordinates* wgs84GeodeticCoordinates = geodeticShiftToWGS84(
1624  sourceIndex, sourceCoordinates );
1625 
1626  long wgs84EllipsoidIndex;
1627  _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
1628  _ellipsoidLibraryImplementation->ellipsoidParameters(
1629  wgs84EllipsoidIndex, &a, &f );
1630 
1631  Geocentric geocentricFromGeodetic( a, f );
1632  CartesianCoordinates* wgs84CartesianCoordinates =
1633  geocentricFromGeodetic.convertFromGeodetic( wgs84GeodeticCoordinates );
1634 
1635  CartesianCoordinates* targetCartesianCoordinates =
1636  geocentricShiftFromWGS84(
1637  wgs84CartesianCoordinates->x(),
1638  wgs84CartesianCoordinates->y(),
1639  wgs84CartesianCoordinates->z(), targetIndex );
1640 
1641  _ellipsoidLibraryImplementation->ellipsoidIndex(
1642  targetDatum->ellipsoidCode(), &E_Index );
1643  _ellipsoidLibraryImplementation->ellipsoidParameters( E_Index, &a, &f );
1644 
1645  Geocentric geocentricToGeodetic( a, f );
1646  GeodeticCoordinates* targetGeodeticCoordinates =
1647  geocentricToGeodetic.convertToGeodetic( targetCartesianCoordinates );
1648 
1649  delete targetCartesianCoordinates;
1650  delete wgs84CartesianCoordinates;
1651  delete wgs84GeodeticCoordinates;
1652 
1653  return targetGeodeticCoordinates;
1654  }
1655  else
1656  { /* Use Molodensky if possible for both stages */
1657  GeodeticCoordinates* wgs84GeodeticCoordinates = geodeticShiftToWGS84(
1658  sourceIndex, sourceCoordinates );
1659 
1660  GeodeticCoordinates* targetGeodeticCoordinates = geodeticShiftFromWGS84(
1661  wgs84GeodeticCoordinates, targetIndex );
1662 
1663  delete wgs84GeodeticCoordinates;
1664 
1665  return targetGeodeticCoordinates;
1666  }
1667  }
1668  else
1669  throw CoordinateConversionException( ErrorMessages::ellipse );
1670  }
1671 }
1672 
1673 
1674 GeodeticCoordinates* DatumLibraryImplementation::geodeticShiftFromWGS84(
1675  const GeodeticCoordinates* sourceCoordinates, const long targetIndex )
1676 {
1677 /*
1678  * The function geodeticShiftFromWGS84 shifts geodetic coordinates relative to WGS84
1679  * to geodetic coordinates relative to a given local datum.
1680  *
1681  * WGS84Longitude : Longitude in radians relative to WGS84 (input)
1682  * WGS84Latitude : Latitude in radians relative to WGS84 (input)
1683  * WGS84Height : Height in meters relative to WGS84 (input)
1684  * targetIndex : Index of destination datum (input)
1685  * targetLongitude : Longitude (rad) relative to destination datum (output)
1686  * targetLatitude : Latitude (rad) relative to destination datum (output)
1687  * targetHeight : Height in meters relative to destination datum (output)
1688  *
1689  */
1690 
1691  double WGS84_a; /* Semi-major axis of WGS84 ellipsoid in meters */
1692  double WGS84_f; /* Flattening of WGS84 ellisoid */
1693  double a; /* Semi-major axis of ellipsoid in meters */
1694  double da; /* Difference in semi-major axes */
1695  double f; /* Flattening of ellipsoid */
1696  double df; /* Difference in flattening */
1697  double dx;
1698  double dy;
1699  double dz;
1700  long E_Index;
1701 
1702  double WGS84Longitude = sourceCoordinates->longitude();
1703  double WGS84Latitude = sourceCoordinates->latitude();
1704  double WGS84Height = sourceCoordinates->height();
1705 
1706  if( ( targetIndex < 0) || (targetIndex >= datumList.size() ) )
1707  throw CoordinateConversionException( ErrorMessages::invalidIndex );
1708  if(( WGS84Latitude < ( -90 * PI_OVER_180 ) ) ||
1709  ( WGS84Latitude > ( 90 * PI_OVER_180 ) ) )
1710  throw CoordinateConversionException( ErrorMessages::latitude );
1711  if( ( WGS84Longitude < ( -PI ) ) || ( WGS84Longitude > TWO_PI ) )
1712  throw CoordinateConversionException( ErrorMessages::longitude );
1713 
1714  Datum* localDatum = datumList[targetIndex];
1715  switch( localDatum->datumType() )
1716  {
1717  case DatumType::wgs72Datum:
1718  {
1719  GeodeticCoordinates* targetGeodeticCoordinates = geodeticShiftWGS84ToWGS72( WGS84Longitude, WGS84Latitude, WGS84Height );
1720  return targetGeodeticCoordinates;
1721  }
1722  case DatumType::wgs84Datum:
1723  {
1724  return new GeodeticCoordinates( CoordinateType::geodetic, WGS84Longitude, WGS84Latitude, WGS84Height );
1725  }
1726  case DatumType::sevenParamDatum:
1727  case DatumType::threeParamDatum:
1728  {
1729  if( _ellipsoidLibraryImplementation )
1730  {
1731  _ellipsoidLibraryImplementation->ellipsoidIndex( localDatum->ellipsoidCode(), &E_Index );
1732  _ellipsoidLibraryImplementation->ellipsoidParameters( E_Index, &a, &f );
1733 
1734  long wgs84EllipsoidIndex;
1735  _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
1736  _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &WGS84_a, &WGS84_f );
1737 
1738  if( ( localDatum->datumType() == DatumType::sevenParamDatum ) ||
1739  ( WGS84Latitude < ( -MOLODENSKY_MAX ) ) ||
1740  ( WGS84Latitude > MOLODENSKY_MAX ) )
1741  { /* Use 3-step method */
1742  Geocentric geocentricFromGeodetic( WGS84_a, WGS84_f );
1743  CartesianCoordinates* wgs84CartesianCoordinates = geocentricFromGeodetic.convertFromGeodetic( sourceCoordinates );
1744 
1745  CartesianCoordinates* localCartesianCoordinates = geocentricShiftFromWGS84( wgs84CartesianCoordinates->x(), wgs84CartesianCoordinates->y(),
1746  wgs84CartesianCoordinates->z(), targetIndex );
1747 
1748  Geocentric geocentricToGeodetic( a, f );
1749  GeodeticCoordinates* targetGeodeticCoordinates = geocentricToGeodetic.convertToGeodetic( localCartesianCoordinates );
1750 
1751  delete localCartesianCoordinates;
1752  delete wgs84CartesianCoordinates;
1753 
1754  return targetGeodeticCoordinates;
1755  }
1756  else
1757  { /* Use Molodensky's method */
1758  da = a - WGS84_a;
1759  df = f - WGS84_f;
1760  dx = -( localDatum->deltaX() );
1761  dy = -( localDatum->deltaY() );
1762  dz = -( localDatum->deltaZ() );
1763 
1764  GeodeticCoordinates* targetGeodeticCoordinates = molodenskyShift( WGS84_a, da, WGS84_f, df, dx, dy, dz,
1765  WGS84Longitude, WGS84Latitude, WGS84Height );
1766 
1767  return targetGeodeticCoordinates;
1768  }
1769  }
1770  }
1771  default:
1772  throw CoordinateConversionException( ErrorMessages::datumType );
1773  }
1774 }
1775 
1776 
1777 GeodeticCoordinates* DatumLibraryImplementation::geodeticShiftToWGS84( const long sourceIndex, const GeodeticCoordinates* sourceCoordinates )
1778 {
1779 /*
1780  * The function geodeticShiftToWGS84 shifts geodetic coordinates relative to a given source datum
1781  * to geodetic coordinates relative to WGS84.
1782  *
1783  * sourceIndex : Index of source datum (input)
1784  * sourceLongitude : Longitude in radians relative to source datum (input)
1785  * sourceLatitude : Latitude in radians relative to source datum (input)
1786  * sourceHeight : Height in meters relative to source datum (input)
1787  * WGS84Longitude : Longitude in radians relative to WGS84 (output)
1788  * WGS84Latitude : Latitude in radians relative to WGS84 (output)
1789  * WGS84Height : Height in meters relative to WGS84 (output)
1790  *
1791  */
1792 
1793  double WGS84_a; /* Semi-major axis of WGS84 ellipsoid in meters */
1794  double WGS84_f; /* Flattening of WGS84 ellisoid */
1795  double a; /* Semi-major axis of ellipsoid in meters */
1796  double da; /* Difference in semi-major axes */
1797  double f; /* Flattening of ellipsoid */
1798  double df; /* Difference in flattening */
1799  double dx;
1800  double dy;
1801  double dz;
1802  long E_Index;
1803 
1804  double sourceLongitude = sourceCoordinates->longitude();
1805  double sourceLatitude = sourceCoordinates->latitude();
1806  double sourceHeight = sourceCoordinates->height();
1807 
1808  if( ( sourceIndex < 0 ) || ( sourceIndex >= datumList.size() ) )
1809  throw CoordinateConversionException( ErrorMessages::invalidIndex );
1810  if(( sourceLatitude < ( -90 * PI_OVER_180 ) ) ||
1811  ( sourceLatitude > ( 90 * PI_OVER_180 ) ) )
1812  throw CoordinateConversionException( ErrorMessages::latitude );
1813  if( ( sourceLongitude < ( -PI ) ) || ( sourceLongitude > TWO_PI ) )
1814  throw CoordinateConversionException( ErrorMessages::longitude );
1815 
1816  Datum* localDatum = datumList[sourceIndex];
1817  switch( localDatum->datumType() )
1818  {
1819  case DatumType::wgs72Datum:
1820  { /* Special case for WGS72 */
1821  return geodeticShiftWGS72ToWGS84( sourceLongitude, sourceLatitude, sourceHeight );
1822  }
1823  case DatumType::wgs84Datum:
1824  { /* Just copy */
1825  return new GeodeticCoordinates(CoordinateType::geodetic, sourceLongitude, sourceLatitude, sourceHeight);
1826  }
1827  case DatumType::sevenParamDatum:
1828  case DatumType::threeParamDatum:
1829  {
1830  if( _ellipsoidLibraryImplementation )
1831  {
1832  _ellipsoidLibraryImplementation->ellipsoidIndex( localDatum->ellipsoidCode(), &E_Index );
1833  _ellipsoidLibraryImplementation->ellipsoidParameters( E_Index, &a, &f );
1834 
1835  if( ( localDatum->datumType() == DatumType::sevenParamDatum ) ||
1836  ( sourceLatitude < (-MOLODENSKY_MAX ) ) ||
1837  ( sourceLatitude > MOLODENSKY_MAX ) )
1838  { /* Use 3-step method */
1839  Geocentric geocentricFromGeodetic( a, f );
1840  CartesianCoordinates* localCartesianCoordinates =
1841  geocentricFromGeodetic.convertFromGeodetic( sourceCoordinates );
1842 
1843  CartesianCoordinates* wgs84CartesianCoordinates = geocentricShiftToWGS84( sourceIndex, localCartesianCoordinates->x(), localCartesianCoordinates->y(), localCartesianCoordinates->z() );
1844 
1845  long wgs84EllipsoidIndex;
1846  _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
1847  _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &WGS84_a, &WGS84_f );
1848 
1849  Geocentric geocentricToGeodetic( WGS84_a, WGS84_f );
1850  GeodeticCoordinates* wgs84GeodeticCoordinates = geocentricToGeodetic.convertToGeodetic( wgs84CartesianCoordinates );
1851 
1852  delete wgs84CartesianCoordinates;
1853  delete localCartesianCoordinates;
1854 
1855  return wgs84GeodeticCoordinates;
1856  }
1857  else
1858  { /* Use Molodensky's method */
1859  long wgs84EllipsoidIndex;
1860  _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
1861  _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &WGS84_a, &WGS84_f );
1862 
1863  da = WGS84_a - a;
1864  df = WGS84_f - f;
1865  dx = localDatum->deltaX();
1866  dy = localDatum->deltaY();
1867  dz = localDatum->deltaZ();
1868 
1869  GeodeticCoordinates* wgs84GeodeticCoordinates = molodenskyShift( a, da, f, df, dx, dy, dz, sourceLongitude, sourceLatitude, sourceHeight );
1870 
1871  return wgs84GeodeticCoordinates;
1872  }
1873  }
1874  }
1875  default:
1876  throw CoordinateConversionException( ErrorMessages::datumType );
1877  }
1878 }
1879 
1880 
1881 void DatumLibraryImplementation::retrieveDatumType(
1882  const long index,
1883  DatumType::Enum *datumType )
1884 {
1885 /*
1886  * The function retrieveDatumType returns the type of the datum referenced by
1887  * index.
1888  *
1889  * index : The index of a given datum in the datum table. (input)
1890  * datumType : The type of datum referenced by index. (output)
1891  *
1892  */
1893 
1894  if( index < 0 || index >= datumList.size() )
1895  throw CoordinateConversionException( ErrorMessages::invalidIndex );
1896  else
1897  *datumType = datumList[index]->datumType();
1898 }
1899 
1900 
1901 void DatumLibraryImplementation::validDatum(
1902  const long index,
1903  double longitude,
1904  double latitude,
1905  long *result )
1906 {
1907 /*
1908  * The function validDatum checks whether or not the specified location is within the
1909  * validity rectangle for the specified datum. It returns zero if the specified
1910  * location is NOT within the validity rectangle, and returns 1 otherwise.
1911  *
1912  * index : The index of a given datum in the datum table (input)
1913  * latitude : Latitude of the location to be checked in radians (input)
1914  * longitude : Longitude of the location to be checked in radians (input)
1915  * result : Indicates whether location is inside (1) or outside (0)
1916  * of the validity rectangle of the specified datum (output)
1917  */
1918  *result = 0;
1919 
1920  if( ( index < 0 ) || ( index >= datumList.size() ) )
1921  throw CoordinateConversionException( ErrorMessages::invalidIndex );
1922  if( ( latitude < MIN_LAT ) || ( latitude > MAX_LAT ) )
1923  throw CoordinateConversionException( ErrorMessages::latitude );
1924  if( ( longitude < MIN_LON ) || ( longitude > MAX_LON ) )
1925  throw CoordinateConversionException( ErrorMessages::longitude );
1926 
1927  Datum* datum = datumList[index];
1928 
1929  double west_longitude = datum->westLongitude();
1930  double east_longitude = datum->eastLongitude();
1931 
1932  /* The west and east longitudes may be in the range 0 to 360
1933  or -180 to 180, longitude is always in the -180 to 180 range
1934  Figure out which range west and east longitudes are in.
1935  If west and east are in the -180 to 180 range and west > east, put them in the 0 to 360 range and adjust longitude if
1936  necessary.
1937  If west and east are in the 0 to 360 range and west > east, put them in the -180 to 180 range. If west < east, adjust
1938  longitude to the 0 to 360 range
1939  */
1940  if( ( west_longitude < 0 ) || ( east_longitude < 0 ) )
1941  {
1942  if( west_longitude > east_longitude )
1943  {
1944  if( west_longitude < 0 )
1945  west_longitude += TWO_PI;
1946  if( east_longitude < 0 )
1947  east_longitude += TWO_PI;
1948  if( longitude < 0 )
1949  longitude += TWO_PI;
1950  }
1951  }
1952  else if( ( west_longitude > PI ) || ( east_longitude > PI ) )
1953  {
1954  if( west_longitude > east_longitude )
1955  {
1956  if( west_longitude > PI )
1957  west_longitude -= TWO_PI;
1958  if( east_longitude > PI )
1959  east_longitude -= TWO_PI;
1960  }
1961  else
1962  {
1963  if( longitude < 0 )
1964  longitude += TWO_PI;
1965  }
1966  }
1967 
1968  if( ( datum->southLatitude() <= latitude ) &&
1969  ( latitude <= datum->northLatitude() ) &&
1970  ( west_longitude <= longitude ) &&
1971  ( longitude <= east_longitude ) )
1972  {
1973  *result = 1;
1974  }
1975  else
1976  {
1977  *result = 0;
1978  }
1979 }
1980 
1981 
1982 void DatumLibraryImplementation::setEllipsoidLibraryImplementation(
1983  EllipsoidLibraryImplementation* __ellipsoidLibraryImplementation )
1984 {
1985 /*
1986  * The function setEllipsoidLibraryImplementation sets the ellipsoid library information
1987  * which is needed to create datums and calculate datum shifts.
1988  *
1989  * __ellipsoidLibraryImplementation : Ellipsoid library implementation(input)
1990  *
1991  */
1992 
1993  _ellipsoidLibraryImplementation = __ellipsoidLibraryImplementation;
1994 }
1995 
1996 
1997 /************************************************************************/
1998 /* PRIVATE FUNCTIONS
1999  *
2000  */
2001 
2002 void DatumLibraryImplementation::loadDatums()
2003 {
2004 /*
2005  * The function loadDatums creates the datum table from two external
2006  * files. If an error occurs, the initialization stops and an error code is
2007  * returned. This function must be called before any of the other functions
2008  * in this component.
2009  */
2010 
2011  long index = 0, i = 0;
2012  char *PathName = NULL;
2013  char* FileName7 = 0;
2014  FILE *fp_7param = NULL;
2015  FILE *fp_3param = NULL;
2016  char* FileName3 = 0;
2017  const int file_name_length = 23;
2018 
2019  CCSThreadLock lock(&mutex);
2020 
2021  /* Check the environment for a user provided path, else current directory; */
2022  /* Build a File Name, including specified or default path: */
2023 
2024 #ifdef NDK_BUILD
2025  PathName = "/data/data/com.baesystems.msp.geotrans/lib/";
2026  FileName7 = new char[ 80 ];
2027  strcpy( FileName7, PathName );
2028  strcat( FileName7, "lib7paramdat.so" );
2029 #else
2030  PathName = getenv( "MSPCCS_DATA" );
2031  if (PathName != NULL)
2032  {
2033  FileName7 = new char[ strlen( PathName ) + 13 ];
2034  strcpy( FileName7, PathName );
2035  strcat( FileName7, "/" );
2036  }
2037  else
2038  {
2039  FileName7 = new char[ file_name_length ];
2040  strcpy( FileName7, "../../data/" );
2041  }
2042  strcat( FileName7, "7_param.dat" );
2043 #endif
2044 
2045  /* Open the File READONLY, or Return Error Condition: */
2046 
2047  if (( fp_7param = fopen( FileName7, "r" ) ) == NULL)
2048  {
2049  delete [] FileName7;
2050  FileName7 = 0;
2051 
2052  char message[256] = "";
2053  if (NULL == PathName)
2054  {
2055  strcpy( message, "Environment variable undefined: MSPCCS_DATA." );
2056  }
2057  else
2058  {
2059  strcpy( message, ErrorMessages::datumFileOpenError );
2060  strcat( message, ": 7_param.dat\n" );
2061  }
2062  throw CoordinateConversionException( message );
2063  }
2064 
2065  /* Open the File READONLY, or Return Error Condition: */
2066 
2067  /* WGS84 datum entry */
2068  datumList.push_back( new Datum(
2069  index, "WGE", "WE", "World Geodetic System 1984", DatumType::wgs84Datum,
2070  0.0, 0.0, 0.0, -PI, +PI, -PI / 2.0, +PI / 2.0, false ) );
2071  index ++;
2072 
2073 
2074  /* WGS72 datum entry */
2075  datumList.push_back( new Datum(
2076  index, "WGC", "WD", "World Geodetic System 1972", DatumType::wgs72Datum,
2077  0.0, 0.0, 0.0, -PI, +PI, -PI / 2.0, +PI / 2.0, false ) );
2078 
2079  index ++;
2080 
2081  datum7ParamCount = 0;
2082  /* build 7-parameter datum table entries */
2083  while ( !feof(fp_7param ) )
2084  {
2085  char code[DATUM_CODE_LENGTH];
2086 
2087  bool userDefined = false;
2088 
2089  if (fscanf(fp_7param, "%s ", code) <= 0)
2090  {
2091  fclose( fp_7param );
2092  throw CoordinateConversionException( ErrorMessages::datumFileParseError );
2093  }
2094  else
2095  {
2096  if( code[0] == '*' )
2097  {
2098  userDefined = true;
2099  for( int i = 0; i < DATUM_CODE_LENGTH; i++ )
2100  code[i] = code[i+1];
2101  }
2102  }
2103 
2104  char name[DATUM_NAME_LENGTH];
2105  if (fscanf(fp_7param, "\"%32[^\"]\"", name) <= 0)
2106  name[0] = '\0';
2107 
2108  char ellipsoidCode[ELLIPSOID_CODE_LENGTH];
2109  double deltaX;
2110  double deltaY;
2111  double deltaZ;
2112  double rotationX;
2113  double rotationY;
2114  double rotationZ;
2115  double scaleFactor;
2116 
2117  if( fscanf( fp_7param, " %s %lf %lf %lf %lf %lf %lf %lf ",
2118  ellipsoidCode, &deltaX, &deltaY, &deltaZ,
2119  &rotationX, &rotationY, &rotationZ, &scaleFactor ) <= 0 )
2120  {
2121  fclose( fp_7param );
2122  throw CoordinateConversionException( ErrorMessages::datumFileParseError );
2123  }
2124  else
2125  { /* convert from degrees to radians */
2126 
2127  rotationX /= SECONDS_PER_RADIAN;
2128  rotationY /= SECONDS_PER_RADIAN;
2129  rotationZ /= SECONDS_PER_RADIAN;
2130 
2131  datumList.push_back( new SevenParameterDatum(
2132  index, code, ellipsoidCode, name, DatumType::sevenParamDatum,
2133  deltaX, deltaY, deltaZ, -PI, +PI, -PI / 2.0, +PI / 2.0,
2134  rotationX, rotationY, rotationZ, scaleFactor, userDefined ) );
2135  }
2136  index++;
2137  datum7ParamCount++;
2138  }
2139  fclose( fp_7param );
2140 
2141 #ifdef NDK_BUILD
2142  PathName = "/data/data/com.baesystems.msp.geotrans/lib/";
2143  FileName3 = new char[ 80 ];
2144  strcpy( FileName3, PathName );
2145  strcat( FileName3, "lib3paramdat.so" );
2146 #else
2147  if (PathName != NULL)
2148  {
2149  FileName3 = new char[ strlen( PathName ) + 13 ];
2150  strcpy( FileName3, PathName );
2151  strcat( FileName3, "/" );
2152  }
2153  else
2154  {
2155  FileName3 = new char[ file_name_length ];
2156  strcpy( FileName3, "../../data/" );
2157  }
2158  strcat( FileName3, "3_param.dat" );
2159 #endif
2160 
2161  if (( fp_3param = fopen( FileName3, "r" ) ) == NULL)
2162  {
2163  delete [] FileName7;
2164  FileName7 = 0;
2165  delete [] FileName3;
2166  FileName3 = 0;
2167 
2168  char message[256] = "";
2169  strcpy( message, ErrorMessages::datumFileOpenError );
2170  strcat( message, ": 3_param.dat\n" );
2171  throw CoordinateConversionException( message );
2172  }
2173 
2174  datum3ParamCount = 0;
2175 
2176  /* build 3-parameter datum table entries */
2177  while( !feof( fp_3param ) )
2178  {
2179  char code[DATUM_CODE_LENGTH];
2180 
2181  bool userDefined = false;
2182 
2183  if( fscanf( fp_3param, "%s ", code ) <= 0 )
2184  {
2185  fclose( fp_3param );
2186  throw CoordinateConversionException( ErrorMessages::datumFileParseError );
2187  }
2188  else
2189  {
2190  if( code[0] == '*' )
2191  {
2192  userDefined = true;
2193  for( int i = 0; i < DATUM_CODE_LENGTH; i++ )
2194  code[i] = code[i+1];
2195  }
2196  }
2197 
2198  char name[DATUM_NAME_LENGTH];
2199  if( fscanf( fp_3param, "\"%32[^\"]\"", name) <= 0 )
2200  name[0] = '\0';
2201 
2202  char ellipsoidCode[ELLIPSOID_CODE_LENGTH];
2203  double deltaX;
2204  double deltaY;
2205  double deltaZ;
2206  double sigmaX;
2207  double sigmaY;
2208  double sigmaZ;
2209  double eastLongitude;
2210  double westLongitude;
2211  double northLatitude;
2212  double southLatitude;
2213 
2214  if( fscanf( fp_3param, " %s %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf ",
2215  ellipsoidCode, &deltaX, &sigmaX, &deltaY, &sigmaY, &deltaZ, &sigmaZ,
2216  &southLatitude, &northLatitude, &westLongitude, &eastLongitude ) <= 0 )
2217  {
2218  fclose( fp_3param );
2219  throw CoordinateConversionException( ErrorMessages::datumFileParseError );
2220  }
2221  else
2222  {
2223  southLatitude *= PI_OVER_180;
2224  northLatitude *= PI_OVER_180;
2225  westLongitude *= PI_OVER_180;
2226  eastLongitude *= PI_OVER_180;
2227 
2228  datumList.push_back( new ThreeParameterDatum(
2229  index, code, ellipsoidCode, name, DatumType::threeParamDatum,
2230  deltaX, deltaY, deltaZ, westLongitude, eastLongitude,
2231  southLatitude, northLatitude, sigmaX, sigmaY, sigmaZ,
2232  userDefined ) );
2233  }
2234 
2235  index++;
2236  datum3ParamCount++;
2237  }
2238  fclose( fp_3param );
2239 
2240  delete [] FileName7;
2241  FileName7 = 0;
2242  delete [] FileName3;
2243  FileName3 = 0;
2244 }
2245 
2246 
2247 void DatumLibraryImplementation::write3ParamFile()
2248 {
2249 /*
2250  * The function write3ParamFile writes the 3 parameter datums in the datum list
2251  * to the 3_param.dat file.
2252  */
2253 
2254  char datum_name[DATUM_NAME_LENGTH+2];
2255  char *PathName = NULL;
2256  char FileName[FILENAME_LENGTH];
2257  FILE *fp_3param = NULL;
2258 
2259  CCSThreadLock lock(&mutex);
2260 
2261  /*output updated 3-parameter datum table*/
2262  PathName = getenv( "MSPCCS_DATA" );
2263  if (PathName != NULL)
2264  {
2265  strcpy( FileName, PathName );
2266  strcat( FileName, "/" );
2267  }
2268  else
2269  {
2270  strcpy( FileName, "../../data/" );
2271  }
2272  strcat( FileName, "3_param.dat" );
2273 
2274  if ((fp_3param = fopen(FileName, "w")) == NULL)
2275  { /* fatal error */
2276  char message[256] = "";
2277  if (NULL == PathName)
2278  {
2279  strcpy( message, "Environment variable undefined: MSPCCS_DATA." );
2280  }
2281  else
2282  {
2283  strcpy( message, ErrorMessages::datumFileOpenError );
2284  strcat( message, ": 3_param.dat\n" );
2285  }
2286  throw CoordinateConversionException( message );
2287  }
2288 
2289  /* write file */
2290  long index = MAX_WGS + datum7ParamCount;
2291  int size = datumList.size();
2292  while( index < size )
2293  {
2294  ThreeParameterDatum* _3parameterDatum = ( ThreeParameterDatum* )datumList[index];
2295  if( _3parameterDatum )
2296  {
2297  strcpy( datum_name, "\"" );
2298  strcat( datum_name, datumList[index]->name());
2299  strcat( datum_name, "\"" );
2300  if( _3parameterDatum->userDefined() )
2301  fprintf( fp_3param, "*");
2302 
2303  fprintf( fp_3param, "%-6s %-33s%-2s %4.0f %4.0f %4.0f %4.0f %5.0f %4.0f %4.0f %4.0f %4.0f %4.0f \n",
2304  _3parameterDatum->code(),
2305  datum_name,
2306  _3parameterDatum->ellipsoidCode(),
2307  _3parameterDatum->deltaX(),
2308  _3parameterDatum->sigmaX(),
2309  _3parameterDatum->deltaY(),
2310  _3parameterDatum->sigmaY(),
2311  _3parameterDatum->deltaZ(),
2312  _3parameterDatum->sigmaZ(),
2313  ( _3parameterDatum->southLatitude() * _180_OVER_PI ),
2314  ( _3parameterDatum->northLatitude() * _180_OVER_PI ),
2315  ( _3parameterDatum->westLongitude() * _180_OVER_PI ),
2316  ( _3parameterDatum->eastLongitude() * _180_OVER_PI ) );
2317  }
2318 
2319  index++;
2320  }
2321 
2322  fclose( fp_3param );
2323 
2324 }
2325 
2326 
2327 void DatumLibraryImplementation::write7ParamFile()
2328 {
2329 /*
2330  * The function write3ParamFile writes the 7 parameter datums in the datum list
2331  * to the 7_param.dat file.
2332  */
2333 
2334  char datum_name[DATUM_NAME_LENGTH+2];
2335  char *PathName = NULL;
2336  char FileName[FILENAME_LENGTH];
2337  FILE *fp_7param = NULL;
2338 
2339  CCSThreadLock lock(&mutex);
2340 
2341  /*output updated 7-parameter datum table*/
2342  PathName = getenv( "MSPCCS_DATA" );
2343  if (PathName != NULL)
2344  {
2345  strcpy( FileName, PathName );
2346  strcat( FileName, "/" );
2347  }
2348  else
2349  {
2350  strcpy( FileName, "../../data/" );
2351  }
2352  strcat( FileName, "7_param.dat" );
2353 
2354  if ((fp_7param = fopen(FileName, "w")) == NULL)
2355  { /* fatal error */
2356  char message[256] = "";
2357  if (NULL == PathName)
2358  {
2359  strcpy( message, "Environment variable undefined: MSPCCS_DATA." );
2360  }
2361  else
2362  {
2363  strcpy( message, ErrorMessages::datumFileOpenError );
2364  strcat( message, ": 7_param.dat\n" );
2365  }
2366  throw CoordinateConversionException( message );
2367  }
2368 
2369  /* write file */
2370  long index = MAX_WGS;
2371  int endIndex = datum7ParamCount + MAX_WGS;
2372  while( index < endIndex )
2373  {
2374  SevenParameterDatum* _7parameterDatum = ( SevenParameterDatum* )datumList[index];
2375  if( _7parameterDatum )
2376  {
2377  strcpy( datum_name, "\"" );
2378  strcat( datum_name, datumList[index]->name());
2379  strcat( datum_name, "\"" );
2380  if( _7parameterDatum->userDefined() )
2381  fprintf( fp_7param, "*");
2382 
2383  fprintf( fp_7param, "%-6s %-33s%-2s %4.0f %4.0f %4.0f % 4.3f % 4.3f % 4.3f % 4.10f \n",
2384  _7parameterDatum->code(),
2385  datum_name,
2386  _7parameterDatum->ellipsoidCode(),
2387  _7parameterDatum->deltaX(),
2388  _7parameterDatum->deltaY(),
2389  _7parameterDatum->deltaZ(),
2390  _7parameterDatum->rotationX() * SECONDS_PER_RADIAN,
2391  _7parameterDatum->rotationY() * SECONDS_PER_RADIAN,
2392  _7parameterDatum->rotationZ() * SECONDS_PER_RADIAN,
2393  _7parameterDatum->scaleFactor() );
2394  }
2395 
2396  index++;
2397  }
2398 
2399  fclose( fp_7param );
2400 }
2401 
2402 
2403 GeodeticCoordinates* DatumLibraryImplementation::geodeticShiftWGS84ToWGS72( const double WGS84Longitude, const double WGS84Latitude, const double WGS84Height )
2404 {
2405 /*
2406  * The function geodeticShiftWGS84ToWGS72 shifts a geodetic coordinate (latitude, longitude in radians
2407  * and height in meters) relative to WGS84 to a geodetic coordinate
2408  * (latitude, longitude in radians and height in meters) relative to WGS72.
2409  *
2410  * WGS84Longitude : Longitude in radians relative to WGS84 (input)
2411  * WGS84Latitude : Latitude in radians relative to WGS84 (input)
2412  * WGS84Height : Height in meters relative to WGS84 (input)
2413  * WGS72Longitude : Longitude in radians relative to WGS72 (output)
2414  * WGS72Latitude : Latitude in radians relative to WGS72 (output)
2415  * WGS72Height : Height in meters relative to WGS72 (output)
2416  */
2417 
2418  double Delta_Lat;
2419  double Delta_Lon;
2420  double Delta_Hgt;
2421  double WGS84_a; /* Semi-major axis of WGS84 ellipsoid */
2422  double WGS84_f; /* Flattening of WGS84 ellipsoid */
2423  double WGS72_a; /* Semi-major axis of WGS72 ellipsoid */
2424  double WGS72_f; /* Flattening of WGS72 ellipsoid */
2425  double da; /* WGS72_a - WGS84_a */
2426  double df; /* WGS72_f - WGS84_f */
2427  double Q;
2428  double sin_Lat;
2429  double sin2_Lat;
2430 
2431  long wgs84EllipsoidIndex;
2432  _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
2433  _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &WGS84_a, &WGS84_f );
2434 
2435  long wgs72EllipsoidIndex;
2436  _ellipsoidLibraryImplementation->ellipsoidIndex( "WD", &wgs72EllipsoidIndex );
2437  _ellipsoidLibraryImplementation->ellipsoidParameters( wgs72EllipsoidIndex, &WGS72_a, &WGS72_f );
2438 
2439  da = WGS72_a - WGS84_a;
2440  df = WGS72_f - WGS84_f;
2441  Q = PI / 648000;
2442  sin_Lat = sin(WGS84Latitude);
2443  sin2_Lat = sin_Lat * sin_Lat;
2444 
2445  Delta_Lat = (-4.5 * cos(WGS84Latitude)) / (WGS84_a*Q)
2446  + (df * sin(2*WGS84Latitude)) / Q;
2447  Delta_Lat /= SECONDS_PER_RADIAN;
2448  Delta_Lon = -0.554 / SECONDS_PER_RADIAN;
2449  Delta_Hgt = -4.5 * sin_Lat + WGS84_a * df * sin2_Lat - da - 1.4;
2450 
2451  double wgs72Latitude = WGS84Latitude + Delta_Lat;
2452  double wgs72Longitude = WGS84Longitude + Delta_Lon;
2453  double wgs72Height = WGS84Height + Delta_Hgt;
2454 
2455  if (wgs72Latitude > PI_OVER_2)
2456  wgs72Latitude = PI_OVER_2 - (wgs72Latitude - PI_OVER_2);
2457  else if (wgs72Latitude < -PI_OVER_2)
2458  wgs72Latitude = -PI_OVER_2 - (wgs72Latitude + PI_OVER_2);
2459 
2460  if (wgs72Longitude > PI)
2461  wgs72Longitude -= TWO_PI;
2462  if (wgs72Longitude < -PI)
2463  wgs72Longitude += TWO_PI;
2464 
2465  return new GeodeticCoordinates(CoordinateType::geodetic, wgs72Longitude, wgs72Latitude, wgs72Height);
2466 }
2467 
2468 
2469 GeodeticCoordinates* DatumLibraryImplementation::geodeticShiftWGS72ToWGS84( const double WGS72Longitude, const double WGS72Latitude, const double WGS72Height )
2470 {
2471 /*
2472  * The function geodeticShiftWGS72ToWGS84 shifts a geodetic coordinate (latitude, longitude in radians
2473  * and height in meters) relative to WGS72 to a geodetic coordinate
2474  * (latitude, longitude in radians and height in meters) relative to WGS84.
2475  *
2476  * WGS72Longitude : Longitude in radians relative to WGS72 (input)
2477  * WGS72Latitude : Latitude in radians relative to WGS72 (input)
2478  * WGS72Height : Height in meters relative to WGS72 (input)
2479  * WGS84Longitude : Longitude in radians relative to WGS84 (output)
2480  * WGS84Latitude : Latitude in radians relative to WGS84 (output)
2481  * WGS84Height : Height in meters relative to WGS84 (output)
2482  */
2483 
2484  double Delta_Lat;
2485  double Delta_Lon;
2486  double Delta_Hgt;
2487  double WGS84_a; /* Semi-major axis of WGS84 ellipsoid */
2488  double WGS84_f; /* Flattening of WGS84 ellipsoid */
2489  double WGS72_a; /* Semi-major axis of WGS72 ellipsoid */
2490  double WGS72_f; /* Flattening of WGS72 ellipsoid */
2491  double da; /* WGS84_a - WGS72_a */
2492  double df; /* WGS84_f - WGS72_f */
2493  double Q;
2494  double sin_Lat;
2495  double sin2_Lat;
2496 
2497  long wgs84EllipsoidIndex;
2498  _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
2499  _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &WGS84_a, &WGS84_f );
2500 
2501  long wgs72EllipsoidIndex;
2502  _ellipsoidLibraryImplementation->ellipsoidIndex( "WD", &wgs72EllipsoidIndex );
2503  _ellipsoidLibraryImplementation->ellipsoidParameters( wgs72EllipsoidIndex, &WGS72_a, &WGS72_f );
2504 
2505  da = WGS84_a - WGS72_a;
2506  df = WGS84_f - WGS72_f;
2507  Q = PI / 648000;
2508  sin_Lat = sin(WGS72Latitude);
2509  sin2_Lat = sin_Lat * sin_Lat;
2510 
2511  Delta_Lat = (4.5 * cos(WGS72Latitude)) / (WGS72_a*Q) + (df * sin(2*WGS72Latitude)) / Q;
2512  Delta_Lat /= SECONDS_PER_RADIAN;
2513  Delta_Lon = 0.554 / SECONDS_PER_RADIAN;
2514  Delta_Hgt = 4.5 * sin_Lat + WGS72_a * df * sin2_Lat - da + 1.4;
2515 
2516  double wgs84Latitude = WGS72Latitude + Delta_Lat;
2517  double wgs84Longitude = WGS72Longitude + Delta_Lon;
2518  double wgs84Height = WGS72Height + Delta_Hgt;
2519 
2520  if (wgs84Latitude > PI_OVER_2)
2521  wgs84Latitude = PI_OVER_2 - (wgs84Latitude - PI_OVER_2);
2522  else if (wgs84Latitude < -PI_OVER_2)
2523  wgs84Latitude = -PI_OVER_2 - (wgs84Latitude + PI_OVER_2);
2524 
2525  if (wgs84Longitude > PI)
2526  wgs84Longitude -= TWO_PI;
2527  if (wgs84Longitude < -PI)
2528  wgs84Longitude += TWO_PI;
2529 
2530  return new GeodeticCoordinates(CoordinateType::geodetic, wgs84Longitude, wgs84Latitude, wgs84Height);
2531 }
2532 
2533 
2534 CartesianCoordinates* DatumLibraryImplementation::geocentricShiftWGS84ToWGS72( const double X_WGS84, const double Y_WGS84, const double Z_WGS84 )
2535 {
2536 /*
2537  * The function geocentricShiftWGS84ToWGS72 shifts a geocentric coordinate (X, Y, Z in meters) relative
2538  * to WGS84 to a geocentric coordinate (X, Y, Z in meters) relative to WGS72.
2539  *
2540  * X_WGS84 : X coordinate relative to WGS84 (input)
2541  * Y_WGS84 : Y coordinate relative to WGS84 (input)
2542  * Z_WGS84 : Z coordinate relative to WGS84 (input)
2543  * X : X coordinate relative to WGS72 (output)
2544  * Y : Y coordinate relative to WGS72 (output)
2545  * Z : Z coordinate relative to WGS72 (output)
2546  */
2547 
2548  double a_72; /* Semi-major axis in meters of WGS72 ellipsoid */
2549  double f_72; /* Flattening of WGS72 ellipsoid */
2550  double a_84; /* Semi-major axis in meters of WGS84 ellipsoid */
2551  double f_84; /* Flattening of WGS84 ellipsoid */
2552 
2553  /* Set WGS84 ellipsoid params */
2554  long wgs84EllipsoidIndex;
2555  _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
2556  _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &a_84, &f_84 );
2557 
2558  Geocentric geocentric84( a_84, f_84 );
2559 
2560  GeodeticCoordinates* wgs84GeodeticCoordinates = geocentric84.convertToGeodetic( new CartesianCoordinates( CoordinateType::geocentric, X_WGS84, Y_WGS84, Z_WGS84 ) );
2561 
2562  GeodeticCoordinates* wgs72GeodeticCoordinates = geodeticShiftWGS84ToWGS72( wgs84GeodeticCoordinates->longitude(), wgs84GeodeticCoordinates->latitude(), wgs84GeodeticCoordinates->height() );
2563 
2564  /* Set WGS72 ellipsoid params */
2565  long wgs72EllipsoidIndex;
2566  _ellipsoidLibraryImplementation->ellipsoidIndex( "WD", &wgs72EllipsoidIndex );
2567  _ellipsoidLibraryImplementation->ellipsoidParameters( wgs72EllipsoidIndex, &a_72, &f_72 );
2568 
2569  Geocentric geocentric72( a_72, f_72 );
2570 
2571  CartesianCoordinates* wgs72GeocentricCoordinates = geocentric72.convertFromGeodetic( wgs72GeodeticCoordinates );
2572 
2573  delete wgs72GeodeticCoordinates;
2574  delete wgs84GeodeticCoordinates;
2575 
2576  return wgs72GeocentricCoordinates;
2577 }
2578 
2579 
2580 CartesianCoordinates* DatumLibraryImplementation::geocentricShiftWGS72ToWGS84( const double X, const double Y, const double Z )
2581 {
2582 /*
2583  * The function geocentricShiftWGS72ToWGS84 shifts a geocentric coordinate (X, Y, Z in meters) relative
2584  * to WGS72 to a geocentric coordinate (X, Y, Z in meters) relative to WGS84.
2585  *
2586  * X : X coordinate relative to WGS72 (input)
2587  * Y : Y coordinate relative to WGS72 (input)
2588  * Z : Z coordinate relative to WGS72 (input)
2589  * X_WGS84 : X coordinate relative to WGS84 (output)
2590  * Y_WGS84 : Y coordinate relative to WGS84 (output)
2591  * Z_WGS84 : Z coordinate relative to WGS84 (output)
2592  */
2593 
2594  double a_72; /* Semi-major axis in meters of WGS72 ellipsoid */
2595  double f_72; /* Flattening of WGS72 ellipsoid */
2596  double a_84; /* Semi-major axis in meters of WGS84 ellipsoid */
2597  double f_84; /* Flattening of WGS84 ellipsoid */
2598 
2599  /* Set WGS72 ellipsoid params */
2600  long wgs72EllipsoidIndex;
2601  _ellipsoidLibraryImplementation->ellipsoidIndex( "WD", &wgs72EllipsoidIndex );
2602  _ellipsoidLibraryImplementation->ellipsoidParameters( wgs72EllipsoidIndex, &a_72, &f_72 );
2603 
2604  Geocentric geocentric72( a_72, f_72 );
2605  GeodeticCoordinates* wgs72GeodeticCoordinates = geocentric72.convertToGeodetic( new CartesianCoordinates( CoordinateType::geocentric, X, Y, Z ) );
2606 
2607  GeodeticCoordinates* wgs84GeodeticCoordinates = geodeticShiftWGS72ToWGS84( wgs72GeodeticCoordinates->longitude(), wgs72GeodeticCoordinates->latitude(), wgs72GeodeticCoordinates->height() );
2608 
2609  /* Set WGS84 ellipsoid params */
2610  long wgs84EllipsoidIndex;
2611  _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
2612  _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &a_84, &f_84 );
2613 
2614  Geocentric geocentric84( a_84, f_84 );
2615  CartesianCoordinates* wgs84GeocentricCoordinates = geocentric84.convertFromGeodetic( wgs84GeodeticCoordinates );
2616 
2617  delete wgs84GeodeticCoordinates;
2618  delete wgs72GeodeticCoordinates;
2619 
2620  return wgs84GeocentricCoordinates;
2621 }
2622 
2623 
2624 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:01 for GeographicTranslator by doxygen 1.8.2