温度平流分析求助

liuna

请问输出的温度平流图冷暖混杂严重是什么原因噻，求大佬掌眼，拜托
import xarray as xr
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import cartopy.mpl.ticker as cticker


class T_advection(object):
    def __init__(self, file_u, file_v, file_t):
        self.f_u = file_u
        self.f_v = file_v
        self.f_t = file_t

    def cal_advection(self):
        a = 6370000  # Earth radius in meters
        print("\n==== 开始计算温度平流 ====")

        # Use the first time step of the selected data
        u0 = self.f_u.isel(time=0)
        v0 = self.f_v.isel(time=0)
        t = self.f_t.isel(time=0)

        # 提取坐标（确保lon0/lat0在打印前定义）
        lon0 = u0.lon
        lat0 = u0.lat

        # 检查数据维度和坐标
        print("\n--- 1. 数据基本信息 ---")
        print(f"u0 维度: {u0.dims}, 形状: {u0.shape}")
        print(f"v0 维度: {v0.dims}, 形状: {v0.shape}")
        print(f"t 维度: {t.dims}, 形状: {t.shape}")

        # 检查经纬度坐标
        print("\n--- 2. 经纬度坐标检查 ---")
        print(f"经度范围: {lon0.min().values:.2f}°E 到 {lon0.max().values:.2f}°E")
        print(f"纬度范围: {lat0.min().values:.2f}°N 到 {lat0.max().values:.2f}°N")

        # 检查纬度方向（是否南→北或北→南）
        if lat0.values[0] < lat0.values[1]:
            print("纬度方向: 南→北（从小到大）")
        else:
            print("纬度方向: 北→南（从大到小）")

        # 计算网格间距
        print("\n--- 4. 计算网格间距 ---")
        dlon = (np.gradient(lon0, axis=1) * np.pi / 180).reshape((1, -1))
        dlat = (np.gradient(lat0, axis=0) * np.pi / 180).reshape((-1, 1))
        print(f"经度网格间距范围: {dlon.min():.6f} 到 {dlon.max():.6f} 弧度")
        print(f"纬度网格间距范围: {dlat.min():.6f} 到 {dlat.max():.6f} 弧度")

        coslat = (np.cos(lat0 * np.pi / 180)).reshape((-1, 1))
        dx = a * coslat * dlon
        dy = a * dlat
        print(f"dx 范围（米）: {dx.min():.2f} 到 {dx.max():.2f}")
        print(f"dy 范围（米）: {dy.min():.2f} 到 {dy.max():.2f}")

        # 计算温度梯度
        print("\n--- 5. 计算温度梯度 ---")
        t_np = np.array(t)
        dT_dx = np.gradient(np.array(t), axis=1) / dx
        dT_dy = np.gradient(np.array(t), axis=0) / dy
        print(f"dT_dx 范围 (K/m): {dT_dx.min():.6e} 到 {dT_dx.max():.6e}")
        print(f"dT_dy 范围 (K/m): {dT_dy.min():.6e} 到 {dT_dy.max():.6e}")

        # 计算温度平流
        print("\n--- 6. 计算温度平流 ---")
        temadvection0 = -(u0 * dT_dx + v0 * dT_dy)
        print(f"温度平流范围 (K/s): {temadvection0.min():.6e} 到 {temadvection0.max():.6e}")
        print(f"温度平流平均值 (K/s): {temadvection0.mean():.6e}")

        return lon0, lat0, u0, v0, temadvection0

    def draw_cmap_advection(self, lon, lat, u, v, temadvection):
        # Set projection extent
        proj = ccrs.PlateCarree(central_longitude=90)
        leftlon, rightlon, lowerlat, upperlat = (0, 180, 0, 90)
        img_extent = [leftlon, rightlon, lowerlat, upperlat]

        # Add basic features
        fig = plt.figure(figsize=(12, 8))
        ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=proj)
        ax.set_extent(img_extent, crs=ccrs.PlateCarree())
        ax.set_title('Nov 2014 Temperature Advection', loc='center', fontsize=18)
        ax.set_title(r'units: 10$^{-5}$ K/s', loc='right', fontsize=12)
        ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
        ax.add_feature(cfeature.LAKES, alpha=0.5)
        ax.set_xticks(np.array([0, 60, 120, 180]), crs=ccrs.PlateCarree())
        ax.set_yticks(np.array([0, 30, 60, 90]), crs=ccrs.PlateCarree())
        ax.xaxis.set_major_formatter(cticker.LongitudeFormatter())
        ax.yaxis.set_major_formatter(cticker.LatitudeFormatter())

        # Add data (quiver, contourf, colorbar)
        q1 = ax.quiver(lon[::8], lat[::8], u.values[::8, ::8], v.values[::8, ::8],
                       pivot='mid', scale=150, color='g', width=0.0038,
                       headwidth=3, transform=ccrs.PlateCarree())
        ax.quiverkey(q1, X=0.9, Y=-0.12, U=15, label='15 m/s', labelpos='E')

        c1 = ax.contourf(lon, lat, temadvection / 1e-05, zorder=0,
                         levels=np.arange(-5, 5.5, 0.5), extend='both',
                         transform=ccrs.PlateCarree(), cmap=plt.cm.bwr)

        position = fig.add_axes([0.12, 0.12, 0.4, 0.02])
        fig.colorbar(c1, cax=position, orientation='horizontal', format='%.2f')

        plt.show()
        print("\n==== 温度平流图已绘制 ====")


# Load and prepare data
def load_and_prepare_data():
    # Load data (replace with your actual file paths)
    print("\n==== 开始加载数据 ====")
    udata = xr.open_dataset(r'F:\BNU\data load\Arctic\wind\uwnd.mon.mean.nc')['uwnd']
    vdata = xr.open_dataset(r'F:\BNU\data load\Arctic\wind\vwnd.mon.mean.nc')['vwnd']
    tdata = xr.open_dataset(r'F:\BNU\data load\Arctic\hgt\daily\air.mon.mean-NCEP20.nc')['air']
    print("数据加载完成")

    # Select 850 hPa level
    print("\n--- 选择 850 hPa 层级 ---")
    udata_850 = udata.sel(level=850)
    vdata_850 = vdata.sel(level=850)
    tdata_850 = tdata.sel(level=850)
    print(udata_850)
    print(vdata_850)
    print(tdata_850)

    # Select November 2014 data
    print("\n--- 选择 2014 年 11 月数据 ---")
    udata_2014_11 = udata_850.sel(time=slice('2014-11-01', '2014-11-30'))
    vdata_2014_11 = vdata_850.sel(time=slice('2014-11-01', '2014-11-30'))
    tdata_2014_11 = tdata_850.sel(time=slice('2014-11-01', '2014-11-30'))

    print(udata_2014_11)
    print(vdata_2014_11)
    print(tdata_2014_11)

    return udata_2014_11, vdata_2014_11, tdata_2014_11


# Main execution
if __name__ == "__main__":
    try:
        # Load and prepare data
        u_data, v_data, t_data = load_and_prepare_data()

        # Create T_advection instance and calculate
        print("\n==== 开始创建 T_advection 实例并计算 ====")
        t_adv = T_advection(u_data, v_data, t_data)
        lon, lat, u, v, temadvection = t_adv.cal_advection()

        # Plot the results
        print("\n==== 开始绘制温度平流图 ====")
        t_adv.draw_cmap_advection(lon, lat, u, v, temadvection)

    except Exception as e:
        print(f"\n==== 错误 ====")
        print(f"错误信息: {str(e)}")
        print("请检查数据或计算步骤")