气象统计方法实习6

xiaoxiao啊

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


import numpy as np
import xarray as xr
from eofs.standard import Eof
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
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

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)

coslat = np.cos(np.deg2rad(tpsst.lat.values))
solver = Eof(tpsst.values, weights=coslat[..., np.newaxis])
eof = solver.eofs(neofs=10)
pc = solver.pcs(npcs=10)
var = solver.varianceFraction(neigs=10)

def eof_plot():
    fig = plt.figure(figsize=(12, 10), facecolor='w')
    clevels = np.arange(-.25, .3, .05)
    for j in range(2):
        ax = fig.add_axes([0, 1 - j * 0.28, 1, 0.35], projection=ccrs.PlateCarree(210))
        ax.set_title(f'EOF{j + 1}', loc='left')
        ax.set_title(f'{round(var[j] * 100, 2)}%', loc='right')
        ax.add_feature(cfeature.LAND.with_scale('110m'), facecolor='w', edgecolor='dimgrey', linewidth=0.3, zorder=2)
        ax.set_extent([120, 300, -20, 20], ccrs.PlateCarree())
        ax.set_xticks(np.arange(120, 300 + 40, 40), crs=ccrs.PlateCarree())
        ax.set_yticks(np.arange(-20, 20 + 10, 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(2)
        cf = ax.contourf(tpsst.lon, tpsst.lat, eof[j], cmap=cmaps.NCV_blue_red, levels=clevels, transform=ccrs.PlateCarree(), zorder=1)
        
        ax2 = make_axes_locatable(ax).new_horizontal(size='50%', pad=-1.2, axes_class=plt.Axes)
        fig.add_axes(ax2)
        ax2.set_title(f'PC{j + 1}', loc='left')
        bars = ax2.bar(range(1978, 2008), pc[:, j], facecolor='#ff9999', edgecolor='r', width=1)
        for bar, height in zip(bars, pc[:, j]):
            if height < 0:
                bar.set(facecolor='#9999ff', edgecolor='blue')
        for i in ['bottom', 'top', 'left', 'right']:
            ax2.spines【i】.set_linewidth(0.5)     # 记得改成英文的中括号，气象家园无法显示英文中括号加i
        ax2.tick_params(which='major', width=0.5, length=5)
        ax2.set_xlabel('year')
        ax2.set_xlim(1977, 2008)
   
    cax = fig.add_axes([0.19, 0.7, 0.35, 0.025])
    cb = fig.colorbar(cf, cax=cax, orientation='horizontal', ticks=clevels)
    cb.outline.set_linewidth(0.1)
    cb.ax.tick_params(length=0)

    plt.savefig('eof.png', dpi=300, bbox_inches='tight')
eof_plot()

pc1 = pc[:, 0]
pc1_norm = (pc1 - pc1.mean()) / pc1.std()
year = np.arange(1978, 1978 + 30)
print('El Nino years: ', year[pc1_norm >= 0.5])
print('La Nina years: ', year[pc1_norm <= -0.5])

nino_yr = np.array([1982, 1986, 1987, 1991, 1994, 1997, 2002, 2006])
nina_yr = np.array([1984, 1985, 1988, 1995, 1996, 1998, 1999, 2000, 2005, 2007])
tpsst_nino = tpsst.loc[tpsst.time.isin(nino_yr)]
tpsst_nina = tpsst.loc[tpsst.time.isin(nina_yr)]
tpsst_ano = tpsst_nino.mean('time') - tpsst_nina.mean('time')

def tpsst_ano_plot():
    fig = plt.figure(figsize=(6, 4), facecolor='w')
    ax = fig.add_subplot(projection=ccrs.PlateCarree(210))
    ax.set_title('Winter El Nino & La Nina TPSST Ano', loc='left')
    ax.add_feature(cfeature.LAND.with_scale('110m'), facecolor='w', edgecolor='dimgrey', linewidth=0.3, zorder=2)
    ax.set_extent([120, 300, -20, 20], ccrs.PlateCarree())
    ax.set_xticks(np.arange(120, 300 + 40, 40), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(-20, 20 + 10, 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(2)
    clevels = np.arange(-3.5, 4, 0.5)
    cf = ax.contourf(tpsst.lon, tpsst.lat, tpsst_ano, cmap=cmaps.NCV_blue_red, levels=clevels, transform=ccrs.PlateCarree(), zorder=1)
    cb = fig.colorbar(cf, orientation='horizontal', ticks=clevels, shrink=0.9)
    cb.outline.set_linewidth(0.1)
    cb.ax.tick_params(length=0)
    plt.savefig('tpsst_ano.png', dpi=300, bbox_inches='tight')
tpsst_ano_plot()

# 重建
tpsst_re = np.dot(pc, eof.reshape(10, 7 * 32)).reshape(30, 7, 32)
tpsst_re = xr.DataArray(tpsst_re, coords=[time, lat, lon], dims=['time', 'lat', 'lon'])
tpsst_re = tpsst_re.where(tpsst_re!=1e+33, np.nan)

def tpsst_98re_plot():
    fig = plt.figure(figsize=(6, 4), facecolor='w')
    ax = fig.add_subplot(projection=ccrs.PlateCarree(210))
    ax.set_title('1998 Winter TPSST Rebuild', loc='left')
    ax.add_feature(cfeature.LAND.with_scale('110m'), facecolor='w', edgecolor='dimgrey', linewidth=0.3, zorder=2)
    ax.set_extent([120, 300, -20, 20], ccrs.PlateCarree())
    ax.set_xticks(np.arange(120, 300 + 40, 40), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(-20, 20 + 10, 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(2)
    clevels = np.arange(-2.5, 3, 0.5)
    cf = ax.contourf(tpsst.lon, tpsst.lat, tpsst_re.loc[1998], cmap=cmaps.NCV_blue_red, levels=clevels, transform=ccrs.PlateCarree(), zorder=1)
    cb = fig.colorbar(cf, orientation='horizontal', ticks=clevels, shrink=0.9)
    cb.outline.set_linewidth(0.1)
    cb.ax.tick_params(length=0)
    plt.savefig('tpsst_98re.png', dpi=300, bbox_inches='tight')
tpsst_98re_plot()

def tpsst_98_plot():
    fig = plt.figure(figsize=(6, 4), facecolor='w')
    ax = fig.add_subplot(projection=ccrs.PlateCarree(210))
    ax.set_title('1998 Winter TPSST ', loc='left')
    ax.add_feature(cfeature.LAND.with_scale('110m'), facecolor='w', edgecolor='dimgrey', linewidth=0.3, zorder=2)
    ax.set_extent([120, 300, -20, 20], ccrs.PlateCarree())
    ax.set_xticks(np.arange(120, 300 + 40, 40), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(-20, 20 + 10, 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(2)
    clevels = np.arange(-2.5, 3, 0.5)
    cf = ax.contourf(tpsst.lon, tpsst.lat, (tpsst - tpsst.mean('time')).loc[1998], cmap=cmaps.NCV_blue_red,
                     levels=clevels, transform=ccrs.PlateCarree(), zorder=1)
    cb = fig.colorbar(cf, orientation='horizontal', ticks=clevels, shrink=0.9)
    cb.outline.set_linewidth(0.1)
    cb.ax.tick_params(length=0)
    plt.savefig('tpsst_98.png', dpi=300, bbox_inches='tight')
tpsst_98_plot()

lxy287131416

牛啊牛啊 学习一下

_丢丢_

非常感谢！