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

Question

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

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

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

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

Answer 1

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

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

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

3) Код ниже приведенного ниже кода преобразует набор данных, похожий на ваш (насколько я могу судить) на 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()