00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082 #include <math.h>
00083 #include "UTM.h"
00084 #include "TransverseMercator.h"
00085 #include "UTMParameters.h"
00086 #include "MapProjectionCoordinates.h"
00087 #include "UTMCoordinates.h"
00088 #include "GeodeticCoordinates.h"
00089 #include "CoordinateConversionException.h"
00090 #include "ErrorMessages.h"
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103 using namespace MSP::CCS;
00104
00105
00106
00107
00108
00109
00110
00111 #define EPSILON 1.75e-7
00112
00113 const double PI = 3.14159265358979323e0;
00114 const double PI_OVER_180 = 3.14159265358979323e0 / 180.0;
00115 const double MIN_LAT = ( (-80.5 * PI) / 180.0 );
00116 const double MAX_LAT = ( (84.5 * PI) / 180.0 );
00117 const double MIN_EASTING = 100000.0;
00118 const double MAX_EASTING = 900000.0;
00119 const double MIN_NORTHING = 0.0;
00120 const double MAX_NORTHING = 10000000.0;
00121
00122
00123
00124
00125
00126
00127
00128 UTM::UTM( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, long override ) :
00129 CoordinateSystem(),
00130 UTM_Override( 0 )
00131 {
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144 double inv_f = 1 / ellipsoidFlattening;
00145 char errorStatus[500] = "";
00146
00147 if (ellipsoidSemiMajorAxis <= 0.0)
00148 {
00149 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00150 }
00151 if ((inv_f < 250) || (inv_f > 350))
00152 {
00153 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00154 }
00155 if ((override < 0) || (override > 60))
00156 {
00157 strcat( errorStatus, ErrorMessages::zoneOverride );
00158 }
00159
00160 if( strlen( errorStatus ) > 0)
00161 throw CoordinateConversionException( errorStatus );
00162
00163 semiMajorAxis = ellipsoidSemiMajorAxis;
00164 flattening = ellipsoidFlattening;
00165
00166 UTM_Override = override;
00167
00168 double centralMeridian;
00169 double originLatitude = 0;
00170 double falseEasting = 500000;
00171 double falseNorthing = 0;
00172 double scale = 0.9996;
00173
00174 for(int zone = 1; zone <= 60; zone++)
00175 {
00176 if (zone >= 31)
00177 centralMeridian = ((6 * zone - 183) * PI_OVER_180);
00178 else
00179 centralMeridian = ((6 * zone + 177) * PI_OVER_180);
00180
00181 transverseMercatorMap[zone] = new TransverseMercator(semiMajorAxis, flattening, centralMeridian, originLatitude, falseEasting, falseNorthing, scale);
00182 }
00183 }
00184
00185
00186 UTM::UTM( const UTM &u )
00187 {
00188 int zone = 1;
00189 std::map< int, TransverseMercator* > tempTransverseMercatorMap = u.transverseMercatorMap;
00190 std::map< int, TransverseMercator* >::iterator _iterator = tempTransverseMercatorMap.begin();
00191 while( _iterator != tempTransverseMercatorMap.end() && zone <= 60 )
00192 {
00193 transverseMercatorMap[zone] = new TransverseMercator( *_iterator->second );
00194 zone++;
00195 }
00196
00197 semiMajorAxis = u.semiMajorAxis;
00198 flattening = u.flattening;
00199 UTM_Override = u.UTM_Override;
00200 }
00201
00202
00203 UTM::~UTM()
00204 {
00205 while( transverseMercatorMap.begin() != transverseMercatorMap.end() )
00206 {
00207 delete ( ( *transverseMercatorMap.begin() ).second );
00208 transverseMercatorMap.erase( transverseMercatorMap.begin() );
00209 }
00210 }
00211
00212
00213 UTM& UTM::operator=( const UTM &u )
00214 {
00215 if( this != &u )
00216 {
00217 int zone = 1;
00218 std::map< int, TransverseMercator* > tempTransverseMercatorMap = u.transverseMercatorMap;
00219 std::map< int, TransverseMercator* >::iterator _iterator = tempTransverseMercatorMap.begin();
00220 while( _iterator != tempTransverseMercatorMap.end() && zone <= 60 )
00221 {
00222 transverseMercatorMap[zone]->operator=( *_iterator->second );
00223 zone++;
00224 }
00225
00227 semiMajorAxis = u.semiMajorAxis;
00228 flattening = u.flattening;
00229 UTM_Override = u.UTM_Override;
00230 }
00231
00232 return *this;
00233 }
00234
00235
00236 UTMParameters* UTM::getParameters() const
00237 {
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247 return new UTMParameters( CoordinateType::universalTransverseMercator, UTM_Override );
00248 }
00249
00250
00251 MSP::CCS::UTMCoordinates* UTM::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00252 {
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268 long Lat_Degrees;
00269 long Long_Degrees;
00270 long temp_zone;
00271 char hemisphere;
00272 double False_Northing = 0;
00273 char errorStatus[50] = "";
00274
00275 double longitude = geodeticCoordinates->longitude();
00276 double latitude = geodeticCoordinates->latitude();
00277
00278 if ((latitude < (MIN_LAT - EPSILON)) || (latitude >= (MAX_LAT + EPSILON)))
00279 {
00280 strcat( errorStatus, ErrorMessages::latitude );
00281 }
00282 if ((longitude < -PI) || (longitude > (2*PI)))
00283 {
00284 strcat( errorStatus, ErrorMessages::longitude );
00285 }
00286
00287 if( strlen( errorStatus ) > 0)
00288 throw CoordinateConversionException( errorStatus );
00289
00290 if((latitude > -1.0e-9) && (latitude < 0))
00291 latitude = 0.0;
00292
00293 if (longitude < 0)
00294 longitude += (2*PI) + 1.0e-10;
00295
00296 Lat_Degrees = (long)(latitude * 180.0 / PI);
00297 Long_Degrees = (long)(longitude * 180.0 / PI);
00298
00299 if (longitude < PI)
00300 temp_zone = (long)(31 + ((longitude * 180.0 / PI) / 6.0));
00301 else
00302 temp_zone = (long)(((longitude * 180.0 / PI) / 6.0) - 29);
00303
00304 if (temp_zone > 60)
00305 temp_zone = 1;
00306
00307
00308 if (UTM_Override)
00309 {
00310 if ((temp_zone == 1) && (UTM_Override == 60))
00311 temp_zone = UTM_Override;
00312 else if ((temp_zone == 60) && (UTM_Override == 1))
00313 temp_zone = UTM_Override;
00314 else if (((temp_zone-1) <= UTM_Override) && (UTM_Override <= (temp_zone+1)))
00315 temp_zone = UTM_Override;
00316 else
00317 throw CoordinateConversionException( ErrorMessages::zoneOverride );
00318 }
00319 else
00320 {
00321
00322 if ((Lat_Degrees > 55) && (Lat_Degrees < 64) && (Long_Degrees > -1)
00323 && (Long_Degrees < 3))
00324 temp_zone = 31;
00325 if ((Lat_Degrees > 55) && (Lat_Degrees < 64) && (Long_Degrees > 2)
00326 && (Long_Degrees < 12))
00327 temp_zone = 32;
00328 if ((Lat_Degrees > 71) && (Long_Degrees > -1) && (Long_Degrees < 9))
00329 temp_zone = 31;
00330 if ((Lat_Degrees > 71) && (Long_Degrees > 8) && (Long_Degrees < 21))
00331 temp_zone = 33;
00332 if ((Lat_Degrees > 71) && (Long_Degrees > 20) && (Long_Degrees < 33))
00333 temp_zone = 35;
00334 if ((Lat_Degrees > 71) && (Long_Degrees > 32) && (Long_Degrees < 42))
00335 temp_zone = 37;
00336 }
00337
00338 TransverseMercator transverseMercator = *transverseMercatorMap[temp_zone];
00339
00340 if (latitude < 0)
00341 {
00342 False_Northing = 10000000;
00343 hemisphere = 'S';
00344 }
00345 else
00346 hemisphere = 'N';
00347
00348 GeodeticCoordinates tempGeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00349 MapProjectionCoordinates* transverseMercatorCoordinates = transverseMercator.convertFromGeodetic( &tempGeodeticCoordinates );
00350 double easting = transverseMercatorCoordinates->easting();
00351 double northing = transverseMercatorCoordinates->northing() + False_Northing;
00352
00353 if ((easting < MIN_EASTING) || (easting > MAX_EASTING))
00354 {
00355 delete transverseMercatorCoordinates;
00356 throw CoordinateConversionException( ErrorMessages::easting );
00357 }
00358
00359 if ((northing < MIN_NORTHING) || (northing > MAX_NORTHING))
00360 {
00361 delete transverseMercatorCoordinates;
00362 throw CoordinateConversionException( ErrorMessages::northing );
00363 }
00364
00365 delete transverseMercatorCoordinates;
00366
00367 return new UTMCoordinates( CoordinateType::universalTransverseMercator, temp_zone, hemisphere, easting, northing );
00368 }
00369
00370
00371 MSP::CCS::GeodeticCoordinates* UTM::convertToGeodetic( MSP::CCS::UTMCoordinates* utmCoordinates )
00372 {
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 double False_Northing = 0;
00389 char errorStatus[50] = "";
00390
00391 long zone = utmCoordinates->zone();
00392 char hemisphere = utmCoordinates->hemisphere();
00393 double easting = utmCoordinates->easting();
00394 double northing = utmCoordinates->northing();
00395
00396 if ((zone < 1) || (zone > 60))
00397 strcat( errorStatus, ErrorMessages::zone );
00398 if ((hemisphere != 'S') && (hemisphere != 'N'))
00399 strcat( errorStatus, ErrorMessages::hemisphere );
00400 if ((easting < MIN_EASTING) || (easting > MAX_EASTING))
00401 strcat( errorStatus, ErrorMessages::easting );
00402 if ((northing < MIN_NORTHING) || (northing > MAX_NORTHING))
00403 strcat( errorStatus, ErrorMessages::northing );
00404
00405 if( strlen( errorStatus ) > 0)
00406 throw CoordinateConversionException( errorStatus );
00407
00408 TransverseMercator transverseMercator = *transverseMercatorMap[zone];
00409
00410 if (hemisphere == 'S')
00411 False_Northing = 10000000;
00412
00413 MapProjectionCoordinates transverseMercatorCoordinates( CoordinateType::transverseMercator, easting, northing - False_Northing );
00414 GeodeticCoordinates* geodeticCoordinates = transverseMercator.convertToGeodetic( &transverseMercatorCoordinates );
00415 geodeticCoordinates->setWarningMessage("");
00416
00417 double latitude = geodeticCoordinates->latitude();
00418
00419 if ((latitude < (MIN_LAT - EPSILON)) || (latitude >= (MAX_LAT + EPSILON)))
00420 {
00421 delete geodeticCoordinates;
00422 throw CoordinateConversionException( ErrorMessages::northing );
00423 }
00424
00425 return geodeticCoordinates;
00426 }
00427
00428