[C++] 좌표변환 유틸 (WGS to KATEC , KATEC to WGS)

갑자기 ASP에서 좌표변환 유틸이 필요하다는 요청에 dll 로 적업하게 됨. 4년만에  Visual studio 를 사용하는 느낌은 너무 생소했다 ㅋ

뭐 이런저런 알고리즘이 있는거 같기는 하지만.. 솔직히 이해는 못하고. 기존에 작업된 자바코드를 토대로 그냥 C++로 변환(??) 했음. 알고리즘을 물어본다면.. 죄송하다는 말밖에는 ㅡㅡ;;

Coordinate.h

#include

typedef struct{
    double x;   //x좌표
    double y;   //y좌표
} DPoint;

typedef struct{
    double lon;                  
    double lat;                  
    double h;
} LonAndLat;

typedef struct{
    double x;                  
    double y;                  
    double h;
} TM;

typedef struct{
    int ind;
    double r_major, r_minor;
    double scale_factor;
    double lon_center, lat_origin;
    double false_northing, false_easting;
    double e0, e1, e2, e3;
    double e, es, esp, ml0;
} Adjustment;

class CCoordinate {
public:
    DPoint convertToKATEC(double dx, double dy);
    TM convLLToTM(LonAndLat inputL, int coordinates, int ellipse);
    LonAndLat convTMToLL(TM inputT, int inputCoordinates, int targetEllipse);
    LonAndLat datumTrans(LonAndLat inputL, int inputEllipse, int outputEllipse);
    Adjustment setAdjustment(int ellipse, int coordinates);
    TM LLtoTM(LonAndLat L, Adjustment A);
    LonAndLat TMToLL(TM T, Adjustment A);
    LonAndLat convertToWGS(int x, int y);

    double mlfn(double e0, double e1, double e2, double e3, double phi);
    double e0fn(double x);
    double e1fn(double x);
    double e2fn(double x);
    double e3fn(double x);
    double asinz(double con);
    double dabs(double x);
  
    BSTR m_x, m_y;
}

Coordinate.cpp

double DFACT = 36000.000000;
int KATEC = 3;
int WGS84 = 1;
double PI = 3.14159265358979323846;
double minor[]={6356078.96325f, 6356752.3142f};
double major[]={6377397.155f, 6378137.0f};

double dX_W2B = 128;
double dY_W2B = -481;
double dZ_W2B = -664;
int Bessel = 0; //TM은 Bessel 만 사용..
double SF[] = {1, 1, 1, 0.9999, 0.9996, 0.9996};
double FE[] = {200000, 200000, 200000, 400000, 500000, 500000};
double FN[] = {500000, 500000,  500000, 600000, 0.0, 0.0};
double LatCen[] = {0.663225115757845, 0.663225115757845, 0.663225115757845,  0.663225115757845,0, 0};
double LonCen[] = {2.18171200985643, 2.21661859489632, 2.2515251799362,  2.23402144255274, 2.25147473507269,  2.14675497995303};
double EPSLN = 0.0000000001;

DPoint CCoordinate::convertToKATEC(double dx, double dy)
{
    DPoint point;
    LonAndLat ll;

    memset(&point, 0x00, sizeof(DPoint));  
    memset(&ll, 0x00, sizeof(LonAndLat));  

    point.x = (double)dx*DFACT;
    point.y = (double)dy*DFACT;

    double dX = point.x/DFACT;
    double dY = point.y/DFACT;

    ll.lon = dX;
    ll.lat = dY;
    ll.h   = 0;
  
    TM tm = convLLToTM(ll, KATEC, WGS84);
  
    point.x = (int)tm.x;
    point.y = (int)tm.y;

    return point;
}

LonAndLat CCoordinate::convertToWGS(int x, int y)
{
    // 좌표계 초기화
    LonAndLat ll;
    TM tm;

    memset(&tm, 0x00, sizeof(TM));  
    memset(&ll, 0x00, sizeof(LonAndLat));  

    tm.x = (double)x;
    tm.y = (double)y;

    // 좌표 변환
    ll = convTMToLL(tm, KATEC, WGS84);
  
    return ll;

}

