短期气候预测实习2

xiaoxiao啊

基于jupyter notebook


# 导入所需库
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


# 更改字体
plt.rcParams['font.sans-serif'] = ['Times New Roman']
plt.rcParams['axes.unicode_minus'] = False


# 处理数据
f = xr.open_dataset('hgt.mon.mean.nc')  # 读取数据集
hgt_Eu= f.hgt.loc[f.time.dt.month.isin([1]), 500, 70:20, 40:140].loc['1948':'2008']  # 欧亚地区一月高度场


# EOF正交分解
coslat = np.cos(np.deg2rad(hgt_Eu.lat.values))  # 纬度权重
solver = Eof(hgt_Eu.values, weights=coslat[..., np.newaxis])  # EOF分解器
eof = solver.eofs(neofs=3)  # 模态场
pc = solver.pcs(npcs=3)  # 时间序列
var = solver.varianceFraction(neigs=3)  # 解释方差


# 绘画EOF三个模态及其时间序列
def eof_plot():
    fig = plt.figure(figsize=(12, 10), facecolor='w')
    for j in range(3):
        ax1 = fig.add_axes([0.1, 1 - 0.3 * j, 0.8, 0.20], projection=ccrs.PlateCarree())
        ax1.set_title(f'EOF{j + 1}', loc='left')
        ax1.set_title(f'{round(var[j] * 100, 1)}%', loc='right')
        ax1.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.6)
        ax1.add_feature(cfeature.LAKES.with_scale('110m'), facecolor='none', edgecolor='k', linewidth=0.6)
        ax1.set_extent([40, 140, 20, 70], crs=ccrs.PlateCarree())
        ax1.set_xticks(np.arange(40, 160, 20), crs=ccrs.PlateCarree())
        ax1.set_yticks(np.arange(20, 80, 10), crs=ccrs.PlateCarree())
        ax1.xaxis.set_major_formatter(LongitudeFormatter())
        ax1.yaxis.set_major_formatter(LatitudeFormatter())
        ax1.tick_params(which='major', width=0.5, length=5)
        ax1.spines['geo'].set_linewidth(0.5)
        ax1.set_aspect(1.2)
        m1 = ax1.contourf(hgt_Eu.lon, hgt_Eu.lat, eof[j], extend='both', cmap='bwr', transform=ccrs.PlateCarree())

        ax2 = make_axes_locatable(ax1).new_horizontal(size='80%', pad=0.5, axes_class=plt.Axes)
        fig.add_axes(ax2)
        ax2.set_title(f'PC{j + 1}', loc='left')
        bars = ax2.bar(range(1948, 2009), 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_xlim(1947, 2009)

eof_plot()