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 "VanDerGrinten.h"
00087 #include "MapProjection3Parameters.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 MAX_LAT = ( 90 * (PI / 180.0) );
00114 const double TWO_PI = (2.0 * PI);
00115 const double TWO_OVER_PI = (2.0 / PI);
00116 const double PI_OVER_3 = (PI / 3.0);
00117 const double ONE_THIRD = (1.0 / 3.0);
00118
00119
00120
00121
00122
00123
00124
00125
00126 VanDerGrinten::VanDerGrinten( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double falseEasting, double falseNorthing ) :
00127 CoordinateSystem(),
00128 es2( 0.0066943799901413800 ),
00129 es4( 4.4814723452405e-005 ),
00130 es6( 3.0000678794350e-007 ),
00131 Ra( 6371007.1810824 ),
00132 PI_Ra( 20015109.356056 ),
00133 Grin_Origin_Long( 0.0 ),
00134 Grin_False_Easting( 0.0 ),
00135 Grin_False_Northing( 0.0 )
00136 {
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153 double inv_f = 1 / ellipsoidFlattening;
00154 char errorStatus[500] = "";
00155
00156 if (ellipsoidSemiMajorAxis <= 0.0)
00157 {
00158 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00159 }
00160 if ((inv_f < 250) || (inv_f > 350))
00161 {
00162 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00163 }
00164 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00165 {
00166 strcat( errorStatus, ErrorMessages::centralMeridian );
00167 }
00168
00169 if( strlen( errorStatus ) > 0)
00170 throw CoordinateConversionException( errorStatus );
00171
00172 semiMajorAxis = ellipsoidSemiMajorAxis;
00173 flattening = ellipsoidFlattening;
00174
00175 es2 = 2 * flattening - flattening * flattening;
00176 es4 = es2 * es2;
00177 es6 = es4 * es2;
00178
00179 Ra = semiMajorAxis * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 /3024.0);
00180 PI_Ra = PI * Ra;
00181 if (centralMeridian > PI)
00182 centralMeridian -= TWO_PI;
00183 Grin_Origin_Long = centralMeridian;
00184 Grin_False_Easting = falseEasting;
00185 Grin_False_Northing = falseNorthing;
00186 }
00187
00188
00189 VanDerGrinten::VanDerGrinten( const VanDerGrinten &v )
00190 {
00191 semiMajorAxis = v.semiMajorAxis;
00192 flattening = v.flattening;
00193 es2 = v.es2;
00194 es4 = v.es4;
00195 es6 = v.es6;
00196 Ra = v.Ra;
00197 PI_Ra = v.PI_Ra;
00198 Grin_Origin_Long = v.Grin_Origin_Long;
00199 Grin_False_Easting = v.Grin_False_Easting;
00200 Grin_False_Northing = v.Grin_False_Northing;
00201 }
00202
00203
00204 VanDerGrinten::~VanDerGrinten()
00205 {
00206 }
00207
00208
00209 VanDerGrinten& VanDerGrinten::operator=( const VanDerGrinten &v )
00210 {
00211 if( this != &v )
00212 {
00213 semiMajorAxis = v.semiMajorAxis;
00214 flattening = v.flattening;
00215 es2 = v.es2;
00216 es4 = v.es4;
00217 es6 = v.es6;
00218 Ra = v.Ra;
00219 PI_Ra = v.PI_Ra;
00220 Grin_Origin_Long = v.Grin_Origin_Long;
00221 Grin_False_Easting = v.Grin_False_Easting;
00222 Grin_False_Northing = v.Grin_False_Northing;
00223 }
00224
00225 return *this;
00226 }
00227
00228
00229 MapProjection3Parameters* VanDerGrinten::getParameters() const
00230 {
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245 return new MapProjection3Parameters( CoordinateType::vanDerGrinten, Grin_Origin_Long, Grin_False_Easting, Grin_False_Northing );
00246 }
00247
00248
00249 MSP::CCS::MapProjectionCoordinates* VanDerGrinten::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00250 {
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264 double dlam;
00265 double aa, aasqr;
00266 double gg;
00267 double pp, ppsqr;
00268 double gg_MINUS_ppsqr, ppsqr_PLUS_aasqr;
00269 double in_theta;
00270 double theta;
00271 double sin_theta, cos_theta;
00272 double qq;
00273 double easting, northing;
00274 char errorStatus[50] = "";
00275
00276 double longitude = geodeticCoordinates->longitude();
00277 double latitude = geodeticCoordinates->latitude();
00278
00279 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00280 {
00281 strcat( errorStatus, ErrorMessages::latitude );
00282 }
00283 if ((longitude < -PI) || (longitude > TWO_PI))
00284 {
00285 strcat( errorStatus, ErrorMessages::longitude );
00286 }
00287
00288 if( strlen( errorStatus ) > 0)
00289 throw CoordinateConversionException( errorStatus );
00290
00291 dlam = longitude - Grin_Origin_Long;
00292 if (dlam > PI)
00293 {
00294 dlam -= TWO_PI;
00295 }
00296 if (dlam < -PI)
00297 {
00298 dlam += TWO_PI;
00299 }
00300
00301 if (latitude == 0.0)
00302 {
00303 easting = Ra * dlam + Grin_False_Easting;
00304 northing = 0.0;
00305 }
00306 else if (dlam == 0.0 || floatEq(latitude,MAX_LAT,.00001) || floatEq(latitude,-MAX_LAT,.00001))
00307 {
00308 in_theta = fabs(TWO_OVER_PI * latitude);
00309
00310 if (in_theta > 1.0)
00311 in_theta = 1.0;
00312 else if (in_theta < -1.0)
00313 in_theta = -1.0;
00314
00315 theta = asin(in_theta);
00316 easting = 0.0;
00317 northing = PI_Ra * tan(theta / 2) + Grin_False_Northing;
00318 if (latitude < 0.0)
00319 northing *= -1.0;
00320 }
00321 else
00322 {
00323 aa = 0.5 * fabs(PI / dlam - dlam / PI);
00324 in_theta = fabs(TWO_OVER_PI * latitude);
00325
00326 if (in_theta > 1.0)
00327 in_theta = 1.0;
00328 else if (in_theta < -1.0)
00329 in_theta = -1.0;
00330
00331 theta = asin(in_theta);
00332 sin_theta = sin(theta);
00333 cos_theta = cos(theta);
00334 gg = cos_theta / (sin_theta + cos_theta - 1);
00335 pp = gg * (2 / sin_theta - 1);
00336 aasqr = aa * aa;
00337 ppsqr = pp * pp;
00338 gg_MINUS_ppsqr = gg - ppsqr;
00339 ppsqr_PLUS_aasqr = ppsqr + aasqr;
00340 qq = aasqr + gg;
00341 easting = PI_Ra * (aa * (gg_MINUS_ppsqr) +
00342 sqrt(aasqr * (gg_MINUS_ppsqr) * (gg_MINUS_ppsqr) -
00343 (ppsqr_PLUS_aasqr) * (gg * gg - ppsqr))) /
00344 (ppsqr_PLUS_aasqr) + Grin_False_Easting;
00345 if (dlam < 0.0)
00346 easting *= -1.0;
00347 northing = PI_Ra * (pp * qq - aa * sqrt ((aasqr + 1) * (ppsqr_PLUS_aasqr) - qq * qq)) /
00348 (ppsqr_PLUS_aasqr) + Grin_False_Northing;
00349 if (latitude < 0.0)
00350 northing *= -1.0;
00351 }
00352
00353 return new MapProjectionCoordinates( CoordinateType::vanDerGrinten, easting, northing );
00354 }
00355
00356
00357 MSP::CCS::GeodeticCoordinates* VanDerGrinten::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00358 {
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372 double dx, dy;
00373 double xx, xxsqr;
00374 double yy, yysqr, two_yysqr;
00375 double xxsqr_PLUS_yysqr;
00376 double c1;
00377 double c2;
00378 double c3, c3sqr;
00379 double c2_OVER_3c3;
00380 double dd;
00381 double a1;
00382 double m1;
00383 double i;
00384 double theta1;
00385 double temp;
00386 const double epsilon = 1.0e-2;
00387 double delta = PI_Ra + epsilon;
00388 double longitude, latitude;
00389 char errorStatus[50] = "";
00390
00391 double easting = mapProjectionCoordinates->easting();
00392 double northing = mapProjectionCoordinates->northing();
00393
00394 if ((easting > (Grin_False_Easting + delta)) ||
00395 (easting < (Grin_False_Easting - delta)))
00396 {
00397 strcat( errorStatus, ErrorMessages::easting );
00398 }
00399 if ((northing > (Grin_False_Northing + delta)) ||
00400 (northing < (Grin_False_Northing - delta)))
00401 {
00402 strcat( errorStatus, ErrorMessages::northing );
00403 }
00404
00405 if( strlen( errorStatus ) > 0)
00406 throw CoordinateConversionException( errorStatus );
00407
00408 temp = sqrt(easting * easting + northing * northing);
00409
00410 if ((temp > (Grin_False_Easting + PI_Ra + epsilon)) ||
00411 (temp > (Grin_False_Northing + PI_Ra + epsilon)) ||
00412 (temp < (Grin_False_Easting - PI_Ra - epsilon)) ||
00413 (temp < (Grin_False_Northing - PI_Ra - epsilon)))
00414 {
00415 throw CoordinateConversionException( ErrorMessages::radius );
00416 }
00417
00418 dy = northing - Grin_False_Northing;
00419 dx = easting - Grin_False_Easting;
00420 xx = dx / PI_Ra;
00421 yy = dy / PI_Ra;
00422 xxsqr = xx * xx;
00423 yysqr = yy * yy;
00424 xxsqr_PLUS_yysqr = xxsqr + yysqr;
00425 two_yysqr = 2 * yysqr;
00426
00427 if (northing == 0.0)
00428 latitude = 0.0;
00429 else
00430 {
00431 c1 = - fabs(yy) * (1 + xxsqr_PLUS_yysqr);
00432 c2 = c1 - two_yysqr + xxsqr;
00433 c3 = - 2 * c1 + 1 + two_yysqr + (xxsqr_PLUS_yysqr) * (xxsqr_PLUS_yysqr);
00434 c2_OVER_3c3 = c2 / (3.0 * c3);
00435 c3sqr = c3 * c3;
00436 dd = yysqr / c3 + ((2 * c2 * c2 * c2) / (c3sqr * c3) - (9 * c1 * c2) / (c3sqr)) / 27;
00437 a1 = (c1 - c2 * c2_OVER_3c3) /c3;
00438 m1 = 2 * sqrt(-ONE_THIRD * a1);
00439 i = 3 * dd/ (a1 * m1);
00440 if ((i > 1.0)||(i < -1.0))
00441 latitude = MAX_LAT;
00442 else
00443 {
00444 theta1 = ONE_THIRD * acos(3 * dd / (a1 * m1));
00445 latitude = PI * (-m1 * cos(theta1 + PI_OVER_3) - c2_OVER_3c3);
00446 }
00447 }
00448 if (northing < 0.0)
00449 latitude *= -1.0;
00450
00451 if (xx == 0.0)
00452 longitude = Grin_Origin_Long;
00453 else
00454 {
00455 longitude = PI * (xxsqr_PLUS_yysqr - 1 +
00456 sqrt(1 + (2 * xxsqr - two_yysqr) + (xxsqr_PLUS_yysqr) * (xxsqr_PLUS_yysqr))) /
00457 (2 * xx) + Grin_Origin_Long;
00458 }
00459 if (latitude > PI_OVER_2)
00460 latitude = PI_OVER_2;
00461 else if (latitude < -PI_OVER_2)
00462 latitude = -PI_OVER_2;
00463
00464 if (longitude > PI)
00465 longitude -= TWO_PI;
00466 if (longitude < -PI)
00467 longitude += TWO_PI;
00468
00469 if (longitude > PI)
00470 longitude = PI;
00471 else if (longitude < -PI)
00472 longitude = -PI;
00473
00474 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00475 }
00476
00477
00478 double VanDerGrinten::floatEq( double x, double v, double epsilon )
00479 {
00480 return (((v - epsilon) < x) && (x < (v + epsilon)));
00481 }
00482
00483
00484