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