Хотелось бы сделать карту плотности неба в формате healpix с помощью healpy

Question

У меня есть набор больших>10M объектов, файлов с RA и склонениями. Я хотел бы сделать из них карты плотности для всего неба, используя, я полагаю, healpix/healpy. Мой текущий код выглядит так:

 m = hp.ang2pix(512, ra, dec, lonlat=True)
 NSIDE = 512
 np.arange(hp.nside2npix(NSIDE))
 hp.visufunc.mollview(m) 

и я получаю ошибку:

 ValueError: Wrong pixel number (it is not 12*nside**2)

Что я делаю неправильно??

Спасибо ник

Answer 1

Здесь m - массив длины ra, (и dec). Сначала вам нужно конвертировать m в карту healpix, [или массив] длиной 12*NSIDE^2.

Для этого вы можете использовать numpy.bincount [очень быстро, и дает вам количество объектов в каждом пикселе], или scipy.stats.binned_statistic, [очень медленно, но позволяет вычислять любую "статистику", например, np.std. и т. д. вы хотите, данных, которые находятся в каждом пикселе]

def gen_fast_map(ip_, nside=512):
    npixel  = hp.nside2npix(nside)
    map_ = np.bincount(ip_,minlength=npixel)
    return map_

map = gen_fast_map(m)
hp.visufunc.mollview(map)