UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
VanDerGrinten.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: VAN DER GRINTEN
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates
9  * (latitude and longitude in radians) and Van Der Grinten projection
10  * coordinates (easting and northing in meters). The Van Der Grinten
11  * projection employs a spherical Earth model. The Spherical Radius
12  * used is the the radius of the sphere having the same area as the
13  * ellipsoid.
14  *
15  * ERROR HANDLING
16  *
17  * This component checks parameters for valid values. If an invalid value
18  * is found, the error code is combined with the current error code using
19  * the bitwise or. This combining allows multiple error codes to be
20  * returned. The possible error codes are:
21  *
22  * GRIN_NO_ERROR : No errors occurred in function
23  * GRIN_LAT_ERROR : Latitude outside of valid range
24  * (-90 to 90 degrees)
25  * GRIN_LON_ERROR : Longitude outside of valid range
26  * (-180 to 360 degrees)
27  * GRIN_EASTING_ERROR : Easting outside of valid range
28  * (False_Easting +/- ~20,000,000 m,
29  * depending on ellipsoid parameters)
30  * GRIN_NORTHING_ERROR : Northing outside of valid range
31  * (False_Northing +/- ~20,000,000 m,
32  * depending on ellipsoid parameters)
33  * GRIN_RADIUS_ERROR : Coordinates too far from pole,
34  * depending on ellipsoid and
35  * projection parameters
36  * GRIN_CENT_MER_ERROR : Central meridian outside of valid range
37  * (-180 to 360 degrees)
38  * GRIN_A_ERROR : Semi-major axis less than or equal to zero
39  * GRIN_INV_F_ERROR : Inverse flattening outside of valid range
40  * (250 to 350)
41  *
42  * REUSE NOTES
43  *
44  * VAN DER GRINTEN is intended for reuse by any application that performs a
45  * Van Der Grinten projection or its inverse.
46  *
47  * REFERENCES
48  *
49  * Further information on VAN DER GRINTEN can be found in the Reuse Manual.
50  *
51  * VAN DER GRINTEN originated from : U.S. Army Topographic Engineering Center
52  * Geospatial Information Division
53  * 7701 Telegraph Road
54  * Alexandria, VA 22310-3864
55  *
56  * LICENSES
57  *
58  * None apply to this component.
59  *
60  * RESTRICTIONS
61  *
62  * VAN DER GRINTEN has no restrictions.
63  *
64  * ENVIRONMENT
65  *
66  * VAN DER GRINTEN was tested and certified in the following environments:
67  *
68  * 1. Solaris 2.5 with GCC, version 2.8.1
69  * 2. Windows 95 with MS Visual C++, version 6
70  *
71  * MODIFICATIONS
72  *
73  * Date Description
74  * ---- -----------
75  * 3-1-07 Original Code
76  *
77  */
78 
79 
80 /***************************************************************************/
81 /*
82  * INCLUDES
83  */
84 
85 #include <math.h>
86 #include "VanDerGrinten.h"
87 #include "MapProjection3Parameters.h"
88 #include "MapProjectionCoordinates.h"
89 #include "GeodeticCoordinates.h"
90 #include "CoordinateConversionException.h"
91 #include "ErrorMessages.h"
92 
93 /*
94  * math.h - Standard C math library
95  * VanDerGrinten.h - Is for prototype error checking
96  * MapProjectionCoordinates.h - defines map projection coordinates
97  * GeodeticCoordinates.h - defines geodetic coordinates
98  * CoordinateConversionException.h - Exception handler
99  * ErrorMessages.h - Contains exception messages
100  */
101 
102 
103 using namespace MSP::CCS;
104 
105 
106 /***************************************************************************/
107 /* DEFINES
108  *
109  */
110 
111 const double PI = 3.14159265358979323e0; /* PI */
112 const double PI_OVER_2 = ( PI / 2.0);
113 const double MAX_LAT = ( 90 * (PI / 180.0) ); /* 90 degrees in radians */
114 const double TWO_PI = (2.0 * PI);
115 const double TWO_OVER_PI = (2.0 / PI);
116 const double PI_OVER_3 = (PI / 3.0);
117 const double ONE_THIRD = (1.0 / 3.0);
118 
119 
120 
121 /************************************************************************/
122 /* FUNCTIONS
123  *
124  */
125 
126 VanDerGrinten::VanDerGrinten( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double falseEasting, double falseNorthing ) :
127  CoordinateSystem(),
128  es2( 0.0066943799901413800 ),
129  es4( 4.4814723452405e-005 ),
130  es6( 3.0000678794350e-007 ),
131  Ra( 6371007.1810824 ),
132  PI_Ra( 20015109.356056 ),
133  Grin_Origin_Long( 0.0 ),
134  Grin_False_Easting( 0.0 ),
135  Grin_False_Northing( 0.0 )
136 {
137 /*
138  * The constructor receives the ellipsoid parameters and
139  * projection parameters as inputs, and sets the corresponding state
140  * variables. If any errors occur, an exception is thrown with a description
141  * of the error.
142  *
143  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
144  * ellipsoidFlattening : Flattening of ellipsoid (input)
145  * centralMeridian : Longitude in radians at the center of (input)
146  * the projection
147  * falseEasting : A coordinate value in meters assigned to the
148  * central meridian of the projection. (input)
149  * falseNorthing : A coordinate value in meters assigned to the
150  * origin latitude of the projection (input)
151  */
152 
153  double inv_f = 1 / ellipsoidFlattening;
154 
155  if (ellipsoidSemiMajorAxis <= 0.0)
156  { /* Semi-major axis must be greater than zero */
157  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
158  }
159  if ((inv_f < 250) || (inv_f > 350))
160  { /* Inverse flattening must be between 250 and 350 */
161  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
162  }
163  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
164  { /* origin longitude out of range */
165  throw CoordinateConversionException( ErrorMessages::centralMeridian );
166  }
167 
168  semiMajorAxis = ellipsoidSemiMajorAxis;
169  flattening = ellipsoidFlattening;
170 
171  es2 = 2 * flattening - flattening * flattening;
172  es4 = es2 * es2;
173  es6 = es4 * es2;
174  /* spherical radius */
175  Ra = semiMajorAxis * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 /3024.0);
176  PI_Ra = PI * Ra;
177  if (centralMeridian > PI)
178  centralMeridian -= TWO_PI;
179  Grin_Origin_Long = centralMeridian;
180  Grin_False_Easting = falseEasting;
181  Grin_False_Northing = falseNorthing;
182 }
183 
184 
185 VanDerGrinten::VanDerGrinten( const VanDerGrinten &v )
186 {
187  semiMajorAxis = v.semiMajorAxis;
188  flattening = v.flattening;
189  es2 = v.es2;
190  es4 = v.es4;
191  es6 = v.es6;
192  Ra = v.Ra;
193  PI_Ra = v.PI_Ra;
194  Grin_Origin_Long = v.Grin_Origin_Long;
195  Grin_False_Easting = v.Grin_False_Easting;
196  Grin_False_Northing = v.Grin_False_Northing;
197 }
198 
199 
200 VanDerGrinten::~VanDerGrinten()
201 {
202 }
203 
204 
205 VanDerGrinten& VanDerGrinten::operator=( const VanDerGrinten &v )
206 {
207  if( this != &v )
208  {
209  semiMajorAxis = v.semiMajorAxis;
210  flattening = v.flattening;
211  es2 = v.es2;
212  es4 = v.es4;
213  es6 = v.es6;
214  Ra = v.Ra;
215  PI_Ra = v.PI_Ra;
216  Grin_Origin_Long = v.Grin_Origin_Long;
217  Grin_False_Easting = v.Grin_False_Easting;
218  Grin_False_Northing = v.Grin_False_Northing;
219  }
220 
221  return *this;
222 }
223 
224 
225 MapProjection3Parameters* VanDerGrinten::getParameters() const
226 {
227 /*
228  * The function getParameters returns the current ellipsoid
229  * parameters, and Van Der Grinten projection parameters.
230  *
231  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
232  * ellipsoidFlattening : Flattening of ellipsoid (output)
233  * centralMeridian : Longitude in radians at the center of (output)
234  * the projection
235  * falseEasting : A coordinate value in meters assigned to the
236  * central meridian of the projection. (output)
237  * falseNorthing : A coordinate value in meters assigned to the
238  * origin latitude of the projection (output)
239  */
240 
241  return new MapProjection3Parameters( CoordinateType::vanDerGrinten, Grin_Origin_Long, Grin_False_Easting, Grin_False_Northing );
242 }
243 
244 
245 MSP::CCS::MapProjectionCoordinates* VanDerGrinten::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
246 {
247 /*
248  * The function convertFromGeodetic converts geodetic (latitude and
249  * longitude) coordinates to Van Der Grinten projection (easting and northing)
250  * coordinates, according to the current ellipsoid and Van Der Grinten projection
251  * parameters. If any errors occur, an exception is thrown with a description
252  * of the error.
253  *
254  * longitude : Longitude (lambda) in radians (input)
255  * latitude : Latitude (phi) in radians (input)
256  * easting : Easting (X) in meters (output)
257  * northing : Northing (Y) in meters (output)
258  */
259 
260  double dlam; /* Longitude - Central Meridan */
261  double aa, aasqr;
262  double gg;
263  double pp, ppsqr;
264  double gg_MINUS_ppsqr, ppsqr_PLUS_aasqr;
265  double in_theta;
266  double theta;
267  double sin_theta, cos_theta;
268  double qq;
269  double easting, northing;
270 
271  double longitude = geodeticCoordinates->longitude();
272  double latitude = geodeticCoordinates->latitude();
273 
274  if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
275  { /* Latitude out of range */
276  throw CoordinateConversionException( ErrorMessages::latitude );
277  }
278  if ((longitude < -PI) || (longitude > TWO_PI))
279  { /* Longitude out of range */
280  throw CoordinateConversionException( ErrorMessages::longitude );
281  }
282 
283  dlam = longitude - Grin_Origin_Long;
284  if (dlam > PI)
285  {
286  dlam -= TWO_PI;
287  }
288  if (dlam < -PI)
289  {
290  dlam += TWO_PI;
291  }
292 
293  if (latitude == 0.0)
294  {
295  easting = Ra * dlam + Grin_False_Easting;
296  northing = 0.0;
297  }
298  else if (dlam == 0.0 || floatEq(latitude,MAX_LAT,.00001) || floatEq(latitude,-MAX_LAT,.00001))
299  {
300  in_theta = fabs(TWO_OVER_PI * latitude);
301 
302  if (in_theta > 1.0)
303  in_theta = 1.0;
304  else if (in_theta < -1.0)
305  in_theta = -1.0;
306 
307  theta = asin(in_theta);
308  easting = 0.0;
309  northing = PI_Ra * tan(theta / 2) + Grin_False_Northing;
310  if (latitude < 0.0)
311  northing *= -1.0;
312  }
313  else
314  {
315  aa = 0.5 * fabs(PI / dlam - dlam / PI);
316  in_theta = fabs(TWO_OVER_PI * latitude);
317 
318  if (in_theta > 1.0)
319  in_theta = 1.0;
320  else if (in_theta < -1.0)
321  in_theta = -1.0;
322 
323  theta = asin(in_theta);
324  sin_theta = sin(theta);
325  cos_theta = cos(theta);
326  gg = cos_theta / (sin_theta + cos_theta - 1);
327  pp = gg * (2 / sin_theta - 1);
328  aasqr = aa * aa;
329  ppsqr = pp * pp;
330  gg_MINUS_ppsqr = gg - ppsqr;
331  ppsqr_PLUS_aasqr = ppsqr + aasqr;
332  qq = aasqr + gg;
333  easting = PI_Ra * (aa * (gg_MINUS_ppsqr) +
334  sqrt(aasqr * (gg_MINUS_ppsqr) * (gg_MINUS_ppsqr) -
335  (ppsqr_PLUS_aasqr) * (gg * gg - ppsqr))) /
336  (ppsqr_PLUS_aasqr) + Grin_False_Easting;
337  if (dlam < 0.0)
338  easting *= -1.0;
339  northing = PI_Ra * (pp * qq - aa * sqrt ((aasqr + 1) * (ppsqr_PLUS_aasqr) - qq * qq)) /
340  (ppsqr_PLUS_aasqr) + Grin_False_Northing;
341  if (latitude < 0.0)
342  northing *= -1.0;
343  }
344 
345  return new MapProjectionCoordinates( CoordinateType::vanDerGrinten, easting, northing );
346 }
347 
348 
349 MSP::CCS::GeodeticCoordinates* VanDerGrinten::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
350 {
351 /*
352  * The function convertToGeodetic converts Grinten projection
353  * (easting and northing) coordinates to geodetic (latitude and longitude)
354  * coordinates, according to the current ellipsoid and Grinten projection
355  * coordinates. If any errors occur, an exception is thrown with a description
356  * of the error.
357  *
358  * easting : Easting (X) in meters (input)
359  * northing : Northing (Y) in meters (input)
360  * longitude : Longitude (lambda) in radians (output)
361  * latitude : Latitude (phi) in radians (output)
362  */
363 
364  double dx, dy;
365  double xx, xxsqr;
366  double yy, yysqr, two_yysqr;
367  double xxsqr_PLUS_yysqr;
368  double c1;
369  double c2;
370  double c3, c3sqr;
371  double c2_OVER_3c3;
372  double dd;
373  double a1;
374  double m1;
375  double i;
376  double theta1;
377  double temp;
378  const double epsilon = 1.0e-2;
379  double delta = PI_Ra + epsilon;
380  double longitude, latitude;
381 
382  double easting = mapProjectionCoordinates->easting();
383  double northing = mapProjectionCoordinates->northing();
384 
385  if ((easting > (Grin_False_Easting + delta)) ||
386  (easting < (Grin_False_Easting - delta)))
387  { /* Easting out of range */
388  throw CoordinateConversionException( ErrorMessages::easting );
389  }
390  if ((northing > (Grin_False_Northing + delta)) ||
391  (northing < (Grin_False_Northing - delta)))
392  { /* Northing out of range */
393  throw CoordinateConversionException( ErrorMessages::northing );
394  }
395 
396  temp = sqrt(easting * easting + northing * northing);
397 
398  if ((temp > (Grin_False_Easting + PI_Ra + epsilon)) ||
399  (temp > (Grin_False_Northing + PI_Ra + epsilon)) ||
400  (temp < (Grin_False_Easting - PI_Ra - epsilon)) ||
401  (temp < (Grin_False_Northing - PI_Ra - epsilon)))
402  { /* Point is outside of projection area */
403  throw CoordinateConversionException( ErrorMessages::radius );
404  }
405 
406  dy = northing - Grin_False_Northing;
407  dx = easting - Grin_False_Easting;
408  xx = dx / PI_Ra;
409  yy = dy / PI_Ra;
410  xxsqr = xx * xx;
411  yysqr = yy * yy;
412  xxsqr_PLUS_yysqr = xxsqr + yysqr;
413  two_yysqr = 2 * yysqr;
414 
415  if (northing == 0.0)
416  latitude = 0.0;
417  else
418  {
419  c1 = - fabs(yy) * (1 + xxsqr_PLUS_yysqr);
420  c2 = c1 - two_yysqr + xxsqr;
421  c3 = - 2 * c1 + 1 + two_yysqr + (xxsqr_PLUS_yysqr) * (xxsqr_PLUS_yysqr);
422  c2_OVER_3c3 = c2 / (3.0 * c3);
423  c3sqr = c3 * c3;
424  dd = yysqr / c3 + ((2 * c2 * c2 * c2) / (c3sqr * c3) - (9 * c1 * c2) / (c3sqr)) / 27;
425  a1 = (c1 - c2 * c2_OVER_3c3) /c3;
426  m1 = 2 * sqrt(-ONE_THIRD * a1);
427  i = 3 * dd/ (a1 * m1);
428  if ((i > 1.0)||(i < -1.0))
429  latitude = MAX_LAT;
430  else
431  {
432  theta1 = ONE_THIRD * acos(3 * dd / (a1 * m1));
433  latitude = PI * (-m1 * cos(theta1 + PI_OVER_3) - c2_OVER_3c3);
434  }
435  }
436  if (northing < 0.0)
437  latitude *= -1.0;
438 
439  if (xx == 0.0)
440  longitude = Grin_Origin_Long;
441  else
442  {
443  longitude = PI * (xxsqr_PLUS_yysqr - 1 +
444  sqrt(1 + (2 * xxsqr - two_yysqr) + (xxsqr_PLUS_yysqr) * (xxsqr_PLUS_yysqr))) /
445  (2 * xx) + Grin_Origin_Long;
446  }
447  if (latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
448  latitude = PI_OVER_2;
449  else if (latitude < -PI_OVER_2)
450  latitude = -PI_OVER_2;
451 
452  if (longitude > PI)
453  longitude -= TWO_PI;
454  if (longitude < -PI)
455  longitude += TWO_PI;
456 
457  if (longitude > PI) /* force distorted values to 180, -180 degrees */
458  longitude = PI;
459  else if (longitude < -PI)
460  longitude = -PI;
461 
462  return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
463 }
464 
465 
466 double VanDerGrinten::floatEq( double x, double v, double epsilon )
467 {
468  return (((v - epsilon) < x) && (x < (v + epsilon)));
469 }
470 
471 
472 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:02 for GeographicTranslator by doxygen 1.8.2