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