一、要求
计算EU（欧亚）遥相关指数，输出1月份该指数年际变化的时间序列；
计算EU遥相关指数与同期环流场（500hPa高度场或海平面气压场）的相关系数；
计算EU遥相关指数与同期我国气温的相关系数。
二、资料说明
（1）格点资料
NCEP/NCAR 1948－2012年（65年）的500百帕月平均高度场资料
资料范围为（900S-900N,00-3600E)
网格距为2.50×2.50,纬向格点数为144，经向格点数为73
资料为GRD格式，资料从南到北、自西向东排列，每月为一个记录，按年逐月排放。
（2）站点资料
我国气候中心整编的160站月平均气温资料；
全国160个台站；
所给的资料是1月份的（1951.1-2013.6）；
资料为txt格式，详见资料说明。
三、方法说明
EU遥相关指数的定义:
相关系数计算：
四、步骤说明
- 编程计算1月份EU遥相关指数（标准化），画出指数图；
- 编程计算1月份EU遥相关指数与500hPa高度场的相关系数分布图，并绘制图形；
- 编程计算EU遥相关指数与1月份我国气温的相关系数分布图，并绘制图形。
五、代码
本次代码分为两个程序。
（1）EU遥相关指数（标准化）与EU遥相关指数与500hPa高度场的相关系数
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
import scipy.stats as stats
from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 绘图
def linedraw(y,syear,endyear,title):
fig, ax1 = plt.subplots(figsize=(8, 5))
x = np.arange(syear, endyear+1)
ax1.axis('auto') # 参数 'scale','auto','square'
ax1.plot(x, y, c='navy', linewidth=2, linestyle='-', marker='s', markerfacecolor='w',
markeredgecolor='navy', label='EU', alpha=0.7)
ax1.set_xlabel('年份', fontsize=12)
ax1.set_ylabel('欧亚遥相关指数', fontsize=12)
xticks = np.arange(syear, endyear+1, 4)
ax1.set_xticks(xticks)
ax1.set_ylim(-3, 2)
ax1.legend(loc='upper left', edgecolor='black',fontsize=12)
ax1.set_title(title, fontsize=14)
ax1.spines['bottom'].set_linewidth(1.5)
ax1.spines['left'].set_linewidth(1.5)
ax1.spines['top'].set_linewidth(1.5)
ax1.spines['right'].set_linewidth(1.5)
ax1.grid(axis='both', ls='--',linewidth=1.5) # axis = 'both','x','y'
plt.show()
def counterdraw(z, parea, lon, lat, resolution):
# 绘制地图底图
fig = plt.figure(figsize=(10, 6))
proj = ccrs.PlateCarree(central_longitude=180)
leftlon, rightlon, lowerlat, upperlat = (0, 360, -90, 90)
img_extent = [leftlon, rightlon, lowerlat, upperlat]
m = np.min(z)
n = np.max(z)
ax = fig.add_axes([0.06, 0., 0.98, 0.98], projection=proj)
c = ax.contourf(lon, lat, z, levels=np.arange(-1, 1.4, 0.4),transform=ccrs.PlateCarree(), cmap='bwr')
c1 = ax.contour(lon, lat, z, levels=np.arange(-1, 1, 0.2),transform=ccrs.PlateCarree(), colors='k')
nx = -np.ones(parea[0, :].shape) * 90 + parea[0, :] * resolution
ny = parea[1, :] * resolution
sig1 = ax.scatter(ny, nx, marker='.', s=1, c='k', alpha=0.9, transform=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.add_feature(cfeature.LAKES)
ax.set_extent(img_extent, crs=ccrs.PlateCarree())
ax.set_xticks(np.arange(leftlon, rightlon + 20, 20), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(lowerlat, upperlat + 20, 20), crs=ccrs.PlateCarree())
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
#ax.set_title(name, loc='center', fontsize=16)
xlocs = np.arange(leftlon, rightlon+20, 20)
ylocs = np.arange(lowerlat, upperlat, 20)
plt.colorbar(c, shrink = 0.4, pad = 0.06)
plt.clabel(c1, fontsize=8)
#plt.savefig(r'C:\Users\HUAWEI\Desktop\700hPa露点温度差20时.png')
plt.show()
# 数据读取
f1 = np.fromfile(r'hgt500.grd',dtype=np.float32) #读取高度场数据
f1 = np.array(f1)
hgt = np.reshape(f1,(780,73,144))
lon = np.arange(0,360,2.5)
lat = np.arange(-90,92.5,2.5)
janhgt = hgt[0:780:12,:,:] # 只取1月份的数据
#print(janhgt.shape)
# 计算EU遥相关指数
EUs = np.zeros(65) # 每年一个结果共65年
for i in range(65):
EUs[i] = -(1./4.) * janhgt[i,58,8] + (1./2.) * janhgt[i,58,30] - (1./4.) * janhgt[i,52,58]
std = np.std(EUs)
for i in range(65):
EUs[i] = (EUs[i] - np.mean(EUs)) / std
#print(EUs.shape)
# 计算相关系数
Corr = np.zeros((73,144)) # 高度场为格点数据，有经度，纬度，时间，三个维度，经纬度确定一个空间位置，每个空间位置上有长度为65（时间数）数据向量，用此向量与EU指数构成的向量求取相关系数。
pvalue = np.zeros((73,144))
for i in range(73):
for j in range(144):
Corr[i,j],pvalue[i,j] = stats.pearsonr(EUs[:],janhgt[:,i,j]) # 相关系数矩阵
# t检验 pvalue小于显著性水平（0.01）时，拒绝原假设，认为两个变量之间具有显著相关性。
#print(np.min(pvalue))
area = np.array(np.where(pvalue<0.01)) # 获得显著区域的格点坐标，并存储在area中
linedraw(EUs,1948,2012,'1948-2012年1月份欧亚遥相关指数') # 绘制EU指数折线图
counterdraw(Corr, area, lon, lat, 2.5) # 绘制相关系数分布图
（2）EU遥相关指数与1月份我国气温的相关系数
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
import scipy.stats as stats
from cartopy.io.shapereader import Reader
import xarray as xr
import geopandas as gpd
import salem
# 绘图
def counterdraw(z, lon, lat, nx, ny, leftlon, rightlon, lowerlat, upperlat, map):
# 绘制地图底图
fig = plt.figure(figsize=(10, 6))
proj = ccrs.PlateCarree(central_longitude=(leftlon+rightlon)/2.)
img_extent = [leftlon-3, rightlon-1, lowerlat-2.5, upperlat-1]
ax = fig.add_axes([0.1, 0.01, 0.9, 1], projection=proj)
c = ax.contourf(lon, lat, z, levels=np.arange(-0.5, 0.6, 0.1),transform=ccrs.PlateCarree(), cmap='bwr')
#c = ax.contour(lon, lat, z,levels=np.arange(-1, 1, 0.2),transform=ccrs.PlateCarree(), cmap='bwr')
sig1 = ax.scatter(nx, ny, marker='.', s=1, c='k', alpha=0.9, transform=ccrs.PlateCarree())
ax.add_geometries(Reader(map).geometries(),crs=ccrs.PlateCarree(),
facecolor='none', edgecolor='k', linewidth=1)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.add_feature(cfeature.LAKES)
ax.set_extent(img_extent, crs=ccrs.PlateCarree())
ax.set_xticks(np.arange(leftlon, rightlon, 5), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(lowerlat, upperlat, 5), crs=ccrs.PlateCarree())
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
#ax.set_title(name, loc='center', fontsize=16)
plt.colorbar(c, shrink = 0.4, pad = 0.06)
#plt.clabel(c, fontsize=8)
#plt.savefig(r'C:\Users\HUAWEI\Desktop\700hPa露点温度差20时.png')
plt.show()
# 温度数据读取
filename=r'E:\CL03\sx3\data\t1601.txt' # 存储温度站点数据的地址
t = pd.read_csv(filename,header=None,sep='\s+')
jant = np.array(t)/10.
jant = np.reshape(jant,(63,160))
#print(jant[:,120])
#averjant = jant.mean((0)) # 按0（行）轴求平均（列方向），得到列的平均
#print(averjant.shape)
# 高度数据读取
f1 = np.fromfile(r'hgt500.grd',dtype=np.float32)
f1 = np.array(f1)
hgt = np.reshape(f1,(780,73,144))
janhgt = hgt[0:780:12,:,:]
# 站点与经纬度读取
filename2 = r'E:\CL03\sx3\data\lat_lon.txt' # 站点经纬度数据地址
f2 = pd.read_csv(filename2,header=None,sep='\s+')
locs = np.array(f2)
# 计算EU遥相关指数
EUs = np.zeros(65)
for i in range(65):
EUs[i] = -(1./4.) * janhgt[i,58,8] + (1./2.) * janhgt[i,58,30] - (1./4.) * janhgt[i,52,58]
std = np.std(EUs)
for i in range(65):
EUs[i] = (EUs[i] - np.mean(EUs)) / std
Corr = np.zeros((160)) # 相关系数
pvalue = np.zeros((160)) # p值，用于判断显著性，后续也会用黑点绘制在图中，显示显著区域。
for i in range(160):
Corr[i],pvalue[i] = stats.pearsonr(EUs[3:],jant[:-1,i])
# 将站点数据插值到格点上，方便绘图，Python的填色图绘制都是以格点数据为基础的。
from scipy.interpolate import griddata
x = locs[:,2] # 经度
y = locs[:,1] # 纬度
#z = jant[0,:] # 观测值
z = Corr
x_min, x_max = int(min(x)), int(max(x))+5 # 设置插值的经度范围
y_min, y_max = int(min(y)), int(max(y))+6 # 设置插值的纬度范围
resolution = 1 # 空间分辨率
xi = np.arange(x_min, x_max, resolution)
yi = np.arange(y_min, y_max, resolution)
xi, yi = np.meshgrid(xi, yi) # 格点化
# 使用克里金插值
Corr_i = griddata((x, y), z, (xi, yi), method='nearest') # method='cubic'三次样条插值方法，method='nearest','linear'等， 这里'nearest'效果最好
pvalue_i = griddata((x, y), pvalue, (xi, yi), method='nearest')
# zi纬度为横坐标，经度为纵坐标
# 构造Dataset
lon = np.arange(x_min, x_max, resolution)
lat = np.arange(y_min, y_max, resolution)
data = xr.DataArray(Corr_i, coords=[lat, lon], dims=["latitude","longitude"])
ds = data.to_dataset(name='corrs')
corrs = ds.corrs
datap = xr.DataArray(pvalue_i, coords=[lat, lon], dims=["latitude","longitude"])
dsp = datap.to_dataset(name='pvalues')
pvalues = dsp.pvalues
# 绘图
shp = r"E:\map\china.shp" # 存储中国地图.shp文件等的地址
# 创建掩膜数据，即将超出中国以外的地区掩盖，不绘制出来。
tp_shp = gpd.read_file(shp)
corrs_tp = corrs.salem.roi(shape=tp_shp)
pvalues_tp = pvalues.salem.roi(shape=tp_shp)
# 寻找显著区位置
#print(pvalue_i.shape)
#area = np.array(np.where(pvalue_i<0.5))
#print(area)
area1 = np.array(np.where(pvalues_tp<0.05)) # 当p值小于0.05时，认为是显著的，并以此为依据确定显著区域的个点坐标，存储在area1中。
#print(area1)
# 显著区域格点坐标转换为经纬度坐标，用以绘图
ny = np.ones(area1[0, :].shape) * y_min + area1[0, :] * resolution
nx = x_min*np.ones(area1[1, :].shape) + area1[1, :] * resolution
#print(nx,ny)
counterdraw(corrs_tp.values,lon,lat,nx,ny,x_min,x_max,y_min,y_max,shp) # 绘图
六、结果
图 1 1948-2021年1月份EU遥相关指数
图 2 1月份EU遥相关指数与500hPa高度场相关系数分布图
图 3 EU遥相关指数与1月份我国气温相关系数分布图
七、结语
从图1可以看到EU指数随时间的变化，图2、图3可以看到EU指数与环流形势、我国气温的相关系数分布，显著区域也用黑点显示出来，据此可以很方便的进行进一步分析讨论。
由于本人水平限制，多有错误，欢迎批评指正！
扫一扫
752
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
