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