Develope
GPS로부터 받은 WGS84 좌표계를 Bessel Datum을 이용한 TM 128 (KATEC) 평면좌표계로 변환하는 소스 코드
친절한 웬디양~ㅎㅎ
2009. 9. 22. 21:34
GPS로부터 받은 WGS84 좌표계를 사용하지말고 평면 좌표계로 변환하여 활용할 것!!!
============================================================
<2007년 강의자료에 있음>
// GPS로부터 받은 좌표를 wx, wy에 넣어서 호출한 후에 TM 평면 좌표계인 xx, yy 사용할 것
void ConvertCoordinates(double wx, double wy, double *xx, double *yy)
{
wx -= 12800.;
wy -= 3600.;
wx = (128. + wx / 60.);
wy = (36. + wy / 60.);
double tx, ty;
double bx, by;
wgs2bessel(wx, wy, &bx, &by);
bessel2tm(bx, by, xx, yy, 129.0, 38.0);
}
=============================================================================
#include "coord.h"
void detorgbl(double *pbx, double *orglongitude)
{
if (*pbx >= 124 && *pbx < 126) {
*orglongitude = 125.0;
} else if (*pbx >= 126 && *pbx < 128) {
*orglongitude = 127.0;
} else if (*pbx >= 128 && *pbx < 130) {
*orglongitude = 129.0;
} else {
*orglongitude = 0;
}
*orglongitude += 10.405/3600.0;
*pbx -= *orglongitude;
}
void deg2rad(double *ptx, double *pty)
{
(*ptx) *= (M_PI/180.0);
(*pty) *= (M_PI/180.0);
}
void rad2deg(double *ptx, double *pty)
{
(*ptx) *= (180.0/M_PI);
(*pty) *= (180.0/M_PI);
}
double Bfnxco(double lat)
{
double fnxco;
fnxco = BESSEL_A * (1 - BESSEL_EE) * ((BA * lat) - 0.5 * BB * sin(2 * lat) +
0.25 * BC * sin(4 * lat) - 1. / 6. * BD * sin(6 * lat));
return fnxco;
}
/*
* 124 <= bx < 126: orglon = 125.0 + 10.405 / 3600.0; (서부원점)
* 126 <= bx < 128: orglon = 127.0 + 10.405 / 3600.0; (중부원점)
* 128 <= bx < 130: orglon = 129.0 + 10.405 / 3600.0; (동부원점)
*
* orglat = 38.0;
*
*/
void bessel2tm(double bx, double by, double *ptx, double *pty,
double orglongitude, double orglatitude)
{
double ptblx, ptbly;
double orgscalefactor;
double b, N, e_prime, nn, t, e, n;
ptblx = bx;
ptbly = by;
orgscalefactor = 1.0;
orglongitude += 10.405 / 3600.0;
ptblx -= orglongitude;
deg2rad(&ptblx, &ptbly);
orglatitude *= (M_PI / 180.0);
b = Bfnxco(ptbly) - Bfnxco(orglatitude);
N = BESSEL_A / sqrt(1 - BESSEL_EE*pow(sin(ptbly), 2.));
e_prime = BESSEL_EE / (1 - BESSEL_EE);
nn = e_prime * pow(cos(ptbly), 2.);
t = tan(ptbly);
e = orgscalefactor * (ptblx * N * cos(ptbly) + pow(ptblx, 3.) /
6. * N * pow(cos(ptbly), 3.) * (1 - t * t + nn) +
pow(ptblx, 5.) / 120. * N * pow(cos(ptbly), 5.) *
(5 - 18 * t * t + pow(t, 4.)));
n = orgscalefactor * (b + pow(ptblx, 2.) /
2. * N * sin(ptbly) * cos(ptbly) + pow(ptblx, 4.) /
24. * N * sin(ptbly) * pow(cos(ptbly), 3.) * (5 - t * t + 9 *nn));
*ptx = e + FALSE_EAST;
*pty = n + FALSE_NORTH;
}
void wgs2bessel(double wx, double wy, double *pbx, double *pby)
{
double rn, rm, d_pi, d_lamda, d_h;
double h;
*pbx = wx;
*pby = wy;
h = 0.0;
deg2rad(&wx, &wy);
rn = WGS84_A / sqrt(1 - WGS84_EE*pow(sin(wy), 2.));
rm = WGS84_A*(1-WGS84_EE) / pow(sqrt(1 - WGS84_EE*pow(sin(wy), 2.)), 3.);
d_pi = (-W2B_DELTAX * sin(wy) * cos(wx) -
W2B_DELTAY * sin(wy) * sin(wx) +
W2B_DELTAZ * cos(wy) +
W2B_DELTAA * (rn * WGS84_EE * sin(wy) * cos(wy)) / WGS84_A +
W2B_DELTAF * (rm * WGS84_A / WGS84_B + rn * WGS84_B / WGS84_A) *
sin(wy) * cos(wy)) /
((rm + h) * sin(M_PI/180. * 1/3600.));
d_lamda = (-W2B_DELTAX * sin(wx) + W2B_DELTAY * cos(wx))
/ ((rn + h) * cos(wy) * sin(M_PI/180. * 1/ 3600.));
d_h = W2B_DELTAX * cos(wy) * cos(wx) +
W2B_DELTAY * cos(wy) * sin(wx) +
W2B_DELTAZ * sin(wy) - W2B_DELTAA * WGS84_A / rn +
W2B_DELTAF * WGS84_B / WGS84_A * rn * pow(sin(wy), 2.);
*pbx += d_lamda / 3600.0;
*pby += d_pi / 3600.0;
}
|
|
==============================================================================
coord.h (2.4 KB), Download : 31
#ifndef _COORD_H_
#define _COORD_H_
#include "math.h"
#ifndef M_PI
#define M_PI ((double)3.14159265358979323846)
#endif
/*
* TM(제주도일 경우 FALSE_NORTH = 550000)
*/
#define FALSE_EAST (200000)
#define FALSE_NORTH (500000)
/*
* BESSEL
*/
#define BESSEL_A 6377397.155
#define BESSEL_RF 299.1528128
#define BESSEL_B (BESSEL_A - BESSEL_A / BESSEL_RF)
#define BESSEL_EE (0.006674372231315) // (a*a - b*b)/(a*a)
#define BA (1.005037306048555)
#define BB (0.005047849240300)
#define BC (0.000010563786831)
#define BD (0.000000020633322)
/*
* W2B
*/
#define W2B_DELTAX 128
#define W2B_DELTAY -481
#define W2B_DELTAZ -664
#define W2B_DELTAA -739.845
#define W2B_DELTAF -0.000010037483
/*
* WGS84
*/
#define WGS84_A 6378137.0
#define WGS84_F (1.0 / 298.257223563)
#define WGS84_RF 298.257223563
#define WGS84_B (WGS84_A - WGS84_A / WGS84_RF)
#define WGS84_EE (2.0 * WGS84_F - WGS84_F * WGS84_F)
#define WA (1.00503730604896035494)
#define WB (0.00504784924129776918)
#define WC (0.00001056378689647333)
#define WD (0.00000002063334976042)
/*
* B2W
*/
#define B2W_DELTAX -147
#define B2W_DELTAY 506
#define B2W_DELTAZ 687
#ifdef __cplusplus
extern "C" {
#endif
/*
* MISC
*/
void detorgbl(double *pbx, double *orglongitude);
void deg2rad(double *ptx, double *pty);
void rad2deg(double *ptx, double *pty);
double Bfnxco(double lat);
double Wfnxco(double lat);
/*
* -- UTM
* 51 Zone: 120 - 126 --> orglongitude = 123.0;
* 52 Zone: 126 - 132 --> orglongitude = 129.0; [default in COREA]
*
* -- TM
* 124 <= bx < 126: orglon = 125.0 + 10.405"(초) = 125.00289; (서부원점)
* 126 <= bx < 128: orglon = 127.0 + 10.405"(초) = 127.00289; (중부원점)
* 128 <= bx < 130: orglon = 129.0 + 10.405"(초) = 129.00289; (동부원점)
*
* orglat = 38.0;
*
*/
/*
* COORD CONVERSION
*/
void bessel2tm(double bx, double by, double *ptx, double *pty,
double orglongitude, double orglatitude);
void tm2bessel(double tx, double ty, double *pbx, double *pby,
double orglongitude);
void bessel2utm(double bx, double by, double *pux, double *puy,
double orglongitude);
void utm2bessel(double ux, double uy, double *pbx, double *pby,
double orglongitude);
void bessel2wgs(double bx, double by, double *pwx, double *pwy);
void wgs2bessel(double wx, double wy, double *pbx, double *pby);
void wgs2tmII(double lon, double lat, double Orlon, double Orlat,
double *ptx, double *pty);
#ifdef __cplusplus
}
#endif
#endif /* _COORD_H_ */
