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 #include <stdio.h>
00085 #include <math.h>
00086 #include "Mercator.h"
00087 #include "MercatorStandardParallelParameters.h"
00088 #include "MercatorScaleFactorParameters.h"
00089 #include "MapProjectionCoordinates.h"
00090 #include "GeodeticCoordinates.h"
00091 #include "CoordinateConversionException.h"
00092 #include "ErrorMessages.h"
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104 using namespace MSP::CCS;
00105
00106
00107
00108
00109
00110
00111
00112 const double PI = 3.14159265358979323e0;
00113 const double PI_OVER_2 = ( PI / 2.0e0);
00114 const double TWO_PI = (2.0 * PI);
00115 const double MAX_LAT = ( (PI * 89.5) / 180.0 );
00116 const double MIN_SCALE_FACTOR = 0.3;
00117 const double MAX_SCALE_FACTOR = 3.0;
00118
00119
00120
00121
00122
00123
00124
00125 Mercator::Mercator( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double standardParallel, double falseEasting, double falseNorthing, double* scaleFactor ) :
00126 CoordinateSystem(),
00127 coordinateType( CoordinateType::mercatorStandardParallel ),
00128 Merc_e( 0.08181919084262188000 ),
00129 Merc_es( 0.0066943799901413800 ),
00130 Merc_Cent_Mer( 0.0 ),
00131 Merc_Standard_Parallel( 0.0 ),
00132 Merc_False_Easting( 0.0 ),
00133 Merc_False_Northing( 0.0 ),
00134 Merc_Scale_Factor( 1.0 ),
00135 Merc_ab( 0.00335655146887969400 ),
00136 Merc_bb( 0.00000657187271079536 ),
00137 Merc_cb( 0.00000001764564338702 ),
00138 Merc_db( 0.00000000005328478445 ),
00139 Merc_Delta_Easting( 20237883.0 ),
00140 Merc_Delta_Northing( 23421740.0 )
00141 {
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163 double es2;
00164 double es3;
00165 double es4;
00166 double sin_olat;
00167 double inv_f = 1 / ellipsoidFlattening;
00168 char errorStatus[500] = "";
00169
00170 if (ellipsoidSemiMajorAxis <= 0.0)
00171 {
00172 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00173 }
00174 if ((inv_f < 250) || (inv_f > 350))
00175 {
00176 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00177 }
00178 if ((standardParallel < -MAX_LAT) || (standardParallel > MAX_LAT))
00179 {
00180 strcat( errorStatus, ErrorMessages::originLatitude );
00181 }
00182 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00183 {
00184 strcat( errorStatus, ErrorMessages::centralMeridian );
00185 }
00186
00187 if( strlen( errorStatus ) > 0)
00188 throw CoordinateConversionException( errorStatus );
00189
00190 semiMajorAxis = ellipsoidSemiMajorAxis;
00191 flattening = ellipsoidFlattening;
00192
00193 Merc_Standard_Parallel = standardParallel;
00194 if (centralMeridian > PI)
00195 centralMeridian -= TWO_PI;
00196 Merc_Cent_Mer = centralMeridian;
00197 Merc_False_Northing = falseNorthing;
00198 Merc_False_Easting = falseEasting;
00199
00200 Merc_es = 2 * flattening - flattening * flattening;
00201 Merc_e = sqrt(Merc_es);
00202 sin_olat = sin(Merc_Standard_Parallel);
00203 Merc_Scale_Factor = cos(Merc_Standard_Parallel) / sqrt(1.e0 - Merc_es * sin_olat * sin_olat );
00204
00205
00206 es2 = Merc_es * Merc_es;
00207 es3 = es2 * Merc_es;
00208 es4 = es3 * Merc_es;
00209 Merc_ab = Merc_es / 2.e0 + 5.e0 * es2 / 24.e0 + es3 / 12.e0
00210 + 13.e0 * es4 / 360.e0;
00211 Merc_bb = 7.e0 * es2 / 48.e0 + 29.e0 * es3 / 240.e0
00212 + 811.e0 * es4 / 11520.e0;
00213 Merc_cb = 7.e0 * es3 / 120.e0 + 81.e0 * es4 / 1120.e0;
00214 Merc_db = 4279.e0 * es4 / 161280.e0;
00215 *scaleFactor = Merc_Scale_Factor;
00216
00217
00218
00219
00220
00221
00222 GeodeticCoordinates gcTemp( CoordinateType::geodetic, ( Merc_Cent_Mer + PI ), MAX_LAT );
00223 MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcTemp );
00224 Merc_Delta_Easting = tempCoordinates->easting();
00225 Merc_Delta_Northing = tempCoordinates->northing();
00226 delete tempCoordinates;
00227
00228 if(Merc_False_Easting)
00229 Merc_Delta_Easting -= Merc_False_Easting;
00230 if (Merc_Delta_Easting < 0)
00231 Merc_Delta_Easting = -Merc_Delta_Easting;
00232 Merc_Delta_Easting *= 1.01;
00233
00234 if(Merc_False_Northing)
00235 Merc_Delta_Northing -= Merc_False_Northing;
00236 if (Merc_Delta_Northing < 0)
00237 Merc_Delta_Northing = -Merc_Delta_Northing;
00238 Merc_Delta_Northing *= 1.01;
00239 }
00240
00241
00242 Mercator::Mercator( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double falseEasting, double falseNorthing, double scaleFactor ) :
00243 CoordinateSystem(),
00244 coordinateType( CoordinateType::mercatorScaleFactor ),
00245 Merc_e( 0.08181919084262188000 ),
00246 Merc_es( 0.0066943799901413800 ),
00247 Merc_Cent_Mer( 0.0 ),
00248 Merc_Standard_Parallel( 0.0 ),
00249 Merc_False_Easting( 0.0 ),
00250 Merc_False_Northing( 0.0 ),
00251 Merc_Scale_Factor( 1.0 ),
00252 Merc_ab( 0.00335655146887969400 ),
00253 Merc_bb( 0.00000657187271079536 ),
00254 Merc_cb( 0.00000001764564338702 ),
00255 Merc_db( 0.00000000005328478445 ),
00256 Merc_Delta_Easting( 20237883.0 ),
00257 Merc_Delta_Northing( 23421740.0 )
00258 {
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278 double es2;
00279 double es3;
00280 double es4;
00281 double Merc_Scale_Factor2;
00282 double inv_f = 1 / ellipsoidFlattening;
00283 char errorStatus[500] = "";
00284
00285 if (ellipsoidSemiMajorAxis <= 0.0)
00286 {
00287 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00288 }
00289 if ((inv_f < 250) || (inv_f > 350))
00290 {
00291 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00292 }
00293 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00294 {
00295 strcat( errorStatus, ErrorMessages::centralMeridian );
00296 }
00297 if ((scaleFactor < MIN_SCALE_FACTOR) || (scaleFactor > MAX_SCALE_FACTOR))
00298 {
00299 strcat( errorStatus, ErrorMessages::scaleFactor );
00300 }
00301
00302 if( strlen( errorStatus ) > 0)
00303 throw CoordinateConversionException( errorStatus );
00304
00305 semiMajorAxis = ellipsoidSemiMajorAxis;
00306 flattening = ellipsoidFlattening;
00307
00309 if (centralMeridian > PI)
00310 centralMeridian -= TWO_PI;
00311 Merc_Cent_Mer = centralMeridian;
00312 Merc_False_Northing = falseNorthing;
00313 Merc_False_Easting = falseEasting;
00314 Merc_Scale_Factor = scaleFactor;
00315
00316 Merc_es = 2 * flattening - flattening * flattening;
00317 Merc_e = sqrt(Merc_es);
00318
00319
00320
00321
00322 es2 = Merc_es * Merc_es;
00323 es3 = es2 * Merc_es;
00324 es4 = es3 * Merc_es;
00325 Merc_ab = Merc_es / 2.e0 + 5.e0 * es2 / 24.e0 + es3 / 12.e0
00326 + 13.e0 * es4 / 360.e0;
00327 Merc_bb = 7.e0 * es2 / 48.e0 + 29.e0 * es3 / 240.e0
00328 + 811.e0 * es4 / 11520.e0;
00329 Merc_cb = 7.e0 * es3 / 120.e0 + 81.e0 * es4 / 1120.e0;
00330 Merc_db = 4279.e0 * es4 / 161280.e0;
00331
00332
00333 Merc_Scale_Factor2 = Merc_Scale_Factor * Merc_Scale_Factor;
00334 Merc_Standard_Parallel = asin(sqrt((1 - Merc_Scale_Factor2)/(1 - Merc_Scale_Factor2 * Merc_es)));
00335
00336
00337
00338
00339
00340
00341 GeodeticCoordinates gcTemp( CoordinateType::geodetic, ( Merc_Cent_Mer + PI ), MAX_LAT );
00342 MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcTemp );
00343 Merc_Delta_Easting = tempCoordinates->easting();
00344 Merc_Delta_Northing = tempCoordinates->northing();
00345 delete tempCoordinates;
00346
00347 if(Merc_False_Easting)
00348 Merc_Delta_Easting -= Merc_False_Easting;
00349 if (Merc_Delta_Easting < 0)
00350 Merc_Delta_Easting = -Merc_Delta_Easting;
00351 Merc_Delta_Easting *= 1.01;
00352
00353 if(Merc_False_Northing)
00354 Merc_Delta_Northing -= Merc_False_Northing;
00355 if (Merc_Delta_Northing < 0)
00356 Merc_Delta_Northing = -Merc_Delta_Northing;
00357 Merc_Delta_Northing *= 1.01;
00358 }
00359
00360
00361 Mercator::Mercator( const Mercator &m )
00362 {
00363 coordinateType = m.coordinateType;
00364 semiMajorAxis = m.semiMajorAxis;
00365 flattening = m.flattening;
00366 Merc_e = m.Merc_e;
00367 Merc_es = m.Merc_es;
00368 Merc_Cent_Mer = m.Merc_Cent_Mer;
00369 Merc_Standard_Parallel = m.Merc_Standard_Parallel;
00370 Merc_False_Easting = m.Merc_False_Easting;
00371 Merc_False_Northing = m.Merc_False_Northing;
00372 Merc_Scale_Factor = m.Merc_Scale_Factor;
00373 Merc_ab = m.Merc_ab;
00374 Merc_bb = m.Merc_bb;
00375 Merc_cb = m.Merc_cb;
00376 Merc_db = m.Merc_db;
00377 Merc_Delta_Easting = m.Merc_Delta_Easting;
00378 Merc_Delta_Northing = m.Merc_Delta_Northing;
00379 }
00380
00381
00382 Mercator::~Mercator()
00383 {
00384 }
00385
00386
00387 Mercator& Mercator::operator=( const Mercator &m )
00388 {
00389 if( this != &m )
00390 {
00391 coordinateType = m.coordinateType;
00392 semiMajorAxis = m.semiMajorAxis;
00393 flattening = m.flattening;
00394 Merc_e = m.Merc_e;
00395 Merc_es = m.Merc_es;
00396 Merc_Cent_Mer = m.Merc_Cent_Mer;
00397 Merc_Standard_Parallel = m.Merc_Standard_Parallel;
00398 Merc_False_Easting = m.Merc_False_Easting;
00399 Merc_False_Northing = m.Merc_False_Northing;
00400 Merc_Scale_Factor = m.Merc_Scale_Factor;
00401 Merc_ab = m.Merc_ab;
00402 Merc_bb = m.Merc_bb;
00403 Merc_cb = m.Merc_cb;
00404 Merc_db = m.Merc_db;
00405 Merc_Delta_Easting = m.Merc_Delta_Easting;
00406 Merc_Delta_Northing = m.Merc_Delta_Northing;
00407 }
00408
00409 return *this;
00410 }
00411
00412
00413 MercatorStandardParallelParameters* Mercator::getStandardParallelParameters() const
00414 {
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431 return new MercatorStandardParallelParameters(CoordinateType::mercatorStandardParallel, Merc_Cent_Mer, Merc_Standard_Parallel, Merc_Scale_Factor, Merc_False_Easting, Merc_False_Northing);
00432 }
00433
00434
00435 MercatorScaleFactorParameters* Mercator::getScaleFactorParameters() const
00436 {
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454 return new MercatorScaleFactorParameters( CoordinateType::mercatorScaleFactor, Merc_Cent_Mer, Merc_Scale_Factor, Merc_False_Easting, Merc_False_Northing );
00455 }
00456
00457
00458 MSP::CCS::MapProjectionCoordinates* Mercator::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00459 {
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473 double ctanz2;
00474 double e_x_sinlat;
00475 double Delta_Long;
00476 double tan_temp;
00477 double pow_temp;
00478 char errorStatus[50] = "";
00479
00480 double longitude = geodeticCoordinates->longitude();
00481 double latitude = geodeticCoordinates->latitude();
00482
00483 if ((latitude < -MAX_LAT) || (latitude > MAX_LAT))
00484 {
00485 strcat( errorStatus, ErrorMessages::latitude );
00486 }
00487 if ((longitude < -PI) || (longitude > TWO_PI))
00488 {
00489 strcat( errorStatus, ErrorMessages::longitude );
00490 }
00491
00492 if( strlen( errorStatus ) > 0)
00493 throw CoordinateConversionException( errorStatus );
00494
00495 if (longitude > PI)
00496 longitude -= TWO_PI;
00497 e_x_sinlat = Merc_e * sin(latitude);
00498 tan_temp = tan(PI / 4.e0 + latitude / 2.e0);
00499 pow_temp = pow( ((1.e0 - e_x_sinlat) / (1.e0 + e_x_sinlat)),
00500 (Merc_e / 2.e0) );
00501 ctanz2 = tan_temp * pow_temp;
00502 double northing = Merc_Scale_Factor * semiMajorAxis * log(ctanz2) + Merc_False_Northing;
00503 Delta_Long = longitude - Merc_Cent_Mer;
00504 if (Delta_Long > PI)
00505 Delta_Long -= TWO_PI;
00506 if (Delta_Long < -PI)
00507 Delta_Long += TWO_PI;
00508 double easting = Merc_Scale_Factor * semiMajorAxis * Delta_Long
00509 + Merc_False_Easting;
00510
00511 return new MapProjectionCoordinates( coordinateType, easting, northing );
00512 }
00513
00514
00515 MSP::CCS::GeodeticCoordinates* Mercator::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00516 {
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530 double dx;
00531 double dy;
00532 double xphi;
00533 char errorStatus[50] = "";
00534
00535 double easting = mapProjectionCoordinates->easting();
00536 double northing = mapProjectionCoordinates->northing();
00537
00538 if ((easting < (Merc_False_Easting - Merc_Delta_Easting))
00539 || (easting > (Merc_False_Easting + Merc_Delta_Easting)))
00540 {
00541 strcat( errorStatus, ErrorMessages::easting );
00542 }
00543 if ((northing < (Merc_False_Northing - Merc_Delta_Northing))
00544 || (northing > (Merc_False_Northing + Merc_Delta_Northing)))
00545 {
00546 strcat( errorStatus, ErrorMessages::northing );
00547 }
00548
00549 if( strlen( errorStatus ) > 0)
00550 throw CoordinateConversionException( errorStatus );
00551
00552 dy = northing - Merc_False_Northing;
00553 dx = easting - Merc_False_Easting;
00554 double longitude = Merc_Cent_Mer + dx / (Merc_Scale_Factor * semiMajorAxis);
00555 xphi = PI_OVER_2
00556 - 2.e0 * atan(1.e0 / exp(dy / (Merc_Scale_Factor * semiMajorAxis)));
00557 double latitude = xphi + Merc_ab * sin(2.e0 * xphi) + Merc_bb * sin(4.e0 * xphi)
00558 + Merc_cb * sin(6.e0 * xphi) + Merc_db * sin(8.e0 * xphi);
00559
00560 if (longitude > PI)
00561 longitude -= TWO_PI;
00562 if (longitude < -PI)
00563 longitude += TWO_PI;
00564
00565 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00566 }
00567
00568
00569
00570