Cosmic Background Noise

작년 10월부터 OSM 간보기를 하다가 올 3월초부터 OSM 작업을 시작했는데, 4개월만 더 있으면  OSM 작업하겠다고 생각한지 만 1년이 된다.

저번 주와 마찬가지로 이번 주에도 작업을 그리 많이 하지 못했다. 이번 주는 화요일, 수요일 합쳐 약 7시간 정도 작업. 이 기사 쓰는데 쓰는 시간이 더 걸리는 것 같다. 국가교통정보센터의 지도가 내년에 업데이트 되어도 그것을 손쉽게 업데이트할 방법을 찾다보니... 도로 작업에 너무 많은 시간을 까먹고 있다.

국가교통정보센터의 지도에서 추출한 도로로 다음 네 단계로 작업으로 나눴다. 라우팅은 미심쩍은 것들이 많아 이번 주 작업에서 뺐다.

1. 노이즈 제거

이런 노이즈의 '크기'는 대략 1-4m 정도. 2점/3점까지 수평/연직이 발견되면 한 점이나 두 점을 삭제하고 2점/3점의 평균 거리로 첫 번째 위치를 이동하여 보정했다.

 

2. 선로 합치기

 

국도 및 지방도는 고속도로처럼 뚜렷한 경계선에 의해 구분되지 않으므로 2개의 라인이 서로 다른 방향을 가진 벡터로 진행할 필요가 없다. 물론 국도나 지방도 중에서 4/8차선이 없는 것은 아니지만, 2개의 라인으로 표현되면 파일 크기가 커지고, 굳이 그렇게 도로 폭을 크게 묘사할 필요가 없다고 판단해 국도와 지방도의 선로를 합치기로 했다.

 

선로를 합치기 위한 조건은, 두 벡터의 시작점/끝점이 일정 거리 내에 인접해 있고 벡터의 방향만 반대일 때이다. 아울러 polyline의 이름과 타잎이 일치해야 한다.

 

시작점, 끝점은 경위도좌표로 나타나므로 두 좌표의 거리를 재어 파라미터로 정한 거리 안에 들어오면 인접했다고 판단할 수 있다. 평면 상의 두 점의 거리는 sqrt( (x2 - x1)^2 + (y2 - y1)^2 ). 물론 경위도로 표현되는 두 좌표는 구체인 지구 구면의 거리가 되어야 하므로 평면 거리는 쓸모없다. 구면 상의 두 점의 거리는 구면기하 삼각함수를 이용한 Haversine 공식을 사용해 계산:

 

d = acos(sin(lat1) * sin(lat2) +  cos(lat1) * cos(lat2) * cos(lon2 - lon1)) * R
R(지구반지름) = 6371000m

 

자전 및 공전을 하는 지구가 완전한 구체는 아니라서(WGS-84 ellipsoid 모델에서 적도 반지름은 6378km, 극 반지름은 6357km) Haversine식은 0.3% 가량의 오차를 가진다(즉, 1km 거리에서 3m 가량의 오차). 오차를 줄이려면Vincenty 공식을 사용하면 되지만(오차가 무려 1mm), 좌표간 거리가 대단한 정밀도를 가질 필요가 없고 계산량도 많아 Haversine 공식을 사용하기로 했다. 지구면의 두 지점에 관한 여러 종류의 수식은 여기를 참조.

 

왼쪽: 작업 전, 오른쪽: 작업 후

국도 및 지방도에 좌표 거리를 25m 이내로 했다. 선로 합치기는 비교적 잘 되었다. 예상대로 파일 크기는 1/2로 줄었다.

 

3. 인접 차선 결합

 

국가교통정보센터의 도로는 대부분 잘게 토막나 있다. 이것들을 하나의 긴 도로로 묶어 놓는다. 선로 합치기와 마찬가지로 벡터의 시작점/끝점이 일정 거리 내에 인접해 있으면 합친다. 경우의 수는 네 가지다.

 

vector1: B1 ---->---->---->---- E1
vector2: B2 ---->---->---->---- E2

 

