REMOTE SENSING

QGIS 파이썬 콘솔: UAV 기반 RGB 영상의 식생지수 계산

유병혁 2021. 5. 7. 00:43

이번 글에서는 QGIS 파이썬 콘솔에서 UAV 기반 RGB 영상의 식생지수를 계산해 보겠습니다.

 

현재 널리 쓰이고 있는 저가 UAV는 거의 대부분 일반적인 RGB 카메라를 탑재하고 있는데요, 전세계 70%를 점유하고 있는 중국 DJI 드론도 대부분 sRGB 스타일의 카메라를 탑재하고 있습니다. sRGB는 1996년에 HP와 Microsoft과 협력하여 만든 표준 RGB 색 공간(standard RGB color space)입니다.

 

보통 다중분광 위성영상의 경우 RGB 밴드와 NIR 밴드를 함께 취득하여 NDVI와 같은 식생지수를 계산했었는데요, UAV 기반 RGB영상은 공간해상도는 매우 우수한 반면, 분광해상도는 원격탐사 목적에 제한적이라고 볼 수 있습니다.

 

최근에 UAV에 탑재 가능한 적외선 카메라들 상당수는 기존 카메라를 NIR 필터로 개조한 형태입니다.

사진 출처: http://diydrones.com/profiles/blogs/diy-camera-filter-swap-on-a-canon-powershot

 

이런 경우에는 natural-color composite와 false-color composite 영상을 동시에 취득할 수가 없는데요, 여기서는 RGB 영상만으로 식생지수를 계산하고, 이것을 QGIS에서 처리하는 방법을 정리해 보겠습니다. 먼저, RGB 영상의 식생지수에 대해 알아보겠습니다. 다양한 식생지수가 있지만, 가장 단순한 지수를 계산해 보겠습니다.

 

Woebbecke et al. (1995)은 식생과 비식생 영역을 구분할 수 있는 강도(intensity) 영상인 '초과녹색식생지수(excess green vegetation index, 줄여서 ExG)'를 제안하였습니다.

 

Journal Article:

Woebbecke, D.M., Meyer, G.E., B.K., V., Mortensen, D.A., 1995. Color indices for weed

identification under various soil, residue and lighting conditions. Transactions of the ASAE 38, 259–269.

 

ExG는 아래와 같이 정의되는데요,

 

 

여기서 r, g, b는,

 

 

 

여기서 R*, G*, B*는 아래와 같이 정의됩니다.

 

 

Rm, Gm, Bm은 각 원색(primary color)에 대한 최대 색조값(maximum tonal value)을 나타내는데요, 따라서 R*, G*, B*는 기존 값 범위가 0~1 범위로 일반화된 RGB 값으로 정의됩니다. UAV 기반 RGB 24bit 영상에서는 Rm = Gm = Bm = 255이기 때문에, 위 두 식은 아래와 같이 간소화할 수 있습니다.

 

 

 

따라서,

 

 

r, g, b는 각 원색의 색 좌표(chromatic coordinate)로 정의될 수 있는데요, 아래 그림은 r좌표와 g좌표로 정의된 색도분포표(chromaticity diagram)입니다. 1931년 국제조명위원회(CIE)가 제정한 것으로, x축은 적색 좌표, y축은 녹색 좌표입니다.

 

그림 출처:

https://ecse.rpi.edu/~schubert/Light-Emitting-Diodes-dot-org/chap17/F17-03%20Chromaticity%20diagram%20(Gage).jpg

 

자, 그럼 이제 UAV 영상에서 초과녹색생지수 ExG를 계산해 볼까요?!

위와 같이 R, G, B값을 r, g, b 색 좌표로 변환한 후 2g - r - b를 적용하면 되겠는데요,

 

가장 쉬운 방법은 QGIS 상단 메뉴 '래스터 > 래스터 계산기'를 이용하는 것인데, 저는 파이썬 코드로 처리하겠습니다.

QGIS에서 상단 메뉴 '플러그인 > 파이썬 콘솔'을 실행하시면, 파이썬 명령어를 실행하실 수 있는 쉘이 보입니다.

 

아래 코드는 일단 캔버스로 래스터 레이어를 추가하고 각 R, G, B 밴드를 변수로 지정하는 과정입니다.

from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
# 래스터 레이어 추가
rasterpath = "D:/GEODATA/odm_orthophoto.tif"
bohLayer = QgsRasterLayer(rasterpath, "ODM orthophoto")
# 캔버스로 래스터 레이어 추가
QgsProject.instance().addMapLayer(bohLayer)
entries = []
#1번 밴드 정의(Red)
boh1 = QgsRasterCalculatorEntry()
boh1.ref = 'boh@1'
boh1.raster = bohLayer
boh1.bandNumber = 1
entries.append(boh1)
#2번 밴드 정의(Green)
boh2 = QgsRasterCalculatorEntry()
boh2.ref = 'boh@2'
boh2.raster = bohLayer
boh2.bandNumber = 2
entries.append(boh2)
#3번 밴드 정의(Blue)
boh3 = QgsRasterCalculatorEntry()
boh3.ref = 'boh@3'
boh3.raster = bohLayer
boh3.bandNumber = 3
entries.append(boh3)

위의 코드를 QGIS 파이썬 콘솔에서 실행한 결과 화면은 아래와 같습니다.

속리산국립공원 문장대 훼손지 일원을 DJI 인스파이어 1 젠뮤즈 X3로 촬영한 영상입니다.

아래는 앞서 소개드린 초과녹색식생지수(ExG)를 계산하기 위한 코드입니다.

#ExG 계산
outputLayer = "D:/GEODATA/odm_orthophoto_exg.tif"
calc = QgsRasterCalculator('2 * (boh@2/(boh@1+boh@2+boh@3))-(boh@1/(boh@1+boh@2+boh@3)) - (boh@3/(boh@1+boh@2+boh@3))', outputLayer, 'GTiff', bohLayer.extent(), bohLayer.width(), bohLayer.height(), entries)
calc.processCalculation()
addLayer = QgsRasterLayer(outputLayer, "ExG")
QgsProject.instance().addMapLayer(addLayer)

 

위 코드를 실행한 ExG 계산 결과는 아래 화면과 같은데요,

 

식생 영역과 비식생 영역(노출된 암석, 훼손지, 인공물)이 대비되고 있음을 볼 수 있습니다.

RGB 식생지수는 이외에도 많은 연구가 있는데요, 같은 방식으로 공식을 적용해보면 되겠습니다.