2

I am trying to compute angular power spectrum from a masked map using anafast in healpy. How does the python version of anafast take care of the effects of masking, compared to e.g., the F90 version, which explicitly takes the mask file as an optional input?

Thanks!

Progman
  • 16,827
  • 6
  • 33
  • 48
Erwin Lau
  • 21
  • 2

2 Answers2

2

It's easiest to simply set all masked values to hp.UNSEEN or to 0 (internally, all UNSEEN values will be set to 0 anyways during hp.anafast()). In the example below, I use numpy.where().

import numpy as np
import healpy as hp

masked_map = np.where(mask, raw_map, hp.UNSEEN)
cl = hp.anafast(masked_map)

You could also pass the mask by turning the input map into a masked array, which does effectively the same as the approach described above.

Internally, all masked values are set to 0 for the computation of the power spectrum, so be aware of maps that do not have zero mean. In any case, you might want to subtract the monopole before computing the power spectrum.

Daniel Lenz
  • 3,334
  • 17
  • 36
1

In my codes, I usually use the masked array to define masks. An example code can be like:

import healpy as hp
import numpy as np

mask = hp.read_map("mask.fits")
map = hp.read_map("map.fits")
map_masked = hp.ma(map)
map_masked.mask = np.logical_not(mask)

C_l = hp.anafast(map_masked)

There is also the pymaster library in python, which is the implementation of NaMaster library. It is specifically used to compute various power spectra with mapped maps and binning. You define an nmt field object with a list of maps and the mask. Then use various built-in functions to compute various spectra like compute_full_master. A sample:

import pymaster as nmt

f_0 = nmt.NmtFiels(mask, [map])
b =  nmt.NmtBin.from_lmax_linear(lmax, nlb=21)
C_l = nmt.compute_full_master(f_0, f_0, b)

Just note that using the pymaster library is a bit computationally expensive.

Jyothis Chandran
  • 11
  • 4