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 #include <math.h>
00087 #include "Cassini.h"
00088 #include "MapProjection4Parameters.h"
00089 #include "MapProjectionCoordinates.h"
00090 #include "GeodeticCoordinates.h"
00091 #include "CoordinateConversionException.h"
00092 #include "ErrorMessages.h"
00093 #include "WarningMessages.h"
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106 using namespace MSP::CCS;
00107
00108
00109
00110
00111
00112
00113
00114 const double PI = 3.14159265358979323e0;
00115 const double PI_OVER_2 = ( PI / 2.0);
00116 const double TWO_PI = (2.0 * PI);
00117 const double THIRTY_ONE = (31.0 * PI / 180);
00118
00119 double cassCoeffTimesSine( double coeff, double x, double latit )
00120 {
00121 return coeff * (sin (x * latit));
00122 }
00123
00124
00125
00126
00127
00128
00129
00130 Cassini::Cassini( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double originLatitude, double falseEasting, double falseNorthing ) :
00131 CoordinateSystem(),
00132 es2( 0.0066943799901413800 ),
00133 es4( 4.4814723452405e-005 ),
00134 es6( 3.0000678794350e-007 ),
00135 M0( 0.0 ),
00136 c0( .99832429845280 ),
00137 c1( .0025146070605187 ),
00138 c2( 2.6390465943377e-006 ),
00139 c3( 3.4180460865959e-009 ),
00140 One_Minus_es2( .99330562000986 ),
00141 a0( .0025188265843907 ),
00142 a1( 3.7009490356205e-006 ),
00143 a2( 7.4478137675038e-009 ),
00144 a3( 1.7035993238596e-011 ),
00145 Cass_Origin_Long( 0.0 ),
00146 Cass_Origin_Lat( 0.0 ),
00147 Cass_False_Easting( 0.0 ),
00148 Cass_False_Northing( 0.0 ),
00149 Cass_Min_Easting( -20037508.4 ),
00150 Cass_Max_Easting( 20037508.4 ),
00151 Cass_Min_Northing( -56575846.0 ),
00152 Cass_Max_Northing( 56575846.0 )
00153 {
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172 double j,three_es4;
00173 double x, e1, e2, e3, e4;
00174 double lat, sin2lat, sin4lat, sin6lat;
00175 double inv_f = 1 / ellipsoidFlattening;
00176 char errorStatus[500] = "";
00177
00178 if (ellipsoidSemiMajorAxis <= 0.0)
00179 {
00180 strcat( errorStatus, MSP::CCS::ErrorMessages::semiMajorAxis );
00181 }
00182 if ((inv_f < 250) || (inv_f > 350))
00183 {
00184 strcat( errorStatus, MSP::CCS::ErrorMessages::ellipsoidFlattening );
00185 }
00186 if ((originLatitude < -PI_OVER_2) || (originLatitude > PI_OVER_2))
00187 {
00188 strcat( errorStatus, MSP::CCS::ErrorMessages::originLatitude );
00189 }
00190 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00191 {
00192 strcat( errorStatus, MSP::CCS::ErrorMessages::centralMeridian );
00193 }
00194
00195 if( strlen( errorStatus ) > 0)
00196 throw CoordinateConversionException( errorStatus );
00197
00198 semiMajorAxis = ellipsoidSemiMajorAxis;
00199 flattening = ellipsoidFlattening;
00200
00201 Cass_Origin_Lat = originLatitude;
00202 if (centralMeridian > PI)
00203 centralMeridian -= TWO_PI;
00204 Cass_Origin_Long = centralMeridian;
00205 Cass_False_Northing = falseNorthing;
00206 Cass_False_Easting = falseEasting;
00207 es2 = 2 * flattening - flattening * flattening;
00208 es4 = es2 * es2;
00209 es6 = es4 * es2;
00210 j = 45.0 * es6 / 1024.0;
00211 three_es4 = 3.0 * es4;
00212 c0 = 1 - es2 / 4.0 - three_es4 / 64.0 - 5.0 * es6 / 256.0;
00213 c1 = 3.0 * es2 /8.0 + three_es4 / 32.0 + j;
00214 c2 = 15.0 * es4 / 256.0 + j;
00215 c3 = 35.0 * es6 / 3072.0;
00216 lat = c0 * Cass_Origin_Lat;
00217 sin2lat = cassCoeffTimesSine(c1, 2.0, Cass_Origin_Lat);
00218 sin4lat = cassCoeffTimesSine(c2, 4.0, Cass_Origin_Lat);
00219 sin6lat = cassCoeffTimesSine(c3, 6.0, Cass_Origin_Lat);
00220 M0 = cassM( lat, sin2lat, sin4lat, sin6lat );
00221
00222 One_Minus_es2 = 1.0 - es2;
00223 x = sqrt (One_Minus_es2);
00224 e1 = (1 - x) / (1 + x);
00225 e2 = e1 * e1;
00226 e3 = e2 * e1;
00227 e4 = e3 * e1;
00228 a0 = 3.0 * e1 / 2.0 - 27.0 * e3 / 32.0;
00229 a1 = 21.0 * e2 / 16.0 - 55.0 * e4 / 32.0;
00230 a2 = 151.0 * e3 / 96.0;
00231 a3 = 1097.0 * e4 /512.0;
00232
00233 if (Cass_Origin_Long > 0)
00234 {
00235 GeodeticCoordinates gcMax( CoordinateType::geodetic, Cass_Origin_Long - PI, THIRTY_ONE );
00236 MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcMax );
00237 Cass_Max_Northing = tempCoordinates->northing();
00238 delete tempCoordinates;
00239
00240 GeodeticCoordinates gcMin( CoordinateType::geodetic, Cass_Origin_Long - PI, -THIRTY_ONE );
00241 tempCoordinates = convertFromGeodetic( &gcMin );
00242 Cass_Min_Northing = tempCoordinates->northing();
00243 delete tempCoordinates;
00244
00245 Cass_Max_Easting = 19926188.9;
00246 Cass_Min_Easting = -20037508.4;
00247 }
00248 else if (Cass_Origin_Long < 0)
00249 {
00250 GeodeticCoordinates gcMax( CoordinateType::geodetic, PI + Cass_Origin_Long, THIRTY_ONE );
00251 MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcMax );
00252 Cass_Max_Northing = tempCoordinates->northing();
00253 delete tempCoordinates;
00254
00255 GeodeticCoordinates gcMin( CoordinateType::geodetic, PI + Cass_Origin_Long, -THIRTY_ONE );
00256 tempCoordinates = convertFromGeodetic( &gcMin );
00257 Cass_Min_Northing = tempCoordinates->northing();
00258 delete tempCoordinates;
00259
00260 Cass_Max_Easting = 20037508.4;
00261 Cass_Min_Easting = -19926188.9;
00262 }
00263 else
00264 {
00265 GeodeticCoordinates gcMax( CoordinateType::geodetic, PI, THIRTY_ONE );
00266 MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcMax );
00267 Cass_Max_Northing = tempCoordinates->northing();
00268 delete tempCoordinates;
00269
00270 GeodeticCoordinates gcMin( CoordinateType::geodetic, PI, -THIRTY_ONE );
00271 tempCoordinates = convertFromGeodetic( &gcMin );
00272 Cass_Min_Northing = tempCoordinates->northing();
00273 delete tempCoordinates;
00274
00275 Cass_Max_Easting = 20037508.4;
00276 Cass_Min_Easting = -20037508.4;
00277 }
00278
00279 if(Cass_False_Northing)
00280 {
00281 Cass_Min_Northing -= Cass_False_Northing;
00282 Cass_Max_Northing -= Cass_False_Northing;
00283 }
00284 }
00285
00286
00287 Cassini::Cassini( const Cassini &c )
00288 {
00289 semiMajorAxis = c.semiMajorAxis;
00290 flattening = c.flattening;
00291 es2 = c.es2;
00292 es4 = c.es4;
00293 es6 = c.es6;
00294 M0 = c.M0;
00295 c0 = c.c0;
00296 c1 = c.c1;
00297 c2 = c.c2;
00298 c3 = c.c3;
00299 One_Minus_es2 = c.One_Minus_es2;
00300 a0 = c.a0;
00301 a1 = c.a1;
00302 a2 = c.a2;
00303 a3 = c.a3;
00304 Cass_Origin_Long = c.Cass_Origin_Long;
00305 Cass_Origin_Lat = c.Cass_Origin_Lat;
00306 Cass_False_Easting = c.Cass_False_Easting;
00307 Cass_False_Northing = c.Cass_False_Northing;
00308 Cass_Min_Easting = c.Cass_Min_Easting;
00309 Cass_Max_Easting = c.Cass_Max_Easting;
00310 Cass_Min_Northing = c.Cass_Min_Northing;
00311 Cass_Max_Northing = c.Cass_Max_Northing;
00312 }
00313
00314
00315 Cassini::~Cassini()
00316 {
00317 }
00318
00319
00320 Cassini& Cassini::operator=( const Cassini &c )
00321 {
00322 if( this != &c )
00323 {
00324 semiMajorAxis = c.semiMajorAxis;
00325 flattening = c.flattening;
00326 es2 = c.es2;
00327 es4 = c.es4;
00328 es6 = c.es6;
00329 M0 = c.M0;
00330 c0 = c.c0;
00331 c1 = c.c1;
00332 c2 = c.c2;
00333 c3 = c.c3;
00334 One_Minus_es2 = c.One_Minus_es2;
00335 a0 = c.a0;
00336 a1 = c.a1;
00337 a2 = c.a2;
00338 a3 = c.a3;
00339 Cass_Origin_Long = c.Cass_Origin_Long;
00340 Cass_Origin_Lat = c.Cass_Origin_Lat;
00341 Cass_False_Easting = c.Cass_False_Easting;
00342 Cass_False_Northing = c.Cass_False_Northing;
00343 Cass_Min_Easting = c.Cass_Min_Easting;
00344 Cass_Max_Easting = c.Cass_Max_Easting;
00345 Cass_Min_Northing = c.Cass_Min_Northing;
00346 Cass_Max_Northing = c.Cass_Max_Northing;
00347 }
00348
00349 return *this;
00350 }
00351
00352
00353 MapProjection4Parameters* Cassini::getParameters() const
00354 {
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371 return new MapProjection4Parameters( CoordinateType::cassini, Cass_Origin_Long, Cass_Origin_Lat, Cass_False_Easting, Cass_False_Northing );
00372 }
00373
00374
00375 MSP::CCS::MapProjectionCoordinates* Cassini::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00376 {
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390 double lat, sin2lat, sin4lat, sin6lat;
00391 double RD;
00392 double dlam;
00393 double NN;
00394 double TT;
00395 double AA, A2, A3, A4, A5;
00396 double CC;
00397 double MM;
00398 char errorStatus[256] = "";
00399
00400 double longitude = geodeticCoordinates->longitude();
00401 double latitude = geodeticCoordinates->latitude();
00402 double tlat = tan(latitude);
00403 double clat = cos(latitude);
00404 double slat = sin(latitude);
00405
00406 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00407 {
00408 strcat( errorStatus, MSP::CCS::ErrorMessages::latitude );
00409 }
00410 if ((longitude < -PI) || (longitude > TWO_PI))
00411 {
00412 strcat( errorStatus, MSP::CCS::ErrorMessages::longitude );
00413 }
00414
00415 if( strlen( errorStatus ) > 0)
00416 throw CoordinateConversionException( errorStatus );
00417
00418 dlam = longitude - Cass_Origin_Long;
00419
00420 if (fabs(dlam) > (4.0 * PI / 180))
00421 {
00422 strcat( errorStatus, MSP::CCS::WarningMessages::longitude );
00423 }
00424
00425 if (dlam > PI)
00426 {
00427 dlam -= TWO_PI;
00428 }
00429 if (dlam < -PI)
00430 {
00431 dlam += TWO_PI;
00432 }
00433 RD = cassRd( slat );
00434 NN = semiMajorAxis / RD;
00435 TT = tlat * tlat;
00436 AA = dlam * clat;
00437 A2 = AA * AA;
00438 A3 = AA * A2;
00439 A4 = AA * A3;
00440 A5 = AA * A4;
00441 CC = es2 * clat * clat / One_Minus_es2;
00442 lat = c0 * latitude;
00443 sin2lat = cassCoeffTimesSine(c1, 2.0, latitude);
00444 sin4lat = cassCoeffTimesSine(c2, 4.0, latitude);
00445 sin6lat = cassCoeffTimesSine(c3, 6.0, latitude);
00446 MM = cassM( lat, sin2lat, sin4lat, sin6lat );
00447
00448 double easting = NN * (AA - (TT * A3 / 6.0) - (8.0 - TT + 8.0 * CC) *
00449 (TT * A5 / 120.0)) + Cass_False_Easting;
00450 double northing = MM - M0 + NN * tlat * ((A2 / 2.0) + (5.0 - TT +
00451 6.0 * CC) * A4 / 24.0) + Cass_False_Northing;
00452
00453 return new MapProjectionCoordinates( CoordinateType::cassini, errorStatus, easting, northing );
00454 }
00455
00456
00457 MSP::CCS::GeodeticCoordinates* Cassini::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00458 {
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472 double dx;
00473 double dy;
00474 double mu1;
00475 double sin2mu, sin4mu, sin6mu, sin8mu;
00476 double M1;
00477 double phi1;
00478 double tanphi1, sinphi1, cosphi1;
00479 double T1, T;
00480 double N1;
00481 double RD, R1;
00482 double DD, D2, D3, D4, D5;
00483 const double epsilon = 1.0e-1;
00484 double longitude, latitude;
00485 char errorStatus[256] = "";
00486
00487 double easting = mapProjectionCoordinates->easting();
00488 double northing = mapProjectionCoordinates->northing();
00489
00490 if ((easting < (Cass_False_Easting + Cass_Min_Easting))
00491 || (easting > (Cass_False_Easting + Cass_Max_Easting)))
00492 {
00493 strcat( errorStatus, MSP::CCS::ErrorMessages::easting );
00494 }
00495 if ((northing < (Cass_False_Northing + Cass_Min_Northing - epsilon))
00496 || (northing > (Cass_False_Northing + Cass_Max_Northing + epsilon)))
00497 {
00498 strcat( errorStatus, MSP::CCS::ErrorMessages::northing );
00499 }
00500
00501 if( strlen( errorStatus ) > 0)
00502 throw CoordinateConversionException( errorStatus );
00503
00504 dy = northing - Cass_False_Northing;
00505 dx = easting - Cass_False_Easting;
00506 M1 = M0 + dy;
00507 mu1 = M1 / (semiMajorAxis * c0);
00508
00509 sin2mu = cassCoeffTimesSine(a0, 2.0, mu1);
00510 sin4mu = cassCoeffTimesSine(a1, 4.0, mu1);
00511 sin6mu = cassCoeffTimesSine(a2, 6.0, mu1);
00512 sin8mu = cassCoeffTimesSine(a3, 8.0, mu1);
00513 phi1 = mu1 + sin2mu + sin4mu + sin6mu + sin8mu;
00514
00515 if (floatEq(phi1,PI_OVER_2,.00001))
00516 {
00517 latitude = PI_OVER_2;
00518 longitude = Cass_Origin_Long;
00519 }
00520 else if (floatEq(phi1,-PI_OVER_2,.00001))
00521 {
00522 latitude = -PI_OVER_2;
00523 longitude = Cass_Origin_Long;
00524 }
00525 else
00526 {
00527 tanphi1 = tan(phi1);
00528 sinphi1 = sin(phi1);
00529 cosphi1 = cos(phi1);
00530 T1 = tanphi1 * tanphi1;
00531 RD = cassRd( sinphi1 );
00532 N1 = semiMajorAxis / RD;
00533 R1 = N1 * One_Minus_es2 / (RD * RD);
00534 DD = dx / N1;
00535 D2 = DD * DD;
00536 D3 = D2 * DD;
00537 D4 = D3 * DD;
00538 D5 = D4 * DD;
00539 T = (1.0 + 3.0 * T1);
00540 latitude = phi1 - (N1 * tanphi1 / R1) * (D2 / 2.0 - T * D4 / 24.0);
00541
00542 longitude = Cass_Origin_Long + (DD - T1 * D3 / 3.0 + T * T1 * D5 / 15.0) / cosphi1;
00543
00544 if (latitude > PI_OVER_2)
00545 latitude = PI_OVER_2;
00546 else if (latitude < -PI_OVER_2)
00547 latitude = -PI_OVER_2;
00548
00549 if (longitude > PI)
00550 longitude -= TWO_PI;
00551 if (longitude < -PI)
00552 longitude += TWO_PI;
00553
00554 if (longitude > PI)
00555 longitude = PI;
00556 else if (longitude < -PI)
00557 longitude = -PI;
00558 }
00559 if (fabs(longitude - Cass_Origin_Long) > (4.0 * PI / 180))
00560 {
00561 strcat( errorStatus, MSP::CCS::WarningMessages::longitude );
00562 }
00563
00564 return new GeodeticCoordinates( CoordinateType::geodetic, errorStatus, longitude, latitude );
00565 }
00566
00567
00568 double Cassini::cassM( double c0lat, double c1s2lat, double c2s4lat, double c3s6lat )
00569 {
00570 return semiMajorAxis*(c0lat-c1s2lat+c2s4lat-c3s6lat);
00571 }
00572
00573
00574 double Cassini::cassRd( double sinlat )
00575 {
00576 return sqrt(1.0 - es2 * (sinlat * sinlat));
00577 }
00578
00579
00580 double Cassini::floatEq( double x, double v, double epsilon )
00581 {
00582 return ((v - epsilon) < x) && (x < (v + epsilon));
00583 }
00584
00585
00586
00587