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 "Eckert4.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 using namespace MSP::CCS;
00101
00102
00103
00104
00105
00106
00107
00108 const double PI = 3.14159265358979323e0;
00109 const double PI_OVER_2 = ( PI / 2.0);
00110 const double TWO_PI = (2.0 * PI) ;
00111 const double two_PLUS_PI_OVER_2 = (2.0 + PI / 2.0);
00112
00113
00114 double calculateNum( double theta, double sinTheta, double cosTheta )
00115 {
00116 return theta + sinTheta * cosTheta + 2.0 * sinTheta;
00117 }
00118
00119
00120
00121
00122
00123
00124
00125 Eckert4::Eckert4( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double falseEasting, double falseNorthing ) :
00126 CoordinateSystem(),
00127 es2( 0.0066943799901413800 ),
00128 es4( 4.4814723452405e-005 ),
00129 es6( 3.0000678794350e-007 ),
00130 Ra0( 2690082.6043273 ),
00131 Ra1( 8451143.5741087 ),
00132 Eck4_Origin_Long( 0.0 ),
00133 Eck4_False_Easting( 0.0 ),
00134 Eck4_False_Northing( 0.0 ),
00135 Eck4_Delta_Northing( 8451144.0 ),
00136 Eck4_Max_Easting( 16902288.0 ),
00137 Eck4_Min_Easting( -16902288.0 )
00138 {
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 double Ra;
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 ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00168 {
00169 strcat( errorStatus, MSP::CCS::ErrorMessages::centralMeridian );
00170 }
00171
00172 if( strlen( errorStatus ) > 0)
00173 throw CoordinateConversionException( errorStatus );
00174
00175 semiMajorAxis = ellipsoidSemiMajorAxis;
00176 flattening = ellipsoidFlattening;
00177
00178 es2 = 2 * flattening - flattening * flattening;
00179 es4 = es2 * es2;
00180 es6 = es4 * es2;
00181
00182 Ra = semiMajorAxis * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 / 3024.0);
00183 Ra0 = 0.4222382 * Ra;
00184 Ra1 = 1.3265004 * Ra;
00185 if (centralMeridian > PI)
00186 centralMeridian -= TWO_PI;
00187 Eck4_Origin_Long = centralMeridian;
00188 Eck4_False_Easting = falseEasting;
00189 Eck4_False_Northing = falseNorthing;
00190 if (Eck4_Origin_Long > 0)
00191 {
00192 Eck4_Max_Easting = 16808386.0;
00193 Eck4_Min_Easting = -16902288.0;
00194 }
00195 else if (Eck4_Origin_Long < 0)
00196 {
00197 Eck4_Max_Easting = 16902288.0;
00198 Eck4_Min_Easting = -16808386.0;
00199 }
00200 else
00201 {
00202 Eck4_Max_Easting = 16902288.0;
00203 Eck4_Min_Easting = -16902288.0;
00204 }
00205 }
00206
00207
00208 Eckert4::Eckert4( const Eckert4 &e )
00209 {
00210 semiMajorAxis = e.semiMajorAxis;
00211 flattening = e.flattening;
00212 es2 = e.es2;
00213 es4 = e.es4;
00214 es6 = e.es6;
00215 Ra0 = e.Ra0;
00216 Ra1 = e.Ra1;
00217 Eck4_Origin_Long = e.Eck4_Origin_Long;
00218 Eck4_False_Easting = e.Eck4_False_Easting;
00219 Eck4_False_Northing = e.Eck4_False_Northing;
00220 Eck4_Delta_Northing = e.Eck4_Delta_Northing;
00221 Eck4_Max_Easting = e.Eck4_Max_Easting;
00222 Eck4_Min_Easting = e.Eck4_Min_Easting;
00223 }
00224
00225
00226 Eckert4::~Eckert4()
00227 {
00228 }
00229
00230
00231 Eckert4& Eckert4::operator=( const Eckert4 &e )
00232 {
00233 if( this != &e )
00234 {
00235 semiMajorAxis = e.semiMajorAxis;
00236 flattening = e.flattening;
00237 es2 = e.es2;
00238 es4 = e.es4;
00239 es6 = e.es6;
00240 Ra0 = e.Ra0;
00241 Ra1 = e.Ra1;
00242 Eck4_Origin_Long = e.Eck4_Origin_Long;
00243 Eck4_False_Easting = e.Eck4_False_Easting;
00244 Eck4_False_Northing = e.Eck4_False_Northing;
00245 Eck4_Delta_Northing = e.Eck4_Delta_Northing;
00246 Eck4_Max_Easting = e.Eck4_Max_Easting;
00247 Eck4_Min_Easting = e.Eck4_Min_Easting;
00248 }
00249
00250 return *this;
00251 }
00252
00253
00254 MapProjection3Parameters* Eckert4::getParameters() const
00255 {
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270 return new MapProjection3Parameters( CoordinateType::eckert4, Eck4_Origin_Long, Eck4_False_Easting, Eck4_False_Northing );
00271 }
00272
00273
00274 MSP::CCS::MapProjectionCoordinates* Eckert4::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00275 {
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290 double sin_theta, cos_theta;
00291 double num;
00292 double dlam;
00293 double delta_theta = 1.0;
00294 double dt_tolerance = 4.85e-10;
00295
00296 int count = 200;
00297 char errorStatus[50] = "";
00298
00299 double longitude = geodeticCoordinates->longitude();
00300 double latitude = geodeticCoordinates->latitude();
00301 double slat = sin(latitude);
00302 double theta = latitude / 2.0;
00303
00304 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00305 {
00306 strcat( errorStatus, MSP::CCS::ErrorMessages::latitude );
00307 }
00308 if ((longitude < -PI) || (longitude > TWO_PI))
00309 {
00310 strcat( errorStatus, MSP::CCS::ErrorMessages::longitude );
00311 }
00312
00313 if( strlen( errorStatus ) > 0)
00314 throw CoordinateConversionException( errorStatus );
00315
00316 dlam = longitude - Eck4_Origin_Long;
00317 if (dlam > PI)
00318 {
00319 dlam -= TWO_PI;
00320 }
00321 if (dlam < -PI)
00322 {
00323 dlam += TWO_PI;
00324 }
00325 while (fabs(delta_theta) > dt_tolerance && count)
00326 {
00327 sin_theta = sin(theta);
00328 cos_theta = cos(theta);
00329 num = calculateNum(theta, sin_theta, cos_theta);
00330 delta_theta = -(num - two_PLUS_PI_OVER_2 * slat) /
00331 (2.0 * cos_theta * (1.0 + cos_theta));
00332 theta += delta_theta;
00333 count --;
00334 }
00335
00336 if(!count)
00337 throw CoordinateConversionException( ErrorMessages::northing );
00338
00339 double easting = Ra0 * dlam * (1.0 + cos(theta)) + Eck4_False_Easting;
00340 double northing = Ra1 * sin(theta) + Eck4_False_Northing;
00341
00342 return new MapProjectionCoordinates( CoordinateType::eckert4, easting, northing );
00343 }
00344
00345
00346 MSP::CCS::GeodeticCoordinates* Eckert4::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00347 {
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362 double theta;
00363 double sin_theta, cos_theta;
00364 double num;
00365 double dx, dy;
00366 double i;
00367 char errorStatus[50] = "";
00368
00369 double easting = mapProjectionCoordinates->easting();
00370 double northing = mapProjectionCoordinates->northing();
00371
00372 if ((easting < (Eck4_False_Easting + Eck4_Min_Easting))
00373 || (easting > (Eck4_False_Easting + Eck4_Max_Easting)))
00374 {
00375 strcat( errorStatus, MSP::CCS::ErrorMessages::easting );
00376 }
00377 if ((northing < (Eck4_False_Northing - Eck4_Delta_Northing))
00378 || (northing > (Eck4_False_Northing + Eck4_Delta_Northing)))
00379 {
00380 strcat( errorStatus, MSP::CCS::ErrorMessages::northing );
00381 }
00382
00383 if( strlen( errorStatus ) > 0)
00384 throw CoordinateConversionException( errorStatus );
00385
00386 dy = northing - Eck4_False_Northing;
00387 dx = easting - Eck4_False_Easting;
00388 i = dy / Ra1;
00389 if (i > 1.0)
00390 i = 1.0;
00391 else if (i < -1.0)
00392 i = -1.0;
00393
00394 theta = asin(i);
00395 sin_theta = sin(theta);
00396 cos_theta = cos(theta);
00397 num = calculateNum(theta, sin_theta, cos_theta);
00398
00399 double latitude = asin(num / two_PLUS_PI_OVER_2);
00400 double longitude = Eck4_Origin_Long + dx / (Ra0 * (1 + cos_theta));
00401
00402 if (latitude > PI_OVER_2)
00403 latitude = PI_OVER_2;
00404 else if (latitude < -PI_OVER_2)
00405 latitude = -PI_OVER_2;
00406
00407 if (longitude > PI)
00408 longitude -= TWO_PI;
00409 if (longitude < -PI)
00410 longitude += TWO_PI;
00411
00412 if (longitude > PI)
00413 longitude = PI;
00414 else if (longitude < -PI)
00415 longitude = -PI;
00416
00417 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00418 }
00419
00420
00421
00422