E1, B2 인접: vector2의 B2를 지우고(공통) vector2를 vector1뒤에 연결한 다음 vector2 삭제(공통)
B1, E2 인접: 위와 비슷
B1, B2 인접: vector2의 방향을 바꿔 vector1에 연결
E1, E2 인접: 위와 비슷.

왼쪽: 떨어져서 별개로 선택되는 도로들, 오른쪽: 연결된 도로

4. 벡터 평탄화

 

벡터가 완만하게 변화하는 구간에서만 평탄화를 실시한다. 평탄화를 하면 지점의 갯수가 줄어든다. 변화폭은 일정 갯수의 지점에서 벡터 방향(bearing)의 각도 변화의 총합을 측정하여 파라미터로 결정한다. 예: 10점 동안 각도 변화가 10도 미만이면 평탄화.

평탄화는 급격한 변화 구간(예를 들면 도로의 인터체인지의 원형 진입로)은 그대로 놔 두고 완만한 변화 구간을 압축하는 효과가 있어 파일 크기를 줄일 수 있다. 3점 각도 변화 15도에서 원래 데이터를 크게 손상시키지 않으며 그 크기를 대략 1/2로 줄였다.

잘 뒤져보면 위에 언급한 네 가지 작업에 적합한 좋은 알고리즘이 있을 법하지만, 굳이 속도나 효율을 요구하는게 아니라서 일단 이대로 프로그래밍하기로 했다.

네 가지 작업의 효과(또는 취약점)을 확인할 수 있도록 파라미터를 선택하여 처리할 수 있도록 GUI 프로그램을 짰다.

이전의 해안선 처리에 사용하던 프로그램과 통합하고, polygon 처리할 때도 smoothing을 넣으니까 vertex 수가 10% 가량 줄었다. 빙고.

프로그램을 돌려보니 출력 파일의 크기가 원본의 24~30% 정도 밖에 되지 않는다. 원래 계획했던, 350MB를 100MB 이내로 줄이는 것이 가능하다. 한동안은 한가하게 출력 파일을 검증하는 작업을 할 것이다. 더디긴 하지만 작업은 계획대로 진행되고 있다.

 

 

OSM의 coastline은 PGS에서 자동 변환한 것으로, polygon이 아닌 polyline이라서 OSM 파일을 Garmin GPS용 .img로 변환하면 바다와 육지 경계만 표시될 뿐, 바다는 파랗고 육지는 육지색으로 표시되지 않는 문제가 있다. 여기에다가 transparent 옵션을 주고 .img 파일을 만들면 해안선은 엉성한 basemap에 묻혀 사실상 GPS에서 보이지 않는다. basemap에 묻히면 해안선 뿐만 아니라 시퍼런 바다에 등고선이 돋아있는 희안한 현상도 볼 수 있다.


남해안. 아래쪽은 새파란 바다여야 한다.

OSM의 랜더러들은 해안선을 판단해 바다와 육지 경계를 표시하므로 문제가 없지만, mkgmap, cgpsmapper 등에서 coastline을 적절히 처리할 방법을 찾아야 한다. mkgmap 제작자가 해안선을 처리할 방법을 나름대로 강구중인 것 같으나, 언제될 지 알 수 없다. 주욱 OSM 지도 작업을 하면서 느낀 거지만, 일단은 자력갱생 해야 한다.

현재로선 OSM 지도의 주 활용처가 Garmin GPS 디바이스인 관계로 이 문제를 해결하기 위해 임시변통 격으로(mkgmap에서 해안선을 다루게 될 때까지) polygon으로 구성된 해안선 데이터를 얻어와 바다를 구성하는 수 밖에 없다.

세계자전거오지여행의 운영자가 국가교통정보센터에서 공개한 한반도 도로 지도를 가공한 mp 파일에는 해안선 polygon 데이터가 있다(최근 데이터에는 없음). 이것을 뜯어와 img를 만들어 봤더니 드디어 바다와 육지가 구분되었다. 그런데 문제가 좀 있다.

