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 "MillerCylindrical.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 TWO_PI = (2.0 * PI);
00112
00113
00114
00115
00116
00117
00118
00119 MillerCylindrical::MillerCylindrical( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double falseEasting, double falseNorthing ) :
00120 CoordinateSystem(),
00121 es2( 0.0066943799901413800 ),
00122 es4( 4.4814723452405e-005 ),
00123 es6( 3.0000678794350e-007 ),
00124 Ra( 6371007.1810824 ),
00125 Mill_Origin_Long( 0.0 ),
00126 Mill_False_Easting( 0.0 ),
00127 Mill_False_Northing( 0.0 ),
00128 Mill_Delta_Northing( 14675058.0 ),
00129 Mill_Max_Easting( 20015110.0 ),
00130 Mill_Min_Easting( -20015110.0 )
00131 {
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148 double inv_f = 1 / ellipsoidFlattening;
00149 char errorStatus[500] = "";
00150
00151 if (ellipsoidSemiMajorAxis <= 0.0)
00152 {
00153 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00154 }
00155 if ((inv_f < 250) || (inv_f > 350))
00156 {
00157 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00158 }
00159 if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
00160 {
00161 strcat( errorStatus, ErrorMessages::centralMeridian );
00162 }
00163
00164 if( strlen( errorStatus ) > 0)
00165 throw CoordinateConversionException( errorStatus );
00166
00167 semiMajorAxis = ellipsoidSemiMajorAxis;
00168 flattening = ellipsoidFlattening;
00169
00170 es2 = 2 * flattening - flattening * flattening;
00171 es4 = es2 * es2;
00172 es6 = es4 * es2;
00173
00174 Ra = semiMajorAxis * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 /3024.0);
00175 if (centralMeridian > PI)
00176 centralMeridian -= TWO_PI;
00177 Mill_Origin_Long = centralMeridian;
00178 Mill_False_Easting = falseEasting;
00179 Mill_False_Northing = falseNorthing;
00180 if (Mill_Origin_Long > 0)
00181 {
00182 Mill_Max_Easting = 19903915.0;
00183 Mill_Min_Easting = -20015110.0;
00184 }
00185 else if (Mill_Origin_Long < 0)
00186 {
00187 Mill_Max_Easting = 20015110.0;
00188 Mill_Min_Easting = -19903915.0;
00189 }
00190 else
00191 {
00192 Mill_Max_Easting = 20015110.0;
00193 Mill_Min_Easting = -20015110.0;
00194 }
00195 }
00196
00197
00198 MillerCylindrical::MillerCylindrical( const MillerCylindrical &mc )
00199 {
00200 semiMajorAxis = mc.semiMajorAxis;
00201 flattening = mc.flattening;
00202 es2 = mc.es2;
00203 es4 = mc.es4;
00204 es6 = mc.es6;
00205 Ra = mc.Ra;
00206 Mill_Origin_Long = mc.Mill_Origin_Long;
00207 Mill_False_Easting = mc.Mill_False_Easting;
00208 Mill_False_Northing = mc.Mill_False_Northing;
00209 Mill_Delta_Northing = mc.Mill_Delta_Northing;
00210 Mill_Max_Easting = mc.Mill_Max_Easting;
00211 Mill_Min_Easting = mc.Mill_Min_Easting;
00212 }
00213
00214
00215 MillerCylindrical::~MillerCylindrical()
00216 {
00217 }
00218
00219
00220 MillerCylindrical& MillerCylindrical::operator=( const MillerCylindrical &mc )
00221 {
00222 if( this != &mc )
00223 {
00224 semiMajorAxis = mc.semiMajorAxis;
00225 flattening = mc.flattening;
00226 es2 = mc.es2;
00227 es4 = mc.es4;
00228 es6 = mc.es6;
00229 Ra = mc.Ra;
00230 Mill_Origin_Long = mc.Mill_Origin_Long;
00231 Mill_False_Easting = mc.Mill_False_Easting;
00232 Mill_False_Northing = mc.Mill_False_Northing;
00233 Mill_Delta_Northing = mc.Mill_Delta_Northing;
00234 Mill_Max_Easting = mc.Mill_Max_Easting;
00235 Mill_Min_Easting = mc.Mill_Min_Easting;
00236 }
00237
00238 return *this;
00239 }
00240
00241
00242 MapProjection3Parameters* MillerCylindrical::getParameters() const
00243 {
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258 return new MapProjection3Parameters( CoordinateType::millerCylindrical, Mill_Origin_Long, Mill_False_Easting, Mill_False_Northing );
00259 }
00260
00261
00262 MSP::CCS::MapProjectionCoordinates* MillerCylindrical::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00263 {
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277 double dlam;
00278 char errorStatus[50] = "";
00279
00280 double longitude = geodeticCoordinates->longitude();
00281 double latitude = geodeticCoordinates->latitude();
00282 double slat = sin(0.8 * latitude);
00283
00284 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00285 {
00286 strcat( errorStatus, ErrorMessages::latitude );
00287 }
00288 if ((longitude < -PI) || (longitude > TWO_PI))
00289 {
00290 strcat( errorStatus, ErrorMessages::longitude );
00291 }
00292
00293 if( strlen( errorStatus ) > 0)
00294 throw CoordinateConversionException( errorStatus );
00295
00296 dlam = longitude - Mill_Origin_Long;
00297 if (dlam > PI)
00298 {
00299 dlam -= TWO_PI;
00300 }
00301 if (dlam < -PI)
00302 {
00303 dlam += TWO_PI;
00304 }
00305 double easting = Ra * dlam + Mill_False_Easting;
00306 double northing = (Ra / 1.6) * log((1.0 + slat) /
00307 (1.0 - slat)) + Mill_False_Northing;
00308
00309 return new MapProjectionCoordinates( CoordinateType::millerCylindrical, easting, northing );
00310 }
00311
00312
00313 MSP::CCS::GeodeticCoordinates* MillerCylindrical::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00314 {
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328 double dx, dy;
00329 char errorStatus[50] = "";
00330
00331 double easting = mapProjectionCoordinates->easting();
00332 double northing = mapProjectionCoordinates->northing();
00333
00334 if ((easting < (Mill_False_Easting + Mill_Min_Easting))
00335 || (easting > (Mill_False_Easting + Mill_Max_Easting)))
00336 {
00337 strcat( errorStatus, ErrorMessages::easting );
00338 }
00339 if ((northing < (Mill_False_Northing - Mill_Delta_Northing)) ||
00340 (northing > (Mill_False_Northing + Mill_Delta_Northing) ))
00341 {
00342 strcat( errorStatus, ErrorMessages::northing );
00343 }
00344
00345 if( strlen( errorStatus ) > 0)
00346 throw CoordinateConversionException( errorStatus );
00347
00348 dy = northing - Mill_False_Northing;
00349 dx = easting - Mill_False_Easting;
00350 double latitude = atan(sinh(0.8 * dy / Ra)) / 0.8;
00351 double longitude = Mill_Origin_Long + dx / Ra;
00352
00353 if (latitude > PI_OVER_2)
00354 latitude = PI_OVER_2;
00355 else if (latitude < -PI_OVER_2)
00356 latitude = -PI_OVER_2;
00357
00358 if (longitude > PI)
00359 longitude -= TWO_PI;
00360 if (longitude < -PI)
00361 longitude += TWO_PI;
00362
00363 if (longitude > PI)
00364 longitude = PI;
00365 else if (longitude < -PI)
00366 longitude = -PI;
00367
00368 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00369 }
00370
00371
00372
00373