PyQGIS: KOMPSAT-3 위성영상 분광지수 이미지 생성 코드 예제

안녕하세요? 이번 글은 PyQGIS를 통해 KOMPSAT-3 위성영상의 분광지수 이미지를 생성하는 코드 예제를 제공해 봅니다. QGIS 메뉴 바에서 "플러그인 > 파이썬 콘솔"에서 실행해보시면 되겠습니다.

 

먼저, 필요한 라이브러리를 호출하겠습니다.

import os
from qgis.core import (
    QgsRasterLayer, 
    QgsProject,
    QgsSingleBandPseudoColorRenderer,
    QgsStyle
)
from qgis.analysis import (
    QgsRasterCalculator,
    QgsRasterCalculatorEntry,
)
from osgeo import gdal

 

KOMPSAT-3 위성영상의 밴드별 경로, 생성할 분광지수 이미지의 경로를 설정합니다. 여기서는 NDVI(정규식생지수)와 NDWI(정규수분지수)를 계산해 보겠습니다.

# 경로 설정
folder_name = "K3_20240114043957_62211_09341264_L1G"
raster_dir = "D:GEODATA/{folder_name}")

b_band_path = os.path.join(raster_dir, f"{folder_name}_B.tif")
g_band_path = os.path.join(raster_dir, f"{folder_name}_G.tif")
r_band_path = os.path.join(raster_dir, f"{folder_name}_R.tif")
n_band_path = os.path.join(raster_dir, f"{folder_name}_N.tif")

output_dir = "D:GEODATA/Output"
ndvi_output_path = os.path.join(output_dir, "NDVI.tif")
ndwi_output_path = os.path.join(output_dir, "NDWI.tif")

 

분광지수 계산을 위한 함수를 정의합니다.

def calculate_index(expression, output_path, layer1, layer2):
    # 분광지수 계산
    entries = []

    layer1_entry = QgsRasterCalculatorEntry()
    layer1_entry.ref = f'{layer1.name()}@1'
    layer1_entry.raster = layer1
    layer1_entry.bandNumber = 1
    entries.append(layer1_entry)

    layer2_entry = QgsRasterCalculatorEntry()
    layer2_entry.ref = f'{layer2.name()}@1'
    layer2_entry.raster = layer2
    layer2_entry.bandNumber = 1
    entries.append(layer2_entry)

    calculator = QgsRasterCalculator(
        expression,
        output_path,
        "GTiff",
        layer1.extent(),
        layer1.width(),
        layer1.height(),
        entries
    )
    calculator.processCalculation()

 

NDVI와 NDWI를 계산합니다.

# NDVI 계산
nir_layer = QgsRasterLayer(n_band_path, "NIR")
red_layer = QgsRasterLayer(r_band_path, "Red")
calculate_index('(NIR@1 - Red@1) / (NIR@1 + Red@1)', ndvi_output_path, nir_layer, red_layer)

# NDWI 계산
green_layer = QgsRasterLayer(g_band_path, "Green")
calculate_index('(Green@1 - NIR@1) / (Green@1 + NIR@1)', ndwi_output_path, green_layer, nir_layer)

 

NDVI와 NDWI 이미지에 색상표를 적용해서 시각화하기 위한 함수를 정의해 줍니다. 해당 함수의 출처는 아래와 같습니다. 기존 코드에서 이산을 선형으로 변경하여 작성했습니다.

Setting the visual style of a raster from python console in QGIS

def render_raster(layer, band, spectrum):
    prov = layer.dataProvider()
    src_ds = gdal.Open(layer.source())
    src_band = src_ds.GetRasterBand(band)

    if src_band.GetMinimum() is None or src_band.GetMaximum() is None:
        src_band.ComputeStatistics(0)
    band_min = src_band.GetMinimum()
    band_max = src_band.GetMaximum()

    fcn = QgsColorRampShader()
    fcn.setColorRampType(QgsColorRampShader.Interpolated)

    item_list = [
        QgsColorRampShader.ColorRampItem(
            band_min + (n / (len(spectrum) - 1)) * (band_max - band_min),
            QColor(color), lbl=f"{band_min + (n / (len(spectrum) - 1)) * (band_max - band_min):.2f}"
        )
        for n, color in enumerate(spectrum)
    ]

    fcn.setColorRampItemList(item_list)
    shader = QgsRasterShader()
    shader.setRasterShaderFunction(fcn)

    renderer = QgsSingleBandPseudoColorRenderer(prov, band, shader)
    renderer.setClassificationMin(band_min)
    renderer.setClassificationMax(band_max)

    layer.setRenderer(renderer)
    layer.triggerRepaint()

 

이어서 레이어를 추가하고 색상표를 업데이트하는 함수를 정의합니다.

def add_layer_with_rendering(file_path, name, spectrum):
    layer = QgsRasterLayer(file_path, name)
    if layer.isValid():
        QgsProject.instance().addMapLayer(layer)
        render_raster(layer, 1, spectrum)
    else:
        print(f"{name} layer is not valid!")

 

이제 해당 함수들을 활용해서 NDVI와 NDWI 이미지를 레이어 추가 및 렌더링 합니다.

# NDVI 레이어 추가 및 렌더링
ndvi_color_list = [
    '#a50026', '#d73027', '#f46d43', '#fdae61', '#fee08b', 
    '#d9ef8b', '#a6d96a', '#66bd63', '#1a9850', '#006837'
]
add_layer_with_rendering(ndvi_output_path, "NDVI", ndvi_color_list)

# NDWI 레이어 추가 및 렌더링
ndwi_color_list = [
    '#f7fbff', '#deebf7', '#c6dbef', '#9ecae1', '#6baed6', 
    '#4292c6', '#2171b5', '#08519c', '#08306b'
]
add_layer_with_rendering(ndwi_output_path, "NDWI", ndwi_color_list)

print("NDVI and NDWI layers successfully added to the QGIS project.")

 

결과는 다음과 같습니다.