import glob from astropy.io import fits from astropy.table import Table from astropy.time import Time import ntpath import os import numpy as np from astropy.coordinates import SkyCoord from astropy import units as u from parse_log import get_teltime import pyvo as vo def mk_table(listfiles): path_arr = [] obj_arr = [] ra_arr = [] dec_arr = [] time_arr = [] timeall = [] dat_arr = [] am_arr = [] exptime_arr = [] for i, fl in enumerate(listfiles): frame = fits.open(fl) hdr = frame[0].header if (type(hdr['OBJECT']) is str): path_arr.append(fl) obj_arr.append(hdr['OBJECT']) ra_arr.append(hdr['RA-D']) dec_arr.append(hdr['DEC-D']) time_arr.append(hdr['UT-TIME']) dat_arr.append(hdr['UT-DATE']) exptime_arr.append(hdr['EXPTIME']) am_arr.append(hdr['AIRMASS']) timeall.append(Time(hdr['UT-DATE']+'T'+hdr['UT-TIME'], format='isot', scale='utc').jd) frame.close() t = Table([path_arr, obj_arr, ra_arr, dec_arr, timeall, time_arr, dat_arr, exptime_arr, am_arr], names=('PATH', 'OBJECT','RA', 'DEC','TIMEUTC', 'TIME', 'DATE', 'EXPTIME', 'AIRMASS')) t.sort(['DATE', 'TIME']) return t service = vo.dal.TAPService("http://tapvizier.u-strasbg.fr/TAPVizieR/tap") indir = '/data1/Data/MagE/ut181214_15/' flatsdir = '/data1/Data/MagE/ut181214_15/' logfile = '/data1/Data/MagE/ut181214_15/obs_log_20181214_mage.txt' telt = get_teltime(logfile) outdir = 'processed_'+indir.split('/')[-2] confdir = 'proc_'+indir.split('/')[-2] proc_dir = '/u/kirg/mage_proc/' full_path_conf = os.path.join(proc_dir, confdir) if confdir+'/' not in glob.glob('*/'): os.system('mkdir '+confdir) all_obs_root = os.path.abspath(os.path.join(os.path.dirname(indir), '..')) infiles = glob.glob(indir+'mage????.fits') flatfiles = glob.glob(flatsdir+'mage????.fits') all_dirs = glob.glob(os.path.abspath(os.path.join(all_obs_root, '*'))) os.chdir(all_obs_root) if outdir+'/' not in glob.glob('*/'): os.system('mkdir '+outdir) outpath = os.path.join(all_obs_root, outdir) os.chdir(proc_dir) t = mk_table(infiles) t_flats = mk_table(flatfiles) xe_flash_files = Table(dtype=t_flats.dtype) flat_files = Table(dtype=t_flats.dtype) obj_files = Table(dtype=t.dtype) sky_files = Table(dtype=t.dtype) thar_files = Table(dtype=t.dtype) tel_files = Table(dtype=t.dtype) for row in t_flats: if (('xe' in row['OBJECT'].lower()) & ('fl' in row['OBJECT'].lower())): xe_flash_files.add_row(row) if ((('flat' in row['OBJECT'].lower()) | ('quartz' in row['OBJECT'].lower())) & (('skyflat' not in row['OBJECT']) & ('SkyFlat' not in row['OBJECT']) & ('sky' not in row['OBJECT']))): flat_files.add_row(row) for row in t: if row['EXPTIME']>800: obj_files.add_row(row) if (('_sky' in row['OBJECT']) | ('-sky' in row['OBJECT']) | (' sky' in row['OBJECT']) | (row['EXPTIME'] == 300.0)): sky_files.add_row(row) if (('ThAr' in row['OBJECT']) | ('arc' in row['OBJECT'])): thar_files.add_row(row) if ((('HIP' in row['OBJECT']) | ('HR' in row['OBJECT'])) & (row['EXPTIME']>10)): tel_files.add_row(row) tel_stars = tel_files.copy() tel_files_c = SkyCoord(ra=tel_files['RA']*u.degree, dec=tel_files['DEC']*u.degree) sky_files_c = SkyCoord(ra=sky_files['RA']*u.degree, dec=sky_files['DEC']*u.degree) tel_stars.sort('RA') telstars = Table(dtype=t.dtype) for i,frame in enumerate(tel_stars): if i == 0: sc1 = SkyCoord(frame['RA']*u.deg,frame['DEC']*u.deg) telstars.add_row(frame) else: sc2 = SkyCoord(frame['RA']*u.deg,frame['DEC']*u.deg) sep = sc2.separation(sc1) if (sep.value > 2./60.): telstars.add_row(frame) sc1 = sc2 nexp_objs = {} xefl0 = ntpath.basename(xe_flash_files[0]['PATH']) print telstars tst_coords = SkyCoord(ra=telstars['RA']*u.degree, dec=telstars['DEC']*u.degree) xefl0 = int(xefl0.replace('mage','').replace('.fits', '')) for obj_frame in obj_files: time_frame = obj_frame['TIMEUTC'] c = SkyCoord(ra=obj_frame['RA']*u.degree, dec=obj_frame['DEC']*u.degree) idx, d2d, d3d = c.match_to_catalog_sky(tst_coords) frms_sep = tst_coords[idx].separation(tel_files_c) alltelfrms = tel_files[(frms_sep.value < 2.0/60.0)] alltelfrms['dt'] = np.abs(alltelfrms['TIMEUTC']-time_frame) alltelfrms.sort('dt') selectedtel = alltelfrms[0] query_hip = "SELECT * FROM \"I/239/hip_main\" WHERE CONTAINS(POINT('ICRS',RAICRS,DEICRS),CIRCLE('ICRS',%.5f,%.5f,%.5f))=1" % (selectedtel['RA'], selectedtel['DEC'], 1./60.) hip_cat = service.search(query_hip) Vmag = hip_cat[0]['Vmag'] idx_sky, d2d_sky, d3d = c.match_to_catalog_sky(sky_files_c) if d2d_sky.value < 1./60.: sub_sky = 1 else: sub_sky = 0 obj_name = obj_frame['OBJECT'] am = obj_frame['AIRMASS'] thar_obj = thar_files[(thar_files['TIMEUTC'] - time_frame)>0][0:3] if sub_sky == 1: skyflsarr = [] all_obj_sky = Table(dtype=t.dtype) for j in range(len(sky_files['OBJECT'])): if c.separation(sky_files_c[j]).value<4./60.: skyflsarr.append(ntpath.basename(sky_files[j]['PATH'])) all_obj_sky.add_row(sky_files[j]) skyflsstr = str(skyflsarr) sf_0 = all_obj_sky[(all_obj_sky['TIMEUTC'] - time_frame)<0] sf_1 = all_obj_sky[(all_obj_sky['TIMEUTC'] - time_frame)>0] if len(skyflsarr)>=2: if len(sf_0) == 0: sky1 = all_obj_sky[0] sky2 = all_obj_sky[1] if len(sf_1) == 0: sky1 = all_obj_sky[-1] sky2 = all_obj_sky[-2] if ((len(sf_0)>0) & (len(sf_1)>0)): sky1 = sf_0[-1] sky2 = sf_1[0] skyflsstr = str([ntpath.basename(sky1['PATH']), ntpath.basename(sky2['PATH'])]) else: skyflsstr = str(skyflsarr) tel_files['Dam'] = np.abs(tel_files['AIRMASS'] - am) #tel_files.sort('Dam') tellurics = tel_files[0:2].copy() tellurics.sort('TIMEUTC') tel_time = selectedtel['TIMEUTC'] tel_thar_obj = thar_files[(thar_files['TIMEUTC'] - tel_time)>0][0:3] if obj_name not in nexp_objs.keys(): nexp_objs[obj_name] = 1 else: nexp_objs[obj_name] = nexp_objs[obj_name] + 1 confname = obj_name + '_exp'+str(nexp_objs[obj_name])+'.pro' outname = "%s_%03d" % (obj_name, nexp_objs[obj_name]) flat0 = ntpath.basename(flat_files[0]['PATH']) flat0 = int(flat0.replace('mage','').replace('.fits', '')) f = open(os.path.join(full_path_conf, confname),"w+") f.write("data_path='%s'\r\n" % indir) f.write("calib_path='%s'\r\n" % flatsdir) f.write("outp_path='%s/'\r\n" % outpath) f.write("outname = '%s'\r\n" % outname) f.write("telname = '%s'\r\n" % selectedtel['OBJECT']) f.write("\r\n") f.write("conf_str={ $\r\n") f.write(" flat_file:calib_path+'mage'+string(%i+findgen(%i),format='(i4.4)')+'.fits' ,$\r\n" % (flat0, len(flat_files))) f.write(" xe_file:calib_path+'mage'+string(%i+findgen(%i),format='(i4.4)')+'.fits',$\r\n" % (xefl0,len(xe_flash_files))) f.write(" arc_file:data_path+['%s','%s','%s'],$\r\n" % (ntpath.basename(thar_obj[0]['PATH']), ntpath.basename(thar_obj[1]['PATH']),ntpath.basename(thar_obj[2]['PATH']))) f.write(" obj_file:data_path+'%s',$\r\n" % (ntpath.basename(obj_frame['PATH']))) f.write(" cr_reject:1,$ ;;; cosmic ray hit rejection using LACosmic\r\n") f.write(" sky_mask:[0,dblarr(10)+1,0],$ ;;; sky slitmask bottom-to-top: 0=exclude; 1=include\r\n") if sub_sky == 1: f.write(" sky_file:data_path+%s,$ ;;; offset sky exposure (here the same as obj_file, initialized later)\r\n" % (skyflsstr)) f.write(" sky_normexp:-1.0,$ ;;; sky exposure normalization factor\r\n") f.write(" sub_sky:%i,$ ;;; sky subtraction (1=Yes/0=No)\r\n" % sub_sky) f.write(" corr_tel:1,$ ;;; perform telluric correction (1=Yes/0=No)\r\n") f.write(" flux_corr:1,$ ;;; perform flux calibration (1=Yes/0=No, corr_tel must be Yes for Yes)\r\n") f.write(" outfile_merged: outp_path+outname+'_merged.fits',$ ;;; merged flux calibrated output (no telluric correction)\r\n") f.write(" outfile_merged_corr: outp_path+outname+'_merged_corr.fits',$ ;;; merged telluric corrected flux calibrated output\r\n") f.write(" raw_2d_output:0, $ ;;; output raw counts (1=Yes/0=No)\r\n") f.write(" outfile_orders: outp_path+outname+'_test_2d.fits' $ ;;; unmerged flux calibrated output (no telluric correction)\r\n") f.write("}\r\n") f.write("tel_conf_str={ $\r\n") f.write(" tel_file:data_path+'%s',$ ;; telluric star %s\r\n" % (ntpath.basename(selectedtel['PATH']),selectedtel['OBJECT'])) f.write(" tel_arc_file:data_path+['%s','%s','%s'],$ ;;; telluric specific arc frames\r\n" % (ntpath.basename(tel_thar_obj[0]['PATH']), ntpath.basename(tel_thar_obj[1]['PATH']),ntpath.basename(tel_thar_obj[2]['PATH']))) f.write(" tel_arc_sci:0,$ ;;; use science arc frames (1=Yes/0=No)\r\n") f.write(" tel_sci_frame:0,$ ;;; use science frame as telluric (1=Yes/0=No)\r\n") f.write(" sky_mask_tel:-1,$ ;;; telluric sky subtraction mask, same as on science frame\r\n") f.write(" sub_sky_tel:1,$ ;;; subtract sky from telluric (1=Yes/0=No)\r\n") f.write(" tel_slit_reg:3+indgen(32),$ ;;; slit goes from 0 to 36 pix, bottom-to-top\r\n") f.write(" tel_texp:%.2f,$ ;;; telluric exposure time on the slit\r\n" % telt[ntpath.basename(selectedtel['PATH'])]) f.write(" tel_magv:%.2f,$ ;;; telluric V band magnitude (assuming an A0V spectum)\r\n" % Vmag) f.write(" tel_vr0: 0.0,$ ;;; radial velocity of the telluric star (initial guess)\r\n") f.write(" tel_vsini0:100.0,$ ;;; rotational velocity v*sin(i) of the telluric star (initial guess)\r\n") f.write(" tel_spec_list:file_search(getenv('MAGE_PIPELINE_PATH')+'calib_MagE/stellar_templates/phoenix/lte09*-0.0.P*'),$ ;;; telluric star template(s)\r\n") f.write(" tel_novsini:0,$ ;;; skip fitting for v*sin(i) of the telluric star (0=Fitting/1=Fix)\r\n") f.write(" tel_gausslosvd:0,$ ;;; LOSVD shape for the telluric star (0=Rot/1=Gaussian); for Gaussian vsini becomes sigma\r\n") f.write(" tel_output:1,$\r\n") f.write(" outfile_merged_telluric: outp_path+telname+'_%s_merged_corr.fits'$\r\n" % (obj_frame['DATE'].replace('-',''))) f.write("}\r\nrun_mage_pipeline,conf_str,tel_conf_str=tel_conf_str\r\nend\r\n") f.close()