Commit ee74c132 authored by LUSTIG Peter's avatar LUSTIG Peter

further improvements, added explications

parent a2c391a9
......@@ -2,6 +2,12 @@ from nikamap import NikaMap
from astropy.table import Table
import numpy as np
import astropy.units as u
from nikamap.utils import pos_cat
import matplotlib.pyplot as plt
from astropy.wcs import WCS
#%load_ext autoreload
#%matplotlib tk
def GetInRange(array, low=None, high=None, extra=None):
......@@ -42,23 +48,69 @@ def Shift(shiftrange, relcoord):
class CatalogShifter:
def __init__(self, catalog, mapshape, wcs):
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
maplimits = wcs.wcs_pix2world(*(np.array(((0, 0), mapshape)).T), 0)
maplimits = np.sort(maplimits, axis=-1)
ra, dec = catalog['ra'], catalog['dec']
sourcelimits = np.array([[np.amax(ra), np.amin(ra)],
# borders of source area
sourcelimits = np.array([[np.amin(ra), np.amax(ra)],
[np.amin(dec), np.amax(dec)]])
shiftlimits = maplimits - sourcelimits
# 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))
self.catalog = catalog
self._maplmits = maplimits
self._shiftlimits = shiftlimits
def __call__(self, relra=None, reldec=None):
'''
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
'''
ra_shiftlimits, dec_shiftlimits = self._shiftlimits
maplimits = self._maplmits
outcat = self.catalog # this should create a copy
if relra is None:
......@@ -66,20 +118,59 @@ class CatalogShifter:
if reldec is None:
reldec = np.random.random_sample()
rashift = Shift(np.sort(ra_shiftlimits), relra)
decshift = Shift(np.sort(dec_shiftlimits), reldec)
rashift = Shift(ra_shiftlimits, relra)
decshift = Shift(dec_shiftlimits, reldec)
outcat['ra'] = outcat['ra'] + rashift
outcat['dec'] = outcat['dec'] + decshift
_, outcat = GetInRange(outcat['ra'], *np.sort(maplimits[0]), outcat)
_, outcat = GetInRange(outcat['dec'], *np.sort(maplimits[1]), outcat)
return outcat
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))
if __name__ == "__main__":
# test_catalogshifter()
nmfile = "/home/peter/Dokumente/Uni/Paris/Stage/data/map.fits"
nmfile = ("/home/peter/Dokumente/Uni/Paris/Stage/N2CLS/ForMarseille/"
"HLS091828_common_mode_kids_out/map_JK.fits")
tablefile = "/home/peter/Dokumente/Uni/Paris/Stage/N2CLS/SIDES_NIKA.fits"
nm = NikaMap.read(nmfile)
wcs = nm.wcs
......@@ -87,7 +178,12 @@ if __name__ == "__main__":
catalog = Table.read(tablefile)
catalog.keys()
sh = CatalogShifter(catalog, shape, wcs)
sh()
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)
# %matplotlib tk
nm.plot()
plt.show(block=True)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment