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")

QGIS Raster to Vector use Python console step1

處理步驟
1.從顏色分類的 RGB 影像中取得所有像素

2. 過濾白色背景區域

3. 使用 K-Means RGB Centroid 進行色彩匹配

4. 將像素映射為指定的 數值分類值

5. 產生可用於 QGIS、ArcGIS、Python 的 GeoTIFF 數值影像

6. 完整設定 NoData、座標系、像元大小

 import numpy as np

from osgeo import gdal, gdalconst

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

# 1. 路徑設定

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

input_path = r"D:\download\87ad4ea7-158c-4f27-9911-57fa9c1ebb472_modified.tif"

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

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

# 2. 讀取 RGB 影像

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

ds = gdal.Open(input_path, gdalconst.GA_ReadOnly)

if ds is None:

    raise Exception("❌ 無法開啟輸入 TIFF")

R = ds.GetRasterBand(1).ReadAsArray().astype(np.float32)

G = ds.GetRasterBand(2).ReadAsArray().astype(np.float32)

B = ds.GetRasterBand(3).ReadAsArray().astype(np.float32)

rows, cols = R.shape

pixels = np.stack([R.flatten(), G.flatten(), B.flatten()], axis=1)

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

# 3. 白色背景過濾（非常重要）

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

# 白底 = RGB 皆 > 225

bg_mask = (pixels[:,0] > 225) & (pixels[:,1] > 225) & (pixels[:,2] > 225)

# 非背景像素

non_bg_pixels = pixels[~bg_mask]

print("原始像素：", len(pixels))

print("非背景像素（有效區域）：", len(non_bg_pixels))

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

# 4. 使用 K-Means 得到的 RGB centroid

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

centroids = np.array([

    [117, 249,  20],   # 群 0 → 32.5

    [237,  28,  64],   # 群 1 → 45

    [148, 202, 235],   # 群 2 → 10

    [231, 224,  30],   # 群 3 → 37.5

    [ 19,  24,  20],   # 群 4 → NoData

    [ 56, 245, 161],   # 群 5 → 27.5

    [102,  97, 227],   # 群 6 → 22.5

], dtype=np.float32)

mapping = {

    0: 32.5,

    1: 45,

    2: 10,

    3: 37.5,

    4: -9999,

    5: 27.5,

    6: 22.5

}

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

# 5. 計算 pixel → 最近 centroid

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

dist = ((pixels[:, None, :] - centroids[None, :, :]) ** 2).sum(axis=2)

cluster_idx = np.argmin(dist, axis=1)

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

# 6. 對背景像素直接指定為 -9999

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

value_arr = np.vectorize(mapping.get)(cluster_idx)

value_arr[bg_mask] = -9999

value_arr = value_arr.reshape(rows, cols).astype(np.float32)

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

# 7. 輸出 GeoTIFF

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

driver = gdal.GetDriverByName("GTiff")

out_ds = driver.Create(output_path, cols, rows, 1, gdalconst.GDT_Float32)

out_ds.SetGeoTransform(ds.GetGeoTransform())

out_ds.SetProjection(ds.GetProjection())

band = out_ds.GetRasterBand(1)

band.WriteArray(value_arr)

band.SetNoDataValue(-9999)

band.FlushCache()

del out_ds

print("🎉 已完成：輸出數值影像 →", output_path)