[GIS] IDW(Inverse Distance Weighting)

IDW는 이미 알고 있는 값으로부터 알고자 하는 값을 보간하는 방법입니다. IDW를 사용하여 주어진 점 x에 대한 보간된 값 u를 결정하는 일반화된 형태의 보간 함수는 다음과 같습니다.

사용자 삽입 이미지
N은 이미 알고 있는 값의 개수, w는 가중치의 값, u는 앞서 말한 계산되어 나온 보간된 값입니다. IDW에서 중요한것은 가중치의 값에 해당되는 w에 대한 함수가 여러개 존재하며, Spepard 방식과 Liszka 방식 그리고 이들의 변종이 존재합니다. 여기서는 Spepard 방식에 대한 w값에 대해서만 살펴 보겠습니다.

사용자 삽입 이미지
위의 식이 Spepard님이 정의한 가중치 w의 값입니다. d는 보간하고자 하는 지점(x)와 이미 알고 있는 지점(xk) 사이의 공간적 거리입니다. 바로 거리(Distance)의 역(Inverse)에 대한 가중치(Weighting)라는 의미로 IDW가 된 것입니다.

p는 0보다 큰 실수값입니다. 이 p값의 범위에 따라 전체적인 보간된 양상이 다양하게 결정됩니다. p의 범위가 0~1이면 전체적인 양상이 좁고 날카로우며 1보다 크면 넓고 부드럽게 퍼져서 보간이 됩니다. 눈에 보이는 보간된 양상을 글로써 표현하려니 한계가 있는데… 이 부분에 대해서 실제 구현을 통해 살펴보도록 하겠습니다.

끝으로 IDW처럼, 이미 알고 있는 값을 통해 다른 값을 추정(보간)하는 방법 중 Kriging 기법이 있습니다. 기회가 닿는 다면 이 Kriging 기법에 대해서도 논의해 보고 싶습니다.

[GIS] 개발중인 맵 엔진으로 표현한 통계 지도

간단하게 서울시 행정구역도와 가상의 데이터를 이용해서 구성해본 통계 지도입니다. 먼저 첫번째와 두번째는 남자, 여자, 노인이라는 항목에 대한 인구수를 가상의 데이터로 하여 Bar Chart 형식으로 표현한 내용이며 두번째는 동일한 주제로 하여 Pie Chart 형태로 표현한 화면입니다.

사용자 삽입 이미지
사용자 삽입 이미지
그리고 아래, 세번째는 실제 아파트의 밀집 정도를 색상 밀도도로 표현한 것입니다. 그리고 네번째는 표현한 밀집도를 등가선을 추출해 표현한 것입니다. 밀도도와 등가선을 추출하는 기능은 개발중인 맵 엔진에 Plug-In 할 수 있는 별도의 확장 기능 개념으로 개발하였습니다.

사용자 삽입 이미지
사용자 삽입 이미지
아래는 좀더 좁은 간격으로 등가선을 딴 결과 화면입니다.

사용자 삽입 이미지
이외에도 개발중인 맵 엔진를 이용하여 지도 위에 다양한 형태의 통계치를 표현할 수 있습니다.

보간 – Catmull-Rom Spline

많은 스플라인의 종류 중에 하나인 큐빅 스플라인을 1차원의 보간에 적용하는 것에 대해 살펴보겠습니다. Catmull-Rom 스플라인을 구성하는 구분된 부드러운 곡선들을 나타내는 키 프레임 집합을 가지며 모든 키는 곡선 상에 위치합니다. 이 루틴을 사용하기 위해서 4개의 키 프레임 값이 필요합니다. 이 4개의 키 값을 v0, v1, v2, v3라고 하고, 여기에 보간을 위하여 v1에서 v2 사이의 지정된 0~1까지의 범위를 가지는 실수값 x가 존재합니다. 아래의 f(x)의 반환값은 x값에 의해 결정이 됩니다.

사용자 삽입 이미지
여기서 M은 다음처럼 정의됩니다.

사용자 삽입 이미지
아래의 이미지는 v1에서 v2 사이의 곡선의 한 예를 나타난 것입니다. 이 곡선은 위의 수식에서 x 값을 0에서 1.0 사이의 값을 이용해 얻을 수 있습니다.

