목적

배경

• 공간 정보 데이터는 맞는데, 그렇다고 이를 바로 지도에 표출하려니 데이터가 다른 분야에서 다루는 데이터라서 공간 정보 데이터의 형식에 맞춰줘야 할때
• 기상 데이터들은 공간에서 일어나는 현상을 다룬 데이터라 위치 정보들이 있지만, 우리가 생각하는 EPSG:4326 혹은 EPSG:3857 같은 GIS 계에서 흔히 쓰는 좌표계를 쓰지 않는다
  • lambert conformal coordinate
  • 각 픽셀마다 좌표 정보가 있다지만, 다음의 두 가지 형식으로 저장된다
    • 픽셀 마다 좌표 정보가 있다. 각 픽셀들은 lcc 좌표계 내에서 등간격으로 존재하고 표기는 위경도로 표기된다
    • 상대적인 범위 내에 격자 형식으로 좌표를 기술할 수도 있다. 아래 그림은 기상청에서 쓰는 격자 좌표계를 기술한 문서의 일부
      

      

• 기상 데이터들은 정해진 데이터 모델에 위치 정보와 데이터를 넣어 그에 맞춰 시각화를 한다. 지구 위에서 보려면 별도의 처리 과정이 필요하다.
  • 데이터 모델을 통해 가시화 하면 이런 형태로 가시화 한 것을 볼 수 있다. 위성으로 측정한 지표면 절대온도 데이터
    

    

  • lcc에서는 분명 직사각형으로 보였는데 EPSG:4326에서는 이렇게 부채꼴로 보인다.
    

    

대상

• 기상 데이터와 같이 위치 정보는 있는데, 자주 쓰는 좌표계로 데이터를 georeference 해주고 싶을 때
• gcp를 이용해서(특히 네 구석의 점) 어떻게 georeference를 하는가? 코드로 하고 싶다
• npy로 출력한 이미지를 어떻게 지구상에 표출할까?

내용

상황

• 주어진 건, 모델의 예측 결과(npy)와 각 array의 요소들의 위경도 좌표 리스트
  • lat, lon 정보가 따로 따로 2차원 배열로 들어옴
• lcc 좌표계 내에서 정보는 아래와 같은 모양새로 분포하고 있다.
  

  

코드

작업할 좌표계에 맞춰 기준 좌표 준비하기

좌표는 위경도로 서술되어 있지만 이미지를 가장 쉽게 Georeference하는 방법은 직사각형의 이미지의 상하좌우 꼭지를 그대로 좌표계에 붙여넣는 것이다.

이미지의 회전, w-e/n-s 간의 gsd 차이 등을 신경쓰지 않아도 되기 때문이다.

좌표는 위경도 좌표계에서 부채꼴의 양상을 띄기 때문에, 만약 이를 그대로 Georeference하려면 이미지의 테두리에 있는 상당한 점들을 써야 한다.

그래서 위경도 좌표계로 서술한 점을 lcc 좌표계에 맞게 변환한다.

import pyproj

proj_102012 = "+proj=lcc +lat_0=0 +lon_0=105 +lat_1=30 +lat_2=62 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs"

# 1. 기준으로 삼을 네 구석의 점을 준비
corner_points = [(45.55741, 118.57613),
 (45.55741, 133.8462),
 (32.30801, 118.57613),
 (32.30801, 133.8462)]

# 2. 위경도 좌표계로 기술된 점을 원본 lcc 좌표계의 점으로 변환
lcc_corner_points = [pyproj.Proj(proj_102012)(x[1],x[0]) for x in corner_points]

gdal gcp 생성하기

from osgeo import gdal

# 1. 생성할 이미지의 사이즈를 지정
# width가 col, height가 row에 해당함
row = lat_array.shape[0]
col = lat_array.shape[1]

# 2. 맞춰줄 픽셀 좌표 정의
# 픽셀 좌표와 실 좌표를 어떻게 맞춰줄 지 정하려면 이미지의 어느 부분이 실제 어느 위치에 존재하는 지 이어줘야 한다
# 이미지의 원점은 좌상단. upper_left, upper_right, lower_right, lower_left
pixel_index = [(0,0), (0,col),(row,col),(row, 0)]
gcp_list = []

# 3. 픽셀 좌표와 맞는 실 좌표의 쌍을 만든다. 
for x in range(4) : 
    gcp_list.append((pixel_index[x][0],pixel_index[x][1],lcc_corner_points[x][0],lcc_corner_points[x][1]))
'''
gcp_list = 
[(0, 0, 1013848.356912874, 5479167.3746448895),
 (0, 576, 2116635.162965722, 5784374.250365057),
 (720, 576, 2629740.4931766586, 4450362.210115054),
 (720, 0, 1259620.9893719808, 4071168.434106591)]
'''  
# 4. gdal 내의 gcp class를 이용해서 gcp 객체들을 생성
gdal_gcp_ins_list = []

