PyQGIS: 출현/배경 GPKG 파일 생성 코드 예제

안녕하세요? 이번 글은 출현 데이터와 배경 데이터를 병합하여 하나의 GPKG 파일을 만들고 이를 QGIS 레이어로 추가하는 코드 예제를 소개해 보겠습니다. 해당 코드는 Python용 종 분포 모델링(SDM: Species Distribution Modeling) 도구, elapid(엘라피드)의 일부 소스 코드를 QGIS Python 콘솔에서 동작하도록 수정한 것입니다.

elapid(엘라피드): Python용 종 분포 모델링 도구 소개

 

먼저 필요한 라이브러리를 호출합니다.

import os
import numpy as np
import pandas as pd
import geopandas as gpd

from osgeo import gdal
from shapely.geometry import Point

출현 데이터

출현 데이터를 읽어오겠습니다.

# 출현 데이터 열기
presence = gpd.read_file(gpkg_path)
presence.head()

배경 데이터

래스터 범위에서 무작위로 지리적 샘플 포인트를 생성하는 함수를 정의한 후, 이를 통해 배경 데이터를 생성해 보겠습니다.

def sample_raster(raster_path, count, nodata=None):
    """래스터(extent)를 기반으로 무작위 지리적 샘플 포인트 생성

    Args:
        raster_path (str): 샘플링할 래스터 파일 경로
        count (int): 생성할 샘플 수
        nodata (float): 이 값이 포함된 픽셀을 nodata 마스크에 추가

    Returns:
        샘플링된 포인트 좌표를 포함한 GeoSeries
    """

    # 파일 존재 여부 확인
    if not os.path.exists(raster_path):
        raise FileNotFoundError(f"파일을 찾을 수 없습니다: {raster_path}")

    # GDAL을 사용하여 래스터 파일 열기
    raster = gdal.Open(raster_path)
    if raster is None:
        raise ValueError(f"래스터 파일을 열 수 없습니다: {raster_path}")

    band = raster.GetRasterBand(1)  # 첫 번째 밴드 가져오기
    transform = raster.GetGeoTransform()  # 지리적 변환 정보 가져오기
    cols, rows = raster.RasterXSize, raster.RasterYSize  # 래스터의 크기

    # NoData 마스크 설정
    if nodata is not None:
        nodata_value = band.GetNoDataValue() or nodata
        mask = band.ReadAsArray() != nodata_value
    else:
        mask = np.ones((rows, cols), dtype=bool)  # 전체를 유효한 영역으로 설정

    # 유효한 픽셀의 인덱스 가져오기
    valid_indices = np.column_stack(np.where(mask))

    if len(valid_indices) < count:
        raise ValueError("유효한 픽셀 수가 요청된 샘플 수보다 적습니다.")

    # 무작위 샘플링
    sampled_indices = valid_indices[
        np.random.choice(valid_indices.shape[0], count, replace=False)
    ]

    # 픽셀 인덱스를 지리 좌표로 변환
    sampled_points = [
        Point(transform[0] + col * transform[1] + row * transform[2],
              transform[3] + col * transform[4] + row * transform[5])
        for row, col in sampled_indices
    ]

    # 좌표계 설정
    spatial_ref = raster.GetSpatialRef()
    crs = spatial_ref.GetAuthorityCode(None) if spatial_ref else None
    points = gpd.GeoSeries(sampled_points, crs=crs)

    return gpd.GeoDataFrame(geometry=points)

# 래스터에서 배경 데이터 생성
background = sample_raster(raster_path, count=10_000, nodata=0)
background.head()

출현-배경 레이어 병합 후 추가

이제 출현 데이터와 배경 데이터를 병합한 후, QGIS에 스타일 적용 후 추가하는 함수를 정의합니다. QGIS 파이썬 콘솔에서 실행한 결과는 다음과 같습니다.

def merge_presence_background(presence, background, output_path, opacity=0.5, size=1.0):
    """Presence와 Absence 데이터를 병합하여 QGIS에 스타일 적용 후 추가

    Args:
        presence (gpd.GeoDataFrame): 존재(presence) 데이터
        background (gpd.GeoDataFrame): 배경(background) 데이터
        output_path (str): 저장할 GeoPackage 파일 경로
        opacity (float): 투명도 (기본값: 0.5)
        size (float): 점 크기 (기본값: 1.0)
    """
    # presence의 CRS로 통일
    background = background.to_crs(presence.crs)
    
    # Presence (1) 및 Background (0) 클래스 추가
    presence["class"] = 1
    background["class"] = 0

    # 병합하여 하나의 GeoDataFrame 생성
    combined_df = gpd.GeoDataFrame(pd.concat([presence, background], ignore_index=True))

    # GeoPackage로 저장
    os.makedirs(os.path.dirname(output_path), exist_ok=True)
    layer_name = "Presence-Absence Layer"
    combined_df.to_file(output_path, layer=layer_name, driver="GPKG")

    print(f"파일 저장 완료: {output_path} (레이어: {layer_name})")

    # QGIS에 레이어 추가 및 스타일 설정
    vl = QgsVectorLayer(output_path, layer_name, "ogr")
    if not vl.isValid():
        print(f"레이어를 불러올 수 없습니다: {output_path}")
        return

    # 클래스별 스타일 설정
    categories = []

    # Presence (blue)
    symbol = QgsSymbol.defaultSymbol(vl.geometryType())
    symbol.setColor(QColor("blue"))
    symbol.setOpacity(opacity)
    symbol.setSize(size)
    categories.append(QgsRendererCategory(1, symbol, "Presence"))

    # Absence (red)
    symbol = QgsSymbol.defaultSymbol(vl.geometryType())
    symbol.setColor(QColor("red"))
    symbol.setOpacity(opacity)
    symbol.setSize(size)
    categories.append(QgsRendererCategory(0, symbol, "Absence"))

    # 카테고리 렌더러 적용
    renderer = QgsCategorizedSymbolRenderer("class", categories)
    vl.setRenderer(renderer)

    # QGIS 프로젝트에 레이어 추가
    QgsProject.instance().addMapLayer(vl)
    print(f"QGIS에 레이어 추가 및 스타일 적용 완료: {layer_name}")

# Presence-Absence 데이터 병합 후 추가
merge_presence_background(presence, background, output_path)