add_catalog_sources.py 10.9 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
252
253
254
255
256
257
        '''
        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`
            and an angle that are used as coord_in_center and rotation.
        Returns
        -------
        catalog : :class:``astropy.table.Table``
            catalog containing new source positions
        '''
258
259
260
261
        catalog = self.catalog
        catcoords = self.catcoords
        centerframe = self.centerframe

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

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

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

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

281

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

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

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

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

302
303
304
305
306
307
    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()
308
    print('done')
309

310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
    #%%
    '''
    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)
333

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

    nm.plot()

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