0

我datapoint[37][19]在phi-theta空间中有这种格式的数据。但是因为我的数据不能覆盖整个天空，所以datapoint数组中有一些NaN。总共大约有一半 NaN datapoint。大约 9/10 的非 NaN 值datapoint是负数，其中大约 1/10 是正数。

我试过这个插值函数：

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

但它返回了错误。我在问如何将包含 NaN、正值和负值的数据插入到 Healpix 可以使用它来制作地图的级别。我有一张由底图制作的未平滑地图。

在此处输入图像描述

4

1 回答 1

0

我不知道您的问题是关于 1）处理丢失的数据，2）处理 NaN，还是 3）将球体上的任意数据转换为 Healpix 地图。

1）在大面积的天空中插入缺失的数据至少需要对数据有一些统计知识，以便对缺失的内容进行有约束的实现。但只有在进行非局部操作（如卷积或梯度）时才需要填补这些空白，因此这取决于您打算如何处理数据。

2) 将缺失数据设置为 NaN 肯定会搞砸所有可用的插值方案。

3) 下面的 python 代码将与您的数据集相似的数据集（据我所知）转换为 2 张 Healpix 地图，一张使用最近夹点 (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()
于 2016-11-02T11:26:29.677 回答