UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
UTM.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: UTM
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between geodetic coordinates
9  * (latitude and longitudes) and Universal Transverse Mercator (UTM)
10  * projection (zone, hemisphere, easting, and northing) coordinates.
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  * UTM_NO_ERROR : No errors occurred in function
20  * UTM_LAT_ERROR : Latitude outside of valid range
21  * (-80.5 to 84.5 degrees)
22  * UTM_LON_ERROR : Longitude outside of valid range
23  * (-180 to 360 degrees)
24  * UTM_EASTING_ERROR : Easting outside of valid range
25  * (100,000 to 900,000 meters)
26  * UTM_NORTHING_ERROR : Northing outside of valid range
27  * (0 to 10,000,000 meters)
28  * UTM_ZONE_ERROR : Zone outside of valid range (1 to 60)
29  * UTM_HEMISPHERE_ERROR : Invalid hemisphere ('N' or 'S')
30  * UTM_ZONE_OVERRIDE_ERROR: Zone outside of valid range
31  * (1 to 60) and within 1 of 'natural' zone
32  * UTM_A_ERROR : Semi-major axis less than or equal to zero
33  * UTM_INV_F_ERROR : Inverse flattening outside of valid range
34  * (250 to 350)
35  *
36  * REUSE NOTES
37  *
38  * UTM is intended for reuse by any application that performs a Universal
39  * Transverse Mercator (UTM) projection or its inverse.
40  *
41  * REFERENCES
42  *
43  * Further information on UTM can be found in the Reuse Manual.
44  *
45  * UTM originated from : U.S. Army Topographic Engineering Center
46  * Geospatial Information Division
47  * 7701 Telegraph Road
48  * Alexandria, VA 22310-3864
49  *
50  * LICENSES
51  *
52  * None apply to this component.
53  *
54  * RESTRICTIONS
55  *
56  * UTM has no restrictions.
57  *
58  * ENVIRONMENT
59  *
60  * UTM was tested and certified in the following environments:
61  *
62  * 1. Solaris 2.5 with GCC, version 2.8.1
63  * 2. MSDOS with MS Visual C++, version 6
64  *
65  * MODIFICATIONS
66  *
67  * Date Description
68  * ---- -----------
69  * 2-27-07 Original C++ Code
70  * 7/18/2010 NGL BAEts27181 Add logic for special zones over southern
71  * Norway and Svalbard. This was removed with BAEts24468
72  * and we are being asked to put it back in.
73  * 5-09-11 KNL BAEts28908, add default constructor.
74  * 1-16-16 AL MSP_DR30125 added ellipsoid code into constructor
75  * and passes ellipsoid code to TransMercator
76  * 1-21-16 KC BAE_MSP00030211, removed the shift from longitude.
77  * Shift is applied when determining the zone.
78  *
79  */
80 
81 
82 /***************************************************************************/
83 /*
84  * INCLUDES
85  */
86 
87 #include <math.h>
88 #include "UTM.h"
89 #include "TransverseMercator.h"
90 #include "UTMParameters.h"
91 #include "MapProjectionCoordinates.h"
92 #include "UTMCoordinates.h"
93 #include "GeodeticCoordinates.h"
94 #include "CoordinateConversionException.h"
95 #include "ErrorMessages.h"
96 
97 /*
98  * UTM.h - Defines the function prototypes for the utm module.
99  * TransverseMercator.h - Is used to convert transverse mercator coordinates
100  * MapProjectionCoordinates.h - defines map projection coordinates
101  * UTMCoordinates.h - defines UTM coordinates
102  * GeodeticCoordinates.h - defines geodetic coordinates
103  * CoordinateConversionException.h - Exception handler
104  * ErrorMessages.h - Contains exception messages
105  */
106 
107 using namespace MSP::CCS;
108 
109 /***************************************************************************/
110 /*
111  * DEFINES
112  */
113 
114 #define PI 3.14159265358979323e0
115 #define PI_OVER_180 (3.14159265358979323e0 / 180.0)
116 #define MIN_LAT ((-80.5 * PI) / 180.0) /* -80.5 degrees in radians */
117 #define MAX_LAT (( 84.5 * PI) / 180.0) /* 84.5 degrees in radians */
118 #define MIN_EASTING 100000.0
119 #define MAX_EASTING 900000.0
120 #define MIN_NORTHING 0.0
121 #define MAX_NORTHING 10000000.0
122 
123 #define EPSILON 1.75e-7 /* approx 1.0e-5 deg (~1 meter) in radians */
124 
125 /************************************************************************/
126 /* FUNCTIONS
127  *
128  */
129 
130 UTM::UTM()
131 {
132 /*
133  * The default constructor
134  */
135  // DR30125 default WGS-84
136  strcpy( ellipsCode, "WE" );
137  double ellipsoidSemiMajorAxis = 6378137.0;
138  double ellipsoidFlattening = 1 / 298.257223563;
139  double inv_f = 1 / ellipsoidFlattening;
140 
141  semiMajorAxis = ellipsoidSemiMajorAxis;
142  flattening = ellipsoidFlattening;
143  UTM_Override = 0;
144 
145  double centralMeridian;
146  double originLatitude = 0;
147  double falseEasting = 500000;
148  double falseNorthing = 0;
149  double scale = 0.9996;
150 
151  for(int zone = 1; zone <= 60; zone++)
152  {
153  if (zone >= 31)
154  centralMeridian = ((6 * zone - 183) * PI_OVER_180);
155  else
156  centralMeridian = ((6 * zone + 177) * PI_OVER_180);
157 
158  transverseMercatorMap[zone] = new TransverseMercator(
159  semiMajorAxis, flattening, centralMeridian, originLatitude,
160  falseEasting, falseNorthing, scale, ellipsCode);
161  }
162 }
163 
164 UTM::UTM(
165  double ellipsoidSemiMajorAxis,
166  double ellipsoidFlattening,
167  char *ellipsoidCode,
168  long override
169  ) :
170  UTM_Override( 0 )
171 {
172 /*
173  * The constructor receives the ellipsoid parameters and
174  * UTM zone override parameter as inputs, and sets the corresponding state
175  * variables. If any errors occur, an exception is thrown with a description
176  * of the error.
177  *
178  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
179  * ellipsoidFlattening : Flattening of ellipsoid (input)
180  * override : UTM override zone, 0 indicates no override (input)
181  */
182  /* get ellipsoid Code */
183  strcpy( ellipsCode, ellipsoidCode );
184  double inv_f = 1 / ellipsoidFlattening;
185 
186  if (ellipsoidSemiMajorAxis <= 0.0)
187  { /* Semi-major axis must be greater than zero */
188  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
189  }
190  if ((inv_f < 250) || (inv_f > 350))
191  { /* Inverse flattening must be between 250 and 350 */
192  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
193  }
194  if ((override < 0) || (override > 60))
195  {
196  throw CoordinateConversionException( ErrorMessages::zoneOverride );
197  }
198 
199  semiMajorAxis = ellipsoidSemiMajorAxis;
200  flattening = ellipsoidFlattening;
201 
202  UTM_Override = override;
203 
204  double centralMeridian;
205  double originLatitude = 0;
206  double falseEasting = 500000;
207  double falseNorthing = 0;
208  double scale = 0.9996;
209 
210  for(int zone = 1; zone <= 60; zone++)
211  {
212  if (zone >= 31)
213  centralMeridian = ((6 * zone - 183) * PI_OVER_180);
214  else
215  centralMeridian = ((6 * zone + 177) * PI_OVER_180);
216 
217  transverseMercatorMap[zone] = new TransverseMercator(
218  semiMajorAxis, flattening, centralMeridian, originLatitude,
219  falseEasting, falseNorthing, scale, ellipsCode);
220  }
221 }
222 
223 
224 UTM::UTM( const UTM &u )
225 {
226  int zone = 1;
227  std::map< int, TransverseMercator* > tempTransverseMercatorMap =
228  u.transverseMercatorMap;
229  std::map< int, TransverseMercator* >::iterator _iterator =
230  tempTransverseMercatorMap.begin();
231  while( _iterator != tempTransverseMercatorMap.end() && zone <= 60 )
232  {
233  transverseMercatorMap[zone] = new TransverseMercator( *_iterator->second );
234  zone++;
235  }
236 
237  semiMajorAxis = u.semiMajorAxis;
238  flattening = u.flattening;
239  UTM_Override = u.UTM_Override;
240 }
241 
242 
243 UTM::~UTM()
244 {
245  while( transverseMercatorMap.begin() != transverseMercatorMap.end() )
246  {
247  delete ( ( *transverseMercatorMap.begin() ).second );
248  transverseMercatorMap.erase( transverseMercatorMap.begin() );
249  }
250 }
251 
252 
253 UTM& UTM::operator=( const UTM &u )
254 {
255  if( this != &u )
256  {
257  int zone = 1;
258  std::map< int, TransverseMercator* > tempTransverseMercatorMap =
259  u.transverseMercatorMap;
260  std::map< int, TransverseMercator* >::iterator _iterator =
261  tempTransverseMercatorMap.begin();
262  while( _iterator != tempTransverseMercatorMap.end() && zone <= 60 )
263  {
264  transverseMercatorMap[zone]->operator=( *_iterator->second );
265  zone++;
266  }
267 
269  semiMajorAxis = u.semiMajorAxis;
270  flattening = u.flattening;
271  UTM_Override = u.UTM_Override;
272  }
273 
274  return *this;
275 }
276 
277 
278 UTMParameters* UTM::getParameters() const
279 {
280 /*
281  * The function getParameters returns the current ellipsoid
282  * parameters and UTM zone override parameter.
283  *
284  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
285  * ellipsoidFlattening : Flattening of ellipsoid (output)
286  * override : UTM override zone, zero indicates no override (output)
287  */
288 
289  return new UTMParameters(
290  CoordinateType::universalTransverseMercator, UTM_Override );
291 }
292 
293 
294 MSP::CCS::UTMCoordinates* UTM::convertFromGeodetic(
295  MSP::CCS::GeodeticCoordinates* geodeticCoordinates,
296  int utmZoneOverride )
297 {
298 /*
299  * The function convertFromGeodetic converts geodetic (latitude and
300  * longitude) coordinates to UTM projection (zone, hemisphere, easting and
301  * northing) coordinates according to the current ellipsoid and UTM zone
302  * override parameters. If any errors occur, an exception is thrown
303  * with a description of the error.
304  *
305  * longitude : Longitude in radians (input)
306  * latitude : Latitude in radians (input)
307  * zone : UTM zone (output)
308  * hemisphere : North or South hemisphere (output)
309  * easting : Easting (X) in meters (output)
310  * northing : Northing (Y) in meters (output)
311  */
312 
313  long Lat_Degrees;
314  long Long_Degrees;
315  long temp_zone;
316  char hemisphere;
317  double False_Northing = 0;
318 
319  double longitude = geodeticCoordinates->longitude();
320  double latitude = geodeticCoordinates->latitude();
321 
322  if ((latitude < (MIN_LAT - EPSILON)) || (latitude >= (MAX_LAT + EPSILON)))
323  { /* latitude out of range */
324  throw CoordinateConversionException( ErrorMessages::latitude );
325  }
326  if ((longitude < (-PI - EPSILON)) || (longitude > (2*PI + EPSILON)))
327  { /* longitude out of range */
328  throw CoordinateConversionException( ErrorMessages::longitude );
329  }
330 
331  if((latitude > -1.0e-9) && (latitude < 0))
332  latitude = 0.0;
333 
334  if (longitude < 0)
335  longitude += (2*PI);
336 
337  Lat_Degrees = (long)(latitude * 180.0 / PI);
338  Long_Degrees = (long)(longitude * 180.0 / PI);
339 
340  if (longitude < PI)
341  temp_zone = (long)(31 + (((longitude+1.0e-10) * 180.0 / PI) / 6.0));
342  else
343  temp_zone = (long)((((longitude+1.0e-10) * 180.0 / PI) / 6.0) - 29);
344 
345  if (temp_zone > 60)
346  temp_zone = 1;
347 
348  /* allow UTM zone override up to +/- one zone of the calculated zone */
349  if( utmZoneOverride )
350  {
351  if ((temp_zone == 1) && (utmZoneOverride == 60))
352  temp_zone = utmZoneOverride;
353  else if ((temp_zone == 60) && (utmZoneOverride == 1))
354  temp_zone = utmZoneOverride;
355  else if (((temp_zone-1) <= utmZoneOverride) &&
356  (utmZoneOverride <= (temp_zone+1)))
357  temp_zone = utmZoneOverride;
358  else
359  throw CoordinateConversionException( ErrorMessages::zoneOverride );
360  }
361  else if( UTM_Override )
362  {
363  if ((temp_zone == 1) && (UTM_Override == 60))
364  temp_zone = UTM_Override;
365  else if ((temp_zone == 60) && (UTM_Override == 1))
366  temp_zone = UTM_Override;
367  else if (((temp_zone-1) <= UTM_Override) &&
368  (UTM_Override <= (temp_zone+1)))
369  temp_zone = UTM_Override;
370  else
371  throw CoordinateConversionException( ErrorMessages::zoneOverride );
372  }
373  else /* not UTM zone override */
374  {
375  /* check for special zone cases over southern Norway and Svalbard */
376  if ((Lat_Degrees > 55) && (Lat_Degrees < 64) && (Long_Degrees > -1)
377  && (Long_Degrees < 3))
378  temp_zone = 31;
379  if ((Lat_Degrees > 55) && (Lat_Degrees < 64) && (Long_Degrees > 2)
380  && (Long_Degrees < 12))
381  temp_zone = 32;
382  if ((Lat_Degrees > 71) && (Long_Degrees > -1) && (Long_Degrees < 9))
383  temp_zone = 31;
384  if ((Lat_Degrees > 71) && (Long_Degrees > 8) && (Long_Degrees < 21))
385  temp_zone = 33;
386  if ((Lat_Degrees > 71) && (Long_Degrees > 20) && (Long_Degrees < 33))
387  temp_zone = 35;
388  if ((Lat_Degrees > 71) && (Long_Degrees > 32) && (Long_Degrees < 42))
389  temp_zone = 37;
390  }
391 
392  TransverseMercator *transverseMercator = transverseMercatorMap[temp_zone];
393 
394  if (latitude < 0)
395  {
396  False_Northing = 10000000;
397  hemisphere = 'S';
398  }
399  else
400  hemisphere = 'N';
401 
402  GeodeticCoordinates tempGeodeticCoordinates(
403  CoordinateType::geodetic, longitude, latitude );
404  MapProjectionCoordinates* transverseMercatorCoordinates =
405  transverseMercator->convertFromGeodetic( &tempGeodeticCoordinates );
406  double easting = transverseMercatorCoordinates->easting();
407  double northing = transverseMercatorCoordinates->northing() + False_Northing;
408 
409  if ((easting < MIN_EASTING) || (easting > MAX_EASTING))
410  {
411  delete transverseMercatorCoordinates;
412  throw CoordinateConversionException( ErrorMessages::easting );
413  }
414 
415  if ((northing < MIN_NORTHING) || (northing > MAX_NORTHING))
416  {
417  delete transverseMercatorCoordinates;
418  throw CoordinateConversionException( ErrorMessages::northing );
419  }
420 
421  delete transverseMercatorCoordinates;
422 
423  return new UTMCoordinates(
424  CoordinateType::universalTransverseMercator,
425  temp_zone, hemisphere, easting, northing );
426 }
427 
428 
429 MSP::CCS::GeodeticCoordinates* UTM::convertToGeodetic(
430  MSP::CCS::UTMCoordinates* utmCoordinates )
431 {
432 /*
433  * The function convertToGeodetic converts UTM projection (zone,
434  * hemisphere, easting and northing) coordinates to geodetic(latitude
435  * and longitude) coordinates, according to the current ellipsoid
436  * parameters. If any errors occur, an exception is thrown with a description
437  * of the error.
438  *
439  * zone : UTM zone (input)
440  * hemisphere : North or South hemisphere (input)
441  * easting : Easting (X) in meters (input)
442  * northing : Northing (Y) in meters (input)
443  * longitude : Longitude in radians (output)
444  * latitude : Latitude in radians (output)
445  */
446 
447  double False_Northing = 0;
448 
449  long zone = utmCoordinates->zone();
450  char hemisphere = utmCoordinates->hemisphere();
451  double easting = utmCoordinates->easting();
452  double northing = utmCoordinates->northing();
453 
454  if ((zone < 1) || (zone > 60))
455  throw CoordinateConversionException( ErrorMessages::zone );
456  if ((hemisphere != 'S') && (hemisphere != 'N'))
457  throw CoordinateConversionException( ErrorMessages::hemisphere );
458  if ((easting < MIN_EASTING) || (easting > MAX_EASTING))
459  throw CoordinateConversionException( ErrorMessages::easting );
460  if ((northing < MIN_NORTHING) || (northing > MAX_NORTHING))
461  throw CoordinateConversionException( ErrorMessages::northing );
462 
463  TransverseMercator *transverseMercator = transverseMercatorMap[zone];
464 
465  if (hemisphere == 'S')
466  False_Northing = 10000000;
467 
468  MapProjectionCoordinates transverseMercatorCoordinates(
469  CoordinateType::transverseMercator, easting, northing - False_Northing );
470  GeodeticCoordinates* geodeticCoordinates =
471  transverseMercator->convertToGeodetic( &transverseMercatorCoordinates );
472  geodeticCoordinates->setWarningMessage("");
473 
474  double latitude = geodeticCoordinates->latitude();
475 
476  if ((latitude < (MIN_LAT - EPSILON)) || (latitude >= (MAX_LAT + EPSILON)))
477  { /* latitude out of range */
478  delete geodeticCoordinates;
479  throw CoordinateConversionException( ErrorMessages::northing );
480  }
481 
482  return geodeticCoordinates;
483 }
484 
485 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:02 for GeographicTranslator by doxygen 1.8.2