如何有效地在 healpix 地图上执行顶帽(磁盘状)平滑?
How to efficiently perform a top-hat (disk-like) smoothing on a healpix map?
我有一个高分辨率的 healpix 贴图 (nside = 4096),我想在给定半径的磁盘中对其进行平滑处理,比如说 10 arcmin。
作为 healpy 的新手并阅读了文档,我发现一个 - 不太好的 - 方法是执行 "cone search",即在每个像素周围找到磁盘内的像素,对它们进行平均并将这个新值赋予中心的像素。然而,这非常耗时。
import numpy as np
import healpy as hp
kappa = hp.read_map("zs_1.0334.fits") #Reading my file
NSIDE = 4096
t = 0.00290888 #10 arcmin
new_array = []
n = len(kappa)
for i in range(n):
a = hp.query_disc(NSIDE,hp.pix2vec(NSIDE,i),t)
new_array.append(np.mean(kappa[a]))
我认为 healpy.sphtfunc.smoothing 函数可能会有一些帮助,因为它声明您可以输入任何自定义光束 window 函数,但我完全不明白它是如何工作的...
非常感谢您的帮助!
按照建议,我可以通过指定自定义(圆形)梁 window.
轻松使用 healpy.sphtfunc.smoothing 函数
要计算光束 window,这是我的问题,healpy.sphtfunc.beam2bl 在 top-hat.
的情况下非常有用且简单
适当的l_max大概是2*Nside,但根据具体的地图可能会更小。例如,可以计算 angular power-spectra(Cls)并检查它是否因小于 l_max 的 l 而衰减,这可能有助于获得更多时间。
非常感谢在评论区给予帮助的大家!
因为我花了一定时间试图弄清楚函数平滑是如何工作的。有一些代码可以让您进行 top_hat 平滑。
干杯,
import healpy as hp
import numpy as np
import matplotlib.pyplot as plt
def top_hat(b, radius):
return np.where(abs(b)<=radius, 1, 0)
nside = 128
npix = hp.nside2npix(nside)
#create a empy map
tst_map = np.zeros(npix)
#put a source in the middle of the map with value = 100
pix = hp.ang2pix(nside, np.pi/2, 0)
tst_map[pix] = 100
#Compute the window function in the harmonic spherical space which will smooth the map.
b = np.linspace(0,np.pi,10000)
bw = top_hat(b, np.radians(45)) #top_hat function of radius 45°
beam = hp.sphtfunc.beam2bl(bw, b, nside*3)
#Smooth map
tst_map_smoothed = hp.smoothing(tst_map, beam_window=beam)
hp.mollview(tst_map_smoothed)
plt.show()
我有一个高分辨率的 healpix 贴图 (nside = 4096),我想在给定半径的磁盘中对其进行平滑处理,比如说 10 arcmin。
作为 healpy 的新手并阅读了文档,我发现一个 - 不太好的 - 方法是执行 "cone search",即在每个像素周围找到磁盘内的像素,对它们进行平均并将这个新值赋予中心的像素。然而,这非常耗时。
import numpy as np
import healpy as hp
kappa = hp.read_map("zs_1.0334.fits") #Reading my file
NSIDE = 4096
t = 0.00290888 #10 arcmin
new_array = []
n = len(kappa)
for i in range(n):
a = hp.query_disc(NSIDE,hp.pix2vec(NSIDE,i),t)
new_array.append(np.mean(kappa[a]))
我认为 healpy.sphtfunc.smoothing 函数可能会有一些帮助,因为它声明您可以输入任何自定义光束 window 函数,但我完全不明白它是如何工作的...
非常感谢您的帮助!
按照建议,我可以通过指定自定义(圆形)梁 window.
轻松使用 healpy.sphtfunc.smoothing 函数要计算光束 window,这是我的问题,healpy.sphtfunc.beam2bl 在 top-hat.
的情况下非常有用且简单适当的l_max大概是2*Nside,但根据具体的地图可能会更小。例如,可以计算 angular power-spectra(Cls)并检查它是否因小于 l_max 的 l 而衰减,这可能有助于获得更多时间。
非常感谢在评论区给予帮助的大家!
因为我花了一定时间试图弄清楚函数平滑是如何工作的。有一些代码可以让您进行 top_hat 平滑。
干杯,
import healpy as hp
import numpy as np
import matplotlib.pyplot as plt
def top_hat(b, radius):
return np.where(abs(b)<=radius, 1, 0)
nside = 128
npix = hp.nside2npix(nside)
#create a empy map
tst_map = np.zeros(npix)
#put a source in the middle of the map with value = 100
pix = hp.ang2pix(nside, np.pi/2, 0)
tst_map[pix] = 100
#Compute the window function in the harmonic spherical space which will smooth the map.
b = np.linspace(0,np.pi,10000)
bw = top_hat(b, np.radians(45)) #top_hat function of radius 45°
beam = hp.sphtfunc.beam2bl(bw, b, nside*3)
#Smooth map
tst_map_smoothed = hp.smoothing(tst_map, beam_window=beam)
hp.mollview(tst_map_smoothed)
plt.show()