使用Cartopy绘制密度图时在不连续区域出现问题

问题背景 #

在Cartopy的NorthPolarStereo投影上绘制气旋密度差值时，数据缺失造成的不连续区域附近出现明显密度偏移，推测由网格转换过程中的变形导致。

生成的图像：

原实现代码：

import numpy as np
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt

lat1 = df1['xlat']
lon1 = df1['xlon']

lat2 = df2['xlat']
lon2 = df2['xlon']


bins1 = [np.linspace(lon1.min(), lon1.max(), 100), np.linspace(lat1.min(), lat1.max(), 100),]
hist1, xedges1, yedges1 = np.histogram2d(lon1, lat1, bins=bins1)

bins2 = [np.linspace(lon2.min(), lon2.max(), 100), np.linspace(lat2.min(), lat2.max(), 100),]
hist2, xedges2, yedges2 = np.histogram2d(lon2, lat2, bins=bins2)
    
hist_smooth1 = gaussian_filter(hist1, sigma=2)
hist_smooth2 = gaussian_filter(hist2, sigma=2)

fig, axs_total_difference = plt.subplots(1, 1, subplot_kw={'projection': ccrs.NorthPolarStereo(central_longitude=0.0, globe=None)}, figsize=(10, 7))

# Plot - Total Difference
total_difference = hist_smooth1 - hist_smooth2
density_total_difference = axs_total_difference.contourf(xedges1[:-1], yedges1[:-1], total_difference.T, transform=ccrs.PlateCarree(), cmap='bwr')
axs_total_difference.set_title('Total Difference (ERA5 - ICON)')

axs_total_difference.add_feature(cartopy.feature.OCEAN, color='white', zorder=0)
axs_total_difference.add_feature(cartopy.feature.LAND, color='lightgray', zorder=0, linewidth=0.5, edgecolor='black')
axs_total_difference.gridlines(draw_labels=True, linewidth=0.5, color='gray', xlocs=range(-180, 180, 15), ylocs=range(-90, 90, 15), x_inline=False, y_inline=False)
axs_total_difference.coastlines(resolution='50m', linewidth=0.3, color='black')

# Colorbar for total difference
cbar_ax_total = plt.gcf().add_axes([0.92, 0.15, 0.02, 0.7])  
cbar_total = plt.colorbar(density_total_difference, cax=cbar_ax_total)
cbar_total.set_label('Total Difference')

# Save or show the plot
plt.savefig('Cyclone Track Total Difference', dpi='figure', format='png')
plt.show()
plt.close()

问题根源 #

1. 网格投影顺序错误：先在经纬度（PlateCarree）网格统计密度，再投影到极射投影，经纬网格到极射投影的非线性变形会导致边界区域的密度单元位置偏移，叠加高斯平滑后放大了偏移效果。
2. 网格不统一：两个数据集使用各自的经纬度范围生成bins，导致直方图网格未对齐，差值计算本身存在系统误差。
3. 无掩码平滑：高斯平滑未区分有效数据和缺失区域，将有效区域的数值扩散到无数据的边界区域，加剧偏移。

解决方案 #

改进思路 #

• 先将所有气旋坐标转换到NorthPolarStereo投影的笛卡尔坐标系，再在该投影网格上统计密度，避免后续投影变形。
• 使用统一的全局网格覆盖北极区域，确保两个数据集的直方图完全对齐。
• 添加掩码处理，禁止高斯平滑向无数据区域扩散数值。

修改后的代码 #

import numpy as np
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# 定义极射投影与经纬度投影
stereo_proj = ccrs.NorthPolarStereo(central_longitude=0.0)
plate_carree = ccrs.PlateCarree()

# 将经纬度坐标转换为极射投影的笛卡尔坐标
x1, y1 = stereo_proj.transform_points(plate_carree, lon1.values, lat1.values).T
x2, y2 = stereo_proj.transform_points(plate_carree, lon2.values, lat2.values).T

# 过滤转换后的无效坐标
valid_mask1 = ~(np.isnan(x1) | np.isnan(y1))
x1, y1 = x1[valid_mask1], y1[valid_mask1]
valid_mask2 = ~(np.isnan(x2) | np.isnan(y2))
x2, y2 = x2[valid_mask2], y2[valid_mask2]

# 设置覆盖北极区域的统一网格（单位：米）
x_bounds = [-5e6, 5e6]
y_bounds = [-5e6, 5e6]
bins = [np.linspace(*x_bounds, 100), np.linspace(*y_bounds, 100)]

# 在极射投影网格上计算直方图
hist1, xedges, yedges = np.histogram2d(x1, y1, bins=bins)
hist2, _, _ = np.histogram2d(x2, y2, bins=bins)

# 创建数据掩码，标记无数据的网格单元
mask1 = hist1 == 0
mask2 = hist2 == 0

# 带掩码的高斯平滑，避免无数据区域被扩散填充
hist_smooth1 = gaussian_filter(hist1.astype(float), sigma=2)
hist_smooth1[mask1] = 0

hist_smooth2 = gaussian_filter(hist2.astype(float), sigma=2)
hist_smooth2[mask2] = 0

# 绘图
fig, ax = plt.subplots(1, 1, subplot_kw={'projection': stereo_proj}, figsize=(10, 7))

total_difference = hist_smooth1 - hist_smooth2
# 直接绘制极射投影网格数据，无需额外transform
density_diff = ax.contourf(xedges[:-1], yedges[:-1], total_difference.T, cmap='bwr')

ax.set_title('Total Difference (ERA5 - ICON)')
ax.add_feature(cfeature.OCEAN, color='white', zorder=0)
ax.add_feature(cfeature.LAND, color='lightgray', zorder=0, linewidth=0.5, edgecolor='black')
ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', xlocs=range(-180, 180, 15), ylocs=range(60, 90, 5), x_inline=False, y_inline=False)
ax.coastlines(resolution='50m', linewidth=0.3, color='black')

# 添加色条
cbar_ax = plt.gcf().add_axes([0.92, 0.15, 0.02, 0.7])
cbar = plt.colorbar(density_diff, cax=cbar_ax)
cbar.set_label('Total Difference')

plt.savefig('Cyclone Track Total Difference.png', dpi=300, bbox_inches='tight')
plt.show()
plt.close()

关键改进点说明 #

• 坐标先转换：提前将经纬度转为极射投影坐标，密度统计直接在目标投影网格上进行，彻底避免投影变形带来的偏移。
• 统一网格：使用固定的北极区域网格范围，确保两个数据集的直方图单元完全对齐，消除差值计算的系统误差。
• 掩码平滑：通过掩码限制高斯平滑的作用范围，防止无数据区域影响边界的密度数值，减少偏移现象。

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