add_catalog_sources.py 11 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
from nikamap.utils import pos_cat
import matplotlib.pyplot as plt
from astropy.wcs import WCS
8
from astropy.coordinates import SkyCoord
9
'''
10 11 12
%load_ext autoreload
%autoreload 2
%matplotlib tk
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 51

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


52
class OldCatalogShifter:
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
    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
78
        maplimits = wcs.wcs_pix2world(*(np.array(((0, 0), mapshape)).T), 0)
79
        maplimits = np.sort(maplimits, axis=-1)
80 81

        ra, dec = catalog['ra'], catalog['dec']
82 83 84

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

87 88 89 90 91 92 93
        # 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))
94 95 96 97 98

        self.catalog = catalog
        self._shiftlimits = shiftlimits

    def __call__(self, relra=None, reldec=None):
99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
        '''
        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
        '''
114 115 116 117 118 119 120 121 122

        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()

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

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

        return outcat


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 168 169
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))


170 171 172 173 174 175 176 177 178 179 180 181 182 183
class transUniformCoords():
    def __init__(self, catalog):

        ra, dec = u.Quantity(catalog['ra']), u.Quantity(catalog['dec'])

        ralimits = u.Quantity([np.amin(ra), np.amax(ra)]).to_value(u.deg)
        declimits = u.Quantity([np.amin(dec), np.amax(dec)]).to_value(u.rad)

        decuniformlimits = np.sort((1 - np.sin(declimits)) / 2)

        self.ralimits = ralimits
        self.decuniformlimits = decuniformlimits


LUSTIG Peter's avatar
LUSTIG Peter committed
184
    def __call__(self, relcoordinate, angle):
185 186 187 188 189 190 191 192 193
        '''
        Parameters
        ----------
        Coordinate: float between 0 and 1
        '''

        ralimits = self.ralimits
        decuniformlimits = self.decuniformlimits

LUSTIG Peter's avatar
LUSTIG Peter committed
194
        relra, reldec = relcoordinate
195 196
        ra = rel_to_coord(*ralimits, relra) * u.deg
        angle = angle * 360 * u.deg
LUSTIG Peter's avatar
LUSTIG Peter committed
197
        # decuniform_inrange = rel_to_coord(*decuniformlimits, reldec)
198 199 200 201 202 203 204 205 206 207

        # rescale uniform coordinate to hold in range
        reldec = rel_to_coord(*decuniformlimits, reldec)
        dec = - np.arcsin(2 * reldec - 1)
        # dec = - np.arcsin(2 * decuniformlimits - 1)
        dec = u.Quantity(dec, unit=u.rad).to(u.deg)

        return SkyCoord(ra, dec), angle


208 209
class CatalogShifter(object):
    def __init__(self, catalog, mapshape, mapwcs):
LUSTIG Peter's avatar
LUSTIG Peter committed
210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228
        '''
        Class that takes source 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

        Returns
        -------
        CatalogShifter object that can be called to create shifted catalogs
        '''
229 230 231

        # create centerframe
        centerpixel = (np.array(nm.shape)/2).astype(int)
232 233
        centerworld = nm.wcs.wcs_pix2world(*centerpixel, 0)
        centerframe = SkyCoord(*(centerworld*u.deg)).skyoffset_frame()
234 235 236 237 238 239 240

        catcoords = SkyCoord(catalog['ra'], catalog['dec'], unit=u.deg)

        self.catalog = catalog
        self.catcoords = catcoords
        self.centerframe = centerframe

241
    def __call__(self, coord_in_center, rotation=0*u.deg, transform=None):
LUSTIG Peter's avatar
LUSTIG Peter committed
242 243 244 245 246 247 248 249 250 251
        '''
        Parameters
        ----------
        coord_in_center : :class:`astropy.coordinates.SkyCoord`
            coordinate that is shifted in the center of the map
        rotation : :class:`astropy.units.Quantity`
            catalog is rotated through this angle
        transform: callable or None.
            If not None, takes coord_in_center and rotation as argument
            and transforms them to a :class:`astropy.coordinates.SkyCoord`
LUSTIG Peter's avatar
LUSTIG Peter committed
252 253
            and an angle that are used as coord_in_center and rotation to
            perform transformation.
