Handle white noise with healpy 3 not-uniform and partial sky coverage

import healpy as hp
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
import astropy.units as u

In this third notebook, we will handle a case of not-uniform and partial sky coverage.

# Number based on Simons Observatory SAT UHF1 array of detectors

net = 10. * u.Unit("uK * sqrt(s)")

5 years with a efficiency of 20%:

integration_time_total = 5 * u.year * .2

Download a hitmap

We can download a simulated hitmap for a Simons Observatory band, for now however, we assume a uniform coverage.

hitmap_url = "https://portal.nersc.gov/project/sobs/so_mapsims_data/v0.2/healpix/ST0_UHF1_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz"

!wget $hitmap_url

--2020-06-20 14:07:02--  https://portal.nersc.gov/project/sobs/so_mapsims_data/v0.2/healpix/ST0_UHF1_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz
Resolving portal.nersc.gov (portal.nersc.gov)... 128.55.201.128
Connecting to portal.nersc.gov (portal.nersc.gov)|128.55.201.128|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 3212515 (3.1M) [application/x-gzip]
Saving to: ‘ST0_UHF1_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz.5’

ST0_UHF1_01_of_20.n 100%[===================>]   3.06M  --.-KB/s    in 0.1s    

2020-06-20 14:07:02 (29.0 MB/s) - ‘ST0_UHF1_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz.5’ saved [3212515/3212515]

hitmap = hp.read_map("ST0_UHF1_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz")

/home/zonca/zonca/p/software/healpy/healpyvenv/lib/python3.7/site-packages/healpy/fitsfunc.py:352: UserWarning: If you are not specifying the input dtype and using the default np.float64 dtype of read_map(), please consider that it will change in a future version to None as to keep the same dtype of the input file: please explicitly set the dtype if it is important to you.
  "If you are not specifying the input dtype and using the default "

NSIDE = 512
ORDERING = NESTED in fits file
INDXSCHM = IMPLICIT
Ordering converted to RING

hitmap /= hitmap.max()

hp.mollview(hitmap)

Generic partial sky survey

We have now a sky coverage which not uniform

nside = 512
npix = hp.nside2npix(nside)
npix_hit = hitmap.sum()
ave_hit = hitmap.mean()

hitmap = hitmap / hitmap.sum() * integration_time_total.to(u.s)

hp.mollview(hitmap, unit=hitmap.unit)

variance_per_pixel = \
    (net**2 / hitmap).decompose()

/home/zonca/zonca/p/software/healpy/healpyvenv/lib/python3.7/site-packages/astropy/units/quantity.py:477: RuntimeWarning: divide by zero encountered in true_divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)

variance_per_pixel[np.isinf(variance_per_pixel)] = 0

m = np.random.normal(scale = np.sqrt(variance_per_pixel),
                     size=len(variance_per_pixel)) * np.sqrt(variance_per_pixel).unit

variance_per_pixel.max()

$4.51983 \times 10^{-7} \; \mathrm{K^{2}}$

m.value[hitmap==0] = hp.UNSEEN

m = m.to(u.uK)

m.value[hitmap==0] = hp.UNSEEN

hp.mollview(m, unit=m.unit, min=-10, max=10, title="White noise map")

Power spectrum

sky_fraction = hp.mask_good(m).sum()/len(m)

cl = hp.anafast(m) / sky_fraction

cl[100:].mean()

0.004372360023653951

pixel_area = hp.nside2pixarea(nside)

white_noise_cl = (variance_per_pixel[variance_per_pixel>0].mean() * pixel_area).to(u.uK**2)

white_noise_cl

$0.0044656828 \; \mathrm{\mu K^{2}}$

white_noise_cl_uniform = 1.5341266e-5 * u.uK**2

plt.figure(figsize=(6,4))
plt.loglog(cl, label="Map power spectrum", alpha=.7)
plt.hlines(white_noise_cl.value, 0, len(cl), color="blue",
           label="White noise level")
plt.hlines(white_noise_cl_uniform.value, 0, len(cl),
           label="White noise level uniform")
plt.title("Fullsky white noise spectrum")
plt.xlabel("$\ell$")
plt.ylabel("$C_\ell [\mu K ^ 2]$");

white_noise_cl / white_noise_cl_full_sky

$43.205324 \; \mathrm{}$

sky_fraction

0.38526121775309247

cl_apodized = hp.anafast(m * hitmap) / np.mean(hitmap**2) / sky_fraction

plt.figure(figsize=(6,4))
plt.loglog(cl, label="Map power spectrum", alpha=.7)
plt.loglog(cl_apodized, label="Apodized map", alpha=.7)

plt.hlines(white_noise_cl.value, 0, len(cl), color="blue",
           label="White noise level")
plt.hlines(white_noise_cl_uniform.value, 0, len(cl),
           label="White noise level uniform")
plt.title("Fullsky white noise spectrum")
plt.legend()
plt.xlabel("$\ell$")
plt.ylabel("$C_\ell [\mu K ^ 2]$");