Wie behandelt Healpix NaN, um Daten auf der Himmelskarte zu interpolieren?

Question

Ich habe meine Daten in diesem Format datapoint[37][19] im phi-theta Raum. Aber weil meine Daten nicht den gesamten Himmel abdecken können, gibt es im Array datapoint einige NaN. Es gibt ungefähr die Hälfte NaN im ganzen datapoint. Etwa 9/10 Nicht-NaN-Werte in datapoint sind negativ, etwa 1/10 von ihnen sind positiv.

habe ich versucht, diese Interpolation Funktion:

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

Aber es kehrte erros. Ich frage, wie interpoliere ich die Daten, die NaN enthalten, positive und negative Werte auf die Ebene, die Healpix verwenden kann, um Karten zu erstellen. Ich habe eine ungeglättete Karte aus Basemap erstellt.

Answer 1

Ich weiß nicht, ob Ihre Frage zu 1 ist) mit fehlenden Daten handelt, 2) mit NaNs tun hat, oder 3) auf der Kugel in eine Healpix Karte beliebige Daten drehen.

1) Die Interpolation fehlender Daten in einem großen Bereich des Himmels erfordert mindestens ein gewisses statistisches Wissen über Ihre Daten, um eine eingeschränkte Realisierung dessen zu ermöglichen, was fehlt. Das Füllen dieser Lücken ist jedoch nur erforderlich, wenn Sie nicht-lokale Operationen ausführen (wie Faltungen oder Verläufe). Es hängt also davon ab, was Sie mit den Daten tun möchten.

2) Das Setzen der fehlenden Daten auf NaN wird sicherlich alle verfügbaren Interpolationsschemata vermasseln.

3) Der folgende Python-Code wandelt einen Datensatz in 2 Healpix-Maps um, der eine NGP-Abtastung (Next Grip Point) und eine BSpline-Interpolation verwendet. Beachten Sie, dass die zweite wahrscheinlich nicht in Anwesenheit von NaNs, arbeiten wird, während die erste sehr robust ist.

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()