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