UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
LambertConformalConic.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: LAMBERT
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates
9  * (latitude and longitude in radians) and Lambert Conformal Conic
10  * (1 or 2 Standard Parallel) projection coordinates (easting and northing in meters) defined
11  * by one standard parallel and specified scale true along that parallel, or two standard parallels.
12  * When both standard parallel parameters
13  * are set to the same latitude value, the result is a Lambert
14  * Conformal Conic projection with one standard parallel at the
15  * specified latitude.
16  *
17  * ERROR HANDLING
18  *
19  * This component checks parameters for valid values. If an invalid value
20  * is found the error code is combined with the current error code using
21  * the bitwise or. This combining allows multiple error codes to be
22  * returned. The possible error codes are:
23  *
24  * LAMBERT_NO_ERROR : No errors occurred in function
25  * LAMBERT_LAT_ERROR : Latitude outside of valid range
26  * (-90 to 90 degrees)
27  * LAMBERT_LON_ERROR : Longitude outside of valid range
28  * (-180 to 360 degrees)
29  * LAMBERT_EASTING_ERROR : Easting outside of valid range
30  * (depends on ellipsoid and projection
31  * parameters)
32  * LAMBERT_NORTHING_ERROR : Northing outside of valid range
33  * (depends on ellipsoid and projection
34  * parameters)
35  * LAMBERT_ORIGIN_LAT_ERROR : Origin latitude outside of valid
36  * range (-89 59 59.0 to 89 59 59.0 degrees)
37  * or within one second of 0 degrees
38  * LAMBERT_CENT_MER_ERROR : Central meridian outside of valid range
39  * (-180 to 360 degrees)
40  * LAMBERT_SCALE_FACTOR_ERROR : Scale factor greater than minimum scale
41  * factor (1.0e-9)
42  * LAMBERT_A_ERROR : Semi-major axis less than or equal to zero
43  * LAMBERT_INV_F_ERROR : Inverse flattening outside of valid range
44  * (250 to 350)
45  *
46  *
47  * REUSE NOTES
48  *
49  * LAMBERT is intended for reuse by any application that performs a Lambert
50  * Conformal Conic (1 or 2 Standard Parallel) projection or its inverse.
51  *
52  * REFERENCES
53  *
54  * Further information on LAMBERT can be found in the Reuse Manual.
55  *
56  * LAMBERT originated from:
57  * Information Technologoy - Spatial Reference Model(SRM)
58  * ISO/IEC FDIS 18026
59  *
60  * LICENSES
61  *
62  * None apply to this component.
63  *
64  * RESTRICTIONS
65  *
66  * LAMBERT has no restrictions.
67  *
68  * ENVIRONMENT
69  *
70  * LAMBERT was tested and certified in the following environments:
71  *
72  * 1. Solaris 2.5 with GCC, version 2.8.1
73  * 2. Windows 98/2000/XP with MS Visual C++, version 6
74  *
75  * MODIFICATIONS
76  *
77  * Date Description
78  * ---- -----------
79  * 03-05-2005 Original Code
80  * 03-02-2007 Original C++ Code
81  * 02-25-2009 Merged Lambert 1 and 2
82  *
83  */
84 
85 
86 /***************************************************************************/
87 /*
88  * INCLUDES
89  */
90 
91 #include <math.h>
92 #include "LambertConformalConic.h"
93 #include "MapProjection5Parameters.h"
94 #include "MapProjection6Parameters.h"
95 #include "MapProjectionCoordinates.h"
96 #include "GeodeticCoordinates.h"
97 #include "CoordinateConversionException.h"
98 #include "ErrorMessages.h"
99 
100 /*
101  * math.h - Standard C math library
102  * LambertConformalConic.h - Is for prototype error checking
103  * CoordinateConversionException.h - Exception handler
104  * MapProjectionCoordinates.h - defines map projection coordinates
105  * GeodeticCoordinates.h - defines geodetic coordinates
106  * ErrorMessages.h - Contains exception messages
107  */
108 
109 
110 using namespace MSP::CCS;
111 
112 
113 /***************************************************************************/
114 /* DEFINES
115  *
116  */
117 
118 const double PI = 3.14159265358979323e0; /* PI */
119 const double PI_OVER_2 = (PI / 2.0);
120 const double PI_OVER_4 = (PI / 4.0);
121 const double PI_OVER_180 = (PI / 180.0);
122 const double MAX_LAT = (( PI * 89.99972222222222) / 180.0); /* 89 59 59.0 degrees in radians */
123 const double TWO_PI = (2.0 * PI);
124 const double MIN_SCALE_FACTOR = 1.0e-9;
125 const double ONE_SECOND = 4.89e-6;
126 
127 
128 /************************************************************************/
129 /* FUNCTIONS
130  *
131  */
132 
133 LambertConformalConic::LambertConformalConic( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double originLatitude, double falseEasting, double falseNorthing, double scaleFactor ) :
134  CoordinateSystem(),
135  coordinateType( CoordinateType::lambertConformalConic1Parallel ),
136  es( 0.08181919084262188000 ),
137  es_OVER_2( .040909595421311 ),
138  Lambert_1_n( 0.70710678118655 ),
139  Lambert_1_rho0( 6388838.2901212 ),
140  Lambert_1_rho_olat( 6388838.2901211 ),
141  Lambert_1_t0( 0.41618115138974 ),
142  Lambert_Origin_Long( 0.0 ),
143  Lambert_Origin_Latitude( (45.0 * PI / 180.0) ),
144  Lambert_False_Easting( 0.0 ),
145  Lambert_False_Northing( 0.0 ),
146  Lambert_Scale_Factor( 1.0 ),
147  Lambert_2_Origin_Lat( (45 * PI / 180) ),
148  Lambert_2_Std_Parallel_1( (40 * PI / 180) ),
149  Lambert_2_Std_Parallel_2( (50 * PI / 180) ),
150  Lambert_Delta_Easting( 400000000.0 ),
151  Lambert_Delta_Northing( 400000000.0 )
152 {
153 /*
154  * The constructor receives the ellipsoid parameters and
155  * Lambert Conformal Conic (1 Standard Parallel) projection parameters as
156  * inputs, and sets the corresponding state variables.
157  * If any errors occur, an exception is thrown with a description of the error.
158  *
159  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
160  * ellipsoidFlattening : Flattening of ellipsoid (input)
161  * centralMeridian : Longitude of origin, in radians (input)
162  * originLatitude : Latitude of origin, in radians (input)
163  * falseEasting : False easting, in meters (input)
164  * falseNorthing : False northing, in meters (input)
165  * scaleFactor : Projection scale factor (input)
166  *
167  */
168 
169  double es2;
170  double inv_f = 1 / ellipsoidFlattening;
171 
172  if (ellipsoidSemiMajorAxis <= 0.0)
173  { /* Semi-major axis must be greater than zero */
174  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
175  }
176  if ((inv_f < 250) || (inv_f > 350))
177  { /* Inverse flattening must be between 250 and 350 */
178  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
179  }
180  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
181  { /* Origin Longitude out of range */
182  throw CoordinateConversionException( ErrorMessages::centralMeridian );
183  }
184 
185  semiMajorAxis = ellipsoidSemiMajorAxis;
186  flattening = ellipsoidFlattening;
187 
188  if (centralMeridian > PI)
189  centralMeridian -= TWO_PI;
190  Lambert_Origin_Long = centralMeridian;
191  Lambert_False_Easting = falseEasting;
192 
193  es2 = 2.0 * flattening - flattening * flattening;
194  es = sqrt(es2);
195  es_OVER_2 = es / 2.0;
196 
197  CommonParameters* parameters = setCommonLambert1StandardParallelParameters(
198  originLatitude, falseNorthing, scaleFactor);
199 
200  Lambert_1_n = parameters->_lambertN;
201  Lambert_1_rho0 = parameters->_lambertRho0;
202  Lambert_1_rho_olat = parameters->_lambertRhoOlat;
203  Lambert_1_t0 = parameters->_lambertT0;
204  Lambert_Origin_Latitude = parameters->_lambertOriginLatitude;
205  Lambert_False_Northing = parameters->_lambertFalseNorthing;
206  Lambert_Scale_Factor = parameters->_lambertScaleFactor;
207 
208  Lambert_2_Origin_Lat = parameters->_lambertOriginLatitude;
209 
210  double sinOlat = sin(Lambert_Origin_Latitude);
211  double esSinOlat = es * sinOlat;//Math.sin(lambertOriginLat);
212  double w0 = sqrt(1 - es2 * sinOlat * sinOlat);
213  double f0 = cos(Lambert_Origin_Latitude) / (w0 * pow(Lambert_1_t0, Lambert_1_n));
214  double c = Lambert_Scale_Factor * f0;
215 
216  double phi = 0.0;
217  double tempPhi = 1.0;
218  Lambert_2_Std_Parallel_1 = calculateLambert2StandardParallel(es2, phi, tempPhi, c);
219 
220  phi = 89.0 * PI_OVER_180;
221  tempPhi = 0.0;
222  Lambert_2_Std_Parallel_2 = calculateLambert2StandardParallel(es2, phi, tempPhi, c);
223 
224  delete parameters;
225  parameters = 0;
226 }
227 
228 
229 LambertConformalConic::LambertConformalConic( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double originLatitude, double standardParallel1, double standardParallel2, double falseEasting, double falseNorthing ) :
230  CoordinateSystem(),
231  coordinateType( CoordinateType::lambertConformalConic2Parallels ),
232  es( 0.081819190842621 ),
233  es_OVER_2( 0.040909595421311 ),
234  Lambert_1_n( 0.70710678118655 ),
235  Lambert_1_rho0( 6388838.2901212 ),
236  Lambert_1_rho_olat( 6388838.2901211 ),
237  Lambert_1_t0( 0.41618115138974 ),
238  Lambert_Origin_Long( 0.0 ),
239  Lambert_Origin_Latitude( (45 * PI / 180) ),
240  Lambert_False_Easting( 0.0 ),
241  Lambert_False_Northing( 0.0 ),
242  Lambert_Scale_Factor( 1.0 ),
243  Lambert_2_Origin_Lat( (45 * PI / 180) ),
244  Lambert_2_Std_Parallel_1( (40 * PI / 180) ),
245  Lambert_2_Std_Parallel_2( (50 * PI / 180) ),
246  Lambert_Delta_Easting( 400000000.0 ),
247  Lambert_Delta_Northing( 400000000.0 )
248 {
249 /*
250  * The constructor receives the ellipsoid parameters and
251  * Lambert Conformal Conic (2 Standard Parallel) projection parameters as inputs, and sets the
252  * corresponding state variables. If any errors occur, an exception is thrown
253  * with a description of the error.
254  *
255  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
256  * ellipsoidFlattening : Flattening of ellipsoid (input)
257  * centralMeridian : Longitude of origin, in radians (input)
258  * originLatitude : Latitude of origin, in radians (input)
259  * standardParallel1 : First standard parallel, in radians (input)
260  * standardParallel2 : Second standard parallel, in radians (input)
261  * falseEasting : False easting, in meters (input)
262  * falseNorthing : False northing, in meters (input)
263  *
264  * Note that when the two standard parallel parameters are both set to the
265  * same latitude value, the result is a Lambert Conformal Conic projection
266  * with one standard parallel at the specified latitude.
267  */
268 
269  double es2;
270  double es_sin;
271  double t0;
272  double t1;
273  double t2;
274  double t_olat;
275  double m0;
276  double m1;
277  double m2;
278  double m_olat;
279  double n; /* Ratio of angle between meridians */
280  double const_value;
281  double Lambert_lat0;
282  double Lambert_k0;
283  double Lambert_false_northing;
284  double inv_f = 1 / ellipsoidFlattening;
285 
286  if (ellipsoidSemiMajorAxis <= 0.0)
287  { /* Semi-major axis must be greater than zero */
288  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
289  }
290  if ((inv_f < 250) || (inv_f > 350))
291  { /* Inverse flattening must be between 250 and 350 */
292  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
293  }
294  if ((originLatitude < -MAX_LAT) || (originLatitude > MAX_LAT))
295  { /* Origin Latitude out of range */
296  throw CoordinateConversionException( ErrorMessages::originLatitude );
297  }
298  if ((standardParallel1 < -MAX_LAT) || (standardParallel1 > MAX_LAT))
299  { /* First Standard Parallel out of range */
300  throw CoordinateConversionException( ErrorMessages::standardParallel1 );
301  }
302  if ((standardParallel2 < -MAX_LAT) || (standardParallel2 > MAX_LAT))
303  { /* Second Standard Parallel out of range */
304  throw CoordinateConversionException( ErrorMessages::standardParallel2 );
305  }
306  if ((standardParallel1 == 0) && (standardParallel2 == 0))
307  { /* First & Second Standard Parallels are both 0 */
308  throw CoordinateConversionException( ErrorMessages::standardParallel1_2 );
309  }
310  if (standardParallel1 == -standardParallel2)
311  { /* Parallels are the negation of each other */
312  throw CoordinateConversionException(
313  ErrorMessages::standardParallelHemisphere );
314  }
315  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
316  { /* Central meridian out of range */
317  throw CoordinateConversionException( ErrorMessages::centralMeridian );
318  }
319 
320  semiMajorAxis = ellipsoidSemiMajorAxis;
321  flattening = ellipsoidFlattening;
322 
323  Lambert_2_Origin_Lat = originLatitude;
324  Lambert_2_Std_Parallel_1 = standardParallel1;
325  Lambert_2_Std_Parallel_2 = standardParallel2;
326 
327  if (centralMeridian > PI)
328  centralMeridian -= TWO_PI;
329  Lambert_Origin_Long = centralMeridian;
330  Lambert_False_Easting = falseEasting;
331 
332  es2 = 2 * flattening - flattening * flattening;
333  es = sqrt(es2);
334  es_OVER_2 = es / 2.0;
335 
336  if (fabs(Lambert_2_Std_Parallel_1 - Lambert_2_Std_Parallel_2) > 1.0e-10)
337  {
338  es_sin = esSin(sin(originLatitude));
339  m_olat = lambertM(cos(originLatitude), es_sin);
340  t_olat = lambertT(originLatitude, es_sin);
341 
342  es_sin = esSin(sin(Lambert_2_Std_Parallel_1));
343  m1 = lambertM(cos(Lambert_2_Std_Parallel_1), es_sin);
344  t1 = lambertT(Lambert_2_Std_Parallel_1, es_sin);
345 
346  es_sin = esSin(sin(Lambert_2_Std_Parallel_2));
347  m2 = lambertM(cos(Lambert_2_Std_Parallel_2), es_sin);
348  t2 = lambertT(Lambert_2_Std_Parallel_2, es_sin);
349 
350  n = log(m1 / m2) / log(t1 / t2);
351 
352  Lambert_lat0 = asin(n);
353 
354  es_sin = esSin(sin(Lambert_lat0));
355  m0 = lambertM(cos(Lambert_lat0), es_sin);
356  t0 = lambertT(Lambert_lat0, es_sin);
357 
358  Lambert_k0 = (m1 / m0) * (pow(t0 / t1, n));
359 
360  const_value = ((semiMajorAxis * m2) / (n * pow(t2, n)));
361 
362  Lambert_false_northing =
363  (const_value * pow(t_olat, n)) - (const_value * pow(t0, n))
364  + falseNorthing;
365  }
366  else
367  {
368  Lambert_lat0 = Lambert_2_Std_Parallel_1;
369  Lambert_k0 = 1.0;
370  Lambert_false_northing = falseNorthing;
371  }
372 
373  CommonParameters* parameters = setCommonLambert1StandardParallelParameters(
374  Lambert_lat0, Lambert_false_northing, Lambert_k0);
375 
376  Lambert_1_n = parameters->_lambertN;
377  Lambert_1_rho0 = parameters->_lambertRho0;
378  Lambert_1_rho_olat = parameters->_lambertRhoOlat;
379  Lambert_1_t0 = parameters->_lambertT0;
380  Lambert_Origin_Latitude = parameters->_lambertOriginLatitude;
381  Lambert_False_Northing = parameters->_lambertFalseNorthing;
382  Lambert_Scale_Factor = parameters->_lambertScaleFactor;
383 
384  delete parameters;
385  parameters = 0;
386 }
387 
388 LambertConformalConic::LambertConformalConic( const LambertConformalConic &lcc )
389 {
390  coordinateType = lcc.coordinateType;
391  semiMajorAxis = lcc.semiMajorAxis;
392  flattening = lcc.flattening;
393  es = lcc.es;
394  es_OVER_2 = lcc.es_OVER_2;
395  Lambert_1_n = lcc.Lambert_1_n;
396  Lambert_1_rho0 = lcc.Lambert_1_rho0;
397  Lambert_1_rho_olat = lcc.Lambert_1_rho_olat;
398  Lambert_1_t0 = lcc.Lambert_1_t0;
399  Lambert_Origin_Long = lcc.Lambert_Origin_Long;
400  Lambert_Origin_Latitude = lcc.Lambert_Origin_Latitude;
401  Lambert_False_Easting = lcc.Lambert_False_Easting;
402  Lambert_False_Northing = lcc.Lambert_False_Northing;
403  Lambert_Scale_Factor = lcc.Lambert_Scale_Factor;
404  Lambert_2_Origin_Lat = lcc.Lambert_2_Origin_Lat;
405  Lambert_2_Std_Parallel_1 = lcc.Lambert_2_Std_Parallel_1;
406  Lambert_2_Std_Parallel_2 = lcc.Lambert_2_Std_Parallel_2;
407  Lambert_Delta_Easting = lcc.Lambert_Delta_Easting;
408  Lambert_Delta_Northing = lcc.Lambert_Delta_Northing;
409 }
410 
411 
412 LambertConformalConic::~LambertConformalConic()
413 {
414 }
415 
416 
417 LambertConformalConic& LambertConformalConic::operator=(
418  const LambertConformalConic &lcc )
419 {
420  if( this != &lcc )
421  {
422  coordinateType = lcc.coordinateType;
423  semiMajorAxis = lcc.semiMajorAxis;
424  flattening = lcc.flattening;
425  es = lcc.es;
426  es_OVER_2 = lcc.es_OVER_2;
427  Lambert_1_n = lcc.Lambert_1_n;
428  Lambert_1_rho0 = lcc.Lambert_1_rho0;
429  Lambert_1_rho_olat = lcc.Lambert_1_rho_olat;
430  Lambert_1_t0 = lcc.Lambert_1_t0;
431  Lambert_Origin_Long = lcc.Lambert_Origin_Long;
432  Lambert_Origin_Latitude = lcc.Lambert_Origin_Latitude;
433  Lambert_False_Easting = lcc.Lambert_False_Easting;
434  Lambert_False_Northing = lcc.Lambert_False_Northing;
435  Lambert_Scale_Factor = lcc.Lambert_Scale_Factor;
436  Lambert_2_Origin_Lat = lcc.Lambert_2_Origin_Lat;
437  Lambert_2_Std_Parallel_1 = lcc.Lambert_2_Std_Parallel_1;
438  Lambert_2_Std_Parallel_2 = lcc.Lambert_2_Std_Parallel_2;
439  Lambert_Delta_Easting = lcc.Lambert_Delta_Easting;
440  Lambert_Delta_Northing = lcc.Lambert_Delta_Northing;
441  }
442 
443  return *this;
444 }
445 
446 
447 MapProjection5Parameters* LambertConformalConic::get1StandardParallelParameters() const
448 {
449 /*
450  * The function get1StandardParallelParameters returns the current ellipsoid
451  * parameters and Lambert Conformal Conic (1 Standard Parallel) projection parameters.
452  *
453  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
454  * ellipsoidFlattening : Flattening of ellipsoid (output)
455  * centralMeridian : Longitude of origin, in radians (output)
456  * originLatitude : Latitude of origin, in radians (output)
457  * falseEasting : False easting, in meters (output)
458  * falseNorthing : False northing, in meters (output)
459  * scaleFactor : Projection scale factor (output)
460  */
461 
462  return new MapProjection5Parameters( CoordinateType::lambertConformalConic1Parallel, Lambert_Origin_Long, Lambert_Origin_Latitude, Lambert_Scale_Factor, Lambert_False_Easting, Lambert_False_Northing );
463 }
464 
465 
466 MapProjection6Parameters* LambertConformalConic::get2StandardParallelParameters() const
467 {
468 /*
469  * The function get2StandardParallelParameters returns the current ellipsoid
470  * parameters and Lambert Conformal Conic (2 Standard Parallel) projection parameters.
471  *
472  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
473  * ellipsoidFlattening : Flattening of ellipsoid (output)
474  * centralMeridian : Longitude of origin, in radians (output)
475  * originLatitude : Latitude of origin, in radians (output)
476  * standardParallel1 : First standard parallel, in radians (output)
477  * standardParallel2 : Second standard parallel, in radians (output)
478  * falseEasting : False easting, in meters (output)
479  * falseNorthing : False northing, in meters (output)
480  */
481 
482  return new MapProjection6Parameters(
483  CoordinateType::lambertConformalConic2Parallels,
484  Lambert_Origin_Long, Lambert_2_Origin_Lat,
485  Lambert_2_Std_Parallel_1, Lambert_2_Std_Parallel_2,
486  Lambert_False_Easting, Lambert_False_Northing );
487 }
488 
489 
490 MSP::CCS::MapProjectionCoordinates* LambertConformalConic::convertFromGeodetic(
491  MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
492 {
493 /*
494  * The function convertFromGeodetic converts Geodetic (latitude and
495  * longitude) coordinates to Lambert Conformal Conic (1 or 2 Standard Parallel)
496  * projection (easting and northing) coordinates, according to the current
497  * ellipsoid and Lambert Conformal Conic (1 or 2 Standard Parallel) projection
498  * parameters. If any errors occur, an exception is thrown with a
499  * description of the error.
500  *
501  * longitude : Longitude, in radians (input)
502  * latitude : Latitude, in radians (input)
503  * easting : Easting (X), in meters (output)
504  * northing : Northing (Y), in meters (output)
505  */
506 
507  double t;
508  double rho;
509  double dlam;
510  double theta;
511 
512  double longitude = geodeticCoordinates->longitude();
513  double latitude = geodeticCoordinates->latitude();
514 
515  if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
516  { /* Latitude out of range */
517  throw CoordinateConversionException( ErrorMessages::latitude );
518  }
519  if ((longitude < -PI) || (longitude > TWO_PI))
520  { /* Longitude out of range */
521  throw CoordinateConversionException( ErrorMessages::longitude );
522  }
523 
524  if (fabs(fabs(latitude) - PI_OVER_2) > 1.0e-10)
525  {
526  t = lambertT(latitude, esSin(sin(latitude)));
527  rho = Lambert_1_rho0 * pow(t / Lambert_1_t0, Lambert_1_n);
528  }
529  else
530  {
531  if ((latitude * Lambert_1_n) <= 0)
532  { /* Point can not be projected */
533  throw CoordinateConversionException( ErrorMessages::latitude );
534  }
535  rho = 0.0;
536  }
537 
538  dlam = longitude - Lambert_Origin_Long;
539 
540  if (dlam > PI)
541  {
542  dlam -= TWO_PI;
543  }
544  if (dlam < -PI)
545  {
546  dlam += TWO_PI;
547  }
548 
549  theta = Lambert_1_n * dlam;
550 
551  double easting = rho * sin(theta) + Lambert_False_Easting;
552  double northing = Lambert_1_rho_olat - rho * cos(theta) + Lambert_False_Northing;
553 
554  return new MapProjectionCoordinates( coordinateType, easting, northing );
555 }
556 
557 
558 MSP::CCS::GeodeticCoordinates* LambertConformalConic::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
559 {
560 /*
561  * The function convertToGeodetic converts Lambert Conformal
562  * Conic (1 or 2 Standard Parallel) projection (easting and northing) coordinates to Geodetic
563  * (latitude and longitude) coordinates, according to the current ellipsoid
564  * and Lambert Conformal Conic (1 or 2 Standard Parallel) projection parameters. If any
565  * errors occur, an exception is thrown with a description of the error.
566  *
567  * easting : Easting (X), in meters (input)
568  * northing : Northing (Y), in meters (input)
569  * longitude : Longitude, in radians (output)
570  * latitude : Latitude, in radians (output)
571  */
572 
573  double dx;
574  double dy;
575  double rho;
576  double rho_olat_MINUS_dy;
577  double t;
578  double PHI;
579  double es_sin;
580  double tempPHI = 0.0;
581  double theta = 0.0;
582  double tolerance = 4.85e-10;
583  int count = 30;
584  double longitude, latitude;
585 
586  double easting = mapProjectionCoordinates->easting();
587  double northing = mapProjectionCoordinates->northing();
588 
589  if ((easting < (Lambert_False_Easting - Lambert_Delta_Easting))
590  ||(easting > (Lambert_False_Easting + Lambert_Delta_Easting)))
591  { /* Easting out of range */
592  throw CoordinateConversionException( ErrorMessages::easting );
593  }
594  if ((northing < (Lambert_False_Northing - Lambert_Delta_Northing))
595  || (northing > (Lambert_False_Northing + Lambert_Delta_Northing)))
596  { /* Northing out of range */
597  throw CoordinateConversionException( ErrorMessages::northing );
598  }
599 
600  dy = northing - Lambert_False_Northing;
601  dx = easting - Lambert_False_Easting;
602  rho_olat_MINUS_dy = Lambert_1_rho_olat - dy;
603  rho = sqrt(dx * dx + (rho_olat_MINUS_dy) * (rho_olat_MINUS_dy));
604 
605  if (Lambert_1_n < 0.0)
606  {
607  rho *= -1.0;
608  dx *= -1.0;
609  rho_olat_MINUS_dy *= -1.0;
610  }
611 
612  if (rho != 0.0)
613  {
614  theta = atan2(dx, rho_olat_MINUS_dy) / Lambert_1_n;
615  t = Lambert_1_t0 * pow(rho / Lambert_1_rho0, 1 / Lambert_1_n);
616  PHI = PI_OVER_2 - 2.0 * atan(t);
617  while (fabs(PHI - tempPHI) > tolerance && count)
618  {
619  tempPHI = PHI;
620  es_sin = esSin(sin(PHI));
621  PHI = PI_OVER_2 - 2.0 * atan(t * pow((1.0 - es_sin) / (1.0 + es_sin), es_OVER_2));
622  count --;
623  }
624 
625  if(!count)
626  throw CoordinateConversionException( ErrorMessages::northing );
627 
628  latitude = PHI;
629  longitude = theta + Lambert_Origin_Long;
630 
631  if (fabs(latitude) < 2.0e-7) /* force lat to 0 to avoid -0 degrees */
632  latitude = 0.0;
633  if (latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
634  latitude = PI_OVER_2;
635  else if (latitude < -PI_OVER_2)
636  latitude = -PI_OVER_2;
637 
638  if (longitude > PI)
639  {
640  if (longitude - PI < 3.5e-6)
641  longitude = PI;
642  else
643  longitude -= TWO_PI;
644  }
645  if (longitude < -PI)
646  {
647  if (fabs(longitude + PI) < 3.5e-6)
648  longitude = -PI;
649  else
650  longitude += TWO_PI;
651  }
652 
653  if (fabs(longitude) < 2.0e-7) /* force lon to 0 to avoid -0 degrees */
654  longitude = 0.0;
655  if (longitude > PI) /* force distorted values to 180, -180 degrees */
656  longitude = PI;
657  else if (longitude < -PI)
658  longitude = -PI;
659  }
660  else
661  {
662  if (Lambert_1_n > 0.0)
663  latitude = PI_OVER_2;
664  else
665  latitude = -PI_OVER_2;
666  longitude = Lambert_Origin_Long;
667  }
668 
669  return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
670 }
671 
672 
673 LambertConformalConic::CommonParameters*
674 LambertConformalConic::setCommonLambert1StandardParallelParameters(
675  double originLatitude, double falseNorthing, double scaleFactor )
676 {
689  double _esSin;
690  double m0;
691 
692  if (((originLatitude < -MAX_LAT) || (originLatitude > MAX_LAT)) ||
693  (originLatitude > -ONE_SECOND) && (originLatitude < ONE_SECOND))
694  { /* Origin Latitude out of range */
695  throw CoordinateConversionException( ErrorMessages::originLatitude );
696  }
697  if (scaleFactor < MIN_SCALE_FACTOR)
698  {
699  throw CoordinateConversionException( ErrorMessages::scaleFactor );
700  }
701 
702  CommonParameters* parameters = new CommonParameters();
703 
704  parameters->_lambertOriginLatitude = originLatitude;
705  parameters->_lambertFalseNorthing = falseNorthing;
706  parameters->_lambertScaleFactor = scaleFactor;
707 
708  parameters->_lambertN = sin(originLatitude);
709 
710  _esSin = esSin(sin(originLatitude));
711 
712  m0 = cos(originLatitude) / sqrt(1.0 - _esSin * _esSin);
713  parameters->_lambertT0 = lambertT(originLatitude, _esSin);
714 
715  parameters->_lambertRho0 =
716  semiMajorAxis * scaleFactor * m0 / parameters->_lambertN;
717 
718  parameters->_lambertRhoOlat = parameters->_lambertRho0;
719 
720  return parameters;
721 }
722 
723 
724 double LambertConformalConic::calculateLambert2StandardParallel(
725  double es2, double phi, double tempPhi, double c)
726 {
727 /*
728  * Calculates the
729  * Lambert Conformal Conic (2 Standard Parallel) standard parallel value.
730  *
731  */
732  double tolerance = 1.0e-11;
733  int count = 30;
734  while (fabs(phi - tempPhi) > tolerance && count > 0)
735  {
736  tempPhi = phi;
737  double sinPhi = sin(phi);
738  double esSinPhi = es * sinPhi;
739  double tPhi = lambertT(phi, esSinPhi);
740  double wPhi = sqrt(1 - es2 * sinPhi * sinPhi);
741  double fPhi = cos(phi) / (wPhi * pow(tPhi, Lambert_1_n));
742  double fprPhi = ((Lambert_1_n - sinPhi) * (1 - es2)) / (pow(wPhi, 3) * pow(tPhi, Lambert_1_n));
743 
744  phi = phi + (c - fPhi) / fprPhi;
745 
746  count--;
747  }
748 
749  return phi;
750 }
751 
752 
753 double LambertConformalConic::lambertM( double clat, double essin )
754 {
755  return clat / sqrt(1.0 - essin * essin);
756 }
757 
758 
759 double LambertConformalConic::lambertT( double lat, double essin )
760 {
761  return tan(PI_OVER_4 - lat / 2) / pow((1.0 - essin) / (1.0 + essin), es_OVER_2);
762 }
763 
764 
765 double LambertConformalConic::esSin( double sinlat )
766 {
767  return es * sinlat;
768 }
769 
770 
771 
772 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:02 for GeographicTranslator by doxygen 1.8.2