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 #include <math.h>
00088 #include "UPS.h"
00089 #include "PolarStereographic.h"
00090 #include "PolarStereographicScaleFactorParameters.h"
00091 #include "MapProjectionCoordinates.h"
00092 #include "UPSCoordinates.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 using namespace MSP::CCS;
00109
00110
00111
00112
00113
00114
00115
00116 #define EPSILON 1.75e-7
00117
00118 const double PI = 3.14159265358979323e0;
00119 const double PI_OVER = (PI/2.0e0);
00120 const double MAX_LAT = 90.0 * (PI / 180.0);
00121 const double MAX_ORIGIN_LAT = 81.114528 * (PI / 180.0);
00122 const double MIN_NORTH_LAT = 83.5 * (PI / 180.0);
00123 const double MAX_SOUTH_LAT = -79.5 * (PI / 180.0);
00124 const double MIN_EAST_NORTH = 0.0;
00125 const double MAX_EAST_NORTH = 4000000.0;
00126
00127 const double UPS_False_Easting = 2000000;
00128 const double UPS_False_Northing = 2000000;
00129 const double UPS_Origin_Longitude = 0.0;
00130
00131
00132
00133
00134
00135
00136
00137 UPS::UPS( double ellipsoidSemiMajorAxis, double ellipsoidFlattening ) :
00138 CoordinateSystem(),
00139 UPS_Origin_Latitude( MAX_ORIGIN_LAT )
00140 {
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150 double inv_f = 1 / ellipsoidFlattening;
00151 char errorStatus[500] = "";
00152
00153 if (ellipsoidSemiMajorAxis <= 0.0)
00154 {
00155 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00156 }
00157 if ((inv_f < 250) || (inv_f > 350))
00158 {
00159 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00160 }
00161
00162 if( strlen( errorStatus ) > 0)
00163 throw CoordinateConversionException( errorStatus );
00164
00165 semiMajorAxis = ellipsoidSemiMajorAxis;
00166 flattening = ellipsoidFlattening;
00167
00168 polarStereographicMap['N'] = new PolarStereographic(semiMajorAxis, flattening, UPS_Origin_Longitude, .994, 'N', UPS_False_Easting, UPS_False_Northing);
00169 polarStereographicMap['S'] = new PolarStereographic(semiMajorAxis, flattening, UPS_Origin_Longitude, .994, 'S', UPS_False_Easting, UPS_False_Northing);
00172 }
00173
00174
00175 UPS::UPS( const UPS &u )
00176 {
00177 std::map< char, PolarStereographic* > tempPolarStereographicMap = u.polarStereographicMap;
00178 polarStereographicMap['N'] = new PolarStereographic( *tempPolarStereographicMap['N'] );
00179 polarStereographicMap['S'] = new PolarStereographic( *tempPolarStereographicMap['S'] );
00180 semiMajorAxis = u.semiMajorAxis;
00181 flattening = u.flattening;
00182 UPS_Origin_Latitude = u.UPS_Origin_Latitude;
00183 }
00184
00185
00186 UPS::~UPS()
00187 {
00188 while( polarStereographicMap.begin() != polarStereographicMap.end() )
00189 {
00190 delete ( ( *polarStereographicMap.begin() ).second );
00191 polarStereographicMap.erase( polarStereographicMap.begin() );
00192 }
00193 }
00194
00195
00196 UPS& UPS::operator=( const UPS &u )
00197 {
00198 if( this != &u )
00199 {
00200 std::map< char, PolarStereographic* > tempPolarStereographicMap = u.polarStereographicMap;
00201 polarStereographicMap['N']->operator=( *tempPolarStereographicMap['N'] );
00202 polarStereographicMap['S']->operator=( *tempPolarStereographicMap['S'] );
00203 semiMajorAxis = u.semiMajorAxis;
00204 flattening = u.flattening;
00205 UPS_Origin_Latitude = u.UPS_Origin_Latitude;
00206 }
00207
00208 return *this;
00209 }
00210
00211
00212 MSP::CCS::UPSCoordinates* UPS::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00213 {
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228 char hemisphere;
00229 char errorStatus[50] = "";
00230
00231 double longitude = geodeticCoordinates->longitude();
00232 double latitude = geodeticCoordinates->latitude();
00233
00234 if ((latitude < -MAX_LAT) || (latitude > MAX_LAT))
00235 {
00236 strcat( errorStatus, ErrorMessages::latitude );
00237 }
00238 else if ((latitude < 0) && (latitude >= (MAX_SOUTH_LAT + EPSILON)))
00239 strcat( errorStatus, ErrorMessages::latitude );
00240 else if ((latitude >= 0) && (latitude < (MIN_NORTH_LAT - EPSILON)))
00241 strcat( errorStatus, ErrorMessages::latitude );
00242 if ((longitude < -PI) || (longitude > (2 * PI)))
00243 {
00244 strcat( errorStatus, ErrorMessages::longitude );
00245 }
00246
00247 if( strlen( errorStatus ) > 0)
00248 throw CoordinateConversionException( errorStatus );
00249
00250 if (latitude < 0)
00251 {
00252 UPS_Origin_Latitude = -MAX_ORIGIN_LAT;
00253 hemisphere = 'S';
00254 }
00255 else
00256 {
00257 UPS_Origin_Latitude = MAX_ORIGIN_LAT;
00258 hemisphere = 'N';
00259 }
00260
00261 PolarStereographic polarStereographic = *polarStereographicMap[hemisphere];
00262 MapProjectionCoordinates* polarStereographicCoordinates = polarStereographic.convertFromGeodetic( geodeticCoordinates );
00263
00264 double easting = polarStereographicCoordinates->easting();
00265 double northing = polarStereographicCoordinates->northing();
00266 delete polarStereographicCoordinates;
00267
00268 return new UPSCoordinates( CoordinateType::universalPolarStereographic, hemisphere, easting, northing );
00269 }
00270
00271
00272 MSP::CCS::GeodeticCoordinates* UPS::convertToGeodetic( MSP::CCS::UPSCoordinates* upsCoordinates )
00273 {
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287 char errorStatus[50] = "";
00288
00289 char hemisphere = upsCoordinates->hemisphere();
00290 double easting = upsCoordinates->easting();
00291 double northing = upsCoordinates->northing();
00292
00293 if ((hemisphere != 'N') && (hemisphere != 'S'))
00294 strcat( errorStatus, ErrorMessages::hemisphere );
00295 if ((easting < MIN_EAST_NORTH) || (easting > MAX_EAST_NORTH))
00296 strcat( errorStatus, ErrorMessages::easting );
00297 if ((northing < MIN_EAST_NORTH) || (northing > MAX_EAST_NORTH))
00298 strcat( errorStatus, ErrorMessages::northing );
00299
00300 if (hemisphere =='N')
00301 {
00302 UPS_Origin_Latitude = MAX_ORIGIN_LAT;
00303 }
00304 else if (hemisphere =='S')
00305 {
00306 UPS_Origin_Latitude = -MAX_ORIGIN_LAT;
00307 }
00308
00309 if( strlen( errorStatus ) > 0)
00310 throw CoordinateConversionException( errorStatus );
00311
00312 MapProjectionCoordinates polarStereographicCoordinates( CoordinateType::polarStereographicStandardParallel, easting, northing );
00313 PolarStereographic polarStereographic = *polarStereographicMap[hemisphere];
00314 GeodeticCoordinates* geodeticCoordinates = polarStereographic.convertToGeodetic( &polarStereographicCoordinates );
00315
00316 double latitude = geodeticCoordinates->latitude();
00317
00318 if ((latitude < 0) && (latitude >= (MAX_SOUTH_LAT + EPSILON)))
00319 {
00320 delete geodeticCoordinates;
00321 throw CoordinateConversionException( ErrorMessages::latitude );
00322 }
00323 if ((latitude >= 0) && (latitude < (MIN_NORTH_LAT - EPSILON)))
00324 {
00325 delete geodeticCoordinates;
00326 throw CoordinateConversionException( ErrorMessages::latitude );
00327 }
00328
00329 return geodeticCoordinates;
00330 }
00331
00332