add_catalog_sources.py 5.49 KB
Newer Older
1 2 3 4
from nikamap import NikaMap
from astropy.table import Table
import numpy as np
import astropy.units as u
5 6 7 8 9 10
from nikamap.utils import pos_cat
import matplotlib.pyplot as plt
from astropy.wcs import WCS

#%load_ext autoreload
#%matplotlib tk
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50


def GetInRange(array, low=None, high=None, extra=None):

    lmask = np.zeros(len(array), dtype=bool)
    hmask = np.zeros(len(array), dtype=bool)

    if low is not None:
        lmask = array < low

    if high is not None:
        hmask = array > high

    tmask = lmask | hmask

    if extra is not None:
        return array[~tmask], extra[~tmask]
    else:
        return array[~tmask]


def uniform(a, b, size=None):
    return (b - a) * np.random.random_sample(size) + a


def rel_to_coord(a, b, relcoord):
    return (b - a) * relcoord + a

'''
def Shift(biglimits, smalllimits, relcoord):
    shiftrange = np.sort(biglimits-smalllimits)
    shift = rel_to_coord(*shiftrange)
    return shift
'''
def Shift(shiftrange, relcoord):
    shift = rel_to_coord(*shiftrange, relcoord)
    return shift


class CatalogShifter:
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
    def __init__(self, catalog, mapshape, wcs, coverage=1):
        '''
        Class that takes sources catalog and shifts source positions to area
        that is covered by a map with shape ``mapshape`` and `astropy.wcs.WCS`
        ``wcs``

        Parameters
        ----------
        catalog : :class:`astropy.table.Table`
            Table containing the source positions. Must provide columns
            'ra' and 'dec'
        mapshape : tuple
            shape of the map that shall contain shifted sources
        wcs : :class:`astropy.wcs.WCS`
            wcs of the map
        coverage : float
            value between [0, 1]. At least (1-coverage)**2 * mapsurface
            is inside area of catalog sources.

        Returns
        -------
        CatalogShifter object that can be used to create shifted catalogs
        '''

        # borders of empty map
76
        maplimits = wcs.wcs_pix2world(*(np.array(((0, 0), mapshape)).T), 0)
77
        maplimits = np.sort(maplimits, axis=-1)
78 79

        ra, dec = catalog['ra'], catalog['dec']
80 81 82

        # borders of source area
        sourcelimits = np.array([[np.amin(ra), np.amax(ra)],
83 84
                                 [np.amin(dec), np.amax(dec)]])

85 86 87 88 89 90 91
        # shiftlimits for full coverage
        shiftlimits = np.sort(maplimits - sourcelimits, axis=-1)

        # modify limits if partial coverage allowed
        mapsize = maplimits[:, 1] - maplimits[:, 0]
        shiftlimits[:, 0] -= (mapsize * (1-coverage))
        shiftlimits[:, 1] += (mapsize * (1-coverage))
92 93 94 95 96

        self.catalog = catalog
        self._shiftlimits = shiftlimits

    def __call__(self, relra=None, reldec=None):
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
        '''
        Parameters
        ----------
        relra : float between [0, 1]
            Optional. Relative RA position of shifted sources in map. If
            ``None`` it is chosen randomly.
        reldec : float between [0, 1]
            Optional. Relative DEC position of shifted sources in map. If
            ``None`` it is chosen randomly.

        Returns
        -------
        catalog : :class:``astropy.table.Table``
            catalog containing new source positions
        '''
112 113 114 115 116 117 118 119 120

        ra_shiftlimits, dec_shiftlimits = self._shiftlimits
        outcat = self.catalog   # this should create a copy

        if relra is None:
            relra = np.random.random_sample()
        if reldec is None:
            reldec = np.random.random_sample()

121 122
        rashift = Shift(ra_shiftlimits, relra)
        decshift = Shift(dec_shiftlimits, reldec)
123 124 125 126 127 128 129

        outcat['ra'] = outcat['ra'] + rashift
        outcat['dec'] = outcat['dec'] + decshift

        return outcat


130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
def CatalogFor_pos_cat(catalog, fluxkey='flux', fluxunit=None, minflux=None):
    catkeys = catalog.keys()

    if '_ra' not in catkeys:
        assert 'ra' in catkeys, 'you must provide RA for source'
        catalog.rename_column('ra', '_ra')

    if '_dec' not in catkeys:
        assert 'dec' in catkeys, 'you must provide DEC for source'
        catalog.rename_column('dec', '_dec')

    assert fluxkey in catkeys, 'fluxkey not matching with any column name'
    catalog.rename_column(fluxkey, 'flux')

    if catalog['flux'].unit is None:
        catalog['flux'].unit = fluxunit
    if minflux is not None:
        catalog = catalog[u.Quantity(catalog['flux']) > minflux]

    return catalog


def test_catalogshifter():
    wcs = WCS()
    data = np.zeros((2, 2))
    ra, dec = np.meshgrid(np.arange(4), np.arange(4))
    ra = ra.reshape(ra.size)
    dec = dec.reshape(dec.size)
    flux = 10 * ra + dec
    catalog = Table()
    catalog['flux'] = flux
    catalog['ra'] = ra
    catalog['dec'] = dec
    sh = CatalogShifter(catalog, (2, 2), wcs)
    sh(0,0, coverage=0)
    #print(sh(0,0, coverage=1))


168
if __name__ == "__main__":
169 170
    # test_catalogshifter()

171
    nmfile = "/home/peter/Dokumente/Uni/Paris/Stage/data/map.fits"
172 173
    nmfile = ("/home/peter/Dokumente/Uni/Paris/Stage/N2CLS/ForMarseille/"
              "HLS091828_common_mode_kids_out/map_JK.fits")
174 175 176 177 178 179 180
    tablefile = "/home/peter/Dokumente/Uni/Paris/Stage/N2CLS/SIDES_NIKA.fits"
    nm = NikaMap.read(nmfile)
    wcs = nm.wcs
    shape = nm.shape
    catalog = Table.read(tablefile)
    catalog.keys()

181 182 183 184 185 186
    sh = CatalogShifter(catalog, shape, wcs, coverage=.5)
    cat = CatalogFor_pos_cat(sh(1, 1), fluxkey='SNIKA1200',
                             fluxunit=u.Jy, minflux=1*u.mJy)
    print(np.amin(cat['_ra']))
    nm.add_gaussian_sources(cat_gen=pos_cat, catalog=cat, wcs=nm.wcs)

187 188
    # %matplotlib tk
    nm.plot()
189
    plt.show(block=True)