온맵 pdf 를 geotiff 로 변환하기

[날틀_김상용]

순서는 다음과 같습니다.

1. 국토지리정보원에서 온맵 다운받기

2. 파일명에서 한글 삭제(프로그램 renamer 추천)

3. 온맵 .pdf  불필요한 레이어 및 워터마크 삭제 (Acrobat pro 추천)

4. pdf 를 래스터 이미지로 변환 (pdftoolbox 추천)

5. 맵매칭에 필요한 Ground Control Point 만들기 (python)

6. 맵매칭으로 gcp 포인트 수정 (qgis 추천)

7. 래스터 이미지와 gcp 포인트를 이용해서 geotiff 만들기 (python)

아시다시피 글로벌 매퍼를 사용하면 작업단계가 줄어들고 시간도 절약됩니다.

14일 동안 사용할 수 있는 임시라이센스를 받을 수도 있습니다.

Acrobat pro 와 pdttoolbox 도 상용프로그램입니다만, 시험버전으로도 30일, 14일동안

모든 기능을 사용할 수 있습니다.

3. 온맵 .pdf  불필요한 레이어 및 워터마크 삭제 (Acrobat pro 추천)

어크로뱃 프로의 일괄작업 중에 워터마크 삭제 기능이 있습니다.

항공사진의 경우 아래 자바스크립트를 일괄작업에 만들어 넣으면 됩니다,.

var oLayers = this.getOCGs();
for (var i=0; i<oLayers.length; i++) {
     if (oLayers[i].name == "항공사진") oLayers[i].initState = false;
}

4. pdf 를 래스터 이미지로 변환 (pdftoolbox 추천)

어크로뱃 프로에서 일괄로 이미지로 내보내기 해도 되지만 안티얼리어싱 기능이 없어

이미지 품질이 떨어집니다.

pdftoolbox cil 버전을 사용해서 배치파일을 만들면 고품질의 이미지로 변환할 수 있습니다.

pdftoolbox --saveasimg --imgformat=png --resolution=350 33602.pdf

5. 맵매칭에 필요한 Ground Control Point 만들기 (python)

맵매칭은 여러가지 방법이 있지만 여기서는 qgis 를 이용합니다.

qgis 설치 때 같이 깔린  파이썬을 이용해서 맵매칭에 필요한 대략적인 gcp 포인트를 만들어줍니다.

아래가 5만 축척 코드 예제입니다.

50k_coord.py

import glob
import sys

file_list = glob.glob("M:\\onmap_mod\\50k_350dpi\\*.png")

for file in file_list:
 print file
 code_no = file[-9:-4]
 y = float(code_no[:2])
 x1 = float(code_no[2])
 x2 = float(code_no[3:])
 point_file = "M:\\onmap_mod\\50k_350dpi\\" + code_no + ".points"
 f = open(point_file, 'w')
    
 if x2 <= 4:
  North = y + 0.25 * 3
 elif 4 < x2 <= 8:
  North = y + 0.25 * 2
 elif 8 < x2 <= 12:
  North = y + 0.25
 else:
  North = y
  
 if (x2 % 4) == 1:
  East = 120 + x1
 elif (x2 % 4) == 2:
  East = 120 + x1 + 0.25
 elif (x2 % 4) == 3:
  East = 120 + x1 + 0.25 * 2
 else:
  East = 120 + x1 + 0.25 * 3
  
 East_list = [East, (East+1/3.0*0.25), (East+2/3.0*0.25), (East + 0.25)]
 North_list = [North, (North+1/3.0*0.25), (North+2/3.0*0.25), (North + 0.25)]
  
 f.write("mapX,mapY,pixelX,pixelY,enable\n")
  
 j = 0
 while j <= 3:
  i = 0
  while i <= 3:
   East_x = East_list[i]
   North_y = North_list[j]
   pnt_x = 670.0 + 2050.0 * i
   pnt_y = -8500.0 + 2550 * j
   pnt = [str(East_x), str(North_y), str(pnt_x), str(pnt_y), str(1)]
   data = ",".join(pnt)
   f.write(data + "\n")
   i = i + 1
  j = j + 1
 f.close()
 

6. 맵매칭으로 gcp 포인트 수정 (qgis 추천)

qgis 의 지오레퍼런스를 실행시킨 후 도곽별로 래스터 이미지와 gcp 포인트를 불러온 후

gcp 포인트를 수정후 gcp 포인트를 저장해 줍니다.

7. 래스터 이미지와 gcp 포인트를 이용해서 geotiff 만들기 (python)

5만 축척 도엽 코딩예제입니다.

onmap50k_match.py

import os
import glob

file_list = glob.glob("d:/50k/*.points")

for file in file_list:
 
 fdir = file[:-13]
 fname = file[-12:-7]
 
 print (file + ',' + fdir + ',' + fname)
 
 f = open(file,'r')
 lines = f.readlines()
 
 f2 = open(fname + ".bat", 'w')
 f2.write("gdal_translate")
 j = 0
 for i in lines:
  if j == 0:
   pass
  else:
   gcp = i.split(',')
   gcps = gcp[2] + ' ' + str(abs(float(gcp[3]))) + ' ' + gcp[0] + ' ' + gcp[1]
   f2.write(' -gcp ' + gcps)
  j = j + 1
 tImg1 = fdir + '/' + fname + ".png"
 tImg2 = fdir + '/' + fname + "t.png"
 tImg3 = fdir + '/' + fname + ".tif"
 f2.write(' "' + tImg1 + '" "' + tImg2 + '"' + "\n")
 f2.write('gdalwarp -t_srs EPSG:4326 -r cubic -tps -co COMPRESS=LZW "' + tImg2 + '" "' + tImg3 + '"')
   
 f.close()
 f2.close()

---

테스트는 해봤지만 저도 작업중(4단계까지 완성)입니다.

첨부파일첨부된 파일이 2개 있습니다.

• etc 파일 50k_coord.py 다운로드
• etc 파일 onmap50k_match.py 다운로드

[남정권]

좋은 정보 감사합니다.

[산하]

비슷한 과정을 거쳐본지라 그냥 경험적인 이야깁니다.
4. pdf 원본자체가 고품질 그래픽이 아니기때문에 jpg로 변환시 300DPI만해도 pc용으로도 차고넘치는 것 같습니다.
jpg변환시 품질만보자면 일러스트가 제일 나은듯합니다.
7. -t_srs EPSG:4326 품질저하가 발생하지 않던가요? (직사각형 원본이미지-> 정사각형으로 변환)
EPSG:3395 로 하면 원본 직사각형 그대로 만들 수 있는데....projection 변환을 한번더 해줘야합니다.
-r cubic 품질이 중요하다면 lanczos 와 비교해보고 결정하시는게 좋을듯합니다.
그외에 -te 옵션을 사용해서 도곽바깥쪽을 삭제할 수도 있습니다.

[날틀_김상용]

정보 공유 감사합니다.
350dpi 로 한 이유는 타일화할 때 줌레벨이 바뀌는 해상도가 3xx dpi 이였습니다.
말씀하신 것처럼 EPSG;4326 으로 하면 정사각형으로 변하는데 cubic 랜더링을 하면 품질저하는 못느낄정도입니다.
lanczos 과 비교해봐야겠습니다.
위 코드에는 안썼지만 gdal_translate로 외곽을 잘라주려고 했는데 -te 옵션이 있었네요 ^^

[날틀_김상용]

아래 그림을 보면, 글자에 관한 한 cubic 이 나은 듯합니다.
원본, cubic, lanczos 순서입니다.