Catalogue conversion

Triumvirate reads catalogue data and files assuming they contain the Cartesian coordinate columns 'x', 'y' and 'z'.

If your catalogue is in sky coordinates 'RA', 'DEC' and 'Z' (case-sensitive), you can perform the coordinate transform following the demo Python script below, which can also read in multiple catalogue files using glob patterns. No additional packages need to be installed as the one(s) needed here are dependencies of triumvirate.

"""Construct a :class:`~triumvirate.catalogue.ParticleCatalogue` instance
from (a) catalogue file(s) (format: text, FITS, HDF5 etc.) containing sky
coordinates and redshifts.

"""
from glob import glob

import astropy.units as u
import numpy as np
from astropy.coordinates import SkyCoord
from astropy.cosmology import FlatLambdaCDM, Planck15, Planck18
from astropy.table import Table

from triumvirate.catalogue import ParticleCatalogue


def get_fiducial_cosmology(name='desi'):
    """Get the fiducial cosmological model from a predefined list.

    Parameters
    ----------
    name : {'desi', 'planck15', 'planck18'}, optional
        Name of the fiducial model (default is 'desi').

    Returns
    -------
    :class:`astropy.cosmology.Cosmology`
        Fiducial cosmological model.

    """
    COSMO = {
        'desi': FlatLambdaCDM(
            H0=67.36,
            Om0=0.3138,
            Ob0=0.0493,
            Tcmb0=2.7255,
            Neff=3.045998221453422,
            m_nu=[0.06, 0., 0.],
        ),
        'planck15': Planck15,
        'planck18': Planck18,
    }

    return COSMO[name]


def read_sky_catalogue(catalogue_filepath, cosmo='desi', nfile=None):
    r"""Read (a) catalogue files containing sky coordinates and add
    Cartesian coordinates in units of :math:`h^{-1}\,\mathrm{Mpc}`.

    Parameters
    ----------
    catalogue_filepath : str or Path-like
        File path or glob pattern to the catalogue(s) containing sky
        coordinate columns 'RA', 'DEC' and 'Z'.
    cosmo : {'desi', 'boss', 'planck15', 'planck18'}, optional
        Name of the fiducial cosmology (default is 'desi').
    nfile : int, optional
        Maximum number of files if multiple are read in.  If `None`
        (default), all matched files are read in and concatenated
        after sorting by filename.

    Returns
    -------
    :class:`numpy.ndarray`
        Catalogue as a structured data array containing the Cartesian
        coordinate columns 'x', 'y' and 'z'.

    """
    # Read the catalogue file(s) as a data table.
    datatabs = []
    ifile = 0
    for cfpath in sorted(glob(catalogue_filepath)):
        datatabs.append(Table.read(cfpath))
        ifile += 1
        if nfile is not None and ifile == nfile:
            break
    datatab = np.concatenate(datatabs)

    # Extract sky coordinates.
    ra, dec = datatab['RA'] * u.degree, datatab['DEC'] * u.degree

    # Calculate comoving distances from the given redshifts and fiducial
    # cosmological model.
    cosmo = get_fiducial_cosmology(cosmo)
    dc = cosmo.comoving_distance(datatab['Z'])

    # Convert sky to Cartesian coordinates and add to the data table.
    coords = SkyCoord(ra=ra, dec=dec, distance=dc, frame='icrs')

    datatab['x'] = coords.cartesian.x.value * cosmo.h
    datatab['y'] = coords.cartesian.y.value * cosmo.h
    datatab['z'] = coords.cartesian.z.value * cosmo.h

    return datatab


if __name__ == '__main__':

    catalogue_data = read_sky_catalogue("<catalogue-filepath>")

    catalogue = ParticleCatalogue(
        catalogue_data['x'], catalogue_data['y'], catalogue_data['z'],
        nz=catalogue_data['NX'],
        # Total systematic weight: e.g. WEIGHT_COMP * WEIGHT_SYS * WEIGHT_ZFAIL
        ws=catalogue_data['WEIGHT'],
        # Total clustering (optimal) weight: e.g. WEIGHT_FKP
        wc=catalogue_data['WEIGHT_FKP'],
    )

    # Continue your script after here...