import os import sys import numpy as np from astropy.io import fits from astropy import wcs from astropy.wcs import WCS from astropy.io.fits import update from astropy.io.fits import getdata import glob from fitsio import FITS from pathlib import Path from astropy.table import QTable, Table, Column from astropy import units as u from astropy.coordinates import SkyCoord # High-level coordinates from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames from astropy.coordinates import Angle, Latitude, Longitude # Angles from astropy.time import Time, TimeDelta, TimezoneInfo, TimeFromEpoch import statistics import shutil MJDREF = 51543.875 TZ_UTC = TimezoneInfo(utc_offset=0*u.hour) TZ_MSK = TimezoneInfo(utc_offset=3*u.hour) TIMESTR_FORMAT = '%Y.%m.%d %H:%M:%S' def mission2date_utc(timesec): mjdref = Time(MJDREF, format='mjd') delta = TimeDelta(timesec, format='sec') dtime = mjdref + delta return dtime.to_datetime(timezone=TZ_UTC) def mission2date_msk(timesec): mjdref = Time(MJDREF, format='mjd') delta = TimeDelta(timesec, format='sec') dtime = mjdref + delta + 3*u.hour return dtime.to_datetime(timezone=TZ_MSK) def create_folder(folder): if not (os.path.exists(folder)): os.makedirs(folder) def remove_file(filename): if(os.path.isfile(filename)==True): os.remove(filename) def do_evtool_esass(evfile=None,events=None,outfile=None,evlist=None, gti=None,region=None,emin=None,emax=None, rmlock=False, do_center=False, ra_cen=None, de_cen=None): eventfiles=None if(events): eventfiles="eventfiles=\'{}\'".format((" ").join(events)) if(evlist): eventfiles="eventfiles=@{}".format(evlist) if(evfile): eventfiles="eventfiles={}".format(evfile) if not (eventfiles): print("ERROR: Event files not provided") if(do_center==True and evfile): """ re-center event file """ if not (ra_cen and de_cen): print("Please provide center coordinates") sys.exit() test_exe('radec2xy') cmd=["radec2xy", "{}".format(evfile), "{}".format(ra_cen), "{}".format(de_cen)] print((" ").join(cmd)) os.system((" ").join(cmd)) emin="emin={}".format(emin) if(emin) else '' emax="emax={}".format(emax) if(emax) else '' gti="gti=\'{}\'".format(gti) if(gti) else '' region="region=\'{}\'".format(region) if(region) else '' cmd=["evtool", eventfiles, gti,region,emin,emax, "outfile={}".format(outfile), "image=yes", "flag=0x2000", "pattern=15" ] print((" ").join(cmd)) test_exe('evtool') os.system((" ").join(cmd)) if(rmlock==True): lockfiles=outfile.replace(".fits", "*.lock") for filename in glob.glob(lockfiles): print("Removing {}".format(filename)) if(os.path.isfile(filename)==True): os.remove(filename) pass def set_bit(value, bit): return value | (1<250 creates new file with added _badpix.fits """ nraw=384 skip_after=250 nskip=(nraw - skip_after) nrows0=nraw*nskip fout=filename.replace(".fits", "_badpix.fits") shutil.copyfile(filename, fout) hdul = fits.open(fout) #hdul.info() tstart=hdul[1].header['TSTART'] tstop=hdul[1].header['TSTOP'] nrows=hdul[1].header['NAXIS2'] f = FITS(fout,'rw') data = f['EVENTS'].read() #print(data.dtype) for row in range(nrows): if (data['RAWX'][row] >= skip_after): data['FLAG'][row]=set_bit(data['FLAG'][row],13) #hdr = fits[0].read_header() f[1].write(data) f[1].write_comment("Changed by Roman Krivonos") f[1].write_history("All RAWx>=250 marked as BAD 0x2000") data_badpix = f['BADPIX6'].read() data4 = np.zeros(nrows0,dtype=data_badpix.dtype) n=0 for i in range(skip_after,nraw): for j in range(nraw): data4['RAWX'][n]=i data4['RAWY'][n]=j data4['TIMEMIN'][n]=tstart data4['TIMEMAX'][n]=tstop data4['BADFLAG'][n]=3 data4['TYPE'][n]=3 data4['PHAMAX'][n]=12000 data4['PHAMIN'][n]=0 n=n+1 f[4].append(data4) f.close() def init_events(key=None, eband_selected=[0], eband_index=None, ra_cen=None, de_cen=None,do_init=True,vign=True, emin_kev=None, emax_kev=None, infile_dir=None, outfile_dir=None, do_obsmode=False,do_center=False,do_evtool=False,do_expmap=False): expmaps=[] expmap_all=[] infile_expmap=[] emin=[] emax=[] emin_ev_str=[] emax_ev_str=[] outfile_evtool=[] outfile_expmap=[] outfile_post='.fits' if not (infile_dir): print("infile_dir=?") return if not (outfile_dir): print("outfile_dir=?") return print("init events -- reading key {}".format(key)) infile="{}/{}.fits".format(infile_dir,key) if(do_obsmode==True and do_init==True): """ correct OBS_MODE in files """ lockfile="{}/{}.obsmode.lock".format(infile_dir,key) if not os.path.isfile(lockfile): with fits.open(infile) as hdu: for h in hdu: extname = h.header['EXTNAME'] if "EXTNAME" in h.header else "None" changeto = 'SURVEY' if ('scan' in infile) else 'POINTING' if not (h.header['OBS_MODE'] == changeto): print("*** {} {} --> {}".format(extname,h.header['OBS_MODE'],changeto)) h.header['OBS_MODE'] = changeto hdu.writeto(infile,overwrite=True) Path(lockfile).touch() else: print("Lock file {} is found, skipping OBS_MODE correction.".format(lockfile)) pass if(do_center==True and do_init==True): """ re-center original events files """ if not (ra_cen and de_cen): print("Please provide center coordinates") sys.exit() cmd=["radec2xy", "{}".format(infile), "{}".format(ra_cen), "{}".format(de_cen)] lockfile="{}/{}.radec2xy.lock".format(infile_dir,key) if not os.path.isfile(lockfile): print((" ").join(cmd)) test_exe('radec2xy') os.system((" ").join(cmd)) Path(lockfile).touch() else: print("Lock file {} is found, skipping radec2xy.".format(lockfile)) pass outfile_evtool="{}/{}_EventList_en{}{}".format(outfile_dir, key, eband_index, outfile_post) cmd=["evtool", "eventfiles=%s" %(infile), "outfile=%s" %(outfile_evtool), "emin=%f" %(emin_kev), "emax=%f" %(emax_kev), #"region=%s" %(region), "image=yes", "flag=0x2000", "pattern=15" ] # run the command if(do_evtool==True and do_init==True): #log = subprocess.check_call(cmd) print((" ").join(cmd)) test_exe('evtool') os.system((" ").join(cmd)) """ correct OBS_MODE """ with fits.open(outfile_evtool) as hdu: for h in hdu: extname = h.header['EXTNAME'] if "EXTNAME" in h.header else "None" changeto = 'SURVEY' if ('scan' in infile) else 'POINTING' if not (h.header['OBS_MODE'] == changeto): print("*** {} {} --> {}".format(extname,h.header['OBS_MODE'],changeto)) h.header['OBS_MODE'] = changeto hdu.writeto(outfile_evtool,overwrite=True) if(vign==True): withvignetting='yes' outfile_expmap="{}_ExposureMap_en{}.vign{}".format(os.path.join(outfile_dir,key), eband_index, outfile_post) else: withvignetting='no' outfile_expmap="{}_ExposureMap_en{}.novign{}".format(os.path.join(outfile_dir,key), eband_index, outfile_post) if(do_expmap==True and do_init==True): cmd=["expmap", "inputdatasets=%s" %(outfile_evtool), "emin=%s" %(emin_kev), "emax=%s" %(emax_kev), "templateimage=%s" %(outfile_evtool), "singlemaps=%s" %(outfile_expmap), "withvignetting={}".format(withvignetting), "withsinglemaps=yes","withfilebadpix=yes", "withmergedmaps=no", ] print((" ").join(cmd)) test_exe('expmap') os.system((" ").join(cmd)) return outfile_evtool,outfile_expmap def do_expmap_esass(expmaps=None,outfile=None,emin_kev=None,emax_kev=None,refimage=None): cmd=["expmap", "inputdatasets=\'{}\'".format((" ").join(expmaps)), "emin=%s" %(emin_kev), "emax=%s" %(emax_kev), "templateimage=%s" %(refimage), "singlemaps=%s" %(outfile), "withvignetting=yes", "withsinglemaps=no","withfilebadpix=yes", "withmergedmaps=yes", ] print((" ").join(cmd)) test_exe('expmap') os.system((" ").join(cmd)) def do_fimgmerge_ftools(maps=None,outfile=None): """ takes first map as reference and merges the rest """ tmp="maps.txt" f = open(tmp, "w") for jj in range(1,len(maps)): f.write("{},0,0\n".format(maps[jj])) f.close() cmd=["fimgmerge", "infile={}".format(maps[0]), "list=@{}".format(tmp), "outfile={}".format(outfile), "clobber=yes", ] print((" ").join(cmd)) test_exe('fimgmerge') os.system((" ").join(cmd)) def test_exe(program): """ Tests if executable exists in PATH """ import os def is_exe(fpath): return os.path.isfile(fpath) and os.access(fpath, os.X_OK) fpath, fname = os.path.split(program) if fpath: if is_exe(program): return program else: for path in os.environ["PATH"].split(os.pathsep): exe_file = os.path.join(path, program) if is_exe(exe_file): return exe_file 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,refimage=None,scale=1200,crval=None,devmax=30,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 crval: print("ERROR: Please provide center coordinates array: crval=") print("NOTES: See uds.config.wcslist") return if not refimage: print("ERROR: Please provide reference image: refimage=") return hdulist = fits.open(refimage) w = WCS(hdulist[0].header) w.wcs.crval=crval naxis1=hdulist[0].header['NAXIS1'] naxis2=hdulist[0].header['NAXIS2'] if not catalog: print("ERROR: Please provide reference catalog: 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.fits".format(outkey)) print("Output file: {}".format(fout)) 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.fits", ".src.fits") ref_fout=fout.replace(".cross.fits", ".ref.fits") 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 = [] delta_sep = [] for rec in tbdata: if not (rec['det_like'] > dlmin and rec['det_like'] < dlmax): continue if (rec['id'] in skip.keys()): print("Skip ID={} DET_LIKE={:.2f}".format(rec['id'],rec['det_like'])) else: print("Take ID={} DET_LIKE={:.2f}".format(rec['id'],rec['det_like'])) 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']) sep=src_crd.separation(ref_crd).arcsec delta_sep.append(sep) print("delta [arcsec] RA, Dec, sep {:.4f} {:.4f}".format((rec['src_ra']-rec['ref_ra'])*3600, (rec['src_dec']-rec['ref_dec'])*3600,sep)) tstart, tstop = read_tstart_tstop(infile=refimage) date_obs = mission2date_utc(tstart) date_end = mission2date_utc(tstop) stat_fout=fout.replace(".cross.fits", ".dat") with open(stat_fout, 'w') as writer: writer.write("[DATE-OBS {}, DATE-END {}] nref, mean [arcsec] dRA, dDec, mean/median sep (mean/median): {} & {:.4f} & {:.4f} & {:.4f}/{:.4f}\n".format( date_obs,date_end,len(delta_ra), statistics.mean(delta_ra)*3600, statistics.mean(delta_dec)*3600, statistics.mean(delta_sep), statistics.median(delta_sep))) 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() def do_adapt_ciao(infile=None,outfile=None,expmap=None,function='tophat',expcut=None): if not (infile and expmap and outfile): print("ERROR: Please provide input and output files") return radfile=infile.replace(".fits", ".adapt.scale.fits") sumfile=infile.replace(".fits", ".sum.scale.fits") test_exe('dmimgadapt') cmd=["dmimgadapt", "infile={}".format(infile), "outfile={}".format(outfile), "function={}".format(function), "minrad=1", "maxrad=30", "numrad=30", "radscale=log", "counts=30", "radfile={}".format(radfile), "sumfile={}".format(sumfile), "clobber=yes",] remove_file(outfile) print((" ").join(cmd)) os.system((" ").join(cmd)) cmd=["dmimgadapt", "infile={}".format(infile), "innormfile={}".format(expmap), "outfile={}".format(outfile), "function={}".format(function), "minrad=1", "maxrad=30", "numrad=30", "radscale=log", "counts=30", "inradfile={}".format(radfile), "sumfile={}".format(sumfile), "clobber=yes",] remove_file(outfile) print((" ").join(cmd)) os.system((" ").join(cmd)) if (expcut): cntmap_data, cntmap_hdr = fits.getdata(outfile, ext=0, header=True) expmap_data, expmap_hdr = fits.getdata(expmap, ext=0, header=True) index = np.where(expmap_data < expcut) cntmap_data[index]=0.0 fits.writeto(outfile, cntmap_data, cntmap_hdr, overwrite=True) """ dmimgadapt artmap_ait.fits adapt.fits tophat 1 30 30 log 30 radfile=map.scale.fits sumfile=sum.fits clobber=yes #dmimgadapt artmap_ait.fits adapt.fits func=gaussian min=1.0 max=30 num=30 rad=linear counts=30 sumfile=sum.fits normfile=norm.fits radfile=scale.fits """ def read_tstart_tstop(infile=None): """ looks like evtool does not set START-TSTOP values correctly. Take these times from event list. """ if not infile: print("ERROR: Please provide input file: infile=") return #print("Reading {}".format(infile)) hdul = fits.open(infile) tbdata_ref = hdul[1].data tstart=tbdata_ref['TIME'][0] tstop=tbdata_ref['TIME'][-1] #print("*** TSTART diff {:.1f} ks".format((tstart-hdul[1].header['TSTART'])/1000)) #print("*** TSTOP diff {:.1f} ks".format((hdul[1].header['TSTOP']-tstop)/1000)) return tstart, tstop def make_rate_map(cntmap=None, expmap=None, outfile=None): if not (cntmap and expmap and outfile): print("ERROR: Please provide counts, exposure ane output file") return cntmap_data, cntmap_hdr = fits.getdata(cntmap, ext=0, header=True) expmap_data, expmap_hdr = fits.getdata(expmap, ext=0, header=True) rate = np.zeros(expmap_data.shape) index = np.where(expmap_data > 0.0) rate[index]=cntmap_data[index]/expmap_data[index] fits.writeto(outfile, rate, overwrite=True) def wcs_match_ciao(catalog=None, method='trans',radius=15,residlim=7): xfm=catalog.replace(".fits", ".xfm") src_list=catalog.replace(".fits", ".shu2019.src.fits") if not (os.path.isfile(src_list)==True): print("Not found {}".format(src_list)) print("Did you run cross-match?") sys.exit() ref_list=catalog.replace(".fits", ".shu2019.ref.fits") if not (os.path.isfile(ref_list)==True): print("Not found {}".format(ref_list)) print("Did you run cross-match?") sys.exit() log=catalog.replace(".fits", ".xfm.log") cmd=["wcs_match", "infile={}".format(src_list), "refsrcfile={}".format(ref_list), "outfile={}".format(xfm), "wcsfile={}".format(src_list), "logfile={}".format(log), "radius={}".format(radius), "residlim={}".format(residlim), "verbose=1", "method={}".format(method), "clobber=yes", ] test_exe('wcs_match') os.system((" ").join(cmd)) print((" ").join(cmd)) def wcs_update_ciao(events,ext=2,crval=None,transformfile=None,clean=True): """ Prepares eRosita attitude file for Chandra's wcs_update command """ if not (transformfile): print("ERROR: Please provide transformfile") return if not (os.path.isfile(transformfile)==True): print("Not found {}".format(transformfile)) print("Did you run wcs_update_ciao?") sys.exit() fnaspsol=events.replace(".fits", ".aspsol") fntmp=events.replace(".fits", ".tmp") fnattcor=events.replace(".fits", ".attcorr.fits") hdul = fits.open(events,mode='readonly') corratt=hdul[ext] naxis1=corratt.header['NAXIS1'] naxis2=corratt.header['NAXIS2'] corratt_extname=corratt.header['EXTNAME'] """ create new empty Q_ATT column """ q_att=np.zeros((naxis2,4)) cols = [] cols.append(fits.Column(name='q_att', format='4D', array=q_att)) orig_cols = corratt.columns new_cols = fits.ColDefs(cols) hdu0 = fits.PrimaryHDU() hdu0.header['EXTNAME'] = ('ASPSOL', 'Name of this binary table extension') hdu0.header['HDUNAME'] = ('ASPSOL', 'ASCDM block name') if(crval): hdu0.header['RA_NOM']=(crval[0], '[deg] Pointing RA') hdu0.header['DEC_NOM']=(crval[1], '[deg] Pointing Dec') hdu0.header['ROLL_NOM']=(0.0 , '[deg] Pointing Roll') hdu0.header['RA_PNT']=(crval[0], '[deg] Nominal RA') hdu0.header['DEC_PNT']=(crval[1], '[deg] Nominal Dec') hdu0.header['ROLL_PNT']=(0.0 , '[deg] Nominal Roll') hdu = fits.BinTableHDU.from_columns(orig_cols + new_cols, header=hdu0.header) hdu.writeto(fnaspsol, overwrite=True) """ First, update attitude information in ASPSOL (Chandra notation), which corresponds to CORRATT (eRosita notation) """ cmd=["wcs_update", "infile={}".format(fnaspsol), "outfile={}".format(fntmp), "transformfile={}".format(transformfile), "clobber=yes", ] os.system((" ").join(cmd)) print((" ").join(cmd)) """ Replace attitude extension """ shutil.copyfile(events, fnattcor) data, hdr = getdata(fntmp, 1, header=True) hdr['EXTNAME']=corratt_extname update(fnattcor, data, ext, header=hdr) """ Second, calculate new RA/Dec for events using updated attitude, using evatt from eSASS """ cmd=["evatt", "EVENTFILE={}".format(fnattcor), "ATTFILE={}".format(fnattcor), "GTIFILE={}".format(fnattcor), ] os.system((" ").join(cmd)) print((" ").join(cmd)) if(clean==True): remove_file(fnaspsol) remove_file(fntmp) return fnattcor def create_detmask_merged(expmaps,outfile,minval=None): """ creates mask from each exposure and then stacks each masks. Header is taken from first exposure file. """ tmlist={} for expmap in expmaps: hdul = fits.open(expmap) emap = hdul[0].data ehdr = hdul[0].header instrume = ehdr['INSTRUME'] threshold = minval if(minval) else 0 mask = np.where(emap > threshold, 1, 0) print("--- detmask {} -- {}".format(expmap,instrume)) if not (instrume in tmlist.keys()): tmlist[instrume]=1 else: pass if 'merged_mask' in locals(): merged_mask = np.add(merged_mask, mask) else: merged_mask = mask merged_hdr = ehdr if 'merged_hdr' in locals(): index=1 for tm in tmlist.keys(): merged_hdr["INSTRUM{}".format(index)]=tm index=index+1 merged_hdr['INSTRUME']='merged' merged_hdr['NINST']=len(tmlist) merged_hdr['OBS_MODE']=' ' if 'merged_mask' in locals() and 'merged_hdr' in locals(): fits.writeto(outfile, np.where(merged_mask > 0, 1, 0), header=merged_hdr, overwrite=True) def create_expmap_merged(expmaps,outfile,scale=7.0): """ Adds exposure from TMs. Header is taken from first exposure file. """ tmlist={} for expmap in expmaps: hdul = fits.open(expmap) emap = np.divide(hdul[0].data, scale) ehdr = hdul[0].header instrume = ehdr['INSTRUME'] print("--- expmap {} -- {}".format(expmap,instrume)) if not (instrume in tmlist.keys()): tmlist[instrume]=1 else: pass if 'merged_map' in locals(): merged_map = np.add(merged_map, emap) else: merged_map = emap merged_hdr = ehdr if 'merged_hdr' in locals(): index=1 for tm in tmlist.keys(): merged_hdr["INSTRUM{}".format(index)]=tm index=index+1 merged_hdr['INSTRUME']='merged' merged_hdr['NINST']=len(tmlist) merged_hdr['OBS_MODE']=' ' if 'merged_map' in locals() and 'merged_hdr' in locals(): fits.writeto(outfile, merged_map, header=merged_hdr, overwrite=True) def do_erbackmap_esass(image,expimage,boxlist,detmask,emin,emax,outfile_backmap,cheese_mask): test_exe('erbackmap') cmd=["erbackmap", "image=%s" %(image), "expimage=%s" %(expimage), "boxlist={}".format(boxlist), "detmask=%s" %(detmask), "emin=%s" %(emin), "emax=%s" %(emax), "bkgimage=%s" %(outfile_backmap), "cheesemask=%s" %(cheese_mask), "idband=1", "scut=0.001", "mlmin=6", "maxcut=0.5", "fitmethod=smooth smoothval=15", "snr=40.", ] remove_file(cheese_mask) remove_file(outfile_backmap) os.system((" ").join(cmd)) print((" ").join(cmd)) def filter_catprep(filename, expcut=100,dlmin=6.0,dlmax=10,scale=60*60,ext_like=0.0,outkey='main'): if(os.path.isfile(filename)==False): print("File not found {}".format(filename)) print("Filter catprep {}".format(filename)) fout_selected=filename.replace(".fits", ".{}.selected.reg".format(outkey)) fout_skip=filename.replace(".fits", ".{}.skip.reg".format(outkey)) fout_extended=filename.replace(".fits", ".extended.reg") hdul = fits.open(filename) tbdata = hdul[1].data catsel=[] catskip=[] catext=[] skip_count=0 selected_count=0 keepalive_count=0 for rec in tbdata: if (rec['ext_like'] > ext_like): catext.append({'ra':rec['ra'],'dec':rec['dec'],'radec_err':rec['radec_err'], 'det_like':rec['det_like_0'],'ext_like':rec['ext_like'], 'ext':rec['ext'],'ext_err':rec['ext_err'], 'src_id':rec['id_src']}) if ((rec['det_like_0'] > dlmin and rec['det_like_0'] < dlmax)): catsel.append({'ra':rec['ra'],'dec':rec['dec'],'radec_err':rec['radec_err'],'det_like':rec['det_like_0'], 'src_id':rec['id_src']}) selected_count=selected_count + 1 else: catskip.append({'ra':rec['ra'],'dec':rec['dec'],'radec_err':rec['radec_err'],'det_like':rec['det_like_0'], 'src_id':rec['id_src']}) skip_count = skip_count+1 print("total={} skip_count={} selected={}".format(len(tbdata),skip_count,selected_count)) with open(fout_selected, 'w') as writer: for rec in catsel: writer.write("fk5;circle({}, {}, {}) # color=white text={{{} {:.2f}}}\n".format(rec['ra'],rec['dec'],rec['radec_err']/scale, rec['src_id'],rec['det_like'])) with open(fout_skip, 'w') as writer: for rec in catskip: writer.write("fk5;circle({}, {}, {}) # color=red text={{{} {:.2f}}}\n".format(rec['ra'],rec['dec'],rec['radec_err']/scale, rec['src_id'],rec['det_like'])) with open(fout_extended, 'w') as writer: for rec in catext: writer.write("fk5;circle({}, {}, {}) # color=magenta text={{{} {:.2f} {:.2f}}}\n".format(rec['ra'], rec['dec'], rec['ext']/scale, rec['src_id'], rec['det_like'], rec['ext_like'])) def filter_mllist(filename, expcut=100,dlcut=6.0,dlmin=6.0,dlmax=10,scale=60*60,ext_like=0.0): if(os.path.isfile(filename)==False): print("File not found {}".format(filename)) print("Filter mllist {}".format(filename)) fout_selected=filename.replace(".fits", ".selected.reg") fout_skip=filename.replace(".fits", ".skip.reg") fout_extended=filename.replace(".fits", ".extended.reg") hdul = fits.open(filename) tbdata = hdul[1].data """ create list of unique source IDs """ id_src_set = set(tbdata['id_src']) src_id_list = list(id_src_set) catsel=[] catskip=[] catext=[] skip_count=0 selected_count=0 keepalive_count=0 for src_id in src_id_list: tab_tm0=tbdata[(tbdata['id_src'] == src_id) & (tbdata['id_inst'] == 0) & (tbdata['id_band'] == 0)] """ averaged over all instruments and bands """ ext_tm0=tbdata[(tbdata['id_src'] == src_id) & (tbdata['id_inst'] == 0) & (tbdata['id_band'] == 0) & (tbdata['ext_like'] > ext_like)] """ averaged over all instruments and bands """ if(len(ext_tm0)>0): catext.append({'ra':ext_tm0['ra'][0],'dec':ext_tm0['dec'][0],'radec_err':ext_tm0['radec_err'][0], 'det_like':ext_tm0['det_like'][0],'ext_like':ext_tm0['ext_like'][0], 'ext':ext_tm0['ext'][0],'ext_err':ext_tm0['ext_err'][0], 'src_id':ext_tm0['id_src'][0]}) mask_src_id = (tbdata['id_src'] == src_id) & (tbdata['id_inst'] == 0) """ get all source records except merged one """ tab_src=tbdata[(tbdata['id_src'] == src_id) & (tbdata['id_inst'] != 0) & (tbdata['id_band'] == 1) & (tbdata['ml_exp'] > expcut)] """ only individual instruments, and only first energy band is selected here """ keepalive=False mask = tab_src['det_like'] > dlcut if (any(mask) and tab_tm0['det_like'] < dlmin): keepalive=True keepalive_count=keepalive_count+1 if(keepalive): tab_src=tbdata[(tbdata['id_src'] == src_id) & (tbdata['id_inst'] != 0) & (tbdata['id_band'] == 1)] print("KEEP ALIVE ID={} DL={:.2f} | radec: {:.4f} {:.4f} DL0={:.2f} TEXP={:.1f}".format(src_id, tab_tm0['det_like'][0], tab_tm0['ra'][0], tab_tm0['dec'][0], tab_src['det_like'][0],tab_src['ml_exp'][0])) if ((tab_tm0['det_like'] > dlmin and tab_tm0['det_like'] < dlmax) or keepalive): catsel.append({'ra':tab_tm0['ra'][0],'dec':tab_tm0['dec'][0],'radec_err':tab_tm0['radec_err'][0],'det_like':tab_tm0['det_like'][0], 'src_id':tab_tm0['id_src'][0]}) selected_count=selected_count + 1 else: catskip.append({'ra':tab_tm0['ra'][0],'dec':tab_tm0['dec'][0],'radec_err':tab_tm0['radec_err'][0],'det_like':tab_tm0['det_like'][0], 'src_id':tab_tm0['id_src'][0]}) skip_count = skip_count+1 print("total={} skip_count={} keepalive_count={} selected={}".format(len(src_id_list),skip_count,keepalive_count,selected_count)) with open(fout_selected, 'w') as writer: for rec in catsel: writer.write("fk5;circle({}, {}, {}) # color=white text={{{} {:.2f}}}\n".format(rec['ra'],rec['dec'],rec['radec_err']/scale, rec['src_id'],rec['det_like'])) with open(fout_skip, 'w') as writer: for rec in catskip: writer.write("fk5;circle({}, {}, {}) # color=red text={{{} {:.2f}}}\n".format(rec['ra'],rec['dec'],rec['radec_err']/scale, rec['src_id'],rec['det_like'])) with open(fout_extended, 'w') as writer: for rec in catext: writer.write("fk5;circle({}, {}, {}) # color=magenta text={{{} {:.2f} {:.2f}}}\n".format(rec['ra'],rec['dec'],rec['ext']/scale, rec['src_id'],rec['det_like'],rec['ext_like'])) def create_sensmap(sensmap=None,expmap=None,backmap=None,detmask=None,areatab=None, emin=None,emax=None,ecf=None, detlike=10): test_exe('ersensmap') cmd=['ersensmap', "expimages=\'{}\'".format(expmap), "bkgimages=\'{}\'".format(backmap), "detmasks=\'{}\'".format(detmask), "sensimage={}".format(sensmap), "emin=\'{}\'".format(emin), "emax=\'{}\'".format(emax), "ecf=\'{}\'".format(ecf), "area_table={}".format(areatab), "method=APER", "aper_type=BOX", "aper_size=4.5", "likemin={}".format(detlike), "extlikemin=6.0", "ext=6.0", "extentmodel=beta", "detmask_flag=yes", "shapelet_flag=no", "photon_flag=no", "ext_flag=no", "extlike_flag=no", "compress_flag=no", "area_flag=no",] remove_file(sensmap) remove_file(areatab) os.system((" ").join(cmd)) print((" ").join(cmd)) def group_spectra(path,group_min=25): test_exe('grppha') flist=glob.glob(path) for f in flist: fout=f.replace(".fits", ".fits.grp") print(f) cmd=["grppha", "infile={}".format(f), "outfile={}".format(fout), "comm=\'GROUP MIN {}\'".format(group_min), "tempc=EXIT", "clobber=yes", "chatter=0", # with chatter <= 5 being very quite and chatter >= 20 very verbose. ] print((" ").join(cmd)) os.system((" ").join(cmd))