QGIS Raster to Vector use Python console step2

資料處理步驟：
1.讀取 Raster 並計算 25 等級等值線級距
2.用 GDAL 產生等值線（Contour Lines）
3.Polygonize：Raster 轉等值面多邊形
4.將多邊形圖層載入 QGIS
5.為每個 polygon 計算其 VALUE 所屬級距（LEVEL、MIN、MAX）
6.套用 25 階色彩漸層分類
7.加上級距標籤並加入 QGIS 專案

from osgeo import gdal, ogr, osr

import numpy as np

from PyQt5.QtGui import QColor

from qgis.core import (

    QgsVectorLayer, QgsProject, QgsGraduatedSymbolRenderer, QgsRendererRange,

    QgsSymbol, QgsPalLayerSettings, QgsVectorLayerSimpleLabeling, QgsField

)

from qgis.PyQt.QtCore import QVariant

# ----------------------------

# 0. Path

# ----------------------------

input_raster = r"D:\download\value.tif"

contour_line = r"D:\download\contour25.shp"

contour_poly = r"D:\download\contour25_polygon.shp"

# ----------------------------

# 1. Read raster + get 25-class levels

# ----------------------------

ds = gdal.Open(input_raster)

band = ds.GetRasterBand(1)

arr = band.ReadAsArray()

nodata = band.GetNoDataValue()

valid = arr[arr != nodata]

min_v = float(valid.min())

max_v = float(valid.max())

# 26 contour levels = 25 classes

levels_list = np.linspace(min_v, max_v, 26).tolist()

print("等值線級距：")

for i in range(25):

    print(f"{i+1}: {levels_list[i]:.2f} ~ {levels_list[i+1]:.2f}")

# ----------------------------

# 2. Create contour lines (using raster band, correct API)

# ----------------------------

drv = ogr.GetDriverByName("ESRI Shapefile")

try: drv.DeleteDataSource(contour_line)

except: pass

out_ds = drv.CreateDataSource(contour_line)

srs = osr.SpatialReference()

srs.ImportFromWkt(ds.GetProjection())

line_lyr = out_ds.CreateLayer("contour", srs, ogr.wkbLineString)

line_lyr.CreateField(ogr.FieldDefn("ELEV", ogr.OFTReal))

print("⚙ 產生等值線…")

gdal.ContourGenerate(

    band,                 # <-- 正確：使用 raster band

    ds.GetGeoTransform()[1],

    abs(ds.GetGeoTransform()[5]),

    levels_list,

    0,

    nodata,

    line_lyr,

    0,

    1

)

out_ds = None

print("✔ 等值線完成：", contour_line)

# ----------------------------

# 3. Polygonize (convert raster → polygon)

# ----------------------------

try: drv.DeleteDataSource(contour_poly)

except: pass

poly_ds = drv.CreateDataSource(contour_poly)

poly_lyr = poly_ds.CreateLayer("poly", srs, ogr.wkbPolygon)

poly_lyr.CreateField(ogr.FieldDefn("VALUE", ogr.OFTReal))

print("⚙ Polygonize…")

gdal.Polygonize(band, None, poly_lyr, 0, [], callback=None)

poly_ds = None

print("✔ 等值面完成：", contour_poly)

# ----------------------------

# 4. Load polygon layer to QGIS

# ----------------------------

vl = QgsVectorLayer(contour_poly, "Geothermal_25Levels", "ogr")

prov = vl.dataProvider()

prov.addAttributes([

    QgsField("MIN", QVariant.Double),

    QgsField("MAX", QVariant.Double),

    QgsField("LEVEL", QVariant.Int)

])

vl.updateFields()

# ----------------------------

# 5. Compute class (LEVEL, MIN, MAX)

# ----------------------------

i_min = vl.fields().indexOf("MIN")

i_max = vl.fields().indexOf("MAX")

i_lev = vl.fields().indexOf("LEVEL")

i_val = vl.fields().indexOf("VALUE")

with edit(vl):

    for f in vl.getFeatures():

        v = f["VALUE"]

        if v == nodata:

            f[i_lev] = -1

            f[i_min] = -9999

            f[i_max] = -9999

        else:

            for idx in range(25):

                if levels_list[idx] <= v <= levels_list[idx+1]:

                    f[i_lev] = idx + 1

                    f[i_min] = levels_list[idx]

                    f[i_max] = levels_list[idx+1]

                    break

        vl.updateFeature(f)

print("✔ LEVEL/MIN/MAX 填寫完成")

# ----------------------------

# 6. Style → 25-color gradient

# ----------------------------

ranges = []

for idx in range(25):

    low = levels_list[idx]

    high = levels_list[idx+1]

    label = f"{low:.1f} - {high:.1f}"

    sym = QgsSymbol.defaultSymbol(vl.geometryType())

    sym.setColor(QColor.fromHsv(int(idx * 360/25), 255, 255))

    rng = QgsRendererRange(low, high, sym, label)

    ranges.append(rng)

renderer = QgsGraduatedSymbolRenderer("VALUE", ranges)

renderer.setMode(QgsGraduatedSymbolRenderer.Custom)

vl.setRenderer(renderer)

# ----------------------------

# 7. Labels

# ----------------------------

label = QgsPalLayerSettings()

label.fieldName = "LEVEL"

label.enabled = True

vl.setLabelsEnabled(True)

vl.setLabeling(QgsVectorLayerSimpleLabeling(label))

QgsProject.instance().addMapLayer(vl)

print("🎉 全部完成：25 等級等值面 + 標籤 + 上色 已加入 QGIS")