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