#!/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)
########################################################################################