nezaval/nsd
1
0
Fork 0
forked from erosita/uds
Code Pull Requests Activity

cross match

Browse Source
This commit is contained in:
Roman Krivonos
2023-03-07 16:01:23 +03:00
parent bfb7f7dc21
commit ae61b520cf
56 changed files with 63153 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

157
uds/uds/utils.py
View File

@@ -254,3 +254,160 @@ def test_exe(program):
print("\n*** Command {} not found ***\n".format(program))
sys.exit()
def save_ds9reg(filename,scale=60):
if(os.path.isfile(filename)==True):
fout=filename.replace(".fits", ".reg")
hdul = fits.open(filename)
tbdata = hdul[1].data
with open(fout, 'w') as writer:
for rec in tbdata:
writer.write("fk5;circle({}, {}, {})\n".format(rec['ra'],rec['dec'],rec['radec_err']/scale))
def save_ermldet_ds9reg(filename,scale=60,id_instr=1,id_band=1,ext_like=0.0,label='det_like'):
if(os.path.isfile(filename)==True):
#fout=filename.replace(".fits", ".instr{}.{}.reg".format(id_instr,label))
fout=filename.replace(".fits", ".{}.reg".format(label))
hdul = fits.open(filename)
tbdata = hdul[1].data
with open(fout, 'w') as writer:
for rec in tbdata:
if not (rec['id_inst'] == id_instr and rec['id_band'] == id_band):
continue
if (rec['ext_like'] > ext_like):
continue
if(abs(rec[label])>0.0):
writer.write("fk5;circle({}, {}, {}) # text={{{:.1f}}}\n".format(rec['ra'],rec['dec'],rec['radec_err']/scale,rec[label]))
def save_catprep_ds9reg(filename,scale=60,id_instr=1,id_band=1,ext_like=0.0,label='det_like_0'):
if(os.path.isfile(filename)==True):
#fout=filename.replace(".fits", ".instr{}.{}.reg".format(id_instr,label))
fout=filename.replace(".fits", ".{}.reg".format(label))
hdul = fits.open(filename)
tbdata = hdul[1].data
with open(fout, 'w') as writer:
for rec in tbdata:
if (rec['ext_like'] > ext_like):
continue
if(abs(rec[label])>0.0):
writer.write("fk5;circle({}, {}, {}) # text={{{:.1f}}}\n".format(rec['ra'],rec['dec'],rec['radec_err']/scale,rec[label]))
def crossmatch_shu2019(filename,scale=1200,devmax=30,w=None,naxis1=1,naxis2=1,dlmin=6.0,dlmax=10000,ext_like=0.0,outkey='shu2019', catalog=None):
if(os.path.isfile(filename)==False):
print("File not found {}".format(filename))
print("Start cross-match with Gaia-unWISE")
if not catalog:
print("ERROR: Please provide catalog.")
return
hdul = fits.open(catalog)
tbdata_ref = hdul[1].data
cat=[]
for rec in tbdata_ref:
cat.append({"ra":rec['ra'],"dec":rec['dec'],"gaia":rec["gaia_sourceid"],"unwise":rec["unwise_objid"]})
hdul = fits.open(filename)
tbdata_src = hdul[1].data
srcs=[]
for rec in tbdata_src:
if (rec['ext_like'] > ext_like):
continue
srcs.append({"id":rec['id_src'],
"ra":rec['ra'],
"det_like":rec['det_like_0'],
"ext_like":rec['ext_like'],
"dec":rec['dec'],
"radec_err":rec['radec_err']/scale})
fout=filename.replace(".fits", ".{}.cross".format(outkey))
if(os.path.isfile(fout)==True):
os.remove(fout)
t = Table(names=('ID','CNT',
'DET_LIKE', 'EXT_LIKE',
'SRC_RA', 'SRC_DEC',
'REF_RA', 'REF_DEC',
'SEP','GAIA','UNWISE'),
dtype=('i4','i4','f4', 'f4','f8', 'f8','f8', 'f8','f4','i8','S16'))
d = {}
skip = {}
# do correlation
for src in srcs:
src_crd = SkyCoord(src['ra'], src['dec'], frame=FK5(), unit="deg")
sep=0
for scat in cat:
scat_crd = SkyCoord(scat['ra'], scat['dec'], frame=FK5(), unit="deg")
sep = src_crd.separation(scat_crd).arcsec
if(sep < devmax):
sid=src['id']
if(sid in d.keys()):
d[sid]=d[sid]+1
skip[sid]=1
else:
d[sid]=1
t.add_row((src['id'], d[sid], src['det_like'], src['ext_like'],
src['ra'], src['dec'],
scat['ra'], scat['dec'],
sep, scat['gaia'], scat['unwise']))
t.write(fout, format='fits')
src_tab = Table(names=('ID','RA', 'DEC'), dtype=('i4','f8','f8',))
src_fout=fout.replace(".cross", ".{}.src".format(outkey))
ref_fout=fout.replace(".cross", ".{}.ref".format(outkey))
ref_tab = Table(names=('GAIA','RA', 'DEC'), dtype=('i8','f8','f8',))
hdul = fits.open(fout)
tbdata = hdul[1].data
srcs=[]
src_map=np.zeros((naxis1, naxis2))
ref_map=np.zeros((naxis1, naxis2))
delta_ra = []
delta_dec = []
for rec in tbdata:
#if (rec['det_like'] < 10.0):
# continue
if (rec['id'] in skip.keys()):
print("Skip ID={} DET_LIKE={:.2f}".format(rec['id'],rec['det_like']))
else:
src_crd = SkyCoord(rec['src_ra'], rec['src_dec'], frame=FK5(), unit="deg")
ref_crd = SkyCoord(rec['ref_ra'], rec['ref_dec'], frame=FK5(), unit="deg")
src_x, src_y = np.round(w.world_to_pixel(src_crd))
src_map[int(src_x), int(src_y)]=1
ref_x, ref_y = np.round(w.world_to_pixel(ref_crd))
ref_map[int(ref_x), int(ref_y)]=1
src_tab.add_row((rec['id'],rec['src_ra'],rec['src_dec']))
ref_tab.add_row((rec['gaia'],rec['ref_ra'],rec['ref_dec']))
delta_ra.append(rec['src_ra']-rec['ref_ra'])
delta_dec.append(rec['src_dec']-rec['ref_dec'])
print("delta RA {:.6f} delta DEC {:.6f}".format(rec['src_ra']-rec['ref_ra'],
rec['src_dec']-rec['ref_dec']))
print("Mean delta RA: {}, delta Dec: {} (arcsec)".format(statistics.mean(delta_ra)*3600,statistics.mean(delta_dec)*3600))
src_tab.write(src_fout, format='fits', overwrite=True)
with fits.open(src_fout) as f:
f[0].verify('fix')
f[1].verify('fix')
with fits.open(src_fout, 'update') as f:
f[0].header=f[0].header+w.to_header()
f[0].data=src_map
#f[1].header=f[1].header+w.to_header()
ref_tab.write(ref_fout, format='fits', overwrite=True)
with fits.open(ref_fout) as f:
f[0].verify('fix')
f[1].verify('fix')
with fits.open(ref_fout, 'update') as f:
f[0].header=f[0].header+w.to_header()
f[0].data=ref_map
#f[1].header=f[1].header+w.to_header()