如何在Healpy地图上叠加实测与建模数据的等高线？

我手头有两类数据：

• 实测数据：范艾伦带、南大西洋异常等地球捕获带电粒子数据
• 建模数据：AE8/AP8或AE9/AP9模型输出数据

目标是将两类数据绘制在同一张Healpy地图中，在单图内对比实测与建模结果：将实测数据绘制成黑白底图，叠加建模数据的等高线（区分质子和电子）。

最终成功实现方案 #

我找到了Healpix专用绘图函数库中的等高线绘制函数，使用效果良好。按照处理实测数据生成log_mymap的流程，先将建模数据处理为对数地图log_mymap_mod，再通过以下代码实现叠加：

MyMap = hp.newvisufunc.projview(log_mymap,                                
                                title='map_and_contour', 
                                coord=["G"], 
                                graticule=True,
                                graticule_labels=True,
                                flip='geo',
                                cmap=cmap_custom,
                                projection_type = 'hammer',
                                phi_convention="symmetrical",
                                cbar=False,
                                hold=True,
                                alpha=0.3
                               )

healpix_contour(log_mymap_mod, levels=3, colors='firebrick', alpha=1.0, linewidths=1.0)

plt.show()

最初尝试（未成功） #

实测数据包含经度（0-360°）、纬度和每秒计数，我先将经度转换为-180°到180°范围，再生成对数地图。建模数据（标记为_mod）的经度同样是0-360°，采用相同的经度转换方式尝试叠加等高线，但未成功，代码如下：

lonAll2 = np.zeros_like(lonAll)

for i in range(len(lonAll)):
    if lonAll[i] >= 180.0:
        lonAll2[i] = lonAll[i] - 360.0
    else:
        lonAll2[i] = lonAll[i]

theta = (90.0 - latAll) * (np.pi / 180.0) 
phi = lonAll2 * (np.pi / 180.0)  

k = 4
nside = 2 ** k

npix = hp.nside2npix(nside)
healpix_sum = np.zeros(npix)
healpix_count = np.zeros(npix)

for i in range(len(lonAll)):
    lon = lonAll2[i]
    lat = latAll[i]
    pixel_idx = hp.ang2pix(nside, np.radians(90 - lat), np.radians(lon))
    count_rate = rateAll[i]
    
    healpix_sum[pixel_idx] += count_rate
    if count_rate > 0:
        healpix_count[pixel_idx] += 1

mask = healpix_count > 0
healpix_mean = np.where(mask, healpix_sum / healpix_count, 0)

log_mymap = np.log10(healpix_mean)

data_mod=['a.txt']
    
lonAll_mod = np.array([])
latAll_mod = np.array([])
rateAll_mod = np.array([])

for filename_mod in data_mod:
    if exists(filename_mod):
        df = pd.read_table(filename_mod, sep=',', skiprows=14)
        
        lonAll_mod = np.concatenate((lonAll_mod, df['lon(deg)'].values))
        latAll_mod = np.concatenate((latAll_mod, df['lat(deg)'].values))
        rateAll_mod = np.concatenate((rateAll_mod, df['flux'].values))

    else:
        print("File not found:", filename_mod)

cmap_custom2 = plt.cm.get_cmap('Greys')

MyMap = hp.newvisufunc.projview(log_mymap,
                                title='VZLU_rch2_k4_mean_geq0_contour', 
                                coord=["G"], 
                                graticule=True,
                                graticule_labels=True,
                                flip='geo',
                                cmap=cmap_custom2,
                                projection_type = 'hammer',
                                phi_convention="symmetrical",
                                cbar=False
                               )

ax = plt.gca()

lonAll_mod2 = np.zeros_like(lonAll_mod)
for i in range(len(lonAll_mod)):
    if lonAll_mod[i] >= 180.0:
        lonAll_mod2[i] = lonAll_mod[i] - 360.0
    else:
        lonAll_mod2[i] = lonAll_mod[i]

lon_grid, lat_grid = np.meshgrid(np.linspace(-180, 180, int(np.sqrt(rateAll_mod.size))), 
                                 np.linspace(-90, 90, int(np.sqrt(rateAll_mod.size))))

rateAll_mod_grid = griddata((lonAll_mod2, latAll_mod), rateAll_mod, (lon_grid, lat_grid), method='cubic')

theta_grid, phi_grid = hp.pix2ang(nside=nside, ipix=hp.ang2pix(nside, lon_grid, lat_grid, lonlat=True), lonlat=True)

contour = ax.contour(theta_grid, phi_grid, rateAll_mod_grid, levels=5, colors='red')

plt.show()

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