사용자 삽입 이미지

아래의 코드는 위에서 설명한 내용을 C언어로 구현한 것입니다.

/* Coefficients for Matrix M */
#define M11  0.0 
#define M12  1.0
#define M13  0.0
#define M14  0.0
#define M21 -0.5
#define M22  0.0
#define M23  0.5
#define M24  0.0
#define M31  1.0
#define M32 -2.5
#define M33  2.0
#define M34 -0.5
#define M41 -0.5
#define M42  1.5
#define M43 -1.5
#define M44  0.5

double catmullRomSpline(float x, float v0,float v1, float v2,float v3) {
    double c1,c2,c3,c4;

    c1 = M12*v1;
    c2 = M21*v0 + M23*v2;
    c3 = M31*v0 + M32*v1 + M33*v2 + M34*v3;
    c4 = M41*v0 + M42*v1 + M43*v2 + M44*v3;

    return(((c4*x + c3)*x +c2)*x + c1);
}

이 글의 원문은 http://www.lighthouse3d.com/opengl/maths/index.php?catmullrom 입니다.

삼각형 폴리곤과 편선(Ray)의 교차점

이 글은 http://www.lighthouse3d.com/opengl/maths/index.php?raytriint의 글을 통해 재구성한 글로, 변역하면서 쉽게 이해할 수 있도록 내용을 약간 수정했습니다.

사용자 삽입 이미지
매개변수방정식의 형태로 편선(시작점에서 어떤 방향으로 무한하게 진행하는 선, Ray)과 삼각형의 폴리곤이 주어질때.. 이 둘이 교차하는가를 판단하는 방법에 대해 증명은 생략하고 그 방법에 대해서 알아 보겠습니다.

삼각형을 구성하는 점은 다음처럼 정의할 수 있습니다.

삼각형의  구성 점(u,v) = (1-u-v)*p0 + u*p1 + v*p2
여기서, p0, p1, p2는 삼각형 위의 이미 알고 있는 점, u >= 0, v >= 0, u + v <= 1.0

그리고 편선에 대한 매개변수방정식을 다음처럼 정의할 수 있습니다.

편선의 구성 점(t) = p + t * d
여기서, p는 편선의 시작점으로 이미 알고 있는 점, d는 편선의 방향 벡터

이제 편선과 삼각형의 교점은 삼각형의 구성 점(u,v) = 편선의 구성 점(t) 로부터 다음과 같습니다.

p + t * d = (1-u-v) * p0 + u*p1 + v*p2

결국… 교차 여부는 세 값(t, u, v)가 위의 공식을 만족하고 u와 v가 앞서 정의한 조건을 만족하면 교차하는 경우이다. 이런 수학적인 논의는 다음으로 미루기로 하고…. 여러분이 원하는 C 코드로 대신합니다.

/* a = b - c */
#define vector(a,b,c) \
        (a)[0] = (b)[0] - (c)[0]; \
        (a)[1] = (b)[1] - (c)[1]; \
        (a)[2] = (b)[2] - (c)[2];

int rayIntersectsTriangle(float *p, float *d, float *v0, float *v1, float *v2) {
    float e1[3],e2[3],h[3],s[3],q[3];
    float a,f,u,v;
 
    vector(e1,v1,v0);
    vector(e2,v2,v0);
    crossProduct(h,d,e2);
    a = innerProduct(e1,h);
 
    if (a > -0.00001 && a < 0.00001) return(false);
 
    f = 1/a;
    vector(s,p,v0);
    u = f * (innerProduct(s,h));
 
    if (u < 0.0 || u > 1.0) return (false);
 
    crossProduct(q,s,e1);
    v = f * innerProduct(d,q);
    if (v < 0.0 || u + v > 1.0) return (false);

    // at this stage we can compute t to find out where 
    // the intersection point is on the line
    t = f * innerProduct(e2,q);
    if (t > 0.00001) // ray intersection
        return (true);
    else // this means that there is line intersection but not ray intersection
        return (false);
}

vector 함수는 두개의 좌표로부터 백터를 만드는 함수이고 innerProduct와 crossProduct는 내적과 외적을 구하는 함수입니다.

