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
00084
00085
00086
00087
00088
00089 #include <math.h>
00090 #include "TransverseMercator.h"
00091 #include "MapProjection5Parameters.h"
00092 #include "MapProjectionCoordinates.h"
00093 #include "GeodeticCoordinates.h"
00094 #include "CoordinateConversionException.h"
00095 #include "ErrorMessages.h"
00096 #include "WarningMessages.h"
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108 using namespace MSP::CCS;
00109
00110
00111
00112
00113
00114
00115
00116 const double PI = 3.14159265358979323e0;
00117 const double PI_OVER_2 = (PI/2.0e0);
00118 const double MAX_LAT = ((PI * 89.99)/180.0);
00119 const double MAX_DELTA_LONG = ((PI * 90)/180.0);
00120 const double MIN_SCALE_FACTOR = 0.3;
00121 const double MAX_SCALE_FACTOR = 3.0;
00122
00123
00124
00125
00126
00127
00128
00129 TransverseMercator::TransverseMercator( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double latitudeOfTrueScale, double falseEasting, double falseNorthing, double scaleFactor ) :
00130 CoordinateSystem(),
00131 TranMerc_es( 0.0066943799901413800 ),
00132 TranMerc_ebs( 0.0067394967565869 ),
00133 TranMerc_Origin_Long( 0.0 ),
00134 TranMerc_Origin_Lat( 0.0 ),
00135 TranMerc_False_Easting( 0.0 ),
00136 TranMerc_False_Northing( 0.0 ),
00137 TranMerc_Scale_Factor( 1.0 ),
00138 TranMerc_ap( 6367449.1458008 ),
00139 TranMerc_bp( 16038.508696861 ),
00140 TranMerc_cp( 16.832613334334 ),
00141 TranMerc_dp( 0.021984404273757 ),
00142 TranMerc_ep( 3.1148371319283e-005 ),
00143 TranMerc_Delta_Easting( 40000000.0 ),
00144 TranMerc_Delta_Northing( 40000000.0 )
00145 {
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163 double tn;
00164 double tn2;
00165 double tn3;
00166 double tn4;
00167 double tn5;
00168 double TranMerc_b;
00169 double inv_f = 1 / ellipsoidFlattening;
00170 char errorStatus[500] = "";
00171
00172 if (ellipsoidSemiMajorAxis <= 0.0)
00173 {
00174 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00175 }
00176 if ((inv_f < 250) || (inv_f > 350))
00177 {
00178 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00179 }
00180 if ((latitudeOfTrueScale < -PI_OVER_2) || (latitudeOfTrueScale > PI_OVER_2))
00181 {
00182 strcat( errorStatus, ErrorMessages::originLatitude );
00183 }
00184 if ((centralMeridian < -PI) || (centralMeridian > (2*PI)))
00185 {
00186 strcat( errorStatus, ErrorMessages::centralMeridian );
00187 }
00188 if ((scaleFactor < MIN_SCALE_FACTOR) || (scaleFactor > MAX_SCALE_FACTOR))
00189 {
00190 strcat( errorStatus, ErrorMessages::scaleFactor );
00191 }
00192
00193 if( strlen( errorStatus ) > 0)
00194 throw CoordinateConversionException( errorStatus );
00195
00196 semiMajorAxis = ellipsoidSemiMajorAxis;
00197 flattening = ellipsoidFlattening;
00198
00199 TranMerc_Origin_Lat = latitudeOfTrueScale;
00200 if (centralMeridian > PI)
00201 centralMeridian -= (2*PI);
00202 TranMerc_Origin_Long = centralMeridian;
00203 TranMerc_False_Northing = falseNorthing;
00204 TranMerc_False_Easting = falseEasting;
00205 TranMerc_Scale_Factor = scaleFactor;
00206
00207
00208 TranMerc_es = 2 * flattening - flattening * flattening;
00209
00210 TranMerc_ebs = (1 / (1 - TranMerc_es)) - 1;
00211
00212 TranMerc_b = semiMajorAxis * (1 - flattening);
00213
00214 tn = (semiMajorAxis - TranMerc_b) / (semiMajorAxis + TranMerc_b);
00215 tn2 = tn * tn;
00216 tn3 = tn2 * tn;
00217 tn4 = tn3 * tn;
00218 tn5 = tn4 * tn;
00219
00220 TranMerc_ap = semiMajorAxis * (1.e0 - tn + 5.e0 * (tn2 - tn3)/4.e0
00221 + 81.e0 * (tn4 - tn5)/64.e0 );
00222 TranMerc_bp = 3.e0 * semiMajorAxis * (tn - tn2 + 7.e0 * (tn3 - tn4)
00223 /8.e0 + 55.e0 * tn5/64.e0 )/2.e0;
00224 TranMerc_cp = 15.e0 * semiMajorAxis * (tn2 - tn3 + 3.e0 * (tn4 - tn5 )/4.e0) /16.0;
00225 TranMerc_dp = 35.e0 * semiMajorAxis * (tn3 - tn4 + 11.e0 * tn5 / 16.e0) / 48.e0;
00226 TranMerc_ep = 315.e0 * semiMajorAxis * (tn4 - tn5) / 512.e0;
00227
00228 MapProjectionCoordinates* deltaCoordinates = 0;
00229 try
00230 {
00231 GeodeticCoordinates gcDeltaNorthing( CoordinateType::geodetic, MAX_DELTA_LONG + TranMerc_Origin_Long, MAX_LAT );
00232 deltaCoordinates = convertFromGeodetic( &gcDeltaNorthing );
00233 TranMerc_Delta_Northing = deltaCoordinates->northing();
00234 TranMerc_Delta_Northing ++;
00235
00236 delete deltaCoordinates;
00237
00238 GeodeticCoordinates gcDeltaEasting( CoordinateType::geodetic, MAX_DELTA_LONG + TranMerc_Origin_Long, 0 );
00239 deltaCoordinates = convertFromGeodetic( &gcDeltaEasting );
00240 TranMerc_Delta_Easting = deltaCoordinates->easting();
00241 TranMerc_Delta_Easting++;
00242
00243 delete deltaCoordinates;
00244 }
00245 catch(CoordinateConversionException e)
00246 {
00247 if( deltaCoordinates )
00248 delete deltaCoordinates;
00249
00250 TranMerc_Delta_Easting = 40000000.0;
00251 TranMerc_Delta_Northing = 40000000.0;
00252 }
00253 }
00254
00255
00256 TransverseMercator::TransverseMercator( const TransverseMercator &tm )
00257 {
00258 semiMajorAxis = tm.semiMajorAxis;
00259 flattening = tm.flattening;
00260 TranMerc_es = tm.TranMerc_es;
00261 TranMerc_ebs = tm.TranMerc_ebs;
00262 TranMerc_Origin_Long = tm.TranMerc_Origin_Long;
00263 TranMerc_Origin_Lat = tm.TranMerc_Origin_Lat;
00264 TranMerc_False_Easting = tm.TranMerc_False_Easting;
00265 TranMerc_False_Northing = tm.TranMerc_False_Northing;
00266 TranMerc_Scale_Factor = tm.TranMerc_Scale_Factor;
00267 TranMerc_ap = tm.TranMerc_ap;
00268 TranMerc_bp = tm.TranMerc_bp;
00269 TranMerc_cp = tm.TranMerc_cp;
00270 TranMerc_dp = tm.TranMerc_dp;
00271 TranMerc_ep = tm.TranMerc_ep;
00272 TranMerc_Delta_Easting = tm.TranMerc_Delta_Easting;
00273 TranMerc_Delta_Northing = tm.TranMerc_Delta_Northing;
00274 }
00275
00276
00277 TransverseMercator::~TransverseMercator()
00278 {
00279 }
00280
00281
00282 TransverseMercator& TransverseMercator::operator=( const TransverseMercator &tm )
00283 {
00284 if( this != &tm )
00285 {
00286 semiMajorAxis = tm.semiMajorAxis;
00287 flattening = tm.flattening;
00288 TranMerc_es = tm.TranMerc_es;
00289 TranMerc_ebs = tm.TranMerc_ebs;
00290 TranMerc_Origin_Long = tm.TranMerc_Origin_Long;
00291 TranMerc_Origin_Lat = tm.TranMerc_Origin_Lat;
00292 TranMerc_False_Easting = tm.TranMerc_False_Easting;
00293 TranMerc_False_Northing = tm.TranMerc_False_Northing;
00294 TranMerc_Scale_Factor = tm.TranMerc_Scale_Factor;
00295 TranMerc_ap = tm.TranMerc_ap;
00296 TranMerc_bp = tm.TranMerc_bp;
00297 TranMerc_cp = tm.TranMerc_cp;
00298 TranMerc_dp = tm.TranMerc_dp;
00299 TranMerc_ep = tm.TranMerc_ep;
00300 TranMerc_Delta_Easting = tm.TranMerc_Delta_Easting;
00301 TranMerc_Delta_Northing = tm.TranMerc_Delta_Northing;
00302 }
00303
00304 return *this;
00305 }
00306
00307
00308 MapProjection5Parameters* TransverseMercator::getParameters() const
00309 {
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325 return new MapProjection5Parameters( CoordinateType::transverseMercator, TranMerc_Origin_Long, TranMerc_Origin_Lat, TranMerc_Scale_Factor, TranMerc_False_Easting, TranMerc_False_Northing );
00326 }
00327
00328
00329 MSP::CCS::MapProjectionCoordinates* TransverseMercator::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00330 {
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344 double c;
00345 double c2;
00346 double c3;
00347 double c5;
00348 double c7;
00349 double dlam;
00350 double eta;
00351 double eta2;
00352 double eta3;
00353 double eta4;
00354 double s;
00355 double sn;
00356 double t;
00357 double tan2;
00358 double tan3;
00359 double tan4;
00360 double tan5;
00361 double tan6;
00362 double t1;
00363 double t2;
00364 double t3;
00365 double t4;
00366 double t5;
00367 double t6;
00368 double t7;
00369 double t8;
00370 double t9;
00371 double tmd;
00372 double tmdo;
00373 double temp_Origin;
00374 double temp_Long;
00375 char errorStatus[256] = "";
00376
00377 double longitude = geodeticCoordinates->longitude();
00378 double latitude = geodeticCoordinates->latitude();
00379
00380 if ((latitude < -MAX_LAT) || (latitude > MAX_LAT))
00381 {
00382 strcat( errorStatus, ErrorMessages::latitude );
00383 }
00384 if (longitude > PI)
00385 longitude -= (2 * PI);
00386 if ((longitude < (TranMerc_Origin_Long - MAX_DELTA_LONG))
00387 || (longitude > (TranMerc_Origin_Long + MAX_DELTA_LONG)))
00388 {
00389 if (longitude < 0)
00390 temp_Long = longitude + 2 * PI;
00391 else
00392 temp_Long = longitude;
00393 if (TranMerc_Origin_Long < 0)
00394 temp_Origin = TranMerc_Origin_Long + 2 * PI;
00395 else
00396 temp_Origin = TranMerc_Origin_Long;
00397 if ((temp_Long < (temp_Origin - MAX_DELTA_LONG))
00398 || (temp_Long > (temp_Origin + MAX_DELTA_LONG)))
00399 strcat( errorStatus, ErrorMessages::longitude );
00400 }
00401
00402 if( strlen( errorStatus ) > 0)
00403 throw CoordinateConversionException( errorStatus );
00404
00405
00406
00407
00408 dlam = longitude - TranMerc_Origin_Long;
00409
00410 if (fabs(dlam) > (9.0 * PI / 180))
00411 {
00412 strcat( errorStatus, MSP::CCS::WarningMessages::longitude );
00413 }
00414
00415 if (dlam > PI)
00416 dlam -= (2 * PI);
00417 if (dlam < -PI)
00418 dlam += (2 * PI);
00419 if (fabs(dlam) < 2.e-10)
00420 dlam = 0.0;
00421
00422 s = sin(latitude);
00423 c = cos(latitude);
00424 c2 = c * c;
00425 c3 = c2 * c;
00426 c5 = c3 * c2;
00427 c7 = c5 * c2;
00428 t = tan (latitude);
00429 tan2 = t * t;
00430 tan3 = tan2 * t;
00431 tan4 = tan3 * t;
00432 tan5 = tan4 * t;
00433 tan6 = tan5 * t;
00434 eta = TranMerc_ebs * c2;
00435 eta2 = eta * eta;
00436 eta3 = eta2 * eta;
00437 eta4 = eta3 * eta;
00438
00439
00440 sn = sphsn(latitude);
00441
00442
00443 tmd = sphtmd(latitude);
00444
00445
00446 tmdo = sphtmd (TranMerc_Origin_Lat);
00447
00448
00449 t1 = (tmd - tmdo) * TranMerc_Scale_Factor;
00450 t2 = sn * s * c * TranMerc_Scale_Factor/ 2.e0;
00451 t3 = sn * s * c3 * TranMerc_Scale_Factor * (5.e0 - tan2 + 9.e0 * eta
00452 + 4.e0 * eta2) /24.e0;
00453
00454 t4 = sn * s * c5 * TranMerc_Scale_Factor * (61.e0 - 58.e0 * tan2
00455 + tan4 + 270.e0 * eta - 330.e0 * tan2 * eta + 445.e0 * eta2
00456 + 324.e0 * eta3 -680.e0 * tan2 * eta2 + 88.e0 * eta4
00457 -600.e0 * tan2 * eta3 - 192.e0 * tan2 * eta4) / 720.e0;
00458
00459 t5 = sn * s * c7 * TranMerc_Scale_Factor * (1385.e0 - 3111.e0 *
00460 tan2 + 543.e0 * tan4 - tan6) / 40320.e0;
00461
00462 double northing = TranMerc_False_Northing + t1 + pow(dlam, 2.e0) * t2
00463 + pow(dlam,4.e0) * t3 + pow(dlam,6.e0) * t4
00464 + pow(dlam,8.e0) * t5;
00465
00466
00467 t6 = sn * c * TranMerc_Scale_Factor;
00468 t7 = sn * c3 * TranMerc_Scale_Factor * (1.e0 - tan2 + eta ) /6.e0;
00469 t8 = sn * c5 * TranMerc_Scale_Factor * (5.e0 - 18.e0 * tan2 + tan4
00470 + 14.e0 * eta - 58.e0 * tan2 * eta + 13.e0 * eta2 + 4.e0 * eta3
00471 - 64.e0 * tan2 * eta2 - 24.e0 * tan2 * eta3 )/ 120.e0;
00472 t9 = sn * c7 * TranMerc_Scale_Factor * ( 61.e0 - 479.e0 * tan2
00473 + 179.e0 * tan4 - tan6 ) /5040.e0;
00474
00475 double easting = TranMerc_False_Easting + dlam * t6 + pow(dlam,3.e0) * t7
00476 + pow(dlam,5.e0) * t8 + pow(dlam,7.e0) * t9;
00477
00478 return new MapProjectionCoordinates( CoordinateType::transverseMercator, errorStatus, easting, northing );
00479 }
00480
00481
00482 MSP::CCS::GeodeticCoordinates* TransverseMercator::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00483 {
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497 double c;
00498 double de;
00499 double dlam;
00500 double eta;
00501 double eta2;
00502 double eta3;
00503 double eta4;
00504 double ftphi;
00505 int i;
00506 double s;
00507 double sn;
00508 double sr;
00509 double t;
00510 double tan2;
00511 double tan4;
00512 double t10;
00513 double t11;
00514 double t12;
00515 double t13;
00516 double t14;
00517 double t15;
00518 double t16;
00519 double t17;
00520 double tmd;
00521 double tmdo;
00522 char errorStatus[256] = "";
00523
00524 double easting = mapProjectionCoordinates->easting();
00525 double northing = mapProjectionCoordinates->northing();
00526
00527 if ((easting < (TranMerc_False_Easting - TranMerc_Delta_Easting))
00528 ||(easting > (TranMerc_False_Easting + TranMerc_Delta_Easting)))
00529 {
00530 strcat( errorStatus, ErrorMessages::easting );
00531 }
00532 if ((northing < (TranMerc_False_Northing - TranMerc_Delta_Northing))
00533 || (northing > (TranMerc_False_Northing + TranMerc_Delta_Northing)))
00534 {
00535 strcat( errorStatus, ErrorMessages::northing );
00536 }
00537
00538 if( strlen( errorStatus ) > 0)
00539 throw CoordinateConversionException( errorStatus );
00540
00541
00542 tmdo = sphtmd(TranMerc_Origin_Lat);
00543
00544
00545 tmd = tmdo + (northing - TranMerc_False_Northing) / TranMerc_Scale_Factor;
00546
00547
00548 sr = sphsr(0.e0);
00549 ftphi = tmd/sr;
00550
00551 for (i = 0; i < 5 ; i++)
00552 {
00553 t10 = sphtmd (ftphi);
00554 sr = sphsr(ftphi);
00555 ftphi = ftphi + (tmd - t10) / sr;
00556 }
00557
00558
00559 sr = sphsr(ftphi);
00560
00561
00562 sn = sphsn(ftphi);
00563
00564
00565 s = sin(ftphi);
00566 c = cos(ftphi);
00567
00568
00569 t = tan(ftphi);
00570 tan2 = t * t;
00571 tan4 = tan2 * tan2;
00572 eta = TranMerc_ebs * pow(c, 2.0);
00573 eta2 = eta * eta;
00574 eta3 = eta2 * eta;
00575 eta4 = eta3 * eta;
00576 de = easting - TranMerc_False_Easting;
00577 if (fabs(de) < 0.0001)
00578 de = 0.0;
00579
00580
00581 t10 = t / (2.e0 * sr * sn * pow(TranMerc_Scale_Factor, 2.0));
00582 t11 = t * (5.e0 + 3.e0 * tan2 + eta - 4.e0 * pow(eta, 2.0)
00583 - 9.e0 * tan2 * eta) / (24.e0 * sr * pow(sn, 3.0)
00584 * pow(TranMerc_Scale_Factor, 4.0));
00585 t12 = t * (61.e0 + 90.e0 * tan2 + 46.e0 * eta + 45.E0 * tan4
00586 - 252.e0 * tan2 * eta - 3.e0 * eta2 + 100.e0
00587 * eta3 - 66.e0 * tan2 * eta2 - 90.e0 * tan4
00588 * eta + 88.e0 * eta4 + 225.e0 * tan4 * eta2
00589 + 84.e0 * tan2* eta3 - 192.e0 * tan2 * eta4)
00590 / ( 720.e0 * sr * pow(sn,5.0) * pow(TranMerc_Scale_Factor, 6.0) );
00591 t13 = t * ( 1385.e0 + 3633.e0 * tan2 + 4095.e0 * tan4 + 1575.e0
00592 * pow(t, 6.0))/ (40320.e0 * sr * pow(sn, 7.0) * pow(TranMerc_Scale_Factor, 8.0));
00593 double latitude = ftphi - pow(de, 2.0) * t10 + pow(de, 4.0) * t11 - pow(de, 6.0) * t12
00594 + pow(de, 8.0) * t13;
00595
00596 t14 = 1.e0 / (sn * c * TranMerc_Scale_Factor);
00597
00598 t15 = (1.e0 + 2.e0 * tan2 + eta) / (6.e0 * pow(sn, 3.0) * c *
00599 pow(TranMerc_Scale_Factor, 3.0));
00600
00601 t16 = (5.e0 + 6.e0 * eta + 28.e0 * tan2 - 3.e0 * eta2
00602 + 8.e0 * tan2 * eta + 24.e0 * tan4 - 4.e0
00603 * eta3 + 4.e0 * tan2 * eta2 + 24.e0
00604 * tan2 * eta3) / (120.e0 * pow(sn, 5.0) * c
00605 * pow(TranMerc_Scale_Factor, 5.0));
00606
00607 t17 = (61.e0 + 662.e0 * tan2 + 1320.e0 * tan4 + 720.e0
00608 * pow(t, 6.0)) / (5040.e0 * pow(sn, 7.0) * c
00609 * pow(TranMerc_Scale_Factor, 7.0));
00610
00611
00612 dlam = de * t14 - pow(de, 3.0) * t15 + pow(de, 5.0) * t16 - pow(de, 7.0) * t17;
00613
00614
00615 double longitude = TranMerc_Origin_Long + dlam;
00616
00617 if(fabs(latitude) > (90.0 * PI / 180.0))
00618 throw CoordinateConversionException( ErrorMessages::northing );
00619
00620 if((longitude) > (PI))
00621 {
00622 longitude -= (2 * PI);
00623 if(fabs(longitude) > PI)
00624 throw CoordinateConversionException( ErrorMessages::easting );
00625 }
00626 else if((longitude) < (-PI))
00627 {
00628 longitude += (2 * PI);
00629 if(fabs(longitude) > PI)
00630 throw CoordinateConversionException( ErrorMessages::easting );
00631 }
00632
00633 if (fabs(dlam) > (9.0 * PI / 180))
00634 {
00635
00636
00637 strcat( errorStatus, MSP::CCS::WarningMessages::longitude );
00638 }
00639
00640 return new GeodeticCoordinates( CoordinateType::geodetic, errorStatus, longitude, latitude );
00641 }
00642
00643
00644
00645
00646
00647
00648
00649 double TransverseMercator::sphtmd( double latitude )
00650 {
00651 return TranMerc_ap * latitude
00652 - TranMerc_bp * sin(2.e0 * latitude) + TranMerc_cp * sin(4.e0 * latitude)
00653 - TranMerc_dp * sin(6.e0 * latitude) + TranMerc_ep * sin(8.e0 * latitude);
00654 }
00655
00656
00657 double TransverseMercator::sphsn( double latitude )
00658 {
00659 return semiMajorAxis / sqrt( 1.e0 - TranMerc_es * pow(sin(latitude), 2.0));
00660 }
00661
00662
00663 double TransverseMercator::sphsr( double latitude )
00664 {
00665 double denom = sqrt(1.e0 - TranMerc_es * pow(sin(latitude), 2.0));
00666 return semiMajorAxis * (1.e0 - TranMerc_es) / pow(denom, 3.0);
00667 }
00668
00669
00670
00671