guanlan-src/scripts/sim_hub/wake_anim.py

"""动态尾流仿真: 风向按实测风玫瑰扫过, 看尾流在真实机位上横扫。
布局 = 真实山地风场 90 台机位(已脱敏刚性变换); 风向频率 = 该场 1,218,810 条实测 10min 记录。
引擎 = PyWake (DTU, MIT) · Bastankhah-Porte-Agel 2014。"""
import os, csv, math, numpy as np, matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib import font_manager as fm
from pyproj import Transformer
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
from py_wake.flow_map import XYGrid

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


t=Transformer.from_crs("EPSG:4545","EPSG:4326",always_xy=True)
lat=[];lon=[]
for r in csv.DictReader(open("/Users/yuanying/wind-analytics/app_ETL/configs/contracts/mawang_machine_map.csv",encoding="utf-8-sig")):
    try: a=float(r["lat"]); o=float(r["lon"])
    except Exception: continue
    if 73<o<136 and 18<a<54: lat.append(a); lon.append(o)
    elif 1e5<o<1e6 and 1e6<a<6e6:
        lo,la=t.transform(o,a); lat.append(la); lon.append(lo)
lat=np.array(lat); lon=np.array(lon)
R=6371000.
x=R*np.radians(lon-lon.mean())*math.cos(math.radians(lat.mean())); y=R*np.radians(lat-lat.mean())
TH=math.radians(20.0)                      # 脱敏旋转; 风向同步 **减去** 该角 (已验逐台不变)
xr=x*math.cos(TH)-y*math.sin(TH); yr=x*math.sin(TH)+y*math.cos(TH)
xr-=xr.min(); yr-=yr.min()
ROT=math.degrees(TH)

u=np.arange(0,26.); p=np.clip(((u-3)/8)**3,0,1)*2e6; 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)
wt=WindTurbine(name="w",diameter=118.,hub_height=80.,powerCtFunction=PowerCtTabular(u,p,'w',ct))
wfm=Bastankhah_PorteAgel_2014(UniformSite(p_wd=[1],ti=.11,ws=8.),wt,k=0.0324555)

pad=1400.
gx=np.linspace(xr.min()-pad, xr.max()+pad, 460)
gy=np.linspace(yr.min()-pad, yr.max()+pad, 210)
grid=XYGrid(x=gx, y=gy)

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]     # 实测, 30° 一档, 从 0° 起
def rose_pct(wd_true):
    return ROSE[int((wd_true%360)//30)]

WDS=list(range(90, 271, 5))    # 扫过实测主风区 (90-270° 覆盖约 92% 的时间)
frames=os.path.join(os.path.dirname(__file__),"frames"); os.makedirs(frames, exist_ok=True)
for f in os.listdir(frames): os.remove(os.path.join(frames,f))

for k, wd_true in enumerate(WDS):
    WD=(wd_true-ROT)%360
    sim=wfm(xr,yr,wd=WD,ws=8.)
    fmap=sim.flow_map(grid=grid, wd=WD, ws=8.)
    Z=fmap.WS_eff.squeeze().values
    assert Z.shape==(len(gy),len(gx)), f"方向变了 {Z.shape}"
    dfc=(8.0-np.asarray(sim.WS_eff.squeeze().values,dtype=float))/8.0*100.
    # 不报 dfc.max(): 高斯尾流多重叠加会把有效风速压到近零(模型饱和), 真实尾流恢复更快,
    # 该极值不可作对外数字。改报稳健量: 受影响台数 + 平均亏损。
    n_hit = int((dfc > 10).sum()); dfc_mean = float(dfc.mean())

    fig=plt.figure(figsize=(14.4,7.2), facecolor=BG)
    ax=fig.add_axes([.052,.10,.70,.70]); ax.set_facecolor(BG)
    ax.contourf(gx,gy,Z,levels=np.linspace(5.6,8.15,48),cmap=WAKE_CMAP,extend="both")
    ax.scatter(xr,yr,s=30,c=dfc,cmap=DEF_CMAP,vmin=0,vmax=60,edgecolors="#ffffff",linewidths=.6,zorder=6)
    fu=-math.sin(math.radians(WD)); fv=-math.cos(math.radians(WD)); L=.09
    ax.annotate("",xy=(.06+fu*L,.12+fv*L),xytext=(.06,.12),xycoords="axes fraction",
                arrowprops=dict(arrowstyle="-|>",color=ACC,lw=2.4))
    ax.set_aspect("equal"); ax.set_xlim(gx[0],gx[-1]); ax.set_ylim(gy[0],gy[-1])
    ax.tick_params(colors=MUT,labelsize=8)
    for s in ax.spines.values(): s.set_color(LINE)
    ax.set_xlabel("场内相对坐标 (m) — 已脱敏",color=MUT,fontsize=9)

    axr=fig.add_axes([.782,.10,.175,.40],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=.55,edgecolor=LINE)
    axr.bar(np.radians((wd_true//30)*30+15),ROSE[int((wd_true%360)//30)],width=np.radians(29),
            color=AMB,edgecolor="#8a5510")
    axr.plot([np.radians(wd_true)]*2,[0,36],color=INK,lw=2.0)
    axr.set_theta_zero_location("N"); axr.set_theta_direction(-1)
    axr.set_yticklabels([]); axr.tick_params(colors=MUT,labelsize=7)
    axr.grid(color=LINE,lw=.5); axr.spines["polar"].set_color(LINE)
    axr.set_title("实测风向分布",color=INK2,fontsize=9.5,pad=10)

    fig.text(.052,.945,"风向扫过时，尾流在整个场里横扫",color=INK,fontsize=19,weight="bold")
    fig.text(.052,.900,"90 台真实机位（已脱敏）· 风向频率来自 1,218,810 条实测 10 分钟记录 · 引擎 PyWake",
             color=MUT,fontsize=10)
    fig.text(.775,.745,"来流方向",color=INK2,fontsize=13)
    fig.text(.775,.665,f"{wd_true:d}°",color=INK,fontsize=34,weight="bold")   # 数字用默认字体, 不指定 monospace
    fig.text(.775,.615,f"该风向占全年 {rose_pct(wd_true):.1f}%",color=AMB,fontsize=12.5)
    fig.text(.775,.585,f"受尾流影响 {n_hit} 台 · 平均亏损 {dfc_mean:.1f}%",color=INK2,fontsize=12.5)
    fig.savefig(os.path.join(frames,f"f{k:03d}.png"),dpi=110,facecolor=BG); plt.close(fig)
    if k%8==0: print("frame",k,"/",len(WDS),"wd",wd_true)
print("FRAMES DONE", len(WDS))