请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 3608|回复: 1

[经验总结] 气象统计方法实习2

[复制链接]

新浪微博达人勋

发表于 2022-7-6 19:24:38 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 xiaoxiao啊 于 2022-7-7 17:08 编辑

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


import numpy as np
import xarray as xr
from scipy.stats.mstats import ttest_ind
from scipy.interpolate import interp2d
from xskillscore import pearson_r, pearson_r_p_value
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

nino_yr = np.array([1982, 1986, 1991, 1997, 2002])
nina_yr = np.array([1984, 1988, 1998, 1999, 2007])

lon = np.linspace(5, 357.5, 48)
lat = np.linspace(-87.5, 85, 24)
time = np.arange(1978, 1978 + 30)
slp = np.fromfile('NCEP_slp_30y_Wt.dat', dtype='<f4').reshape(30, 24, 48)
slp = xr.DataArray(slp, coords=[time, lat, lon], dims=['time', 'lat', 'lon'])

slp_nino = slp.loc[slp.time.isin(nino_yr)]
slp_nina = slp.loc[slp.time.isin(nina_yr)]
slp_ano = slp_nino.mean('time') - slp_nina.mean('time')
_, slp_p = ttest_ind(slp_nino, slp_nina, axis=0, equal_var=False)

new_lon = np.linspace(0, 360, 144)
new_lat = np.linspace(-90, 90, 72)
mesh_new_lon, mesh_new_lat = np.meshgrid(new_lon, new_lat)

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

new_slp_ano = my_interploate(var=slp_ano)
new_slp_p = my_interploate(var=slp_p)

def ano_plot():
    fig = plt.figure(figsize=(8, 5), facecolor='w')
    ax = fig.add_subplot(111, projection=ccrs.PlateCarree(180))
    ax.set_title('Winter El Nino & La Nina Slp Ano', loc='left')
    ax.add_feature(cfeature.COASTLINE.with_scale('110m'), edgecolor='dimgrey', linewidth=0.5)
    ax.set_extent([0, 360, -90, 90], ccrs.PlateCarree())
    ax.set_xticks(np.arange(0, 360, 60), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(-90, 120, 30), 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)
    clevels = np.arange(-15, 18, 3)
    cf = ax.contourf(new_lon, new_lat, new_slp_ano, cmap=cmaps.NCV_blue_red, levels=clevels, transform=ccrs.PlateCarree())
    cb = fig.colorbar(cf, orientation='horizontal', shrink=0.6, aspect=20, pad=0.12, ticks=clevels)
    cb.outline.set_linewidth(0.1)
    cb.ax.tick_params(length=0)
    dot_area = np.where(new_slp_p < 0.05)
    dot = ax.scatter(mesh_new_lon[dot_area], mesh_new_lat[dot_area], color='k', s=1, linewidths=0, transform=ccrs.PlateCarree())
    plt.savefig('slp_ano.png', dpi=300, bbox_inches='tight')
ano_plot()

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)
nino34 = tpsst.loc[:, -5:5, 190:240].mean(['lat', 'lon'])  # nino3.4区海温指数

lon = np.linspace(5, 357.5, 48)
lat = np.linspace(-87.5, 85, 24)
time = np.arange(1978, 1978 + 30)

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

nino34_3dims = to_3dims(var=nino34)
r = pearson_r(nino34_3dims, slp, dim='time')
p = pearson_r_p_value(nino34_3dims, slp, dim='time')
new_r = my_interploate(var=r)
new_p = my_interploate(var=p)

