127 lines
4.6 KiB
Python
127 lines
4.6 KiB
Python
#!/usr/bin/env python3
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
样本数据离线验证 / 渲染基准生成器。
|
||
|
||
用途：在未搭好 Qt/VTK 环境前，用 Python(matplotlib) 复现剖面散点(#17)、网格等值面(#18)
|
||
渲染效果，验证「样本解析 + colorBar 色阶映射 + 异常叠加」逻辑正确；并量化两条剖面的几何
|
||
关系，为 dd_voxel 可信度(设计 §10/§14 K-10)提供依据。
|
||
|
||
产出的 PNG 作为后续 VTK 渲染的「地面真值」对照图。
|
||
|
||
运行：
|
||
python tools/validate_samples.py
|
||
输出：docs/_validate/ref_17_scatter.png, ref_18_grid.png（该目录已 .gitignore）
|
||
|
||
依赖：numpy, matplotlib
|
||
"""
|
||
import json
|
||
import os
|
||
import re
|
||
import sys
|
||
|
||
import numpy as np
|
||
import matplotlib
|
||
matplotlib.use("Agg")
|
||
import matplotlib.pyplot as plt
|
||
from matplotlib.colors import ListedColormap, BoundaryNorm
|
||
|
||
DATA = os.path.join(os.path.dirname(__file__), "..", "docs", "剖面网格数据的色阶数据2等文件")
|
||
OUT = os.path.join(os.path.dirname(__file__), "..", "docs", "_validate")
|
||
|
||
|
||
def load(name):
|
||
with open(os.path.join(DATA, name), encoding="utf-8") as f:
|
||
return json.load(f)["data"]
|
||
|
||
|
||
def parse_color(c):
|
||
"""支持 #RRGGBB 与 rgba(r,g,b,a)；返回 0-1 RGB。"""
|
||
c = c.strip()
|
||
if c.startswith("#"):
|
||
return (int(c[1:3], 16) / 255, int(c[3:5], 16) / 255, int(c[5:7], 16) / 255)
|
||
nums = re.findall(r"[\d.]+", c)
|
||
return (float(nums[0]) / 255, float(nums[1]) / 255, float(nums[2]) / 255)
|
||
|
||
|
||
def colorbar(cb):
|
||
"""colorBar: [[value, color], ...] -> (sorted values, colors)。"""
|
||
pairs = sorted(((float(v), parse_color(c)) for v, c in cb), key=lambda p: p[0])
|
||
return np.array([p[0] for p in pairs]), [p[1] for p in pairs]
|
||
|
||
|
||
def render_scatter():
|
||
raw = load("剖面原数据1.txt")
|
||
vals, cols = colorbar(load("剖面原数据的色阶数据1.txt")["properties"]["colorBar"])
|
||
x, y, v = np.array(raw["xlist"]), np.array(raw["ylist"]), np.array(raw["vlist"])
|
||
cmap = ListedColormap(cols)
|
||
norm = BoundaryNorm(np.append(vals, vals[-1] * 1.2), cmap.N)
|
||
plt.figure(figsize=(11, 3.2))
|
||
plt.scatter(x, y, c=v, cmap=cmap, norm=norm, s=6, marker="s")
|
||
plt.gca().invert_yaxis()
|
||
plt.title("scatter (raw profile) ~#17")
|
||
plt.colorbar(shrink=0.7)
|
||
plt.tight_layout()
|
||
plt.savefig(os.path.join(OUT, "ref_17_scatter.png"), dpi=90)
|
||
plt.close()
|
||
print(f"[#17] pts={len(x)} x[{x.min():.1f},{x.max():.1f}] "
|
||
f"y[{y.min():.1f},{y.max():.1f}] v[{v.min():.1f},{v.max():.1f}]")
|
||
|
||
|
||
def render_grid():
|
||
g = load("剖面网格数据1.txt")
|
||
vals, cols = colorbar(load("剖面网格数据的色阶数据1.txt")["properties"]["colorBar"])
|
||
gx, gy, gv = np.array(g["x"]), np.array(g["y"]), np.array(g["v"]) # gv: [22][100]
|
||
X, Y = np.meshgrid(gx, gy)
|
||
plt.figure(figsize=(11, 3.0))
|
||
cf = plt.contourf(X, Y, gv, levels=vals, colors=cols[:len(vals) - 1])
|
||
plt.contour(X, Y, gv, levels=vals, colors="k", linewidths=0.4)
|
||
for a in load("剖面网格数据1——对应的异常圈定数据.txt"):
|
||
coord = a.get("location", {}).get("coordinate", [])
|
||
if coord and isinstance(coord[0], dict):
|
||
plt.plot([p["x"] for p in coord], [p["y"] for p in coord], "k--", lw=1.5)
|
||
plt.gca().invert_yaxis()
|
||
plt.title("grid bands + contour + anomaly ~#18")
|
||
plt.colorbar(cf, shrink=0.7)
|
||
plt.tight_layout()
|
||
plt.savefig(os.path.join(OUT, "ref_18_grid.png"), dpi=90)
|
||
plt.close()
|
||
print(f"[#18] grid V{gv.shape} vmin/max {g['vmin']:.1f}/{g['vmax']:.1f}")
|
||
|
||
|
||
def analyze_voxel_geometry():
|
||
"""量化两条剖面几何关系：判断 dd_voxel 三维插值的数据支撑充分性。"""
|
||
def line(name):
|
||
r = load(name)
|
||
return np.column_stack([r["projectXList"], r["projectYList"]])
|
||
|
||
def principal_dir(p):
|
||
c = p - p.mean(0)
|
||
_, _, vt = np.linalg.svd(c, full_matrices=False)
|
||
return vt[0], p.mean(0)
|
||
|
||
p1, p2 = line("剖面原数据1.txt"), line("剖面原数据2.txt")
|
||
d1, c1 = principal_dir(p1)
|
||
d2, c2 = principal_dir(p2)
|
||
ang = np.degrees(np.arccos(abs(np.clip(d1 @ d2, -1, 1))))
|
||
n1 = np.array([-d1[1], d1[0]])
|
||
perp_gap = abs((c2 - c1) @ n1)
|
||
print(f"[voxel] line1 pts={len(p1)} line2 pts={len(p2)} | "
|
||
f"angle={ang:.1f}deg perp_gap={perp_gap:.1f}m "
|
||
f"=> {'crossing' if ang > 30 else 'sub-parallel'}; "
|
||
f"volume between lines is interpolated/extrapolated "
|
||
f"(credible voxel needs >=3 lines or a 3D grid)")
|
||
|
||
|
||
def main():
|
||
os.makedirs(OUT, exist_ok=True)
|
||
render_scatter()
|
||
render_grid()
|
||
analyze_voxel_geometry()
|
||
print("written:", sorted(os.listdir(OUT)))
|
||
|
||
|
||
if __name__ == "__main__":
|
||
sys.stdout.reconfigure(encoding="utf-8")
|
||
main()
|