Как 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()