3

I've got a high-resolution healpix map (nside = 4096) that I want to smooth in disks of a given radius, let's say 10 arcmin.

Being very new to healpy and having read the documentation I found that one - not so good - way to do this was to perform a "cone search", that is to find around each pixels the ones inside the disk, average them and give this new value to the pixel at the center. However this is very time-consuming.

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

I think the healpy.sphtfunc.smoothing function could be of some help as it states that you can enter any custom beam window function but I don't understand how this works at all...

Thanks a lot for your help !

Hermanter
  • 225
  • 1
  • 9
  • You could try a [convolution](https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.signal.convolve2d.html) with a "round" kernel. – warped May 14 '19 at 21:49
  • Smoothing in healpy is usually done in Fourier space, where the window function is the Fourier transform of your smoothing kernel (typically a Gaussian). If you wanted to use a tophat filter, you'll have to calculate the transform of that, which will be a sinc function (https://en.wikipedia.org/wiki/Rectangular_function#Fourier_transform_of_the_rectangular_function). – Daniel Lenz May 14 '19 at 22:10
  • I'm afraid I don't have time to look into this in more detail, but are you 100% sure you need a tophat filter and a Gaussian window does not work for your application? – Daniel Lenz May 14 '19 at 22:11
  • @DanielLenz Unfortunately yes I am sure I need a top-hat... The thing is I'm not sure how to actually perform the convolution in Fourier space since it's not just a "flat plane" nor a regular grid... Can I just treat everything as it was for the convolution ? Shouldn't it be in spherical harmonics space ? – Hermanter May 14 '19 at 22:20
  • Some background on this subject can be found e.g. in White (1992, https://journals.aps.org/prd/pdf/10.1103/PhysRevD.46.4198), his Equations 40-46. The idea is as follow: (i) You calculate your beam window function in Fourier space (this is the tricky part, and is described in the reference for a Gaussian beam (ii) You convert your map to the spherical harmonic coefficients a_lm, then multiply your alm with your window function, and then convert back to the map. hp.smoothing() does all of this for you, you just need to provide the correct window function. – Daniel Lenz May 15 '19 at 16:33
  • What I'd also suggest is to raise an issue on the healpy github (https://github.com/healpy/healpy). That's likely a better place for a more in-depth discussion. Others might have gone through something similar, and perhaps this could be an interesting feature request as well. – Daniel Lenz May 15 '19 at 16:35

2 Answers2

1

As suggested, I can easily make use of the healpy.sphtfunc.smoothing function by specifying a custom (circular) beam window.

To compute the beam window, which was my problem, healpy.sphtfunc.beam2bl is very useful and simple in the case of a top-hat.

The appropriated l_max would roughly be 2*Nside but it can be smaller depending on specific maps. One could for example compute the angular power-spectra (the Cls) and check if it dampens for smaller l than l_max which could help gain some more time.

Thanks a lot to everyone who helped in the comments section!

Hermanter
  • 225
  • 1
  • 9
  • Can I suggest you add a code snippet to illustrate what you did? Maybe just generate a healpy map with random numbers and then smooth it with a circular beam? ` healpy.sphtfunc.beam2bl` is not super clear to the rest of us. – I.P. Freeley Nov 14 '19 at 17:23
1

since I spent a certain amount of time trying to figure out how the function smoothing was working. There is a bit of code that allows you to do a top_hat smoothing.

Cheers,

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()
Sullivan Marafico
  • 11
  • 1