LUSTIG Peter's avatar
LUSTIG Peter committed
254 255 256 257 258
        Returns
        -------
        catalog : :class:``astropy.table.Table``
            catalog containing new source positions
        '''
259 260 261 262
        catalog = self.catalog
        catcoords = self.catcoords
        centerframe = self.centerframe

263 264 265
        if transform is not None:
            coord_in_center, rotation = transform(coord_in_center, rotation)

266 267
        # create new frame with new center and transform soucres into
        sourcecenterframe = coord_in_center.skyoffset_frame(rotation=rotation)
268
        sourcecoordinates = catcoords.transform_to(sourcecenterframe)
269 270

        # change variable in center with variable in mapcenter
271 272
        relativecoordinates = SkyCoord(sourcecoordinates.lon.deg * u.deg,
                                       sourcecoordinates.lat.deg * u.deg,
273 274 275 276
                                       frame=centerframe)
        absolutecoordinates = relativecoordinates.transform_to('icrs')

        # write in table copy
277
        ra, dec = absolutecoordinates.ra.deg, absolutecoordinates.dec.deg
278 279 280 281
        catalog['ra'] = ra
        catalog['dec'] = dec
        return catalog

282

283
if __name__ == "__main__":
284 285
    # test_catalogshifter()

286
    nmfile = "/home/peter/Dokumente/Uni/Paris/Stage/data/map.fits"
287 288
    nmfile = ("/home/peter/Dokumente/Uni/Paris/Stage/N2CLS/ForMarseille/"
              "HLS091828_common_mode_kids_out/map_JK.fits")
289
    tablefile = "/home/peter/Dokumente/Uni/Paris/Stage/N2CLS/SIDES_NIKA.fits"
290 291 292 293 294 295 296
    catalog = Table.read(tablefile)
    catalog['ra'].unit = u.deg
    catalog['dec'].unit = u.deg
    tr = transUniformCoords(catalog)
    tr(0.5,.5,.5)

    #%%
297 298 299
    nm = NikaMap.read(nmfile)
    wcs = nm.wcs
    shape = nm.shape
300

301 302
    catalog = catalog[catalog['SNIKA1200'] > 1E-3]

303 304 305 306 307 308
    sh = CatalogShifter(catalog, (501, 501), nm.wcs)
    catalog = sh(SkyCoord(1.4*u.deg, 1.4*u.deg))
    cat = CatalogFor_pos_cat(catalog, fluxkey='SNIKA1200',
                             fluxunit=u.Jy, minflux=1*u.mJy)
    nm.add_gaussian_sources(cat_gen=pos_cat, catalog=cat, wcs=nm.wcs)
    nm.plot()
309
    print('done')
310

311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333
    #%%
    '''
    nmcenterpixel = (np.array(nm.shape)/2).astype(int)
    nmcenterworld = nm.wcs.wcs_pix2world(*nmcenterpixel, 0)
    nmcenterframe = SkyCoord(*(nmcenterworld*u.deg)).skyoffset_frame()
    catcoords = SkyCoord(catalog['ra']*u.deg, catalog['dec']*u.deg)
    # coord_to_center = catcoords[0]


    rotation=45*u.deg
    coord_in_center = SkyCoord(1.3*u.deg, 1.4*u.deg)
    newcenterframe = coord_in_center.skyoffset_frame(rotation=rotation)
    newcoordinates = catcoords.transform_to(newcenterframe)
    newcoordinates
    newnew = SkyCoord(newcoordinates.lon.deg * u.deg,
                      newcoordinates.lat.deg * u.deg, frame=nmcenterframe)
    newnew = newnew.transform_to('icrs')
    ra, dec = newnew.ra.deg, newnew.dec.deg
    catalog['_ra'] = ra
    catalog['_dec'] = dec
    cat = catalog
    cat = CatalogFor_pos_cat(catalog, fluxkey='SNIKA1200',
                             fluxunit=u.Jy, minflux=1*u.mJy)
334

335 336 337 338 339 340
    nm.add_gaussian_sources(cat_gen=pos_cat, catalog=cat, wcs=nm.wcs)

    nm.plot()

    '''
    '''
341 342 343 344 345
    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)
346
    '''
347
    # %matplotlib tk