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 #include <math.h>
00085 #include "Bonne.h"
00086 #include "Sinusoidal.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
00104
00105 using namespace MSP::CCS;
00106
00107
00108
00109
00110
00111
00112
00113 const double PI = 3.14159265358979323e0;
00114 const double PI_OVER_2 = (PI / 2.0);
00115 const double TWO_PI = (2.0 * PI);
00116
00117
00118 double bonnCoeffTimesSine( double coeff, double x, double latit )
00119 {
00120 return coeff * (sin(x * latit));
00121 }
00122
00123
00124
00125
00126
00127
00128
00129 Bonne::Bonne( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double originLatitude, double falseEasting, double falseNorthing ) :
00130 CoordinateSystem(),
00131 sinusoidal( 0 ),
00132 es2( 0.0066943799901413800 ),
00133 es4( 4.4814723452405e-005 ),
00134 es6( 3.0000678794350e-007 ),
00135 M1( 4984944.3782319 ),
00136 m1( .70829317069372 ),
00137 c0( .99832429845280 ),
00138 c1( .0025146070605187 ),
00139 c2( 2.6390465943377e-006 ),
00140 c3( 3.4180460865959e-009 ),
00141 a0( .0025188265843907 ),
00142 a1( 3.7009490356205e-006 ),
00143 a2( 7.4478137675038e-009 ),
00144 a3( 1.7035993238596e-011 ),
00145 Bonn_Origin_Long( 0.0 ),
00146 Bonn_Origin_Lat( ((45 * PI) / 180.0) ),
00147 Bonn_False_Easting( 0.0 ),
00148 Bonn_False_Northing( 0.0 ),
00149 Sin_Bonn_Origin_Lat( .70710678118655 ),
00150 Bonn_am1sin( 6388838.2901211 ),
00151 Bonn_Max_Easting( 20027474.0 ),
00152 Bonn_Min_Easting( -20027474.0 ),
00153 Bonn_Delta_Northing( 20003932.0 )
00154 {
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173 double j, three_es4;
00174 double x,e1,e2,e3,e4;
00175 double clat;
00176 double sin2lat, sin4lat, sin6lat, lat;
00177 double inv_f = 1 / ellipsoidFlattening;
00178 char errorStatus[500] = "";
00179
00180 if (ellipsoidSemiMajorAxis <= 0.0)
00181 {
00182 strcat( errorStatus, MSP::CCS::ErrorMessages::semiMajorAxis );
00183 }
00184 if ((inv_f < 250) || (inv_f > 350))
00185 {
00186 strcat( errorStatus, MSP::CCS::ErrorMessages::ellipsoidFlattening );
00187 }
00188 if ((originLatitude < -PI_OVER_2) || (originLatitude > PI_OVER_2))
00189 {
00190 strcat( errorStatus, MSP::CCS::ErrorMessages::originLatitude );
00191 }
00192 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00193 {
00194 strcat( errorStatus, MSP::CCS::ErrorMessages::centralMeridian );
00195 }
00196
00197 if( strlen( errorStatus ) > 0)
00198 throw CoordinateConversionException( errorStatus );
00199
00200 semiMajorAxis = ellipsoidSemiMajorAxis;
00201 flattening = ellipsoidFlattening;
00202
00203 Bonn_Origin_Lat = originLatitude;
00204 if (centralMeridian > PI)
00205 centralMeridian -= TWO_PI;
00206 Bonn_Origin_Long = centralMeridian;
00207 Bonn_False_Northing = falseNorthing;
00208 Bonn_False_Easting = falseEasting;
00209 if (Bonn_Origin_Lat == 0.0)
00210 {
00211 if (Bonn_Origin_Long > 0)
00212 {
00213 Bonn_Max_Easting = 19926189.0;
00214 Bonn_Min_Easting = -20037509.0;
00215 }
00216 else if (Bonn_Origin_Long < 0)
00217 {
00218 Bonn_Max_Easting = 20037509.0;
00219 Bonn_Min_Easting = -19926189.0;
00220 }
00221 else
00222 {
00223 Bonn_Max_Easting = 20037509.0;
00224 Bonn_Min_Easting = -20037509.0;
00225 }
00226 Bonn_Delta_Northing = 10001966.0;
00227
00228 sinusoidal = new Sinusoidal( semiMajorAxis, flattening, Bonn_Origin_Long, Bonn_False_Easting, Bonn_False_Northing );
00229 }
00230 else
00231 {
00232 Sin_Bonn_Origin_Lat = sin(Bonn_Origin_Lat);
00233
00234 es2 = 2 * flattening - flattening * flattening;
00235 es4 = es2 * es2;
00236 es6 = es4 * es2;
00237 j = 45.0 * es6 / 1024.0;
00238 three_es4 = 3.0 * es4;
00239 c0 = 1 - es2 / 4.0 - three_es4 / 64.0 - 5.0 * es6 / 256.0;
00240 c1 = 3.0 * es2 / 8.0 + three_es4 / 32.0 + j;
00241 c2 = 15.0 * es4 / 256.0 + j;
00242 c3 = 35.0 * es6 / 3072.0;
00243
00244 clat = cos(Bonn_Origin_Lat);
00245 m1 = bonnm(clat, Sin_Bonn_Origin_Lat);
00246
00247 lat = c0 * Bonn_Origin_Lat;
00248 sin2lat = bonnCoeffTimesSine(c1, 2.0, Bonn_Origin_Lat);
00249 sin4lat = bonnCoeffTimesSine(c2, 4.0, Bonn_Origin_Lat);
00250 sin6lat = bonnCoeffTimesSine(c3, 6.0, Bonn_Origin_Lat);
00251 M1 = bonnM(lat, sin2lat, sin4lat, sin6lat);
00252
00253 x = sqrt (1.0 - es2);
00254 e1 = (1.0 - x) / (1.0 + x);
00255 e2 = e1 * e1;
00256 e3 = e2 * e1;
00257 e4 = e3 * e1;
00258 a0 = 3.0 * e1 / 2.0 - 27.0 * e3 / 32.0;
00259 a1 = 21.0 * e2 / 16.0 - 55.0 * e4 / 32.0;
00260 a2 = 151.0 * e3 / 96.0;
00261 a3 = 1097.0 * e4 / 512.0;
00262 if (Sin_Bonn_Origin_Lat == 0.0)
00263 Bonn_am1sin = 0.0;
00264 else
00265 Bonn_am1sin = semiMajorAxis * m1 / Sin_Bonn_Origin_Lat;
00266
00267 Bonn_Max_Easting = 20027474.0;
00268 Bonn_Min_Easting = -20027474.0;
00269 Bonn_Delta_Northing = 20003932.0;
00270 }
00271 }
00272
00273
00274 Bonne::Bonne( const Bonne &b )
00275 {
00276 sinusoidal = new Sinusoidal( *( b.sinusoidal ) );
00277 semiMajorAxis = b.semiMajorAxis;
00278 flattening = b.flattening;
00279 es2 = b.es2;
00280 es4 = b.es4;
00281 es6 = b.es6;
00282 M1 = b.M1;
00283 m1 = b.m1;
00284 c0 = b.c0;
00285 c1 = b.c1;
00286 c2 = b.c2;
00287 c3 = b.c3;
00288 a0 = b.a0;
00289 a1 = b.a1;
00290 a2 = b.a2;
00291 a3 = b.a3;
00292 Bonn_Origin_Long = b.Bonn_Origin_Long;
00293 Bonn_Origin_Lat = b.Bonn_Origin_Lat;
00294 Bonn_False_Easting = b.Bonn_False_Easting;
00295 Bonn_False_Northing = b.Bonn_False_Northing;
00296 Sin_Bonn_Origin_Lat = b.Sin_Bonn_Origin_Lat;
00297 Bonn_am1sin = b.Bonn_am1sin;
00298 Bonn_Max_Easting = b.Bonn_Max_Easting;
00299 Bonn_Min_Easting = b.Bonn_Min_Easting;
00300 Bonn_Delta_Northing = b.Bonn_Delta_Northing;
00301 }
00302
00303
00304 Bonne::~Bonne()
00305 {
00306 if( sinusoidal )
00307 {
00308 delete sinusoidal;
00309 sinusoidal = 0;
00310 }
00311 }
00312
00313
00314 Bonne& Bonne::operator=( const Bonne &b )
00315 {
00316 if( this != &b )
00317 {
00318 sinusoidal->operator=( *b.sinusoidal );
00319 semiMajorAxis = b.semiMajorAxis;
00320 flattening = b.flattening;
00321 es2 = b.es2;
00322 es4 = b.es4;
00323 es6 = b.es6;
00324 M1 = b.M1;
00325 m1 = b.m1;
00326 c0 = b.c0;
00327 c1 = b.c1;
00328 c2 = b.c2;
00329 c3 = b.c3;
00330 a0 = b.a0;
00331 a1 = b.a1;
00332 a2 = b.a2;
00333 a3 = b.a3;
00334 Bonn_Origin_Long = b.Bonn_Origin_Long;
00335 Bonn_Origin_Lat = b.Bonn_Origin_Lat;
00336 Bonn_False_Easting = b.Bonn_False_Easting;
00337 Bonn_False_Northing = b.Bonn_False_Northing;
00338 Sin_Bonn_Origin_Lat = b.Sin_Bonn_Origin_Lat;
00339 Bonn_am1sin = b.Bonn_am1sin;
00340 Bonn_Max_Easting = b.Bonn_Max_Easting;
00341 Bonn_Min_Easting = b.Bonn_Min_Easting;
00342 Bonn_Delta_Northing = b.Bonn_Delta_Northing;
00343 }
00344
00345 return *this;
00346 }
00347
00348
00349 MapProjection4Parameters* Bonne::getParameters() const
00350 {
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367 return new MapProjection4Parameters( CoordinateType::bonne, Bonn_Origin_Long, Bonn_Origin_Lat, Bonn_False_Easting, Bonn_False_Northing );
00368 }
00369
00370
00371 MSP::CCS::MapProjectionCoordinates* Bonne::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00372 {
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386 double dlam;
00387 double mm;
00388 double MM;
00389 double rho;
00390 double EE;
00391 double lat, sin2lat, sin4lat, sin6lat;
00392 double easting, northing;
00393 char errorStatus[50] = "";
00394
00395 double longitude = geodeticCoordinates->longitude();
00396 double latitude = geodeticCoordinates->latitude();
00397 double clat = cos(latitude);
00398 double slat = sin(latitude);
00399
00400 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00401 {
00402 strcat( errorStatus, MSP::CCS::ErrorMessages::latitude );
00403 }
00404 if ((longitude < -PI) || (longitude > TWO_PI))
00405 {
00406 strcat( errorStatus, MSP::CCS::ErrorMessages::longitude );
00407 }
00408
00409 if( strlen( errorStatus ) > 0)
00410 throw CoordinateConversionException( errorStatus );
00411
00412 if (Bonn_Origin_Lat == 0.0)
00413 return sinusoidal->convertFromGeodetic( geodeticCoordinates );
00414 else
00415 {
00416 dlam = longitude - Bonn_Origin_Long;
00417 if (dlam > PI)
00418 {
00419 dlam -= TWO_PI;
00420 }
00421 if (dlam < -PI)
00422 {
00423 dlam += TWO_PI;
00424 }
00425 if ((latitude - Bonn_Origin_Lat) == 0.0 && floatEq(fabs(latitude),PI_OVER_2,.00001))
00426 {
00427 easting = 0.0;
00428 northing = 0.0;
00429 }
00430 else
00431 {
00432 mm = bonnm(clat, slat);
00433 lat = c0 * latitude;
00434 sin2lat = bonnCoeffTimesSine(c1, 2.0, latitude);
00435 sin4lat = bonnCoeffTimesSine(c2, 4.0, latitude);
00436 sin6lat = bonnCoeffTimesSine(c3, 6.0, latitude);
00437 MM = bonnM(lat, sin2lat, sin4lat, sin6lat);
00438
00439 rho = Bonn_am1sin + M1 - MM;
00440 if (rho == 0)
00441 EE = 0;
00442 else
00443 EE = semiMajorAxis * mm * dlam / rho;
00444 easting = rho * sin(EE) + Bonn_False_Easting;
00445 northing = Bonn_am1sin - rho * cos(EE) + Bonn_False_Northing;
00446 }
00447
00448 return new MapProjectionCoordinates( CoordinateType::bonne, easting, northing );
00449 }
00450 }
00451
00452
00453 MSP::CCS::GeodeticCoordinates* Bonne::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00454 {
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468 double dx;
00469 double dy;
00470 double mu;
00471 double MM;
00472 double mm;
00473 double am1sin_dy;
00474 double rho;
00475 double sin2mu, sin4mu, sin6mu, sin8mu;
00476 double clat, slat;
00477 double longitude, latitude;
00478 char errorStatus[50] = "";
00479
00480 double easting = mapProjectionCoordinates->easting();
00481 double northing = mapProjectionCoordinates->northing();
00482
00483 if ((easting < (Bonn_False_Easting + Bonn_Min_Easting))
00484 || (easting > (Bonn_False_Easting + Bonn_Max_Easting)))
00485 {
00486 strcat( errorStatus, MSP::CCS::ErrorMessages::easting );
00487 }
00488 if ((northing < (Bonn_False_Northing - Bonn_Delta_Northing))
00489 || (northing > (Bonn_False_Northing + Bonn_Delta_Northing)))
00490 {
00491 strcat( errorStatus, MSP::CCS::ErrorMessages::northing );
00492 }
00493
00494 if( strlen( errorStatus ) > 0)
00495 throw CoordinateConversionException( errorStatus );
00496
00497 if (Bonn_Origin_Lat == 0.0)
00498 return sinusoidal->convertToGeodetic( mapProjectionCoordinates );
00499 else
00500 {
00501 dy = northing - Bonn_False_Northing;
00502 dx = easting - Bonn_False_Easting;
00503 am1sin_dy = Bonn_am1sin - dy;
00504 rho = sqrt(dx * dx + am1sin_dy * am1sin_dy);
00505 if (Bonn_Origin_Lat < 0.0)
00506 rho = -rho;
00507 MM = Bonn_am1sin + M1 - rho;
00508
00509 mu = MM / (semiMajorAxis * c0);
00510 sin2mu = bonnCoeffTimesSine(a0, 2.0, mu);
00511 sin4mu = bonnCoeffTimesSine(a1, 4.0, mu);
00512 sin6mu = bonnCoeffTimesSine(a2, 6.0, mu);
00513 sin8mu = bonnCoeffTimesSine(a3, 8.0, mu);
00514 latitude = mu + sin2mu + sin4mu + sin6mu + sin8mu;
00515
00516 if (floatEq(fabs(latitude),PI_OVER_2,.00001))
00517 {
00518 longitude = Bonn_Origin_Long;
00519 }
00520 else
00521 {
00522 clat = cos(latitude);
00523 slat = sin(latitude);
00524 mm = bonnm(clat, slat);
00525
00526 if (Bonn_Origin_Lat < 0.0)
00527 {
00528 dx = -dx;
00529 am1sin_dy = -am1sin_dy;
00530 }
00531 longitude = Bonn_Origin_Long + rho * (atan2(dx, am1sin_dy)) /
00532 (semiMajorAxis * mm);
00533 }
00534
00535 if (latitude > PI_OVER_2)
00536 latitude = PI_OVER_2;
00537 else if (latitude < -PI_OVER_2)
00538 latitude = -PI_OVER_2;
00539
00540 if (longitude > PI)
00541 longitude -= TWO_PI;
00542 if (longitude < -PI)
00543 longitude += TWO_PI;
00544
00545 if (longitude > PI)
00546 longitude = PI;
00547 else if (longitude < -PI)
00548 longitude = -PI;
00549
00550 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00551 }
00552 }
00553
00554
00555 double Bonne::bonnm( double coslat, double sinlat )
00556 {
00557 return coslat/sqrt(1.0 - es2*sinlat*sinlat);
00558 }
00559
00560
00561 double Bonne::bonnM( double c0lat, double c1s2lat, double c2s4lat, double c3s6lat )
00562 {
00563 return semiMajorAxis*(c0lat-c1s2lat+c2s4lat-c3s6lat);
00564 }
00565
00566
00567 double Bonne::floatEq( double x, double v, double epsilon )
00568 {
00569 return ((v - epsilon) < x) && (x < (v + epsilon));
00570 }
00571
00572
00573
00574