模式与观测对比的散点图

功能:

1. 分析模式模拟数据时，常用散点图来表征模式模拟数据和观测数据的一致性，散点图分布越接近参考线y=x，模拟值越接近观测值

程序清单：

# 导入包:

import xarray as xr

import numpy as np

import matplotlib.pyplot as plt

#读入模式数据

data=xr.open_dataset('F:/python/data/6.24.7_conc_DJF.nc')

lon = data['lon']

lat = data['lat']

conc_DJF = data['conc_DJF']

#读入观测数据

data_obs = np.load('F:/python/data/pm25_15_winter.npy')

lat_obs = data_obs[:,0]

lon_obs = data_obs[:,1]

pm25 = data_obs[:,2]

#去除缺测的站点值

obs = xr.DataArray(data_obs,dims=['x','y'])

obs = obs.dropna(dim='x')

#将模式的网格点数据根据站点位置进行插值

da = xr.DataArray(conc_DJF,[('lat', lat.data),('lon', lon.data)])

mod = []

for i in range(len(obs)):

    latID = np.argmin(abs(obs[i,0]-lat).data)  #取绝对值最近的值

    lonID = np.argmin(abs(obs[i,1]-lon).data)

    mod.append(conc_DJF[latID,lonID])

#计算归一化平均偏差NMB值

a = mod - obs[:,2]

NMB = np.sum(mod - obs[:,2])/np.sum(obs[:,2])

print("%.2f%%" %(NMB*100))

#计算相关系数

from scipy.stats import pearsonr

R,P = pearsonr(mod,obs[:,2])

print(R,P)

#创建画布，画布大小8*8，分辨率300

fig = plt.figure(figsize=(8, 8),dpi=300)

ax=fig.add_subplot(1,1,1)

#设置刻度标签朝内

plt.rcParams['xtick.direction'] = 'in'  

plt.rcParams['ytick.direction'] = 'in'

#绘制散点，obs为观测数据，mod为模式数据，facecolor为填充颜色，edgecolor为边框颜色，linewidth为线宽，marker为散点形状

ax.scatter(obs[:,2], mod, s=20, facecolor='none', edgecolor="forestgreen",linewidths=1.5, marker='^')   

#采用对数坐标系，semilogx或semilogy为半对数坐标系

ax.loglog()   

x = np.linspace(0.1,1000)

#设置x轴和y轴等长

plt.axis('scaled')  

#绘制y=x，y=3x以及y=1/3x参考线，zorder设置显示顺序

ax.plot(x,x,ls="-",c='k',lw=1.5,zorder=2)   

ax.plot(x,3*x,ls="--",c='k',lw=1.5,zorder=2)

ax.plot(x,1/3*x,ls="--",c='k',lw=1.5,zorder=2)

#对数坐标轴下必须>0，<或=均不可

ax.set(xlim=(4,4e2), ylim=(4,4e2))

#设置横纵坐标刻度     

ax.set_xticks([5,10,25,50,100,200,400])

ax.set_xticklabels([5,10,25,50,100,200,400],fontsize=18)

ax.set_yticks([5,10,25,50,100,200,400])

ax.set_yticklabels([5,10,25,50,100,200,400],fontsize=18)

#设置上、左、右、下边框宽度

for spine in ['top','left','right','bottom']:   

    ax.spines[spine].set_linewidth(1.5)

#在图上x=6，y=280的位置显示归一化平均偏差NMB值和相关系数

ax.text(6, 280, 'NMB = '+"%.2f%%" %(NMB*100), fontsize =18)    

ax.text(6, 230, 'R = '+str('%.2f' % R), fontsize =18)

ax.set_xlabel('Observed aerosol conc.',fontsize=20)

ax.set_ylabel('Modeled aerosol conc.',fontsize=20)

plt.show()