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 #include <string.h>
00081 #include <math.h>
00082 #include "NZMG.h"
00083 #include "EllipsoidParameters.h"
00084 #include "MapProjectionCoordinates.h"
00085 #include "GeodeticCoordinates.h"
00086 #include "CoordinateConversionException.h"
00087 #include "ErrorMessages.h"
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100 using namespace MSP::CCS;
00101
00102
00103
00104
00105
00106
00107
00108 const double PI = 3.14159265358979323e0;
00109 const double PI_OVER_2 = (PI / 2.0);
00110 const double TWO_PI = (2.0 * PI);
00111 const double MAX_LAT = (-33.5 * PI / 180.0);
00112 const double MIN_LAT = (-48.5 * PI / 180.0);
00113 const double MAX_LON = (180.0 * PI / 180.0);
00114 const double MIN_LON = (165.5 * PI / 180.0);
00115
00116 const char* International = "IN";
00117
00118
00119 const double NZMG_Origin_Lat = (-41.0);
00120 const double NZMG_Origin_Long = (173.0 * PI / 180.0);
00121 const double NZMG_False_Northing = 6023150.0;
00122 const double NZMG_False_Easting = 2510000.0;
00123
00124
00125 const double NZMG_Max_Easting = 3170000.0;
00126 const double NZMG_Max_Northing = 6900000.0;
00127 const double NZMG_Min_Easting = 1810000.0;
00128 const double NZMG_Min_Northing = 5160000.0;
00129
00130 struct Complex
00131 {
00132 double real;
00133 double imag;
00134 };
00135
00136 double A[] = { 0.6399175073, -0.1358797613, 0.063294409,
00137 -0.02526853, 0.0117879, -0.0055161,
00138 0.0026906, -0.001333, 0.00067, -0.00034 };
00139
00140 Complex B[] = { { 0.7557853228, 0.0 },
00141 { 0.249204646, 0.003371507 },
00142 { -0.001541739, 0.041058560 },
00143 { -0.10162907, 0.01727609 },
00144 { -0.26623489, -0.36249218 },
00145 { -0.6870983, -1.1651967 } };
00146
00147 Complex C[] = { { 1.3231270439, 0.0 },
00148 { -0.577245789, -0.007809598 },
00149 { 0.508307513, -0.112208952 },
00150 { -0.15094762, 0.18200602 },
00151 { 1.01418179, 1.64497696 },
00152 { 1.9660549, 2.5127645 } };
00153
00154 double D[] = { 1.5627014243, 0.5185406398, -0.03333098,
00155 -0.1052906, -0.0368594, 0.007317,
00156 0.01220, 0.00394, -0.0013 };
00157
00158
00159
00160
00161
00162
00163
00164 Complex add(Complex z1, Complex z2)
00165 {
00166 Complex z;
00167
00168 z.real = z1.real + z2.real;
00169 z.imag = z1.imag + z2.imag;
00170
00171 return z;
00172 }
00173
00174
00175
00176 Complex multiply(Complex z1, Complex z2)
00177 {
00178 Complex z;
00179
00180 z.real = z1.real * z2.real - z1.imag * z2.imag;
00181 z.imag = z1.imag * z2.real + z1.real * z2.imag;
00182
00183 return z;
00184 }
00185
00186
00187
00188 Complex divide(Complex z1, Complex z2)
00189 {
00190 Complex z;
00191 double denom;
00192
00193 denom = z2.real * z2.real + z2.imag * z2.imag;
00194 z.real = (z1.real * z2.real + z1.imag * z2.imag) / denom;
00195 z.imag = (z1.imag * z2.real - z1.real * z2.imag) / denom;
00196
00197 return z;
00198 }
00199
00200
00201 NZMG::NZMG( char* ellipsoidCode ) :
00202 CoordinateSystem( 6378388.0, 1 / 297.0 )
00203 {
00204
00205
00206
00207
00208
00209
00210
00211
00212 char errorStatus[500] = "";
00213
00214 strcpy( NZMGEllipsoidCode, "IN" );
00215
00216 if (strcmp(ellipsoidCode, International) != 0)
00217 {
00218 strcat( errorStatus, ErrorMessages::nzmgEllipsoid );
00219 }
00220
00221 if( strlen( errorStatus ) > 0)
00222 throw CoordinateConversionException( errorStatus );
00223
00224 strcpy( NZMGEllipsoidCode, ellipsoidCode );
00225 }
00226
00227
00228 NZMG::NZMG( const NZMG &n )
00229 {
00230 semiMajorAxis = n.semiMajorAxis;
00231 flattening = n.flattening;
00232 }
00233
00234
00235 NZMG::~NZMG()
00236 {
00237 }
00238
00239
00240 NZMG& NZMG::operator=( const NZMG &n )
00241 {
00242 if( this != &n )
00243 {
00244 semiMajorAxis = n.semiMajorAxis;
00245 flattening = n.flattening;
00246 }
00247
00248 return *this;
00249 }
00250
00251
00252 EllipsoidParameters* NZMG::getParameters() const
00253 {
00254
00255
00256
00257
00258
00259
00260
00261 return new EllipsoidParameters( semiMajorAxis, flattening, (char*)NZMGEllipsoidCode );
00262 }
00263
00264
00265 MSP::CCS::MapProjectionCoordinates* NZMG::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
00266 {
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280 Complex Zeta, z;
00281 int n;
00282 double dphi;
00283 double du, dlam;
00284 char errorStatus[50] = "";
00285
00286 double longitude = geodeticCoordinates->longitude();
00287 double latitude = geodeticCoordinates->latitude();
00288
00289 if ((latitude < MIN_LAT) || (latitude > MAX_LAT))
00290 {
00291 strcat( errorStatus, ErrorMessages::latitude );
00292 }
00293 if ((longitude < MIN_LON) || (longitude > MAX_LON))
00294 {
00295 strcat( errorStatus, ErrorMessages::longitude );
00296 }
00297
00298 if( strlen( errorStatus ) > 0)
00299 throw CoordinateConversionException( errorStatus );
00300
00301 dphi = (latitude * (180.0 / PI) - NZMG_Origin_Lat) * 3600.0 * 1.0e-5;
00302 du = A[9];
00303 for (n = 8; n >= 0; n--)
00304 du = du * dphi + A[n];
00305 du *= dphi;
00306
00307 dlam = longitude - NZMG_Origin_Long;
00308
00309 Zeta.real = du;
00310 Zeta.imag = dlam;
00311
00312 z.real = B[5].real;
00313 z.imag = B[5].imag;
00314 for (n = 4; n >= 0; n--)
00315 {
00316 z = multiply(z, Zeta);
00317 z = add(B[n], z);
00318 }
00319 z = multiply(z, Zeta);
00320
00321 double easting = (z.imag * semiMajorAxis) + NZMG_False_Easting;
00322 double northing = (z.real * semiMajorAxis) + NZMG_False_Northing;
00323
00324 if ((easting < NZMG_Min_Easting) || (easting > NZMG_Max_Easting))
00325 throw CoordinateConversionException( ErrorMessages::easting );
00326 if ((northing < NZMG_Min_Northing) || (northing > NZMG_Max_Northing))
00327 throw CoordinateConversionException( ErrorMessages::northing );
00328
00329 return new MapProjectionCoordinates( CoordinateType::newZealandMapGrid, easting, northing );
00330 }
00331
00332
00333 MSP::CCS::GeodeticCoordinates* NZMG::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
00334 {
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348 int i, n;
00349 Complex coeff;
00350 Complex z, Zeta, Zeta_Numer, Zeta_Denom, Zeta_sqr;
00351 double dphi;
00352 char errorStatus[50] = "";
00353
00354 double easting = mapProjectionCoordinates->easting();
00355 double northing = mapProjectionCoordinates->northing();
00356
00357 if ((easting < NZMG_Min_Easting) || (easting > NZMG_Max_Easting))
00358 {
00359 strcat( errorStatus, ErrorMessages::easting );
00360 }
00361 if ((northing < NZMG_Min_Northing) || (northing > NZMG_Max_Northing))
00362 {
00363 strcat( errorStatus, ErrorMessages::northing );
00364 }
00365
00366 if( strlen( errorStatus ) > 0)
00367 throw CoordinateConversionException( errorStatus );
00368
00369 z.real = (northing - NZMG_False_Northing) / semiMajorAxis;
00370 z.imag = (easting - NZMG_False_Easting) / semiMajorAxis;
00371
00372 Zeta.real = C[5].real;
00373 Zeta.imag = C[5].imag;
00374 for (n = 4; n >= 0; n--)
00375 {
00376 Zeta = multiply(Zeta, z);
00377 Zeta = add(C[n], Zeta);
00378 }
00379 Zeta = multiply(Zeta, z);
00380
00381 for (i = 0; i < 2; i++)
00382 {
00383 Zeta_Numer.real = 5.0 * B[5].real;
00384 Zeta_Numer.imag = 5.0 * B[5].imag;
00385 Zeta_Denom.real = 6.0 * B[5].real;
00386 Zeta_Denom.imag = 6.0 * B[5].imag;
00387 for (n = 4; n >= 1; n--)
00388 {
00389 Zeta_Numer = multiply(Zeta_Numer, Zeta);
00390 coeff.real = n * B[n].real;
00391 coeff.imag = n * B[n].imag;
00392 Zeta_Numer = add(coeff, Zeta_Numer);
00393
00394 Zeta_Denom = multiply(Zeta_Denom, Zeta);
00395 coeff.real = (n+1) * B[n].real;
00396 coeff.imag = (n+1) * B[n].imag;
00397 Zeta_Denom = add(coeff, Zeta_Denom);
00398 }
00399 Zeta_sqr = multiply(Zeta, Zeta);
00400
00401 Zeta_Numer = multiply(Zeta_Numer, Zeta_sqr);
00402 Zeta_Numer = add(z, Zeta_Numer);
00403
00404 Zeta_Denom = multiply(Zeta_Denom, Zeta);
00405 Zeta_Denom = add(B[0], Zeta_Denom);
00406
00407 Zeta = divide(Zeta_Numer, Zeta_Denom);
00408 }
00409 dphi = D[8];
00410 for (n = 7; n >= 0; n--)
00411 dphi = dphi * Zeta.real + D[n];
00412 dphi *= Zeta.real;
00413
00414 double latitude = NZMG_Origin_Lat + (dphi * 1.0e5 / 3600.0);
00415 latitude *= PI / 180.0;
00416 double longitude = NZMG_Origin_Long + Zeta.imag;
00417
00418 if ((longitude > PI) && (longitude - PI < 1.0e-6))
00419 longitude = PI;
00420
00421 if ((latitude < MIN_LAT) || (latitude > MAX_LAT))
00422 throw CoordinateConversionException( ErrorMessages::latitude );
00423 if ((longitude < MIN_LON) || (longitude > MAX_LON))
00424 throw CoordinateConversionException( ErrorMessages::longitude );
00425
00426 return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
00427 }
00428
00429
00430
00431