#!/usr/bin/env python # """ This file is part of ESA's XMM-Newton Scientific Analysis System (SAS). SAS is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. SAS is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with SAS. If not, see . Name = 'epicSED' VERSION = '1.4'; author = 'Jose Vicente Perea'; date = '2021-03-02'; ChangeLog ========= Version 1.4 - 2021-03-02 (JVP) ------------ + Skip plots if no FLUX columns in OMsrcList Usually when all OM images are WHITE filter Version 1.3 - 2021-02-19 (JVP) ------------ + Importing modules moved after input parameters parser. Version output string fixed Version 1.2 - 2020-11-30 (JVP) ------------ + OM images: Skip OM images plotting if no OM image files or there are not images of the proper filters Version 1.1 - 2020-11-30 (JVP) ------------ + OM images: If more than one image mode per filter discard mode != favorite ( = 'L' ) Version 1.0 - 2020-09-28 (JVP) ------------ + Initial version Task to plot the EPIC and OM combined Spectral Energy Distribution """ __version__ = "1.4" __author__ = "Jose V. Perea" __copyright__ = "ESA (C) 2000-2020" __license__ = "GPL" version_num = __version__ srcPositionsPlot = True verbose = False h_only_one_row = 0.25 s_two_images = 0.10 favoriteOMmode = 'L' ######################################################################################## # def plotOMimages( imageFiles, src_num_coords, srcid_coords, src_num=0, fig='figure', y0_sed=0, font=12, ecol=['k'] ): #print('\tplotOMimages ...') #t0 = time.time() from astropy.visualization import simple_norm ########################### #imageFiles = { 'V': imageFiles['V'] } #imageFiles = { 'V': imageFiles['V'], 'B': imageFiles['B'] } #imageFiles = { 'V': imageFiles['V'], 'B': imageFiles['B'], 'U': imageFiles['U'] } #imageFiles = { 'V': imageFiles['V'], 'B': imageFiles['B'], 'U': imageFiles['U'], 'L': imageFiles['L'] } #imageFiles = { 'V': imageFiles['V'], 'B': imageFiles['B'], 'U': imageFiles['U'], # 'L': imageFiles['L'], 'M': imageFiles['M'] } #imageFiles = { 'V': imageFiles['V'], 'B': imageFiles['B'], 'U': imageFiles['U'], # 'L': imageFiles['L'], 'M': imageFiles['M'], 'S': imageFiles['S'] } # 7 images: #imageFiles = { 'V': imageFiles['V'], 'B': imageFiles['B'], 'U': imageFiles['U'], # 'L': imageFiles['L'], 'M': imageFiles['M'], 'S': imageFiles['S'] } ## print('=====================================') ## for key, val in imageFiles.items(): ## print(f'\t\tkey: {key}\t\tval: {val}') ## print('=====================================') ########################### Nimages = len( [ y for x in imageFiles.values() for y in x ] ) N = Nimages # Return if no OMimages with a standard circle radius if Nimages == 0: return 0, 0.005 * u.degree twoRows=False rows = 1 h = h_only_one_row #y00 = 0.3 y00 = 0.0 if Nimages > 3: # ===> Or > 2 . CHOOSE THE BEST OPTION rows = 2 twoRows=True N = 3 #h = 0.18 h = 0.17 #y00 = y0_sed - 0.23 y00 = 0.175 #print(f'N = {N}') files = imageFiles.values() idx=0 c=0 for filter in filtersSort: try: file = imageFiles[filter] except: continue if len(file) == 0: #print(f'Skipping: {filter} ----------------> CONTINUING') continue elif len(file) == 1: file = file[0] else: print(f'More than 1 image per filter {filter} ==> More tha 1 obs. mode') #------------------------------------------------------------------ idx +=1 with fits.open( file, memmap=False) as OMimage: OMima = OMimage[0].data OMhdr = OMimage[0].header #print(f'OFFSET Filter: {filter} ========= DEOFFSET: ', OMhdr['DEOFFSET'] ) #print(f'OFFSET Filter: {filter} ========= RAOFFSET: ', OMhdr['RAOFFSET'] ) with warnings.catch_warnings(): warnings.simplefilter("ignore", category=Warning) wcsOM = WCS( OMhdr, naxis=2 ) omFilter = OMhdr['FILTER'] # Start ploting #print(f'\tidx: {idx}\t{file}\t{omFilter}') if idx == 4: c=0 center = 0.5 s = 0.03 if Nimages == 2: s = s_two_images w = h x0 = center - ( s* (N-1)/2 + w* (N/2) ) ## y0 = y0_sed - 0.24 ## test y0 = 0.41 y0 = y00 if idx > 3: y0 = y00 - h/2 - s *2.5 if verbose: print(f'\t\trows: {rows}, idx: {idx}, c: {c} ; x0: {x0} y0: {y0}\n') position = [ x0 + w*(c) + s*(c), y0, w, h ] ax_ = fig.add_axes( [ x0 + w*(c) + s*(c), y0, w, h ], projection = wcsOM ) ax_.set_aspect('equal', 'box') if verbose: print(f'\t\tax_ position after set_aspect\t: {ax_.get_position()}') ax_om_pos = ax_.get_position() #ax_.annotate('OM Images:', xy=(0.05, 0.60), xycoords='figure fraction', fontsize=font-2) if idx == 1: ax_.annotate('OM Images:', xy=(0.05, ax_om_pos.y1 + 0.01), xycoords='figure fraction', fontsize=font-2) font_om = font-6 ax_.set_title(omFilter, fontsize=font_om) ax_.coords[0].set_axislabel('RA', fontsize=font_om) ax_.coords[1].set_axislabel('DEC', fontsize=font_om) ax_.tick_params(labelsize=font_om) if c > 0: #ax_.coords[0].set_ticks_visible(False) ax_.coords[1].set_ticklabel_visible(False) ax_.coords[1].set_axislabel('') if twoRows and idx < 4: ax_.coords[0].set_ticklabel_visible(False) #ax_.coords[0].set_axislabel('') c +=1 norm = simple_norm(OMima, stretch='log', percent=99.9) im = ax_.imshow( OMima, norm=norm ) # ...... EPIC Source mark src_num_px, src_num_py = wcsOM.world_to_pixel( src_num_coords['radec_corr'][src_num] ) # position in pixel ####ax_.plot_coord(src_num_coords['radec'][src_num], '+', c='k', ms=80, mew=4) ax_.plot_coord(src_num_coords['radec_corr'][src_num], 'o', c='k', ms=18, mew=1, fillstyle='none') # Green circle ( 20 pix radius ) circ_radius = 20 circ = plt.Circle((src_num_px, src_num_py), radius=circ_radius, ec='green', fill=False, linewidth=1) ax_.add_patch(circ) Aworld= wcsOM.pixel_to_world( src_num_px, src_num_py ) Bworld= wcsOM.pixel_to_world( src_num_px + circ_radius, src_num_py ) circ_radius_world = Aworld.separation(Bworld) # Image zoom around the EPIC source #print(f'src_num_px, src_num_py = {src_num_px, src_num_py}') boxSize = 20 # pixels xlo = src_num_px - boxSize xhi = src_num_px + boxSize ylo = src_num_py - boxSize yhi = src_num_py + boxSize #print(f'--- DEBUG Limits from wcsOM (xlo, xhi), (ylo,yhi) = ({xlo},{xhi}), ({ylo},{yhi})') #print(f'ax_.get_xlim(), ax_.get_ylim() = {ax_.get_xlim()}, {ax_.get_ylim()}') ax_.set_xlim(xlo, xhi) ax_.set_ylim(ylo, yhi) #print('=============== DEBUG ====================') col=0 ## print( f'OM SRCID : {srcid_coords}') for om_id, coords in srcid_coords.items(): ##print(f'DEBUG --- om_id, coords, color = {om_id}, {coords}\t{ecol[col]}') ax_.plot_coord(coords, "o", c=ecol[col], ms=12, mew=1, fillstyle='none') col+=1 if verbose: print() del OMima del OMimage #print('=============== DEBUG ====================') #t1=time.time() #print(f'\tplotOMimages\t: {t1-t0:.4f} seconds') return y0, circ_radius_world ######################################################################################## # def plotEPICimage( epicimage, src_num_coords, srcid_coords, src_num=0, fig='figure', y0_om=0, font=12, ecol=['k'], green_circ_radius=0.1 ): #print(f'\tplotEPICimage ... DEBUG ... font: {font}') #t0=time.time() epicImageFile = epicimage ## EPIC 3COlor image ---------------------------------------------- with fits.open(epicImageFile, memmap=False) as image: ima = image[0].data hdr = image[0].header with warnings.catch_warnings(): warnings.simplefilter("ignore", category=Warning) wcs = WCS( hdr, naxis=2 ) instr = hdr['INSTRUME'] instr = instr.replace(' ',' ').replace(' ','+') #print(f'\t{epicImageFile}\t{instr}') import ModuleUtil convertedImage = ModuleUtil.convertImage( ima ) ima4 = convertedImage['image'] xmin = convertedImage['limits']['xmin'] xmax = convertedImage['limits']['xmax'] ymin = convertedImage['limits']['ymin'] ymax = convertedImage['limits']['ymax'] wcsCrop = wcs[ xmin:xmax, ymin:ymax ] # WCS crop ## EPIC 3COlor image ---------------------------------------------- ###y0 = 0.35 # TEMPORAL VALUE center = 0.5 s = s_two_images #h_epic = h_only_one_row h_epic = 0.28 w_epic = h_epic #x0_epic = center - w_epic/2 x0_epic = center - ( s/2 + w_epic ) #y0_epic = y0_om - 0.3 #y0_epic = 0 y0_epic = 0.335 #ax_epic = fig.add_axes([x0_epic, y0_epic, w_epic, h_epic], projection = wcsCrop) ax_epic_1 = fig.add_axes([x0_epic, y0_epic, w_epic, h_epic], projection = wcsCrop) #print(f'\tax_epic_1 pos : {ax_epic_1.get_position()}') #print(f'\tax_epic_1 pos : {ax_epic_1.get_position().bounds}') # ax_om_pos.y1 + 0.02 #ax_epic_1.annotate(f'EPIC 3Colour combined image ({instr})', xy=(0.05, 0.59), xycoords='figure fraction', fontsize=font-2) ax_epic_1.annotate(f'EPIC 3Colour combined image', xy=(0.05, 0.59), xycoords='figure fraction', fontsize=font-2) #ax_epic_1.annotate(f'EPIC 3Colour image ({instr})', xy=(0.05, ax_epic_1.get_position().y1), xycoords='figure fraction', fontsize=font-2) ax_epic_2 = fig.add_axes([x0_epic + w_epic + s, y0_epic, w_epic, h_epic], projection = wcsCrop) #print(f'\tax_epic_2 pos : {ax_epic_2.get_position()}') #print(f'\tax_epic_2 pos : {ax_epic_2.get_position().bounds}') ax_epic_1.imshow( ima4[ xmin:xmax, ymin:ymax, : ] ) ax_epic_2.imshow( ima4[ xmin:xmax, ymin:ymax, : ] ) font_epic = font-6 ax_epic_1.coords[0].set_axislabel('RA' , fontsize=font_epic) ax_epic_1.coords[1].set_axislabel('DEC', fontsize=font_epic) ax_epic_2.coords[0].set_axislabel('RA' , fontsize=font_epic) ax_epic_2.coords[1].set_axislabel('DEC', fontsize=font_epic) #ax_epic_2.coords[1].set_ticklabel_visible(False) #ax_epic_2.coords[1].set_axislabel('') ax_epic_1.tick_params(labelsize=font_epic) ax_epic_2.tick_params(labelsize=font_epic) # Image zoom around the EPIC source px, py = wcsCrop.world_to_pixel( src_num_coords['radec'][src_num] ) # EPIC source position in pixel boxSize = 30 # pixels xlo = px - boxSize xhi = px + boxSize ylo = py - boxSize yhi = py + boxSize #print(f'--- DEBUG Limits from wcsCrop EPIC (xlo, xhi), (ylo,yhi) = ({xlo},{xhi}), ({ylo},{yhi})') ax_epic_2.set_xlim(xlo, xhi) ax_epic_2.set_ylim(ylo, yhi) # ...... EPIC Source mark #### ax_epic_1.plot_coord(src_num_coords['radec'][src_num], '+', c='k', ms=20) ax_epic_2.plot_coord(src_num_coords['radec'][src_num], '+', c='k', ms=20) # ...... OM Source marks c=0 ##### srcid_coords = src_num_data[src_num]['coords']['om'] #print( f'OM SRCID : {srcid_coords}') for om_id, coords in srcid_coords.items(): ####ax_epic_1.plot_coord(coords, '.', c=ecol[c], ms=10) ax_epic_2.plot_coord(coords, '.', c=ecol[c], ms=6) c+=1 #ax_epic.set_title( f'EPIC 3Colour image\n{instr}', fontsize=font-2 ) # ax_epic_1.set_title( f'EPIC 3Colour image ({instr})', fontsize=font-2 ) from astropy import units as u from astropy.visualization.wcsaxes import SphericalCircle scale = 2 r_1 = SphericalCircle((src_num_coords['radec'][src_num].ra, src_num_coords['radec'][src_num].dec), green_circ_radius * scale, edgecolor='green', facecolor='none', linewidth=1, transform=ax_epic_1.get_transform('fk5')) ax_epic_1.add_patch(r_1) r_2 = SphericalCircle((src_num_coords['radec'][src_num].ra, src_num_coords['radec'][src_num].dec), green_circ_radius, edgecolor='green', facecolor='none', linewidth=1, transform=ax_epic_2.get_transform('fk5')) ax_epic_2.add_patch(r_2) del ima del ima4 del image #t1=time.time() #print(f'\tplotEPICimage\t: {t1-t0:.4f} seconds') return y0_epic ####################################################################################### # def plotSED( epicfile, omfile, epicimage='EPICimage', omimages='omImageFiles', outSrcFiles=None): #t0 = time.time() with warnings.catch_warnings(): warnings.simplefilter("ignore", category=Warning) # To ignore: WARNING: UnitsWarning: 'cts' did not parse as fits unit epicSrcList = Table.read( epicfile, hdu='SRCLIST' ) omSrcList = Table.read( omfile , hdu='SRCLIST' ) omSrcListHeader = fits.getheader(omfile, 'SRCLIST') obsid = fits.getval(omfile, 'OBS_ID', 0) epicSrcList = epicSrcList[ epicSrcList['N_OPT_CO'] > 0 ] epicSrcList = epicSrcList[ 'SRC_NUM', 'RA', 'DEC', 'RA_CORR', 'DEC_CORR', 'EP_1_FLUX','ERR_EP_1_FLUX', 'EP_2_FLUX','ERR_EP_2_FLUX', 'EP_3_FLUX','ERR_EP_3_FLUX', 'EP_4_FLUX','ERR_EP_4_FLUX', 'EP_5_FLUX','ERR_EP_5_FLUX', 'OM_ID1', 'OM_ID1_DIST', 'OM_ID2', 'OM_ID2_DIST', 'OM_ID3', 'OM_ID3_DIST', 'OM_ID4', 'OM_ID4_DIST', 'OM_ID5', 'OM_ID5_DIST', 'N_OPT_CO', 'EP_EXTENT', 'EP_EXT_ERR', 'EP_EXT_ML', 'OMFLAG' ] total_epic_sources = len(epicSrcList) #print(epicSrcList['SRC_NUM','OM_ID1','OM_ID2','OM_ID3','OM_ID4','OM_ID5', 'N_OPT_CO']) #t00 = time.time() #print(f'\tRead SrcLists\t: {t00-t0:.4f} seconds') print( f'Number of EPIC sources with OM counterpart: {total_epic_sources}') elo=[0.2, 0.5, 1.0, 2.0, 4.] # EPIC bands lo ehi=[0.5, 1.0, 2.0, 4.0, 12.] # EPIC bands hi ene = (np.array(ehi)+np.array(elo))/2.0 # Energy position dne = (np.array(ehi)-np.array(elo))/2.0 # Energy "error" = band width #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #print() # ... OM Flux Conversion Factors ...... om_flux_cols = [s for s in omSrcList.colnames if s.endswith('FLUX')] om_filters = [ s.replace('_AB_FLUX','') for s in om_flux_cols ] om_flux_err_cols = [s for s in omSrcList.colnames if s.endswith('FLUX_ERR')] om_all_cols = om_flux_cols + om_flux_err_cols if not len(om_all_cols): print(f'No FLUX columns in omSrcList. Plots not produced.') sys.exit(0) # Conversion factors E2JY = {} for filter in om_filters: # 'E2JY' + filter_name (Ex.: E2JYV value V flux conversion erg/s/cm**2/A) # sqrt( E2JYV * 1E-23 * 3E+18 ) E2JY[filter] = np.sqrt( omSrcListHeader[str('E2JY')+filter] * 1E-23 * 3E+18 ) # ... OM Flux Conversion Factors ...... src_num_data = {} for src_num in epicSrcList['SRC_NUM']: #tloop0 = time.time() #print('------------------------------------------------------------') ###if verbose: print( 'Obsid: {0}. EPIC Source: {1} ...'.format( obsid, src_num ) ) t = epicSrcList[ epicSrcList['SRC_NUM'] == src_num ] # ============ EPIC bands fluxes ============ fep_cols = [s for s in epicSrcList.colnames if (s.endswith('_FLUX') and not s.startswith('ERR_')) ] fep = t[ fep_cols ] efep = t[ [s for s in epicSrcList.colnames if (s.endswith('_FLUX') & s.startswith('ERR_')) ] ] epicT = Table(names = ['fep', 'efep']) for r in fep.colnames: epicT.add_row( [ fep[r], efep['ERR_'+r] ] ) src_num_data[src_num] = {} src_num_data[src_num]['epicData'] = epicT # EPIC SRC_NUM coordinates src_num_coords = {} src_num_coords_corr = {} ecoords = t[ [s for s in epicSrcList.colnames if (s.startswith('RA') or s.startswith('DEC')) ] ] src_num_coords[src_num] = SkyCoord(ecoords['RA'] , ecoords['DEC'] ) src_num_coords_corr[src_num] = SkyCoord(ecoords['RA_CORR'] , ecoords['DEC_CORR'] ) src_num_data[src_num]['coords'] = {} src_num_data[src_num]['coords']['epic'] = {} src_num_data[src_num]['coords']['epic']['radec'] = src_num_coords src_num_data[src_num]['coords']['epic']['radec_corr'] = src_num_coords_corr #print(f'DEBUG------ SRC_NUM: {src_num} Coords: {src_num_coords_corr[src_num]}') #tloop1 = time.time() #print(f'\tCollect EPIC data {src_num}\t\t: {tloop1-tloop0:.4f} seconds') ## DEBUG print(f'\tSRC_NUM : {src_num} ', end='SRCIDs : ') # ============ OM data ============ t_om = t[ [s for s in epicSrcList.colnames if (s.startswith('OM_ID') & s[-1].isdigit()) ] ] srcid_coords = {} srcid_fluxes = {} epic_om_distances = {} for col_id in t_om.colnames: tloop2 = time.time() om_id = t_om[col_id][0] if not om_id: continue om_id_dist = t[ str(col_id)+ '_DIST'][0] epic_om_distances[om_id] = om_id_dist if verbose: print('\t{0}: {1}\tdist to SRC_NUM {3}: {2:.2f} arcsec'.format(col_id, om_id, om_id_dist, src_num)) # OM data: om_data = omSrcList[ omSrcList['SRCID'] == om_id ] if len(om_data) > 1: print('\t* Table size = ', len(om_data)) print('\t* More than 1 OM src with the same SRCID. Skipping OM data') continue om_flux = om_data[ om_all_cols ] omcoords = om_data[ [s for s in om_data.colnames if (s.endswith('RA_CORR') or s.endswith('DEC_CORR'))] ] srcid_coords[om_id] = SkyCoord( omcoords['RA_CORR'], omcoords['DEC_CORR'] ) # - OM wavelength (energy) positions: #om_ene = Column( 12.3984 / np.array( list(E2JY.values()) ) , name='ww') om_ene = Column( np.array( list(E2JY.values()) ) , name='ww') t_om_flux = Table(names = ['om_flux', 'om_flux_err']) for r in om_flux_cols: t_om_flux.add_row( [ om_flux[r], om_flux[r+'_ERR'] ] ) t_om_flux.add_column( om_ene, index=0 ) t_om_flux.add_column( t_om_flux['ww'] * t_om_flux['om_flux'] , name='ff') t_om_flux.add_column( t_om_flux['ww'] * t_om_flux['om_flux_err'] , name='eff') srcid_fluxes[om_id] = t_om_flux #tloop3 = time.time() #print(f'\tCollect OM data, col_id: {col_id}\t: {tloop3-tloop2:.4f} seconds') ## print() ## DEBUG print(f'{om_id} ', end='') #tloop4 = time.time() #print(f'\t{tloop4-tloop0:.4f}') src_num_data[src_num]['omData'] = srcid_fluxes src_num_data[src_num]['distances'] = epic_om_distances src_num_data[src_num]['coords']['om'] = srcid_coords #print('============================================================================') #t1=time.time() #print(f'\tInfo from source lists\t: {t1-t0:.4f} seconds') del epicSrcList del omSrcList #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # # Start PLOT # #............................................. ecol=['blue','red','orange','olive','wheat'] ###font=24 font=14 ms = 6 x0_sed, y0_sed, w_sed, h_sed = [0.10, 0.65, 0.80, 0.30] #............................................. epic_sources = 0 for src_num in src_num_data.keys(): t2 = time.time() epic_sources += 1 # Output plotfile: plotfile = outSrcFiles[str(src_num)] if 'omData' not in src_num_data[src_num].keys(): print('* WARNING: No omData. Skipping plot') continue om_counterparts = list( src_num_data[src_num]['omData'].keys() ) n_om_counterparts = len( om_counterparts ) print(f'* Plot {epic_sources}/{total_epic_sources} - SRC_NUM: {src_num} - Possible OM counterparts {n_om_counterparts}, SRCIDs: {om_counterparts}') fig = plt.figure() #print(f'=====> Orig figure size = {fig.get_size_inches()}') #fig.set_size_inches(15.0,20.0) #fig.set_size_inches(12.0,24.0) #print(f'=====> Orig figure size = {fig.get_size_inches()}') #size_scale = 1.0 # If reduce the size ==> reduce font size and symbols size #fig.set_size_inches( size_scale* 15.0, size_scale* 20.0) fig.set_size_inches(9.0, 12.0) ax = fig.add_axes([x0_sed, y0_sed, w_sed, h_sed]) #print(f'\tax_sed pos : {ax.get_position()}') #print(f'\tax_sed pos : {ax.get_position().bounds}') src_num_fluxes = src_num_data[src_num]['epicData'] ax.errorbar( ene, src_num_fluxes['fep'], xerr = dne, yerr = src_num_fluxes['efep'], fmt='.', c='k', ms=ms, ecolor='k' ) ax.set_xscale('log') ax.set_yscale('log') ax.set_xlabel('Energy (keV)', fontsize=font) ax.set_ylabel('erg/cm$^2$/s', fontsize=font) ax.set_title( f'{obsid} EPIC SRC_NUM# {src_num} ( {str(n_om_counterparts)} possible OM counterparts )\n', fontsize=font ) xlimits = [0.001,20] ax.set_xlim( xlimits ) ax.tick_params(labelsize=font-2) ax.tick_params(direction="in", length=12, which='major') ax.tick_params(direction="in", length=6, which='minor') ax.grid() ax.xaxis.set_major_formatter(FormatStrFormatter("%.3f")) c=0 om_non_NaN = {} srcid_fluxes = src_num_data[src_num]['omData'] epic_om_distances = src_num_data[src_num]['distances'] for om_id, fluxes in srcid_fluxes.items(): ax.errorbar( 12.3984 / fluxes['ww'], fluxes['ff'], yerr = fluxes['eff'], fmt='.', c=ecol[c], ms=ms, ecolor=ecol[c] ) ax.annotate('OM SRCID: {0} (EPIC-OM: {1:.2f} arcsec)'.format(om_id, epic_om_distances[om_id]), xy=(0.25, 0.90 - c*0.05), xycoords='axes fraction', fontsize=font-4, color=ecol[c]) #print('-------------- ', om_id, ' ---------------------------') om_non_NaN[om_id] = np.count_nonzero(~np.isnan( fluxes['ff'] )) c+=1 #t3=time.time() #print(f'\tPlot SED for {src_num}\t: {t3-t2:.4f} seconds') if srcPositionsPlot: pass ##------------------------------------------------------------------------------------------------ ## VISUALIZATION PLOTS ## OM and EPIC 3Colour images ## -------------------------- # ...... OM IMAGE ....................................... # y0_om, om_circ_radius_world = plotOMimages( omimages, src_num_data[src_num]['coords']['epic'], src_num_data[src_num]['coords']['om'], src_num=src_num, fig=fig, y0_sed=y0_sed, font=font, ecol=ecol ) # ...... EPIC 3Color IMAGE ....................................... # y0_epic = plotEPICimage( epicimage, src_num_data[src_num]['coords']['epic'], src_num_data[src_num]['coords']['om'], src_num=src_num, fig=fig, y0_om=y0_om, font=font, ecol=ecol, green_circ_radius=om_circ_radius_world ) # plt.tight_layout() not compatible with Axes with specific positions. It works with ax = plt.subplot(...) print(f'\t{plotfile}') #t4 = time.time() # To ignore: WARNING: invalid value encountered in subtract A_scaled -= a_min with warnings.catch_warnings(): warnings.simplefilter("ignore", category=RuntimeWarning) plt.savefig( plotfile ) plt.close('all') if verbose: print('\tDone') #t5=time.time() #print(f'\tSave png {src_num}\t: {t5-t4:.4f} seconds') #### TO BE COMMENTED <========================== ###if epic_sources == 4: break #if epic_sources == N_sources: break # ######################################################################################## #import sys import argparse #import ModuleUtil # OM filters order ......................... filtersSort = ['V', 'B', 'U', 'L', 'M', 'S'] def input(): # Inut Options parser = argparse.ArgumentParser(description='EPIC OM combined SED', formatter_class=argparse.ArgumentDefaultsHelpFormatter) parser.add_argument('--version', '-v', action='version', version=f'%(prog)s (%(prog)s-{version_num}) [-]') parser.add_argument('--omsrclist' , help='Input EPIC Summary source list file', default='PobsidEPX000OBSMLI0000.FTZ') parser.add_argument('--epicsrclist', help='Input OM Summary source list file' , default='PobsidOMX000OBSMLI0000.FTZ') parser.add_argument('--epic3COLimage', help='Input EPIC 3Colour image FITS file', default='PobsidEPX0003COLIM8000.FTZ') parser.add_argument('--omimages' , help='Input OM image FITS files' , default='PobsidOMX000*SIMAG{filter1}000.FTZ,PobsidOMX000*SIMAG{filter2}000.FTZ') parser.add_argument('--plotformat', help='Plot format (png/pdf)', choices=['png', 'pdf'], default='png') parser.add_argument('--outfiles', help='Output src1:filename1,src2:filename2,...', default='EPIC-OM_SED') args = parser.parse_args() return args def main(args): omSrcList = args.omsrclist epicSrcList = args.epicsrclist epicImage = args.epic3COLimage omImages = args.omimages.split(',') plotformat = args.plotformat outfiles = args.outfiles.split(',') omImageFiles = {} for filter in filtersSort: omImageFiles[filter] = [ x for x in omImages if 'SIMAG'+filter in x ] nOMimages = len( [ y for x in omImageFiles.values() for y in x ] ) # Count all files from all filters discardedOMimages = [] # If the number of OM images > 6 ==> More than 1 mode per filter ==> # Slect most common observing mode = 'L' (Low resolution) = favoriteOMmode # Except if there are more 'H' images (High resolution) # Do not select 'U' because a smaller FoV if nOMimages > 6: OMmodes = [ y[-14] for x in omImageFiles.values() for y in x ] #print(f'OMmodes = {OMmodes}') favorite = favoriteOMmode # 'L' nH = len([ x for x in OMmodes if x == 'H' ]) nL = len([ x for x in OMmodes if x == 'L' ]) if nH > nL: favorite = 'H' for filter,val in omImageFiles.items(): if len(val) > 1: for f in val: OMmode = f[-14] # 'H', 'L', 'U' if OMmode == favorite: omImageFiles[filter] = [f] else: discardedOMimages.append(f) # Look for more than 1 mode per filter # Remove the image != favoriteOMmode = 'L' # for filter, files in omImageFiles.items(): if len(files) > 1: print(f'Number of Modes per filter {filter} > 1 : {files}') for f in files: OMmode = f[-14] # 'H', 'L'(favoriteOMmode), 'U' if OMmode != favoriteOMmode: print(f'Image discarded for Mode {OMmode}: {f}') files.remove(f) omImageFiles[filter] = files discardedOMimages.append(f) # EPIC and OM Source lists: print('\t\tepic\t: {0}'.format(epicSrcList)) print('\t\tom\t: {0}'.format(omSrcList)) # Output plots: src1:filename1, src2:filename2, ... srcFiles = {} for j in outfiles: e = j.split(':') srcFiles[e[0]] = e[1] for src,filename in srcFiles.items(): print(f'\t\t\tsrc_num: {src} file: {filename}') # EPIC and OM images: print(f'\t\t\tepic\t: {epicImage}') if omImageFiles: for filter, files in omImageFiles.items(): print(f'\t\t\tom\t: {filter} : {files}') if 'discardedOMimages' in locals(): print(f'\t\t\tdiscarded : {discardedOMimages}') if not nOMimages: print(f'No OM images in valid filters ({filtersSort})') #t1=time.time() #print(f'\tFiles ...\t: {t1-t0:.4f} seconds') # plotSED and Images......... plotSED( epicSrcList, omSrcList, epicimage=epicImage, omimages=omImageFiles, outSrcFiles=srcFiles ) #################### if __name__=='__main__': argsObj = input() import sys, time time0 = time.time() import warnings import numpy as np from astropy.io import fits from astropy.table import Table, Column from astropy.wcs import WCS from astropy.coordinates import SkyCoord from astropy import units as u #-- -- -- import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt plt.style.use('classic') #%matplotlib inline from matplotlib.ticker import FormatStrFormatter #-- -- -- time1 = time.time() Dt = time1 - time0 ##print(f'\nImport modules time = {Dt:.4f} seconds\n') main(argsObj) sys.exit(0) ########################################################################################