방법
벡터를 래스터로 변경하는 다양한 방법이 존재한다. 그 중 TIN을 사용하기로한다. 사용할 유틸리티에서의 TIN은 입력 벡터 자료가 POLYGON, MULTIPOLOGON, LINE, LINESTRING 등의 관계없이 지오메트리를 구성하는 '점' 들에 대해 Delaunay triangulation method를 적용한다. DTM은 주어진 모든 점들을 이용하여 삼각형의 집합을 만들어 내는데, 내각의 최소값이 최대가 되도록 하는 방법이다.(시각적으로는 찌그러진 삼각형이 최소가 되도록 또는 가장 이쁜?? 삼격형 집합이 되도록) 이렇게 만들어진 각각의 삼각형을 Barycentric Interpolation method로 래스터화 하면 DEM 생성이 완료된다. 어떤 장단점이 있는지 잘 모르겠지만 IDW, Spline, Kriging 등에 비해 연산량은 적을 것으로 보인다. 다른 알고리즘 대비 불연속성(인접 삼각형이 공유하는 각 변에서 두드러지는 경계 현상)이 특징일 것으로 보인다.
입력자료
우선 DEM 생성에 필요한 자료는 해당 지역에 대한 3차원 점들의 집합니다. 3차원의 점은 위경도(x,y)해발고도(z) 값을 가지며, 많으면 많을수록 좋다. 가장 좋을 최소 조건은 각 격자 마다 최소 하나의 점들이 포함되는 것이다. 수평해상도가 작으면 가능할 지도 모르지만 100m 급 이하로 올라가면 앞서 언급한 공간보간법을 적용해야 할 것이다.
국가 공간정보 포털의 '연속수치지형도(http://market.nsdi.go.kr/goods/detail.do?gno=14943)' 항목을 보면, 다음과 같은 자료를 DEM 생성에 이용할 수 있을 것으로 확인된다.
1. 등고선(5,000:1데이터라고 한다. 이 등고선은 어떻게 만드는 걸까?? 이 등고선을 만들때 DEM이나 DSM을 소스로 만들었을 것 같은데, 그렇다면 등고선으로 또다시 DEM을 제작하는 것은 몇 가지 특별한 경우를 제외하면, 옳바르지 않은 방법이다. 가능한 경우 DSM/DEM을 직접 이용하면 좋을텐데..)
2. 표고점('점' 데이터이다. 사람이 직접 돌아다니면서 기록한건가??)
3. 해안선(해안선 데이터가 필요한 이유는, 등고선 데이터의 특이성?? 때문이다. 등고선 데이터는 5m의 수직 해상도를 갖고 있는데, 동해 쪽의 바다와 최인접한 등고선을 보면 0으로 끝나는 것이 아니라 5로 끝난다. 동해 쪽이라 절벽이 많아 0으로 끝나지 않는 경우가 있을 수 있을 것 같다. 하지만 모래사장과 같이 바다가 인접한 곳도 0으로 끝나지 않음은 이해되지 않는다.)
입력자료 보정
입력자료는 몇 가지 보정이 필요하다. 가장 큰 문제는 국가공간정보 포털에서 제공하는 연속수치지형도의 오류이다. 명백히 0또는 마이너스 고도가 아닌 지점임에도 불구하고, 0또는 마이스너 고도값을 갖는 등고선이 존재한다. 이를 제거해야 하는데 답을 정할 수 없다. 왜냐하면 실제로 0이거나 마이너스인 고도가 존재하기 때문이다.(해안 부분근에 지하주차장 공사하는 곳은 마이너스 고도를 갖는다) 따라서 첫째 내륙에는 0이하의 고도가 없음을 가정한다. 둘째는 주변과 특이적으로 고도차이를 보이는 등고선 들이 존재한다. 따라서 둘째, 인접지점과의 고도차이에대해 편차가 1% 이상은 곳은 제거한다. (편차 1%이상은 지점은 인접 30미터에 대해 고도차이가 50미터 이상은 지점들로 측정되었다.)
입력자료로 사용될 벡터자료 중
등고선은 LINESTRING형식이며, CONT라는 필드에 5m 수직해상도의 고도값이 포함되어 있다.
표고점은 POINT형식으며, NUME라는 필드에 1m 보다 더 상세한 수직해상도의 고도값이 포함되어 있다.
해안선은 LINESTRING형식이며, 고도값일 포함하는 필드는 없지만, 해안과 인접해 있으므로 고도를 0으로 가정할 수 있다.(절벽은 어쩌고?? 무시하자.ㅠ.ㅠ)
Rasterization에 사용할 GDAL 유틸리티는 Z필드 헤더가 일치해야한다. 따라서 등고선의 CONT필드를 기준으로, 표고점파일과 해안선 파일을 수정한한다.
표고점 파일은 CONT필드를 새로 만들고, NUME필드의 값을 복사한다.(또는 NUME필드 명을 CONT로 변경한다.)
해안선 파일은 CONT필디를 새로 만들고, 값을 0으로 채운다.
이제 이 3가지 파일의 하나의 GIS파일처럼 보이게 해야하는데, 여러가지 방법 중 GDAL의 VRT(http://www.gdal.org/drv_vrt.html)를 이용한 방법이 손이 덜 가는것 같다.
merge.vrt라는 파일을 아래와 같은 내용을로 만들어 보았다.
<OGRVRTDataSource>
<OGRVRTUnionLayer name="unionLayer">
<!-- 표고점 -->
<OGRVRTLayer name="N3P_F0020000"><SrcDataSource>refined_inputs/elevation_point.add_CONT</SrcDataSource></OGRVRTLayer>
<!-- 해안선 -->
<OGRVRTLayer name="N3L_E0080000"><SrcDataSource>refined_inputs/coastline.add_CONT</SrcDataSource></OGRVRTLayer>
<!-- 등고선 -->
<OGRVRTLayer name="N3L_F0010000_01"><SrcDataSource>refined_inputs/isoline/01</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_02"><SrcDataSource>refined_inputs/isoline/02</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_03"><SrcDataSource>refined_inputs/isoline/03</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_04"><SrcDataSource>refined_inputs/isoline/04</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_05"><SrcDataSource>refined_inputs/isoline/05</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_06"><SrcDataSource>refined_inputs/isoline/06</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_07"><SrcDataSource>refined_inputs/isoline/07</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_08"><SrcDataSource>refined_inputs/isoline/08</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_09"><SrcDataSource>refined_inputs/isoline/09</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_10"><SrcDataSource>refined_inputs/isoline/10</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_11"><SrcDataSource>refined_inputs/isoline/11</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_12"><SrcDataSource>refined_inputs/isoline/12</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_13"><SrcDataSource>refined_inputs/isoline/13</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_14"><SrcDataSource>refined_inputs/isoline/14</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_15"><SrcDataSource>refined_inputs/isoline/15</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_16"><SrcDataSource>refined_inputs/isoline/16</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_17"><SrcDataSource>refined_inputs/isoline/17</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_18"><SrcDataSource>refined_inputs/isoline/18</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_19"><SrcDataSource>refined_inputs/isoline/19</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_20"><SrcDataSource>refined_inputs/isoline/20</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_21"><SrcDataSource>refined_inputs/isoline/21</SrcDataSource></OGRVRTLayer>
<OGRVRTLayer name="N3L_F0010000_22"><SrcDataSource>refined_inputs/isoline/22</SrcDataSource></OGRVRTLayer>
</OGRVRTUnionLayer>
</OGRVRTDataSource>
중요한 점은 POINT형식의 자료인 '표고점'이 맨 앞에 기술되어야 한다는 것이다.
아래의 명령으로 대한민국 전체에 대한 DEM파일을 생성 할 수 있다고 생각했으나...
점들이 너무 많은 이유로 오류가 발생(512GB 램을 가진 서버에서도 실패, 실제 메모리 사용량을 훨씬 적었지만..)했다.
특정 크기로 타일링하여, 만들어진 타일들을 이어 붙이는 방식을 선택해야 했다. 여기서 중요한 점은 타일의 경계면에서 불연속성이 발생하지 않도록 하는 점이다.
gdal_grid -ot Float32 -a_srs EPSG:5179 -a linear -zfield CONT -outsize 300 300 -spat 884450 1467570 893510 1476630 -txe 884480 893480 -tye 1467600 1476600 -co tiled=yes -co compress=deflate -co BIGTIFF=YES -where 'CONT>=0' ./prepared_inputs/merge.vrt ./output/012.001.tiff
위 명령은 가로영역(884480 ~ 89348), 세로 영역(1467600 ~ 1476600)에 대해 300x300크기의 타일을 생성한다.(영역의 좌표는 EPSG:5179이다.) 그런데 -spat로 입력된 영역을 보면 타일의 영역보다 조금 더 크게 버퍼를 두었다. 이렇게 버퍼를 두면 타일을 합칠 때 발생하는 불연속을 제거 할 수 있다.
이 타일을 전국에 대해 제작 후 병합하면 되겠다.
다수의 코어를 이용하여 타일일 생성하는 코드이다.
32코어 64스레드를 이용하여 마스킹 처리 포함 30미터 해상도 기준 일 주일 정도 소요 되었고, 용량을 Float32, 256x256 deflate compression으로 350MB이다.
10미터 해상도로(Float32, 256x256 Deflate compression) 3.2GB 이다.
1미터 해상도로 312GB이다.(마스킹 작업 중에 있으나 2주일 째 0%에서 올라가지 못하고 있다. 5,000:1 원본자료로는 5m 격자 해상도가 의미 있다고 한다.-출처제시필요)
import math import multiprocessing import os import subprocess RESOLUTION = 1 TILE_SIZE = 300 EXTENT_SIZE = TILE_SIZE * RESOLUTION BUF_SIZE = 5 OUTPUT_DIR = './output/1m' EXTENT = (776480, 1458600, 1387950, 2044840) LEFT, BOTTOM, RIGHT, TOP = EXTENT X_TILES = int(math.ceil(float(RIGHT - LEFT) / EXTENT_SIZE)) Y_TILES = int(math.ceil(float(TOP - BOTTOM) / EXTENT_SIZE)) W = X_TILES * EXTENT_SIZE H = Y_TILES * EXTENT_SIZE def make_grid(x_tile_index, y_tile_index): extent_x = x_tile_index * EXTENT_SIZE + LEFT extent_y = y_tile_index * EXTENT_SIZE + BOTTOM l_ = extent_x r_ = extent_x + EXTENT_SIZE b_ = extent_y t_ = extent_y + EXTENT_SIZE sl_ = l_ - BUF_SIZE sr_ = r_ + BUF_SIZE sb_ = b_ - BUF_SIZE st_ = t_ + BUF_SIZE cmd = [ 'gdal_grid', '-q', '-ot', 'Float32', '-a', 'linear', # '-a', 'linear:radius=5000', '-a_srs', 'EPSG:5179', '-zfield', 'CONT', '-outsize', str(TILE_SIZE), str(TILE_SIZE), '-spat', str(sl_), str(st_), str(sr_), str(sb_), '-txe', str(l_), str(r_), '-tye', str(t_), str(b_), '-co', 'tiled=yes', '-co', 'compress=deflate', '-co', 'BIGTIFF=YES', # '-where', 'CONT>0 or DIVI in (\'CLD000\', \'CLD001\', \'CLD002\')',
'-where', "CONT>=0", './prepared_inputs/merge.vrt', os.path.join(OUTPUT_DIR, '{:03d}.{:03d}.tiff'.format(x_tile_index, y_tile_index)), ] print(' '.join(cmd)) subprocess.call(cmd) def main(): if not os.path.exists(OUTPUT_DIR): os.makedirs(OUTPUT_DIR) tasks = [] for y_tile_index in range(Y_TILES): for x_tile_index in range(X_TILES): tasks.append( (x_tile_index, y_tile_index) ) print('task count: {}'.format(len(tasks))) pool = multiprocessing.Pool() for i, task in enumerate(tasks): pool.apply_async( make_grid, args=task ) print('{}/{}'.format(i + 1, len(tasks))) pool.close() pool.join() if __name__ == '__main__': main()
이렇게 만들어진 타일을 gdal_merge.py 유틸리티로 합치면 되겠다.