UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
ObliqueMercator.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: OBLIQUE MERCATOR
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates
9  * (latitude and longitude in radians) and Oblique Mercator
10  * projection coordinates (easting and northing in meters).
11  *
12  * ERROR HANDLING
13  *
14  * This component checks parameters for valid values. If an invalid value
15  * is found the error code is combined with the current error code using
16  * the bitwise or. This combining allows multiple error codes to be
17  * returned. The possible error codes are:
18  *
19  * OMERC_NO_ERROR : No errors occurred in function
20  * OMERC_LAT_ERROR : Latitude outside of valid range
21  * (-90 to 90 degrees)
22  * OMERC_LON_ERROR : Longitude outside of valid range
23  * (-180 to 360 degrees)
24  * OMERC_ORIGIN_LAT_ERROR : Origin latitude outside of valid range
25  * (-89 to 89 degrees)
26  * OMERC_LAT1_ERROR : First latitude outside of valid range
27  * (-89 to 89 degrees, excluding 0)
28  * OMERC_LAT2_ERROR : First latitude outside of valid range
29  * (-89 to 89 degrees)
30  * OMERC_LON1_ERROR : First longitude outside of valid range
31  * (-180 to 360 degrees)
32  * OMERC_LON2_ERROR : Second longitude outside of valid range
33  * (-180 to 360 degrees)
34  * OMERC_LAT1_LAT2_ERROR : First and second latitudes can not be equal
35  * OMERC_DIFF_HEMISPHERE_ERROR: First and second latitudes can not be
36  * in different hemispheres
37  * OMERC_EASTING_ERROR : Easting outside of valid range
38  * (depends on ellipsoid and projection
39  * parameters)
40  * OMERC_NORTHING_ERROR : Northing outside of valid range
41  * (depends on ellipsoid and projection
42  * parameters)
43  * OMERC_A_ERROR : Semi-major axis less than or equal to zero
44  * OMERC_INV_F_ERROR : Inverse flattening outside of valid range
45  * (250 to 350)
46  * OMERC_SCALE_FACTOR_ERROR : Scale factor outside of valid
47  * range (0.3 to 3.0)
48  * OMERC_LON_WARNING : Distortion will result if longitude is 90 degrees or more
49  * from the Central Meridian
50  *
51  * REUSE NOTES
52  *
53  * OBLIQUE MERCATOR is intended for reuse by any application that
54  * performs an Oblique Mercator projection or its inverse.
55  *
56  * REFERENCES
57  *
58  * Further information on OBLIQUE MERCATOR can be found in the Reuse Manual.
59  *
60  * OBLIQUE MERCATOR originated from: U.S. Army Topographic Engineering Center
61  * Geospatial Information Division
62  * 7701 Telegraph Road
63  * Alexandria, VA 22310-3864
64  *
65  * LICENSES
66  *
67  * None apply to this component.
68  *
69  * RESTRICTIONS
70  *
71  * OBLIQUE MERCATOR has no restrictions.
72  *
73  * ENVIRONMENT
74  *
75  * OBLIQUE MERCATOR was tested and certified in the following environments:
76  *
77  * 1. Solaris 2.5 with GCC, version 2.8.1
78  * 2. MSDOS with MS Visual C++, version 6
79  *
80  * MODIFICATIONS
81  *
82  * Date Description
83  * ---- -----------
84  * 06-07-00 Original Code
85  * 03-02-07 Original C++ Code
86  * 05-11-11 BAEts28017 - Fix Oblique Mercator near poles
87  *
88  */
89 
90 
91 /***************************************************************************/
92 /*
93  * INCLUDES
94  */
95 
96 #include <math.h>
97 #include "ObliqueMercator.h"
98 #include "ObliqueMercatorParameters.h"
99 #include "MapProjectionCoordinates.h"
100 #include "GeodeticCoordinates.h"
101 #include "CoordinateConversionException.h"
102 #include "ErrorMessages.h"
103 #include "WarningMessages.h"
104 
105 /*
106  * math.h - Standard C math library
107  * ObliqueMercator.h - Is for prototype error checking
108  * MapProjectionCoordinates.h - defines map projection coordinates
109  * GeodeticCoordinates.h - defines geodetic coordinates
110  * CoordinateConversionException.h - Exception handler
111  * ErrorMessages.h - Contains exception messages
112  * WarningMessages.h - Contains warning messages
113  */
114 
115 
116 using namespace MSP::CCS;
117 
118 
119 /***************************************************************************/
120 /* DEFINES
121  *
122  */
123 
124 const double PI = 3.14159265358979323e0; /* PI */
125 const double PI_OVER_2 = ( PI / 2.0);
126 const double PI_OVER_4 = ( PI / 4.0);
127 const double TWO_PI = ( 2.0 * PI);
128 const double MIN_SCALE_FACTOR = 0.3;
129 const double MAX_SCALE_FACTOR = 3.0;
130 
131 
132 /************************************************************************/
133 /* FUNCTIONS
134  *
135  */
136 
137 ObliqueMercator::ObliqueMercator( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double originLatitude, double longitude1, double latitude1, double longitude2, double latitude2, double falseEasting, double falseNorthing, double scaleFactor ) :
138  CoordinateSystem(),
139  es( 0.08181919084262188000 ),
140  es_OVER_2( .040909595421311 ),
141  OMerc_A( 6383471.9177251 ),
142  OMerc_B( 1.0008420825413 ),
143  OMerc_E( 1.0028158089754 ),
144  OMerc_gamma( .41705894983580 ),
145  OMerc_azimuth( .60940407333533 ),
146  OMerc_Origin_Long( -.46732023406900 ),
147  cos_gamma( .91428423352628 ),
148  sin_gamma( .40507325303611 ),
149  sin_azimuth( .57237890829911 ),
150  cos_azimuth( .81998925927985 ),
151  A_over_B( 6378101.0302010 ),
152  B_over_A( 1.5678647849335e-7 ),
153  OMerc_u( 5632885.2272051 ),
154  OMerc_Origin_Lat( ((45.0 * PI) / 180.0) ),
155  OMerc_Lon_1( ((-5.0 * PI) / 180.0) ),
156  OMerc_Lat_1( ((40.0 * PI) / 180.0) ),
157  OMerc_Lon_2( ((5.0 * PI) / 180.0) ),
158  OMerc_Lat_2( ((50.0 * PI) / 180.0) ),
159  OMerc_False_Easting( 0.0 ),
160  OMerc_False_Northing( 0.0 ),
161  OMerc_Scale_Factor( 1.0 ),
162  OMerc_Delta_Northing( 40000000.0 ),
163  OMerc_Delta_Easting( 40000000.0 )
164 {
165 /*
166  * The constructor receives the ellipsoid parameters and
167  * projection parameters as inputs, and sets the corresponding state
168  * variables. If any errors occur, an exception is thrown with a description
169  * of the error.
170  *
171  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
172  * ellipsoidFlattening : Flattening of ellipsoid (input)
173  * originLatitude : Latitude, in radians, at which the (input)
174  * point scale factor is 1.0
175  * longitude1 : Longitude, in radians, of first point lying on
176  * central line (input)
177  * latitude1 : Latitude, in radians, of first point lying on
178  * central line (input)
179  * longitude2 : Longitude, in radians, of second point lying on
180  * central line (input)
181  * latitude2 : Latitude, in radians, of second point lying on
182  * central line (input)
183  * falseEasting : A coordinate value, in meters, assigned to the
184  * central meridian of the projection (input)
185  * falseNorthing : A coordinate value, in meters, assigned to the
186  * origin latitude of the projection (input)
187  * scaleFactor : Multiplier which reduces distances in the
188  * projection to the actual distance on the
189  * ellipsoid (input)
190  */
191 
192  double inv_f = 1 / ellipsoidFlattening;
193  double es2, one_MINUS_es2;
194  double cos_olat, cos_olat2;
195  double sin_olat, sin_olat2, es2_sin_olat2;
196  double t0, t1, t2;
197  double D, D2, D2_MINUS_1, sqrt_D2_MINUS_1;
198  double H, L, LH;
199  double E2;
200  double F, G, J, P;
201  double dlon;
202 
203  if (ellipsoidSemiMajorAxis <= 0.0)
204  { /* Semi-major axis must be greater than zero */
205  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
206  }
207  if ((inv_f < 250) || (inv_f > 350))
208  { /* Inverse flattening must be between 250 and 350 */
209  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
210  }
211  if ((originLatitude <= -PI_OVER_2) || (originLatitude >= PI_OVER_2))
212  { /* origin latitude out of range - can not be at a pole */
213  throw CoordinateConversionException( ErrorMessages::originLatitude );
214  }
215  if ((latitude1 <= -PI_OVER_2) || (latitude1 >= PI_OVER_2))
216  { /* first latitude out of range - can not be at a pole */
217  throw CoordinateConversionException( ErrorMessages::latitude1 );
218  }
219  if ((latitude2 <= -PI_OVER_2) || (latitude2 >= PI_OVER_2))
220  { /* second latitude out of range - can not be at a pole */
221  throw CoordinateConversionException( ErrorMessages::latitude2 );
222  }
223  if (latitude1 == 0.0)
224  { /* first latitude can not be at the equator */
225  throw CoordinateConversionException( ErrorMessages::latitude1 );
226  }
227  if (latitude1 == latitude2)
228  { /* first and second latitudes can not be equal */
229  throw CoordinateConversionException( ErrorMessages::latitude2 );
230  }
231  if (((latitude1 < 0.0) && (latitude2 > 0.0)) ||
232  ((latitude1 > 0.0) && (latitude2 < 0.0)))
233  { /*first and second points can not be in different hemispheres */
234  throw CoordinateConversionException( ErrorMessages::omercHemisphere );
235  }
236  if ((longitude1 < -PI) || (longitude1 > TWO_PI))
237  { /* first longitude out of range */
238  throw CoordinateConversionException( ErrorMessages::longitude1 );
239  }
240  if ((longitude2 < -PI) || (longitude2 > TWO_PI))
241  { /* first longitude out of range */
242  throw CoordinateConversionException( ErrorMessages::longitude2 );
243  }
244  if ((scaleFactor < MIN_SCALE_FACTOR) || (scaleFactor > MAX_SCALE_FACTOR))
245  { /* scale factor out of range */
246  throw CoordinateConversionException( ErrorMessages::scaleFactor );
247  }
248 
249  semiMajorAxis = ellipsoidSemiMajorAxis;
250  flattening = ellipsoidFlattening;
251 
252  OMerc_Origin_Lat = originLatitude;
253  OMerc_Lon_1 = longitude1;
254  OMerc_Lat_1 = latitude1;
255  OMerc_Lon_2 = longitude2;
256  OMerc_Lat_2 = latitude2;
257  OMerc_False_Northing = falseNorthing;
258  OMerc_False_Easting = falseEasting;
259  OMerc_Scale_Factor = scaleFactor;
260 
261  es2 = 2 * flattening - flattening * flattening;
262  es = sqrt(es2);
263  one_MINUS_es2 = 1 - es2;
264  es_OVER_2 = es / 2.0;
265 
266  cos_olat = cos(OMerc_Origin_Lat);
267  cos_olat2 = cos_olat * cos_olat;
268  sin_olat = sin(OMerc_Origin_Lat);
269  sin_olat2 = sin_olat * sin_olat;
270  es2_sin_olat2 = es2 * sin_olat2;
271 
272  OMerc_B = sqrt(1 + (es2 * cos_olat2 * cos_olat2) / one_MINUS_es2);
273  OMerc_A = (semiMajorAxis * OMerc_B * OMerc_Scale_Factor * sqrt(one_MINUS_es2)) / (1.0 - es2_sin_olat2);
274  A_over_B = OMerc_A / OMerc_B;
275  B_over_A = OMerc_B / OMerc_A;
276 
277  t0 = omercT(OMerc_Origin_Lat, es * sin_olat, es_OVER_2);
278  t1 = omercT(OMerc_Lat_1, es * sin(OMerc_Lat_1), es_OVER_2);
279  t2 = omercT(OMerc_Lat_2, es * sin(OMerc_Lat_2), es_OVER_2);
280 
281  D = (OMerc_B * sqrt(one_MINUS_es2)) / (cos_olat * sqrt(1.0 - es2_sin_olat2));
282  D2 = D * D;
283  if (D2 < 1.0)
284  D2 = 1.0;
285  D2_MINUS_1 = D2 - 1.0;
286  sqrt_D2_MINUS_1 = sqrt(D2_MINUS_1);
287  if (D2_MINUS_1 > 1.0e-10)
288  {
289  if (OMerc_Origin_Lat >= 0.0)
290  OMerc_E = (D + sqrt_D2_MINUS_1) * pow(t0, OMerc_B);
291  else
292  OMerc_E = (D - sqrt_D2_MINUS_1) * pow(t0, OMerc_B);
293  }
294  else
295  OMerc_E = D * pow(t0, OMerc_B);
296  H = pow(t1, OMerc_B);
297  L = pow(t2, OMerc_B);
298  F = OMerc_E / H;
299  G = (F - 1.0 / F) / 2.0;
300  E2 = OMerc_E * OMerc_E;
301  LH = L * H;
302  J = (E2 - LH) / (E2 + LH);
303  P = (L - H) / (L + H);
304 
305  dlon = OMerc_Lon_1 - OMerc_Lon_2;
306  if (dlon < -PI )
307  OMerc_Lon_2 -= TWO_PI;
308  if (dlon > PI)
309  OMerc_Lon_2 += TWO_PI;
310  dlon = OMerc_Lon_1 - OMerc_Lon_2;
311  OMerc_Origin_Long = (OMerc_Lon_1 + OMerc_Lon_2) / 2.0 - (atan(J * tan(OMerc_B * dlon / 2.0) / P)) / OMerc_B;
312 
313  dlon = OMerc_Lon_1 - OMerc_Origin_Long;
314  if (dlon < -PI )
315  OMerc_Origin_Long -= TWO_PI;
316  if (dlon > PI)
317  OMerc_Origin_Long += TWO_PI;
318 
319  dlon = OMerc_Lon_1 - OMerc_Origin_Long;
320  OMerc_gamma = atan(sin(OMerc_B * dlon) / G);
321  cos_gamma = cos(OMerc_gamma);
322  sin_gamma = sin(OMerc_gamma);
323 
324  OMerc_azimuth = asin(D * sin_gamma);
325  cos_azimuth = cos(OMerc_azimuth);
326  sin_azimuth = sin(OMerc_azimuth);
327 
328  if (OMerc_Origin_Lat >= 0)
329  OMerc_u = A_over_B * atan(sqrt_D2_MINUS_1/cos_azimuth);
330  else
331  OMerc_u = -A_over_B * atan(sqrt_D2_MINUS_1/cos_azimuth);
332 }
333 
334 
335 ObliqueMercator::ObliqueMercator( const ObliqueMercator &om )
336 {
337  semiMajorAxis = om.semiMajorAxis;
338  flattening = om.flattening;
339  es = om.es;
340  es_OVER_2 = om.es_OVER_2;
341  OMerc_A = om.OMerc_A;
342  OMerc_B = om.OMerc_B;
343  OMerc_E = om.OMerc_E;
344  OMerc_gamma = om.OMerc_gamma;
345  OMerc_azimuth = om.OMerc_azimuth;
346  OMerc_Origin_Long = om.OMerc_Origin_Long;
347  cos_gamma = om.cos_gamma;
348  sin_gamma = om.sin_gamma;
349  sin_azimuth = om.sin_azimuth;
350  cos_azimuth = om.cos_azimuth;
351  A_over_B = om.A_over_B;
352  B_over_A = om.B_over_A;
353  OMerc_u = om.OMerc_u;
354  OMerc_Origin_Lat = om.OMerc_Origin_Lat;
355  OMerc_Lon_1 = om.OMerc_Lon_1;
356  OMerc_Lat_1 = om.OMerc_Lat_1;
357  OMerc_Lon_2 = om.OMerc_Lon_2;
358  OMerc_Lat_2 = om.OMerc_Lat_2;
359  OMerc_False_Easting = om.OMerc_False_Easting;
360  OMerc_False_Northing = om.OMerc_False_Northing;
361  OMerc_Scale_Factor = om.OMerc_Scale_Factor;
362  OMerc_Delta_Northing = om.OMerc_Delta_Northing;
363  OMerc_Delta_Easting = om.OMerc_Delta_Easting;
364 }
365 
366 
367 ObliqueMercator::~ObliqueMercator()
368 {
369 }
370 
371 
372 ObliqueMercator& ObliqueMercator::operator=( const ObliqueMercator &om )
373 {
374  if( this != &om )
375  {
376  semiMajorAxis = om.semiMajorAxis;
377  flattening = om.flattening;
378  es = om.es;
379  es_OVER_2 = om.es_OVER_2;
380  OMerc_A = om.OMerc_A;
381  OMerc_B = om.OMerc_B;
382  OMerc_E = om.OMerc_E;
383  OMerc_gamma = om.OMerc_gamma;
384  OMerc_azimuth = om.OMerc_azimuth;
385  OMerc_Origin_Long = om.OMerc_Origin_Long;
386  cos_gamma = om.cos_gamma;
387  sin_gamma = om.sin_gamma;
388  sin_azimuth = om.sin_azimuth;
389  cos_azimuth = om.cos_azimuth;
390  A_over_B = om.A_over_B;
391  B_over_A = om.B_over_A;
392  OMerc_u = om.OMerc_u;
393  OMerc_Origin_Lat = om.OMerc_Origin_Lat;
394  OMerc_Lon_1 = om.OMerc_Lon_1;
395  OMerc_Lat_1 = om.OMerc_Lat_1;
396  OMerc_Lon_2 = om.OMerc_Lon_2;
397  OMerc_Lat_2 = om.OMerc_Lat_2;
398  OMerc_False_Easting = om.OMerc_False_Easting;
399  OMerc_False_Northing = om.OMerc_False_Northing;
400  OMerc_Scale_Factor = om.OMerc_Scale_Factor;
401  OMerc_Delta_Northing = om.OMerc_Delta_Northing;
402  OMerc_Delta_Easting = om.OMerc_Delta_Easting;
403  }
404 
405  return *this;
406 }
407 
408 
409 ObliqueMercatorParameters* ObliqueMercator::getParameters() const
410 {
411 /*
412  * The function getParameters returns the current ellipsoid
413  * parameters and Oblique Mercator projection parameters.
414  *
415  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
416  * ellipsoidFlattening : Flattening of ellipsoid (output)
417  * originLatitude : Latitude, in radians, at which the (output)
418  * point scale factor is 1.0
419  * longitude1 : Longitude, in radians, of first point lying on
420  * central line (output)
421  * latitude1 : Latitude, in radians, of first point lying on
422  * central line (output)
423  * longitude2 : Longitude, in radians, of second point lying on
424  * central line (output)
425  * latitude2 : Latitude, in radians, of second point lying on
426  * central line (output)
427  * falseEasting : A coordinate value, in meters, assigned to the
428  * central meridian of the projection (output)
429  * falseNorthing : A coordinate value, in meters, assigned to the
430  * origin latitude of the projection (output)
431  * scaleFactor : Multiplier which reduces distances in the
432  * projection to the actual distance on the
433  * ellipsoid (output)
434  */
435 
436  return new ObliqueMercatorParameters( CoordinateType::obliqueMercator, OMerc_Origin_Lat, OMerc_Lon_1, OMerc_Lat_1, OMerc_Lon_2, OMerc_Lat_2, OMerc_False_Easting, OMerc_False_Northing, OMerc_Scale_Factor );
437 }
438 
439 
440 MSP::CCS::MapProjectionCoordinates* ObliqueMercator::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
441 {
442 /*
443  * The function convertFromGeodetic converts geodetic (latitude and
444  * longitude) coordinates to Oblique Mercator projection (easting and
445  * northing) coordinates, according to the current ellipsoid and Oblique Mercator
446  * projection parameters. If any errors occur, an exception is thrown with a description
447  * of the error.
448  *
449  * longitude : Longitude (lambda), in radians (input)
450  * latitude : Latitude (phi), in radians (input)
451  * easting : Easting (X), in meters (output)
452  * northing : Northing (Y), in meters (output)
453  */
454 
455  double dlam, B_dlam, cos_B_dlam;
456  double t, S, T, V, U;
457  double Q, Q_inv;
458  /* Coordinate axes defined with respect to the azimuth of the center line */
459  /* Natural origin*/
460  double v = 0;
461  double u = 0;
462 
463  double longitude = geodeticCoordinates->longitude();
464  double latitude = geodeticCoordinates->latitude();
465 
466  if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
467  { /* Latitude out of range */
468  throw CoordinateConversionException( ErrorMessages::latitude );
469  }
470  if ((longitude < -PI) || (longitude > TWO_PI))
471  { /* Longitude out of range */
472  throw CoordinateConversionException( ErrorMessages::longitude );
473  }
474 
475  dlam = longitude - OMerc_Origin_Long;
476 
477  char warning[256];
478  warning[0] = '\0';
479  if (fabs(dlam) >= PI_OVER_2)
480  { /* Distortion will result if Longitude is 90 degrees or more from the Central Meridian */
481  strcat( warning, MSP::CCS::WarningMessages::longitude );
482  }
483 
484  if (dlam > PI)
485  {
486  dlam -= TWO_PI;
487  }
488  if (dlam < -PI)
489  {
490  dlam += TWO_PI;
491  }
492 
493  if (fabs(fabs(latitude) - PI_OVER_2) > 1.0e-10)
494  {
495  t = omercT(latitude, es * sin(latitude), es_OVER_2);
496  Q = OMerc_E / pow(t, OMerc_B);
497  Q_inv = 1.0 / Q;
498  S = (Q - Q_inv) / 2.0;
499  T = (Q + Q_inv) / 2.0;
500  B_dlam = OMerc_B * dlam;
501  V = sin(B_dlam);
502  U = ((-1.0 * V * cos_gamma) + (S * sin_gamma)) / T;
503  if (fabs(fabs(U) - 1.0) < 1.0e-10)
504  { /* Point projects into infinity */
505  throw CoordinateConversionException( ErrorMessages::longitude );
506  }
507  else
508  {
509  v = A_over_B * log((1.0 - U) / (1.0 + U)) / 2.0;
510  cos_B_dlam = cos(B_dlam);
511  if (fabs(cos_B_dlam) < 1.0e-10)
512  u = OMerc_A * B_dlam;
513  // Check for longitude span > 90 degrees
514  else if (fabs(B_dlam) > PI_OVER_2)
515  {
516  double temp = atan(((S * cos_gamma) + (V * sin_gamma)) / cos_B_dlam);
517  if (temp < 0.0)
518  u = A_over_B * (temp + PI);
519  else
520  u = A_over_B * (temp - PI);
521  }
522  else
523  u = A_over_B * atan(((S * cos_gamma) + (V * sin_gamma)) / cos_B_dlam);
524  }
525  }
526  else
527  {
528  if (latitude > 0.0)
529  v = A_over_B * log(tan(PI_OVER_4 - (OMerc_gamma / 2.0)));
530  else
531  v = A_over_B * log(tan(PI_OVER_4 + (OMerc_gamma / 2.0)));
532  u = A_over_B * latitude;
533  }
534 
535 
536  u = u - OMerc_u;
537 
538  double easting = OMerc_False_Easting + v * cos_azimuth + u * sin_azimuth;
539  double northing = OMerc_False_Northing + u * cos_azimuth - v * sin_azimuth;
540 
541  return new MapProjectionCoordinates(
542  CoordinateType::obliqueMercator, warning, easting, northing );
543 }
544 
545 
546 MSP::CCS::GeodeticCoordinates* ObliqueMercator::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
547 {
548 /*
549  * The function convertToGeodetic converts Oblique Mercator projection
550  * (easting and northing) coordinates to geodetic (latitude and longitude)
551  * coordinates, according to the current ellipsoid and Oblique Mercator projection
552  * coordinates. If any errors occur, an exception is thrown with a description
553  * of the error.
554  *
555  * easting : Easting (X), in meters (input)
556  * northing : Northing (Y), in meters (input)
557  * longitude : Longitude (lambda), in radians (output)
558  * latitude : Latitude (phi), in radians (output)
559  */
560 
561  double dx, dy;
562  /* Coordinate axes defined with respect to the azimuth of the center line */
563  /* Natural origin*/
564  double u, v;
565  double Q_prime, Q_prime_inv;
566  double S_prime, T_prime, V_prime, U_prime;
567  double t;
568  double es_sin;
569  double u_B_over_A;
570  double phi;
571  double temp_phi = 0.0;
572  int count = 60;
573  double longitude, latitude;
574 
575  double easting = mapProjectionCoordinates->easting();
576  double northing = mapProjectionCoordinates->northing();
577 
578  if ((easting < (OMerc_False_Easting - OMerc_Delta_Easting))
579  || (easting > (OMerc_False_Easting + OMerc_Delta_Easting)))
580  { /* Easting out of range */
581  throw CoordinateConversionException( ErrorMessages::easting );
582  }
583  if ((northing < (OMerc_False_Northing - OMerc_Delta_Northing))
584  || (northing > (OMerc_False_Northing + OMerc_Delta_Northing)))
585  { /* Northing out of range */
586  throw CoordinateConversionException( ErrorMessages::northing );
587  }
588 
589  dy = northing - OMerc_False_Northing;
590  dx = easting - OMerc_False_Easting;
591  v = dx * cos_azimuth - dy * sin_azimuth;
592  u = dy * cos_azimuth + dx * sin_azimuth;
593  u = u + OMerc_u;
594  Q_prime = exp(-1.0 * (v * B_over_A ));
595  Q_prime_inv = 1.0 / Q_prime;
596  S_prime = (Q_prime - Q_prime_inv) / 2.0;
597  T_prime = (Q_prime + Q_prime_inv) / 2.0;
598  u_B_over_A = u * B_over_A;
599  V_prime = sin(u_B_over_A);
600  U_prime = (V_prime * cos_gamma + S_prime * sin_gamma) / T_prime;
601  if (fabs(fabs(U_prime) - 1.0) < 1.0e-10)
602  {
603  if (U_prime > 0)
604  latitude = PI_OVER_2;
605  else
606  latitude = -PI_OVER_2;
607  longitude = OMerc_Origin_Long;
608  }
609  else
610  {
611  t = pow(OMerc_E / sqrt((1.0 + U_prime) / (1.0 - U_prime)), 1.0 / OMerc_B);
612  phi = PI_OVER_2 - 2.0 * atan(t);
613  while (fabs(phi - temp_phi) > 1.0e-10 && count)
614  {
615  temp_phi = phi;
616  es_sin = es * sin(phi);
617  phi = PI_OVER_2 - 2.0 * atan(t * pow((1.0 - es_sin) / (1.0 + es_sin), es_OVER_2));
618  count --;
619  }
620 
621  if(!count)
622  throw CoordinateConversionException( ErrorMessages::northing );
623 
624  latitude = phi;
625  longitude = OMerc_Origin_Long - atan2((S_prime * cos_gamma - V_prime * sin_gamma), cos(u_B_over_A)) / OMerc_B;
626  }
627 
628  if (fabs(latitude) < 2.0e-7) /* force lat to 0 to avoid -0 degrees */
629  latitude = 0.0;
630  if (latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
631  latitude = PI_OVER_2;
632  else if (latitude < -PI_OVER_2)
633  latitude = -PI_OVER_2;
634 
635  if (longitude > PI)
636  longitude -= TWO_PI;
637  if (longitude < -PI)
638  longitude += TWO_PI;
639 
640  if (fabs(longitude) < 2.0e-7) /* force lon to 0 to avoid -0 degrees */
641  longitude = 0.0;
642  if (longitude > PI) /* force distorted values to 180, -180 degrees */
643  longitude = PI;
644  else if (longitude < -PI)
645  longitude = -PI;
646 
647  char warning[256];
648  warning[0] = '\0';
649  if (fabs(longitude - OMerc_Origin_Long) >= PI_OVER_2)
650  { /* Distortion results if Longitude > 90 degrees from the Central Meridian */
651  strcat( warning, MSP::CCS::WarningMessages::longitude );
652  }
653 
654  return new GeodeticCoordinates(
655  CoordinateType::geodetic, warning, longitude, latitude );
656 }
657 
658 
659 double ObliqueMercator::omercT( double lat, double e_sinlat, double e_over_2 )
660 {
661  return (tan(PI_OVER_4 - lat / 2.0)) / (pow((1 - e_sinlat) / (1 + e_sinlat), e_over_2));
662 }
663 
664 
665 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:02 for GeographicTranslator by doxygen 1.8.2