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