气象统计方法实习1

xiaoxiao啊

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


import numpy as np
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cmaps

import warnings
warnings.filterwarnings('ignore')

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)

tpsst_clim = tpsst.mean('time')  # 气候态
tpsst_std = tpsst.std('time') # 均方差
tpsst_82 = tpsst.loc[1982]
tpsst_98 = tpsst.loc[1998]

def tpsst_plot():
    fig = plt.figure(figsize=(10, 9), facecolor='w')
    for j in range(3):
        for i in range(2):
            ax = fig.add_axes([0 + 0.5 * i, 1 - 0.33 * j, 0.4, 0.3], projection=ccrs.PlateCarree(210))
            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)

            if i == 0 and j == 0:
                clevels = np.arange(20, 31, 1)
                ax.set_title('Winter TPSST Clim', loc='left')
                cf = ax.contourf(lon, lat, tpsst_clim, cmap='jet', levels=clevels, transform=ccrs.PlateCarree(), zorder=1)
            if i == 1 and j == 0:
                clevels = np.arange(0, 1.8, 0.2)
                ax.set_title('Winter TPSST Std', loc='left')
                cf = ax.contourf(lon, lat, tpsst_std, cmap='Reds', levels=clevels, transform=ccrs.PlateCarree(), zorder=1)
            if i == 0 and j == 1:
                clevels = np.arange(20, 31, 1)
                ax.set_title('1982 Winter TPSST', loc='left')
                cf = ax.contourf(lon, lat, tpsst_82, cmap='jet', levels=clevels, transform=ccrs.PlateCarree(), zorder=1)
            if i ==1 and j ==1:
                clevels = np.arange(-5, 6, 1)
                ax.set_title('1982 Winter TPSST Ano', loc='left')
                cf = ax.contourf(lon, lat, tpsst_82 - tpsst_clim, cmap=cmaps.NCV_blue_red, levels=clevels, transform=ccrs.PlateCarree(), zorder=1)
            if i == 0 and j == 2:
                clevels = np.arange(20, 31, 1)
                ax.set_title('1998 Winter TPSST', loc='left')
                cf = ax.contourf(lon, lat, tpsst_98, cmap='jet', levels=clevels, transform=ccrs.PlateCarree(), zorder=1)
            if i ==1 and j ==2:
                clevels = np.arange(-2.5, 3, 0.5)
                ax.set_title('1998 Winter TPSST Ano', loc='left')
                cf = ax.contourf(lon, lat, tpsst_98 - tpsst_clim, cmap=cmaps.NCV_blue_red, levels=clevels, transform=ccrs.PlateCarree(), zorder=1)

            cb = fig.colorbar(cf, orientation='horizontal', shrink=0.9, pad=0.17, ticks=clevels)
            cb.outline.set_linewidth(0.1)
            cb.ax.tick_params(length=0)
    plt.savefig('tpsst.png', dpi=300, bbox_inches='tight')
tpsst_plot()

nino34 = tpsst.loc[:, -5:5, 190:240].mean(['lat', 'lon'])  # nino3.4区海温指数
nino34_ano = nino34 - nino34.mean('time')  # 距平值
nino34_norm = nino34_ano / nino34.std('time')  # 标准化

year = np.arange(1978, 1978 + 30)
print('El Nino years: ', year[nino34_norm.values >= 1])
print('La Nina years: ', year[nino34_norm.values <= -1])

def nino34_plot():
    fig = plt.figure(figsize=(5, 10), facecolor='w')

    ax1 = fig.add_axes([0.1, 0.75, 0.8, 0.2])
    ax1.plot(year, nino34, linewidth=1)
    ax1.set_title('Nino3.4 SST Index', loc='left')
    ax1.tick_params(which='major', width=0.5, length=5)
    for i in ['bottom', 'top', 'left', 'right']:
        ax1.spines【i】.set_linewidth(0.5)     # 记得改成英文的中括号，气象家园无法显示英文中括号加i
    ax1.set_xlabel('year')
    ax1.set_xlim(1977, 2008)
    ax1.set_ylim(24.5, 29.5)

    ax2 = fig.add_axes([0.1, 0.45, 0.8, 0.2])
    bars = ax2.bar(year, nino34_ano, facecolor='#ff9999', edgecolor='r', width=1)
    for bar, height in zip(bars, nino34_ano):
        if height < 0:
            bar.set(facecolor='#9999ff', edgecolor='blue')
    ax2.set_title('Nino3.4 SST Index Ano', loc='left')
    ax2.tick_params(which='major', width=0.5, length=5)   
    for i in ['bottom', 'top', 'left', 'right']:
        ax2.spines【i】.set_linewidth(0.5)   # 记得改成英文的中括号，气象家园无法显示中括号加i
    ax2.set_xlabel('year')
    ax2.set_xlim(1977, 2008)
    ax2.set_ylim(-2.5, 2.5)

    ax3 = fig.add_axes([0.1, 0.15, 0.8, 0.2])
    bars = ax3.bar(year, nino34_norm, facecolor='#ff9999', edgecolor='r', width=1)
    for bar, height in zip(bars, nino34_norm):
        if height < 0:
            bar.set(facecolor='#9999ff', edgecolor='blue')
    ax3.set_title('Nino3.4 SST Index Norm', loc='left')
    ax3.tick_params(which='major', width=0.5, length=5)   
    for i in ['bottom', 'top', 'left', 'right']:
        ax3.spines【i】.set_linewidth(0.5)   # 记得改成英文的中括号，气象家园无法显示中括号加i
    ax3.set_xlabel('year')
    ax3.set_xlim(1977, 2008)
    ax3.set_ylim(-2.5, 2.5)
    plt.axhline(1, lw=1, c='red', ls=':')
    plt.axhline(-1, lw=1, c='blue', ls=':')
    plt.savefig('nino34.png', dpi=300, bbox_inches='tight')
nino34_plot()

cookie-o-o

谢谢楼主提供脚本及图例，试着运行了脚本，似乎存在一点点小问题。如ax1.spines.set_linewidth(0.5)报错，似乎应该改为ax1.spines【i】.set_linewidth(0.5)。因为这样改后前面的报错信息就没了。其它几个脚本都存在类似问题。另外，脚本运行时提示有字体问题（不过不影响图形输出）。提醒，前面的【i】实际应为英文格式，但英文格式【i】的帖子保存后并不能显示，故用中文格式【i】替代。也可能楼主帖子里面本来就是正确写法，但由于帖子无法正常显示英文格式的【i】造成的。

puck66

楼主好物，谢谢分享，用到再来取，先踩为敬

xiaoxiao啊

cookie-o-o 发表于 2022-7-7 15:33
谢谢楼主提供脚本及图例，试着运行了脚本，似乎存在一点点小问题。如ax1.spines.set_linewidth(0. ...

对，会自动消失，qaq；感谢指出问题，已经修改并加上注释了；出现字体问题的话是你没有安装该字体；

cookie-o-o

出现字体问题的话是你没有安装该字体；

谢谢指教！