Применить вращение к карте HEALPix в healpy

Я понимаю, что об этом спрашивали много лет назад, но на этой неделе у меня была такая же проблема, и я нашел ваш пост. Я нашел несколько потенциальных решений, поэтому я поделюсь, если кто-то еще наткнется на это и найдет его полезным.

Решение 1. Это зависит от формата, в котором поступают ваши данные. Мои пришли в сетке (тета, фи).

import numpy as np
import healpy as H
map = <your original map>
nside = <your map resolution, mine=256>
npix = H.nside2npix(nside)
pix = N.arange(npix)
t,p = H.pix2ang(nside,pix) #theta, phi

r = H.Rotator(deg=True, rot=[<THETA ROTATION>, <PHI ROTATION>])

map_rot = np.zeros(npix)

for i in pix:
    trot, prot = r(t[i],p[i])
    tpix = int(trot*180./np.pi) #my data came in a theta, phi grid -- this finds its location there
    ppix = int(prot*180./np.pi)
    map_rot[i] = map[ppix,tpix] #this being the rright way round may need double-checking

Решение 2. Еще не закончили тестирование, но наткнулись на него ПОСЛЕ выполнения раздражающей работы, описанной выше...

map_rot = H.mollview(map,deg=True,rot=[<THETA>,<PHI>], return_projected_map=True)

который дает массив 2D numpy. Мне интересно узнать, как преобразовать это обратно в карту healpix...

Я нашел еще одно потенциальное решение после нескольких месяцев поиска. Я еще не проверял это очень много, поэтому, пожалуйста, будьте осторожны!

Решение 2 Саула, приведенное выше, является ключевым (отличное предложение!)

По сути, вы комбинируете функциональность healpy.mollview (также работают gnomview, cartview и orthview) с функциональностью reproject_to_healpix в пакете reproject (http://reproject.readthedocs.org/en/stable/).

Полученная карта подходит для моих угловых масштабов, но я не могу сказать, насколько точно преобразование по сравнению с другими методами.

-----Основная схема----------

Шаг 1. Считайте карту и создайте прямоугольный массив с помощью cartview. Как указал Саул выше, это также один из способов вращения. Если вы просто выполняете стандартное преобразование поворота/координаты, то все, что вам нужно, это ключевое слово coord. От небесных до галактических координат установите coord = ['C','G']

map_Gal = hp.cartview(map_Cel, coord=['C','G'], return_projected_map=True, xsize=desired_xsize, norm='hist',nest=False)

Шаг 2. Напишите шаблон заголовка all-sky FITS (как в примере ниже). Я написал, что моя карта имеет тот же средний размер в пикселях, что и моя желаемая карта HEALPix.

Шаг 3. Используйте reproject.transform_to_healpix

reproject включает функцию для сопоставления "обычного" массива (или файла FITS) с проекцией HEALPix. Объедините это с возможностью возврата массива, созданного healpy.mollview/cartview/orthview/gnomview, и вы сможете повернуть карту HEALPix из одной системы координат (небесной) в другую систему координат (галактическую).

map_Gal_HP, footprint_Gal_HP = rp.reproject_to_healpix((map_Gal, target_header), coord_system_out= 'GALACTIC', nside=nside, nested=False)

По сути, это сводится к этим двум командам. Однако вам придется создать заголовок шаблона, указав масштаб и размер в пикселях, соответствующие промежуточной карте всего неба, которую вы хотите создать.

-----Полный рабочий пример (формат блокнота iPython + пример данных FITS)------

https://github.com/aaroncnb/healpix_coordtrans_example/tree/master

Код там должен работать очень быстро, но это потому, что карты сильно деградировали. Я сделал то же самое для своих карт NSIDE 1024 и 2048, и это заняло около часа.

------Изображения до и после------

Эта функция, кажется, делает свое дело (достаточно медленно, но должно быть лучше, чем цикл for):

def rotate_map(hmap, rot_theta, rot_phi):
    """
    Take hmap (a healpix map array) and return another healpix map array 
    which is ordered such that it has been rotated in (theta, phi) by the 
    amounts given.
    """
    nside = hp.npix2nside(len(hmap))

    # Get theta, phi for non-rotated map
    t,p = hp.pix2ang(nside, np.arange(hp.nside2npix(nside))) #theta, phi

    # Define a rotator
    r = hp.Rotator(deg=False, rot=[rot_phi,rot_theta])

    # Get theta, phi under rotated co-ordinates
    trot, prot = r(t,p)

    # Interpolate map onto these co-ordinates
    rot_map = hp.get_interp_val(hmap, trot, prot)

    return rot_map

Использование этого для данных из PyGSM дает следующее:

hp.mollview(np.log(rotate_map(gsm.generated_map_data, 0,0)))

При вращении phi:

hp.mollview(np.log(rotate_map(gsm.generated_map_data, 0,np.pi)))

Или вращая theta:

hp.mollview(np.log(rotate_map(gsm.generated_map_data, np.pi/4,0)))