Generate point source maps with pixell

cosmology
pysm
Published

May 2, 2024

Testing the pixell sim_objects functionality to create maps of point sources pre-smoothed with a gaussian beam. The purpose is to include this functionality in PySM to be able to generate on the fly maps of source starting from a catalog.

from pixell import enmap, utils, resample, curvedsky as cs, reproject, pointsrcs
import numpy as np
import healpy as hp
fwhm = 5 * utils.degree
shape, wcs = enmap.fullsky_geometry(res=fwhm / 3, proj="car")
shape
(109, 216)
def fwhm2sigma(fwhm):
    return fwhm / (2.0 * np.sqrt(2.0 * np.log(2.0)))
def flux2amp(flux, fwhm):
    sigma = fwhm2sigma(fwhm)
    return flux / (2 * np.pi * sigma**2)
assert flux2amp((2 * np.pi * fwhm2sigma(5) ** 2), 5) == 1
n_sources = 1
flux_sources = np.arange(n_sources) + 10
amplitude_sources = flux2amp(flux_sources, fwhm)
source_pos = np.array([[0], [np.pi / 3]])
r, p = pointsrcs.expand_beam(fwhm2sigma(fwhm))
source_map = pointsrcs.sim_objects(shape, wcs, source_pos, amplitude_sources, ((r, p)))
import matplotlib.pyplot as plt
source_pos
array([[0.        ],
       [1.04719755]])
plt.imshow(source_map)

source_map.argmax(unit="coord")
array([0.        , 1.04719755])
source_map.argmax(unit="pix")
array([54, 72])
source_pos[:, -1]
array([0.        , 1.04719755])
-source_map.argmax(unit="coord") + source_pos[:, -1]
array([ 0.00000000e+00, -6.66133815e-16])
source_map.max()
array(1158.8864, dtype=float32)
source_map.min()
array(0., dtype=float32)
def aperture_photometry(
    thumbs, aperture_radius, annulus_width=None, modrmap=None, pixsizemap=None
):
    """
    Flux from aperture photometry.

    from https://github.com/msyriac/orphics/blob/master/orphics/maps.py

    Parameters
    ----------
    thumb : ndmap
        An (...,Ny,Nx) ndmap (i.e. a pixell enmap) containing the thumbnails.
    aperture_radius : float
        Aperture inner radius in radians
    annulus_width : float
        Annulus width for mean subtraction in radians.
        Defaults to sqrt(2)-1 times the aperture inner radius.
    modrmap : ndmap, optional
        An (Ny,Nx) ndmap containing distances of each pixel from the center in radians.
    modrmap : ndmap, optional
        An (Ny,Nx) ndmap containing pixel areas in steradians.

    Returns
    -------
    flux : ndarray
        (...,) array of aperture photometry fluxes.

    """
    if modrmap is None:
        modrmap = thumbs.modrmap()
    if annulus_width is None:
        annulus_width = (np.sqrt(2.0) - 1.0) * aperture_radius
    # Get the mean background level from the annulus
    mean = thumbs[
        ...,
        np.logical_and(
            modrmap > aperture_radius, modrmap < (aperture_radius + annulus_width)
        ),
    ].mean()
    if pixsizemap is None:
        pixsizemap = thumbs.pixsizemap()
    # Subtract the mean, multiply by pixel areas and sum
    return (((thumbs - mean) * pixsizemap)[..., modrmap <= aperture_radius]).sum(
        axis=-1
    )
from astropy import units as u
box_half_size_rad = 2 * fwhm
box_center = [source_pos[0, -1], source_pos[1, -1]]
box = np.array(
    [
        [box_center[0] - box_half_size_rad, box_center[1] - box_half_size_rad],
        [box_center[0] + box_half_size_rad, box_center[1] + box_half_size_rad],
    ]
)  # in radians
box_center
[0.0, 1.0471975511965976]
cutout = source_map.submap(box)
cutout.max()
array(1158.8864, dtype=float32)
cutout.min()
array(0., dtype=float32)
plt.imshow(cutout)

aperture_photometry(cutout, 2 * fwhm)
9.99300220192394