"""
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 = 'slewFootCalc'
VERSION = '1.1';
author = 'Ole Koenig, Jose Vicente Perea';
date = '2018-07-30';
ChangeLog
=========
Version 1.1 - 2018-07-30 (JVP)
------------
+ Adaptation for Pipeline
Version 1.0 - 2018-06-26 (Ole Koenig)
------------
+ Initial version
Task for Slew footprints
This code aims on determining the footprints of the XMM exposure maps
Author: Ole Koenig ole.koenig@fau.de 26/06/2018
"""
__version__ = "1.1"
__author__ = "Jose V. Perea"
__copyright__ = "ESA (C) 2000-2018"
__license__ = "GPL"
__progname__ = "slewFootCalc"
def usage():
'''
Help.
'''
print("""
usage: slewFootCalc.py --expMap=<...> --ds9footprint=<...>
This code aims on determining the footprints of the XMM exposure maps
Author: Ole Koenig ole.koenig@fau.de 26/06/2018
Basic idea of the code:
Axioms:
The XMM footprint consists of a minimum of 4 and a maximum of 6 points.
It has a clear rhombohedral shape with straight edges.
The four cornerpoints of the image frame are (clock-wise, in pixels):
C1 = bottom left (1,1)
C2 = upper left (1,max_y)
C3 = upper right (max_x,max_y)
C4 = bottom right (max_x,1)
1) Load the fits file
2) Detect NON-illuminated corner points (there are in total 4 corner points).
Reason: By walking the edges from the NON-illuminated corners the code avoids
landing in "gaps" or "stripes of 0 value" in the exposure map.
3) From the corner point walk in x and y direction through the arrays until
finding a value >0. This is the footprint point.
Repeat this in a clock-wise rotation starting in the lower left corner.
4) In total there should be a maximum of 6 footprint points (and minimum of 4).
If the code gives out less then 4 or more then 6 corner points it appends the
filename to "list_of_weird_images.txt"
5) Tranform pixels to WCS world coordinates and write to file
"xmm_expmaps_footprints.txt" in format filename | Polygon(ra1,dec1,ra2,dec2...)
###### DS9 output ######
Additional one can set the option --ds9 which creates a .reg file for ds9
displaying the footprint points as polygon.
########################
The code can be either used with an input file (list of filenames, delimiter=\n)
python3 slewFootCalc.py --file filelist.txt [--ds9]
or with
python3 slewFootCalc.py --file name.fits [--ds9]
which promts the output to stdout.
What the code is NOT capable of:
- Find footprint points in the center of the image (it only walks on the edges)
- Calculate if the shape is a correct parallelogonal hexagon
- handle wavy lines within the image: it assumes a strict rhombohedral shape
- handle gaps (only walks on the edges of the image!)
Note that the x,y pixel definition in fits files runs from 1..n (actually, the
only really reasonable choice!) and not from 0..n-1 as it is in python.
Modified by:
- if either of the input image dimensions is 1 (or 0), exits without producing footprint.
- if the code cannot find the footprints by stepping just use the four corners of the
input image
===================================================================================================
""")
def error(string):
'''
prints the string in red.
'''
sys.stdout.write( string )
def version(version_num):
'''
prints version number.
'''
versionString = '{0} ({0}-{1}) [-]'.format(sys.argv[0].split('/')[-1], version_num)
print(versionString)
def getFootprint(filename, outFootprint):
thr = 1
color_red = '\x1b[6;15;41m'
color_normal = '\x1b[0m'
from astropy.io import fits
from astropy.wcs import WCS
output, damaged = "",False
with fits.open (filename) as hdulist:
NAXIS1 = hdulist[0].header['NAXIS1']
NAXIS2 = hdulist[0].header['NAXIS2']
wcs = WCS(hdulist[0].header)
image = hdulist[0].data
if NAXIS1<2 or NAXIS2<2:
print("File " + filename + "has too few pixels in at least one dimension to "
+ "calculate a footprint")
damaged = True
return " ",damaged
C1, C2, C3, C4 = (0,0), (0,NAXIS2-1), (NAXIS1-1,NAXIS2-1), (NAXIS1-1,0)
startp = [C1,C2,C3,C4]
x_pixel_max = len (image[0]) -1
y_pixel_max = len (image) -1
non_ill, ill = get_start(image, x_pixel_max, y_pixel_max, C1, C2, C3, C4)
fpp=[]
for ii, p in enumerate(non_ill,start=1):
if p=="C1":
for i in range(0, x_pixel_max, thr):
if image[0][i] > 0.0:
fpp.append( (i,0) )
break
if C1 in ill and C1 not in fpp and ii==1:
fpp.append( C1 )
if p=="C1":
for i in range(0, y_pixel_max, thr):
if image[i][0] > 0.0:
fpp.append( (0,i) )
break
if p=="C2":
for i in reversed( range(0, y_pixel_max, thr)):
if image[i][0] > 0.0:
fpp.append( (0,i) )
break
if C2 in ill and C2 not in fpp and ii==2:
fpp.append( C2 )
if p=="C2":
for i in range(0, x_pixel_max, thr):
if image[-1][i] > 0.0:
fpp.append( (i,NAXIS2-1) )
break
if p=="C3":
for i in reversed( range(0, x_pixel_max, thr)):
if image[-1][i] > 0.0:
fpp.append( (i,NAXIS2-1) )
break
if C3 in ill and C3 not in fpp and ii==3:
fpp.append( C3 )
if p=="C3":
for i in reversed( range(0, y_pixel_max, thr)):
if image[i][-1] > 0.0:
fpp.append( (NAXIS1-1,i) )
break
if p=="C4":
for i in range(0, y_pixel_max, thr):
if image[i][-1] > 0.0:
fpp.append( (NAXIS1-1,i) )
break
if C4 in ill and C4 not in fpp and ii==4:
fpp.append( C4 )
if p=="C4":
for i in reversed( range(0, x_pixel_max, thr)):
if image[0][i] > 0.0:
fpp.append( (i,0) )
break
if len(fpp) >= 4:
world = wcs.all_pix2world(fpp,0)
else:
if len(fpp)==0:
print(" " + color_red + "WARNING" + color_normal + ": No corner is illuminated! using image corners")
world = wcs.all_pix2world(startp,0)
elif 06
print(" " + color_red + "WARNING" + color_normal + ": Number of footprints smaller then 4, using image corners.")
world = wcs.all_pix2world(startp,0)
writeds9output( outFootprint, world )
return output, damaged
def get_start (image, x_pixel_max, y_pixel_max, C1, C2, C3, C4):
footprintpoints = [C1,C2,C3,C4]
non_ill = ["","","",""]
if image[0][0]==0.0:
footprintpoints.remove( C1 )
non_ill[0] = "C1"
if image[y_pixel_max][0]==0.0:
footprintpoints.remove( C2 )
non_ill[1] = "C2"
if image[y_pixel_max][x_pixel_max]==0.0:
footprintpoints.remove( C3 )
non_ill[2] = "C3"
if image[0][x_pixel_max]==0.0:
footprintpoints.remove( C4 )
non_ill[3] = "C4"
return non_ill, footprintpoints
def writeds9output (outFootprint, fpp):
import os
print('{0}: Slew footprint ds9 file : {1}'.format(__progname__, outFootprint))
fp = open (outFootprint, 'w')
fp.write("# Region file format: DS9 version 4.0\n# Filename: %s\n" % outFootprint)
fp.write("# global color=green font='helvetica 12 bold' select=1 highlite=1 edit=1 move=1 delete=1 include=1 fixed=0 source\n")
fp.write("linear\npolygon(")
for i in range (len(fpp)):
fp.write ("%f, %f," % (fpp[i][0], fpp[i][1]))
cur_pos = fp.tell()
fp.seek (cur_pos-1, os.SEEK_SET)
fp.truncate()
fp.write(") # color=green width=2")
"""
# Circles
fp.write("FK5\n")
for i in range( len(fpp)):
fp.write ("Circle(%f, %f, 0.005)\n" % (fpp[i][0], fpp[i][1])) # x,y (pixels)
"""
fp.close()
def main(arguments):
'''
This is the main function doing the full loop.
'''
import getopt, os
import ModuleUtil
version_num = __version__
verbose = False
expMap = ''
outFootprint = 'expMapfootprint.ds9'
arguments = arguments.split()
progName = arguments[0]
argsList = arguments[1:]
try:
opts, args = getopt.getopt(argsList, "vo:Vh", ["expMap=", "outFootprint=", \
"version", "verbose=", "help"])
except getopt.GetoptError as err:
usage()
error("ERROR: Could not understand the options!\n\n")
print(err)
print("================================================================================================")
return
for o, a in opts:
if o in ("-v", "--version"):
version(version_num)
return 0
if o in ("-V", "--verbose"):
verbose = ModuleUtil.varCheck(client = progName, inVar = a).checkBool()
if o in ("-h", "--help"):
usage()
return 0
if o in ("--expMap"):
expMap = a
if o in ("-o", "--outFootprint"):
outFootprint = a
if os.path.isfile( expMap ) and os.access(expMap, os.R_OK):
print('{0}: Creating footPrint for Slew exp map : {1}'.format(__progname__, expMap))
output, damaged = getFootprint( expMap, outFootprint )
else:
usage()
error("ERROR: Input exposure map does not exist\n\n")
return
if __name__=='__main__':
import sys
main(" ".join(sys.argv))