Python Scipy二维核密度估计(KDE)出现异常高估孔洞问题排查

研究背景 #

• 研究目标：按风暴模态对美国本土（CONUS）的龙卷风事件开展核密度估计，复现Smith等人2012年公开论文中基于40km网格、按对流模态统计的逐年代龙卷风事件KDE图
• 所用数据集：2003-2011年龙卷风事件的经纬度散点记录
• 原方法参考：经咨询论文作者得知，原研究结果基于ArcGIS平台的Density函数生成，采用四次核（quartic kernel），可自动输出单位面积事件数

复现问题说明 #

• 现有实现与原方法的差异：目前仅找到高斯核（gaussian kernel）的实现方案，未匹配原方法使用的四次核；自行设计缩放系数完成单位换算
• 异常表现：首张KDE图件的复现效果基本符合预期，但后续图件出现异常高估问题：密度估计结果中部出现数值远高于合理水平的大范围空洞区域
• 待排查方向：由于对核密度估计相关统计方法熟悉度有限，无法确定异常诱因是带宽（bandwidth）设置错误，还是估计计算流程存在疏漏
• 当前缩放系数计算逻辑：单网格像元面积为1600 km²，乘以系数10/9，对应0.9个十年的时长换算

复现代码 #

# 加载40km RAP网格
f = np.load("/Users/andrewlyons/Downloads/pperf_grid_template.npz")
lon,lat = f['lon'], f['lat']
f.close()

proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(1, 1, 1, projection=proj)

# 设置绘图空间范围
ax.set_extent([-125,-61,22, 49])

states_provinces = cfeature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='50m',
        facecolor='none')

data = qlcs03

k = kde.gaussian_kde([data['slon'],data['slat']])

# 设置40km计算网格
xi,yi =lon,lat
zi = k(np.vstack([xi.flatten(), yi.flatten()]))

# 应用缩放系数
zi=(zi*(1600*(10/9)))

# 绘制结果
c =ax.contourf(xi, yi, zi.reshape(xi.shape),colors='k',levels=[0.5,1,2,3,4,5,6,7,8],alpha=0.17,transform=proj)
cs = ax.contour(xi, yi, zi.reshape(xi.shape),colors='k',levels=[0.5,1,2,3,4,5,6,7,8],linewidths=2.5,transform=proj,zorder=9)
ax.clabel(cs, fontsize=16, inline=False, colors ='r')
ax.scatter(data['slon'],data['slat'],color ='k',marker = '.',s=0.1)
ax.add_feature(cfeature.BORDERS)
ax.add_feature(states_provinces)
ax.add_feature(cfeature.LAND,facecolor ='white')
ax.add_feature(cfeature.COASTLINE,zorder=11)
ax.add_feature(cfeature.OCEAN,facecolor='white',zorder=10)

plt.title("Convective Mode KDE Test QLCS Tor 2003-2011 events per decade 40km grid                      N="+str(len(data)))
plt.show()

内容的提问来源于stack exchange，提问作者Twisterkid34