插值后的Healpy坐标误差:bis的出现

2024-05-29 07:39:26 发布

您现在位置:Python中文网/ 问答频道 /正文

我有一个由128个点组成的粗略的天空地图,我想制作一个平滑的healpix地图(见附图,左手侧)。文中引用的数字:

here

我加载我的数据,然后为最终的地图制作一个新的经度和纬度数组(例如nside=32)。在

我的输入数据是:

lats = pi/2 + ths   # theta from 0, pi, size 8
lons = phs          # phi from 0, 2pi, size 16
data = sky_data[0]  # shape (8,16)

新的lon/lat数组大小基于nside的像素数:

^{pr2}$

然后通过插值找到这些像素的新数据值,然后从角度转换回像素。在

# new lon/lat
new_lats = hp.pix2ang(nside, pixIdx)[0] # thetas I need to populate with interpolated theta values
new_lons = hp.pix2ang(nside, pixIdx)[1] # phis, same

# interpolation
lut = RectSphereBivariateSpline(lats, lons, data, pole_values=4e-14)
data_interp = lut.ev(new_lats.ravel(), new_lons.ravel()) #interpolate the data
pix = hp.ang2pix(nside, new_lats, new_lons) # convert latitudes and longitudes back to pixels

然后,我用插值值构造healpy映射:

healpix_map = np.zeros(hp.nside2npix(nside), dtype=np.double) # create empty map
healpix_map[pix] = data_interp # assign pixels to new interpolated values
testmap = hp.mollview(healpix_map)

地图的结果是附图的右上方。在

(请原谅使用jet——viridis没有“白色”零,因此使用颜色贴图会添加蓝色背景。)

地图看起来不对:从图中的粗略地图可以看出,右下方应该有一个“热点”,但这里它出现在左上方。在

作为一个健全性检查,我使用matplotlib绘制了一个mollview投影中插值点的散点图,如图2所示,我去掉了标记的边缘,使其看起来像一张地图;)

ax = plt.subplot(111, projection='astro mollweide')
ax.grid()
colors = data_interp
sky=plt.scatter(new_lons, new_lats-pi/2, c = colors, edgecolors='none', cmap ='jet')
plt.colorbar(sky, orientation = 'horizontal')

你可以看到,这张地图,在附图的右下方,产生的正是我所期望的!所以坐标没问题,我完全搞不懂。在

以前有人遇到过这种情况吗?我能做什么?我想在这个和未来的地图上使用healpy函数,所以仅仅使用matplotlib不是一个选择。在

谢谢!在


Tags: to数据mapnewdata地图pihp
2条回答

我发现我必须将pi/2添加到θ中才能使插值生效,因此最终需要应用以下变换才能正确渲染图像:

newnew_lats = pi - new_lats
newnew_lons = pi + new_lons

似乎仍然有一点关于插值的问题,虽然似乎不是那么明显了。我可以试试另一个来比较。在

我不是healpix方面的专家(事实上,我以前从未使用过它——我是一个粒子物理学家),但据我所知,这只是一个惯例问题:在莫尔韦德投影中,healpy将北极(正纬度)放在地图的底部,出于某种原因。我不确定它为什么会这样做,或者这是否是故意的行为,但似乎很清楚这就是发生的事情。如果我遮住赤道以下的一切,即只保留正纬度点

mask = new_lats - pi/2 > 0
pix = hp.ang2pix(nside, new_lats[mask], new_lons[mask])
healpix_map = np.zeros(hp.nside2npix(nside), dtype=np.double)
healpix_map[pix] = data_interp[mask]
testmap = hp.mollview(healpix_map)

它得出一个中心线以上没有数据的绘图:

masked plot from healpy

至少它很容易修复。mollview允许一个rot参数,该参数将在投影之前有效地围绕观察轴旋转球体,而flip参数可以设置为'astro'(默认)或{}来设置东方是显示在左侧还是右侧。一个小实验表明,你得到了你想要的坐标系

^{pr2}$

在元组中,前两个元素是要设置在绘图中心的点的经度和纬度,第三个元素是旋转。所有的单位都是度。没有面具的话,它会给你:

unmasked rotated image

我相信这正是你想要的。在

相关问题 更多 >

    热门问题