for item in gcp_list:
    pixel, line, x, y = map(float, item)
    z = 0
        # pixel, line이라는 단어가 등장하는데, col과 row에 해당한다고 생각하면 된다. 
    gcp = gdal.GCP(x, y, z, pixel, line)
    gdal_gcp_ins_list.append(gcp)

새로운 geotiff 생성

먼저 좌표계 정보를 알려주기 위해서 wkt string으로 표기된 ESRI:101210 좌표계의 정보를 가져왔다.

#lcc_wkt_str = 'PROJCS["Asia_Lambert_Conformal_Conic",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",105],PARAMETER["standard_parallel_1",30],PARAMETER["standard_parallel_2",62],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["ESRI","102012"]]'

이 부분은 gdal의 geotiff 드라이버의 영상 생성 기능을 이용해 geotiff로 영상을 저장하는 함수이다.

from osgeo import gdal, osr
'''
인자 설명
- gcps : gdal의 gcp 클래스 객체가 들어간다. 최소 세 개는 들어가야 한다(평면을 구성하는 최소 꼭지점의 갯수)
- img_path : 이미지를 저장할 경로
- band_num : rgb 저장이 필요하면 3
- radar_array : 저장하려는 numpy array
'''
def radar_to_geotiff(gcps, img_path, band_num, radar_array) : 
    format = "GTiff"
    # 1. 이미지 크기 정하기
        col = radar_array.shape[0]
    row = radar_array.shape[1]    
    driver = gdal.GetDriverByName( format )
    metadata = driver.GetMetadata()
    # 주의해야 할 점 1 : col, row의 순서다 보니 numpy array를 넣으려면 transpose가 필요하다. 
    # 2. geotiff를 생성하기 위해 data source 생성
        '''
        인자들을 설명한다. 
        - img_path : 저장경로
        - col : 영상의 width
        - row : 영상의 height
        - band_num : 당시 RGB로 저장했기에 3이 들어감
        - gdal.GDT_Float32 : 데이터 타입
        - options : 별도 옵션. 이 옵션을 정해주지 않으면 자꾸 qgis에서 gray scale로 파일을 읽어 성가셨음
        '''
        dst_ds=driver.Create(img_path, col, row, band_num ,gdal.GDT_Float32,options = [ 'PHOTOMETRIC=RGB' ])
    for i in range(1, band_num+1) : 
        dst_ds.GetRasterBand(i).WriteArray( np.transpose(radar_array[:,:,i-1]) )
                # nodatavalue를 (255,255,255)로 설정한다    
        dst_ds.GetRasterBand(i).SetNoDataValue(255)

        # 3. 영상의 좌표계를 설정한다. 우리는 ESRI:102012를 쓸 건데, osr의 ImportFromEPSG에 ESRI 코드도 들어감
        # 다만 GDAL 3.x 버전을 사용할 것
    srs_lcc = osr.SpatialReference()
    srs_lcc.ImportFromEPSG(102012)
        # 3-1. data source에 projection 설정
    dst_ds.SetProjection(srs_lcc.ExportToWkt())
        # 4. gcp 설정 
    dst_ds.SetGCPs(gcps, srs_lcc.ExportToWkt())    
        # 5. 4에서 입력한 GCP에서 geotransform을 계산해서 data source에 세팅하기
        '''
        gdal의 c++ 코드를 까면 대부분 GeoTransform이 data source에 설정되어 있다는 것을 가정하고 계산을 한다.
        gdal 함수 중 입력한 gcp에서 역으로 GeoTransform을 계산해서 설정하는 함수가 있다. 꼭 쓰자. 
        '''
    geotransform = gdal.GCPsToGeoTransform(dst_ds.GetGCPs())
    dst_ds.SetGeoTransform(list(geotransform))
    # 6. 설정 저장. data source의 포인터를 저장한 변수에 None을 넣으면 여태 입력한 정보를 영상에 넣어 영상 생성
        dst_ds = None

reproject

이 부분은 정말 아쉽게도 아직 코드로 구현하는 방법을 찾지 못해서, cli를 통해 ESRI:102012에서 EPSG:4326으로 reproject 해줬다. 이 부분은 추후에 코드로 다루는 방법을 보완할 예정

• nodata설정을 해줘야 나중에 시각화 할 때 nodata 부분은 투명하게 볼 수 있다.

**gdalwarp -t_srs EPSG:4326 -dstnodata 255 ${src_path} ${dst_path}**

결론

• 이 문서는 GDAL 3.x 버전을 기준으로 작성되었다. 2.4 버전에서는 lcc좌표계(ESRI:102012)가 인식이 안되더라
• GDAL vrt로 reproject하는 cli를 코드로 대체할 수 있을 것 같다. GDAL vrt에 대한 학습이 필요

전제 및 가정

지오레퍼런싱(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 {파일 이름}