[GIS] 당신이 겪을 수 있는 georeference 문제들과 내 해결법

전제 및 가정

지오레퍼런싱(georeference)을 하기 위해서는 다음 세 가지 정보가 꼭 필요하다.

1. 좌표계 정보(EPSG 코드가 명시되어 있으면 가장 좋다.)
   1. 한국에서 쓰는 좌표계는 몇 가지가 있는데, 중부 원점(5186)을 가장 많이 쓰고 조선소의 경우 동부 원점으로 기술된 경우가 있다.
   2. 한국 내에서 영상의 지리 정보를 기술할 때 표준과 상관없이 기술하는 경우가 대다수이므로, 이럴 때는 다음을 확인해 보자. 한국에서 정의된 좌표계는 좌표계의 원점을 (0,0)으로 뒀을 때, 이를 좌표계 상에서 동쪽으로 얼마, 북쪽으로 얼마 이동한 가상 원점을 사용한다. 예를 들어 동쪽으로 200,000m 북쪽으로 600,000m 이동했다고 적혀있으면 epsg:5186에 해당하는 중부 원점이다.
2. 좌표계에서의 영상의 좌상단 좌표
   1. gsd : pixel 2 geo 할 때 1pixel이 좌표계 상에서는 얼마에 해당하는지 그 비율을 서술함

1. 이미지를 georeference 하기

목적

• 아예 지리 정보가 없는 이미지를 지도 위에 놓기
  • 예 : inference 생성 결과는 단순 이미지일 뿐, 좌표계 정보가 없다.
• 파이썬으로 처리할 수 있는 방법이 없을까?

원리

• GCP를 지정해서 지리정보가 없는 이미지와 그 이미지를 놓을 위치 간의 참조를 만든다.
• 평면 대 평면을 매핑하면 된다.
• 특정 지역을 매핑하는 것이 아니라면, 이미지의 가장자리 4개 꼭지점과 이미지를 놓을 지도 위 영역 가장자리 4개 꼭지점을 매핑한다.
• 매칭 순서
  • 최소 점 3개가 필요하다(평면을 이루는 최소 단위)
  • 픽셀 좌표계 x좌표 ⇒ longitude
  • 픽셀 좌표계 y좌표 ⇒ latitude
  • ccw인지 cw인지는 상관없고 하나로 통일해서 매칭하면 된다
• 주의할 점
  • 두 좌표계 간 원점이 다르다 : 픽셀 좌표계는 원점이 좌상단, 위경도 좌표계는 원점이 중앙에 위치한다. (원점은 0,0)
  • 경도 값은 중간을 0으로, 북반구가 양의 값으로 90까지 증가하고 남반구가 음의 값으로 -90까지 감소한다.
  • 픽셀 좌표계는 상단에서 하단으로 값이 증가한다.
  • 위 두 사항 때문에 좌표값을 매칭할 때 y값과 경도 값을 매칭하는 일이 헷갈린다. 아래 그림에서도 볼 수 있듯 좌상단 꼭지점의 y값은 '값'으로 봤을 때 max값이다.

