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