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