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 #include <math.h>
00088 #include "AzimuthalEquidistant.h"
00089 #include "MapProjection4Parameters.h"
00090 #include "MapProjectionCoordinates.h"
00091 #include "GeodeticCoordinates.h"
00092 #include "CoordinateConversionException.h"
00093 #include "ErrorMessages.h"
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 using namespace MSP::CCS;
00106
00107
00108
00109
00110
00111
00112
00113 const double PI = 3.14159265358979323e0;
00114 const double PI_OVER_2 = ( PI / 2.0);
00115 const double TWO_PI = ( 2.0 * PI);
00116 const double ONE = (1.0 * PI / 180);
00117
00118
00119
00120
00121
00122
00123
00124 AzimuthalEquidistant::AzimuthalEquidistant( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double originLatitude, double falseEasting, double falseNorthing ) :
00125 CoordinateSystem(),
00126 Ra( 6371007.1810824 ),
00127 Sin_Azeq_Origin_Lat( 0.0 ),
00128 Cos_Azeq_Origin_Lat( 1.0 ),
00129 Azeq_Origin_Long( 0.0 ),
00130 Azeq_Origin_Lat( 0.0 ),
00131 Azeq_False_Easting( 0.0 ),
00132 Azeq_False_Northing( 0.0 ),
00133 abs_Azeq_Origin_Lat( 0.0 ),
00134 Azeq_Delta_Northing( 19903915.0 ),
00135 Azeq_Delta_Easting( 19903915.0 )
00136 {
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 double es2, es4, es6;
00156 double inv_f = 1 / ellipsoidFlattening;
00157 char errorStatus[500] = "";
00158
00159 if (ellipsoidSemiMajorAxis <= 0.0)
00160 {
00161 strcat( errorStatus, MSP::CCS::ErrorMessages::semiMajorAxis );
00162 }
00163 if ((inv_f < 250) || (inv_f > 350))
00164 {
00165 strcat( errorStatus, MSP::CCS::ErrorMessages::ellipsoidFlattening );
00166 }
00167 if ((originLatitude < -PI_OVER_2) || (originLatitude > PI_OVER_2))
00168 {
00169 strcat( errorStatus, MSP::CCS::ErrorMessages::originLatitude );
00170 }
00171 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00172 {
00173 strcat( errorStatus, MSP::CCS::ErrorMessages::centralMeridian );
00174 }
00175
00176 if( strlen( errorStatus ) > 0)
00177 throw CoordinateConversionException( errorStatus );
00178
00179 semiMajorAxis = ellipsoidSemiMajorAxis;
00180 flattening = ellipsoidFlattening;
00181
00182 es2 = 2 * flattening - flattening * flattening;
00183 es4 = es2 * es2;
00184 es6 = es4 * es2;
00185
00186 Ra = semiMajorAxis * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 / 3024.0);
00187 Azeq_Origin_Lat = originLatitude;
00188 Sin_Azeq_Origin_Lat = sin(Azeq_Origin_Lat);
00189 Cos_Azeq_Origin_Lat = cos(Azeq_Origin_Lat);
00190 abs_Azeq_Origin_Lat = fabs(Azeq_Origin_Lat);
00191 if (centralMeridian > PI)
00192 centralMeridian -= TWO_PI;
00193 Azeq_Origin_Long = centralMeridian;
00194 Azeq_False_Northing = falseNorthing;
00195 Azeq_False_Easting = falseEasting;
00196
00197 if (fabs(abs_Azeq_Origin_Lat - PI_OVER_2) < 1.0e-10)
00198 {
00199 Azeq_Delta_Northing = 20015110.0;
00200 Azeq_Delta_Easting = 20015110.0;
00201 }
00202 else if (abs_Azeq_Origin_Lat >= 1.0e-10)
00203 {
00204 if (Azeq_Origin_Long > 0.0)
00205 {
00206 GeodeticCoordinates gcDelta( CoordinateType::geodetic, Azeq_Origin_Long - PI + ONE, -Azeq_Origin_Lat );
00207 MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcDelta );
00208 Azeq_Delta_Easting = tempCoordinates->easting();
00209 delete tempCoordinates;
00210 }
00211 else
00212 {
00213 GeodeticCoordinates gcDelta(CoordinateType::geodetic, Azeq_Origin_Long + PI - ONE, -Azeq_Origin_Lat );
00214 MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcDelta );
00215 Azeq_Delta_Easting = tempCoordinates->easting();
00216 delete tempCoordinates;
00217 }
00218
00219 if(Azeq_False_Easting)
00220 Azeq_Delta_Easting -= Azeq_False_Easting;
00221 if (Azeq_Delta_Easting < 0)
00222 Azeq_Delta_Easting = -Azeq_Delta_Easting;
00223
00224 Azeq_Delta_Northing = 19903915.0;
00225 }
00226 else
00227 {
00228 Azeq_Delta_Northing = 19903915.0;
00229 Azeq_Delta_Easting = 19903915.0;
00230 }
00231 }
00232
00233
00234 AzimuthalEquidistant::AzimuthalEquidistant( const AzimuthalEquidistant &ae )
00235 {
00236 semiMajorAxis = ae.semiMajorAxis;
00237 flattening = ae.flattening;
00238 Ra = ae.Ra;
00239 Sin_Azeq_Origin_Lat = ae.Sin_Azeq_Origin_Lat;
00240 Cos_Azeq_Origin_Lat = ae.Cos_Azeq_Origin_Lat;
00241 Azeq_Origin_Long = ae.Azeq_Origin_Long;
00242 Azeq_Origin_Lat = ae.Azeq_Origin_Lat;
00243 Azeq_False_Easting = ae.Azeq_False_Easting;
00244 Azeq_False_Northing = ae.Azeq_False_Northing;
00245 abs_Azeq_Origin_Lat = ae.abs_Azeq_Origin_Lat;
00246 Azeq_Delta_Northing = ae.Azeq_Delta_Northing;
00247 Azeq_Delta_Easting = ae.Azeq_Delta_Easting;
00248 }
00249
00250
00251 AzimuthalEquidistant::~AzimuthalEquidistant()
00252 {
00253 }
00254
00255
00256 AzimuthalEquidistant& AzimuthalEquidistant::operator=( const AzimuthalEquidistant &ae )
00257 {
00258 if( this != &ae )
00259 {
00260 semiMajorAxis = ae.semiMajorAxis;
00261 flattening = ae.flattening;
00262 Ra = ae.Ra;
00263 Sin_Azeq_Origin_Lat = ae.Sin_Azeq_Origin_Lat;
00264 Cos_Azeq_Origin_Lat = ae.Cos_Azeq_Origin_Lat;
00265 Azeq_Origin_Long = ae.Azeq_Origin_Long;
00266 Azeq_Origin_Lat = ae.Azeq_Origin_Lat;
00267 Azeq_False_Easting = ae.Azeq_False_Easting;
00268 Azeq_False_Northing = ae.Azeq_False_Northing;
00269 abs_Azeq_Origin_Lat = ae.abs_Azeq_Origin_Lat;
00270 Azeq_Delta_Northing = ae.Azeq_Delta_Northing;
00271 Azeq_Delta_Easting = ae.Azeq_Delta_Easting;
00272 }
00273
00274 return *this;
00275 }
00276
00277
00278 MapProjection4Parameters* AzimuthalEquidistant::getParameters() const
00279 {
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296 return new MapProjection4Parameters( CoordinateType::azimuthalEquidistant, Azeq_Origin_Long, Azeq_Origin_Lat, Azeq_False_Easting, Azeq_False_Northing );
00297 }
00298
00299
00300 MSP::CCS::MapProjectionCoordinates* AzimuthalEquidistant::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00301 {
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315 double dlam;
00316 double k_prime;
00317 double c;
00318 double cos_c;
00319 double sin_dlam, cos_dlam;
00320 double Ra_kprime;
00321 double Ra_PI_OVER_2_Lat;
00322 double easting, northing;
00323 char errorStatus[50] = "";
00324
00325 double longitude = geodeticCoordinates->longitude();
00326 double latitude = geodeticCoordinates->latitude();
00327 double slat = sin(latitude);
00328 double clat = cos(latitude);
00329
00330 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00331 {
00332 strcat( errorStatus, MSP::CCS::ErrorMessages::latitude );
00333 }
00334 if ((longitude < -PI) || (longitude > TWO_PI))
00335 {
00336 strcat( errorStatus, MSP::CCS::ErrorMessages::longitude );
00337 }
00338
00339 if( strlen( errorStatus ) > 0)
00340 throw CoordinateConversionException( errorStatus );
00341
00342 dlam = longitude - Azeq_Origin_Long;
00343 if (dlam > PI)
00344 {
00345 dlam -= TWO_PI;
00346 }
00347 if (dlam < -PI)
00348 {
00349 dlam += TWO_PI;
00350 }
00351
00352 sin_dlam = sin(dlam);
00353 cos_dlam = cos(dlam);
00354 if (fabs(abs_Azeq_Origin_Lat - PI_OVER_2) < 1.0e-10)
00355 {
00356 if (Azeq_Origin_Lat >= 0.0)
00357 {
00358 Ra_PI_OVER_2_Lat = Ra * (PI_OVER_2 - latitude);
00359 easting = Ra_PI_OVER_2_Lat * sin_dlam + Azeq_False_Easting;
00360 northing = -1.0 * (Ra_PI_OVER_2_Lat * cos_dlam) + Azeq_False_Northing;
00361 }
00362 else
00363 {
00364 Ra_PI_OVER_2_Lat = Ra * (PI_OVER_2 + latitude);
00365 easting = Ra_PI_OVER_2_Lat * sin_dlam + Azeq_False_Easting;
00366 northing = Ra_PI_OVER_2_Lat * cos_dlam + Azeq_False_Northing;
00367 }
00368 }
00369 else if (abs_Azeq_Origin_Lat <= 1.0e-10)
00370 {
00371 cos_c = clat * cos_dlam;
00372 if (fabs(fabs(cos_c) - 1.0) < 1.0e-14)
00373 {
00374 if (cos_c >= 0.0)
00375 {
00376 easting = Azeq_False_Easting;
00377 northing = Azeq_False_Northing;
00378 }
00379 else
00380 {
00381
00382 throw CoordinateConversionException( MSP::CCS::ErrorMessages::projection );
00383 }
00384 }
00385 else
00386 {
00387 c = acos(cos_c);
00388 k_prime = c / sin(c);
00389 Ra_kprime = Ra * k_prime;
00390 easting = Ra_kprime * clat * sin_dlam + Azeq_False_Easting;
00391 northing = Ra_kprime * slat + Azeq_False_Northing;
00392 }
00393 }
00394 else
00395 {
00396 cos_c = (Sin_Azeq_Origin_Lat * slat) + (Cos_Azeq_Origin_Lat * clat * cos_dlam);
00397 if (fabs(fabs(cos_c) - 1.0) < 1.0e-14)
00398 {
00399 if (cos_c >= 0.0)
00400 {
00401 easting = Azeq_False_Easting;
00402 northing = Azeq_False_Northing;
00403 }
00404 else
00405 {
00406
00407 throw CoordinateConversionException( MSP::CCS::ErrorMessages::projection );
00408 }
00409 }
00410 else
00411 {
00412 c = acos(cos_c);
00413 k_prime = c / sin(c);
00414 Ra_kprime = Ra * k_prime;
00415 easting = Ra_kprime * clat * sin_dlam + Azeq_False_Easting;
00416 northing = Ra_kprime * (Cos_Azeq_Origin_Lat * slat - Sin_Azeq_Origin_Lat * clat * cos_dlam) + Azeq_False_Northing;
00417 }
00418 }
00419
00420 return new MapProjectionCoordinates( CoordinateType::azimuthalEquidistant, easting, northing );
00421 }
00422
00423
00424 MSP::CCS::GeodeticCoordinates* AzimuthalEquidistant::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00425 {
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439 double dx, dy;
00440 double rho;
00441 double c;
00442 double sin_c, cos_c, dy_sinc;
00443 double longitude, latitude;
00444 char errorStatus[50] = "";
00445
00446 double easting = mapProjectionCoordinates->easting();
00447 double northing = mapProjectionCoordinates->northing();
00448
00449 if ((easting < (Azeq_False_Easting - Azeq_Delta_Easting))
00450 || (easting > (Azeq_False_Easting + Azeq_Delta_Easting)))
00451 {
00452 strcat( errorStatus, MSP::CCS::ErrorMessages::easting );
00453 }
00454 if ((northing < (Azeq_False_Northing - Azeq_Delta_Northing))
00455 || (northing > (Azeq_False_Northing + Azeq_Delta_Northing)))
00456 {
00457 strcat( errorStatus, MSP::CCS::ErrorMessages::northing );
00458 }
00459
00460 if( strlen( errorStatus ) > 0)
00461 throw CoordinateConversionException( errorStatus );
00462
00463 dy = northing - Azeq_False_Northing;
00464 dx = easting - Azeq_False_Easting;
00465 rho = sqrt(dx * dx + dy * dy);
00466 if (fabs(rho) <= 1.0e-10)
00467 {
00468 latitude = Azeq_Origin_Lat;
00469 longitude = Azeq_Origin_Long;
00470 }
00471 else
00472 {
00473 c = rho / Ra;
00474 sin_c = sin(c);
00475 cos_c = cos(c);
00476 dy_sinc = dy * sin_c;
00477 latitude = asin((cos_c * Sin_Azeq_Origin_Lat) + ((dy_sinc * Cos_Azeq_Origin_Lat) / rho));
00478 if (fabs(abs_Azeq_Origin_Lat - PI_OVER_2) < 1.0e-10)
00479 {
00480 if (Azeq_Origin_Lat >= 0.0)
00481 longitude = Azeq_Origin_Long + atan2(dx, -dy);
00482 else
00483 longitude = Azeq_Origin_Long + atan2(dx, dy);
00484 }
00485 else
00486 longitude = Azeq_Origin_Long + atan2((dx * sin_c), ((rho * Cos_Azeq_Origin_Lat * cos_c) - (dy_sinc * Sin_Azeq_Origin_Lat)));
00487 }
00488
00489 if (latitude > PI_OVER_2)
00490 latitude = PI_OVER_2;
00491 else if (latitude < -PI_OVER_2)
00492 latitude = -PI_OVER_2;
00493
00494 if (longitude > PI)
00495 longitude -= TWO_PI;
00496 if (longitude < -PI)
00497 longitude += TWO_PI;
00498
00499 if (longitude > PI)
00500 longitude = PI;
00501 else if (longitude < -PI)
00502 longitude = -PI;
00503
00504 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00505 }
00506
00507
00508
00509