def corr_plot():
    fig = plt.figure(figsize=(8, 5), facecolor='w')
    ax = fig.add_subplot(111, projection=ccrs.PlateCarree(180))
    ax.set_title('Winter Nino34 & Slp Corr', loc='left')
    ax.add_feature(cfeature.COASTLINE.with_scale('110m'), edgecolor='dimgrey', linewidth=0.5)
    ax.set_extent([0, 360, -90, 90], ccrs.PlateCarree())
    ax.set_xticks(np.arange(0, 360, 60), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(-90, 120, 30), 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)
    clevels = np.arange(-1, 1.2, 0.2)
    cf = ax.contourf(new_lon, new_lat, new_r, cmap=cmaps.NCV_blue_red, levels=clevels, transform=ccrs.PlateCarree())
    cb = fig.colorbar(cf, orientation='horizontal', shrink=0.6, aspect=20, pad=0.12, ticks=clevels)
    cb.outline.set_linewidth(0.1)
    cb.ax.tick_params(length=0)
    dot_area = np.where(new_p < 0.05)
    dot = ax.scatter(mesh_new_lon[dot_area], mesh_new_lat[dot_area], color='k', s=1, linewidths=0, transform=ccrs.PlateCarree())
    plt.savefig('nino34_slp_corr.png', dpi=300, bbox_inches='tight')
corr_plot()

max_place = r.argmax(dim=['lat', 'lon'])
min_place = r.argmin(dim=['lat', 'lon'])
max_lon_id, max_lat_id = max_place['lon'], max_place['lat']
min_lon_id, min_lat_id = min_place['lon'], min_place['lat']
SOI = slp[:, max_lat_id, max_lon_id] - slp[:, min_lat_id, min_lon_id]  # 计算南方涛动指数

nino34_norm = (nino34 - nino34.mean('time')) / nino34.std('time')
SOI_norm = (SOI - SOI.mean('time')) / SOI.std('time')

# 计算自相关系数
def calc_autoCorr(*, n, j, x):
    autoCorr = 0
    for t in range(n - j):
        autoCorr += (x[t] * x[t + j]) / (n - j)
    return float(autoCorr)

# 计算超前滞后相关系数函数
def calc_lagCorr(*, n, j, x, y):
    lagCorr = 0
    for t in range(n - j):
        lagCorr += (x[t] * y[t + j]) / (n - j)
    return float(lagCorr)

nino34_autoCorrs = []
SOI_autoCorrs = []
nino34_SOI_lagCorrs = []

for j in range(30):
    nino34_autoCorr = calc_autoCorr(n=30, j=j, x=nino34_norm)
    nino34_autoCorrs.append(calc_autoCorr(n=30, j=j, x=nino34_norm))

    SOI_autoCorr = calc_autoCorr(n=30, j=j, x=SOI_norm)
    SOI_autoCorrs.append(calc_autoCorr(n=30, j=j, x=SOI_norm))

    nino34_SOI_lagCorr = calc_lagCorr(n=30, j=j, x=nino34_norm, y=SOI_norm)
    nino34_SOI_lagCorrs.append(nino34_SOI_lagCorr)

def autoCorr_plot():
    fig = plt.figure(figsize=(6, 4), facecolor='w')
    ax = fig.add_subplot()
    ax.set_title('autoCorr & lagCorr', loc='left')
    ax.plot(range(30), nino34_autoCorrs, label='nino34 autoCorr', lw=1)
    ax.plot(range(30), SOI_autoCorrs, label='SOI autoCorr', lw=1)
    ax.plot(range(30), nino34_SOI_lagCorrs, label='nino34 & SOI lagCorr', lw=1)
    for i in ['bottom', 'top', 'left', 'right']:
        ax.spines【i】.set_linewidth(0.5)     # 记得改成英文的中括号,气象家园无法显示英文中括号加i
    ax.tick_params(which='major', width=0.5, length=5)
    ax.set_xlabel('Year Intervel')
    ax.set_ylabel('Corr')
    ax.set_xlim(-1, 30)
    ax.set_ylim(-0.7, 1.2)
    ax.legend()
    plt.savefig('autoCorr.png', dpi=300, bbox_inches='tight')
autoCorr_plot()


slp_ano.png
nino34_slp_corr.png
autoCorr.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2022-7-7 15:28:57 | 显示全部楼层
谢谢楼主提供脚本及图例
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表