- 积分
- 1428
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2021-7-21
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
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()
|
|