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;