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 "Polyconic.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 FOURTY_ONE (41.0 * PI / 180);
00118
00119
00120 double polyCoeffTimesSine( double coeff, double x, double latit )
00121 {
00122 return coeff * (sin (x * latit));
00123 }
00124
00125
00126
00127
00128
00129
00130
00131 Polyconic::Polyconic( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double originLatitude, double falseEasting, double falseNorthing ) :
00132 CoordinateSystem(),
00133 es2( 0.0066943799901413800 ),
00134 es4( 4.4814723452405e-005 ),
00135 es6( 3.0000678794350e-007 ),
00136 M0( 0.0 ),
00137 c0( .99832429845280 ),
00138 c1( .0025146070605187 ),
00139 c2( 2.6390465943377e-006 ),
00140 c3( 3.4180460865959e-009 ),
00141 Poly_Origin_Long( 0.0 ),
00142 Poly_Origin_Lat( 0.0 ),
00143 Poly_False_Easting( 0.0 ),
00144 Poly_False_Northing( 0.0 ),
00145 Poly_Max_Easting( 20037509.0 ),
00146 Poly_Max_Northing( 15348215.0 ),
00147 Poly_Min_Easting( -20037509.0 ),
00148 Poly_Min_Northing( -15348215.0 )
00149 {
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168 double j, three_es4;
00169 double lat, sin2lat, sin4lat, sin6lat;
00170 double inv_f = 1 / ellipsoidFlattening;
00171 char errorStatus[500] = "";
00172
00173 if (ellipsoidSemiMajorAxis <= 0.0)
00174 {
00175 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00176 }
00177 if ((inv_f < 250) || (inv_f > 350))
00178 {
00179 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00180 }
00181 if ((originLatitude < -PI_OVER_2) || (originLatitude > PI_OVER_2))
00182 {
00183 strcat( errorStatus, ErrorMessages::originLatitude );
00184 }
00185 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00186 {
00187 strcat( errorStatus, ErrorMessages::centralMeridian );
00188 }
00189
00190 if( strlen( errorStatus ) > 0)
00191 throw CoordinateConversionException( errorStatus );
00192
00193 semiMajorAxis = ellipsoidSemiMajorAxis;
00194 flattening = ellipsoidFlattening;
00195
00196 Poly_Origin_Lat = originLatitude;
00197 if (centralMeridian > PI)
00198 centralMeridian -= TWO_PI;
00199 Poly_Origin_Long = centralMeridian;
00200 Poly_False_Northing = falseNorthing;
00201 Poly_False_Easting = falseEasting;
00202 es2 = 2 * flattening - flattening * flattening;
00203 es4 = es2 * es2;
00204 es6 = es4 * es2;
00205
00206 j = 45.0 * es6 / 1024.0;
00207 three_es4 = 3.0 * es4;
00208 c0 = 1.0 - es2 / 4.0 - three_es4 / 64.0 - 5.0 * es6 / 256.0;
00209 c1 = 3.0 * es2 / 8.0 + three_es4 / 32.0 + j;
00210 c2 = 15.0 * es4 / 256.0 + j;
00211 c3 = 35.0 * es6 / 3072.0;
00212
00213 lat = c0 * Poly_Origin_Lat;
00214 sin2lat = polyCoeffTimesSine(c1, 2.0, Poly_Origin_Lat);
00215 sin4lat = polyCoeffTimesSine(c2, 4.0, Poly_Origin_Lat);
00216 sin6lat = polyCoeffTimesSine(c3, 6.0, Poly_Origin_Lat);
00217 M0 = polyM(lat, sin2lat, sin4lat, sin6lat);
00218
00219 if (Poly_Origin_Long > 0)
00220 {
00221 GeodeticCoordinates gcMax( CoordinateType::geodetic, Poly_Origin_Long - PI, FOURTY_ONE );
00222 MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcMax );
00223 Poly_Max_Northing = tempCoordinates->northing();
00224 delete tempCoordinates;
00225
00226 GeodeticCoordinates gcMin( CoordinateType::geodetic, Poly_Origin_Long - PI, -FOURTY_ONE );
00227 tempCoordinates = convertFromGeodetic( &gcMin );
00228 Poly_Min_Northing = tempCoordinates->northing();
00229 delete tempCoordinates;
00230
00231 Poly_Max_Easting = 19926189.0;
00232 Poly_Min_Easting = -20037509.0;
00233 }
00234 else if (Poly_Origin_Long < 0)
00235 {
00236 GeodeticCoordinates gcMax( CoordinateType::geodetic, Poly_Origin_Long + PI, FOURTY_ONE );
00237 MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcMax );
00238 Poly_Max_Northing = tempCoordinates->northing();
00239 delete tempCoordinates;
00240
00241 GeodeticCoordinates gcMin( CoordinateType::geodetic, Poly_Origin_Long + PI, -FOURTY_ONE );
00242 tempCoordinates = convertFromGeodetic( &gcMin );
00243 Poly_Min_Northing = tempCoordinates->northing();
00244 delete tempCoordinates;
00245
00246 Poly_Max_Easting = 20037509.0;
00247 Poly_Min_Easting = -19926189.0;
00248 }
00249 else
00250 {
00251 GeodeticCoordinates gcMax( CoordinateType::geodetic, PI, FOURTY_ONE );
00252 MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcMax );
00253 Poly_Max_Northing = tempCoordinates->northing();
00254 delete tempCoordinates;
00255
00256 GeodeticCoordinates gcMin( CoordinateType::geodetic, PI, -FOURTY_ONE );
00257 tempCoordinates = convertFromGeodetic( &gcMin ); Poly_Min_Northing = tempCoordinates->northing();
00258 delete tempCoordinates;
00259
00260 Poly_Max_Easting = 20037509.0;
00261 Poly_Min_Easting = -20037509.0;
00262 }
00263
00264 if(Poly_False_Northing)
00265 {
00266 Poly_Max_Northing -= Poly_False_Northing;
00267 Poly_Min_Northing -= Poly_False_Northing;
00268 }
00269 }
00270
00271
00272 Polyconic::Polyconic( const Polyconic &p )
00273 {
00274 semiMajorAxis = p.semiMajorAxis;
00275 flattening = p.flattening;
00276 es2 = p.es2;
00277 es4 = p.es4;
00278 es6 = p.es6;
00279 M0 = p.M0;
00280 c0 = p.c0;
00281 c1 = p.c1;
00282 c2 = p.c2;
00283 c3 = p.c3;
00284 Poly_Origin_Long = p.Poly_Origin_Long;
00285 Poly_Origin_Lat = p.Poly_Origin_Lat;
00286 Poly_False_Easting = p.Poly_False_Easting;
00287 Poly_False_Northing = p.Poly_False_Northing;
00288 Poly_Max_Easting = p.Poly_Max_Easting;
00289 Poly_Max_Northing = p.Poly_Max_Northing;
00290 Poly_Min_Easting = p.Poly_Min_Easting;
00291 Poly_Min_Northing = p.Poly_Min_Northing;
00292 }
00293
00294
00295 Polyconic::~Polyconic()
00296 {
00297 }
00298
00299
00300 Polyconic& Polyconic::operator=( const Polyconic &p )
00301 {
00302 if( this != &p )
00303 {
00304 semiMajorAxis = p.semiMajorAxis;
00305 flattening = p.flattening;
00306 es2 = p.es2;
00307 es4 = p.es4;
00308 es6 = p.es6;
00309 M0 = p.M0;
00310 c0 = p.c0;
00311 c1 = p.c1;
00312 c2 = p.c2;
00313 c3 = p.c3;
00314 Poly_Origin_Long = p.Poly_Origin_Long;
00315 Poly_Origin_Lat = p.Poly_Origin_Lat;
00316 Poly_False_Easting = p.Poly_False_Easting;
00317 Poly_False_Northing = p.Poly_False_Northing;
00318 Poly_Max_Easting = p.Poly_Max_Easting;
00319 Poly_Max_Northing = p.Poly_Max_Northing;
00320 Poly_Min_Easting = p.Poly_Min_Easting;
00321 Poly_Min_Northing = p.Poly_Min_Northing;
00322 }
00323
00324 return *this;
00325 }
00326
00327
00328 MapProjection4Parameters* Polyconic::getParameters() const
00329 {
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346 return new MapProjection4Parameters( CoordinateType::polyconic, Poly_Origin_Long, Poly_Origin_Lat, Poly_False_Easting, Poly_False_Northing );
00347 }
00348
00349
00350 MSP::CCS::MapProjectionCoordinates* Polyconic::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00351 {
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365 double lat, sin2lat, sin4lat, sin6lat;
00366 double dlam;
00367 double NN;
00368 double NN_OVER_tlat;
00369 double MM;
00370 double EE;
00371 double easting, northing;
00372 char errorStatus[256] = "";
00373
00374 double longitude = geodeticCoordinates->longitude();
00375 double latitude = geodeticCoordinates->latitude();
00376 double slat = sin(latitude);
00377
00378 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00379 {
00380 strcat( errorStatus, ErrorMessages::latitude );
00381 }
00382 if ((longitude < -PI) || (longitude > TWO_PI))
00383 {
00384 strcat( errorStatus, ErrorMessages::longitude );
00385 }
00386
00387 if( strlen( errorStatus ) > 0)
00388 throw CoordinateConversionException( errorStatus );
00389
00390 dlam = longitude - Poly_Origin_Long;
00391 if (fabs(dlam) > (PI / 2))
00392 {
00393 strcat( errorStatus, MSP::CCS::WarningMessages::longitude );
00394 }
00395 if (dlam > PI)
00396 {
00397 dlam -= TWO_PI;
00398 }
00399 if (dlam < -PI)
00400 {
00401 dlam += TWO_PI;
00402 }
00403 if (latitude == 0.0)
00404 {
00405 easting = semiMajorAxis * dlam + Poly_False_Easting;
00406 northing = -M0 + Poly_False_Northing;
00407 }
00408 else
00409 {
00410 NN = semiMajorAxis / sqrt(1.0 - es2 * (slat * slat));
00411 NN_OVER_tlat = NN / tan(latitude);
00412 lat = c0 * latitude;
00413 sin2lat = polyCoeffTimesSine(c1, 2.0, latitude);
00414 sin4lat = polyCoeffTimesSine(c2, 4.0, latitude);
00415 sin6lat = polyCoeffTimesSine(c3, 6.0, latitude);
00416 MM = polyM(lat, sin2lat, sin4lat, sin6lat);
00417 EE = dlam * slat;
00418 easting = NN_OVER_tlat * sin(EE) + Poly_False_Easting;
00419 northing = MM - M0 + NN_OVER_tlat * (1.0 - cos(EE)) +
00420 Poly_False_Northing;
00421 }
00422
00423 return new MapProjectionCoordinates( CoordinateType::polyconic, errorStatus, easting, northing );
00424 }
00425
00426
00427 MSP::CCS::GeodeticCoordinates* Polyconic::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00428 {
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442 double dx;
00443 double dy;
00444 double dx_OVER_Poly_a;
00445 double AA;
00446 double BB;
00447 double CC = 0.0;
00448 double PHIn, Delta_PHI = 1.0;
00449 double sin_PHIn;
00450 double PHI, sin2PHI,sin4PHI, sin6PHI;
00451 double Mn, Mn_prime, Ma;
00452 double AA_Ma;
00453 double Ma2_PLUS_BB;
00454 double AA_MINUS_Ma;
00455 double tolerance = 1.0e-12;
00456
00457 double longitude, latitude;
00458 int count = 45000;
00459 char errorStatus[50] = "";
00460
00461 double easting = mapProjectionCoordinates->easting();
00462 double northing = mapProjectionCoordinates->northing();
00463
00464 if ((easting < (Poly_False_Easting + Poly_Min_Easting))
00465 || (easting > (Poly_False_Easting + Poly_Max_Easting)))
00466 {
00467 strcat( errorStatus, ErrorMessages::easting );
00468 }
00469 if ((northing < (Poly_False_Northing + Poly_Min_Northing))
00470 || (northing > (Poly_False_Northing + Poly_Max_Northing)))
00471 {
00472 strcat( errorStatus, ErrorMessages::northing );
00473 }
00474
00475 if( strlen( errorStatus ) > 0)
00476 throw CoordinateConversionException( errorStatus );
00477
00478 dy = northing - Poly_False_Northing;
00479 dx = easting - Poly_False_Easting;
00480 dx_OVER_Poly_a = dx / semiMajorAxis;
00481 if (floatEq(dy,-M0,1))
00482 {
00483 latitude = 0.0;
00484 longitude = dx_OVER_Poly_a + Poly_Origin_Long;
00485 }
00486 else
00487 {
00488 AA = (M0 + dy) / semiMajorAxis;
00489 BB = dx_OVER_Poly_a * dx_OVER_Poly_a + (AA * AA);
00490 PHIn = AA;
00491
00492 while (fabs(Delta_PHI) > tolerance && count)
00493 {
00494 sin_PHIn = sin(PHIn);
00495 CC = sqrt(1.0 - es2 * sin_PHIn * sin_PHIn) * tan(PHIn);
00496 PHI = c0 * PHIn;
00497 sin2PHI = polyCoeffTimesSine(c1, 2.0, PHIn);
00498 sin4PHI = polyCoeffTimesSine(c2, 4.0, PHIn);
00499 sin6PHI = polyCoeffTimesSine(c3, 6.0, PHIn);
00500 Mn = polyM(PHI, sin2PHI, sin4PHI, sin6PHI);
00501 Mn_prime = c0 - 2.0 * c1 * cos(2.0 * PHIn) + 4.0 * c2 * cos(4.0 * PHIn) -
00502 6.0 * c3 * cos(6.0 * PHIn);
00503 Ma = Mn / semiMajorAxis;
00504 AA_Ma = AA * Ma;
00505 Ma2_PLUS_BB = Ma * Ma + BB;
00506 AA_MINUS_Ma = AA - Ma;
00507 Delta_PHI = (AA_Ma * CC + AA_MINUS_Ma - 0.5 * (Ma2_PLUS_BB) * CC) /
00508 (es2 * sin2PHI * (Ma2_PLUS_BB - 2.0 * AA_Ma) /
00509 4.0 * CC + (AA_MINUS_Ma) * (CC * Mn_prime - 2.0 / sin2PHI) - Mn_prime);
00510 PHIn -= Delta_PHI;
00511 count --;
00512 }
00513
00514 if(!count)
00515 throw CoordinateConversionException( ErrorMessages::northing );
00516
00517 latitude = PHIn;
00518
00519 if (latitude > PI_OVER_2)
00520 latitude = PI_OVER_2;
00521 else if (latitude < -PI_OVER_2)
00522 latitude = -PI_OVER_2;
00523
00524 if (floatEq(fabs(latitude),PI_OVER_2,.00001) || (latitude == 0))
00525 longitude = Poly_Origin_Long;
00526
00527 else
00528 {
00529 longitude = (asin(dx_OVER_Poly_a * CC)) / sin(latitude) +
00530 Poly_Origin_Long;
00531 }
00532 }
00533 if (longitude > PI)
00534 longitude -= TWO_PI;
00535 if (longitude < -PI)
00536 longitude += TWO_PI;
00537
00538 if (longitude > PI)
00539 longitude = PI;
00540 else if (longitude < -PI)
00541 longitude = -PI;
00542
00543 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00544 }
00545
00546
00547 double Polyconic::polyM( double c0lat, double c1s2lat, double c2s4lat, double c3s6lat )
00548 {
00549 return semiMajorAxis * (c0lat - c1s2lat + c2s4lat - c3s6lat);
00550 }
00551
00552
00553 double Polyconic::floatEq( double x, double v, double epsilon )
00554 {
00555 return ((v - epsilon) < x) && (x < (v + epsilon));
00556 }
00557
00558
00559
00560