LonAndLat CCoordinate::convTMToLL(TM inputT, int inputCoordinates, int targetEllipse) {
    TM tm;
    memset(&tm, 0x00, sizeof(TM));  

    tm.x = inputT.x;
    tm.y = inputT.y;

    // TM좌표 -> 경위도로 변환
    //Adjustment A = new Adjustment();
    Adjustment A = setAdjustment(Bessel, inputCoordinates);

    LonAndLat inputL = TMToLL(tm, A);
    //System.out.println("좌표변환후 :"+inputL.lon +"  "+inputL.lat);
    //System.out.println(inputL.lon*180.0f/Math.PI+"   "+inputL.lat*180.0f/Math.PI);
    // 타원체를 변경한다.
    LonAndLat outputL = datumTrans(inputL, Bessel, targetEllipse);
    outputL.lon = outputL.lon*180.0f/PI;
    outputL.lat = outputL.lat*180.0f/PI;
    return outputL;
}

TM CCoordinate::convLLToTM(LonAndLat inputL, int coordinates, int ellipse) {

    LonAndLat LL;
    memset(&LL, 0x00, sizeof(LonAndLat));  

    LL.lon = inputL.lon;
    LL.lat = inputL.lat;
    LL.h   = 0;

    //convert degree to radian
    LL.lon = inputL.lon * PI / 180.0f;
    LL.lat = inputL.lat * PI / 180.0f;

    // 타원체를 먼저 변경한다.
    LonAndLat outputL  =  datumTrans(LL, ellipse,  Bessel);
    //System.out.println("타원체 변경후:"+outputL.lon*180/Math.PI+"  "+outputL.lat*180/Math.PI);
    //Adjustment A = new Adjustment();
    // 경위도 좌표를 TM좌표로 변경한다.  
    Adjustment A = setAdjustment(Bessel, coordinates);
    TM T = LLtoTM(outputL, A);
    T.h = outputL.h;
    return T;

}

LonAndLat CCoordinate::TMToLL(TM T, Adjustment A) {
    //System.out.println("A.e0"+A.e0);
    LonAndLat L;
    memset(&L, 0x00, sizeof(LonAndLat));  

    double con, phi, delta_phi, sin_phi, cos_phi, tan_phi;
    double c, cs, t, ts, n, r, d, ds, f, h, g, temp;
    long max_iter = 6;
    long i;
    if (A.ind != 0) {
        f = exp(T.x / (A.r_major * A.scale_factor));
        g = 0.5 * (f - 1.0f / f);
        temp = A.lat_origin + T.y / (A.r_major * A.scale_factor);
        h = cos(temp);
        con = sqrt((1.0f - h * h) / (1.0f + g * g));
        L.lat= asinz(con);
        if (temp < 0)
            L.lat = -L.lat;
        if ((g == 0) && (h == 0))
            L.lon = A.lon_center;
        else
            L.lon = atan(g / h) + A.lon_center;
    }
    // TM to LL inverse equations from here
    T.x = T.x - A.false_easting;
    T.y = T.y - A.false_northing;
    con = (A.ml0 + T.y / A.scale_factor) / A.r_major;
    phi = con;
    i = 0;
    while (true) {
        delta_phi = ((con + A.e1 * sin(2.0f * phi) - A.e2 * sin(4.0f * phi)
                + A.e3 * sin(6.0f * phi)) / A.e0) - phi;
        phi = phi + delta_phi;
        if (dabs(delta_phi) <= EPSLN)
            break;
        if (i >= max_iter) {
            //ShowMessage("Latitude failed to convert...");
        //    System.out.println("Latitude failed to convert...");              
        }
        i++;
    }
    if (fabs(phi) < (PI / 2.0f)) {
        sin_phi = sin(phi);
        cos_phi = cos(phi);
        tan_phi = tan(phi);
        c = A.esp * (cos_phi * cos_phi);
        cs = c * c;
        t = tan_phi * tan_phi;
        ts = t * t;
        con = 1.0f - A.es * sin_phi * sin_phi;
        n = A.r_major / sqrt(con);
        r = n * (1.0f - A.es) / con;
        d = T.x / (n * A.scale_factor);
        ds = d * d;
        L.lat = phi - (n * tan_phi * ds / r) * (0.5f - ds / 24.0f *
            (5.0f + 3.0f * t + 10.0f * c - 4.0f * cs - 9.0f * A.esp - ds / 30.0f *
            (61.0f + 90.0f * t + 298.0f * c + 45.0f * ts - 252.0f * A.esp - 3.0f * cs)));
        L.lon = A.lon_center + (d *
            (1.0f - ds / 6.0f *
            (1.0f + 2.0f * t + c - ds / 20.0f *
            (5.0f - 2.0f * c + 28.0f * t - 3.0f * cs + 8.0f * A.esp + 24.0f * ts))) / cos_phi);
//        *lat = 127.2866484932;
//        *lon = 37.4402108468;
    } else {
        L.lat = PI/2.0f * sin(T.y);
        L.lon = A.lon_center;
    }
    return L;
}

