import xarray as xr import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs def readdat(pth,var): ds = xr.open_dataset(pth) dat = ds[var] lon=ds['lon'] lat=ds['lat'] # 只读取每年1月的数据 dat = dat.sel(time=dat['time.month'].isin([1])) dat = dat.loc['1981-01-01':'2018-12-01','500',:,:] return lon,lat,dat if __name__=='__main__': pth=r'D:\6termsx\qhyc\sx1\hgt.mon.mean.nc' lon,lat,hgt=readdat(pth,'hgt') #print(hgt.shape) #年 纬 经 #print(hgt) wp = np.zeros(38) # 初始化空数组 for i in range(38): wp[i] =0.5*(hgt[i,12,62]-hgt[i,24,62]) wpm=np.mean(wp) #平均 wpa=wp[:]-np.tile(wpm,(wp.shape))#中心化 #print(wpa.shape) wps=wpa/np.tile(np.std(wpa),(wpa.shape))#标准化 #print(wps) #计算WP遥相关指数与同期环流场（500hPa高度场）的相关系数 for i in range(73): for j in range(144): r=np.corrcoef(hgt[:,i,j],wps)[1,0]画出r的全球分布

以下是代码实现：

import xarray as xr 
import numpy as np 
import matplotlib.pyplot as plt 
import cartopy.crs as ccrs 

def readdat(pth,var):
    ds = xr.open_dataset(pth) 
    dat = ds[var] 
    lon=ds['lon'] 
    lat=ds['lat'] 
    # 只读取每年1月的数据 
    dat = dat.sel(time=dat['time.month'].isin([1])) 
    dat = dat.loc['1981-01-01':'2018-12-01','500',:,:] 
    return lon,lat,dat 

if __name__=='__main__':
    pth=r'D:\6termsx\qhyc\sx1\hgt.mon.mean.nc' 
    lon,lat,hgt=readdat(pth,'hgt') 
    wp = np.zeros(38) # 初始化空数组 
    for i in range(38): 
        wp[i] =0.5*(hgt[i,12,62]-hgt[i,24,62]) 
    wpm=np.mean(wp) #平均 
    wpa=wp[:]-np.tile(wpm,(wp.shape))#中心化 
    wps=wpa/np.tile(np.std(wpa),(wpa.shape))#标准化 
    
    #计算WP遥相关指数与同期环流场（500hPa高度场）的相关系数 
    r_map = np.zeros((73,144))
    for i in range(73):
        for j in range(144):
            r_map[i,j] = np.corrcoef(hgt[:,i,j],wps)[1,0]
    
    #画出r的全球分布
    fig = plt.figure(figsize=(8,6))
    ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree())
    ax.coastlines()
    im = ax.contourf(lon,lat,r_map,transform=ccrs.PlateCarree(),cmap='RdBu_r',extend='both')
    cax = plt.colorbar(im,ax=ax)
    cax.ax.set_ylabel('Correlation Coefficient')
    ax.set_title('WP Teleconnection Pattern')
    plt.show()

输出结果：