Как Healpix справляется с NaN для интерполяции данных на карте неба?

У меня есть данные в этом формате datapoint[37][19] в phi-theta пространство. Но поскольку мои данные не могут охватить все небо, поэтому в datapoint массив. В целом около половины NaN datapoint, Около 9/10 не-NaN значений в datapoint отрицательные, около 1/10 из них являются положительными.

Я попробовал эту функцию интерполяции:

scipy.interpolate.RectSphereBivariateSpline(theta,phi,datapoint.T)

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

введите описание изображения здесь

1 ответ

Я не знаю, касается ли ваш вопрос 1) работы с недостающими данными, 2) работы с NaN или 3) преобразования произвольных данных на сфере в карту Healpix.

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

2) Установка недостающих данных в NaN, безусловно, испортит все доступные схемы интерполяции.

3) Код Python ниже превращает набор данных, похожий на ваш (насколько я могу судить) в 2 карты Healpix, одна из которых использует выборку Nearest Grip Point (NGP), а другая - интерполяцию BSpline. Обратите внимание, что второе, скорее всего, не будет работать в присутствии NaN, в то время как первое чрезвычайно устойчиво.

import healpy as hp
import numpy as np
import pylab as pl

datapoint = np.zeros((37,19), dtype=np.float)
datapoint[18,9] = 1.0
datapoint[0,9] = -1.0

nside      = 64
npix       = hp.nside2npix(nside)

# location of Healpix pixels center
ip         = np.arange(npix)
theta_rad, phi_rad = hp.pix2ang(nside, ip)

# map0 : NGP sampling
theta_deg = np.rad2deg(theta_rad)
phi_deg   = np.rad2deg(phi_rad)
hp_0      = datapoint[np.rint(phi_deg/10.).astype(int), \
                      np.rint(theta_deg/10.).astype(int)]
hp.mollview(hp_0,title='NGP map')

# map1: BSpline interpolation
from scipy.interpolate import RectSphereBivariateSpline
epsilon = 1.e-12
th_in = np.linspace(epsilon,  np.pi-epsilon, 19)
ph_in = np.linspace(epsilon,2*np.pi-epsilon, 37)
lut   = RectSphereBivariateSpline(th_in, ph_in, datapoint.T, s=1)
hp_1  = lut.ev(theta_rad, phi_rad)
hp.mollview(hp_1,title='BSpline map')

pl.show()
Другие вопросы по тегам