- 참고 자료 : [https://www.earthdatascience.org/courses/use-data-open-source-python/intro-raster-data-python/fundamentals-raster-data/raster-metadata-in-python/](https://www.earthdatascience.org/courses/use-data-open-source-python/intro-raster-data-python/fundamentals-raster-data/raster-metadata-in-python/)

구현 코드

from osgeo import gdal, ogr, osr

def georef_inference(raster_path, infer_path) :
    '''
    raster_path : 대상이 되는 원본 영상의 경로
    infer_path : 생성한 inference의 경로
    '''

    '''
    원본 영상에서 inference에 매칭할 기준점들 추출(가장자리 4개 점)-> epsg코드 파싱으로 얻기 -> 원본 영상이랑 inference랑 대응되는 점들 gdal.GCP 생성으로 mapping
    -> SetGCP로 inference georeference 하기
    '''
    origin_ds = gdal.Open(raster_path)
    geoTransform = origin_ds.GetGeoTransform()
    minx = geoTransform[0]
    maxy = geoTransform[3]
    maxx = minx + geoTransform[1] * origin_ds.RasterXSize
    miny = maxy + geoTransform[5] * origin_ds.RasterYSize
    origin_points = [minx, miny, maxx, maxy]

    crs_str = origin_ds.GetProjectionRef()
    res = crs_str.replace("[",',').replace("]",',').replace('\"',',').split(',')
    num_list = [x for x in res if x.isdecimal() == True]
    crs_num = num_list[-1]

    target_ds = gdal.Open(infer_path)
    sr = osr.SpatialReference()
    sr.ImportFromEPSG(int(crs_num))
    target_maxx = target_ds.RasterXSize
    target_maxy = target_ds.RasterYSize
    target_points = [0, 0, target_maxx, target_maxy]

        '''
        geotransform[0] = top left x
        geotransform[1] = w-e pixel resolution
        geotransform[2] = 0
        geotransform[3] = top left y
        geotransform[4] = 0
        geotransform[5] = n-s pixel resolution (negative value)
        '''

    gcps = [gdal.GCP(origin_points[0],origin_points[3],0,0,0),
            gdal.GCP(origin_points[0],origin_points[1],0,0,target_points[3]),
            gdal.GCP(origin_points[2],origin_points[1],0,target_points[2],target_points[3]),
            gdal.GCP(origin_points[2],origin_points[3],0,target_points[2],target_points[1])
           ]
    target_ds.SetGCPs(gcps, sr.ExportToWkt())
    target_ds.SetGeoTransform(origin_ds.GetGeoTransform())
    target_ds = None

2. 처음부터 georeference 하기

목적

• 이미지 생성시 georeference하는 방법
• 다른 이미지를 읽어와서 생성하는 게 아니라, array에 담긴 데이터를 georeference 정보를 담아서 저장

원리

코드

'''
    original_img_path : georeference 정보를 가져올 이미지의 경로
    img_write_path : 이미지를 저장할 경로
  im_array : numpy array로 저장된 데이터
'''
def save_img_with_georef(original_img_path, img_write_path, im_array) :
    original_img = gdal.Open(original_img_path)
        # 이미지는 c,h,w의 형태로 들어오는 것을 기대한다. 
    num_channel = im_array.shape[0]
    nrows = im_array.shape[1]
    ncols = im_array.shape[2]

        # gdal에서 파일 생성(create)이 지원되는 드라이버는 geotiff가 대표적이다.
    driver = gdal.GetDriverByName('GTiff')
    output_img = driver.Create(img_write_path, ncols, nrows, num_channel ,gdal.GDT_Int16, ['COMPRESS=LZW'])

    output_img.SetGeoTransform(original_img.GetGeoTransform())
    for i in range(1, num_channel+1) :
        output_img.GetRasterBand(i).WriteArray(im_array[i-1])

    output_img.SetProjection(original_img.GetProjection())
    output_img = None

3. 이미지 좌표 → 지도 좌표

목적

패치 단위로 이미지를 읽어와 얻어낸 이미지 좌표를 지도 좌표로 바꾸고자 할 때

원리

패치의 x, y offset만큼 평행이동하고 (픽셀 : 실제 거리) 간의 비율인 gsd를 곱해 이미지의 크기를 조정한다.

gdal에서 geotransform을 얻을 때 좌상단의 x/y좌표, 이미지의 가로/세로 gsd를 구할 수 있다.

이를 이용해서 변환에 필요한 행렬을 계산한다.

코드

Affine transformation

'''
from shapely import affinity 
from shapely.geometry import Polygon
from shaeply.geometry import Point
'''

def canvas2WorldCoord(geom, geo_transform,xoff=0, yoff=0) :
    tl_x, x_gsd, _, tl_y, _y, y_gsd = geo_transform

    geom = affinity.affine_transform(geom, [1, 0, 0, 1, xoff, yoff])
    return affinity.affine_transform(geom, [x_gsd, 0, 0, y_gsd, tl_x, tl_y])

4. 이미지 + gis 정보가 metadata로 따로 있을 때

목적

이미지를 전달 받았는데, 이미지는 이미지 대로 찍고 georeference 정보는 이미지와 별개로 다른 파일/시스템에 저장되어 있는 경우에 이미지에 georefence 적용하기

원리

이 글의 대전제 부분에도 명시했지만, 이미지의 좌상단 점을 기준으로 원본 좌표계에 먼저 자리를 잡고 gsd에 따라 영상의 비율을 수정하는 작업을 거친다. 좌표계를 모른다면 답이 없다.

코드

from osgeo import gdal, osr
from typing import List

def set_geo_reference(input_path : str , output_path : str , epsg_code : int , geo_transform : List[float]) : 
    input_ds = gdal.Open(input_path, gdal.GA_Update)
    input_ds.SetGeoTransform(geo_transform)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg_code)
    input_ds.SetProjection(srs.ExportToWkt())
    input_ds = None

def get_epsg_code(path : str) -> str :
    ds = gdal.Open(path)
    prj = ds.GetProjection()
    prj = osr.SpatialReference(wkt=prj)
    return int(prj.GetAuthorityCode(None))

def set_reproject(input_path : str, output_path : str , target_epsg_code : int) : 
    cmd_dict = {'input_path' : input_path, 'output_path' : output_path, 'source_srs' : f'{get_epsg_code(input_path)}', 'target_srs' : f'EPSG:{target_epsg_code}'}
    base_command = 'gdalwarp {input_path} {output_path} -s_srs {input_srs} -t_srs {output_srs}'
    process = subprocess.run(base_command.format(input_path=cmd_dict['input_path'],output_path=cmd_dict['output_path'], input_srs=cmd_dict['source_srs'], output_srs=cmd_dict['target_srs']), shell=True, check=True)

설명

1. set_geo_reference
2. 인자로 받아들인 epsg code와 geo_transform을 통해 georeference를 하는 함수
3. get_epsg_code
4. georefence가 된 영상의 epsg code를 얻어내는 함수
5. set_reproject
6. georeference가 된 영상을 다른 좌표계로 reprojection 시키는 함수

환경 설정

혹시 본인이 가지고 있는 docker에서 gdal 2.4.0을 사용하고, pyproj 3.2 미만 버전을 사용한다면

error 6: unable to load proj.4 library (libproj.so), creation of ogrcoordinatetransformation failed.

라는 에러를 볼 수 있다. 이 때 해결법은 더 높은 버전의 pyproj를 깔아보고, 그래도 안되면 다음을 실행한다. proj의 개발 버전을 까는 명령어이다.

sudo apt update
sudo apt install libproj-dev

5. 아예 아무 정보도 없을 때

목적

crs(srs)정보가 누락된 경우, 이를 원본 영상에 할당하고 싶을 때.

코드

gdal_edit.py -a_srs EPSG:4326 {파일 이름}