xuod/castor

A collection of useful stuff.

castor's Issues

Function to convert MANGLE to HEALPIX

Hi!
I am looking for a function which given MANGLE-format maps and translates it into healpix format with given nside. I cannot use mangle software on my 64-bit machine.

I came across you castor package and tried to use function ply2hpx.
I tried to run in as is, copying needed functions as chunk from math module, etc, and found that it does not work for several reasons.

Therefore I made a little re-writing of your code which you can find below (in the comments I've written how I fixed the errors which had made code to stop).
After running it on the SDSS DR12 LOWZ acceptance mask mask_DR12v5_LOWZ_North.ply and mask_DR12v5_LOWZ_South.ply I got the following mask:

This seems good at first glance, but I am confused due to the presence of 'zero-value' strips which are parallel to the coordinate system's lines of equal longitude.

The question is whether it is normal to have such strips due to the peculiarities of MANGLE format, or I did something wrong?
In papers which use LOWZ and/or CMASS SDSS samples, I do not observe such strips.

Thanks!

Appendix: code that does not have errors (on my machine).

import numpy as np
import healpy as hp
from astropy import units as u
import astropy.coordinates as coord
import pymangle

def mangle2healpix(ply: pymangle.Mangle,
                   nside: int, fk52gal: bool = False) -> np.ndarray:
    """
    mangle2healpix translate Mangle-format map (.ply) into healpix map with given nside


    # see https://github.com/xuod/castor/blob/master/castor/cosmo.py package by Cyrille Doux (https://github.com/xuod?tab=repositories). Functions from accirting submodules of 'castor' package.

    Args:
        ply (pymangle.Mangle): pymangle map object
        nside (int): nside for output map
        fk52gal (bool, optional): if transform  galactic coordinates into fk5 coordinatees. Defaults to False.

    Returns:
        np.ndarray: healpix map of the input mangle map
    """

    def ply2hpx(ply, nside, fk52gal=False):
        """
        Fast convertion of a mangle .ply mask into a Healpix map using the weight of
        the polygon where centers of healpix pixels fall.
        Parameters
        ----------
        ply : mangle.Mangle
            mangle.Mangle object.
        nside : int
            `nside` of output healpix map.
        fk52gal : bool
            If true, transform galactic coordinates to fk5 coordinates, useful for galaxy surveys (the default is False).
        Returns
        -------
        array
            Healpix map.
        """

        def chunk(seq, m):
            """
            Divide a list `seq` into m chuncks evenly-sized chuncks.
            Parameters
            ----------
            seq : list
            m : int
            Returns
            -------
            list
                List of sub-lists.
            See also np.array_split
            """
            n = len(seq)
            p = n / m
            r = n % m
            out = []
            last = 0

            for _ in range(m):
                if r > 0:
                    length = p + 1
                else:
                    length = p
                # out.append(seq[last:last + length]) #before

                # after adding int()
                out.append(seq[int(last):int(last + length)])

                last += length
                r -= 1

            return out

        def thetaphi2radec(theta, phi):
            """
            Converts theta and phi (Healpix conventions, in radians) to ra and dec (in degrees)
            """
            ra = np.rad2deg(phi)
            dec = 90. - np.rad2deg(theta)
            return ra, dec

        def _ply2hp_aux(args):
            ply, ra, dec = args
            # return ply.get_polyids(ra, dec) # before

            # after. get_polyids method is no more in current vrsion of pymangle
            return ply.polyid(ra, dec)

        npix = hp.nside2npix(nside)
        theta, phi = hp.pix2ang(nside, np.arange(npix))
        ra, dec = thetaphi2radec(theta, phi)

        if fk52gal:
            coord_gal = coord.SkyCoord(coord.Angle(
                ra*u.degree), coord.Angle(dec*u.degree), frame='galactic')
            coord_fk5 = coord_gal.transform_to('fk5')
            ra = coord_fk5.ra.deg
            dec = coord_fk5.dec.deg

        ral = chunk(ra, 100)
        decl = chunk(dec, 100)

        arglist = [(ply, ral[i], decl[i]) for i in range(100)]

        # ids = np.concatenate(parallel.map(_ply2hp_aux, arglist))  # timesleep=1.0) #before. some strangee error when use parallel.map from 'castor' package: cannot pickle 'Mangle' object

        # after
        ids = np.array([_ply2hp_aux(arg)
                        for arg in arglist])  # not parallelized. should be fine for not very large nside
        ids = np.concatenate(ids)

        pos = np.where(ids > -1)

        hpmap = np.zeros(npix)
        hpmap[pos] = ply.weights[ids[pos]]

        return hpmap

    return ply2hpx(ply=ply, nside=nside, fk52gal=fk52gal)


#masks are here: https://data.sdss.org/sas/dr12/boss/lss/
name = 'LOWZ'
mask1 = pymangle.Mangle(f'/masks/mask_DR12v5_{name}_North.ply')
hpmask1 = mangle2healpix(mask1, nside)
mask2 = pymangle.Mangle(f'/masks/mask_DR12v5_{name}_South.ply')
hpmask2 = mangle2healpix(mask2, nside)

assert np.all(hpmask1*hpmask2 == 0)

lowz_acceptance = hpmask1 + hpmask2

del hpmask1, hpmask2, mask1, mask2, name

hp.mollview(lowz_acceptance, title='LOWZ ACCEPTANCE')