Johan Bregeon / ctools

Commit 1f796dc19ca78de4ad8b9060e49d76df8480efc9

Authored by Jürgen Knödlseder
1 parent 5540ec74

Add tsmerge script written by Michael Mayer

Inline Side-by-side
Showing 2 changed files with 184 additions and 1 deletions
scripts/Makefile.am
  View file @1f796dc
1 1 ## Process this file with automake to produce Makefile.in
2 2  
3 3 # Scripts to be distributed in bindir
4   -binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py
  4 +binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py \
  5 + tsmerge.py
5 6  
6 7 # Scripts to be distributed in pythondir
7 8 pythonscripts = cspull.py cssens.py cstsdist.py cscaldb.py \
... ...
scripts/tsmerge.py 0 → 100644
  View file @1f796dc
  1 +#!/usr/bin/env python
  2 +
  3 +import gammalib
  4 +import sys
  5 +import os
  6 +import glob
  7 +import json
  8 +
  9 +def print_help():
  10 + print ""
  11 + print "tsmerge.py expects two keyword argumens: use via command line as follows"
  12 + print "tsmerge.py files=<list of files> outfile=<desired outfile name>"
  13 + print "-------------------------------------------------------------------------"
  14 + print "*** \"outfile\" should be a path to a fits file"
  15 + print "*** \"files\" should be either be:"
  16 + print " - an ascii-file listing the files your want to merge"
  17 + print " OR - a wildcard string describing the files for merging"
  18 + print "-------------------------------------------------------------------------"
  19 + print ""
  20 + print "Example: python tsmerge,py files=tsmap*.fits outfile=tsmap.fits"
  21 + print ""
  22 + sys.exit(-1)
  23 +
  24 +
  25 +# Maps class holding the different sky maps
  26 +# created by cttsmap
  27 +class maps:
  28 +
  29 + # Constructor
  30 + def __init__(self,fitsfile):
  31 +
  32 + self.filename = fitsfile
  33 + if not os.path.isfile(fitsfile):
  34 + self.error("Maps instance cannot be constructed: File \""+self.filename+"\" does not exist")
  35 +
  36 + # open FITS file
  37 + fits = gammalib.GFits(fitsfile)
  38 +
  39 + # get TS map
  40 + self.tsmap = gammalib.GSkymap()
  41 + self.tsmap.read(fits[0])
  42 +
  43 + # Get status map
  44 + if not fits.contains("STATUS MAP"):
  45 + self.error("No extension \"STATUS MAP\" found in \""+self.filename+"\"")
  46 + self.statusmap = gammalib.GSkymap()
  47 + self.statusmap.read(fits["STATUS MAP"])
  48 +
  49 + # Get other maps
  50 + self.maps = []
  51 + self.mapnames = []
  52 + for hdu in fits:
  53 + if hdu.extname() != "IMAGE" and hdu.extname() != "STATUS MAP":
  54 + skymap = gammalib.GSkymap()
  55 + skymap.read(hdu)
  56 + self.maps.append(skymap)
  57 + self.mapnames.append(hdu.extname())
  58 +
  59 + # len operator
  60 + def __len__(self):
  61 + return len(self.maps)
  62 +
  63 + # operator []
  64 + def __getitem__(self,index):
  65 + return self.maps[index]
  66 +
  67 + # Check the status map for completeness
  68 + def check_status(self):
  69 + check = True
  70 + for pix in self.statusmap:
  71 + if pix < 0.5:
  72 + check = False
  73 + break
  74 + return check
  75 +
  76 + # Save maps to outfile
  77 + def save(self,outfile):
  78 +
  79 + fits = gammalib.GFits()
  80 + self.tsmap.write(fits)
  81 +
  82 + for i in range(len(self.maps)):
  83 + self.maps[i].write(fits)
  84 +
  85 + for i in range(len(self.mapnames)):
  86 + fits[i+1].extname(self.mapnames[i])
  87 +
  88 + if not self.check_status():
  89 + self.statusmap.write(fits)
  90 + fits[fits.size()-1].extname("STATUS MAP")
  91 + fits.saveto(outfile,True)
  92 +
  93 + # Add map instances to each other
  94 + # in this way the TS maps get merged
  95 + def add(self,addmaps):
  96 + if not len(addmaps) == len(self.maps):
  97 + self.error("Maps cannot be added: Maps from \""+self.filename+"\" do not match from \""+addmaps.filename+"\"")
  98 +
  99 + # Loop over bins
  100 + for i in range(self.tsmap.npix()):
  101 + if addmaps.statusmap[i] > 0.5:
  102 + self.tsmap[i] = addmaps.tsmap[i]
  103 + self.statusmap[i] = addmaps.statusmap[i]
  104 + for j in range(len(self.maps)):
  105 + self.maps[j][i] = addmaps[j][i]
  106 + # Unify error messages
  107 + def error(self,message):
  108 + print "*** ERROR ***: "+message
  109 + sys.exit(-1)
  110 +
  111 +
  112 +def pars_from_argv(argv):
  113 +
  114 + if argv is None:
  115 + argv = sys.argv
  116 +
  117 + args = []
  118 + kwargs = {}
  119 + for s in argv[1:]:
  120 + if s.count('=') == 1:
  121 + key, value = s.split('=', 1)
  122 + else:
  123 + key, value = None, s
  124 + try:
  125 + value = json.loads(value)
  126 + except ValueError:
  127 + pass
  128 + if key:
  129 + kwargs[key] = value
  130 + else:
  131 + args.append(value)
  132 + return args, kwargs
  133 +
  134 +
  135 +def get_kwargs(argv):
  136 + kwargs = {}
  137 +
  138 + for s in argv:
  139 + if s.count('=') == 1:
  140 + key, value = s.split('=', 1)
  141 + kwargs[key] = value
  142 +
  143 + else:
  144 + print_help()
  145 +
  146 + return kwargs
  147 +
  148 +if __name__ == "__main__":
  149 + # Initialise list of files
  150 +
  151 + kwargs = get_kwargs(sys.argv[1:])
  152 + if not kwargs.has_key("outfile") or not kwargs.has_key("files"):
  153 + print_help()
  154 + outfile = kwargs["outfile"]
  155 + files = kwargs["files"]
  156 + allfiles = []
  157 + if os.path.isfile(files):
  158 + allfiles = open(files).readlines()
  159 + else:
  160 + allfiles = glob.glob(files)
  161 +
  162 + # Initialise files to be merged
  163 + mergefiles = []
  164 +
  165 + # Test files for the entry status map
  166 + # and add them to the file list
  167 + for fitsfile in allfiles:
  168 + fits = gammalib.GFits(fitsfile)
  169 + if fits.contains("STATUS MAP"):
  170 + mergefiles.append(fitsfile)
  171 + fits.close()
  172 +
  173 + # create a maps instance from the first file
  174 + finalmaps = maps(mergefiles[0])
  175 +
  176 + # add other files to the maps instance
  177 + for fitsfile in mergefiles[1:]:
  178 + finalmaps.add(maps(fitsfile))
  179 +
  180 + # save the filled maps to file
  181 + finalmaps.save(outfile=outfile)
  182 +
... ...