气象统计方法实习3

xiaoxiao啊

气象统计方法所有数据获取方式如下
链接：https://pan.baidu.com/s/1KbxLS8nmvpqkxbVt5Q-iHA
提取码：wumq


import numpy as np
import xarray as xr
from scipy.interpolate import interp2d
from xskillscore import linslope, pearson_r_p_value
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cmaps

plt.rcParams['font.sans-serif'] = ['Times New Roman']
plt.rcParams['axes.unicode_minus'] = False

# TPSST
lon = np.arange(121.88, 121.88 + 32 * 5.65, 5.65)
lat = np.array([-18.0950, -12.3808, -6.6666, -0.9524, 4.7618, 10.4761, 16.1902])
time = np.arange(1978, 1978 + 30)
tpsst = np.fromfile('NCEP_TPSST_30y_Wt.dat', dtype='<f4').reshape(30, 7, 32)
tpsst = xr.DataArray(tpsst, coords=[time, lat, lon], dims=['time', 'lat', 'lon'])
tpsst = tpsst.where(tpsst!=1e+33, np.nan)
nino34 = tpsst.loc[:, -5:5, 190:240].mean(['lat', 'lon'])  # nino3.4区海温指数

# SLP
lon = np.linspace(5, 357.5, 48)
lat = np.linspace(-87.5, 85, 24)
time = np.arange(1978, 1978 + 30)
slp = np.fromfile('NCEP_slp_30y_Wt.dat', dtype='<f4').reshape(30, 24, 48)
slp = xr.DataArray(slp, coords=[time, lat, lon], dims=['time', 'lat', 'lon'])

def to_3dims(*, var):
    var = np.array(var)
    var = np.expand_dims(var, axis=1)
    var = np.repeat(var, 24 * 48, axis=1)
    var = np.reshape(var, [30, 24, 48])
    var = xr.DataArray(var, coords=[time, lat, lon], dims=['time', 'lat', 'lon'])
    return var

new_lon = np.linspace(0, 360, 144)
new_lat = np.linspace(-90, 90, 72)
mesh_new_lon, mesh_new_lat = np.meshgrid(new_lon, new_lat)

def my_interploate(*, var):
    func = interp2d(lon, lat, var, kind='cubic')
    return func(new_lon, new_lat)

nino34_3dims = to_3dims(var=nino34)
k = linslope(slp, nino34_3dims, dim='time')  # slp为预报因子，nino34为预报量
p = pearson_r_p_value(slp, nino34_3dims, dim='time')
new_k = my_interploate(var=k)
new_p = my_interploate(var=p)

def linearRegress_plot():
    fig = plt.figure(figsize=(8, 5), facecolor='w')
    ax = fig.add_subplot(111, projection=ccrs.PlateCarree(180))
    ax.set_title('Winter Nino34 & Slp Linear Regress', loc='left')
    ax.add_feature(cfeature.COASTLINE.with_scale('110m'), edgecolor='dimgrey', linewidth=0.5)
    ax.set_extent([0, 360, -90, 90], ccrs.PlateCarree())
    ax.set_xticks(np.arange(0, 360, 60), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(-90, 120, 30), crs=ccrs.PlateCarree())
    ax.xaxis.set_major_formatter(LongitudeFormatter())
    ax.yaxis.set_major_formatter(LatitudeFormatter())
    ax.tick_params(which='major', width=0.5, length=5)
    ax.spines['geo'].set_linewidth(0.5)
    ax.set_aspect(1.2)
    clevels = np.arange(-2, 2.4, 0.4)
    cf = ax.contourf(new_lon, new_lat, new_k, cmap=cmaps.NCV_blue_red, levels=clevels, transform=ccrs.PlateCarree())
    cb = fig.colorbar(cf, orientation='horizontal', shrink=0.6, aspect=20, pad=0.12, ticks=clevels)
    cb.outline.set_linewidth(0.1)
    cb.ax.tick_params(length=0)
    dot_area = np.where(new_p < 0.05)
    dot = ax.scatter(mesh_new_lon[dot_area], mesh_new_lat[dot_area], color='k', s=1, linewidths=0, transform=ccrs.PlateCarree())

    plt.savefig('linearRegress.png', dpi=300, bbox_inches='tight')
linearRegress_plot()