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 #include <math.h>
00083 #include "Eckert6.h"
00084 #include "MapProjection3Parameters.h"
00085 #include "MapProjectionCoordinates.h"
00086 #include "GeodeticCoordinates.h"
00087 #include "CoordinateConversionException.h"
00088 #include "ErrorMessages.h"
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101 using namespace MSP::CCS;
00102
00103
00104
00105
00106
00107
00108
00109 const double PI = 3.14159265358979323e0;
00110 const double PI_OVER_2 = ( PI / 2.0);
00111 const double MAX_LAT = ( (PI * 90) / 180.0 );
00112 const double TWO_PI = (2.0 * PI);
00113 const double one_PLUS_PI_OVER_2 = (1.0 + PI / 2.0);
00114
00115
00116
00117
00118
00119
00120
00121 Eckert6::Eckert6( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double falseEasting, double falseNorthing ) :
00122 CoordinateSystem(),
00123 es2( 0.0066943799901413800 ),
00124 es4( 4.4814723452405e-005 ),
00125 es6( 3.0000678794350e-007 ),
00126 Ra_Over_Sqrt_Two_Plus_PI( 2809695.5356062 ),
00127 Inv_Ra_Over_Sqrt_Two_Plus_PI( 3.5591044913137e-007 ),
00128 Eck6_Origin_Long( 0.0 ),
00129 Eck6_False_Easting( 0.0 ),
00130 Eck6_False_Northing( 0.0 ),
00131 Eck6_Delta_Northing( 8826919.0 ),
00132 Eck6_Max_Easting( 17653838.0 )
00133 {
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150 double Ra;
00151 double inv_f = 1 / ellipsoidFlattening;
00152 char errorStatus[500] = "";
00153
00154 if (ellipsoidSemiMajorAxis <= 0.0)
00155 {
00156 strcat( errorStatus, MSP::CCS::ErrorMessages::semiMajorAxis );
00157 }
00158 if ((inv_f < 250) || (inv_f > 350))
00159 {
00160 strcat( errorStatus, MSP::CCS::ErrorMessages::ellipsoidFlattening );
00161 }
00162 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00163 {
00164 strcat( errorStatus, MSP::CCS::ErrorMessages::centralMeridian );
00165 }
00166
00167 if( strlen( errorStatus ) > 0)
00168 throw CoordinateConversionException( errorStatus );
00169
00170 semiMajorAxis = ellipsoidSemiMajorAxis;
00171 flattening = ellipsoidFlattening;
00172
00173 es2 = 2 * flattening - flattening * flattening;
00174 es4 = es2 * es2;
00175 es6 = es4 * es2;
00176
00177 Ra = semiMajorAxis * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 /3024.0);
00178 Ra_Over_Sqrt_Two_Plus_PI = Ra / (sqrt(2.0 + PI));
00179 Inv_Ra_Over_Sqrt_Two_Plus_PI = 1 / Ra_Over_Sqrt_Two_Plus_PI;
00180 if (centralMeridian > PI)
00181 centralMeridian -= TWO_PI;
00182 Eck6_Origin_Long = centralMeridian;
00183 Eck6_False_Easting = falseEasting;
00184 Eck6_False_Northing = falseNorthing;
00185 if (Eck6_Origin_Long > 0)
00186 {
00187 Eck6_Max_Easting = 17555761.0;
00188 Eck6_Min_Easting = -17653839.0;
00189 }
00190 else if (Eck6_Origin_Long < 0)
00191 {
00192 Eck6_Max_Easting = 17653838.0;
00193 Eck6_Min_Easting = -17555761.0;
00194 }
00195 else
00196 {
00197 Eck6_Max_Easting = 17653838.0;
00198 Eck6_Min_Easting = -17653838.0;
00199 }
00200 }
00201
00202
00203 Eckert6::Eckert6( const Eckert6 &e )
00204 {
00205 semiMajorAxis = e.semiMajorAxis;
00206 flattening = e.flattening;
00207 es2 = e.es2;
00208 es4 = e.es4;
00209 es6 = e.es6;
00210 Ra_Over_Sqrt_Two_Plus_PI = e.Ra_Over_Sqrt_Two_Plus_PI;
00211 Inv_Ra_Over_Sqrt_Two_Plus_PI = e.Inv_Ra_Over_Sqrt_Two_Plus_PI;
00212 Eck6_Origin_Long = e.Eck6_Origin_Long;
00213 Eck6_False_Easting = e.Eck6_False_Easting;
00214 Eck6_False_Northing = e.Eck6_False_Northing;
00215 Eck6_Delta_Northing = e.Eck6_Delta_Northing;
00216 Eck6_Max_Easting = e.Eck6_Max_Easting;
00217 }
00218
00219
00220 Eckert6::~Eckert6()
00221 {
00222 }
00223
00224
00225 Eckert6& Eckert6::operator=( const Eckert6 &e )
00226 {
00227 if( this != &e )
00228 {
00229 semiMajorAxis = e.semiMajorAxis;
00230 flattening = e.flattening;
00231 es2 = e.es2;
00232 es4 = e.es4;
00233 es6 = e.es6;
00234 Ra_Over_Sqrt_Two_Plus_PI = e.Ra_Over_Sqrt_Two_Plus_PI;
00235 Inv_Ra_Over_Sqrt_Two_Plus_PI = e.Inv_Ra_Over_Sqrt_Two_Plus_PI;
00236 Eck6_Origin_Long = e.Eck6_Origin_Long;
00237 Eck6_False_Easting = e.Eck6_False_Easting;
00238 Eck6_False_Northing = e.Eck6_False_Northing;
00239 Eck6_Delta_Northing = e.Eck6_Delta_Northing;
00240 Eck6_Max_Easting = e.Eck6_Max_Easting;
00241 }
00242
00243 return *this;
00244 }
00245
00246
00247 MapProjection3Parameters* Eckert6::getParameters() const
00248 {
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263 return new MapProjection3Parameters( CoordinateType::eckert6, Eck6_Origin_Long, Eck6_False_Easting, Eck6_False_Northing );
00264 }
00265
00266
00267 MSP::CCS::MapProjectionCoordinates* Eckert6::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00268 {
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282 double dlam;
00283 double delta_theta = 1.0;
00284 double dt_tolerance = 4.85e-10;
00285
00286 int count = 60;
00287 char errorStatus[50] = "";
00288
00289 double longitude = geodeticCoordinates->longitude();
00290 double latitude = geodeticCoordinates->latitude();
00291 double slat = sin(latitude);
00292 double theta = latitude;
00293
00294 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00295 {
00296 strcat( errorStatus, MSP::CCS::ErrorMessages::latitude );
00297 }
00298 if ((longitude < -PI) || (longitude > TWO_PI))
00299 {
00300 strcat( errorStatus, MSP::CCS::ErrorMessages::longitude );
00301 }
00302
00303 if( strlen( errorStatus ) > 0)
00304 throw CoordinateConversionException( errorStatus );
00305
00306 dlam = longitude - Eck6_Origin_Long;
00307 if (dlam > PI)
00308 {
00309 dlam -= TWO_PI;
00310 }
00311 if (dlam < -PI)
00312 {
00313 dlam += TWO_PI;
00314 }
00315 while (fabs(delta_theta) > dt_tolerance && count)
00316 {
00317 delta_theta = -(theta + sin(theta) - one_PLUS_PI_OVER_2 *
00318 slat) / (1.0 + cos(theta));
00319 theta += delta_theta;
00320 count --;
00321 }
00322
00323 if(!count)
00324 throw CoordinateConversionException( ErrorMessages::northing);
00325
00326 double easting = Ra_Over_Sqrt_Two_Plus_PI * dlam * (1.0 + cos(theta)) +
00327 Eck6_False_Easting;
00328 double northing = 2.0 * Ra_Over_Sqrt_Two_Plus_PI * theta + Eck6_False_Northing;
00329
00330 return new MapProjectionCoordinates( CoordinateType::eckert6, easting, northing );
00331 }
00332
00333
00334 MSP::CCS::GeodeticCoordinates* Eckert6::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00335 {
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349 double dx, dy;
00350 double theta;
00351 double i;
00352 double longitude, latitude;
00353 char errorStatus[50] = "";
00354
00355 double easting = mapProjectionCoordinates->easting();
00356 double northing = mapProjectionCoordinates->northing();
00357
00358 if ((easting < (Eck6_False_Easting + Eck6_Min_Easting))
00359 || (easting > (Eck6_False_Easting + Eck6_Max_Easting)))
00360 {
00361 strcat( errorStatus, MSP::CCS::ErrorMessages::easting );
00362 }
00363 if ((northing < (Eck6_False_Northing - Eck6_Delta_Northing))
00364 || (northing > (Eck6_False_Northing + Eck6_Delta_Northing)))
00365 {
00366 strcat( errorStatus, MSP::CCS::ErrorMessages::northing );
00367 }
00368
00369 if( strlen( errorStatus ) > 0)
00370 throw CoordinateConversionException( errorStatus );
00371
00372 dy = northing - Eck6_False_Northing;
00373 dx = easting - Eck6_False_Easting;
00374 theta = Inv_Ra_Over_Sqrt_Two_Plus_PI * dy / 2.0;
00375 i = (theta + sin(theta)) / one_PLUS_PI_OVER_2;
00376 if (i > 1.0)
00377 latitude = MAX_LAT;
00378 else if (i < -1.0)
00379 latitude = -MAX_LAT;
00380 else
00381 latitude = asin(i);
00382 longitude = Eck6_Origin_Long + Inv_Ra_Over_Sqrt_Two_Plus_PI * dx / (1 + cos(theta));
00383
00384 if (latitude > PI_OVER_2)
00385 latitude = PI_OVER_2;
00386 else if (latitude < -PI_OVER_2)
00387 latitude = -PI_OVER_2;
00388
00389 if (longitude > PI)
00390 longitude -= TWO_PI;
00391 if (longitude < -PI)
00392 longitude += TWO_PI;
00393
00394 if (longitude > PI)
00395 longitude = PI;
00396 else if (longitude < -PI)
00397 longitude = -PI;
00398
00399 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00400 }
00401
00402
00403
00404