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
00083 #include <math.h>
00084 #include "UTM.h"
00085 #include "TransverseMercator.h"
00086 #include "UTMParameters.h"
00087 #include "MapProjectionCoordinates.h"
00088 #include "UTMCoordinates.h"
00089 #include "GeodeticCoordinates.h"
00090 #include "CoordinateConversionException.h"
00091 #include "ErrorMessages.h"
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103 using namespace MSP::CCS;
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119 #define PI 3.14159265358979323e0
00120 #define PI_OVER_180 (3.14159265358979323e0 / 180.0)
00121 #define MIN_LAT ((-80.5 * PI) / 180.0)
00122 #define MAX_LAT (( 84.5 * PI) / 180.0)
00123 #define MIN_EASTING 100000.0
00124 #define MAX_EASTING 900000.0
00125 #define MIN_NORTHING 0.0
00126 #define MAX_NORTHING 10000000.0
00127
00128 #define EPSILON 1.75e-7
00129
00130
00131
00132
00133
00134
00135 UTM::UTM()
00136 {
00137
00138
00139
00140
00141 double ellipsoidSemiMajorAxis = 6378137.0;
00142 double ellipsoidFlattening = 1 / 298.257223563;
00143 double inv_f = 1 / ellipsoidFlattening;
00144
00145 semiMajorAxis = ellipsoidSemiMajorAxis;
00146 flattening = ellipsoidFlattening;
00147 UTM_Override = 0;
00148
00149 double centralMeridian;
00150 double originLatitude = 0;
00151 double falseEasting = 500000;
00152 double falseNorthing = 0;
00153 double scale = 0.9996;
00154
00155 for(int zone = 1; zone <= 60; zone++)
00156 {
00157 if (zone >= 31)
00158 centralMeridian = ((6 * zone - 183) * PI_OVER_180);
00159 else
00160 centralMeridian = ((6 * zone + 177) * PI_OVER_180);
00161
00162 transverseMercatorMap[zone] = new TransverseMercator(
00163 semiMajorAxis, flattening, centralMeridian, originLatitude,
00164 falseEasting, falseNorthing, scale);
00165 }
00166 }
00167
00168 UTM::UTM(
00169 double ellipsoidSemiMajorAxis,
00170 double ellipsoidFlattening,
00171 long override ) :
00172 UTM_Override( 0 )
00173 {
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186 double inv_f = 1 / ellipsoidFlattening;
00187 char errorStatus[500] = "";
00188
00189 if (ellipsoidSemiMajorAxis <= 0.0)
00190 {
00191 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00192 }
00193 if ((inv_f < 250) || (inv_f > 350))
00194 {
00195 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00196 }
00197 if ((override < 0) || (override > 60))
00198 {
00199 strcat( errorStatus, ErrorMessages::zoneOverride );
00200 }
00201
00202 if( strlen( errorStatus ) > 0)
00203 throw CoordinateConversionException( errorStatus );
00204
00205 semiMajorAxis = ellipsoidSemiMajorAxis;
00206 flattening = ellipsoidFlattening;
00207
00208 UTM_Override = override;
00209
00210 double centralMeridian;
00211 double originLatitude = 0;
00212 double falseEasting = 500000;
00213 double falseNorthing = 0;
00214 double scale = 0.9996;
00215
00216 for(int zone = 1; zone <= 60; zone++)
00217 {
00218 if (zone >= 31)
00219 centralMeridian = ((6 * zone - 183) * PI_OVER_180);
00220 else
00221 centralMeridian = ((6 * zone + 177) * PI_OVER_180);
00222
00223 transverseMercatorMap[zone] = new TransverseMercator(
00224 semiMajorAxis, flattening, centralMeridian, originLatitude,
00225 falseEasting, falseNorthing, scale);
00226 }
00227 }
00228
00229
00230 UTM::UTM( const UTM &u )
00231 {
00232 int zone = 1;
00233 std::map< int, TransverseMercator* > tempTransverseMercatorMap =
00234 u.transverseMercatorMap;
00235 std::map< int, TransverseMercator* >::iterator _iterator =
00236 tempTransverseMercatorMap.begin();
00237 while( _iterator != tempTransverseMercatorMap.end() && zone <= 60 )
00238 {
00239 transverseMercatorMap[zone] = new TransverseMercator( *_iterator->second );
00240 zone++;
00241 }
00242
00243 semiMajorAxis = u.semiMajorAxis;
00244 flattening = u.flattening;
00245 UTM_Override = u.UTM_Override;
00246 }
00247
00248
00249 UTM::~UTM()
00250 {
00251 while( transverseMercatorMap.begin() != transverseMercatorMap.end() )
00252 {
00253 delete ( ( *transverseMercatorMap.begin() ).second );
00254 transverseMercatorMap.erase( transverseMercatorMap.begin() );
00255 }
00256 }
00257
00258
00259 UTM& UTM::operator=( const UTM &u )
00260 {
00261 if( this != &u )
00262 {
00263 int zone = 1;
00264 std::map< int, TransverseMercator* > tempTransverseMercatorMap =
00265 u.transverseMercatorMap;
00266 std::map< int, TransverseMercator* >::iterator _iterator =
00267 tempTransverseMercatorMap.begin();
00268 while( _iterator != tempTransverseMercatorMap.end() && zone <= 60 )
00269 {
00270 transverseMercatorMap[zone]->operator=( *_iterator->second );
00271 zone++;
00272 }
00273
00275 semiMajorAxis = u.semiMajorAxis;
00276 flattening = u.flattening;
00277 UTM_Override = u.UTM_Override;
00278 }
00279
00280 return *this;
00281 }
00282
00283
00284 UTMParameters* UTM::getParameters() const
00285 {
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295 return new UTMParameters(
00296 CoordinateType::universalTransverseMercator, UTM_Override );
00297 }
00298
00299
00300 MSP::CCS::UTMCoordinates* UTM::convertFromGeodetic(
00301 MSP::CCS::GeodeticCoordinates* geodeticCoordinates,
00302 int utmZoneOverride )
00303 {
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319 long Lat_Degrees;
00320 long Long_Degrees;
00321 long temp_zone;
00322 char hemisphere;
00323 double False_Northing = 0;
00324 char errorStatus[50] = "";
00325
00326 double longitude = geodeticCoordinates->longitude();
00327 double latitude = geodeticCoordinates->latitude();
00328
00329 if ((latitude < (MIN_LAT - EPSILON)) || (latitude >= (MAX_LAT + EPSILON)))
00330 {
00331 strcat( errorStatus, ErrorMessages::latitude );
00332 }
00333 if ((longitude < -PI) || (longitude > (2*PI)))
00334 {
00335 strcat( errorStatus, ErrorMessages::longitude );
00336 }
00337
00338 if( strlen( errorStatus ) > 0)
00339 throw CoordinateConversionException( errorStatus );
00340
00341 if((latitude > -1.0e-9) && (latitude < 0))
00342 latitude = 0.0;
00343
00344 if (longitude < 0)
00345 longitude += (2*PI) + 1.0e-10;
00346
00347 Lat_Degrees = (long)(latitude * 180.0 / PI);
00348 Long_Degrees = (long)(longitude * 180.0 / PI);
00349
00350 if (longitude < PI)
00351 temp_zone = (long)(31 + ((longitude * 180.0 / PI) / 6.0));
00352 else
00353 temp_zone = (long)(((longitude * 180.0 / PI) / 6.0) - 29);
00354
00355 if (temp_zone > 60)
00356 temp_zone = 1;
00357
00358
00359 if( utmZoneOverride )
00360 {
00361 if ((temp_zone == 1) && (utmZoneOverride == 60))
00362 temp_zone = utmZoneOverride;
00363 else if ((temp_zone == 60) && (utmZoneOverride == 1))
00364 temp_zone = utmZoneOverride;
00365 else if (((temp_zone-1) <= utmZoneOverride) &&
00366 (utmZoneOverride <= (temp_zone+1)))
00367 temp_zone = utmZoneOverride;
00368 else
00369 throw CoordinateConversionException( ErrorMessages::zoneOverride );
00370 }
00371 else if( UTM_Override )
00372 {
00373 if ((temp_zone == 1) && (UTM_Override == 60))
00374 temp_zone = UTM_Override;
00375 else if ((temp_zone == 60) && (UTM_Override == 1))
00376 temp_zone = UTM_Override;
00377 else if (((temp_zone-1) <= UTM_Override) &&
00378 (UTM_Override <= (temp_zone+1)))
00379 temp_zone = UTM_Override;
00380 else
00381 throw CoordinateConversionException( ErrorMessages::zoneOverride );
00382 }
00383 else
00384 {
00385
00386 if ((Lat_Degrees > 55) && (Lat_Degrees < 64) && (Long_Degrees > -1)
00387 && (Long_Degrees < 3))
00388 temp_zone = 31;
00389 if ((Lat_Degrees > 55) && (Lat_Degrees < 64) && (Long_Degrees > 2)
00390 && (Long_Degrees < 12))
00391 temp_zone = 32;
00392 if ((Lat_Degrees > 71) && (Long_Degrees > -1) && (Long_Degrees < 9))
00393 temp_zone = 31;
00394 if ((Lat_Degrees > 71) && (Long_Degrees > 8) && (Long_Degrees < 21))
00395 temp_zone = 33;
00396 if ((Lat_Degrees > 71) && (Long_Degrees > 20) && (Long_Degrees < 33))
00397 temp_zone = 35;
00398 if ((Lat_Degrees > 71) && (Long_Degrees > 32) && (Long_Degrees < 42))
00399 temp_zone = 37;
00400 }
00401
00402 TransverseMercator *transverseMercator = transverseMercatorMap[temp_zone];
00403
00404 if (latitude < 0)
00405 {
00406 False_Northing = 10000000;
00407 hemisphere = 'S';
00408 }
00409 else
00410 hemisphere = 'N';
00411
00412 GeodeticCoordinates tempGeodeticCoordinates(
00413 CoordinateType::geodetic, longitude, latitude );
00414 MapProjectionCoordinates* transverseMercatorCoordinates =
00415 transverseMercator->convertFromGeodetic( &tempGeodeticCoordinates );
00416 double easting = transverseMercatorCoordinates->easting();
00417 double northing = transverseMercatorCoordinates->northing() + False_Northing;
00418
00419 if ((easting < MIN_EASTING) || (easting > MAX_EASTING))
00420 {
00421 delete transverseMercatorCoordinates;
00422 throw CoordinateConversionException( ErrorMessages::easting );
00423 }
00424
00425 if ((northing < MIN_NORTHING) || (northing > MAX_NORTHING))
00426 {
00427 delete transverseMercatorCoordinates;
00428 throw CoordinateConversionException( ErrorMessages::northing );
00429 }
00430
00431 delete transverseMercatorCoordinates;
00432
00433 return new UTMCoordinates(
00434 CoordinateType::universalTransverseMercator,
00435 temp_zone, hemisphere, easting, northing );
00436 }
00437
00438
00439 MSP::CCS::GeodeticCoordinates* UTM::convertToGeodetic(
00440 MSP::CCS::UTMCoordinates* utmCoordinates )
00441 {
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457 double False_Northing = 0;
00458 char errorStatus[50] = "";
00459
00460 long zone = utmCoordinates->zone();
00461 char hemisphere = utmCoordinates->hemisphere();
00462 double easting = utmCoordinates->easting();
00463 double northing = utmCoordinates->northing();
00464
00465 if ((zone < 1) || (zone > 60))
00466 strcat( errorStatus, ErrorMessages::zone );
00467 if ((hemisphere != 'S') && (hemisphere != 'N'))
00468 strcat( errorStatus, ErrorMessages::hemisphere );
00469 if ((easting < MIN_EASTING) || (easting > MAX_EASTING))
00470 strcat( errorStatus, ErrorMessages::easting );
00471 if ((northing < MIN_NORTHING) || (northing > MAX_NORTHING))
00472 strcat( errorStatus, ErrorMessages::northing );
00473
00474 if( strlen( errorStatus ) > 0)
00475 throw CoordinateConversionException( errorStatus );
00476
00477 TransverseMercator *transverseMercator = transverseMercatorMap[zone];
00478
00479 if (hemisphere == 'S')
00480 False_Northing = 10000000;
00481
00482 MapProjectionCoordinates transverseMercatorCoordinates(
00483 CoordinateType::transverseMercator, easting, northing - False_Northing );
00484 GeodeticCoordinates* geodeticCoordinates =
00485 transverseMercator->convertToGeodetic( &transverseMercatorCoordinates );
00486 geodeticCoordinates->setWarningMessage("");
00487
00488 double latitude = geodeticCoordinates->latitude();
00489
00490 if ((latitude < (MIN_LAT - EPSILON)) || (latitude >= (MAX_LAT + EPSILON)))
00491 {
00492 delete geodeticCoordinates;
00493 throw CoordinateConversionException( ErrorMessages::northing );
00494 }
00495
00496 return geodeticCoordinates;
00497 }
00498
00499