LonAndLat CCoordinate::datumTrans(LonAndLat inputL, int inputEllipse, int outputEllipse) {
    LonAndLat outputL;
    memset(&outputL, 0x00, sizeof(LonAndLat));  

    double input_a        = major[inputEllipse];
    double input_b        = minor[inputEllipse];
    double output_a        = major[outputEllipse];
    double output_b        = minor[outputEllipse];
    double delta_a, delta_f, delta_Phi, delta_Lamda, delta_H;
    double Rm, Rn;
    double temp, es_temp;
    int gap = inputEllipse - outputEllipse;
    double delta_X = gap*dX_W2B;
    double delta_Y = gap*dY_W2B;
    double delta_Z = gap*dZ_W2B;

    temp = input_b / input_a;
    es_temp = 1.0f - temp*temp;
    delta_a = output_a - input_a;
    delta_f = input_b / input_a - output_b / output_a;

    Rm =(input_a * (1.0f - es_temp)) / pow((1.0f - es_temp * sin(inputL.lat) * sin(inputL.lat)), (3.0f/2.0f));
    Rn = input_a / pow((1.0f - es_temp * sin(inputL.lat) * sin(inputL.lat)), (1.0f/ 2.0f));
    delta_Phi = ((((-delta_X * sin(inputL.lat) * cos(inputL.lon) - delta_Y * sin(inputL.lat) * sin(inputL.lon))
                            + delta_Z * cos(inputL.lat)) + delta_a * Rn * es_temp * sin(inputL.lat) * cos(inputL.lat) / input_a)
                            + delta_f * (Rm / temp + Rn * temp) * sin(inputL.lat) * cos(inputL.lat)) / (Rm + inputL.h);
    //System.out.println(delta_Phi);
    delta_Lamda = ((-delta_X) * sin(inputL.lon) + delta_Y * cos(inputL.lon)) / ((Rn + inputL.h) * cos(inputL.lat));
    delta_H = delta_X * cos(inputL.lat) * cos(inputL.lon) + delta_Y * cos(inputL.lat) * sin(inputL.lon)
            + delta_Z * sin(inputL.lat) - delta_a * input_a / Rn + delta_f * temp * Rn * sin(inputL.lat) * sin(inputL.lat);
    //System.out.println(delta_Lamda);

    outputL.lat = inputL.lat + delta_Phi;
    outputL.lon = inputL.lon + delta_Lamda;
    outputL.h = inputL.h + delta_H;

    return outputL;
}

Adjustment CCoordinate::setAdjustment(int ellipse, int coordinates) {

    Adjustment A;
    memset(&A, 0x00, sizeof(Adjustment));  

    A.r_major = major[ellipse];
    A.r_minor = minor[ellipse];
    A.scale_factor = SF[coordinates];
    A.lon_center = LonCen[coordinates];
    A.lat_origin = LatCen[coordinates];
    A.false_northing = FN[coordinates];
    A.false_easting = FE[coordinates];
    double temp = A.r_minor / A.r_major;
    A.es = 1.0f - temp*temp;
    A.e = sqrt(A.es);
    A.e0 = e0fn(A.es);
    A.e1 = e1fn(A.es);
    A.e2 = e2fn(A.es);
    A.e3 = e3fn(A.es);
    A.ml0 = A.r_major * mlfn(A.e0, A.e1, A.e2, A.e3, A.lat_origin);
    A.esp = A.es / (1.0f - A.es);
    if(A.es < 0.00001f)
            A.ind = 1;
    else
            A.ind = 0;

    return A;
}

