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 #include <math.h>
00089 #include "TransverseCylindricalEqualArea.h"
00090 #include "MapProjection5Parameters.h"
00091 #include "MapProjectionCoordinates.h"
00092 #include "GeodeticCoordinates.h"
00093 #include "CoordinateConversionException.h"
00094 #include "ErrorMessages.h"
00095 #include "WarningMessages.h"
00096
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.0);
00118 const double TWO_PI = (2.0 * PI);
00119 const double MIN_SCALE_FACTOR = 0.3;
00120 const double MAX_SCALE_FACTOR = 3.0;
00121
00122
00123
00124
00125
00126
00127
00128 TransverseCylindricalEqualArea::TransverseCylindricalEqualArea( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double latitudeOfTrueScale, double falseEasting, double falseNorthing, double scaleFactor ) :
00129 CoordinateSystem(),
00130 es2( 0.0066943799901413800 ),
00131 es4( 4.4814723452405e-005 ),
00132 es6( 3.0000678794350e-007 ),
00133 es( .081819190842622 ),
00134 M0( 0.0 ),
00135 qp( 1.9955310875028 ),
00136 One_MINUS_es2( .99330562000986 ),
00137 One_OVER_2es( 6.1110357466348 ),
00138 a0( .0022392088624809 ),
00139 a1( 2.8830839728915e-006 ),
00140 a2( 5.0331826636906e-009 ),
00141 b0( .0025188265843907 ),
00142 b1( 3.7009490356205e-006 ),
00143 b2( 7.4478137675038e-009 ),
00144 b3( 1.7035993238596e-011 ),
00145 c0( .99832429845280 ),
00146 c1( .0025146070605187 ),
00147 c2( 2.6390465943377e-006 ),
00148 c3( 3.4180460865959e-009 ),
00149 Tcea_Origin_Long( 0.0 ),
00150 Tcea_Origin_Lat( 0.0 ),
00151 Tcea_False_Easting( 0.0 ),
00152 Tcea_False_Northing( 0.0 ),
00153 Tcea_Scale_Factor( 1.0 ),
00154 Tcea_Min_Easting( -6398628.0 ),
00155 Tcea_Max_Easting( 6398628.0 ),
00156 Tcea_Min_Northing( -20003931.0 ),
00157 Tcea_Max_Northing( 20003931.0 )
00158 {
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180 double sin_lat_90 = sin(PI_OVER_2);
00181 double x, j, three_es4;
00182 double Sqrt_One_MINUS_es2;
00183 double e1, e2, e3, e4;
00184 double lat, sin2lat, sin4lat, sin6lat;
00185 double temp, temp_northing;
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 ((latitudeOfTrueScale < -PI_OVER_2) || (latitudeOfTrueScale > PI_OVER_2))
00198 {
00199 strcat( errorStatus, ErrorMessages::originLatitude );
00200 }
00201 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00202 {
00203 strcat( errorStatus, ErrorMessages::centralMeridian );
00204 }
00205 if ((scaleFactor < MIN_SCALE_FACTOR) || (scaleFactor > MAX_SCALE_FACTOR))
00206 {
00207 strcat( errorStatus, ErrorMessages::scaleFactor );
00208 }
00209
00210 if( strlen( errorStatus ) > 0)
00211 throw CoordinateConversionException( errorStatus );
00212
00213 semiMajorAxis = ellipsoidSemiMajorAxis;
00214 flattening = ellipsoidFlattening;
00215
00216 Tcea_Origin_Lat = latitudeOfTrueScale;
00217 if (centralMeridian > PI)
00218 centralMeridian -= TWO_PI;
00219 Tcea_Origin_Long = centralMeridian;
00220 Tcea_False_Northing = falseNorthing;
00221 Tcea_False_Easting = falseEasting;
00222 Tcea_Scale_Factor = scaleFactor;
00223
00224 es2 = 2 * flattening - flattening * flattening;
00225 es = sqrt(es2);
00226 One_MINUS_es2 = 1.0 - es2;
00227 Sqrt_One_MINUS_es2 = sqrt(One_MINUS_es2);
00228 One_OVER_2es = 1.0 / (2.0 * es);
00229 es4 = es2 * es2;
00230 es6 = es4 * es2;
00231 x = es * sin_lat_90;
00232 qp = TCEA_Q(sin_lat_90,x);
00233
00234 a0 = es2 / 3.0 + 31.0 * es4 / 180.0 + 517.0 * es6 / 5040.0;
00235 a1 = 23.0 * es4 / 360.0 + 251.0 * es6 / 3780.0;
00236 a2 = 761.0 * es6 / 45360.0;
00237
00238 e1 = (1.0 - Sqrt_One_MINUS_es2) / (1.0 + Sqrt_One_MINUS_es2);
00239 e2 = e1 * e1;
00240 e3 = e2 * e1;
00241 e4 = e3 * e1;
00242 b0 = 3.0 * e1 / 2.0 - 27.0 * e3 / 32.0;
00243 b1 = 21.0 * e2 / 16.0 - 55.0 * e4 / 32.0;
00244 b2 = 151.0 * e3 / 96.0;
00245 b3 = 1097.0 * e4 / 512.0;
00246
00247 j = 45.0 * es6 / 1024.0;
00248 three_es4 = 3.0 * es4;
00249 c0 = 1.0 - es2 / 4.0 - three_es4 / 64.0 - 5.0 * es6 / 256.0;
00250 c1 = 3.0 * es2 / 8.0 + three_es4 / 32.0 + j;
00251 c2 = 15.0 * es4 / 256.0 + j;
00252 c3 = 35.0 * es6 / 3072.0;
00253 lat = c0 * Tcea_Origin_Lat;
00254 sin2lat = TCEA_COEFF_TIMES_SIN(c1, 2.0, Tcea_Origin_Lat);
00255 sin4lat = TCEA_COEFF_TIMES_SIN(c2, 4.0, Tcea_Origin_Lat);
00256 sin6lat = TCEA_COEFF_TIMES_SIN(c3, 6.0, Tcea_Origin_Lat);
00257 M0 = TCEA_M(lat, sin2lat, sin4lat, sin6lat);
00258
00259 GeodeticCoordinates gcTemp( CoordinateType::geodetic, PI, PI_OVER_2 );
00260 MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcTemp );
00261 temp = tempCoordinates->easting();
00262 temp_northing = tempCoordinates->northing();
00263 delete tempCoordinates;
00264
00265 if (temp_northing > 0)
00266 {
00267 Tcea_Min_Northing = temp_northing - 20003931.458986;
00268 Tcea_Max_Northing = temp_northing;
00269
00270 if(Tcea_False_Northing)
00271 {
00272 Tcea_Min_Northing -= Tcea_False_Northing;
00273 Tcea_Max_Northing -= Tcea_False_Northing;
00274 }
00275 }
00276 else if (temp_northing < 0)
00277 {
00278 Tcea_Max_Northing = temp_northing + 20003931.458986;
00279 Tcea_Min_Northing = temp_northing;
00280
00281 if(Tcea_False_Northing)
00282 {
00283 Tcea_Min_Northing -= Tcea_False_Northing;
00284 Tcea_Max_Northing -= Tcea_False_Northing;
00285 }
00286 }
00287 }
00288
00289
00290 TransverseCylindricalEqualArea::TransverseCylindricalEqualArea( const TransverseCylindricalEqualArea &tcea )
00291 {
00292 semiMajorAxis = tcea.semiMajorAxis;
00293 flattening = tcea.flattening;
00294 es2 = tcea.es2;
00295 es4 = tcea.es4;
00296 es6 = tcea.es6;
00297 es = tcea.es;
00298 M0 = tcea.M0;
00299 qp = tcea.qp;
00300 One_MINUS_es2 = tcea.One_MINUS_es2;
00301 One_OVER_2es = tcea.One_OVER_2es;
00302 a0 = tcea.a0;
00303 a1 = tcea.a1;
00304 a2 = tcea.a2;
00305 b0 = tcea.b0;
00306 b1 = tcea.b1;
00307 b2 = tcea.b2;
00308 b3 = tcea.b3;
00309 c0 = tcea.c0;
00310 c1 = tcea.c1;
00311 c2 = tcea.c2;
00312 c3 = tcea.c3;
00313 Tcea_Origin_Long = tcea.Tcea_Origin_Long;
00314 Tcea_Origin_Lat = tcea.Tcea_Origin_Lat;
00315 Tcea_False_Easting = tcea.Tcea_False_Easting;
00316 Tcea_False_Northing = tcea.Tcea_False_Northing;
00317 Tcea_Scale_Factor = tcea.Tcea_Scale_Factor;
00318 Tcea_Min_Easting = tcea.Tcea_Min_Easting;
00319 Tcea_Max_Easting = tcea.Tcea_Max_Easting;
00320 Tcea_Min_Northing = tcea.Tcea_Min_Northing;
00321 Tcea_Max_Northing = tcea.Tcea_Max_Northing;
00322 }
00323
00324
00325 TransverseCylindricalEqualArea::~TransverseCylindricalEqualArea()
00326 {
00327 }
00328
00329
00330 TransverseCylindricalEqualArea& TransverseCylindricalEqualArea::operator=( const TransverseCylindricalEqualArea &tcea )
00331 {
00332 if( this != &tcea )
00333 {
00334 semiMajorAxis = tcea.semiMajorAxis;
00335 flattening = tcea.flattening;
00336 es2 = tcea.es2;
00337 es4 = tcea.es4;
00338 es6 = tcea.es6;
00339 es = tcea.es;
00340 M0 = tcea.M0;
00341 qp = tcea.qp;
00342 One_MINUS_es2 = tcea.One_MINUS_es2;
00343 One_OVER_2es = tcea.One_OVER_2es;
00344 a0 = tcea.a0;
00345 a1 = tcea.a1;
00346 a2 = tcea.a2;
00347 b0 = tcea.b0;
00348 b1 = tcea.b1;
00349 b2 = tcea.b2;
00350 b3 = tcea.b3;
00351 c0 = tcea.c0;
00352 c1 = tcea.c1;
00353 c2 = tcea.c2;
00354 c3 = tcea.c3;
00355 Tcea_Origin_Long = tcea.Tcea_Origin_Long;
00356 Tcea_Origin_Lat = tcea.Tcea_Origin_Lat;
00357 Tcea_False_Easting = tcea.Tcea_False_Easting;
00358 Tcea_False_Northing = tcea.Tcea_False_Northing;
00359 Tcea_Scale_Factor = tcea.Tcea_Scale_Factor;
00360 Tcea_Min_Easting = tcea.Tcea_Min_Easting;
00361 Tcea_Max_Easting = tcea.Tcea_Max_Easting;
00362 Tcea_Min_Northing = tcea.Tcea_Min_Northing;
00363 Tcea_Max_Northing = tcea.Tcea_Max_Northing;
00364 }
00365
00366 return *this;
00367 }
00368
00369
00370 MapProjection5Parameters* TransverseCylindricalEqualArea::getParameters() const
00371 {
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391 return new MapProjection5Parameters( CoordinateType::transverseCylindricalEqualArea, Tcea_Origin_Long, Tcea_Origin_Lat, Tcea_Scale_Factor, Tcea_False_Easting, Tcea_False_Northing );
00392 }
00393
00394
00395 MSP::CCS::MapProjectionCoordinates* TransverseCylindricalEqualArea::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00396 {
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410 double x;
00411 double dlam;
00412 double qq, qq_OVER_qp;
00413 double beta, betac;
00414 double sin2betac, sin4betac, sin6betac;
00415 double PHIc;
00416 double phi, sin2phi, sin4phi, sin6phi;
00417 double sinPHIc;
00418 double Mc;
00419 char errorStatus[256] = "";
00420
00421 double longitude = geodeticCoordinates->longitude();
00422 double latitude = geodeticCoordinates->latitude();
00423 double sin_lat = sin(latitude);
00424
00425 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00426 {
00427 strcat( errorStatus, ErrorMessages::latitude );
00428 }
00429 if ((longitude < -PI) || (longitude > TWO_PI))
00430 {
00431 strcat( errorStatus, ErrorMessages::longitude );
00432 }
00433
00434 if( strlen( errorStatus ) > 0)
00435 throw CoordinateConversionException( errorStatus );
00436
00437 dlam = longitude - Tcea_Origin_Long;
00438 if (fabs(dlam) >= (PI / 2))
00439 {
00440 strcat( errorStatus, MSP::CCS::WarningMessages::longitude );
00441 }
00442
00443 if (dlam > PI)
00444 {
00445 dlam -= TWO_PI;
00446 }
00447 if (dlam < -PI)
00448 {
00449 dlam += TWO_PI;
00450 }
00451 if (latitude == PI_OVER_2)
00452 {
00453 qq = qp;
00454 qq_OVER_qp = 1.0;
00455 }
00456 else
00457 {
00458 x = es * sin_lat;
00459 qq = TCEA_Q(sin_lat, x);
00460 qq_OVER_qp = qq / qp;
00461 }
00462
00463
00464 if (qq_OVER_qp > 1.0)
00465 qq_OVER_qp = 1.0;
00466 else if (qq_OVER_qp < -1.0)
00467 qq_OVER_qp = -1.0;
00468
00469 beta = asin(qq_OVER_qp);
00470 betac = atan(tan(beta) / cos(dlam));
00471
00472 if ((fabs(betac) - PI_OVER_2) > 1.0e-8)
00473 PHIc = betac;
00474 else
00475 {
00476 sin2betac = TCEA_COEFF_TIMES_SIN(a0, 2.0, betac);
00477 sin4betac = TCEA_COEFF_TIMES_SIN(a1, 4.0, betac);
00478 sin6betac = TCEA_COEFF_TIMES_SIN(a2, 6.0, betac);
00479 PHIc = TCEA_L(betac, sin2betac, sin4betac, sin6betac);
00480 }
00481
00482 sinPHIc = sin(PHIc);
00483 double easting = semiMajorAxis * cos(beta) * cos(PHIc) * sin(dlam) /
00484 (Tcea_Scale_Factor * cos(betac) * sqrt(1.0 - es2 *
00485 sinPHIc * sinPHIc)) + Tcea_False_Easting;
00486
00487 phi = c0 * PHIc;
00488 sin2phi = TCEA_COEFF_TIMES_SIN(c1, 2.0, PHIc);
00489 sin4phi = TCEA_COEFF_TIMES_SIN(c2, 4.0, PHIc);
00490 sin6phi = TCEA_COEFF_TIMES_SIN(c3, 6.0, PHIc);
00491 Mc = TCEA_M(phi, sin2phi, sin4phi, sin6phi);
00492
00493 double northing = Tcea_Scale_Factor * (Mc - M0) + Tcea_False_Northing;
00494
00495 return new MapProjectionCoordinates( CoordinateType::transverseCylindricalEqualArea, errorStatus, easting, northing );
00496 }
00497
00498
00499 MSP::CCS::GeodeticCoordinates* TransverseCylindricalEqualArea::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00500 {
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515 double x;
00516 double dx;
00517 double dy;
00518 double Mc;
00519 double MUc;
00520 double sin2mu, sin4mu, sin6mu, sin8mu;
00521 double PHIc;
00522 double Qc;
00523 double sin_lat;
00524 double beta, betac, beta_prime;
00525 double sin2beta, sin4beta, sin6beta;
00526 double cosbetac;
00527 double Qc_OVER_qp;
00528 double temp;
00529 char errorStatus[50] = "";
00530
00531 double easting = mapProjectionCoordinates->easting();
00532 double northing = mapProjectionCoordinates->northing();
00533
00534 if ((easting < (Tcea_False_Easting + Tcea_Min_Easting))
00535 || (easting > (Tcea_False_Easting + Tcea_Max_Easting)))
00536 {
00537 strcat( errorStatus, ErrorMessages::easting );
00538 }
00539 if ((northing < (Tcea_False_Northing + Tcea_Min_Northing))
00540 || (northing > (Tcea_False_Northing + Tcea_Max_Northing)))
00541 {
00542 strcat( errorStatus, ErrorMessages::northing );
00543 }
00544
00545 if( strlen( errorStatus ) > 0)
00546 throw CoordinateConversionException( errorStatus );
00547
00548 dy = northing - Tcea_False_Northing;
00549 dx = easting - Tcea_False_Easting;
00550 Mc = M0 + dy / Tcea_Scale_Factor;
00551 MUc = Mc / (semiMajorAxis * c0);
00552
00553 sin2mu = TCEA_COEFF_TIMES_SIN(b0, 2.0, MUc);
00554 sin4mu = TCEA_COEFF_TIMES_SIN(b1, 4.0, MUc);
00555 sin6mu = TCEA_COEFF_TIMES_SIN(b2, 6.0, MUc);
00556 sin8mu = TCEA_COEFF_TIMES_SIN(b3, 8.0, MUc);
00557 PHIc = MUc + sin2mu + sin4mu + sin6mu + sin8mu;
00558
00559 sin_lat = sin(PHIc);
00560 x = es * sin_lat;
00561 Qc = TCEA_Q(sin_lat, x);
00562 Qc_OVER_qp = Qc / qp;
00563
00564 if (Qc_OVER_qp < -1.0)
00565 Qc_OVER_qp = -1.0;
00566 else if (Qc_OVER_qp > 1.0)
00567 Qc_OVER_qp = 1.0;
00568
00569 betac = asin(Qc_OVER_qp);
00570 cosbetac = cos(betac);
00571 temp = Tcea_Scale_Factor * dx * cosbetac * sqrt(1.0 -
00572 es2 * sin_lat * sin_lat) / (semiMajorAxis * cos(PHIc));
00573 if (temp > 1.0)
00574 temp = 1.0;
00575 else if (temp < -1.0)
00576 temp = -1.0;
00577 beta_prime = -asin(temp);
00578 beta = asin(cos(beta_prime) * sin(betac));
00579
00580 sin2beta = TCEA_COEFF_TIMES_SIN(a0, 2.0, beta);
00581 sin4beta = TCEA_COEFF_TIMES_SIN(a1, 4.0, beta);
00582 sin6beta = TCEA_COEFF_TIMES_SIN(a2, 6.0, beta);
00583 double latitude = TCEA_L(beta, sin2beta, sin4beta, sin6beta);
00584
00585 double longitude = Tcea_Origin_Long - atan(tan(beta_prime) / cosbetac);
00586
00587 if (latitude > PI_OVER_2)
00588 latitude = PI_OVER_2;
00589 else if (latitude < -PI_OVER_2)
00590 latitude = -PI_OVER_2;
00591
00592 if (longitude > PI)
00593 longitude -= TWO_PI;
00594 if (longitude < -PI)
00595 longitude += TWO_PI;
00596
00597 if (longitude > PI)
00598 longitude = PI;
00599 else if (longitude < -PI)
00600 longitude = -PI;
00601
00602 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00603 }
00604
00605
00606 double TransverseCylindricalEqualArea::TCEA_Q( double sinlat, double x )
00607 {
00608 return One_MINUS_es2*(sinlat/(1.0-es2*sinlat*sinlat)-One_OVER_2es*log((1-x)/(1+x)));
00609 }
00610
00611
00612 double TransverseCylindricalEqualArea::TCEA_COEFF_TIMES_SIN( double coeff, double x, double latit )
00613 {
00614 return coeff * sin(x*latit);
00615 }
00616
00617
00618 double TransverseCylindricalEqualArea::TCEA_M( double c0lat, double c1lat, double c2lat, double c3lat )
00619 {
00620 return semiMajorAxis * (c0lat - c1lat + c2lat - c3lat);
00621 }
00622
00623
00624 double TransverseCylindricalEqualArea::TCEA_L( double Beta, double c0lat, double c1lat, double c2lat )
00625 {
00626 return Beta + c0lat + c1lat + c2lat;
00627 }
00628
00629