[GeoTIFF] GeoTIFF 파일의 corner coordinates를 구하자

개요

오늘도 그렇듯 python이랑 gdal을 쓰면서 일을 하고 있었습니다.
제가 구현하려는 기능 중에 '생성된 geotiff 파일의 extent가 좌표계의 유효 범위를 넘어가면 geotiff를 유효 범위에 맞춰 자르기' 가 있었는데요,
평소와 같이 아~ GeoTransform을 이용해서 계산하면 되겠지? 하고 안일하게 생각하다가 생각과 다른 결과물을 얻었습니다.
왜일까요?

오늘 이렇게 쉬워보이는 주제를 들고 왔는데 무심코 쓰다 보면 저처럼 당황하실 수 있습니다.
오늘 다룰 주제는 ApplyGeoTransform입니다.

개념

gdal에서 다루는 GeoTransform의 개념은 다음과 같습니다. gdal geotransform 문서
여기를 보면 이렇게 나옵니다.

GT(0) x-coordinate of the upper-left corner of the upper-left pixel.
GT(1) w-e pixel resolution / pixel width.
GT(2) row rotation (typically zero).
GT(3) y-coordinate of the upper-left corner of the upper-left pixel.
GT(4) column rotation (typically zero).
GT(5) n-s pixel resolution / pixel height (negative value for a north-up image).

그리고 이를 이용해서 corner coordinates를 이렇게 계산하면 된다고 말합니다.

X_geo = GT(0) + X_pixel * GT(1) + Y_line * GT(2)
Y_geo = GT(3) + X_pixel * GT(4) + Y_line * GT(5)

여기까지는 늘 보던 개념과 식입니다.
그런데 말입니다, 단위에 매우 조심해야 합니다.
EPSG:4326은 그 단위가 degree로, 작은 영역의 데이터는 pixel resolution이 소수점 몇몇자리 부터 시작합니다.
이걸 깜빡하고 위 식을 그대로 적용하면 제대로 된 corner coordinates를 얻을 수 없습니다.

해결법은 의외로 단순했습니다.

ApplyGeoTransform

위 함수는 다음과 같이 동작합니다.

ds = gdal.Open(tiff_path)
gt = ds.GetGeoTransform()
xs = ds.RasterXSize
ys = ds.RasterYSize

'''
args : GeoTransform, Pixel_coordinates_x, Pixel_coordinates_y
'''
ulx, uly = gdal.ApplyGeoTransform(gt, 0, 0)
lrx, lry = gdal.ApplyGeoTransform(gt, xs, ys)

힌트는 스텍오버플로우의 How I can read corner coordinates in degrees에서 ApplyGeoTransform 함수를 보고 얻었습니다.

ApplyGeoTransform은 pixel 좌표를 world 좌표계로 바꿔주는 용도로 많이 씁니다.
원래 용도는 GeoTransform이 Affine Transform 계산에 필요한 정보를 담고 있어 여러 변환에 해당 함수를 사용할 수 있습니다.
이 함수와 같이 쓰이는 유용한 함수가 또 하나 있습니다.

InvGeoTransform

Affine Transform이 공간 A에서 B로 행해진다고 할 때, InvGeoTransform은 B에서 A로 변환을 시행할 수 있게 GeoTransform 값을 역으로 바꿔줍니다. 이 함수를 본 지 오랫만이기도 해서 구글로 검색해 봤는데요, 굉장히 좋은 참고 자료를 찾았습니다.The Affine Transform 베를린 홈볼트 대학교의 Dr.Dirk Pflugmacher과 Dr.Matthias Baumann께 감사드립니다. Again, I'm so empressed by your good lecture note. Sincerely appereciate for your works, Dr.Dirk Pflugmacher and Dr.Matthias Baumann.

덕분에 새로운 용례를 알게 되었습니다.
이 두 함수는 이렇게 쓸 수 있습니다. 위 자료의 일부를 발췌했습니다.

gt = ds.GetGeoTransform()
#gt = (-63.86859973338657,0.00026949458523585647,0.0,-24.446662310500248,0.0,-0.00026949458523585647)

invGT = gdal.InvGeoTransform(gt)
#invGT = (236994.0, 3710.6496931091187, 0.0, -90713.0, 0.0, -3710.6496931091187)

# Define some example lat/lon
lat = -25.0
lon = -63.0
# Convert these geo-coordinates into array coordinates
x, y = gdal.ApplyGeoTransform(invGT, lon, lat)
# x, y = (3223.0693341255246, 2053.242327727974)

위 예제에서 볼 수 있듯 역으로 world 좌표계에서 pixel 좌표계로 변환할 때 사용합니다.
참고로 pixel 좌표는 integer여야 하기 때문에 int로 형변환 잊지 마세요!

결론

생각지도 못한 곳에서 정말 오랫만에 ApplyGeoTransform을 만났습니다. 그리고 좋은 학습 자료도 얻었구요.

요즘 챗지피티 유료로 결제해서 비서처럼 쓰고 있는데 GIS 쪽이 희귀해서 그런지 틀리거나 검증되지 않은 코드들을 답변으로 내어놓습니다.
이만큼 여기가 마이너해서 그런거 같은데 참 묘합니다.
gdal의 corner coordinates는 upper left와 lower right로 이루어져 있는데 lower left와 upper right로 주석을 달아놓고,
계산식에서 고려해야 할 사항들을 빼고 계산하는 등의 오류가 있었습니다.
당분간 챗지피티는 초기 진입 장벽을 빠르게 뚫고 나가는 데에 쓰거나 정말 손이 부족할 때 보고서 요약, 생성 등으로 쓰기로 했습니다.