#!/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 = 'freeFWCandOOTimage'
VERSION = '1.7';
author = 'Jose Vicente Perea';
date = '2021-02-17';
ChangeLog
=========
Version 1.7 - 2021-02-17 (JVP)
------------
+ For last versions of matplotlib imshow() the default interpolation = 'antialiased'
and gives unwanted 3Colour images. Back to interpolation = 'none'
Version 1.6 - 2019-03-25 (JVP)
------------
+ Skip free FWC/OoT 3Col final image if input Raw 3Col image is all Null values
Show OoT image in the diagnostic plot if OoT 3Col image is != all Null values
Version 1.5 - 2018-10-08 (JVP)
------------
+ A scale factor applied to the OoT image subtraction. Weigthed by exposure times:
OoTimage / { 1 + [(Sum(emos1) + Sum(emos2)) / Sum(epn)] }
Where Sum(instr) == Sum of exposure times per instrument
+ PNG image size and aspect ratio modified
Version 1.4 - 2018-10-02 (JVP)
------------
+ Scale factors for FWC and ToO image subtraction included in the Primary clean image
Keywords: FWCSCAL and OOTSCAL
Version 1.3 - 2018-08-28 (JVP)
------------
+ Bug fix: Original Prime FITS extension copied to the Prime extension of the final FITS file
+ PNG image width modified
+ Instruments, Filters and Exposure time labels relocated
Version 1.2 - 2018-06-15 (JVP)
------------
+ Filename label in png plot corrected
Extension names modified
Version 1.1 - 2018-06-12 (JVP)
------------
+ FITS file extension numbers re-coded
Version 1.0 - 2018-05-25 (JVP)
------------
+ Initial version
Task for drawing free FWC and OOT EPIC images
"""
__version__ = "1.7"
__author__ = "Jose V. Perea"
__copyright__ = "ESA (C) 2000-2018"
__license__ = "GPL"
#------------------------------------------------------------------------------------------------
def usage():
#------------------------------------------------------------------------------------------------
'''
Help.
'''
print("""
usage: freeFWCandOOTimage.py --input3COLimage=<...> --FWC3COLimage=<...> --OOT3COLimage=<...> [options]
INTRODUCTION
The freeFWCandOOTimage task .....
OPTIONS:
Input:
-i, --input3COLimage= Input EPIC 3Color image
-b, --FWC3COLimage= Filter Wheel Closed 3Color image
--FWCbscale= Scale factor for FWC image
--OOT3COLimage= Out-Of-Time 3Color image
--OOTbscale= Scale factor for OOT image
Output:
-o, --plotfile= Output PNG/PDF image file.
--outformat={png,pdf} Output image format. Default: png
--interactplot Output plot in an interactive graphics window (Y) or in a PNG file (N)
(not implemented yet)
--verbose= Increase verbosity (Yes/No, Y/N)
-v, --version Show version message
Help:
-h, --help Print this help.
DETAILS
...
LICENCE
...
DOCUMENTATION
...
...
===================================================================================================
""")
#------------------------------------------------------------------------------------------------
def error(string):
'''
prints the string in red.
'''
#import sys
#sys.stdout.write("\033[39;31m" + string + "\033[39;29m") # Red color string
sys.stdout.write( string ) # No color string
def version(version_num):
'''
prints version number.
'''
#import sys
#versionString = sys.argv[0].split('/')[-1] + ' ' + version_num
versionString = '{0} ({0}-{1}) [-]'.format(sys.argv[0].split('/')[-1], version_num)
print(versionString)
# ----------------------------------------------------------------------------------------
def convertImage(img):
import numpy as np
import warnings
convertImageOut = {}
ima2 = np.transpose(img, [1,2,0]) # Transpose the 3-D image ( RGB 3Colour image )
#print('Orig image shape : ', img.shape)
#print('Transpose image shape : ', ima2.shape)
with warnings.catch_warnings(): # Logarithmic scape
warnings.simplefilter("ignore", category=RuntimeWarning)
ima3 = (np.log10(ima2)+5.)/(5.-2.52287875)
# Image size (pixels)
ii = np.nansum(ima2, axis=2).nonzero() # Pixels index. X, Y min and max
xmin = ii[0].min()
xmax = ii[0].max()
ymin = ii[1].min()
ymax = ii[1].max()
# jj = list(map(tuple,np.where(np.isnan( ima3 ))))
# if len(jj[0]) > 0:
# ima3[jj]=0.0 # ; 1.0 for white 0 for black
ima3 = np.nan_to_num(ima3)
ima3[ ima3 > 1. ] = 1. # All values > 1 are set to 1 ( ~ normalization )
ima3[ ima3 < 0. ] = 0
#ima4 = ima3
convertImageOut['image'] = ima3
convertImageOut['limits'] = { 'xmin': xmin, 'xmax': xmax, 'ymin': ymin, 'ymax': ymax }
return convertImageOut
# ----------------------------------------------------------------------------------------
def main(arguments):
'''
This is the main function doing the full loop.
'''
import sys, os, getopt
import ModuleUtil
version_num = __version__ # Version number to be updated
verbose = False
diagplot = False
interactplot = False
doallimagesFITS= False
input3COLimage = ''
FWC3COLimage = ''
OOT3COLimage = ''
FWCbscale = 1.0 # Default scale factor for subtraction
OOTbscale = 1.0 # Default scale factor for subtraction
infileCopy = ''
diagPlotFile = 'diagnosticPlot.png'
plotfile = 'outfile.png'
plotformat = 'png'
allimagesFITSfile = 'allimagesFITSfile.fits'
arguments = arguments.split()
progName = arguments[0]
argsList = arguments[1:]
try:
opts, args = getopt.getopt(argsList, "i:b:o:vh", ["input3COLimage=", "FWC3COLimage=", "FWCbscale=", "OOT3COLimage=", "OOTbscale=", \
"plotfile=", "plotformat=", "interactplot=", \
"diagplot=", "diagPlotFile=", \
"allimagesFITSfile=", "doallimagesFITS=", \
"version" , "verbose=", "help"])
except getopt.GetoptError as err:
usage()
error("ERROR: Could not understand the options!\n\n")
print(err) # will print something like "option -a not recognized"
print("================================================================================================")
return
for o, a in opts:
if o in ("-v", "--version"):
version(version_num)
return 0
if o in ("--verbose"):
verbose = ModuleUtil.varCheck(client = progName, inVar = a).checkBool()
if o in ("--interactplot"):
interactplot = ModuleUtil.varCheck(client = progName, inVar = a).checkBool()
if o in ("--diagplot"):
diagplot = ModuleUtil.varCheck(client = progName, inVar = a).checkBool()
if o in ("--diagPlotFile"):
diagPlotFile = a
if o in ("--doallimagesFITS"):
doallimagesFITS = ModuleUtil.varCheck(client = progName, inVar = a).checkBool()
if o in ("--allimagesFITSfile"):
allimagesFITSfile = a
if o in ("-h", "--help"):
usage()
return 0
if o in ("-i", "--input3COLimage"):
input3COLimage = a
if o in ("-b", "--FWC3COLimage"):
FWC3COLimage = a
if o in ("--FWCbscale"):
FWCbscale = float(a)
if o in ("--OOT3COLimage"):
OOT3COLimage = a
if o in ("--OOTbscale"):
OOTbscale = float(a)
if o in ("-o", "--plotfile"):
plotfile = a
if o in ("--plotformat"):
plotformat = a
# If plotfile name already has a valid extension (png or pdf): keep it.
if plotfile.lower().endswith(('.png', '.pdf')):
plotfile = plotfile
else:
plotfile = plotfile + '.' + plotformat
if interactplot: print("Interactive output ...")
infile = input3COLimage
bkgfile = FWC3COLimage
ootfile = OOT3COLimage
# if outfile == '':
# outfile = infile.replace('.FIT','.PNG').replace('.FTZ','.PNG')
# outfile = os.path.basename( outfile )
if verbose: print( 'doallimagesFITS = ', doallimagesFITS )
if verbose: print('Scale factors for FWC and OOT = ', FWCbscale, OOTbscale )
import numpy as np
import warnings
if not interactplot:
import matplotlib
matplotlib.use('Agg') # To avoid problems with DISPLAY variable
# raise RuntimeError('Invalid DISPLAY variable')
# RuntimeError: Invalid DISPLAY variable
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS
with fits.open(infile) as image:
ima = image[0].data
hdr = image[0].header
### writeFITSExtension( EXT_NUM, 'EXT_NAME', ima )
INSTRUMENTS = hdr['INSTRUME'].split()
EXPOSURES = hdr['EXPOSURE'].split()
if os.path.isfile(ootfile) and os.access(ootfile, os.R_OK): # Get OOT factor if OOT background image exists
emos1_exposures = []
emos2_exposures = []
epn_exposures = []
for idx, instr in enumerate(INSTRUMENTS):
#print('\t{0}\t\t{1}\t{2}'.format(idx, instr, EXPOSURES[idx]))
if instr == 'EMOS1':
emos1_exposures.append( int(EXPOSURES[idx]) )
if instr == 'EMOS2':
emos2_exposures.append( int(EXPOSURES[idx]) )
if instr == 'EPN':
epn_exposures.append( int(EXPOSURES[idx]) )
emos1_total_exposure_time = sum(emos1_exposures)
emos2_total_exposure_time = sum(emos2_exposures)
epn_total_exposure_time = sum(epn_exposures)
OoT_factor = ( 1 + ( (emos1_total_exposure_time + emos2_total_exposure_time)/epn_total_exposure_time ) )
if verbose: print( 'OoT_factor { 1 + [(Sum(emos1) + Sum(emos2)) / Sum(epn)] } \t\t: ', OoT_factor )
if np.nan_to_num( ima ).min() == 0.0 and np.nan_to_num( ima ).max() == 0:
print( 'input3COLimage is all null values. EPIC 3Col image not created' )
return 0
OrigImage = convertImage( ima )
wcs = WCS( hdr, naxis=2 )
#----------------------------------------------
# From PR
#
# w = WCS(naxis=2)
# #if verbose: print ('\n', w , '\n')
# w.wcs.crval = [ hdr['crval1'], hdr['crval2'] ]
# w.wcs.crpix = [ hdr['crpix1'], hdr['crpix2'] ]
# w.wcs.ctype = [ hdr['ctype1'], hdr['ctype2'] ]
# w.wcs.cdelt = [ hdr['cdelt1'], hdr['cdelt2'] ]
#----------------------------------------------
# From PR
#
# if bkgfile != '':
# bkg=fits.open(bkgfile)[0].data
# hbk=fits.open(bkgfile)[0].header
# net=ima
# ib=np.where(np.isnan(bkg) is False)
# net=ima
# net[ib]=ima[ib]-bkg[ib]/4.
# ima=net
#
# if ootfile != '':
# oot=fits.open(ootfile)[0].data
# hoo=fits.open(ootfile)[0].header
# io=np.where(np.isnan(oot) is False)
# net=ima
# print(len(io))
# net[io]=ima[io]-oot[io]*2.06
# ima=net
#----------------------------------------------
if os.path.isfile( bkgfile ) and os.access(bkgfile, os.R_OK): # Correct image if FWC background image exists
with fits.open( bkgfile ) as bkgFWC:
bkghdu = bkgFWC[0]
bkg = bkgFWC[0].data
hbk = bkgFWC[0].header
bkg = np.nan_to_num( bkg ) # NaN values to zero to avoid wrong subtraction
net = ima - ( bkg * FWCbscale )
ima = net
if os.path.isfile(ootfile) and os.access(ootfile, os.R_OK): # Correct image if OOT background image exists
with fits.open( ootfile ) as bkgOOT:
oot = bkgOOT[0].data
hoot = bkgOOT[0].header
oot = np.nan_to_num( oot ) # NaN values to zero to avoid wrong subtraction
net = ima - ( (oot/OoT_factor) * OOTbscale )
ima = net
############################################################
#
# -- Write all images FITS output --
#
if doallimagesFITS:
from shutil import copy
infileCopy = allimagesFITSfile
copy( infile, infileCopy )
rawheader = fits.getheader(infileCopy, 0)
nethdu = fits.PrimaryHDU( ima )
with fits.open( infileCopy, mode='update' ) as hdul:
rawHDU = hdul[0]
extid = 0
hdul[extid] = nethdu # Overwrite the 1st HDU (HDUList[0]) with Net image ( Raw - FWC - OoT)
#hdul[extid].name = 'NET image'
hdul.append( rawHDU ) # Append Raw image
extid = extid + 1
hdul[extid].name = 'RAWIMAGE'
if os.path.isfile( bkgfile ) and os.access(bkgfile, os.R_OK):
with fits.open( bkgfile ) as bkghdul:
fwcHDU = bkghdul[0]
hdul.append( fwcHDU ) # Append FWC image
extid = extid + 1
hdul[extid].name = 'FWCIMAGE'
if os.path.isfile(ootfile) and os.access(ootfile, os.R_OK):
with fits.open( ootfile ) as oothdul:
ootHDU = oothdul[0]
hdul.append( ootHDU ) # Append OoT image
extid = extid + 1
hdul[extid].name = 'OOTIMAGE'
rawheader['FWCSCAL'] = ( FWCbscale, 'Scale factor for FWC image subtraction' )
rawheader['OOTSCAL'] = ( OOTbscale, 'Scale factor for OoT image subtraction' )
hdul[0].header = rawheader # Write the 1st extension 'NET image' with the original file primary extension
hdul.flush()
# Write the 1st extension 'NET image' with the original file primary extension
# with fits.open( infileCopy, mode='update' ) as hdulfinal:
#
# hdulfinal[0].header = rawheader
# #hdulfinal[0].name = 'NET image' # Be aware of this. Rename the extension name ensure that the header is correctly overwritten.
# hdulfinal.flush()
if verbose: print( '\n', fits.info( infileCopy ), '\n' )
# if verbose:
# print( '\n', fits.getheader(infileCopy, 0), '\n' )
############################################################
filteredImage = convertImage( ima )
ima4 = filteredImage['image']
xmin = filteredImage['limits']['xmin']
xmax = filteredImage['limits']['xmax']
ymin = filteredImage['limits']['ymin']
ymax = filteredImage['limits']['ymax']
wcsCrop = wcs[ xmin:xmax, ymin:ymax ] # WCS crop
#ra,dec = w2.all_pix2world( [xmin,xmax], [ymin,ymax], 0 )
#### if verbose: print ('\n', wcsCrop , '\n')
#####################################################################################################
facecolor = '0.1'
if interactplot: # If interactive => diagnostic plot
diagplot = True
if diagplot:
print('Plotting a diagnostic plot ...')
# ------ Image plotting --------------------------------------------
#
plt.style.use('dark_background')
plt.rcParams['patch.linewidth'] = 2
if interactplot:
fig = plt.figure()
else:
fig = plt.figure(figsize = [15, 11])
fig.suptitle('Diagnostic plot')
# Filtered 3Col image
ax = fig.add_axes([0.1, 0.55, 0.4, 0.4], projection = wcsCrop )
ax.imshow( ima4[ xmin:xmax, ymin:ymax, : ], origin='lower', interpolation='none' )
ax.set_title( 'Filtered 3Colour image' )
ax.coords[0].set_major_formatter('hh:mm') # Change RA to 'hh:mm'
ax.set_xlabel('Right Ascension')
ax.set_ylabel('Declination')
# Raw 3Col image
ax = fig.add_axes([0.5, 0.55, 0.4, 0.4], projection = wcsCrop )
ax.imshow( OrigImage['image'][ xmin:xmax, ymin:ymax, : ], origin='lower', interpolation='none' )
ax.set_title( 'Raw 3Colour image' )
ax.coords[0].set_major_formatter('hh:mm') # Change RA to 'hh:mm'
ax.set_xlabel('Right Ascension')
ax.set_ylabel('Declination')
if 'bkg' in locals():
# FWC 3Col image
ax = fig.add_axes([0.1, 0.05, 0.4, 0.4], projection = wcsCrop )
ax.imshow( convertImage( bkg )['image'][ xmin:xmax, ymin:ymax, : ] , origin='lower', interpolation='none' )
ax.set_title( 'FWC 3Colour image' )
ax.coords[0].set_major_formatter('hh:mm') # Change RA to 'hh:mm'
ax.set_xlabel('Right Ascension')
ax.set_ylabel('Declination')
if 'oot' in locals() and ( oot.min() != 0.0 or oot.max() != 0.0 ): # and oot != all Null values
# OOT 3Col image
ax = fig.add_axes([0.5, 0.05, 0.4, 0.4], projection = wcsCrop )
ax.imshow( convertImage( oot )['image'][ xmin:xmax, ymin:ymax, : ] , origin='lower', interpolation='none' )
ax.set_title( 'OoT 3Colour image' )
ax.coords[0].set_major_formatter('hh:mm') # Change RA to 'hh:mm'
ax.set_xlabel('Right Ascension')
ax.set_ylabel('Declination')
if interactplot:
plt.show()
sys.exit()
else:
out = diagPlotFile
print("Saving output diagnostic plot file = ", out)
fig.savefig(out, facecolor = facecolor )
if not interactplot:
# ------ Image plotting --------------------------------------------
#
plt.style.use('dark_background')
plt.rcParams['patch.linewidth'] = 2
dpi = 100
fntcolor = 'white'
#width_image_increase = 1.62
#figsize = [(ima.shape[1] + ima.shape[1] / width_image_increase )/dpi, (ima.shape[2]+30)/dpi]
x_size = ymax - ymin
y_size = xmax - xmin
increase = 300
ratio = x_size / y_size
figsize = [ ( 678*ratio + increase )/dpi , 678/dpi ]
if ratio < 1:
figsize = [ ( 678 + increase )/dpi , 678*(1/ratio)/dpi ]
if verbose: print( 'Image shape\t\t: {0}, Visible data x,y: {1}, {2}'.format( ima.shape, x_size, y_size ) )
if verbose: print( 'Output figure size x,y\t: {0:.2f}, {1:.2f}'.format( figsize[0] * dpi, figsize[1] * dpi) )
if verbose: print( 'DEBUG --- x/y ratio, y increase = ', ratio, increase)
fig = plt.figure( figsize = figsize, dpi=dpi )
ax = fig.add_subplot(1, 1, 1, projection = wcsCrop ) # WCS projection
##plt.subplots_adjust(left=000./dpi/figsize[0], right=1.-100./dpi/figsize[0]) # Shuft image to the left
plt.subplots_adjust( left = -0.05 ) # Shuft image to the left
axra = ax.coords[0]
axdec = ax.coords[1]
axra.set_major_formatter('hh:mm') # Change RA to 'hh:mm'
#axdec.set_major_formatter('dd:mm') # 'dd:mm' = Default
# Plot image:
ax.imshow( ima4[ xmin:xmax, ymin:ymax, : ], origin='lower', interpolation='none' )
ax.set_xlabel('Right Ascension')
ax.set_ylabel('Declination')
fntsize = 12
fntcolor = 'white'
plt.title(os.path.basename(infileCopy), fontsize = fntsize, color = fntcolor)
xrange = ax.get_xlim()
yrange = ax.get_ylim()
fscale = ( yrange[1] - yrange[0] )/fig.dpi*4.
fntsize = 10
INSTRUMENTS = hdr['INSTRUME'].split()
FILTERS = hdr['FILTER'].split()
EXPOSURES = hdr['EXPOSURE'].split()
shift = 0.9
if len(INSTRUMENTS) > 3:
INSTRUMENTS.insert(3, '\n')
FILTERS.insert(3, '\n')
EXPOSURES.insert(3, '\n')
shift = 1.8
if len(INSTRUMENTS) > 6:
INSTRUMENTS.insert(7, '\n')
FILTERS.insert(7, '\n')
EXPOSURES.insert(7, '\n')
shift = 2.5
INSTRUMENTS = ' '.join(INSTRUMENTS)
FILTERS = ' '.join(FILTERS)
EXPOSURES = ' '.join(EXPOSURES)
inst_pos = 1.0
filter_pos = 5.0
exp_pos = 9.0
an = ax.text(xrange[1]+fscale, yrange[1]-fscale*inst_pos, 'Instrument:' ,fontsize=fntsize,color=fntcolor)
an = ax.text(xrange[1]+fscale, yrange[1]-fscale*(inst_pos+shift), INSTRUMENTS ,fontsize=fntsize,color=fntcolor)
an = ax.text(xrange[1]+fscale, yrange[1]-fscale*filter_pos, 'Filter:' ,fontsize=fntsize,color=fntcolor)
an = ax.text(xrange[1]+fscale, yrange[1]-fscale*(filter_pos+shift), FILTERS ,fontsize=fntsize,color=fntcolor)
an = ax.text(xrange[1]+fscale, yrange[1]-fscale*exp_pos, 'Exposure (sec):' ,fontsize=fntsize,color=fntcolor)
an = ax.text(xrange[1]+fscale, yrange[1]-fscale*(exp_pos+shift), EXPOSURES ,fontsize=fntsize,color=fntcolor)
an = ax.text(xrange[1]+fscale, yrange[1]-fscale*14.0, 'Object:' ,fontsize=fntsize,color=fntcolor)
an = ax.text(xrange[1]+fscale, yrange[1]-fscale*14.9, hdr['OBJECT'] ,fontsize=fntsize,color=fntcolor)
an = ax.text(xrange[1]+fscale, yrange[1]-fscale*17.0, 'Observer:' ,fontsize=fntsize,color=fntcolor)
an = ax.text(xrange[1]+fscale, yrange[1]-fscale*17.9, hdr['OBSERVER'] ,fontsize=fntsize,color=fntcolor)
an = ax.text(xrange[1]+fscale, yrange[1]-fscale*20.0, 'DATE-OBS:' ,fontsize=fntsize,color=fntcolor)
an = ax.text(xrange[1]+fscale, yrange[1]-fscale*20.9, hdr['DATE-OBS'] ,fontsize=fntsize,color=fntcolor)
an = ax.text(xrange[1]+fscale, yrange[1]-fscale*23.0, 'Image size:' ,fontsize=fntsize,color=fntcolor)
an = ax.text(xrange[1]+fscale, yrange[1]-fscale*23.9, 'pixels ({0}:{1}, {2}:{3})'.format(xmin, xmax, ymin, ymax ) ,fontsize=fntsize, color=fntcolor)
out = plotfile
print("Saving output plot file = ", out)
fig.savefig(out, facecolor = facecolor )
#####################################################################################################
sys.exit()
# ------------------------------------------------------------------------------------------------------------------------
if __name__=='__main__':
import sys
main(" ".join(sys.argv))
# ------------------------------------------------------------------------------------------------------------------------