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 #include <math.h>
00094 #include "AlbersEqualAreaConic.h"
00095 #include "MapProjection6Parameters.h"
00096 #include "MapProjectionCoordinates.h"
00097 #include "GeodeticCoordinates.h"
00098 #include "CoordinateConversionException.h"
00099 #include "ErrorMessages.h"
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111 using namespace MSP::CCS;
00112
00113
00114
00115
00116
00117
00118
00119 const double PI = 3.14159265358979323e0;
00120 const double PI_OVER_2 = ( PI / 2.0);
00121 const double TWO_PI = (2.0 * PI);
00122
00123
00124 double oneMinusSqr( double x )
00125 {
00126 return 1.0 - x * x;
00127 }
00128
00129
00130 double albersM( double clat, double oneminussqressin )
00131 {
00132 return clat / sqrt(oneminussqressin);
00133 }
00134
00135
00136
00137
00138
00139
00140
00141 AlbersEqualAreaConic::AlbersEqualAreaConic( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double originLatitude, double standardParallel1, double standardParallel2, double falseEasting, double falseNorthing ) :
00142 CoordinateSystem(),
00143 es( 0.08181919084262188000 ),
00144 es2( 0.0066943799901413800 ),
00145 C( 1.4896626908850 ),
00146 rho0( 6388749.3391064 ),
00147 n( .70443998701755 ),
00148 Albers_a_OVER_n( 9054194.9882824 ),
00149 one_MINUS_es2( .99330562000986 ),
00150 two_es( .16363838168524 ),
00151 Albers_Origin_Long( 0.0 ),
00152 Albers_Origin_Lat( (45 * PI / 180) ),
00153 Albers_Std_Parallel_1( 40 * PI / 180 ),
00154 Albers_Std_Parallel_2( 50 * PI / 180 ),
00155 Albers_False_Easting( 0.0 ),
00156 Albers_False_Northing( 0.0 ),
00157 Albers_Delta_Northing( 40000000 ),
00158 Albers_Delta_Easting( 40000000 )
00159 {
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180 double sin_lat, sin_lat_1, cos_lat;
00181 double m1, m2, SQRm1;
00182 double q0, q1, q2;
00183 double es_sin, one_MINUS_SQRes_sin;
00184 double nq0;
00185 double inv_f = 1 / ellipsoidFlattening;
00186 char errorStatus[500] = "";
00187
00188 if (ellipsoidSemiMajorAxis <= 0.0)
00189 {
00190 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00191 }
00192 if ((inv_f < 250) || (inv_f > 350))
00193 {
00194 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00195 }
00196 if ((originLatitude < -PI_OVER_2) || (originLatitude > PI_OVER_2))
00197 {
00198 strcat( errorStatus, ErrorMessages::originLatitude );
00199 }
00200 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00201 {
00202 strcat( errorStatus, ErrorMessages::centralMeridian );
00203 }
00204 if ((standardParallel1 < -PI_OVER_2) || (standardParallel1 > PI_OVER_2))
00205 {
00206 strcat( errorStatus, ErrorMessages::standardParallel1 );
00207 }
00208 if ((standardParallel2 < -PI_OVER_2) || (standardParallel2 > PI_OVER_2))
00209 {
00210 strcat( errorStatus, ErrorMessages::standardParallel2 );
00211 }
00212 if ((standardParallel1 == 0.0) && (standardParallel2 == 0.0))
00213 {
00214 strcat( errorStatus, ErrorMessages::standardParallel1_2 );
00215 }
00216 if (standardParallel1 == -standardParallel2)
00217 {
00218 strcat( errorStatus, ErrorMessages::standardParallelHemisphere );
00219 }
00220
00221 if( strlen( errorStatus ) > 0)
00222 throw CoordinateConversionException( errorStatus );
00223
00224 semiMajorAxis = ellipsoidSemiMajorAxis;
00225 flattening = ellipsoidFlattening;
00226
00227 Albers_Origin_Lat = originLatitude;
00228 Albers_Std_Parallel_1 = standardParallel1;
00229 Albers_Std_Parallel_2 = standardParallel2;
00230 if (centralMeridian > PI)
00231 centralMeridian -= TWO_PI;
00232 Albers_Origin_Long = centralMeridian;
00233 Albers_False_Easting = falseEasting;
00234 Albers_False_Northing = falseNorthing;
00235
00236 es2 = 2 * flattening - flattening * flattening;
00237 es = sqrt(es2);
00238 one_MINUS_es2 = 1 - es2;
00239 two_es = 2 * es;
00240
00241 sin_lat = sin(Albers_Origin_Lat);
00242 es_sin = esSine(sin_lat);
00243 one_MINUS_SQRes_sin = oneMinusSqr( es_sin );
00244 q0 = albersQ( sin_lat, one_MINUS_SQRes_sin, es_sin );
00245
00246 sin_lat_1 = sin(Albers_Std_Parallel_1);
00247 cos_lat = cos(Albers_Std_Parallel_1);
00248 es_sin = esSine(sin_lat_1);
00249 one_MINUS_SQRes_sin = oneMinusSqr( es_sin );
00250 m1 = albersM( cos_lat, one_MINUS_SQRes_sin );
00251 q1 = albersQ( sin_lat_1, one_MINUS_SQRes_sin, es_sin );
00252
00253 SQRm1 = m1 * m1;
00254 if (fabs(Albers_Std_Parallel_1 - Albers_Std_Parallel_2) > 1.0e-10)
00255 {
00256 sin_lat = sin(Albers_Std_Parallel_2);
00257 cos_lat = cos(Albers_Std_Parallel_2);
00258 es_sin = esSine(sin_lat);
00259 one_MINUS_SQRes_sin = oneMinusSqr( es_sin );
00260 m2 = albersM( cos_lat, one_MINUS_SQRes_sin );
00261 q2 = albersQ( sin_lat, one_MINUS_SQRes_sin, es_sin );
00262 n = (SQRm1 - m2 * m2) / (q2 - q1);
00263 }
00264 else
00265 n = sin_lat_1;
00266
00267 C = SQRm1 + n * q1;
00268 Albers_a_OVER_n = semiMajorAxis / n;
00269 nq0 = n * q0;
00270 if (C < nq0)
00271 rho0 = 0;
00272 else
00273 rho0 = Albers_a_OVER_n * sqrt(C - nq0);
00274 }
00275
00276
00277 AlbersEqualAreaConic::AlbersEqualAreaConic( const AlbersEqualAreaConic &aeac )
00278 {
00279 semiMajorAxis = aeac.semiMajorAxis;
00280 flattening = aeac.flattening;
00281 es = aeac.es;
00282 es2 = aeac.es2;
00283 C = aeac.C;
00284 rho0 = aeac.rho0;
00285 n = aeac.n;
00286 Albers_a_OVER_n = aeac.Albers_a_OVER_n;
00287 one_MINUS_es2 = aeac.one_MINUS_es2;
00288 two_es = aeac.two_es;
00289 Albers_Origin_Long = aeac.Albers_Origin_Long;
00290 Albers_Origin_Lat = aeac.Albers_Origin_Lat;
00291 Albers_Std_Parallel_1 = aeac.Albers_Std_Parallel_1;
00292 Albers_Std_Parallel_2 = aeac.Albers_Std_Parallel_2;
00293 Albers_False_Easting = aeac.Albers_False_Easting;
00294 Albers_False_Northing = aeac.Albers_False_Northing;
00295 Albers_Delta_Northing = aeac.Albers_Delta_Northing;
00296 Albers_Delta_Easting = aeac.Albers_Delta_Easting;
00297 }
00298
00299
00300 AlbersEqualAreaConic::~AlbersEqualAreaConic()
00301 {
00302 }
00303
00304
00305 AlbersEqualAreaConic& AlbersEqualAreaConic::operator=( const AlbersEqualAreaConic &aeac )
00306 {
00307 if( this != &aeac )
00308 {
00309 semiMajorAxis = aeac.semiMajorAxis;
00310 flattening = aeac.flattening;
00311 es = aeac.es;
00312 es2 = aeac.es2;
00313 C = aeac.C;
00314 rho0 = aeac.rho0;
00315 n = aeac.n;
00316 Albers_a_OVER_n = aeac.Albers_a_OVER_n;
00317 one_MINUS_es2 = aeac.one_MINUS_es2;
00318 two_es = aeac.two_es;
00319 Albers_Origin_Long = aeac.Albers_Origin_Long;
00320 Albers_Origin_Lat = aeac.Albers_Origin_Lat;
00321 Albers_Std_Parallel_1 = aeac.Albers_Std_Parallel_1;
00322 Albers_Std_Parallel_2 = aeac.Albers_Std_Parallel_2;
00323 Albers_False_Easting = aeac.Albers_False_Easting;
00324 Albers_False_Northing = aeac.Albers_False_Northing;
00325 Albers_Delta_Northing = aeac.Albers_Delta_Northing;
00326 Albers_Delta_Easting = aeac.Albers_Delta_Easting;
00327 }
00328
00329 return *this;
00330 }
00331
00332
00333 MapProjection6Parameters* AlbersEqualAreaConic::getParameters() const
00334 {
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353 return new MapProjection6Parameters( CoordinateType::albersEqualAreaConic, Albers_Origin_Long, Albers_Origin_Lat, Albers_Std_Parallel_1, Albers_Std_Parallel_2, Albers_False_Easting, Albers_False_Northing );
00354 }
00355
00356
00357 MSP::CCS::MapProjectionCoordinates* AlbersEqualAreaConic::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00358 {
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372 double dlam;
00373 double sin_lat, cos_lat;
00374 double es_sin, one_MINUS_SQRes_sin;
00375 double q;
00376 double rho;
00377 double theta;
00378 double nq;
00379 char errorStatus[50] = "";
00380
00381 double longitude = geodeticCoordinates->longitude();
00382 double latitude = geodeticCoordinates->latitude();
00383
00384 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00385 {
00386 strcat( errorStatus, ErrorMessages::latitude );
00387 }
00388 if ((longitude < -PI) || (longitude > TWO_PI))
00389 {
00390 strcat( errorStatus, ErrorMessages::longitude );
00391 }
00392
00393 if( strlen( errorStatus ) > 0)
00394 throw CoordinateConversionException( errorStatus );
00395
00396 dlam = longitude - Albers_Origin_Long;
00397 if (dlam > PI)
00398 {
00399 dlam -= TWO_PI;
00400 }
00401 if (dlam < -PI)
00402 {
00403 dlam += TWO_PI;
00404 }
00405 sin_lat = sin(latitude);
00406 cos_lat = cos(latitude);
00407 es_sin = esSine( sin_lat );
00408 one_MINUS_SQRes_sin = oneMinusSqr( es_sin );
00409 q = albersQ( sin_lat, one_MINUS_SQRes_sin, es_sin );
00410 nq = n * q;
00411 if (C < nq)
00412 rho = 0;
00413 else
00414 rho = Albers_a_OVER_n * sqrt(C - nq);
00415
00416
00417 theta = n * dlam;
00418 double easting = rho * sin(theta) + Albers_False_Easting;
00419 double northing = rho0 - rho * cos(theta) + Albers_False_Northing;
00420
00421 return new MapProjectionCoordinates( CoordinateType::albersEqualAreaConic, easting, northing );
00422 }
00423
00424
00425 MSP::CCS::GeodeticCoordinates* AlbersEqualAreaConic::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00426 {
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440 double dy, dx;
00441 double rho0_MINUS_dy;
00442 double q, qconst, q_OVER_2;
00443 double rho, rho_n;
00444 double PHI, Delta_PHI = 1.0;
00445 double sin_phi;
00446 double es_sin, one_MINUS_SQRes_sin;
00447 double theta = 0.0;
00448 int count = 60;
00449 double longitude, latitude;
00450 double tolerance = 4.85e-10;
00451
00452 char errorStatus[50] = "";
00453
00454 double easting = mapProjectionCoordinates->easting();
00455 double northing = mapProjectionCoordinates->northing();
00456
00457 if ((easting < (Albers_False_Easting - Albers_Delta_Easting))
00458 || (easting > Albers_False_Easting + Albers_Delta_Easting))
00459 {
00460 strcat( errorStatus, ErrorMessages::easting );
00461 }
00462 if ((northing < (Albers_False_Northing - Albers_Delta_Northing))
00463 || (northing > Albers_False_Northing + Albers_Delta_Northing))
00464 {
00465 strcat( errorStatus, ErrorMessages::northing );
00466 }
00467
00468 if( strlen( errorStatus ) > 0)
00469 throw CoordinateConversionException( errorStatus );
00470
00471 dy = northing - Albers_False_Northing;
00472 dx = easting - Albers_False_Easting;
00473 rho0_MINUS_dy = rho0 - dy;
00474 rho = sqrt(dx * dx + rho0_MINUS_dy * rho0_MINUS_dy);
00475
00476 if (n < 0)
00477 {
00478 rho *= -1.0;
00479 dy *= -1.0;
00480 dx *= -1.0;
00481 rho0_MINUS_dy *= -1.0;
00482 }
00483
00484 if (rho != 0.0)
00485 theta = atan2(dx, rho0_MINUS_dy);
00486 rho_n = rho * n;
00487 q = (C - (rho_n * rho_n) / (semiMajorAxis * semiMajorAxis)) / n;
00488 qconst = 1 - ((one_MINUS_es2) / (two_es)) * log((1.0 - es) / (1.0 + es));
00489 if (fabs(fabs(qconst) - fabs(q)) > 1.0e-6)
00490 {
00491 q_OVER_2 = q / 2.0;
00492 if (q_OVER_2 > 1.0)
00493 latitude = PI_OVER_2;
00494 else if (q_OVER_2 < -1.0)
00495 latitude = -PI_OVER_2;
00496 else
00497 {
00498 PHI = asin(q_OVER_2);
00499 if (es < 1.0e-10)
00500 latitude = PHI;
00501 else
00502 {
00503 while ((fabs(Delta_PHI) > tolerance) && count)
00504 {
00505 sin_phi = sin(PHI);
00506 es_sin = esSine( sin_phi );
00507 one_MINUS_SQRes_sin = oneMinusSqr( es_sin );
00508 Delta_PHI = (one_MINUS_SQRes_sin * one_MINUS_SQRes_sin) / (2.0 * cos(PHI)) *
00509 (q / (one_MINUS_es2) - sin_phi / one_MINUS_SQRes_sin +
00510 (log((1.0 - es_sin) / (1.0 + es_sin)) / (two_es)));
00511 PHI += Delta_PHI;
00512 count --;
00513 }
00514
00515 if(!count)
00516 throw CoordinateConversionException( ErrorMessages::northing );
00517
00518 latitude = PHI;
00519 }
00520
00521 if (latitude > PI_OVER_2)
00522 latitude = PI_OVER_2;
00523 else if (latitude < -PI_OVER_2)
00524 latitude = -PI_OVER_2;
00525
00526 }
00527 }
00528 else
00529 {
00530 if (q >= 0.0)
00531 latitude = PI_OVER_2;
00532 else
00533 latitude = -PI_OVER_2;
00534 }
00535 longitude = Albers_Origin_Long + theta / n;
00536
00537 if (longitude > PI)
00538 longitude -= TWO_PI;
00539 if (longitude < -PI)
00540 longitude += TWO_PI;
00541
00542 if (longitude > PI)
00543 longitude = PI;
00544 else if (longitude < -PI)
00545 longitude = -PI;
00546
00547 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00548 }
00549
00550
00551 double AlbersEqualAreaConic::esSine( double sinlat )
00552 {
00553 return es * sinlat;
00554 }
00555
00556
00557 double AlbersEqualAreaConic::albersQ( double slat, double oneminussqressin, double essin )
00558 {
00559 return (one_MINUS_es2)*(slat / (oneminussqressin) -
00560 (1 / (two_es)) *log((1 - essin) / (1 + essin)));
00561 }
00562
00563
00564
00565