#! /usr/bin/env python # -*- coding: iso-8859-1 -*- # This program finds the footprint of the illuminated of a fits image on the sky. This # footprint is a polygon in pixel (or WCS) coordinates for each chip. # # The input fits images images can # - be composed of multiple chips # - be concave or convex # - have non-straight outlines # - have bad/masked pixels # - have holes # - holes with islands in them # # There is one main free parameter to the code: # - --smoothing specifies the number of skipped contour points in the footprint. # --smoothing = 1 : all points in each polygon are included # --smoothing = 2 : half of the points in each polygon is removed # # Requirements: # - pyfits (http://www.stsci.edu/resources/software_hardware/pyfits) # - pylab (http://matplotlib.sourceforge.net) # - pywcs (documentation: http://stsdas.stsci.edu/astrolib/pywcs/ # download: https://trac.assembla.com/astrolib) # or # - astropy (http://www.astropy.org/) # # Script inspiration by Felix Stoehr, ST-ECF, 2008 # 2017-08-15: Footprints from the original Felix Stoehr version # 2017-09-04: Footprints from the "contour" function capabilities of "matplotlib" # 2017-09-05: Diagnostic plot Image + Footprint to 'matplotlib graphics window' or 'PNG file' # 2017-09-05: No colors for stdout messages # # Summary of the algorithm used: # 1) Find pixels of the image contour in Image or WCS coordinates # 2) Write the footprint polygons in ds9 region format. Image or WCS coordinates # 3) [optional] Show the Image and the Footprint in a Python matplotlib graphics window (interactplot=Y) or PNG file # # ---------------------------------------------------------------------------------------- def main(arguments): # ---------------------------------------------------------------------------------------- ''' This is the main function doing the full loop. Note that pyfits reads y,x which is why all the plotting routines are using inverted coordinates. For the description of directions in the comments, we use this inverted notation (i.e. as plotted on the screen) ''' import getopt, os version_num = "2.2" # Version number to be updated zerovalue = 0 smoothing = 1 extension = None wcsextension = None verbose = False plotting = "no" interactplot = "no" ds9 = False ds9reg = "yes" skycoord = "yes" infile = "exposureMap.ds" outfile = "footprint_ds9reg.reg" diagplot = "diagPlot.png" date = runtimestring() datestring = getdatetime() try: opts, args = getopt.getopt(arguments.split(), "vwr:z:g:e:E:o:f:Vh", ["verbose", "plotting=", "interactplot=", "ds9reg=", \ "zerovalue=", "smoothing=", "extension=", "wcsextension=", \ "infile=", "outfile=", "diagplot=", "skycoord=", "version", "help"]) except getopt.GetoptError: usage() error("ERROR: Could not understand the options!\n\n") print("================================================================================================") return for o, a in opts: if o in ("-V", "--verbose"): verbose = True alert(whocalls(verbose=verbose) + "Verbose ON"+"\n") if o in ("-z", "--zerovalue"): zerovalue = float(a) alert(whocalls(verbose=verbose) + "Zerovalue = %f " % zerovalue+"\n") if o in ("--smoothing"): smoothing = int(a) alert(whocalls(verbose=verbose) + "Smoothing = %i " % smoothing+"\n") if o in ("-e", "--extension"): extension = int(a) alert(whocalls(verbose=verbose) + "FITS file extension = %d " % extension+"\n") if o in ("-E", "--wcsextension"): wcsextension = int(a) alert(whocalls(verbose=verbose) + "FITS file WCS extension = %d " % wcsextension+"\n") if o in ("--ds9reg"): ds9reg = a if o in ("--skycoord"): skycoord = a if o in ("--plotting"): # Plotting Image + Footprint plotting = a if o in ("--interactplot"): # interactplot = Y : Plotting Image + Footprint in a Python graphics window (for debugging) interactplot = a # interactplot = N : Plotting Image + Footprint in a PNG file if o in ("-i", "--infile"): infile = a alert(whocalls(verbose=verbose) + "Input filename = %s " % infile+"\n") if o in ("-o", "--outfile"): outfile = a alert(whocalls(verbose=verbose) + "Output filename = %s " % outfile+"\n") if o in ("--diagplot"): diagplot = a alert(whocalls(verbose=verbose) + "Diagnostic Output file = %s " % diagplot+"\n") if o in ("-v", "--version"): version(version_num) return 0 if o in ("-h", "--help"): usage() return 0 yes = ('yes', 'y') if ds9reg.lower() in yes: ds9 = True alert(whocalls(verbose=verbose) + "ds9 output ON \n") if skycoord.lower() in yes: skycoord = True alert(whocalls(verbose=verbose) + "ds9 output in Sky coord. \n") else: skycoord = False alert(whocalls(verbose=verbose) + "ds9 output in Image coord. \n") if plotting.lower() in yes: plotting = True alert(whocalls(verbose=verbose) + "Plotting ON\n") else: plotting = False if interactplot.lower() in yes: interactplot = True alert(whocalls(verbose=verbose) + "Interactive plot ON. No diagnostic output file\n") else: interactplot = False import matplotlib matplotlib.use('Agg') # To avoid problems with DISPLAY variable # raise RuntimeError('Invalid DISPLAY variable') # RuntimeError: Invalid DISPLAY variable. if not (os.path.isfile(infile) and os.access(infile, os.R_OK)): usage() error("ERROR: Either file " + infile + " is missing or has not read persmissions!\n\n") print("===================================================================================================") return # #print("===================================================================================================") # #alert("\nfootPrintFinder\n") # #alert("Copyright (C) Felix Stoehr, 2008 Space Telescope European Coordinating Facility\n") # alert("\nfootPrintFinder\n") # alert("Jose Vicente Perea, 2017 XMM-Newton SOC. ESAC - ESA\n") # alert("Based on : Copyright (C) Felix Stoehr, 2008 Space Telescope European Coordinating Facility\n") # alert(date+"\n\n") # import pylab # import matplotlib # matplotlib.use('Agg') # To avoid problems with DISPLAY variable # raise RuntimeError('Invalid DISPLAY variable') # RuntimeError: Invalid DISPLAY variable. filename = infile ### borderpixels, dimshape, wcs = directborder(filename, zerovalue, extension, wcsextension, verbose=verbose) segments, dimshape, wcs, image = getbordersegments(filename, zerovalue, extension, wcsextension, skycoord, verbose=verbose) print(whocalls(verbose=verbose)) print(whocalls(verbose=verbose) + "Number of polygons found: %s" % len(segments) ) # print() # print('-------------------------------------------------------') # # print('segments.shape:\t', len(segments), ' polygons') # print() # print('segments[0].shape:\t', segments[0].shape, '\tlen(): ', len(segments[0])) # print(segments[0]) # print() # print('segments[1].shape:\t', segments[1].shape, '\tlen(): ', len(segments[1])) # print(segments[1]) # print() # #print('segments[47].shape:\t', segments[47].shape) # #print(segments[47]) # #print() # print('segments[0][:,0].shape:\t', segments[0][:,0].shape, '\tlen(): ', len(segments[0][:,0])) # print('segments[0][:,0] : ', segments[0][:,0]) # print() # print('segments[0][:,1].shape:\t', segments[0][:,1].shape, '\tlen(): ', len(segments[0][:,1])) # print('segments[0][:,1] : ', segments[0][:,1]) # # print('-------------------------------------------------------') # writeoutput(filename, diagplot, simplefootprintlist, borderpixels, hierarchy, wcs, centroid, arguments, writepoints, verbose) if ds9: writeds9output(filename, outfile, segments, wcs, skycoord, smoothing, verbose) print(whocalls(verbose=verbose) + "Footprint generation finished successfully.") if plotting: plotfootprints(image, segments, dimshape, wcs, skycoord, smoothing, diagplot, interactplot, verbose) alert(whocalls(verbose=verbose) + runtimestring(date)+"\n") #------------------------------------------------------------------------------------------------ def usage(): #------------------------------------------------------------------------------------------------ ''' Help. ''' print(""" usage: footPrintFinder.py --infile= [options] output: footprint_ds9reg.reg INTRODUCTION The footPrintFinder reads a FITS file, identifies the footprint(s) of the illuminated pixels and writes them into an ascii output file. This program finds the footprint of the illuminated of a fits image on the sky. This footprint is a set of polygons in pixel (or WCS) coordinates for each chip. This output file is written to the local directory and contains the definition of the footprint illuminated part of the image in ds9 region format. That definition is written in X, Y polygons in image coordinates or RA, Dec coordinates if the CTYPE1/CTYPE2 are 'RA---TAN/DEC--TAN'. OPTIONS: Input: -i, --infile= Input exposure map image fits FITS file. Default="exposureMap.ds" -e, --extension FITS file extension to use. Default=first extension with data. -E, --wcsextension FITS file extension to use for the WCS. Default=same extension as the data. Processing: -t, --smoothing Specifies the number of skipped contour points in the footprint which smooths the polygons contour. smoothing = 1 (default value) is equivalent to include all points in the footprint. smoothing = 2 is equivalent to include only half of the points in each polygon. Rest of points are ignored. -z, --zerovalue Value of the pixels that should be identified as 'border'. If a pixel value 'nan', it is replaced with the zerovalue, too. Output: -o, --outfile= Output ds9 region file. Default="footprint_ds9reg.reg" --ds9reg= Create output ds9 region file. Default="yes" --skycoord= Output ds9 region file in projection RA--TAN, DEC-TAN, WCS coordinates. Otherwise the ds9 region file will be produced in Image coordinates. Default="yes" ( RA--TAN, DEC-TAN coord. ) -p, --plotting Plots the footprint. Default=off. --interactplot Output plot in an interactive graphics window (Y) or in a PNG file (N) -v, --verbose Increase verbosity. Help: -h, --help Print this help. DETAILS ... LICENCE ... DOCUMENTATION A description of the footPrintFinder is available in the issue #45 of the ST-ECF newsletter http://www.spacetelescope.org/about/further_information/stecfnewsletters/hst_stecf_0045 =================================================================================================== """) #------------------------------------------------------------------------------------------------ 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 alert(string): #------------------------------------------------------------------------------------------------ ''' prints the string in blue. ''' import sys #sys.stdout.write("\033[39;34m" + string + "\033[39;29m") # Blue color string sys.stdout.write( string ) # No color string #------------------------------------------------------------------------------------------------ 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 warning(string): #------------------------------------------------------------------------------------------------ ''' prints the string in orange. ''' import sys sys.stdout.write("\033[39;33m" + string + "\033[39;29m") #------------------------------------------------------------------------------------------------ def runtimestring(date=None): #------------------------------------------------------------------------------------------------ ''' returns a string of the current time and if a date is given, adds the running time from that date. ''' import time, datetime if date==None: return datetime.datetime.now().strftime("%d/%m/%Y %H:%M:%S") else: thisdate = datetime.datetime.now() difference = thisdate-datetime.datetime(*time.strptime(date, "%d/%m/%Y %H:%M:%S")[:6]) minutes, seconds = divmod(difference.seconds, 60) hours, minutes = divmod(minutes, 60) hours = hours + difference.days * 24 return datetime.datetime.strftime(thisdate, "%d/%m/%Y %H:%M:%S") + " (%dh %02dm %02ds)" %(hours, minutes, seconds) #------------------------------------------------------------------------------------------------ def getdatetime(time=None): #------------------------------------------------------------------------------------------------ ''' returns a SYBSE accepted datetime string. If no datetime object is given, it returns the (local) time. ''' import datetime if time == None: time = datetime.datetime.utcnow() return "%04d-%02d-%02d %02d:%02d:%02d" % (time.year, time.month, time.day, time.hour, time.minute, time.second) #------------------------------------------------------------------------------------------------ def whocalls(level=1, verbose=False): #------------------------------------------------------------------------------------------------ ''' This function returns the caller, the function and the line number. ''' import os, sys adjust = 48 x = sys._getframe(level) if verbose: return ('%s|%s:%s:l=%d'%(getdatetime(), os.path.basename(x.f_code.co_filename), x.f_code.co_name, x.f_lineno)).ljust(adjust+20) else: return ('%s:%s:l=%d'%(os.path.basename(x.f_code.co_filename), x.f_code.co_name, x.f_lineno)).ljust(adjust) #------------------------------------------------------------------------------------------------ def showprogress(i, n, lastprogress): #------------------------------------------------------------------------------------------------ ''' This function shows the progress of an iteration. i=current number, n=total number. ''' import sys progress = int(100.0 * i/n) if i == n-1: sys.stdout.write("100%") else: if progress>lastprogress: lastprogress = progress if (progress % 10)==0: sys.stdout.write("%d%%" % progress) sys.stdout.flush() elif (progress % 2)==0: sys.stdout.write(".") sys.stdout.flush() i = i + 1 return i, lastprogress # ---------------------------------------------------------------------------------------- def getbordersegments(filename, zerovalue, extension, wcsextension, skycoord, verbose=False): # ---------------------------------------------------------------------------------------- ''' This function reads the data in slices (for low memory usage), cleans it an then identifies the borderpixels. Pyfits only allows this kind of access for non-scaled data! The crval, crpix amd cdmatrix are read from the file, too. ''' import matplotlib.pyplot as plt #import pylab, numpy import numpy # Importing pyfits or astroy otherwise try: import pyfits except: warning(whocalls(verbose=verbose) +"WARNING: pyfits is not installed!\n") warning(whocalls(verbose=verbose) +"WARNING: trying astropy...\n") try: import astropy.io.fits as pyfits # (re-using the module name : pyfits) print(whocalls(verbose=verbose) + "astropy found! Continuing ...") except: warning(whocalls(verbose=verbose) +"ERROR: Neither pyfits nor astropy are installed!\n") warning(whocalls(verbose=verbose) +"ERROR: Some of them are mandatory requirement. Exiting...\n") import sys sys.exit(1) # Importing pywcs or astroy.wcs otherwise wcs = {'wcs_is_possible': True} try: import pywcs wcs ['wcs_method'] = 'wcs_pix2sky' # WCS.wcs_pix2sky() : method name in pywcs except: warning(whocalls(verbose=verbose) +"WARNING: pywcs is not installed!\n") warning(whocalls(verbose=verbose) +"WARNING: trying astropy.wcs...\n") try: import astropy.wcs as pywcs # (re-using the module name : pywcs) print(whocalls(verbose=verbose) + "astropy.wcs found! Continuing ...") wcs ['wcs_method'] = 'wcs_pix2world' # WCS.wcs_pix2world() : method name in astropy.wcs except: warning(whocalls(verbose=verbose) +"WARNING: Neither pywcs nor astropy.wcs are installed!\n") warning(whocalls(verbose=verbose) +"WARNING: Sky coordinates will thus not be available.\n") warning(whocalls(verbose=verbose) +"WARNING: Please install pywcs from https://trac.assembla.com/astrolib \n") warning(whocalls(verbose=verbose) +"WARNING: or astropy from http://www.astropy.org \n") wcs = {'wcs_is_possible': False} pass degtorad = 0.017453292519943295 print(whocalls(verbose=verbose) + "Reading fits file %s ..." % filename) try: hdulist = pyfits.open(filename) # pyfits could belong to 'pyfits' or 'astropy' except IOError: error("\n"+whocalls(verbose=verbose) + "ERROR: the file could not be opened!\n") print("===================================================================================================") import sys sys.exit(1) hdu = 0 if extension==None: while True: dshape = [] try: try: dshape = hdulist[hdu].shape except: dshape = hdulist[hdu]._dimShape() except IndexError: error("\n"+whocalls(verbose=verbose) + "ERROR: Cannot read file. Is this really a fits file?\n") print("===================================================================================================") import sys sys.exit(1) if list(dshape) == []: hdu += 1 else: break else: hdu = extension try: n = list(hdulist[hdu].shape) except: n = list(hdulist[hdu]._dimShape()) # This will be needed later for extracting the data nn = n print(whocalls(verbose=verbose) + "Found data with %s in hdu %d ..." % (str(n), hdu)) if len(n)!=2: n = n[-2:] print(whocalls(verbose=verbose) + "Assuming that the image dimensions are the last two in that list (%s)" % (str(n))) if verbose: print(whocalls(verbose=verbose) + "Fits header of this extension:") print(hdulist[hdu].header) # For some FITS files, the WCS information is split into different headers and the EQUINOX/RADE(C)SYS # are in the primary header whereas the rest of the WCS is in later headers if 'EQUINOX' in hdulist[0].header: #if hdulist[0].header.has_key('EQUINOX'): wcs['EQUINOX'] = hdulist[0].header['EQUINOX'] if 'RADECSYS' in hdulist[0].header: #if hdulist[0].header.has_key('RADECSYS'): wcs['RADESYS'] = hdulist[0].header['RADECSYS'] if 'RADESYS' in hdulist[0].header: #if hdulist[0].header.has_key('RADESYS'): wcs['RADESYS'] = hdulist[0].header['RADESYS'] if wcs['wcs_is_possible']: wcshdu = hdu if wcsextension!=None: wcshdu = wcsextension try: wcs['wcs'] = pywcs.WCS(hdulist[wcshdu].header, naxis=2, relax=True) # pywcs could belong to 'pywcs' or 'astropy' except Exception as e: warning(whocalls(verbose=verbose) +"WARNING: It was not possible to retrieve WCS information from the FITS file using extension %d \n" %wcshdu) warning(whocalls(verbose=verbose) +"WARNING: pywcs issued the error message " + str(e) + ".\n") warning(whocalls(verbose=verbose) +"WARNING: Sky coordinates will thus not be available.\n") warning(whocalls(verbose=verbose) +"WARNING: Try -E/--wcsextension if the WCS is in a different header than the data.\n") wcs = {'wcs_is_possible': False} pass print(whocalls(verbose=verbose) + "Cleaning and finding borderpixels ...") print(whocalls(verbose=verbose), end='') # Now find the borders segments contours: nchunks = 2 if skycoord and wcs['wcs_is_possible']: ax = plt.subplot(projection = wcs['wcs']) #cs = ax.contour( hdulist[hdu].data, transform=ax.get_transform( wcs['wcs'] ), levels=[zerovalue], nchunk=nchunks) cs = ax.contour( hdulist[hdu].data, transform=ax.get_transform( wcs['wcs'] ), levels=[zerovalue]) # Contours coordinates in WCS: RA, Dec (degrees) else: ax = plt.subplot() #cs = ax.contour( hdulist[hdu].data, levels=[0], nchunk=nchunks) cs = ax.contour( hdulist[hdu].data, levels=[0]) segments = cs.allsegs[0] image = hdulist[hdu].data return segments, n, wcs, image print() # ---------------------------------------------------------------------------------------- def writeds9output(filename, outfile, segments, wcs, skycoord, smoothing, verbose): # ------------------------------------------------------------------------------------------------------------------------ """ writing ds9 region file(s). One file with the image coordinates, and if possible (i.e. of the projection is RA-TAN/DEC-TAN) also a ds9 WCS file. """ thres = 3 # Lower limit of points numbers in a polygon step = smoothing # Step for polygon smoothing tt = [] for s in segments: t='polygon(' if len(s) >= thres: reduced_s = s[::step] if len(reduced_s) == 1: # Polygons of only 1 point are not permitted in ds9 continue for ss in reduced_s: dd = (ss[0],ss[1],0) if skycoord and wcs['wcs_is_possible']: dd = getattr(wcs['wcs'], wcs['wcs_method'])(ss[0], ss[1], 0) t = '{0}{1},{2},'.format(t, dd[0], dd[1]) tt.append(t[:len(t)-1] + ') # color=green width=2\n') # Write output files outputfilename = outfile try: outputfile = open(outputfilename, "w") except IOError: error("\n"+whocalls(verbose=verbose) + "ERROR: the output file could not be opened for writing!\n") print("===================================================================================================") import sys sys.exit(1) print("# Region file format: DS9 version 4.0", file=outputfile) print("# Filename: %s " % filename, file=outputfile) print('# global color=green font="helvetica 12 bold" select=1 highlite=1 edit=1 move=1 delete=1 include=1 fixed=0 source', file=outputfile) ccord_type = 'image' if skycoord and wcs['wcs_is_possible']: ccord_type = 'linear' print(ccord_type, file=outputfile) outputfile.writelines(tt) outputfile.close() print(whocalls(verbose=verbose) + "Written to %s" % outputfilename) # ------------------------------------------------------------------------------------------------------------------------ def plotfootprints(image, segments, dimshape, wcs, skycoord, smoothing, diagplot, interactplot, verbose): # ------------------------------------------------------------------------------------------------------------------------ """ This function plots the footprints using matplotlib """ import matplotlib.pyplot as plt from astropy.visualization import SqrtStretch # Sqrt scale from astropy.visualization import LogStretch # Log scale from astropy.visualization.mpl_normalize import ImageNormalize norm = ImageNormalize(stretch=LogStretch()) print(whocalls(verbose=verbose) + "Plotting ...") if interactplot: plt.ion() #plt.clf() plt.title('Exposure map + footprint') Xlab = 'X ima' Ylab = 'Y ima' if skycoord and wcs['wcs_is_possible']: Xlab = 'RA' Ylab = 'Dec' plt.xlabel(Xlab) plt.ylabel(Ylab) # Load image try: ima = plt.imshow(image, origin='lower', cmap='hot', norm=norm) except IOError: error("\n"+whocalls(verbose=verbose) + "ERROR: the input image could not be opened!\n") print("===================================================================================================") import sys sys.exit(1) # Plot footprint polygons step = smoothing # Step for polygon smoothing for i in range(len(segments)): plot = plt.plot(segments[i][:,0][::step],segments[i][:,1][::step], 'c') # plot = plt.plot(segments[0][:,0], segments[0][:,1] , 'c') # step = 4 # X_47_steps = segments[47][:,0][::step] # Y_47_steps = segments[47][:,1][::step] # # plot = plt.plot(segments[47][:,0], segments[47][:,1] , 'ro') # plot = plt.plot(X_47_steps, Y_47_steps, 'bx') # # f, ax = plt.subplots(1, 1, sharex=True) # f.suptitle('Exposure map + footprint', fontsize=12) # ax.imshow(image, origin='lower', cmap='hot', norm=norm) # ax.plot(segments[0][:,0],segments[0][:,1]) # ax.plot(segments[47][:,0],segments[47][:,1]) # # Interactive plot in graphics window if interactplot: plt.show() r = '' r = input("Hit return to quit ") # Plot to PNG file else: plt.savefig( diagplot, format='png') # ------------------------------------------------------------------------------------------------------------------------ if __name__=='__main__': import sys main(" ".join(sys.argv[1:])) # ------------------------------------------------------------------------------------------------------------------------