- 积分
- 1428
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2021-7-21
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 xiaoxiao啊 于 2022-7-7 17:09 编辑
基于jupyter notebook
数据获取方式如下
链接:https://pan.baidu.com/s/17TSNCN-XXV5zrMXhIG9iTA
提取码:0un5
# 导入所需库
import numpy as np
import pandas as pd
import xarray as xr
from xskillscore import pearson_r, pearson_r_p_value
from scipy.interpolate import Rbf
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cnmaps import get_map, draw_map, clip_contours_by_map, clip_pcolormesh_by_map
import cmaps
import warnings
warnings.filterwarnings('ignore')
# 更改字体
plt.rcParams['font.sans-serif'] = ['Times New Roman']
plt.rcParams['axes.unicode_minus'] = False
# 数据处理,根据公式计算EU指数
f1 = xr.open_dataset('hgt.mon.mean.nc')
z = f1.hgt.loc[f1.time.dt.month.isin([1]), 500].loc['1948':'2012']
z1 = z.loc[:, 55, 20].drop(['lat', 'lon'])
z2 = z.loc[:, 55, 75].drop(['lat', 'lon'])
z3 = z.loc[:, 40, 145].drop(['lat', 'lon'])
eu = -z1 / 4 + z2 / 2 - z3 / 4
eu_norm = (eu - eu.mean('time')) / eu.std('time') # 标准化后的EU指数
# 绘画EU指数
def eu_plot():
fig = plt.figure(figsize=(6, 3), facecolor='w')
ax = fig.add_subplot(111)
ax.tick_params(which='major', width=0.5, length=5)
bars = ax.bar(range(1948, 2013), eu_norm, facecolor='#ff9999', edgecolor='r', width=1)
for bar, height in zip(bars, eu_norm):
if height < 0:
bar.set(facecolor='#9999ff', edgecolor='blue')
for i in ['bottom', 'top', 'left', 'right']:
ax.spines【i】.set_linewidth(0.5) # 记得改成英文的中括号,气象家园无法显示英文中括号加i
ax.set_xlim(1947, 2013)
ax.set_xlabel('year')
ax.set_title('EU Index(1948-2012)', loc='left')
eu_plot()
# 将一维的EU指数转化为三维的
def to_3dims(*, var):
var = np.array(var)
var = np.expand_dims(var, axis=1)
var = np.repeat(var, 73 * 144, axis=1)
var = np.reshape(var, [65, 73, 144])
var = xr.DataArray(var, coords=[z.time, z.lat, z.lon], dims=['time', 'lat', 'lon'])
return var
eu_norm_3dims = to_3dims(var=eu_norm)
# EU遥相关指数与同期环流场(500hPa高度场)的相关系数
r1 = pearson_r(eu_norm_3dims, z, dim='time')
p1 = pearson_r_p_value(eu_norm_3dims, z, dim='time')
mesh_lon, mesh_lat = np.meshgrid(f1.lon, f1.lat)
# 绘图
def z500_plot():
levels1 = np.arange(-1, 1.2, 0.2)
levels2 = np.arange(-0.8, 1, 0.2)
fig = plt.figure(figsize=(10, 6), facecolor='w')
ax = fig.add_subplot(111, projection=ccrs.PlateCarree(180))
ax.set_title('Corr EU & Z500 in Jan(1948-2012)', loc='left')
ax.add_feature(cfeature.COASTLINE.with_scale('110m'), linewidth=0.5, edgecolor='dimgrey')
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)
cf = ax.contourf(f1.lon, f1.lat, r1, cmap=cmaps.BlueDarkRed18, transform=ccrs.PlateCarree(), levels=levels1)
cb = fig.colorbar(cf, orientation='horizontal', shrink=0.6, aspect=25, pad=0.12, ticks=levels2)
cb.outline.set_linewidth(0.1)
cb.ax.tick_params(length=0)
dot_area = np.where(p1 < 0.05)
dot = ax.scatter(mesh_lon[dot_area], mesh_lat[dot_area], color='k', s=1.2, linewidths=0, transform=ccrs.PlateCarree())
ax.set_aspect(1.2)
z500_plot()
df = pd.read_csv('station.csv', encoding='gbk', names=['station', 'num', 'lat', 'lon'], header=0) # 读取160站站点数据
f2 = np.loadtxt('t1601.txt') # 160站气温数据
t = f2.reshape(63, 160)[:-1]
t = xr.DataArray(t, coords=[z.time[3:], df.station], dims=['time', 'station'])
# 将一维的EU指数转化为二维的
def to_2dims_ch(*, var):
var = np.array(var)
var = np.expand_dims(var, axis=1)
var = np.repeat(var, 160, axis=1)
var = np.reshape(var, [62, 160])
var = xr.DataArray(var, coords=[z.time[3:], df.station], dims=['time', 'station'])
return var
eu_norm_2dims = to_2dims_ch(var=eu_norm[3:])
# EU遥相关指数与同期我国气温的相关系数
r2 = pearson_r(eu_norm_2dims, t, dim='time')
p2 = pearson_r_p_value(eu_norm_2dims, t, dim='time')
lon = df.lon # 站点经纬度
lat = df.lat
grid_lon = np.linspace(73, 135, 120) # 格点经纬度
grid_lat = np.linspace(15, 55, 80)
meshgrid_lon, meshgrid_lat = np.meshgrid(grid_lon, grid_lat)
# 站点数据插值为格点数据,三次样条插值
func = Rbf(df.lon, df.lat, r2, function='cubic')
grid_r2 = func(meshgrid_lon, meshgrid_lat)
func = Rbf(df.lon, df.lat, p2, function='cubic')
grid_p2 = func(meshgrid_lon, meshgrid_lat)
# 绘图
def temp_plot():
levels1 = np.arange(-0.5, 0.6, 0.1)
levels2 = np.arange(-0.4, 0.5, 0.1)
fig = plt.figure(figsize=(6, 5.5), facecolor='w')
ax1 = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=ccrs.PlateCarree(105))
ax1.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5, edgecolor='dimgrey')
ax1.add_feature(cfeature.LAKES.with_scale('110m'), linewidth=0.5, edgecolor='dimgrey', facecolor='none')
ax1.set_extent([70, 140, 15, 55], ccrs.PlateCarree())
ax1.set_title('Corr EU & Temp in Jan(1951-2012)', loc='left')
ax1.set_xticks(np.arange(70, 145, 10), crs=ccrs.PlateCarree())
ax1.set_yticks(np.arange(15, 60, 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)
cf = ax1.contourf(grid_lon, grid_lat, grid_r2, cmap=cmaps.BlueDarkRed18, transform=ccrs.PlateCarree(), levels=levels1)
dot_area = np.where(grid_p2 < 0.05)
dot = ax1.scatter(meshgrid_lon[dot_area], meshgrid_lat[dot_area], color='k', s=1.2, linewidths=0, transform=ccrs.PlateCarree())
china = get_map('中国')
south_sea = get_map('南海')
draw_map(china, color='dimgrey', linewidth=0.5)
draw_map(south_sea, color='k', linewidth=0.6)
clip_contours_by_map(cf, china)
clip_pcolormesh_by_map(dot, china)
ax1.add_geometries(Reader('河流.shp').geometries(), crs=ccrs.PlateCarree(), edgecolor='dimgrey', facecolor='none', lw=0.5)
ax2 = fig.add_axes([0.768, 0.2005, 0.15, 0.15], projection=ccrs.PlateCarree())
ax2.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5, edgecolor='dimgrey')
ax2.set_extent([105, 124, 0, 23], ccrs.PlateCarree())
ax2.spines['geo'].set_linewidth(0.5)
draw_map(china, color='dimgrey', linewidth=0.5)
draw_map(south_sea, color='k', linewidth=0.6)
ax3 = fig.add_axes([0.15, 0.06, 0.7, 0.04])
cb = fig.colorbar(cf, cax=ax3, orientation='horizontal', ticks=levels2)
cb.outline.set_linewidth(0.1)
cb.ax.tick_params(length=0)
temp_plot()
|
-
EU(欧亚)遥相关指数
-
-
EU遥相关指数与同期我国气温的相关系数
|