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 "Orthographic.h"
00090 #include "MapProjection4Parameters.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 MAX_LAT = ( (PI * 90) / 180.0 );
00117 const double TWO_PI = (2.0 * PI);
00118
00119
00120
00121
00122
00123
00124
00125 Orthographic::Orthographic( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double originLatitude, 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 Orth_Origin_Long( 0.0 ),
00132 Orth_Origin_Lat( 0.0 ),
00133 Orth_False_Easting( 0.0 ),
00134 Orth_False_Northing( 0.0 ),
00135 Sin_Orth_Origin_Lat( 0.0 ),
00136 Cos_Orth_Origin_Lat( 1.0 )
00137 {
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156 double inv_f = 1 / ellipsoidFlattening;
00157 char errorStatus[500] = "";
00158
00159 if (ellipsoidSemiMajorAxis <= 0.0)
00160 {
00161 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00162 }
00163 if ((inv_f < 250) || (inv_f > 350))
00164 {
00165 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00166 }
00167 if ((originLatitude < -PI_OVER_2) || (originLatitude > PI_OVER_2))
00168 {
00169 strcat( errorStatus, ErrorMessages::originLatitude );
00170 }
00171 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00172 {
00173 strcat( errorStatus, ErrorMessages::centralMeridian );
00174 }
00175
00176 if( strlen( errorStatus ) > 0)
00177 throw CoordinateConversionException( errorStatus );
00178
00179 semiMajorAxis = ellipsoidSemiMajorAxis;
00180 flattening = ellipsoidFlattening;
00181
00182 es2 = 2 * flattening - flattening * flattening;
00183 es4 = es2 * es2;
00184 es6 = es4 * es2;
00185 Ra = semiMajorAxis * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 /3024.0);
00186 Orth_Origin_Lat = originLatitude;
00187 Sin_Orth_Origin_Lat = sin(Orth_Origin_Lat);
00188 Cos_Orth_Origin_Lat = cos(Orth_Origin_Lat);
00189 if (centralMeridian > PI)
00190 centralMeridian -= TWO_PI;
00191 Orth_Origin_Long = centralMeridian;
00192 Orth_False_Easting = falseEasting;
00193 Orth_False_Northing = falseNorthing;
00194 }
00195
00196
00197 Orthographic::Orthographic( const Orthographic &o )
00198 {
00199 semiMajorAxis = o.semiMajorAxis;
00200 flattening = o.flattening;
00201 es2 = o.es2;
00202 es4 = o.es4;
00203 es6 = o.es6;
00204 Ra = o.Ra;
00205 Orth_Origin_Long = o.Orth_Origin_Long;
00206 Orth_Origin_Lat = o.Orth_Origin_Lat;
00207 Orth_False_Easting = o.Orth_False_Easting;
00208 Orth_False_Northing = o.Orth_False_Northing;
00209 Sin_Orth_Origin_Lat = o.Sin_Orth_Origin_Lat;
00210 Cos_Orth_Origin_Lat = o.Cos_Orth_Origin_Lat;
00211 }
00212
00213
00214 Orthographic::~Orthographic()
00215 {
00216 }
00217
00218
00219 Orthographic& Orthographic::operator=( const Orthographic &o )
00220 {
00221 if( this != &o )
00222 {
00223 semiMajorAxis = o.semiMajorAxis;
00224 flattening = o.flattening;
00225 es2 = o.es2;
00226 es4 = o.es4;
00227 es6 = o.es6;
00228 Ra = o.Ra;
00229 Orth_Origin_Long = o.Orth_Origin_Long;
00230 Orth_Origin_Lat = o.Orth_Origin_Lat;
00231 Orth_False_Easting = o.Orth_False_Easting;
00232 Orth_False_Northing = o.Orth_False_Northing;
00233 Sin_Orth_Origin_Lat = o.Sin_Orth_Origin_Lat;
00234 Cos_Orth_Origin_Lat = o.Cos_Orth_Origin_Lat;
00235 }
00236
00237 return *this;
00238 }
00239
00240
00241 MapProjection4Parameters* Orthographic::getParameters() const
00242 {
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259 return new MapProjection4Parameters( CoordinateType::orthographic, Orth_Origin_Long, Orth_Origin_Lat, Orth_False_Easting, Orth_False_Northing );
00260 }
00261
00262
00263 MSP::CCS::MapProjectionCoordinates* Orthographic::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00264 {
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278 double dlam;
00279 double clat_cdlam;
00280 double cos_c;
00281
00282 char errorStatus[50] = "";
00283
00284 double longitude = geodeticCoordinates->longitude();
00285 double latitude = geodeticCoordinates->latitude();
00286 double slat = sin(latitude);
00287 double clat = cos(latitude);
00288
00289 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00290 {
00291 strcat( errorStatus, ErrorMessages::latitude );
00292 }
00293 if ((longitude < -PI) || (longitude > TWO_PI))
00294 {
00295 strcat( errorStatus, ErrorMessages::longitude );
00296 }
00297
00298 if( strlen( errorStatus ) > 0)
00299 throw CoordinateConversionException( errorStatus );
00300
00301 dlam = longitude - Orth_Origin_Long;
00302 clat_cdlam = clat * cos(dlam);
00303 cos_c = Sin_Orth_Origin_Lat * slat + Cos_Orth_Origin_Lat * clat_cdlam;
00304 if (cos_c < 0.0)
00305 {
00306
00307 throw CoordinateConversionException( ErrorMessages::longitude );
00308 }
00309
00310 if (dlam > PI)
00311 {
00312 dlam -= TWO_PI;
00313 }
00314 if (dlam < -PI)
00315 {
00316 dlam += TWO_PI;
00317 }
00318 double easting = Ra * clat * sin(dlam) + Orth_False_Easting;
00319 double northing = Ra * (Cos_Orth_Origin_Lat * slat - Sin_Orth_Origin_Lat * clat_cdlam) +
00320 Orth_False_Northing;
00321
00322 return new MapProjectionCoordinates( CoordinateType::orthographic, easting, northing );
00323 }
00324
00325
00326 MSP::CCS::GeodeticCoordinates* Orthographic::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00327 {
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341 double cc;
00342 double cos_cc, sin_cc;
00343 double rho;
00344 double dx, dy;
00345 double temp;
00346 double rho_OVER_Ra;
00347 double longitude, latitude;
00348 char errorStatus[50] = "";
00349
00350 double easting = mapProjectionCoordinates->easting();
00351 double northing = mapProjectionCoordinates->northing();
00352
00353 if ((easting > (Orth_False_Easting + Ra)) ||
00354 (easting < (Orth_False_Easting - Ra)))
00355 {
00356 strcat( errorStatus, ErrorMessages::easting );
00357 }
00358 if ((northing > (Orth_False_Northing + Ra)) ||
00359 (northing < (Orth_False_Northing - Ra)))
00360 {
00361 strcat( errorStatus, ErrorMessages::northing );
00362 }
00363
00364 if( strlen( errorStatus ) > 0)
00365 throw CoordinateConversionException( errorStatus );
00366
00367 temp = sqrt(easting * easting + northing * northing);
00368
00369 if ((temp > (Orth_False_Easting + Ra)) || (temp > (Orth_False_Northing + Ra)) ||
00370 (temp < (Orth_False_Easting - Ra)) || (temp < (Orth_False_Northing - Ra)))
00371 {
00372 throw CoordinateConversionException( ErrorMessages::radius );
00373 }
00374
00375 dx = easting - Orth_False_Easting;
00376 dy = northing - Orth_False_Northing;
00377 rho = sqrt(dx * dx + dy * dy);
00378 if (rho == 0.0)
00379 {
00380 latitude = Orth_Origin_Lat;
00381 longitude = Orth_Origin_Long;
00382 }
00383 else
00384 {
00385 rho_OVER_Ra = rho / Ra;
00386
00387 if (rho_OVER_Ra > 1.0)
00388 rho_OVER_Ra = 1.0;
00389 else if (rho_OVER_Ra < -1.0)
00390 rho_OVER_Ra = -1.0;
00391
00392 cc = asin(rho_OVER_Ra);
00393 cos_cc = cos(cc);
00394 sin_cc = sin(cc);
00395 latitude = asin(cos_cc * Sin_Orth_Origin_Lat + (dy * sin_cc * Cos_Orth_Origin_Lat / rho));
00396
00397 if (Orth_Origin_Lat == MAX_LAT)
00398 longitude = Orth_Origin_Long + atan2(dx, -dy);
00399 else if (Orth_Origin_Lat == -MAX_LAT)
00400 longitude = Orth_Origin_Long + atan2(dx, dy);
00401 else
00402 longitude = Orth_Origin_Long + atan2(dx * sin_cc, (rho *
00403 Cos_Orth_Origin_Lat * cos_cc - dy * Sin_Orth_Origin_Lat * sin_cc));
00404
00405 if (latitude > PI_OVER_2)
00406 latitude = PI_OVER_2;
00407 else if (latitude < -PI_OVER_2)
00408 latitude = -PI_OVER_2;
00409
00410 if (longitude > PI)
00411 longitude -= TWO_PI;
00412 if (longitude < -PI)
00413 longitude += TWO_PI;
00414
00415 if (longitude > PI)
00416 longitude = PI;
00417 else if (longitude < -PI)
00418 longitude = -PI;
00419 }
00420
00421 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00422 }
00423
00424
00425
00426