guanlan-src/scripts/sim_hub/wake_real.py

"""真实山地风场尾流图 (已脱敏)。
布局 = 我们做过的一个真实山地风场 90 台机位; 风向分布 = 该场 1,218,810 条实测 10min 记录。
脱敏: 经纬度转本地米制后做刚性变换(旋转+平移) —— 相对几何不变(尾流物理为真), 反查不出位置;
      台号重编, 不出现场名/业主/厂商。
边界: 平坦地形假设(未建地形), 推力曲线为通用示意值 → 只看相对形态与排序, 不报绝对损失。
引擎: PyWake (DTU, MIT) · Bastankhah-Porte-Agel 2014。"""
import os, math, numpy as np, matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib import font_manager as fm
from py_wake.site import UniformSite
from py_wake.wind_turbines import WindTurbine
from py_wake.wind_turbines.power_ct_functions import PowerCtTabular
from py_wake.literature.gaussian_models import Bastankhah_PorteAgel_2014

for cand in ["PingFang SC", "Hiragino Sans GB", "Songti SC", "STHeiti"]:
    try:
        fm.findfont(fm.FontProperties(family=cand), fallback_to_default=False)
        matplotlib.rcParams["font.sans-serif"] = [cand, "DejaVu Sans"]; break
    except Exception: continue
matplotlib.rcParams["axes.unicode_minus"] = False
# ---- 浅蓝主题 (与观澜仿真台页面同一套 token) ----
BG="#edf7f9"; PANEL="#ffffff"; INK="#133047"; INK2="#28536b"; MUT="#61788b"
LINE="#cbe2eb"; ACC="#0e8898"; MEAS="#2077a8"; AMB="#d48818"; BAD="#d94b43"; GRN="#23875b"
from matplotlib.colors import LinearSegmentedColormap
# 自由来流 = 近乎页面底色, 尾流越深越蓝
WAKE_CMAP = LinearSegmentedColormap.from_list("guanlan_wake",
    ["#0d4a63", "#2077a8", "#6cb6d4", "#b6e0ec", "#e8f6fa"])
DEF_CMAP  = LinearSegmentedColormap.from_list("guanlan_deficit",
    ["#e8f6fa", "#f0d9a8", "#d48818", "#d94b43", "#8c2420"])


MAP="/Users/yuanying/wind-analytics/app_ETL/configs/contracts/mawang_machine_map.csv"
# ★ 数据实况: 本表 90 台里混了两套坐标系 —— 50 台是经纬度(WGS84), 40 台(二期,上海电气)是 CGCS2000 3度带 CM108E 投影米制。
#   直接当经纬度读会把比例尺拉爆 1e11 量级。此处逐行判别后统一到 WGS84。
from pyproj import Transformer
utm2ll = Transformer.from_crs("EPSG:4545", "EPSG:4326", always_xy=True)  # CGCS2000 3度带 CM108E -> lon/lat
lat=[]; lon=[]; dia=[]; n_ll=0; n_utm=0
with open(MAP, encoding="utf-8-sig") as fh:
    import csv as _csv
    for r in _csv.DictReader(fh):
        try: a=float(r["lat"]); o=float(r["lon"])
        except (TypeError, ValueError): continue
        if 73 < o < 136 and 18 < a < 54:
            la_, lo_ = a, o; n_ll += 1
        elif 1e5 < o < 1e6 and 1e6 < a < 6e6:  # 投影米制
            lo_, la_ = utm2ll.transform(o, a); n_utm += 1
        else:
            continue
        m = r.get("model","") or ""
        d = 121. if "121" in m else (116. if "116" in m else 118.)
        lat.append(la_); lon.append(lo_); dia.append(d)
