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 #include <math.h>
00086 #include "CylindricalEqualArea.h"
00087 #include "MapProjection4Parameters.h"
00088 #include "MapProjectionCoordinates.h"
00089 #include "GeodeticCoordinates.h"
00090 #include "CoordinateConversionException.h"
00091 #include "ErrorMessages.h"
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103 using namespace MSP::CCS;
00104
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 const double ONE = (1.0 * PI / 180);
00115
00116
00117 double cyeqCoeffTimesSine( double coeff, double c, double Beta )
00118 {
00119 return coeff * sin(c * Beta);
00120 }
00121
00122
00123
00124
00125
00126
00127
00128 CylindricalEqualArea::CylindricalEqualArea( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double originLatitude, double falseEasting, double falseNorthing ) :
00129 CoordinateSystem(),
00130 es( .081819190842622 ),
00131 es2( 0.0066943799901413800 ),
00132 es4( 4.4814723452405e-005 ),
00133 es6( 3.0000678794350e-007 ),
00134 k0( 1.0 ),
00135 Cyeq_a_k0( 6378137.0 ),
00136 two_k0( 2.0 ),
00137 c0( .0022392088624809 ),
00138 c1( 2.8830839728915e-006 ),
00139 c2( 5.0331826636906e-009 ),
00140 Cyeq_Origin_Long( 0.0 ),
00141 Cyeq_Origin_Lat( 0.0 ),
00142 Cyeq_False_Easting( 0.0 ),
00143 Cyeq_False_Northing( 0.0 ),
00144 Cyeq_Max_Easting( 20037509.0 ),
00145 Cyeq_Min_Easting( -20037509.0 ),
00146 Cyeq_Delta_Northing( 6363886.0 )
00147 {
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166 double Sin_Cyeq_Origin_Lat;
00167 double inv_f = 1 / ellipsoidFlattening;
00168 char errorStatus[500] = "";
00169
00170 if (ellipsoidSemiMajorAxis <= 0.0)
00171 {
00172 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00173 }
00174 if ((inv_f < 250) || (inv_f > 350))
00175 {
00176 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00177 }
00178 if ((originLatitude < -PI_OVER_2) || (originLatitude > PI_OVER_2))
00179 {
00180 strcat( errorStatus, ErrorMessages::originLatitude );
00181 }
00182 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00183 {
00184 strcat( errorStatus, ErrorMessages::centralMeridian );
00185 }
00186
00187 if( strlen( errorStatus ) > 0)
00188 throw CoordinateConversionException( errorStatus );
00189
00190 semiMajorAxis = ellipsoidSemiMajorAxis;
00191 flattening = ellipsoidFlattening;
00192
00193 Cyeq_Origin_Lat = originLatitude;
00194 if (centralMeridian > PI)
00195 centralMeridian -= TWO_PI;
00196 Cyeq_Origin_Long = centralMeridian;
00197 Cyeq_False_Northing = falseNorthing;
00198 Cyeq_False_Easting = falseEasting;
00199 es2 = 2 * flattening - flattening * flattening;
00200 es4 = es2 * es2;
00201 es6 = es4 * es2;
00202 es = sqrt(es2);
00203 c0 = es2 / 3.0 + 31.0 * es4 / 180.0 + 517.0 * es6 / 5040.0;
00204 c1 = 23.0 * es4 / 360.0 + 251.0 * es6 / 3780.0;
00205 c2 = 761.0 * es6 / 45360.0;
00206 Sin_Cyeq_Origin_Lat = sin(Cyeq_Origin_Lat);
00207 k0 = cos(Cyeq_Origin_Lat) / sqrt(1.0 - es2 * Sin_Cyeq_Origin_Lat * Sin_Cyeq_Origin_Lat);
00208 Cyeq_a_k0 = semiMajorAxis * k0;
00209 two_k0 = 2.0 * k0;
00210
00211 if (Cyeq_Origin_Long > 0)
00212 {
00213 GeodeticCoordinates gcMax( CoordinateType::geodetic, Cyeq_Origin_Long - PI - ONE, PI_OVER_2 );
00214 MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcMax );
00215 Cyeq_Max_Easting = tempCoordinates->easting();
00216 delete tempCoordinates;
00217
00218 GeodeticCoordinates gcMin( CoordinateType::geodetic, Cyeq_Origin_Long - PI, PI_OVER_2 );
00219 tempCoordinates = convertFromGeodetic( &gcMin );
00220 Cyeq_Min_Easting = tempCoordinates->easting();
00221 delete tempCoordinates;
00222
00223 GeodeticCoordinates gcDelta( CoordinateType::geodetic, PI, PI_OVER_2 );
00224 tempCoordinates = convertFromGeodetic( &gcDelta );
00225 Cyeq_Delta_Northing = tempCoordinates->northing();
00226 delete tempCoordinates;
00227 }
00228 else if (Cyeq_Origin_Long < 0)
00229 {
00230 GeodeticCoordinates gcMax( CoordinateType::geodetic, Cyeq_Origin_Long + PI, PI_OVER_2 );
00231 MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcMax );
00232 Cyeq_Max_Easting = tempCoordinates->easting();
00233 delete tempCoordinates;
00234
00235 GeodeticCoordinates gcMin( CoordinateType::geodetic, Cyeq_Origin_Long + PI + ONE, PI_OVER_2 );
00236 tempCoordinates = convertFromGeodetic( &gcMin );
00237 Cyeq_Min_Easting = tempCoordinates->easting();
00238 delete tempCoordinates;
00239
00240 GeodeticCoordinates gcDelta( CoordinateType::geodetic, PI, PI_OVER_2 );
00241 tempCoordinates = convertFromGeodetic( &gcDelta );
00242 Cyeq_Delta_Northing = tempCoordinates->northing();
00243 delete tempCoordinates;
00244 }
00245 else
00246 {
00247 GeodeticCoordinates gcTemp( CoordinateType::geodetic, PI, PI_OVER_2 );
00248 MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcTemp );
00249 Cyeq_Max_Easting = tempCoordinates->easting();
00250 Cyeq_Delta_Northing = tempCoordinates->northing();
00251 delete tempCoordinates;
00252 Cyeq_Min_Easting = -Cyeq_Max_Easting;
00253 }
00254
00255 if(Cyeq_False_Northing)
00256 Cyeq_Delta_Northing -= Cyeq_False_Northing;
00257 if (Cyeq_Delta_Northing < 0)
00258 Cyeq_Delta_Northing = -Cyeq_Delta_Northing;
00259
00260 if(Cyeq_False_Easting)
00261 {
00262 Cyeq_Min_Easting -= Cyeq_False_Easting;
00263 Cyeq_Max_Easting -= Cyeq_False_Easting;
00264 }
00265 }
00266
00267
00268 CylindricalEqualArea::CylindricalEqualArea( const CylindricalEqualArea &cea )
00269 {
00270 semiMajorAxis = cea.semiMajorAxis;
00271 flattening = cea.flattening;
00272 es = cea.es;
00273 es2 = cea.es2;
00274 es4 = cea.es4;
00275 es6 = cea.es6;
00276 k0 = cea.k0;
00277 Cyeq_a_k0 = cea.Cyeq_a_k0;
00278 two_k0 = cea.two_k0;
00279 c0 = cea.c0;
00280 c1 = cea.c1;
00281 c2 = cea.c2;
00282 Cyeq_Origin_Long = cea.Cyeq_Origin_Long;
00283 Cyeq_Origin_Lat = cea.Cyeq_Origin_Lat;
00284 Cyeq_False_Easting = cea.Cyeq_False_Easting;
00285 Cyeq_False_Northing = cea.Cyeq_False_Northing;
00286 Cyeq_Max_Easting = cea.Cyeq_Max_Easting;
00287 Cyeq_Min_Easting = cea.Cyeq_Min_Easting;
00288 Cyeq_Delta_Northing = cea.Cyeq_Delta_Northing;
00289 }
00290
00291
00292 CylindricalEqualArea::~CylindricalEqualArea()
00293 {
00294 }
00295
00296
00297 CylindricalEqualArea& CylindricalEqualArea::operator=( const CylindricalEqualArea &cea )
00298 {
00299 if( this != &cea )
00300 {
00301 semiMajorAxis = cea.semiMajorAxis;
00302 flattening = cea.flattening;
00303 es = cea.es;
00304 es2 = cea.es2;
00305 es4 = cea.es4;
00306 es6 = cea.es6;
00307 k0 = cea.k0;
00308 Cyeq_a_k0 = cea.Cyeq_a_k0;
00309 two_k0 = cea.two_k0;
00310 c0 = cea.c0;
00311 c1 = cea.c1;
00312 c2 = cea.c2;
00313 Cyeq_Origin_Long = cea.Cyeq_Origin_Long;
00314 Cyeq_Origin_Lat = cea.Cyeq_Origin_Lat;
00315 Cyeq_False_Easting = cea.Cyeq_False_Easting;
00316 Cyeq_False_Northing = cea.Cyeq_False_Northing;
00317 Cyeq_Max_Easting = cea.Cyeq_Max_Easting;
00318 Cyeq_Min_Easting = cea.Cyeq_Min_Easting;
00319 Cyeq_Delta_Northing = cea.Cyeq_Delta_Northing;
00320 }
00321
00322 return *this;
00323 }
00324
00325
00326 MapProjection4Parameters* CylindricalEqualArea::getParameters() const
00327 {
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344 return new MapProjection4Parameters( CoordinateType::cylindricalEqualArea, Cyeq_Origin_Long, Cyeq_Origin_Lat, Cyeq_False_Easting, Cyeq_False_Northing );
00345 }
00346
00347
00348 MSP::CCS::MapProjectionCoordinates* CylindricalEqualArea::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00349 {
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363 double dlam;
00364 double qq;
00365 double x;
00366 char errorStatus[50] = "";
00367
00368 double longitude = geodeticCoordinates->longitude();
00369 double latitude = geodeticCoordinates->latitude();
00370 double sin_lat = sin(latitude);
00371
00372 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00373 {
00374 strcat( errorStatus, ErrorMessages::latitude );
00375 }
00376 if ((longitude < -PI) || (longitude > TWO_PI))
00377 {
00378 strcat( errorStatus, ErrorMessages::longitude );
00379 }
00380
00381 if( strlen( errorStatus ) > 0)
00382 throw CoordinateConversionException( errorStatus );
00383
00384 dlam = longitude - Cyeq_Origin_Long;
00385 if (dlam > PI)
00386 {
00387 dlam -= TWO_PI;
00388 }
00389 if (dlam < -PI)
00390 {
00391 dlam += TWO_PI;
00392 }
00393 x = es * sin_lat;
00394 qq = cyleqarQ( sin_lat, x );
00395
00396 double easting = Cyeq_a_k0 * dlam + Cyeq_False_Easting;
00397 double northing = semiMajorAxis * qq / two_k0 + Cyeq_False_Northing;
00398
00399 return new MapProjectionCoordinates( CoordinateType::cylindricalEqualArea, easting, northing );
00400 }
00401
00402
00403 MSP::CCS::GeodeticCoordinates* CylindricalEqualArea::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00404 {
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419 double sin2beta, sin4beta, sin6beta;
00420 double dx;
00421 double dy;
00422 double qp;
00423 double beta;
00424 double sin_lat = sin(PI_OVER_2);
00425 double i;
00426 double x;
00427 char errorStatus[50] = "";
00428
00429 double easting = mapProjectionCoordinates->easting();
00430 double northing = mapProjectionCoordinates->northing();
00431
00432 if ((easting < (Cyeq_False_Easting + Cyeq_Min_Easting))
00433 || (easting > (Cyeq_False_Easting + Cyeq_Max_Easting)))
00434 {
00435 strcat( errorStatus, ErrorMessages::easting );
00436 }
00437 if ((northing < (Cyeq_False_Northing - fabs(Cyeq_Delta_Northing)))
00438 || (northing > (Cyeq_False_Northing + fabs(Cyeq_Delta_Northing))))
00439 {
00440 strcat( errorStatus, ErrorMessages::northing );
00441 }
00442
00443 if( strlen( errorStatus ) > 0)
00444 throw CoordinateConversionException( errorStatus );
00445
00446 dy = northing - Cyeq_False_Northing;
00447 dx = easting - Cyeq_False_Easting;
00448 x = es * sin_lat;
00449 qp = cyleqarQ( sin_lat, x );
00450 i = two_k0 * dy / (semiMajorAxis * qp);
00451 if (i > 1.0)
00452 i = 1.0;
00453 else if (i < -1.0)
00454 i = -1.0;
00455 beta = asin(i);
00456 sin2beta = cyeqCoeffTimesSine(c0, 2.0, beta);
00457 sin4beta = cyeqCoeffTimesSine(c1, 4.0, beta);
00458 sin6beta = cyeqCoeffTimesSine(c2, 6.0, beta);
00459
00460 double latitude = beta + sin2beta + sin4beta + sin6beta;
00461 double longitude = Cyeq_Origin_Long + dx / Cyeq_a_k0;
00462
00463 if (latitude > PI_OVER_2)
00464 latitude = PI_OVER_2;
00465 else if (latitude < -PI_OVER_2)
00466 latitude = -PI_OVER_2;
00467
00468 if (longitude > PI)
00469 longitude -= TWO_PI;
00470 if (longitude < -PI)
00471 longitude += TWO_PI;
00472
00473 if (longitude > PI)
00474 longitude = PI;
00475 else if (longitude < -PI)
00476 longitude = -PI;
00477
00478 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00479 }
00480
00481
00482 double CylindricalEqualArea::cyleqarQ( double slat, double x )
00483 {
00484 return (1.0-es2)*(slat/(1.0-x*x)-(1.0/(2.0*es)) * log((1.0-x)/(1.0+x)));
00485 }
00486
00487
00488
00489