短期气候预测实习5

xiaoxiao啊

基于jupyter notebook
数据获取方式如下
链接：https://pan.baidu.com/s/17TSNCN-XXV5zrMXhIG9iTA
提取码：0un5


import numpy as np
import pandas as pd
import xarray as xr
from scipy.stats.mstats import ttest_ind
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


time = pd.date_range('1947-01-01', periods=758, freq='M')
lat = np.linspace(-88, 88, 89)
lon = np.linspace(0, 358, 180)


f = np.fromfile('sst.grb', '<f4').reshape(758, 89, 180)
f = np.where(f == 32767, np.nan, f)
sst = xr.DataArray(f, coords=[time, lat, lon], dims=['time', 'lat', 'lon'])
sst = sst.loc[sst.time.dt.month.isin(12), -10:60, 120:300].loc['1950':].assign_coords(time=range(1951, 2011))


f2 = np.loadtxt('ddi.txt', dtype='int32')
_, pr_type = np.where(f2 == 1)
year = np.arange(1951, 2011)
yr1 = year[pr_type == 0]
yr2 = year[pr_type == 1]
yr3 = year[pr_type == 2]


sst1 = sst.loc[sst.time.isin(yr1)]
sst2 = sst.loc[sst.time.isin(yr2)]
sst3 = sst.loc[sst.time.isin(yr3)]


_, p12 = ttest_ind(sst1, sst2)
_, p23 = ttest_ind(sst2, sst3)
_, p13 = ttest_ind(sst1, sst3)


mesh_lon, mesh_lat = np.meshgrid(sst.lon, sst.lat)


def sst_plot():
#     levels = np.arange(-1, 1.2, 0.2)
    levels = np.arange(-0.5, 0.6, 0.1)
    fig = plt.figure(figsize=(8, 6), facecolor='w')
    ax = fig.add_subplot(111, projection=ccrs.PlateCarree(210))
    ax.add_feature(cfeature.LAND.with_scale('110m'), facecolor='grey', zorder=3)
    ax.set_extent([120, 300, -10, 60], ccrs.PlateCarree())
    ax.set_xticks(np.arange(120, 320, 20), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(-10, 70, 10), 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)
    cf = ax.contourf(sst.lon, sst.lat, sst3.mean('time')-sst.mean('time'), transform=ccrs.PlateCarree(),
                     cmap=cmaps.BlueDarkRed18, extend='both', zorder=1, levels=levels)   # 该处变量更具需求自己更换
    cb = fig.colorbar(cf, orientation='horizontal', aspect=35, shrink=0.9, ticks=levels, pad=0.1)
    cb.outline.set_linewidth(0.1)
    cb.ax.tick_params(length=0)
    # 显著性打点
#     dot_area = np.where(p12 < 0.05)
#     dot = ax.scatter(mesh_lon[dot_area], mesh_lat[dot_area], color='k', s=2, linewidths=0, zorder=2, transform=ccrs.PlateCarree())
sst_plot()