使用 healpy 使用 HEALPix 像素化制作二维直方图
Make a 2D histogram with HEALPix pixellization using healpy
数据为天体坐标,举例如下:
import pylab as plt
import numpy as np
l = np.random.uniform(-180, 180, 2000)
b = np.random.uniform(-90, 90, 2000)
我想做一个二维直方图,以便在天空中使用 (l, b)
坐标绘制某个点的密度图,在 Mollweide 投影上使用 HEALPix 像素化。我怎样才能使用 healpy 做到这一点?
教程:
说明了如何绘制一维数组或拟合文件,但我找不到如何使用这种像素化来绘制二维直方图。
我也找到了这个功能,但是没有用,所以卡住了。
hp.projaxes.MollweideAxes.hist2d(l, b, bins=10)
我可以用这种方式在 Mollweide 投影中绘制这些点的图:
l_axis_name ='Latitude l (deg)'
b_axis_name = 'Longitude b (deg)'
fig = plt.figure(figsize=(12,9))
ax = fig.add_subplot(111, projection="mollweide")
ax.grid(True)
ax.scatter(np.array(l)*np.pi/180., np.array(b)*np.pi/180.)
plt.show()
非常感谢您的帮助。
好问题!我编写了一个简短的函数来将目录转换为数字计数的 HEALPix 图:
from astropy.coordinates import SkyCoord
import healpy as hp
import numpy as np
def cat2hpx(lon, lat, nside, radec=True):
"""
Convert a catalogue to a HEALPix map of number counts per resolution
element.
Parameters
----------
lon, lat : (ndarray, ndarray)
Coordinates of the sources in degree. If radec=True, assume input is in the icrs
coordinate system. Otherwise assume input is glon, glat
nside : int
HEALPix nside of the target map
radec : bool
Switch between R.A./Dec and glon/glat as input coordinate system.
Return
------
hpx_map : ndarray
HEALPix map of the catalogue number counts in Galactic coordinates
"""
npix = hp.nside2npix(nside)
if radec:
eq = SkyCoord(lon, lat, 'icrs', unit='deg')
l, b = eq.galactic.l.value, eq.galactic.b.value
else:
l, b = lon, lat
# conver to theta, phi
theta = np.radians(90. - b)
phi = np.radians(l)
# convert to HEALPix indices
indices = hp.ang2pix(nside, theta, phi)
idx, counts = np.unique(indices, return_counts=True)
# fill the fullsky map
hpx_map = np.zeros(npix, dtype=int)
hpx_map[idx] = counts
return hpx_map
然后您可以使用它来填充 HEALPix 地图:
l = np.random.uniform(-180, 180, 20000)
b = np.random.uniform(-90, 90, 20000)
hpx_map = hpx.cat2hpx(l, b, nside=32, radec=False)
此处,nside
决定像素网格的精细程度。
hp.mollview(np.log10(hpx_map+1))
另请注意,通过在银河纬度上均匀采样,您会更喜欢银河两极的数据点。如果你想避免这种情况,你可以用余弦缩小它。
hp.orthview(np.log10(hpx_map+1), rot=[0, 90])
hp.graticule(color='white')
数据为天体坐标,举例如下:
import pylab as plt
import numpy as np
l = np.random.uniform(-180, 180, 2000)
b = np.random.uniform(-90, 90, 2000)
我想做一个二维直方图,以便在天空中使用 (l, b)
坐标绘制某个点的密度图,在 Mollweide 投影上使用 HEALPix 像素化。我怎样才能使用 healpy 做到这一点?
教程:
说明了如何绘制一维数组或拟合文件,但我找不到如何使用这种像素化来绘制二维直方图。
我也找到了这个功能,但是没有用,所以卡住了。
hp.projaxes.MollweideAxes.hist2d(l, b, bins=10)
我可以用这种方式在 Mollweide 投影中绘制这些点的图:
l_axis_name ='Latitude l (deg)'
b_axis_name = 'Longitude b (deg)'
fig = plt.figure(figsize=(12,9))
ax = fig.add_subplot(111, projection="mollweide")
ax.grid(True)
ax.scatter(np.array(l)*np.pi/180., np.array(b)*np.pi/180.)
plt.show()
非常感谢您的帮助。
好问题!我编写了一个简短的函数来将目录转换为数字计数的 HEALPix 图:
from astropy.coordinates import SkyCoord
import healpy as hp
import numpy as np
def cat2hpx(lon, lat, nside, radec=True):
"""
Convert a catalogue to a HEALPix map of number counts per resolution
element.
Parameters
----------
lon, lat : (ndarray, ndarray)
Coordinates of the sources in degree. If radec=True, assume input is in the icrs
coordinate system. Otherwise assume input is glon, glat
nside : int
HEALPix nside of the target map
radec : bool
Switch between R.A./Dec and glon/glat as input coordinate system.
Return
------
hpx_map : ndarray
HEALPix map of the catalogue number counts in Galactic coordinates
"""
npix = hp.nside2npix(nside)
if radec:
eq = SkyCoord(lon, lat, 'icrs', unit='deg')
l, b = eq.galactic.l.value, eq.galactic.b.value
else:
l, b = lon, lat
# conver to theta, phi
theta = np.radians(90. - b)
phi = np.radians(l)
# convert to HEALPix indices
indices = hp.ang2pix(nside, theta, phi)
idx, counts = np.unique(indices, return_counts=True)
# fill the fullsky map
hpx_map = np.zeros(npix, dtype=int)
hpx_map[idx] = counts
return hpx_map
然后您可以使用它来填充 HEALPix 地图:
l = np.random.uniform(-180, 180, 20000)
b = np.random.uniform(-90, 90, 20000)
hpx_map = hpx.cat2hpx(l, b, nside=32, radec=False)
此处,nside
决定像素网格的精细程度。
hp.mollview(np.log10(hpx_map+1))
另请注意,通过在银河纬度上均匀采样,您会更喜欢银河两极的数据点。如果你想避免这种情况,你可以用余弦缩小它。
hp.orthview(np.log10(hpx_map+1), rot=[0, 90])
hp.graticule(color='white')