해안선을 확대하면 디지타이즈된 라인이 아마도 한반도 저해상도 파일을 가지고 작업한 것처럼 라인이 단속적으로 끊어졌다.

성격은 좀 다르지만 국가교통정보센터의 도로에도 비슷한 문제가 있다.

오버줌 상태에서 보면 라인이 매끄럽게 연속적으로 이어진 것이 아니라 단락적으로 지글거린다. 이 때문에 도로만 데이터 크기가 350MB나 되어 OSM 서버에 도로 데이터를 올리는 것을 계속 망설이고 있었다. 해안선 문제 뿐만 아니라, 불필요한 노이즈를 억제하여 도로 데이터의 크기를 줄이는 방법도 역시 필요하다.

이들 문제를 해결하기 위해 벡터를 보간하는 방법을 궁리했다. 첫번째 방법은 벡터의 연속된 인접 3점으로 구성된 삼각형의 무게 중심을 구하는 공식을 이용해 다음 지점을 구하는 것이다. 그림으로 그리면,

좌상단이 첫 포인트, 우상단이 마지막 포인트. 두꺼운 실선은 원래의 좌표 벡터. 푸른선은 3점 포인트를 이용해 구성된 삼각형, 붉은선은 ABC 세 점으로 이루어진 연속적인 삼각형의 중심을 따라 이어진 선.

이렇게 하면 anti alias를 한 것 같은 효과가 나타난다. 원래의 검정 실선보다 붉은 선이 좀 더 부드러워지는 것. 포인트 수는 원래 7개에서 5개로 줄었는데, polyline이 상당히 길 경우 이 알고리즘을 적용하면 포인트의 수가 대략 1/2로 감소한다. 즉 파일 크기가 50% 줄어든다. 식: 각 점 A,B,C가 (x1,y1), (x2,y2),  (x3,y3)로 구성되어 있을 때 삼각형의 중심은 xc = (x1 + x2 + x3) / 3, yc = (y1 + y2 + y3) / 3.

두번째 삼각형과 세번째 삼각형에서 지나친 평균화 때문에 라인이 뭉툭해진다. 이것을 개선하기 위해 알고리즘을 바꿔 마름모의 중심을 이용하되 마지막 3점은 삼각형으로 계산하는 형태로 바꿨다.
 
삼각형보다는 원래 형태가 잘 반영되었다. 사각형(마름모) ABCD의 중심점은 xc = (x1 + x2 + x3 + x4) / 4, yc = (y1 + y2 + y3 + y4) / 4. 이것을 일반화하면 xc = (x1 + x2 + ... + xn ) / n, yc = (y1 + y2 + ... + yn) / n. 어? 정리해 놓고 보니 좌표의 이동 평균을 취한 것이다. 기껏 머리 좀 굴렸더니만 하하. 이때 n이 증가하면 포인트 수는 줄어드는 대신 궤적의 변화는 둔해진다.

과문한 탓에 더 빠르고 훌륭한 알고리즘을 만들진 못했다. 어디 뒤져보면 꽤 괜찮은 알고리즘이 있을 것 같다. 이쯤 해두고 알고리즘을 구현했다. 간단한 프로그램을 작성해 원래의 해안선 Polished file을 변환해 보았다.

원본.

알고리즘 적용.

만들긴 했지만 해안선은 애당초 DEM처럼 정밀한 데이터가 필요하다. 운 좋게 그걸 구하게 되면 이건 폐기할 것이다. 일단은 이동평균법을 이용해 Garmin GPS용 해안선을 예쁘게 만들었고(OSM 서버에는 안 올린다), 주 활용은 도로 데이터의 크기를 줄이는 것이다.

일단 OSM 서버에 도로 업로드하는 것은 당분간 안 하기로 했다. 도로 작업은 앞으로 몇 단계로 나뉘어 진행하게 될 것 같은데... 1. interpolation. 2. merge lanes. 3. connect. 4. add route information.

어쩐지 삽질하는 느낌을 지울 수 없다. 이런 데이터는 어딘가 이미 있다. 단지 그 데이터를 구하지 못해서...