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 #include <math.h>
00084 #include "Mollweide.h"
00085 #include "MapProjection3Parameters.h"
00086 #include "MapProjectionCoordinates.h"
00087 #include "GeodeticCoordinates.h"
00088 #include "CoordinateConversionException.h"
00089 #include "ErrorMessages.h"
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101 using namespace MSP::CCS;
00102
00103
00104
00105
00106
00107
00108
00109 const double PI = 3.14159265358979323e0;
00110 const double PI_OVER_2 = ( PI / 2.0);
00111 const double MAX_LAT = ( (PI * 90) / 180.0 );
00112 const double TWO_PI = (2.0 * PI);
00113
00114
00115
00116
00117
00118
00119
00120 Mollweide::Mollweide( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double falseEasting, double falseNorthing ) :
00121 CoordinateSystem(),
00122 es2( 0.0066943799901413800 ),
00123 es4( 4.4814723452405e-005 ),
00124 es6( 3.0000678794350e-007 ),
00125 Sqrt2_Ra( 9009964.7614632 ),
00126 Sqrt8_Ra( 18019929.522926 ),
00127 Moll_Origin_Long( 0.0 ),
00128 Moll_False_Easting( 0.0 ),
00129 Moll_False_Northing( 0.0 ),
00130 Moll_Delta_Northing( 9009965.0 ),
00131 Moll_Max_Easting( 18019930.0 ),
00132 Moll_Min_Easting( -18019930.0 )
00133 {
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150 double Ra;
00151 double inv_f = 1 / ellipsoidFlattening;
00152 char errorStatus[500] = "";
00153
00154 if (ellipsoidSemiMajorAxis <= 0.0)
00155 {
00156 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00157 }
00158 if ((inv_f < 250) || (inv_f > 350))
00159 {
00160 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00161 }
00162 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00163 {
00164 strcat( errorStatus, ErrorMessages::centralMeridian );
00165 }
00166
00167 if( strlen( errorStatus ) > 0)
00168 throw CoordinateConversionException( errorStatus );
00169
00170 semiMajorAxis = ellipsoidSemiMajorAxis;
00171 flattening = ellipsoidFlattening;
00172
00173 es2 = 2 * flattening - flattening * flattening;
00174 es4 = es2 * es2;
00175 es6 = es4 * es2;
00176
00177 Ra = semiMajorAxis * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 / 3024.0);
00178 Sqrt2_Ra = sqrt(2.0) * Ra;
00179 Sqrt8_Ra = sqrt(8.0) * Ra;
00180 if (centralMeridian > PI)
00181 centralMeridian -= TWO_PI;
00182 Moll_Origin_Long = centralMeridian;
00183 Moll_False_Easting = falseEasting;
00184 Moll_False_Northing = falseNorthing;
00185
00186 if (Moll_Origin_Long > 0)
00187 {
00188 Moll_Max_Easting = 17919819.0;
00189 Moll_Min_Easting = -18019930.0;
00190 }
00191 else if (Moll_Origin_Long < 0)
00192 {
00193 Moll_Max_Easting = 18019930.0;
00194 Moll_Min_Easting = -17919819.0;
00195 }
00196 else
00197 {
00198 Moll_Max_Easting = 18019930.0;
00199 Moll_Min_Easting = -18019930.0;
00200 }
00201 }
00202
00203
00204 Mollweide::Mollweide( const Mollweide &m )
00205 {
00206 semiMajorAxis = m.semiMajorAxis;
00207 flattening = m.flattening;
00208 es2 = m.es2;
00209 es4 = m.es4;
00210 es6 = m.es6;
00211 Sqrt2_Ra = m.Sqrt2_Ra;
00212 Sqrt8_Ra = m.Sqrt8_Ra;
00213 Moll_Origin_Long = m.Moll_Origin_Long;
00214 Moll_False_Easting = m.Moll_False_Easting;
00215 Moll_False_Northing = m.Moll_False_Northing;
00216 Moll_Delta_Northing = m.Moll_Delta_Northing;
00217 Moll_Max_Easting = m.Moll_Max_Easting;
00218 Moll_Min_Easting = m.Moll_Min_Easting;
00219 }
00220
00221
00222 Mollweide::~Mollweide()
00223 {
00224 }
00225
00226
00227 Mollweide& Mollweide::operator=( const Mollweide &m )
00228 {
00229 if( this != &m )
00230 {
00231 semiMajorAxis = m.semiMajorAxis;
00232 flattening = m.flattening;
00233 es2 = m.es2;
00234 es4 = m.es4;
00235 es6 = m.es6;
00236 Sqrt2_Ra = m.Sqrt2_Ra;
00237 Sqrt8_Ra = m.Sqrt8_Ra;
00238 Moll_Origin_Long = m.Moll_Origin_Long;
00239 Moll_False_Easting = m.Moll_False_Easting;
00240 Moll_False_Northing = m.Moll_False_Northing;
00241 Moll_Delta_Northing = m.Moll_Delta_Northing;
00242 Moll_Max_Easting = m.Moll_Max_Easting;
00243 Moll_Min_Easting = m.Moll_Min_Easting;
00244 }
00245
00246 return *this;
00247 }
00248
00249
00250 MapProjection3Parameters* Mollweide::getParameters() const
00251 {
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266 return new MapProjection3Parameters( CoordinateType::mollweide, Moll_Origin_Long, Moll_False_Easting, Moll_False_Northing );
00267 }
00268
00269
00270 MSP::CCS::MapProjectionCoordinates* Mollweide::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00271 {
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285 double dlam;
00286 double theta;
00287 double delta_theta_primed = 0.1745329;
00288 double dtp_tolerance = 4.85e-10;
00289
00290 int count = 60;
00291 char errorStatus[50] = "";
00292
00293 double longitude = geodeticCoordinates->longitude();
00294 double latitude = geodeticCoordinates->latitude();
00295 double PI_Sin_Latitude = PI * sin(latitude);
00296 double theta_primed = latitude;
00297
00298 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00299 {
00300 strcat( errorStatus, ErrorMessages::latitude );
00301 }
00302 if ((longitude < -PI) || (longitude > TWO_PI))
00303 {
00304 strcat( errorStatus, ErrorMessages::longitude );
00305 }
00306
00307 if( strlen( errorStatus ) > 0)
00308 throw CoordinateConversionException( errorStatus );
00309
00310 dlam = longitude - Moll_Origin_Long;
00311 if (dlam > PI)
00312 {
00313 dlam -= TWO_PI;
00314 }
00315 if (dlam < -PI)
00316 {
00317 dlam += TWO_PI;
00318 }
00319 while (fabs(delta_theta_primed) > dtp_tolerance && count)
00320 {
00321 delta_theta_primed = -(theta_primed + sin(theta_primed) -
00322 PI_Sin_Latitude) / (1.0 + cos(theta_primed));
00323 theta_primed += delta_theta_primed;
00324 count --;
00325 }
00326
00327 if(!count)
00328 throw CoordinateConversionException( ErrorMessages::northing );
00329
00330 theta = theta_primed / 2.0;
00331 double easting = (Sqrt8_Ra / PI ) * dlam * cos(theta) +
00332 Moll_False_Easting;
00333 double northing = Sqrt2_Ra * sin(theta) + Moll_False_Northing;
00334
00335 return new MapProjectionCoordinates( CoordinateType::mollweide, easting, northing );
00336 }
00337
00338
00339 MSP::CCS::GeodeticCoordinates* Mollweide::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00340 {
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354 double dx, dy;
00355 double theta = 0.0;
00356 double two_theta;
00357 double i;
00358 double longitude, latitude;
00359 char errorStatus[50] = "";
00360
00361 double easting = mapProjectionCoordinates->easting();
00362 double northing = mapProjectionCoordinates->northing();
00363
00364 if ((easting < (Moll_False_Easting + Moll_Min_Easting))
00365 || (easting > (Moll_False_Easting + Moll_Max_Easting)))
00366 {
00367 strcat( errorStatus, ErrorMessages::easting );
00368 }
00369 if ((northing < (Moll_False_Northing - Moll_Delta_Northing)) ||
00370 (northing >(Moll_False_Northing + Moll_Delta_Northing) ))
00371 {
00372 strcat( errorStatus, ErrorMessages::northing );
00373 }
00374
00375 if( strlen( errorStatus ) > 0)
00376 throw CoordinateConversionException( errorStatus );
00377
00378 dy = northing - Moll_False_Northing;
00379 dx = easting - Moll_False_Easting;
00380 i = dy / Sqrt2_Ra;
00381 if (fabs(i) > 1.0)
00382 {
00383 latitude = MAX_LAT;
00384 if (northing < 0.0)
00385 latitude *= -1.0;
00386 }
00387
00388 else
00389 {
00390 theta = asin(i);
00391 two_theta = 2.0 * theta;
00392 latitude = asin((two_theta + sin(two_theta)) / PI);
00393
00394 if (latitude > PI_OVER_2)
00395 latitude = PI_OVER_2;
00396 else if (latitude < -PI_OVER_2)
00397 latitude = -PI_OVER_2;
00398
00399 }
00400 if (fabs(fabs(latitude) - MAX_LAT) < 1.0e-10)
00401 longitude = Moll_Origin_Long;
00402 else
00403 longitude = Moll_Origin_Long + PI * dx /
00404 (Sqrt8_Ra * cos(theta));
00405
00406 if (longitude > PI)
00407 longitude -= TWO_PI;
00408 if (longitude < -PI)
00409 longitude += TWO_PI;
00410
00411 if (longitude > PI)
00412 longitude = PI;
00413 else if (longitude < -PI)
00414 longitude = -PI;
00415
00416 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00417 }
00418
00419
00420
00421