[GIS] 카텍좌표로부터 GoogleMap 이미지 다운로드

2009년 2월 11일 변경된 구글지도서버의 내용을 반영해, 정상적으로 작동하도록 하였습니다. 소중한 정보를 제공해준 이사돌이님께 감사를 드립니다. E-Mail이나 블로그 링크를 공개하지 않으셔서 이 자리를 통해 감사의 말씀을 드립니다. 또한 소스까지 공개합니다. 아무쪼록 필요하신 분에게 도움이 되시길 바랍니다.

구글에서 제공하는 구글맵 OpenAPI에 대한 관심이 저를 제외하고 매우 뜨거웠던 모양입니다. 관심을 갖게된 배경은 회사에서 개발하고 있는 엔진에 다른 회사에서 제공하는 OpenAPI를 통해 제공하고 있는 지도를 적용해보면 어떨까 하는 생각에서 접근을 하게 되었습니다. 구글맵을 통해 OpenAPI를 접해보는 순간… 머리속에 하얗게 되는것을 느꼈습니다. 웹이라면 단순이 html 정도와 자바스크립트 그리고 근사한 이미지나 플래시로 치장하는게 전부라고 생각했기 때문입니다. OpenAPI는 많은 사람들에게 공개된 API인지라 사용하기가 매우 쉬웠습니다만, 개발자로써 내부적으로 어떻게 구현했느냐를 모르니 참… 우울해지더군요. 바로 자바스크립트 관련 책을 주문하고 OpenAPI와 AJAX와 같은 흐름을 대표하는 Web2.0에 찬찬히… 개발자로써 접근을 시도해 보고자 다짐하게 되었습니다.

뭐 여튼… 처음으로 돌아가, 원래의 목적을 완수하고 위해.. 구글맵을 엔진단에 통합하는 방법을 연구하던 중에, 그 중간 결과로 만든 것을 소개해 드릴겸해서 글을 작성합니다. 구글맵은 경도와 위도값을 주게 되면 그 지점에 대한 이미지를 던져주는데, 이 이미지의 크기는 256 x 256 입니다. 이 이미지를 타일이라고 합니다. 구글맵에는 해상도 레벨에 대한 개념이 있습니다. 1~18 단계로 되어 있고 1단계는 하나의 타일에 전세계가 보입니다. 18단계는 최대로 확대한 타일이미지로 18단계에 해당되는 타일의 개수는 엄청나겠지요. 이런 구글맵에서 제공하는 지도를 활용하기 위해서는 원하는 지점에 대한 타일을 얻는 것이 시발점이 됩니다. 그런데 문제는 이 원하는 지점을 어떤식으로 접근할 것이냐인데… 우리 나라 GIS계에서 사용하는 좌표계는 카텍입니다. 원래 카텍은 CNS 분야에서 사용하는 좌표계인데 원점이 하나인지라… 직각좌표계에서 경위도로 변경할때 원점을 고민할 필요가 없습니다. 바로 이 카텍 좌표를 이용해서 구글이 사용하는 경위도 좌표로 변환하고, 이 변환된 경위도 좌표를 구글맵 서버에게 넘겨줘서 해당 지점의 타일을 받아오면 됩니다.

먼저 구글맵은 WGS84 경위도 좌표계이고 우리나라에서 많이 사용하는 좌표계는 베셀 TM 직각좌표계의 원점을 하나로 한 카텍좌표계입니다. 먼저 경위도 좌표를 카텍과 같은 직각 좌표로 변환하는 방법과 그 반대로 변환하는 방법이 필요한데, 이 부분에 대한 방법은 gTrans라는 Visual Basic으로된 소스를 참고로 하여 작성을 하였습니다. 하지만 이 코드는 카텍좌표에 대한 언급이 전혀 없습니만… 타원체와 투영 그리고 변환공식 등에 대한 기본적인 내용을 숙지 하고 있다면 카텍 좌표계로의 변환은 그리 어렵지 않습니다. 비주얼 베이직으로 된것을 C++으 변환하는게 시간적으로 어려울 뿐이지요. 대학원때 GPS를 분석해보며 살펴봤던 좌표변환이 이제서야 큰 도움이 되더군요. 참고로 카텍 좌표계는 TM128이라고 하는데, 원래 우리나라에서 사용하는 TM 좌표계와 다른 점은 아래와 같고 다른 점은 모두 동일합니다.

