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 <ctype.h>
00085 #include <math.h>
00086 #include <stdio.h>
00087 #include <string.h>
00088 #include "TransverseMercator.h"
00089 #include "BritishNationalGrid.h"
00090 #include "BNGCoordinates.h"
00091 #include "EllipsoidParameters.h"
00092 #include "MapProjectionCoordinates.h"
00093 #include "GeodeticCoordinates.h"
00094 #include "CoordinateConversionException.h"
00095 #include "ErrorMessages.h"
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112 using namespace MSP::CCS;
00113
00114
00115
00116
00117
00118
00119
00120 const double PI = 3.14159265358979323e0;
00121 const double PI_OVER_2 = (PI / 2.0e0);
00122 const double TWO_PI = (2.0e0 * PI);
00123 const double MAX_LAT = (61.5 * PI / 180.0);
00124 const double MIN_LAT = (49.5 * PI / 180.0);
00125 const double MAX_LON = (3.5 * PI / 180.0);
00126 const double MIN_LON = (-10.0 * PI / 180.0);
00127 const char* BNG500GRID = "STNOHJ";
00128 const char* BNG100GRID = "VWXYZQRSTULMNOPFGHJKABCDE";
00129
00130
00131 const double BNG_Origin_Lat = (49.0 * PI / 180.0);
00132 const double BNG_Origin_Long = (-2.0 * PI / 180.0);
00133 const double BNG_False_Northing = -100000.0;
00134 const double BNG_False_Easting = 400000.0;
00135 const double BNG_Scale_Factor = .9996012717;
00136
00137
00138 const double BNG_Max_Easting = 759961.0;
00139 const double BNG_Max_Northing = 1257875.0;
00140 const double BNG_Min_Easting = -133134.0;
00141 const double BNG_Min_Northing = -14829.0;
00142
00143 static const char* Airy = "AA";
00144
00145
00146
00147
00148
00149
00150
00151 void findIndex( char letter, const char* letterArray, long *index )
00152 {
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163 long i = 0;
00164 long not_Found = 1;
00165 long length = strlen(letterArray);
00166
00167 while ((i < length) && (not_Found))
00168 {
00169 if (letter == letterArray[i])
00170 {
00171 *index = i;
00172 not_Found = 0;
00173 }
00174 i++;
00175 }
00176 if (not_Found)
00177 throw CoordinateConversionException( ErrorMessages::bngString );
00178 }
00179
00180
00181 long roundBNG( double value )
00182 {
00183
00184 double ivalue;
00185 long ival;
00186 double fraction = modf (value, &ivalue);
00187 ival = (long)(ivalue);
00188 if ((fraction > 0.5) || ((fraction == 0.5) && (ival%2 == 1)))
00189 ival++;
00190
00191 return (ival);
00192 }
00193
00194
00195 void makeBNGString( char ltrnum[4], long easting, long northing, char* BNGString, long precision )
00196
00197 {
00198
00199 double divisor;
00200 long east;
00201 long north;
00202 long i;
00203 long j;
00204 double unitInterval;
00205
00206 i = 0;
00207 for (j = 0; j < 3; j++)
00208 BNGString[i++] = ltrnum[j];
00209 divisor = pow (10.0, (5.0 - precision));
00210 unitInterval = pow (10.0, (double)precision);
00211
00212 east = roundBNG (easting/divisor);
00213 if (east == unitInterval)
00214 east -= 1;
00215 if ((precision == 0) && (east == 1))
00216 east = 0;
00217 i += sprintf (BNGString + i, "%*.*ld", precision, precision, east);
00218
00219 north = roundBNG (northing/divisor);
00220 if (north == unitInterval)
00221 north -= 1;
00222 if ((precision == 0) && (north == 1))
00223 north = 0;
00224 i += sprintf (BNGString + i, "%*.*ld", precision, precision, north);
00225 }
00226
00227
00228 bool checkOutOfArea( char BNG500, char BNG100 )
00229
00230 {
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241 bool outOfArea = false;
00242
00243 switch (BNG500)
00244 {
00245 case 'S':
00246 switch (BNG100)
00247 {
00248 case 'A':
00249 case 'F':
00250 case 'L':
00251 outOfArea = true;
00252 break;
00253 default:
00254 break;
00255 }
00256 break;
00257 case 'N':
00258 switch (BNG100)
00259 {
00260 case 'V':
00261 outOfArea = true;
00262 break;
00263 default:
00264 break;
00265 }
00266 break;
00267
00268 case 'H':
00269 if (BNG100 < 'L')
00270 outOfArea = true;
00271 break;
00272 case 'T':
00273 switch (BNG100)
00274 {
00275 case 'D':
00276 case 'E':
00277 case 'J':
00278 case 'K':
00279 case 'O':
00280 case 'P':
00281 case 'T':
00282 case 'U':
00283 case 'X':
00284 case 'Y':
00285 case 'Z':
00286 outOfArea = true;
00287 break;
00288 default:
00289 break;
00290 }
00291 break;
00292 case 'O':
00293 switch (BNG100)
00294 {
00295 case 'C':
00296 case 'D':
00297 case 'E':
00298 case 'J':
00299 case 'K':
00300 case 'O':
00301 case 'P':
00302 case 'T':
00303 case 'U':
00304 case 'Y':
00305 case 'Z':
00306 outOfArea = true;
00307 break;
00308 default:
00309 break;
00310 }
00311 break;
00312 case 'J':
00313 switch (BNG100)
00314 {
00315 case 'L':
00316 case 'M':
00317 case 'Q':
00318 case 'R':
00319 case 'V':
00320 case 'W':
00322 break;
00323 default:
00324 outOfArea = true;
00325 break;
00326 }
00327 break;
00328 default:
00329 outOfArea = true;
00330 break;
00331 }
00332
00333 return outOfArea;
00334 }
00335
00336
00337 void breakBNGString( char* BNGString, char letters[3], double* easting, double* northing, long* precision )
00338 {
00339
00340
00341 long i = 0;
00342 long j;
00343 long num_digits = 0;
00344 long num_letters;
00345 long length = strlen(BNGString);
00346
00347 while (BNGString[i] == ' ')
00348 i++;
00349 j = i;
00350 while (isalpha(BNGString[i]))
00351 i++;
00352 num_letters = i - j;
00353 if (num_letters == 2)
00354 {
00355
00356 letters[0] = (char)toupper(BNGString[j]);
00357 letters[1] = (char)toupper(BNGString[j+1]);
00358 letters[2] = 0;
00359 }
00360 else
00361 throw CoordinateConversionException( ErrorMessages::bngString );
00362
00363 checkOutOfArea(letters[0], letters[1]);
00364 while (BNGString[i] == ' ')
00365 i++;
00366 j = i;
00367 if (BNGString[length-1] == ' ')
00368 length --;
00369 while (i < length)
00370 {
00371 if (isdigit(BNGString[i]))
00372 i++;
00373 else
00374 throw CoordinateConversionException( ErrorMessages::bngString );
00375 }
00376
00377 num_digits = i - j;
00378 if ((num_digits <= 10) && (num_digits%2 == 0))
00379 {
00380 long n;
00381 char east_string[6];
00382 char north_string[6];
00383 long east;
00384 long north;
00385 double multiplier;
00386
00387
00388 n = num_digits / 2;
00389 *precision = n;
00390 if (n > 0)
00391 {
00392 strncpy (east_string, BNGString+j, n);
00393 east_string[n] = 0;
00394 sscanf (east_string, "%ld", &east);
00395 strncpy (north_string, BNGString+j+n, n);
00396 north_string[n] = 0;
00397 sscanf (north_string, "%ld", &north);
00398 multiplier = pow (10.0, 5.0 - n);
00399 *easting = east * multiplier;
00400 *northing = north * multiplier;
00401 }
00402 else
00403 {
00404 *easting = 0.0;
00405 *northing = 0.0;
00406 }
00407 }
00408 else
00409 throw CoordinateConversionException( ErrorMessages::bngString );
00410 }
00411
00412
00413
00414
00415
00416
00417
00418 BritishNationalGrid::BritishNationalGrid( char *ellipsoidCode ) :
00419 CoordinateSystem( 6377563.396, 1 / 299.324964600 ),
00420 transverseMercator( 0 ),
00421 BNG_Easting( 0.0 ),
00422 BNG_Northing( 0.0 )
00423 {
00424
00425
00426
00427
00428
00429
00430
00431
00432 strcpy( BNG_Letters, "SV" );
00433 strcpy( BNG_Ellipsoid_Code, "AA");
00434
00435 if ( strcmp( ellipsoidCode, Airy ) != 0 )
00436 {
00437 throw CoordinateConversionException( ErrorMessages::bngEllipsoid );
00438 }
00439
00440 strcpy( BNG_Ellipsoid_Code, ellipsoidCode );
00441
00442 transverseMercator = new TransverseMercator( semiMajorAxis, flattening, BNG_Origin_Long, BNG_Origin_Lat,
00443 BNG_False_Easting, BNG_False_Northing, BNG_Scale_Factor );
00444 }
00445
00446
00447 BritishNationalGrid::BritishNationalGrid( const BritishNationalGrid &bng )
00448 {
00449 transverseMercator = new TransverseMercator( *( bng.transverseMercator ) );
00450 semiMajorAxis = bng.semiMajorAxis;
00451 flattening = bng.flattening;
00452 BNG_Easting = bng.BNG_Easting;
00453 BNG_Northing = bng.BNG_Northing;
00454 }
00455
00456
00457 BritishNationalGrid::~BritishNationalGrid()
00458 {
00459 delete transverseMercator;
00460 transverseMercator = 0;
00461 }
00462
00463
00464 BritishNationalGrid& BritishNationalGrid::operator=( const BritishNationalGrid &bng )
00465 {
00466 if( this != &bng )
00467 {
00468 transverseMercator->operator=( *bng.transverseMercator );
00469 semiMajorAxis = bng.semiMajorAxis;
00470 flattening = bng.flattening;
00471 BNG_Easting = bng.BNG_Easting;
00472 BNG_Northing = bng.BNG_Northing;
00473 }
00474
00475 return *this;
00476 }
00477
00478
00479 EllipsoidParameters* BritishNationalGrid::getParameters() const
00480 {
00481
00482
00483
00484
00485
00486
00487
00488 return new EllipsoidParameters( semiMajorAxis, flattening, (char*)BNG_Ellipsoid_Code );
00489 }
00490
00491
00492 MSP::CCS::BNGCoordinates* BritishNationalGrid::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates, long precision )
00493 {
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507 double TMEasting, TMNorthing;
00508 char errorStatus[50] = "";
00509
00510 double longitude = geodeticCoordinates->longitude();
00511 double latitude = geodeticCoordinates->latitude();
00512
00513 if ((latitude < MIN_LAT) || (latitude > MAX_LAT))
00514 {
00515 strcat( errorStatus, MSP::CCS::ErrorMessages::latitude );
00516 }
00517 if ((longitude < MIN_LON) || (longitude > MAX_LON))
00518 {
00519 strcat( errorStatus, MSP::CCS::ErrorMessages::longitude );
00520 }
00521
00522 if( strlen( errorStatus ) > 0)
00523 throw CoordinateConversionException( errorStatus );
00524
00525 MapProjectionCoordinates* transverseMercatorCoordinates = transverseMercator->convertFromGeodetic( geodeticCoordinates );
00526
00527 TMEasting = transverseMercatorCoordinates->easting();
00528 TMNorthing = transverseMercatorCoordinates->northing();
00529
00530 if ((TMEasting < 0.0) && (TMEasting > -2.0))
00531 transverseMercatorCoordinates->setEasting( 0.0 );
00532 if ((TMNorthing < 0.0) && (TMNorthing > -2.0))
00533 transverseMercatorCoordinates->setNorthing( 0.0 );
00534
00535 TMEasting = transverseMercatorCoordinates->easting();
00536 TMNorthing = transverseMercatorCoordinates->northing();
00537
00538 if ((TMEasting < BNG_Min_Easting) || (TMEasting > BNG_Max_Easting))
00539 {
00540 delete transverseMercatorCoordinates;
00541 throw CoordinateConversionException( ErrorMessages::invalidArea );
00542 }
00543 if ((TMNorthing < BNG_Min_Northing) || (TMNorthing > BNG_Max_Northing))
00544 {
00545 delete transverseMercatorCoordinates;
00546 throw CoordinateConversionException( ErrorMessages::invalidArea );
00547 }
00548
00549 BNGCoordinates* bngCoordinates = 0;
00550 try {
00551 bngCoordinates = convertFromTransverseMercator( transverseMercatorCoordinates, precision );
00552 delete transverseMercatorCoordinates;
00553 }
00554 catch ( CoordinateConversionException e ){
00555 delete transverseMercatorCoordinates;
00556 throw e;
00557 }
00558
00559 return bngCoordinates;
00560 }
00561
00562
00563 MSP::CCS::GeodeticCoordinates* BritishNationalGrid::convertToGeodetic( MSP::CCS::BNGCoordinates* bngCoordinates )
00564 {
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577 double TMEasting, TMNorthing;
00578 long in_Precision;
00579
00580 char* BNGString = bngCoordinates->BNGString();
00581
00582 breakBNGString( BNGString, BNG_Letters, &BNG_Easting, &BNG_Northing, &in_Precision );
00583
00584 MapProjectionCoordinates* transverseMercatorCoordinates = convertToTransverseMercator( bngCoordinates );
00585 TMEasting = transverseMercatorCoordinates->easting();
00586 TMNorthing = transverseMercatorCoordinates->northing();
00587
00588 if ((TMEasting < BNG_Min_Easting) || (TMEasting > BNG_Max_Easting))
00589 {
00590 delete transverseMercatorCoordinates;
00591 throw CoordinateConversionException( ErrorMessages::invalidArea );
00592 }
00593 if ((TMNorthing < BNG_Min_Northing) || (TMNorthing > BNG_Max_Northing))
00594 {
00595 delete transverseMercatorCoordinates;
00596 throw CoordinateConversionException( ErrorMessages::invalidArea );
00597 }
00598
00599 GeodeticCoordinates* geodeticCoordinates = 0;
00600 double latitude;
00601 double longitude;
00602 try {
00603 geodeticCoordinates = transverseMercator->convertToGeodetic( transverseMercatorCoordinates );
00604 latitude = geodeticCoordinates->latitude();
00605 longitude = geodeticCoordinates->longitude();
00606 delete transverseMercatorCoordinates;
00607 }
00608 catch ( CoordinateConversionException e ){
00609 delete transverseMercatorCoordinates;
00610 throw e;
00611 }
00612
00613 if ((latitude < MIN_LAT) || (latitude > MAX_LAT))
00614 {
00615 delete geodeticCoordinates;
00616 throw CoordinateConversionException( ErrorMessages::invalidArea );
00617 }
00618 if ((longitude < MIN_LON) || (longitude > MAX_LON))
00619 {
00620 delete geodeticCoordinates;
00621 throw CoordinateConversionException( ErrorMessages::invalidArea );
00622 }
00623
00624 return geodeticCoordinates;
00625 }
00626
00627
00628 MSP::CCS::BNGCoordinates* BritishNationalGrid::convertFromTransverseMercator( MapProjectionCoordinates* mapProjectionCoordinates, long precision )
00629 {
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642 char letters[4];
00643 long x, y;
00644 long index;
00645 long temp_Easting, temp_Northing;
00646 char BNGString[21];
00647 char errorStatus[50] = "";
00648
00649 double easting = mapProjectionCoordinates->easting();
00650 double northing = mapProjectionCoordinates->northing();
00651
00652 if ((easting < BNG_Min_Easting) || (easting > BNG_Max_Easting))
00653 {
00654 strcat( errorStatus, ErrorMessages::easting );
00655 }
00656 if ((northing < BNG_Min_Northing) || (northing > BNG_Max_Northing))
00657 {
00658 strcat( errorStatus, ErrorMessages::northing );
00659 }
00660
00661 if( strlen( errorStatus ) > 0)
00662 throw CoordinateConversionException( errorStatus );
00663
00664 temp_Easting = roundBNG(easting);
00665 temp_Northing = roundBNG(northing);
00666
00667 temp_Easting += 1000000;
00668 temp_Northing += 500000;
00669
00670 x = temp_Easting / 500000;
00671 y = temp_Northing / 500000;
00672
00673 index = y * 5 + x;
00674 if ((index >= 0) && (index < 25))
00675 {
00676 letters[0] = BNG100GRID[index];
00677 temp_Easting %= 500000;
00678 temp_Northing %= 500000;
00679 x = temp_Easting / 100000;
00680 y = temp_Northing / 100000;
00681 index = y * 5 + x;
00682 if ((index >= 0) && (index < 25))
00683 {
00684 letters[1] = BNG100GRID[index];
00685
00686 if( checkOutOfArea(letters[0], letters[1]) )
00687 throw CoordinateConversionException( ErrorMessages::invalidArea );
00688
00689 letters[2] = ' ';
00690 letters[3] = 0;
00691 temp_Easting %= 100000;
00692 temp_Northing %= 100000;
00693 makeBNGString(letters, temp_Easting, temp_Northing, BNGString, precision);
00694 }
00695 else
00696 {
00697 long five_y = 5 * y;
00698 if ((x >= (25 - five_y)) || (x < -five_y))
00699 throw CoordinateConversionException( ErrorMessages::easting );
00700 if ((five_y >= (25 - x)) || (five_y < -x))
00701 throw CoordinateConversionException( ErrorMessages::northing );
00702 }
00703 }
00704 else
00705 {
00706 long five_y = 5 * y;
00707 if ((x >= (25 - five_y)) || (x < -five_y))
00708 throw CoordinateConversionException( ErrorMessages::easting );
00709 if ((five_y >= (25 - x)) || (five_y < -x))
00710 throw CoordinateConversionException( ErrorMessages::northing );
00711 }
00712
00713 return new BNGCoordinates( CoordinateType::britishNationalGrid, BNGString );
00714 }
00715
00716
00717 MSP::CCS::MapProjectionCoordinates* BritishNationalGrid::convertToTransverseMercator( MSP::CCS::BNGCoordinates* bngCoordinates )
00718 {
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731 long northing_Offset, easting_Offset;
00732 long i = 0;
00733 long j = 0;
00734 long in_Precision;
00735
00736 char* BNGString = bngCoordinates->BNGString();
00737
00738 breakBNGString( BNGString, BNG_Letters, &BNG_Easting, &BNG_Northing, &in_Precision );
00739
00740 findIndex(BNG_Letters[0], BNG500GRID, &i);
00741
00742 northing_Offset = 500000 * (i / 2);
00743 easting_Offset = 500000 * (i % 2);
00744
00745 findIndex(BNG_Letters[1], BNG100GRID, &j);
00746
00747 northing_Offset += 100000 * (j / 5);
00748 easting_Offset += 100000 * (j % 5);
00749
00750 double easting = BNG_Easting + easting_Offset;
00751 double northing = BNG_Northing + northing_Offset;
00752
00753 return new MapProjectionCoordinates( CoordinateType::transverseMercator, easting, northing );
00754 }
00755
00756
00757
00758
00759