print(f"坐标解析: 经纬度(一期) {n_ll} 台 + 3度带投影(二期) {n_utm} 台 = {len(lat)} 台")
lat=np.array(lat); lon=np.array(lon); dia=np.array(dia)
R=6371000.; lat0=lat.mean(); lon0=lon.mean()
x = R*np.radians(lon-lon0)*math.cos(math.radians(lat0))
y = R*np.radians(lat-lat0)
assert np.ptp(x) < 4e4 and np.ptp(y) < 4e4, f"场跨度异常: {np.ptp(x):.0f} x {np.ptp(y):.0f} m"
print(f"场跨度: {np.ptp(x):.0f} m (东西) x {np.ptp(y):.0f} m (南北)")

# --- 脱敏刚性变换: 坐标与风向同步旋转, 尾流几何完全不变 ---
TH = math.radians(20.0)    # 温和旋转脱敏(保持横构图); 风向必须同步 **减去** 该角
#   已验: wd = wd_true - theta 时逐台亏损与不旋转基准逐位相同 (最大差 0.0000);
#   若写成 + theta, 尾流会朝错误方向拖 (曾错, 已修)
xr =  x*math.cos(TH) - y*math.sin(TH)
yr =  x*math.sin(TH) + y*math.cos(TH)
xr -= xr.min(); yr -= yr.min()
WD_TRUE = 165.0                     # 实测主导风向 (33.5%, n=1,218,810)
WD = (WD_TRUE - math.degrees(TH)) % 360   # 已验: 逐台亏损与不旋转基准逐位相同

D_MEAN = float(np.mean(dia))
def mk_wt(D, H, kw):
    u=np.arange(0,26.); p=np.clip(((u-3)/(11-3))**3,0,1)*kw*1e3; p[u<3]=0; p[u>25]=0
    ct=np.where((u>=3)&(u<=25), np.clip(0.8*np.minimum(1,(11/np.maximum(u,3))**2),.05,.8),0)
    return WindTurbine(name="WT", diameter=D, hub_height=H, powerCtFunction=PowerCtTabular(u,p,'w',ct))

site=UniformSite(p_wd=[1], ti=.11, ws=8.)
wfm=Bastankhah_PorteAgel_2014(site, mk_wt(D_MEAN, 80., 2000.), k=0.0324555)
sim=wfm(xr, yr, wd=WD, ws=8.)
fmap=sim.flow_map(wd=WD, ws=8.)
Z=fmap.WS_eff.squeeze().values; X=fmap.x.values; Y=fmap.y.values
# flow_map 的 WS_eff 已是 (len(Y), len(X)); 不得转置。
# (旧代码写 "if Z.shape[0]==len(X): Z=Z.T" —— 网格是 500x500 方阵时每次误触发, 把尾流画成竖向)
assert Z.shape == (len(Y), len(X)), f"flow_map 方向变了: {Z.shape} vs {(len(Y), len(X))}"
deficit = (8.0 - np.asarray(sim.WS_eff.squeeze().values, dtype=float))/8.0*100.

fig=plt.figure(figsize=(14.4, 8.6), facecolor=BG)
ax=fig.add_axes([.055,.085,.70,.735]); ax.set_facecolor(BG)
c=ax.contourf(X, Y, Z, levels=np.linspace(5.6,8.15,64), cmap=WAKE_CMAP, extend="both")
order=np.argsort(-deficit)
sc=ax.scatter(xr, yr, s=34, c=deficit, cmap=DEF_CMAP, vmin=0, vmax=max(12,deficit.max()),
              edgecolors="#ffffff", linewidths=.7, zorder=6)
offs = [(12,12), (12,-18), (-46,14)]
for k, i in enumerate(order[:3]):
    ax.annotate(f"#{i+1}", (xr[i], yr[i]), textcoords="offset points", xytext=offs[k],
                color=AMB, fontsize=12, weight="bold", zorder=8,
                arrowprops=dict(arrowstyle="-", color=AMB, lw=.9, alpha=.75))

