[GIS] gcp로 georeference 하기, lcc 좌표계에서 EPSG:4326으로 reproject 하기

목적

배경

• 공간 정보 데이터는 맞는데, 그렇다고 이를 바로 지도에 표출하려니 데이터가 다른 분야에서 다루는 데이터라서 공간 정보 데이터의 형식에 맞춰줘야 할때
• 기상 데이터들은 공간에서 일어나는 현상을 다룬 데이터라 위치 정보들이 있지만, 우리가 생각하는 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에 대한 학습이 필요