TM CCoordinate::LLtoTM(LonAndLat L, Adjustment A)
{
    TM T;
    memset(&T, 0x00, sizeof(TM));  

    double delta_lon;
    double sin_phi, cos_phi;
    double al, als, b=0, c, t, tq;
    double con, n, ml;

    // LL to TM forward equations from here;
    delta_lon = L.lon - A.lon_center;
    sin_phi = sin(L.lat);
    cos_phi = cos(L.lat);

    if (A.ind != 0) {
            b =  cos_phi * sin(delta_lon);
            if ((fabs(fabs(b) - 1.0f)) < 0.0000000001f) {
                    //ShowMessage("The point is going to the unlimited value...\n");
                    //System.out.println("The point is going to the unlimited value...");
            }
    } else {
            T.x = 0.5f * A.r_major * A.scale_factor * log((1.0f + b) / 1.0f- b);
            con = acos(cos_phi * cos(delta_lon) / sqrt(1.0f - b * b));
            if (L.lat < 0) {
                    con = -con;
                    T.y = A.r_major * A.scale_factor * (con - A.lat_origin);
            }
    }

    al = cos_phi * delta_lon;
    als = al * al;
    c = A.esp * cos_phi * cos_phi;
    tq = tan(L.lat);
    t = tq * tq;
    con = 1.0f - A.es * sin_phi * sin_phi;
    n = A.r_major / sqrt(con);
    ml = A.r_major * mlfn(A.e0, A.e1, A.e2, A.e3, L.lat);

    T.x = A.scale_factor * n * al * (1.0f + als / 6.0f *
            (1.0f - t + c + als / 20.0f *
            (5.0f - 18.0f * t + t * t + 72.0f * c - 58.0f * A.esp)
            )
            ) + A.false_easting;
    T.y = A.scale_factor * (ml - A.ml0 + n * tq *
            (als *
            (0.5f + als / 24.0f *
            (5.0f - t + 9.0f * c + 4.0f * c * c + als / 30.0f *
            (61.0f - 58.0f * t + t * t + 600.0f * c - 330.0f * A.esp)
            )))) + A.false_northing;
    return T;
}

double CCoordinate::mlfn(double e0, double e1, double e2, double e3, double phi)
{
     return (e0 * phi - e1 * sin(2.0f * phi) + e2 * sin(4.0f * phi) - e3 * sin(6.0f * phi));
}

double CCoordinate::e0fn(double x)
{
        return (1.0f - 0.25f*x*(1.0f + x / 16.0f * (3+1.25f*x)));
}
double CCoordinate::e1fn(double x)
{
        return(0.375f * x * (1. + 0.25f*x*(1.0f + 0.46875f*x)));
}
double CCoordinate::e2fn(double x)
{
        return(0.05859375f*x*x*(1.0f+0.75f*x));
}
double CCoordinate::e3fn(double x)
{
        return (x * x * x * (35.0f / 3072.0f));
}

double CCoordinate::asinz(double con)
{
    if (fabs(con) > 1.0f)
        if (con > 1.0f)
            con = 1.0f;
        else
            con = -1.0f;
    return asin(con);
}
double CCoordinate::dabs(double x)
{
    if(x>=0) return x;
    else return -x;      
}

Dll 사용방법

<%
     Set  trans = Server.CreateObject("Interactive.Coordinate")
    
     // WGS 타입을 KATEC 타입으로 변환.
   
     trans.WGS2KATEC 126.83530, 37.62508
     Response.write trans.X
     Response.write " : "
     Response.write trans.Y
   
     // KATEC 타입을 WGS 타입으로 변화
   
     trans.KATEC2WGS 297386, 558722
     Response.write "             "
     Response.write trans.X
     Response.write " : "
     Response.write trans.Y
%>