Commit 3475b41f authored by LUSTIG Peter's avatar LUSTIG Peter

should properly do transformation now

parent f2b534e4
......@@ -6,10 +6,11 @@ from nikamap.utils import pos_cat
import matplotlib.pyplot as plt
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
'''
%load_ext autoreload
%autoreload 2
%matplotlib tk
'''
def GetInRange(array, low=None, high=None, extra=None):
......@@ -166,13 +167,50 @@ def test_catalogshifter():
#print(sh(0,0, coverage=1))
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
def __call__(self, relra, reldec, angle):
'''
Parameters
----------
Coordinate: float between 0 and 1
'''
ralimits = self.ralimits
decuniformlimits = self.decuniformlimits
ra = rel_to_coord(*ralimits, relra) * u.deg
angle = angle * 360 * u.deg
decuniform_inrange = rel_to_coord(*decuniformlimits, reldec)
# 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
class CatalogShifter(object):
def __init__(self, catalog, mapshape, mapwcs):
# create centerframe
centerpixel = (np.array(nm.shape)/2).astype(int)
centerworld = nm.wcs.wcs_pix2world(*nmcenterpixel, 0)
centerframe = SkyCoord(*(nmcenterworld*u.deg)).skyoffset_frame()
centerworld = nm.wcs.wcs_pix2world(*centerpixel, 0)
centerframe = SkyCoord(*(centerworld*u.deg)).skyoffset_frame()
catcoords = SkyCoord(catalog['ra'], catalog['dec'], unit=u.deg)
......@@ -180,27 +218,31 @@ class CatalogShifter(object):
self.catcoords = catcoords
self.centerframe = centerframe
def __call__(self, coord_in_center, rotation=0*u.deg):
def __call__(self, coord_in_center, rotation=0*u.deg, transform=None):
catalog = self.catalog
catcoords = self.catcoords
centerframe = self.centerframe
if transform is not None:
coord_in_center, rotation = transform(coord_in_center, rotation)
# create new frame with new center and transform soucres into
sourcecenterframe = coord_in_center.skyoffset_frame(rotation=rotation)
sourcecoordinates = catcoords.transform_to(newcenterframe)
sourcecoordinates = catcoords.transform_to(sourcecenterframe)
# change variable in center with variable in mapcenter
relativecoordinates = SkyCoord(newcoordinates.lon.deg * u.deg,
newcoordinates.lat.deg * u.deg,
relativecoordinates = SkyCoord(sourcecoordinates.lon.deg * u.deg,
sourcecoordinates.lat.deg * u.deg,
frame=centerframe)
absolutecoordinates = relativecoordinates.transform_to('icrs')
# write in table copy
ra, dec = newnew.ra.deg, newnew.dec.deg
ra, dec = absolutecoordinates.ra.deg, absolutecoordinates.dec.deg
catalog['ra'] = ra
catalog['dec'] = dec
return catalog
if __name__ == "__main__":
# test_catalogshifter()
......@@ -208,16 +250,27 @@ if __name__ == "__main__":
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"
catalog = Table.read(tablefile)
catalog['ra'].unit = u.deg
catalog['dec'].unit = u.deg
tr = transUniformCoords(catalog)
tr(0.5,.5,.5)
#%%
nm = NikaMap.read(nmfile)
wcs = nm.wcs
shape = nm.shape
catalog = Table.read(tablefile)
catalog = catalog[catalog['SNIKA1200'] > 1E-3]
sh = NewCatalogShifter(catalog, (501, 501), nm.wcs)
t = sh(SkyCoord(1.4*u.deg, 1.4*u.deg))
print(t)
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()
print('done')
#%%
'''
nmcenterpixel = (np.array(nm.shape)/2).astype(int)
......
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