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