# 风向箭头
fu = -math.sin(math.radians(WD)); fv = -math.cos(math.radians(WD))   # 顺风单位向量
L = .085
ax.annotate("", xy=(.055+fu*L, .085+fv*L), xytext=(.055, .085), xycoords="axes fraction",
            arrowprops=dict(arrowstyle="-|>", color=ACC, lw=2.4))
ax.text(.055+fu*L+.018, .085+fv*L, "来流方向", transform=ax.transAxes,
        color=ACC, fontsize=11, va="center")
ax.set_aspect("equal"); ax.tick_params(colors=MUT, labelsize=8.5)
for s in ax.spines.values(): s.set_color(LINE)
ax.set_xlabel("场内相对坐标 (m) — 已做刚性变换脱敏", color=MUT, fontsize=9.5)
ax.set_ylabel("场内相对坐标 (m)", color=MUT, fontsize=9.5)

cax=fig.add_axes([.762,.42,.011,.40]); cb=fig.colorbar(c,cax=cax)
cb.set_label("尾流后有效风速 (m/s)", color=INK2, fontsize=9.5)
cb.ax.yaxis.set_tick_params(color=MUT, labelsize=8.5)
plt.setp(plt.getp(cb.ax.axes,'yticklabels'), color=MUT); cb.outline.set_edgecolor(LINE)
cax2=fig.add_axes([.845,.42,.011,.40]); cb2=fig.colorbar(sc,cax=cax2)
cb2.set_label("该机位的风速亏损 (%)", color=INK2, fontsize=9.5)
cb2.ax.yaxis.set_tick_params(color=MUT, labelsize=8.5)
plt.setp(plt.getp(cb2.ax.axes,'yticklabels'), color=MUT); cb2.outline.set_edgecolor(LINE)

# 风玫瑰 (实测)
rose=[0.6,2.5,2.7,4.0,14.2,33.5,26.6,8.1,2.7,2.2,2.4,0.5]
axr=fig.add_axes([.775,.085,.185,.245], projection="polar"); axr.set_facecolor(BG)
th=np.radians(np.arange(0,360,30)+15)
axr.bar(th, rose, width=np.radians(29), color=ACC, alpha=.72, edgecolor=LINE)
axr.bar(np.radians(165), 33.5, width=np.radians(29), color=AMB, edgecolor="#8a5510")
axr.set_theta_zero_location("N"); axr.set_theta_direction(-1)
axr.set_yticklabels([]); axr.tick_params(colors=MUT, labelsize=7.5)
axr.grid(color=LINE, lw=.6); axr.spines["polar"].set_color(LINE)
axr.set_title("实测风向分布 · 主导 165°", color=INK2, fontsize=9.5, pad=6)

fig.text(.055,.955,"一个真实山地风场的尾流图", color=INK, fontsize=20, weight="bold")
fig.text(.055,.912,f"{len(xr)} 台机位为实际坐标（两期、三家整机厂混排）；风向分布来自该场 1,218,810 条实测 10 分钟记录，主导风向占 33.5%。",
         color=INK2, fontsize=12.2)
fig.text(.62,.955,"标注 = 该风向下尾流最深的 3 台", color=AMB, fontsize=10.5)
fig.text(.055,.884,"已脱敏：坐标做刚性变换（相对几何不变，尾流物理为真），台号重编，不含场名与厂商。"
                   "边界：平坦地形假设，推力曲线为通用示意值 → 只看相对形态与排序，不报绝对损失。",
         color=MUT, fontsize=9.3)
fig.text(.055,.018,"引擎：PyWake（丹麦技术大学，MIT 许可）· Bastankhah–Porté-Agel 2014 高斯尾流模型",
         color=MUT, fontsize=8.8)
out=os.path.join(os.path.dirname(__file__),"wake_real.png")
fig.savefig(out, dpi=150, facecolor=BG); print("SAVED", out)
print("n_turbines", len(xr), "| deficit max %.1f%% mean %.1f%%" % (deficit.max(), deficit.mean()))