central scale : 0.9998 or 0.9999
central meridian : E 128 00 00
Origin Latitude : N 38 00 00
False Easting (meters) : 400000
Falsle Northing (meters) : 600000

위의 내용은 http://blog.naver.com/azzuo/100037410143에서 참조했습니다. 위의 차이점을 활용해서 구글의 WGS84 경위도 좌표를 카텍좌표계로 변환하고 또 이 반대로 변환하는 코드를 작성하여 구글서버에 요청을 하여 타일을 받아오는 코드를 작성하였는데… 또 여기서 문제는 원하는 지점(카텍좌표)를 구글이 사용하는 경위도 좌표계로 변환을 하여 구글 서버에게 그 지점의 타일을 요청하는 url을 작성하는 것이 문제입니다. 이에 대한 것은 아래의 GoogleMap, How to Work라는 아티클을 참고하였습니다.

역시 인터넷은 바다입니다. 찾으면 있습니다… 이제 원하는 지점에 대한 타일요청 url을 알았으니 구글맵 서버로부터 쉽게 타일 이미지를 받아올 수 있지요. 이제 실제로 제가 만든 프로그램이 실행되는 모습을 말씀드리겠습니다.

먼저 요청할 지점입니다. 카텍 좌표로 만들어진 지도인데… 제주도와 서울의 여의도 그리고 꼬리 부분에 대한 이미지를 요청해 보겠습니다.

이 세 지점에 대한 좌표가 나타나 있습니다. 이 세 지점의 좌표는 ArcMap에 지도를 올려 놓고 얻었습니다. 이 좌표를 구글맵이 사용하는 WGS84 경위도 좌표계로 변환하고 타일을 요청하는 url을 만들어 구글맵 서버에 요청하면 되는데… 먼저 서울의 여의도 지점을 요청해 보겠습니다.

tstGoogle.exe라는 프로그램을 실행할때 인자로 카텍좌표계로써의 X, Y와 줌레벨(여기서는 13)을 주고 구글서버에게 받아 저장할 타일에 대한 파일명을 주고 실행을 하면 받은 타일에 대한 정보을 출력해주고.. “똥 싸는 중…”라는 메세지와 함께 타일이미지를 저장하게 됩니다. ^^; (최신버전에서는 얌전하게 Downloading으로 변경되었습니다) 아래는 요청해서 받은 타일이미지입니다.


여의도가 담긴 타일을 받았다는 것을 알 수가 있습니다. 이제 나머지 다른 지점에 대해서도 요청을 해서 이미지를 얻어올 수가 있습니다. 아래는 나머지 두 타일에 대한 요청 명령이고 받은 타일 이미지입니다.

tstGoogle.exe 265630 88249 11 d:/제주도.jpg
tstGoogle.exe 541018 385821 5 d:/꼬리.jpg

아래는 tstGoogle 프로그램입니다. 필요하신 분은 사용하시길 바랍니다. 아마도 배치파일을 만들어서 해당 지점의 타일이미지를 긁어 오는 곳에 응용하면 어떨지 싶습니다.

2009년 2월 11일 변경된 구글지도서버의 내용을 반영해, 정상적으로 작동하도록 하였습니다. 소중한 정보를 제공해준 이사돌이님께 다시 한번 더 감사를 드립니다.

아~ 끝으로.. 타일을 보면 요청한 좌표의 지점이 타일의 중앙에 있지 않은데, 구글맵의 타일은 사용자가 요청한 지점이 포함된 타일을 보내준다는 점에 유의하셔야 합니다.

아래의 그림은 제가 만든 프로그램을 이용해서 GoogleMap 서버로부터 타일이미지를 얻어온후 ArcGIS에서 수치지도(Shape)파일과 함께 중첩한 모습입니다. Georeferencing을 위한 2개의 컨트롤포인트(빨강색 십자 심벌)는 제가 만든 프로그램이 계산해준 좌표값입니다.