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
00090
00091
00092
00093
00094
00095
00096 #include <math.h>
00097 #include "ObliqueMercator.h"
00098 #include "ObliqueMercatorParameters.h"
00099 #include "MapProjectionCoordinates.h"
00100 #include "GeodeticCoordinates.h"
00101 #include "CoordinateConversionException.h"
00102 #include "ErrorMessages.h"
00103 #include "WarningMessages.h"
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116 using namespace MSP::CCS;
00117
00118
00119
00120
00121
00122
00123
00124 const double PI = 3.14159265358979323e0;
00125 const double PI_OVER_2 = ( PI / 2.0);
00126 const double PI_OVER_4 = ( PI / 4.0);
00127 const double TWO_PI = ( 2.0 * PI);
00128 const double MIN_SCALE_FACTOR = 0.3;
00129 const double MAX_SCALE_FACTOR = 3.0;
00130
00131
00132
00133
00134
00135
00136
00137 ObliqueMercator::ObliqueMercator( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double originLatitude, double longitude1, double latitude1, double longitude2, double latitude2, double falseEasting, double falseNorthing, double scaleFactor ) :
00138 CoordinateSystem(),
00139 es( 0.08181919084262188000 ),
00140 es_OVER_2( .040909595421311 ),
00141 OMerc_A( 6383471.9177251 ),
00142 OMerc_B( 1.0008420825413 ),
00143 OMerc_E( 1.0028158089754 ),
00144 OMerc_gamma( .41705894983580 ),
00145 OMerc_azimuth( .60940407333533 ),
00146 OMerc_Origin_Long( -.46732023406900 ),
00147 cos_gamma( .91428423352628 ),
00148 sin_gamma( .40507325303611 ),
00149 sin_azimuth( .57237890829911 ),
00150 cos_azimuth( .81998925927985 ),
00151 A_over_B( 6378101.0302010 ),
00152 B_over_A( 1.5678647849335e-7 ),
00153 OMerc_u( 5632885.2272051 ),
00154 OMerc_Origin_Lat( ((45.0 * PI) / 180.0) ),
00155 OMerc_Lon_1( ((-5.0 * PI) / 180.0) ),
00156 OMerc_Lat_1( ((40.0 * PI) / 180.0) ),
00157 OMerc_Lon_2( ((5.0 * PI) / 180.0) ),
00158 OMerc_Lat_2( ((50.0 * PI) / 180.0) ),
00159 OMerc_False_Easting( 0.0 ),
00160 OMerc_False_Northing( 0.0 ),
00161 OMerc_Scale_Factor( 1.0 ),
00162 OMerc_Delta_Northing( 40000000.0 ),
00163 OMerc_Delta_Easting( 40000000.0 )
00164 {
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 double inv_f = 1 / ellipsoidFlattening;
00193 double es2, one_MINUS_es2;
00194 double cos_olat, cos_olat2;
00195 double sin_olat, sin_olat2, es2_sin_olat2;
00196 double t0, t1, t2;
00197 double D, D2, D2_MINUS_1, sqrt_D2_MINUS_1;
00198 double H, L, LH;
00199 double E2;
00200 double F, G, J, P;
00201 double dlon;
00202 char errorStatus[500] = "";
00203
00204 if (ellipsoidSemiMajorAxis <= 0.0)
00205 {
00206 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00207 }
00208 if ((inv_f < 250) || (inv_f > 350))
00209 {
00210 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00211 }
00212 if ((originLatitude <= -PI_OVER_2) || (originLatitude >= PI_OVER_2))
00213 {
00214 strcat( errorStatus, ErrorMessages::originLatitude );
00215 }
00216 if ((latitude1 <= -PI_OVER_2) || (latitude1 >= PI_OVER_2))
00217 {
00218 strcat( errorStatus, ErrorMessages::latitude1 );
00219 }
00220 if ((latitude2 <= -PI_OVER_2) || (latitude2 >= PI_OVER_2))
00221 {
00222 strcat( errorStatus, ErrorMessages::latitude2 );
00223 }
00224 if (latitude1 == 0.0)
00225 {
00226 strcat( errorStatus, ErrorMessages::latitude1 );
00227 }
00228 if (latitude1 == latitude2)
00229 {
00230 strcat( errorStatus, ErrorMessages::latitude2 );
00231 }
00232 if (((latitude1 < 0.0) && (latitude2 > 0.0)) ||
00233 ((latitude1 > 0.0) && (latitude2 < 0.0)))
00234 {
00235 strcat( errorStatus, ErrorMessages::omercHemisphere );
00236 }
00237 if ((longitude1 < -PI) || (longitude1 > TWO_PI))
00238 {
00239 strcat( errorStatus, ErrorMessages::longitude1 );
00240 }
00241 if ((longitude2 < -PI) || (longitude2 > TWO_PI))
00242 {
00243 strcat( errorStatus, ErrorMessages::longitude2 );
00244 }
00245 if ((scaleFactor < MIN_SCALE_FACTOR) || (scaleFactor > MAX_SCALE_FACTOR))
00246 {
00247 strcat( errorStatus, ErrorMessages::scaleFactor );
00248 }
00249
00250 if( strlen( errorStatus ) > 0)
00251 throw CoordinateConversionException( errorStatus );
00252
00253 semiMajorAxis = ellipsoidSemiMajorAxis;
00254 flattening = ellipsoidFlattening;
00255
00256 OMerc_Origin_Lat = originLatitude;
00257 OMerc_Lon_1 = longitude1;
00258 OMerc_Lat_1 = latitude1;
00259 OMerc_Lon_2 = longitude2;
00260 OMerc_Lat_2 = latitude2;
00261 OMerc_False_Northing = falseNorthing;
00262 OMerc_False_Easting = falseEasting;
00263 OMerc_Scale_Factor = scaleFactor;
00264
00265 es2 = 2 * flattening - flattening * flattening;
00266 es = sqrt(es2);
00267 one_MINUS_es2 = 1 - es2;
00268 es_OVER_2 = es / 2.0;
00269
00270 cos_olat = cos(OMerc_Origin_Lat);
00271 cos_olat2 = cos_olat * cos_olat;
00272 sin_olat = sin(OMerc_Origin_Lat);
00273 sin_olat2 = sin_olat * sin_olat;
00274 es2_sin_olat2 = es2 * sin_olat2;
00275
00276 OMerc_B = sqrt(1 + (es2 * cos_olat2 * cos_olat2) / one_MINUS_es2);
00277 OMerc_A = (semiMajorAxis * OMerc_B * OMerc_Scale_Factor * sqrt(one_MINUS_es2)) / (1.0 - es2_sin_olat2);
00278 A_over_B = OMerc_A / OMerc_B;
00279 B_over_A = OMerc_B / OMerc_A;
00280
00281 t0 = omercT(OMerc_Origin_Lat, es * sin_olat, es_OVER_2);
00282 t1 = omercT(OMerc_Lat_1, es * sin(OMerc_Lat_1), es_OVER_2);
00283 t2 = omercT(OMerc_Lat_2, es * sin(OMerc_Lat_2), es_OVER_2);
00284
00285 D = (OMerc_B * sqrt(one_MINUS_es2)) / (cos_olat * sqrt(1.0 - es2_sin_olat2));
00286 D2 = D * D;
00287 if (D2 < 1.0)
00288 D2 = 1.0;
00289 D2_MINUS_1 = D2 - 1.0;
00290 sqrt_D2_MINUS_1 = sqrt(D2_MINUS_1);
00291 if (D2_MINUS_1 > 1.0e-10)
00292 {
00293 if (OMerc_Origin_Lat >= 0.0)
00294 OMerc_E = (D + sqrt_D2_MINUS_1) * pow(t0, OMerc_B);
00295 else
00296 OMerc_E = (D - sqrt_D2_MINUS_1) * pow(t0, OMerc_B);
00297 }
00298 else
00299 OMerc_E = D * pow(t0, OMerc_B);
00300 H = pow(t1, OMerc_B);
00301 L = pow(t2, OMerc_B);
00302 F = OMerc_E / H;
00303 G = (F - 1.0 / F) / 2.0;
00304 E2 = OMerc_E * OMerc_E;
00305 LH = L * H;
00306 J = (E2 - LH) / (E2 + LH);
00307 P = (L - H) / (L + H);
00308
00309 dlon = OMerc_Lon_1 - OMerc_Lon_2;
00310 if (dlon < -PI )
00311 OMerc_Lon_2 -= TWO_PI;
00312 if (dlon > PI)
00313 OMerc_Lon_2 += TWO_PI;
00314 dlon = OMerc_Lon_1 - OMerc_Lon_2;
00315 OMerc_Origin_Long = (OMerc_Lon_1 + OMerc_Lon_2) / 2.0 - (atan(J * tan(OMerc_B * dlon / 2.0) / P)) / OMerc_B;
00316
00317 dlon = OMerc_Lon_1 - OMerc_Origin_Long;
00318 if (dlon < -PI )
00319 OMerc_Origin_Long -= TWO_PI;
00320 if (dlon > PI)
00321 OMerc_Origin_Long += TWO_PI;
00322
00323 dlon = OMerc_Lon_1 - OMerc_Origin_Long;
00324 OMerc_gamma = atan(sin(OMerc_B * dlon) / G);
00325 cos_gamma = cos(OMerc_gamma);
00326 sin_gamma = sin(OMerc_gamma);
00327
00328 OMerc_azimuth = asin(D * sin_gamma);
00329 cos_azimuth = cos(OMerc_azimuth);
00330 sin_azimuth = sin(OMerc_azimuth);
00331
00332 if (OMerc_Origin_Lat >= 0)
00333 OMerc_u = A_over_B * atan(sqrt_D2_MINUS_1/cos_azimuth);
00334 else
00335 OMerc_u = -A_over_B * atan(sqrt_D2_MINUS_1/cos_azimuth);
00336 }
00337
00338
00339 ObliqueMercator::ObliqueMercator( const ObliqueMercator &om )
00340 {
00341 semiMajorAxis = om.semiMajorAxis;
00342 flattening = om.flattening;
00343 es = om.es;
00344 es_OVER_2 = om.es_OVER_2;
00345 OMerc_A = om.OMerc_A;
00346 OMerc_B = om.OMerc_B;
00347 OMerc_E = om.OMerc_E;
00348 OMerc_gamma = om.OMerc_gamma;
00349 OMerc_azimuth = om.OMerc_azimuth;
00350 OMerc_Origin_Long = om.OMerc_Origin_Long;
00351 cos_gamma = om.cos_gamma;
00352 sin_gamma = om.sin_gamma;
00353 sin_azimuth = om.sin_azimuth;
00354 cos_azimuth = om.cos_azimuth;
00355 A_over_B = om.A_over_B;
00356 B_over_A = om.B_over_A;
00357 OMerc_u = om.OMerc_u;
00358 OMerc_Origin_Lat = om.OMerc_Origin_Lat;
00359 OMerc_Lon_1 = om.OMerc_Lon_1;
00360 OMerc_Lat_1 = om.OMerc_Lat_1;
00361 OMerc_Lon_2 = om.OMerc_Lon_2;
00362 OMerc_Lat_2 = om.OMerc_Lat_2;
00363 OMerc_False_Easting = om.OMerc_False_Easting;
00364 OMerc_False_Northing = om.OMerc_False_Northing;
00365 OMerc_Scale_Factor = om.OMerc_Scale_Factor;
00366 OMerc_Delta_Northing = om.OMerc_Delta_Northing;
00367 OMerc_Delta_Easting = om.OMerc_Delta_Easting;
00368 }
00369
00370
00371 ObliqueMercator::~ObliqueMercator()
00372 {
00373 }
00374
00375
00376 ObliqueMercator& ObliqueMercator::operator=( const ObliqueMercator &om )
00377 {
00378 if( this != &om )
00379 {
00380 semiMajorAxis = om.semiMajorAxis;
00381 flattening = om.flattening;
00382 es = om.es;
00383 es_OVER_2 = om.es_OVER_2;
00384 OMerc_A = om.OMerc_A;
00385 OMerc_B = om.OMerc_B;
00386 OMerc_E = om.OMerc_E;
00387 OMerc_gamma = om.OMerc_gamma;
00388 OMerc_azimuth = om.OMerc_azimuth;
00389 OMerc_Origin_Long = om.OMerc_Origin_Long;
00390 cos_gamma = om.cos_gamma;
00391 sin_gamma = om.sin_gamma;
00392 sin_azimuth = om.sin_azimuth;
00393 cos_azimuth = om.cos_azimuth;
00394 A_over_B = om.A_over_B;
00395 B_over_A = om.B_over_A;
00396 OMerc_u = om.OMerc_u;
00397 OMerc_Origin_Lat = om.OMerc_Origin_Lat;
00398 OMerc_Lon_1 = om.OMerc_Lon_1;
00399 OMerc_Lat_1 = om.OMerc_Lat_1;
00400 OMerc_Lon_2 = om.OMerc_Lon_2;
00401 OMerc_Lat_2 = om.OMerc_Lat_2;
00402 OMerc_False_Easting = om.OMerc_False_Easting;
00403 OMerc_False_Northing = om.OMerc_False_Northing;
00404 OMerc_Scale_Factor = om.OMerc_Scale_Factor;
00405 OMerc_Delta_Northing = om.OMerc_Delta_Northing;
00406 OMerc_Delta_Easting = om.OMerc_Delta_Easting;
00407 }
00408
00409 return *this;
00410 }
00411
00412
00413 ObliqueMercatorParameters* ObliqueMercator::getParameters() const
00414 {
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440 return new ObliqueMercatorParameters( CoordinateType::obliqueMercator, OMerc_Origin_Lat, OMerc_Lon_1, OMerc_Lat_1, OMerc_Lon_2, OMerc_Lat_2, OMerc_False_Easting, OMerc_False_Northing, OMerc_Scale_Factor );
00441 }
00442
00443
00444 MSP::CCS::MapProjectionCoordinates* ObliqueMercator::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00445 {
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459 double dlam, B_dlam, cos_B_dlam;
00460 double t, S, T, V, U;
00461 double Q, Q_inv;
00462
00463
00464 double v = 0;
00465 double u = 0;
00466 char errorStatus[256] = "";
00467
00468 double longitude = geodeticCoordinates->longitude();
00469 double latitude = geodeticCoordinates->latitude();
00470
00471 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00472 {
00473 strcat( errorStatus, ErrorMessages::latitude );
00474 }
00475 if ((longitude < -PI) || (longitude > TWO_PI))
00476 {
00477 strcat( errorStatus, ErrorMessages::longitude );
00478 }
00479
00480 if( strlen( errorStatus ) > 0)
00481 throw CoordinateConversionException( errorStatus );
00482
00483 dlam = longitude - OMerc_Origin_Long;
00484
00485 if (fabs(dlam) >= PI_OVER_2)
00486 {
00487 strcat( errorStatus, MSP::CCS::WarningMessages::longitude );
00488 }
00489
00490 if (dlam > PI)
00491 {
00492 dlam -= TWO_PI;
00493 }
00494 if (dlam < -PI)
00495 {
00496 dlam += TWO_PI;
00497 }
00498
00499 if (fabs(fabs(latitude) - PI_OVER_2) > 1.0e-10)
00500 {
00501 t = omercT(latitude, es * sin(latitude), es_OVER_2);
00502 Q = OMerc_E / pow(t, OMerc_B);
00503 Q_inv = 1.0 / Q;
00504 S = (Q - Q_inv) / 2.0;
00505 T = (Q + Q_inv) / 2.0;
00506 B_dlam = OMerc_B * dlam;
00507 V = sin(B_dlam);
00508 U = ((-1.0 * V * cos_gamma) + (S * sin_gamma)) / T;
00509 if (fabs(fabs(U) - 1.0) < 1.0e-10)
00510 {
00511 throw CoordinateConversionException( ErrorMessages::longitude );
00512 }
00513 else
00514 {
00515 v = A_over_B * log((1.0 - U) / (1.0 + U)) / 2.0;
00516 cos_B_dlam = cos(B_dlam);
00517 if (fabs(cos_B_dlam) < 1.0e-10)
00518 u = OMerc_A * B_dlam;
00519 else
00520 u = A_over_B * atan(((S * cos_gamma) + (V * sin_gamma)) / cos_B_dlam);
00521 }
00522 }
00523 else
00524 {
00525 if (latitude > 0.0)
00526 v = A_over_B * log(tan(PI_OVER_4 - (OMerc_gamma / 2.0)));
00527 else
00528 v = A_over_B * log(tan(PI_OVER_4 + (OMerc_gamma / 2.0)));
00529 u = A_over_B * latitude;
00530 }
00531
00532
00533 u = u - OMerc_u;
00534
00535 double easting = OMerc_False_Easting + v * cos_azimuth + u * sin_azimuth;
00536 double northing = OMerc_False_Northing + u * cos_azimuth - v * sin_azimuth;
00537
00538 return new MapProjectionCoordinates( CoordinateType::obliqueMercator, errorStatus, easting, northing );
00539 }
00540
00541
00542 MSP::CCS::GeodeticCoordinates* ObliqueMercator::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00543 {
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557 double dx, dy;
00558
00559
00560 double u, v;
00561 double Q_prime, Q_prime_inv;
00562 double S_prime, T_prime, V_prime, U_prime;
00563 double t;
00564 double es_sin;
00565 double u_B_over_A;
00566 double phi;
00567 double temp_phi = 0.0;
00568 int count = 60;
00569 double longitude, latitude;
00570 char errorStatus[256] = "";
00571
00572 double easting = mapProjectionCoordinates->easting();
00573 double northing = mapProjectionCoordinates->northing();
00574
00575 if ((easting < (OMerc_False_Easting - OMerc_Delta_Easting))
00576 || (easting > (OMerc_False_Easting + OMerc_Delta_Easting)))
00577 {
00578 strcat( errorStatus, ErrorMessages::easting );
00579 }
00580 if ((northing < (OMerc_False_Northing - OMerc_Delta_Northing))
00581 || (northing > (OMerc_False_Northing + OMerc_Delta_Northing)))
00582 {
00583 strcat( errorStatus, ErrorMessages::northing );
00584 }
00585
00586 if( strlen( errorStatus ) > 0)
00587 throw CoordinateConversionException( errorStatus );
00588
00589 dy = northing - OMerc_False_Northing;
00590 dx = easting - OMerc_False_Easting;
00591 v = dx * cos_azimuth - dy * sin_azimuth;
00592 u = dy * cos_azimuth + dx * sin_azimuth;
00593 u = u + OMerc_u;
00594 Q_prime = exp(-1.0 * (v * B_over_A ));
00595 Q_prime_inv = 1.0 / Q_prime;
00596 S_prime = (Q_prime - Q_prime_inv) / 2.0;
00597 T_prime = (Q_prime + Q_prime_inv) / 2.0;
00598 u_B_over_A = u * B_over_A;
00599 V_prime = sin(u_B_over_A);
00600 U_prime = (V_prime * cos_gamma + S_prime * sin_gamma) / T_prime;
00601 if (fabs(fabs(U_prime) - 1.0) < 1.0e-10)
00602 {
00603 if (U_prime > 0)
00604 latitude = PI_OVER_2;
00605 else
00606 latitude = -PI_OVER_2;
00607 longitude = OMerc_Origin_Long;
00608 }
00609 else
00610 {
00611 t = pow(OMerc_E / sqrt((1.0 + U_prime) / (1.0 - U_prime)), 1.0 / OMerc_B);
00612 phi = PI_OVER_2 - 2.0 * atan(t);
00613 while (fabs(phi - temp_phi) > 1.0e-10 && count)
00614 {
00615 temp_phi = phi;
00616 es_sin = es * sin(phi);
00617 phi = PI_OVER_2 - 2.0 * atan(t * pow((1.0 - es_sin) / (1.0 + es_sin), es_OVER_2));
00618 count --;
00619 }
00620
00621 if(!count)
00622 throw CoordinateConversionException( ErrorMessages::northing );
00623
00624 latitude = phi;
00625 longitude = OMerc_Origin_Long - atan2((S_prime * cos_gamma - V_prime * sin_gamma), cos(u_B_over_A)) / OMerc_B;
00626 }
00627
00628 if (fabs(latitude) < 2.0e-7)
00629 latitude = 0.0;
00630 if (latitude > PI_OVER_2)
00631 latitude = PI_OVER_2;
00632 else if (latitude < -PI_OVER_2)
00633 latitude = -PI_OVER_2;
00634
00635 if (longitude > PI)
00636 longitude -= TWO_PI;
00637 if (longitude < -PI)
00638 longitude += TWO_PI;
00639
00640 if (fabs(longitude) < 2.0e-7)
00641 longitude = 0.0;
00642 if (longitude > PI)
00643 longitude = PI;
00644 else if (longitude < -PI)
00645 longitude = -PI;
00646
00647 if (fabs(longitude - OMerc_Origin_Long) >= PI_OVER_2)
00648 {
00649 strcat( errorStatus, MSP::CCS::WarningMessages::longitude );
00650 }
00651
00652 return new GeodeticCoordinates( CoordinateType::geodetic, errorStatus, longitude, latitude );
00653 }
00654
00655
00656 double ObliqueMercator::omercT( double lat, double e_sinlat, double e_over_2 )
00657 {
00658 return (tan(PI_OVER_4 - lat / 2.0)) / (pow((1 - e_sinlat) / (1 + e_sinlat), e_over_2));
00659 }
00660
00661
00662