气象统计方法实习5

xiaoxiao啊

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


import numpy as np
import xarray as xr
from scipy.stats import linregress
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
import cmaps

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

# IOSST
lon = np.arange(37.5, 37.5 + 5.65 * 16, 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)
iosst = np.fromfile('NCEP_IOSST_30y_Wt.dat', dtype='<f4').reshape(30, 7, 16)
iosst = xr.DataArray(iosst, coords=[time, lat, lon], dims=['time', 'lat', 'lon'])
iosst = iosst.where(iosst!=1e+33, np.nan)
IOB = iosst.mean(['lat', 'lon'])
IOBA = IOB - IOB.mean('time')

# 线性趋势
k, b, _, _, _ = linregress(time, IOBA)
y = k * time + b

# 11点滑动平均
IOBA_11 = np.convolve(IOBA, np.ones(11), 'valid') / 11

# 累计距平
x = []
for i in range(30):
    x.append(float(IOBA[:i].sum()))

def my_plot():
    fig = plt.figure(figsize=(4, 8), facecolor='w')
    for j in range(3):
        ax = fig.add_axes([0.1, 1 - 0.3 * j, 0.8, 0.21])
        if j == 0:
            ax.set_title('Linear Detrend', loc='left')
            ax.plot(time, y, c='purple', lw=1)
        elif j == 1:
            ax.set_title('Moveing Average', loc='left')
            ax.plot(range(1983, 1983 + 20), IOBA_11, c='purple', lw=1)
        else:
            ax.set_title('Cumulative Ano', loc='left')
            ax.plot(time, x, c='purple', lw=1)
        bars = ax.bar(time, IOBA, facecolor='#ff9999', width=.9)
        for bar, height in zip(bars, IOBA):
            if height < 0:
                bar.set(facecolor='#9999ff')
        ax.tick_params(which='major', width=0.5, length=5)
        for i in ['bottom', 'top', 'left', 'right']:
            ax.spines【i】.set_linewidth(0.5)     # 记得改成英文的中括号，气象家园无法显示英文中括号加i
        ax.set_xlim(1977, 2008)
        plt.savefig('3_types.png', dpi=300, bbox_inches='tight')
my_plot()

lon = np.arange(1.875, 1.875 + 191 * 1.875, 1.875)
lat = np.array([-59.9986, -58.0939, -56.1893, -54.2846, -52.3799, -50.4752, -48.5705, -46.6658, -44.7611, -42.8564, -40.9517, -39.047, -37.1422, -35.2375, -33.3328, -31.4281, -29.5234, -27.6186, -25.7139, -23.8092, -21.9044, -19.9997, -18.095, -16.1902, -14.2855, -12.3808, -10.476, -8.57131, -6.66657, -4.76184, -2.8571, -0.952368, 0.952368, 2.8571, 4.76184, 6.66657, 8.57131, 10.476, 12.3808, 14.2855, 16.1902, 18.095, 19.9997, 21.9044, 23.8092, 25.7139, 27.6186, 29.5234, 31.4281, 33.3328, 35.2375, 37.1422, 39.047, 40.9517, 42.8564, 44.7611, 46.6658, 48.5705, 50.4752, 52.3799, 54.2846, 56.1893, 58.0939, 59.9986])
time = np.arange(1978, 1978 + 50)
gsst = np.fromfile('NCEP_GSST_50y_Wt.dat', dtype='<f4').reshape(50, 64, 191)
gsst = xr.DataArray(gsst, coords=[time, lat, lon], dims=['time', 'lat', 'lon'])
gsst = gsst.where(gsst!=1e+33, np.nan)

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

time_3dims = to_3dims(var=time)
gsst_k = linslope(time_3dims, gsst, dim='time')
gsst_p = pearson_r_p_value(time_3dims, gsst, dim='time')

new_lon = np.linspace(0, 360, 382)
new_lat = np.linspace(-60, 60, 144)
mesh_lon, mesh_lat = np.meshgrid(lon, lat)

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

new_gsst_k = my_interploate(var=gsst_k.fillna(0))

def linearRegress_plot():
    fig = plt.figure(figsize=(6, 6), facecolor='w')
    ax = fig.add_subplot(111, projection=ccrs.Robinson(180))
    ax.set_title('Winter Global SST Linear Detrend')
    ax.add_feature(cfeature.LAND.with_scale('110m'), facecolor='w', edgecolor='dimgrey', linewidth=0.5, zorder=1)
    ax.set_global()
    ax.spines['geo'].set_linewidth(0.5)
    ax.set_aspect(1.2)
    clevels = np.arange(-.05, .06, .01)
    cf = ax.contourf(new_lon, new_lat, new_gsst_k, cmap=cmaps.NCV_blue_red, levels=clevels, extend='both', transform=ccrs.PlateCarree())
    cb = fig.colorbar(cf, orientation='horizontal', shrink=0.85, aspect=35, pad=0.03, ticks=clevels)
    cb.outline.set_linewidth(0.1)
    cb.ax.tick_params(length=0)
    cb.set_label(r'$[^{\circ}C\ /\ year]$', fontsize=8)
    dot_area = np.where(gsst_p < 0.05)
    dot = ax.scatter(mesh_lon[dot_area], mesh_lat[dot_area], color='k', s=0.8, linewidths=0, transform=ccrs.PlateCarree())
    plt.savefig('gsst_linearRegress.png', dpi=300, bbox_inches='tight')
linearRegress_plot()