diff --git a/ChangeLog b/ChangeLog index c55a3e2..dc8250b 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2016-11-25 +2016-12-01 * Version 1.2.0 released ======================== + Add cssrcdetect script (#1617) Add background subtraction to ctskymap tool (#1882) Add csobsselect script (#1881) Add csmodelselect script (#1879) diff --git a/NEWS b/NEWS index e556c2c..5a73fd4 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.2.0 -25 November 2016 +1 December 2016 Introduction @@ -206,30 +206,33 @@ The class now derives from ctlikelihood (#1828). cscripts -------- -Added the support modules "modutils" and "ioutils" to the cscripts module and -removed the obsolete functions fit(), cterror(), cntmap() and modmap() from the +Add the support modules "modutils" and "ioutils" to the cscripts module and +remove the obsolete functions fit(), cterror(), cntmap() and modmap() from the "obsutils" module (#1830). The csroot2caldb script now takes the zenith and azimuth angles as parameters and writes them into the calibration database (#1835). -Added a csviscube script that computes the visibility cube providing the number +Add csviscube script that computes the visibility cube providing the number of hours that a given position on the sky is observable under a given zenith angle by an IACT array during a given time period. The tool implements so far a Sun constraint, but no Moon constraint (#1835). -A publish() method has been added to cslightcrv and csspec that publishes the -results on a VO hub (#1871). +Add publish() method to cslightcrv and csspec that publishes the results on a +VO hub (#1871). Fix energy boundary handling in cssens that has been corrupted during code refactoring (#1873). -Added csmodelselect script that allows selecting models from a model definition +Add csmodelselect script that allows selecting models from a model definition XML file (#1879). -Added csobsselect script that allows selecting observations from an observation +Add csobsselect script that allows selecting observations from an observation definition XML file (#1881). +Add cssrcdetect script to detect sources in a sky map (#1617). + + Examples -------- diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index 8cc586c..685ec45 100644 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -36,6 +36,7 @@ pythonscripts = $(builddir)/__init__.py \ $(srcdir)/csroot2caldb.py \ $(srcdir)/cssens.py \ $(srcdir)/csspec.py \ + $(srcdir)/cssrcdetect.py \ $(srcdir)/cstsdist.py \ $(srcdir)/cstsmapmerge.py \ $(srcdir)/cstsmapsplit.py \ @@ -65,6 +66,7 @@ parfiles = $(srcdir)/cscaldb.par \ $(srcdir)/csroot2caldb.par \ $(srcdir)/cssens.par \ $(srcdir)/csspec.par \ + $(srcdir)/cssrcdetect.par \ $(srcdir)/cstsdist.par \ $(srcdir)/cstsmapmerge.par \ $(srcdir)/cstsmapsplit.par \ @@ -91,6 +93,7 @@ linkscripts = cscaldb \ csroot2caldb \ cssens \ csspec \ + cssrcdetect \ cstsdist \ cstsmapmerge \ cstsmapsplit \ @@ -119,6 +122,7 @@ testscripts = $(srcdir)/tests/__init__.py \ $(top_srcdir)/test/test_csroot2caldb.py \ $(top_srcdir)/test/test_cssens.py \ $(top_srcdir)/test/test_csspec.py \ + $(top_srcdir)/test/test_cssrcdetect.py \ $(top_srcdir)/test/test_cstsdist.py \ $(top_srcdir)/test/test_cstsmapmerge.py \ $(top_srcdir)/test/test_cstsmapsplit.py \ @@ -160,6 +164,7 @@ testdata = $(top_srcdir)/test/data/crab.xml \ $(top_srcdir)/test/data/pull.dat \ $(top_srcdir)/test/data/sensitivity.dat \ $(top_srcdir)/test/data/sensitivity_integral.dat \ + $(top_srcdir)/test/data/skymap.fits \ $(top_srcdir)/test/data/spectrum.fits \ $(top_srcdir)/test/data/ts.dat \ $(top_srcdir)/test/data/models_spatial.xml \ diff --git a/cscripts/__init__.py.in b/cscripts/__init__.py.in index 9332118..34c92d1 100644 --- a/cscripts/__init__.py.in +++ b/cscripts/__init__.py.in @@ -39,6 +39,7 @@ __all__ = [ "csroot2caldb", "cssens", "csspec", + "cssrcdetect", "cstsdist", "cstsmapmerge", "cstsmapsplit", @@ -71,6 +72,7 @@ from cscripts.csresmap import csresmap from cscripts.csroot2caldb import csroot2caldb from cscripts.cssens import cssens from cscripts.csspec import csspec +from cscripts.cssrcdetect import cssrcdetect from cscripts.cstsdist import cstsdist from cscripts.cstsmapmerge import cstsmapmerge from cscripts.cstsmapsplit import cstsmapsplit diff --git a/cscripts/csmodelinfo.py b/cscripts/csmodelinfo.py index 088e0cc..ed2bf22 100755 --- a/cscripts/csmodelinfo.py +++ b/cscripts/csmodelinfo.py @@ -21,6 +21,7 @@ import sys import gammalib import ctools +from cscripts import modutils # ================= # @@ -86,155 +87,6 @@ class csmodelinfo(ctools.cscript): # Return return - def _ds9attributes(self, model): - """ - Returns attributes for DS9 string - Parameters - ---------- - model : `~gammalib.GModel` - Model - - Returns - ------- - attributes : str - DS9 attribute string - """ - # Determine region color. A model where all parameters are fixed will - # get the color defined by "fixed_color", a model with at least one - # free parameter will get the color defined by "free_color". - color = self['fixed_color'].string() - for par in model: - if par.is_free(): - color = self['free_color'].string() - break - - # Set DS9 attributes - attributes = (' color=%s width=%d font=%s %d %s %s' % - (color, self['width'].integer(), - self['fontfamily'].string(), self['fontsize'].integer(), - self['fontweight'].string(), self['fontslant'].string())) - - # Return attributes - return attributes - - def _model2ds9string(self, model): - """ - Converts model into a DS9 region string - - Parameters - ---------- - model : `~gammalib.GModel` - Model - - Returns - ------- - ds9string : str - DS9 region string - """ - # Initialise DS9 region string - ds9string = '' - - # Retrieve model sky direction. The model is skipped in case that - # it does not provide a sky direction. - is_valid = True - try: - modelpos = model.spatial().dir() - except AttributeError: - msg = ('Skip model "%s" since it has no sky direction.\n' % - model.name()) - self._log(msg) - is_valid = False - - # Continue only if sky direction was found - if is_valid: - - # Retrieve model type and name - modeltype = model.type() - modelname = model.name() - - # Handle point source - if modeltype == 'PointSource': - - # Append point with Right Ascension and Declination to the DS9 - # string. The type of the point is specified by the "pnt_type" - # parameter, the size of the point by the "pnt_mark_size" - # parameter - ds9string += ('point(%.6f,%.6f) # point=%s %d' % - (modelpos.ra_deg(), modelpos.dec_deg(), - self['pnt_type'].string(), - self['pnt_mark_size'].integer())) - - # Handle extended sources - elif modeltype == "ExtendedSource": - - # Retrieve spatial model - spatial = model.spatial() - classname = spatial.classname() - - # Handle radial sources - if 'Radial' in classname: - - # Retrieve short name of model class (e.g. "Gauss" - # or "Disk) - shorttype = classname.split('Radial')[-1] - - # Handle Disk and Shell model - if (classname == 'GModelSpatialRadialDisk' or - classname == 'GModelSpatialRadialShell'): - size = spatial.radius() - - # Handle Gauss Model - elif classname == 'GModelSpatialRadialGauss': - size = spatial.sigma() - - # Skip if source is unknown - else: - msg = ('Skip model "%s" since the radial model "%s" ' - 'is unknown.\n' % (model.name(), classname)) - self._log(msg) - is_valid = False - - # Append circle to DS9 string - if is_valid: - ds9string += ('circle(%.6f,%.6f,%.6f) #' % - (modelpos.ra_deg(), modelpos.dec_deg(), - size*3600.0)) - - # Handle elliptical sources - elif 'Elliptical' in classname: - - # Retrieve short name and source size - shorttype = classname.split('Elliptical')[-1] - size1 = spatial.semimajor() - size2 = spatial.semiminor() - angle = spatial.posangle() - - # Append ellipse to DS9 string - ds9string += ('ellipse(%.6f,%.6f,%.6f,%.6f,%.6f) #' % - (modelpos.ra_deg(), modelpos.dec_deg(), - size1*3600.0, size2*3600.0, angle+90.0)) - - # Skip if source is neither radial nor elliptical - else: - msg = ('Skip model "%s" since the model "%s" is neither ' - 'a point source, a radial source, nor an elliptical ' - 'source.\n' % (model.name(), classname)) - self._log(msg) - is_valid = False - - # Add short model type to modelname - if self['show_ext_type'].boolean(): - modelname +=' ('+shorttype+')' - - # Add DS9 attributes - if is_valid: - ds9string += self._ds9attributes(model) - if self['show_labels'].boolean(): - ds9string += ' text={'+modelname+'}' - - # Return string - return ds9string - # Public methods def run(self): @@ -407,36 +259,34 @@ class csmodelinfo(ctools.cscript): # Log filename self._log_value(gammalib.NORMAL, 'DS9 filename', ds9file.url()) - - # Open file - f = open(ds9file.url(), 'w') - - # Write coordinate system - f.write('fk5\n') - - # Loop over models - for model in self._models: - - # Continue only if point source or extended source model - if (model.type() == 'PointSource' or - model.type() == 'ExtendedSource'): - line = self._model2ds9string(model) - if len(line): - f.write(line+'\n') - - # Logging for diffuse components - elif model.type() == 'DiffuseSource': - self._log('Skipping diffuse model "'+model.name()+'".\n') - - # Logging for background components - else: - if self._logExplicit(): - self._log('Skipping background model "'+model.name()+ - '".\n') - - # Close file - f.close() - + + # Get ds9 parameters + pnt_type = self['pnt_type'].string() + pnt_mark_size = self['pnt_mark_size'].integer() + free_color = self['free_color'].string() + fixed_color = self['fixed_color'].string() + width = self['width'].integer() + fontfamily = self['fontfamily'].string() + fontsize = self['fontsize'].integer() + fontweight = self['fontweight'].string() + fontslant = self['fontslant'].string() + + # Save DS9 region file + errors = modutils.models2ds9file(self._models, ds9file.url(), + pnt_type=pnt_type, + pnt_mark_size=pnt_mark_size, + free_color=free_color, + fixed_color=fixed_color, + width=width, + fontfamily=fontfamily, + fontsize=fontsize, + fontweight=fontweight, + fontslant=fontslant) + + # Log errors + if len(errors): + self._log(errors) + # Return return diff --git a/cscripts/cssrcdetect.par b/cscripts/cssrcdetect.par new file mode 100644 index 0000000..88bda7c --- /dev/null +++ b/cscripts/cssrcdetect.par @@ -0,0 +1,41 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: cssrcdetect.par # +# # +# Author: Juergen Knoedlseder # +# jknodlseder@irap.omp.eu # +# IRAP, Toulouse, France # +# # +# Purpose: Parameter file for the cssrcdetect script # +# # +############################################################################# +# +# Input and output parameters +#============================ +inmap, f, a, skymap.fits,,, "Input sky map file" +outmodel, f, a, models.xml,,, "Output model definition XML file" +outds9file, f, a, ds9.reg,,, "Output DS9 region file" + +# +# Script parameters +#================== +srcmodel, s, a, POINT,POINT,, "Source model type" +bkgmodel, s, a, NONE,NONE|IRF|AEFF|CUBE,, "Background model type" +threshold, r, a, 5.0,,, "Detection threshold (Gaussian sigma)" +maxsrcs, i, h, 20,1,1000, "Maximum number of sources" +exclrad, r, h, 0.2,,, "Exclusion radius (deg)" +fit_pos, b, h, yes,,, "Fit source positions?" +fit_shape, b, h, yes,,, "Fit source shapes?" + +# +# Standard parameters +#==================== +chatter, i, h, 2,0,4, "Chattiness of output" +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" +debug, b, h, no,,, "Debugging mode activated" +mode, s, h, ql,,, "Mode of automatic parameters" +logfile, f, h, cssrcdetect.log,,, "Log filename" diff --git a/cscripts/cssrcdetect.py b/cscripts/cssrcdetect.py new file mode 100755 index 0000000..f053ddd --- /dev/null +++ b/cscripts/cssrcdetect.py @@ -0,0 +1,465 @@ +#!/usr/bin/env python +# ========================================================================== +# Detects sources in a sky map +# +# Copyright (C) 2016 Juergen Knoedlseder +# +# This program 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. +# +# This program 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 this program. If not, see . +# +# ========================================================================== +import sys +import math +import gammalib +import ctools +from cscripts import modutils + + +# ================= # +# cssrcdetect class # +# ================= # +class cssrcdetect(ctools.cscript): + """ + Detects sources in a sky map + """ + + # Constructor + def __init__(self, *argv): + """ + Constructor + """ + # Set name + self._name = 'cssrcdetect' + self._version = '1.2.0' + + # Set protected members + self._models = gammalib.GModels() + + # Initialise application by calling the appropriate class constructor + self._init_cscript(argv) + + # Return + return + + + # Private methods + def _get_parameters(self): + """ + Get parameters from parfile + """ + # Query input parameters + self['inmap'].filename() + + # Query further parameters + self['srcmodel'].string() + self['bkgmodel'].string() + self['threshold'].real() + self['maxsrcs'].integer() + self['exclrad'].real() + self['fit_pos'].boolean() + self['fit_shape'].boolean() + + # Query ahead output model filename + if self._read_ahead(): + self['outmodel'].filename() + self['outds9file'].filename() + + # Write input parameters into logger + self._log_parameters(gammalib.TERSE) + + # Return + return + + def _detect_sources(self, counts): + """ + Detect sources in counts map + + Parameters + ---------- + counts : `~gammalib.GSkyMap()` + Counts map + """ + # Source detection iterations + for i in range(self['maxsrcs'].integer()): + + # Write header + self._log_header3(gammalib.NORMAL, 'Iteration '+str(i+1)) + + # Get map moments + mean, std = self._map_moments(counts) + + # Log map moments + self._log_value(gammalib.NORMAL, 'Map mean', mean) + self._log_value(gammalib.NORMAL, 'Map standard deviation', std) + + # Compute threshold + threshold = mean + self['threshold'].real() * std + + # Log map threshold + self._log_value(gammalib.NORMAL, 'Map threshold', threshold) + + # Get maximum value and corresponding sky direction + value, dir = self._find_maximum(counts, threshold=threshold) + + # If maximum found then log maximum and add model + if dir is not None: + + # Set source name + name = 'Src%3.3d' % (i+1) + + # Log maximum + self._log_value(gammalib.NORMAL, 'Map maximum', str(value)) + self._log_value(gammalib.NORMAL, name+' position', str(dir)) + + # Add model + self._add_model(dir, name) + + # Remove maximum from map + counts = self._remove_maximum(counts, dir, mean, + radius=self['exclrad'].real()) + + # ... otherwise log that no maximum was found and break iterations + else: + + # Log than no maximum was found + self._log_value(gammalib.NORMAL, 'Map maximum', + 'None above threshold') + + # Break + break + + # Return + return + + def _find_maximum(self, map, threshold=0.0): + """ + Find maximum in a sky map + + Parameters + ---------- + map : `~gammalib.GSkyMap()` + Sky map + threshold : float, optional + Threshold for maximum value + + Returns + ------- + value, dir : tuple of float and `~gammalib.GSkyDir()` + Maximum sky map value and corresponding sky direction + """ + # Initialise maximum pixel value and sky direction + value = threshold + dir = None + + # Loop over all pixels and find maximum + for i in range(map.npix()): + if map[i] > value: + value = map[i] + dir = map.inx2dir(i) + + # Return sky direction of maximum + return value, dir + + def _remove_maximum(self, map, dir, value=0.0, radius=0.1): + """ + Remove maximum from sky map by replacing pixels with a given value + + Parameters + ---------- + map : `~gammalib.GSkyMap()` + Sky map + dir : `~gammalib.GSkyDir()` + Sky direction of maximum + value : float, optional + Replacement value + radius : float, optional + Radius within which pixel values are replaced + + Returns + ------- + map : `~gammalib.GSkyMap()` + Sky map with maximum removed + """ + # Copy skymap + map_copy = map.copy() + + # Loop over all pixels and find maximum + for i in range(map_copy.npix()): + mapdir = map_copy.inx2dir(i) + if mapdir.dist_deg(dir) < radius: + map_copy[i] = value + + # Return copy of map + return map_copy + + def _map_moments(self, map): + """ + Determine moments of sky map pixels + + Parameters + ---------- + map : `~gammalib.GSkyMap()` + Sky map + + Returns + ------- + mean, std : tuple of float + Mean and standard deviation of pixel values + """ + # Initialise mean and standard deviation + mean = 0.0 + std = 0.0 + + # Loop over all pixels + for i in range(map.npix()): + mean += map[i] + std += map[i]*map[i] + + # Compute mean and standard deviation + mean /= float(map.npix()) + std /= float(map.npix()) + std -= mean + std = math.sqrt(std) + + # Return mean and standard deviation + return mean, std + + def _add_model(self, dir, name): + """ + Add model to model container + + Parameters + ---------- + dir : `~gammalib.GSkyDir()` + Sky direction of model + name : str + Model name + """ + # Set point source model + model = self._set_ptsrc(dir) + + # Set model name + model.name(name) + + # Append model to container + self._models.append(model) + + # Return + return + + def _set_bkg(self, type): + """ + Set background model + + Parameters + ---------- + type : str + Model type ('IRF', 'AEFF' or 'CUBE') + + Returns + ------- + model : `~gammalib.GModelData()` + Background model + """ + # Set spectral component + spectral = gammalib.GModelSpectralPlaw(1.0, 0.0, + gammalib.GEnergy(1.0, 'TeV')) + + # Set background model + if type == 'IRF': + model = gammalib.GCTAModelIrfBackground(spectral) + elif type == 'AEFF': + model = gammalib.GCTAModelAeffBackground(spectral) + elif type == 'CUBE': + model = gammalib.GCTAModelCubeBackground(spectral) + else: + model = None + + # Set background name + if model is not None: + model.name('Background') + + # Return model + return model + + def _set_ptsrc(self, dir): + """ + Set point source model + + Parameters + ---------- + dir : `~gammalib.GSkyDir()` + Sky direction of model + + Returns + ------- + model : `~gammalib.GModelSky()` + Point source model + """ + # Set spatial component + spatial = gammalib.GModelSpatialPointSource(dir) + + # Get fit parameters + fit_pos = self['fit_pos'].boolean() + fit_shape = self['fit_shape'].boolean() + + # Loop over all parameters + for par in spatial: + + # Handle position parameters + if par.name() == 'RA' or par.name() == 'DEC': + if fit_pos: + par.free() + else: + par.fix() + + # Set spectral component (powerlaw with 1% of Crab) + spectral = gammalib.GModelSpectralPlaw(5.7e-18, -2.48, + gammalib.GEnergy(0.3, 'TeV')) + + # Set sky model + model = gammalib.GModelSky(spatial, spectral) + + # Return model + return model + + + # Public methods + def run(self): + """ + Run the script + """ + # Clear model container + self._models.clear() + + # Switch screen logging on in debug mode + if self._logDebug(): + self._log.cout(True) + + # Get parameters + self._get_parameters() + + # Write header + self._log_header1(gammalib.NORMAL, 'Source detection') + + # Get skymap filename + inmap = self['inmap'].filename() + + # Open sky map file + fits = gammalib.GFits(inmap) + + # Extract primary extension as counts map + counts = gammalib.GSkyMap(fits.image(0)) + + # Detect sources + self._detect_sources(counts) + + # Close sky map file + fits.close() + + # Write detected sources as models into logger + self._log_models(gammalib.NORMAL, self._models, 'Detected source model') + + # Write header + self._log_header1(gammalib.NORMAL, 'Add background model') + + # Add background model if requested + bkgmodel = self._set_bkg(self['bkgmodel'].string()) + if bkgmodel is not None: + + # Append model + self._models.append(bkgmodel) + + # Log model + self._log_string(gammalib.NORMAL, str(bkgmodel)) + + # ... otherwise notify that no background model is added + else: + self._log_value(gammalib.NORMAL, 'Background model', 'None') + + # Return + return + + def save(self): + """ + Save sources into model definition XML file and ds9 file + """ + # Write header + self._log_header1(gammalib.TERSE, 'Save sources') + + # Get output filenames + outmodel = self['outmodel'].filename() + outds9file = self['outds9file'].filename() + + # If model definition XML file exists and clobber flag is false then + # raise an exception + if outmodel.exists() and not self._clobber: + msg = ('Cannot save "'+outmodel.url()+'": File already exists. ' + 'Use parameter clobber=yes to allow overwriting of files.') + raise RuntimeError(msg) + + # If ds9 file exists and clobber flag is false then raise an exception + elif outds9file.exists() and not self._clobber: + msg = ('Cannot save "'+outds9file.url()+'": File already exists. ' + 'Use parameter clobber=yes to allow overwriting of files.') + raise RuntimeError(msg) + + # Otherwise log filename and save file + else: + # Log model definition XML filename + self._log_value(gammalib.NORMAL, 'Model definition XML file', + outmodel.url()) + + # Save models + self._models.save(outmodel) + + # Log ds9 filename + self._log_value(gammalib.NORMAL, 'DS9 region file', + outds9file.url()) + + # Save models as ds9 region file + modutils.models2ds9file(self._models, outds9file.url()) + + # Return + return + + def execute(self): + """ + Execute the script + """ + # Open logfile + self.logFileOpen() + + # Read ahead output parameters + self._read_ahead(True) + + # Run the script + self.run() + + # Save observation definition file + self.save() + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + + # Create instance of application + app = cssrcdetect(sys.argv) + + # Execute application + app.execute() diff --git a/cscripts/modutils.py b/cscripts/modutils.py index 2e05cfe..19b6d95 100755 --- a/cscripts/modutils.py +++ b/cscripts/modutils.py @@ -99,3 +99,298 @@ def test_source(models, srcname, ra=None, dec=None, fitspat=False, # Return model container return outmodels + +# ================================= # +# Returns attributes for ds9 string # +# ================================= # +def ds9attributes(model, free_color='green', fixed_color='magenta', + width=2, fontfamily='helvetica', fontsize=12, + fontweight='normal', fontslant='roman'): + """ + Returns attributes for DS9 string + + Parameters + ---------- + model : `~gammalib.GModel` + Model + free_color : str, optional + Color for sources with free parameters (any DS9 color or hex code) + fixed_color : str, optional + Color for source without free parameters (any DS9 color or hex code) + width : int, optional + Line width for regions + fontfamily : str, optional + Font family for source labels (helvetica,times,courier) + fontsize : int, optional + Font size for source labels + fontweight : str, optional + Font weight for source labels (normal,bold) + fontslant : str, optional + Font slant for source labels (roman,italic) + + Returns + ------- + attributes : str + DS9 attribute string + """ + # Determine region color. A model where all parameters are fixed will + # get the color defined by "fixed_color", a model with at least one + # free parameter will get the color defined by "free_color". + color = fixed_color + for par in model: + if par.is_free(): + color = free_color + break + + # Set DS9 attributes + attributes = (' color=%s width=%d font="%s %d %s %s"' % + (color, width, fontfamily, fontsize, fontweight, fontslant)) + + # Return attributes + return attributes + + +# =========================== # +# Convert model to ds9 string # +# =========================== # +def model2ds9string(model, pnt_type='cross', pnt_mark_size=12, + show_labels=True, show_ext_type=True, + free_color='green', fixed_color='magenta', + width=2, fontfamily='helvetica', fontsize=12, + fontweight='normal', fontslant='roman'): + """ + Converts model into a DS9 region string + + Parameters + ---------- + model : `~gammalib.GModel` + Model + pnt_type : str, optional + Marker type for point sources (circle,box,diamond,cross,x,arrow,boxcircle) + pnt_mark_size : integer, optional + Marker size for point sources + show_labels : boolean, optional + Add source labels? + show_ext_type : boolean, optional + Show type of extended model in source name? + free_color : str, optional + Color for sources with free parameters (any DS9 color or hex code) + fixed_color : str, optional + Color for source without free parameters (any DS9 color or hex code) + width : int, optional + Line width for regions + fontfamily : str, optional + Font family for source labels (helvetica,times,courier) + fontsize : int, optional + Font size for source labels + fontweight : str, optional + Font weight for source labels (normal,bold) + fontslant : str, optional + Font slant for source labels (roman,italic) + + Returns + ------- + ds9string : str + DS9 region string + msg : str + Error message + """ + # Initialise DS9 region string and message string + ds9string = '' + msg = '' + + # Retrieve model sky direction. The model is skipped in case that + # it does not provide a sky direction. + is_valid = True + try: + modelpos = model.spatial().dir() + except AttributeError: + msg = ('Skip model "%s" since it has no sky direction.\n' % + model.name()) + is_valid = False + + # Continue only if sky direction was found + if is_valid: + + # Retrieve model type and name + modeltype = model.type() + modelname = model.name() + + # Handle point source + if modeltype == 'PointSource': + + # Append point with Right Ascension and Declination to the DS9 + # string. The type of the point is specified by the "pnt_type" + # parameter, the size of the point by the "pnt_mark_size" + # parameter + ds9string += ('point(%.6f,%.6f) # point=%s %d' % + (modelpos.ra_deg(), modelpos.dec_deg(), + pnt_type, pnt_mark_size)) + + # Handle extended sources + elif modeltype == "ExtendedSource": + + # Retrieve spatial model + spatial = model.spatial() + classname = spatial.classname() + + # Handle radial sources + if 'Radial' in classname: + + # Retrieve short name of model class (e.g. "Gauss" + # or "Disk) + shorttype = classname.split('Radial')[-1] + + # Handle Disk and Shell model + if (classname == 'GModelSpatialRadialDisk' or + classname == 'GModelSpatialRadialShell'): + size = spatial.radius() + + # Handle Gauss Model + elif classname == 'GModelSpatialRadialGauss': + size = spatial.sigma() + + # Skip if source is unknown + else: + msg = ('Skip model "%s" since the radial model "%s" ' + 'is unknown.\n' % (model.name(), classname)) + is_valid = False + + # Append circle to DS9 string + if is_valid: + ds9string += ('circle(%.6f,%.6f,%.6f) #' % + (modelpos.ra_deg(), modelpos.dec_deg(), + size*3600.0)) + + # Handle elliptical sources + elif 'Elliptical' in classname: + + # Retrieve short name and source size + shorttype = classname.split('Elliptical')[-1] + size1 = spatial.semimajor() + size2 = spatial.semiminor() + angle = spatial.posangle() + + # Append ellipse to DS9 string + ds9string += ('ellipse(%.6f,%.6f,%.6f,%.6f,%.6f) #' % + (modelpos.ra_deg(), modelpos.dec_deg(), + size1*3600.0, size2*3600.0, angle+90.0)) + + # Skip if source is neither radial nor elliptical + else: + msg = ('Skip model "%s" since the model "%s" is neither ' + 'a point source, a radial source, nor an elliptical ' + 'source.\n' % (model.name(), classname)) + is_valid = False + + # Add short model type to modelname + if show_ext_type: + modelname +=' ('+shorttype+')' + + # Add DS9 attributes + if is_valid: + ds9string += ds9attributes(model, free_color=free_color, + fixed_color=fixed_color, + width=width, + fontfamily=fontfamily, + fontsize=fontsize, + fontweight=fontweight, + fontslant=fontslant) + if show_labels: + ds9string += ' text={'+modelname+'}' + else: + ds9string = '' + + # Return ds9 and message strings + return ds9string, msg + + +# ========================= # +# Save models into ds9 file # +# ========================= # +def models2ds9file(models, filename, + pnt_type='cross', pnt_mark_size=12, + show_labels=True, show_ext_type=True, + free_color='green', fixed_color='magenta', + width=2, fontfamily='helvetica', fontsize=12, + fontweight='normal', fontslant='roman'): + """ + Save models in a DS9 region file + + Parameters + ---------- + model : `~gammalib.GModel` + Model + filename : str + DS9 region file name + pnt_type : str, optional + Marker type for point sources (circle,box,diamond,cross,x,arrow,boxcircle) + pnt_mark_size : integer, optional + Marker size for point sources + show_labels : boolean, optional + Add source labels? + show_ext_type : boolean, optional + Show type of extended model in source name? + free_color : str, optional + Color for sources with free parameters (any DS9 color or hex code) + fixed_color : str, optional + Color for source without free parameters (any DS9 color or hex code) + width : int, optional + Line width for regions + fontfamily : str, optional + Font family for source labels (helvetica,times,courier) + fontsize : int, optional + Font size for source labels + fontweight : str, optional + Font weight for source labels (normal,bold) + fontslant : str, optional + Font slant for source labels (roman,italic) + + Returns + ------- + errors : str + Error message + """ + # Initialise error string + errors = '' + + # Open file + f = open(filename, 'w') + + # Write coordinate system + f.write('fk5\n') + + # Loop over models + for model in models: + + # Continue only if point source or extended source model + if (model.type() == 'PointSource' or + model.type() == 'ExtendedSource'): + line, msg = model2ds9string(model, + pnt_type=pnt_type, + pnt_mark_size=pnt_mark_size, + free_color=free_color, + fixed_color=fixed_color, + width=width, + fontfamily=fontfamily, + fontsize=fontsize, + fontweight=fontweight, + fontslant=fontslant) + if len(line): + f.write(line+'\n') + elif len(msg): + errors += msg+'\n' + + # Logging for diffuse components + elif model.type() == 'DiffuseSource': + errors += 'Skipping diffuse model "'+model.name()+'".\n' + + # Logging for background components + else: + errors += 'Skipping background model "'+model.name()+'".\n' + + # Close file + f.close() + + # Return error message string + return errors diff --git a/doc/source/users/reference_manual/cssrcdetect.rst b/doc/source/users/reference_manual/cssrcdetect.rst new file mode 100644 index 0000000..6c3035f --- /dev/null +++ b/doc/source/users/reference_manual/cssrcdetect.rst @@ -0,0 +1,97 @@ +.. _cssrcdetect: + +cssrcdetect +=========== + +Detects sources in sky map and writes results into a model definition XML file. + + +Synopsis +-------- + +This script detects sources in a sky map using a peak detection method and +writes the detected sources into a model definition XML file. In addition, +the script produces also a DS9 region file of the detected sources that can +be overlayed over a sky map displayed with DS9. + +For the moment, the script writes out all sources as point sources. Optionally, +a background model can be added to the model definition XML file so that the +output file can be used directly as input of a :doc:`ctlike` run. + + +General parameters +------------------ + +``inmap [file]`` + Input sky map file. + +``outmodel [file]`` + Output model definition XML file. + +``outds9file [file]`` + Output DS9 region file. + +``srcmodel [string]`` + Source model type. For the moment, all detected sources will be added as + point sources with power law spectral shapes to the model definition + XML file. + +``bkgmodel [string]`` + Background model type. Using ``NONE`` no background model will be added + to the model definition XML file. Using ``IRF`` a background model based + on the template information in the Instrument Response Function will be + added. Using ``AEFF`` a background model based on the shape of the effective + area will be added. Using ``CUBE`` a background model for stacked analysis + will be added. + +``threshold [real]`` + Detection threshold (Gaussian sigma). Only sources above this detection + threshold will be added to the model definiton XML file. + +``(maxsrcs = 20) [integer]`` + Maximum number of sources that should be detected. + +``(exclrad = 0.2) [real]`` + Radius around a detected source that is excluded from further source + detection (degrees). + +``(fit_pos = yes) [boolean]`` + Enable fitting of source positions in model definiton XML file? + +``(fit_shape = yes) [boolean]`` + Enable fitting of source shapes in model definiton XML file? + + +Standard parameters +------------------- + +``(chatter = 2) [integer]`` + Verbosity of the executable: + chatter = 0: no information will be logged + + chatter = 1: only errors will be logged + + chatter = 2: errors and actions will be logged + + chatter = 3: report about the task execution + + chatter = 4: detailed report about the task execution + +``(clobber = yes) [boolean]`` + Specifies whether an existing output model file should be overwritten. + +``(debug = no) [boolean]`` + Enables debug mode. In debug mode the executable will dump any log file + output to the console. + +``(mode = ql) [string]`` + Mode of automatic parameters (default is "ql", i.e. "query and learn"). + +``(logfile = cssrcdetect.log) [filename]`` + Log filename. + + +Related tools or scripts +------------------------ + +None diff --git a/doc/source/users/reference_manual/index.rst b/doc/source/users/reference_manual/index.rst index 8ad2850..94b179e 100644 --- a/doc/source/users/reference_manual/index.rst +++ b/doc/source/users/reference_manual/index.rst @@ -40,7 +40,6 @@ cscripts :maxdepth: 1 cscaldb --- Lists available instrument response functions - csviscube --- Computes visibility cube cslightcrv --- Computes lightcurve csmodelinfo --- Shows model container content csmodelmerge --- Merges several model containers into one file @@ -52,11 +51,12 @@ cscripts csresmap --- Generates residual map cssens --- Computes CTA sensitivity csspec --- Computes spectral points + cssrcdetect --- Detects sources in sky map cstsdist --- Generates Test Statistic distribution cstsmapsplit --- Creates commands to split the Test Statistic map computations cstsmapmerge --- Merges slices from Test Statistic map computations - - + csviscube --- Computes visibility cube + Scripts to manage an IACT database ---------------------------------- diff --git a/test/data/skymap.fits b/test/data/skymap.fits new file mode 100644 index 0000000..7782719 --- /dev/null +++ b/test/data/skymap.fits diff --git a/test/test_cssrcdetect.py b/test/test_cssrcdetect.py new file mode 100755 index 0000000..02e1461 --- /dev/null +++ b/test/test_cssrcdetect.py @@ -0,0 +1,175 @@ +#! /usr/bin/env python +# ========================================================================== +# This scripts performs unit tests for the cssrcdetect script +# +# Copyright (C) 2016 Juergen Knoedlseder +# +# This program 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. +# +# This program 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 this program. If not, see . +# +# ========================================================================== +import os +import gammalib +import cscripts +from testing import test + + +# ================================= # +# Test class for cssrcdetect script # +# ================================= # +class Test(test): + """ + Test class for cssrcdetect script + + This test class makes unit tests for the cssrcdetect script by using it + from the command line and from Python. + """ + + # Constructor + def __init__(self): + """ + Constructor. + """ + # Call base class constructor + test.__init__(self) + + # Set members + self._map = self._datadir + '/skymap.fits' + + # Return + return + + # Set test functions + def set(self): + """ + Set all test functions. + """ + # Set test name + self.name('cssrcdetect') + + # Append tests + self.append(self._test_cmd, 'Test cssrcdetect on command line') + self.append(self._test_python, 'Test cssrcdetect from Python') + + # Return + return + + # Test cssrcdetect on command line + def _test_cmd(self): + """ + Test cssrcdetect on the command line + """ + # Set script name + cssrcdetect = self._script('cssrcdetect') + + # Setup cssrcdetect command + cmd = cssrcdetect+' inmap="'+self._map+'"'+ \ + ' outmodel="cssrcdetect_cmd1.xml"'+ \ + ' outds9file="cssrcdetect_cmd1.reg"'+ \ + ' srcmodel="POINT" bkgmodel="NONE" threshold=10.0'+ \ + ' logfile="cssrcdetect_cmd1.log" chatter=1' + + # Check if execution of wrong command fails + self.test_assert(self._execute('command_that_does_not_exist') != 0, + 'Self test of test script') + + # Check if execution was successful + self.test_assert(self._execute(cmd) == 0, + 'Check successful execution from command line') + + # Check model file + self._check_model_file('cssrcdetect_cmd1.xml', 2) + + # Setup cssrcdetect command + cmd = cssrcdetect+' inmap="skymap_that_does_not_exist"'+ \ + ' outmodel="cssrcdetect_cmd1.xml"'+ \ + ' outds9file="cssrcdetect_cmd1.reg"'+ \ + ' srcmodel="POINT" bkgmodel="NONE" threshold=10.0'+ \ + ' logfile="cssrcdetect_cmd2.log" chatter=2' + + # Check if execution failed + self.test_assert(self._execute(cmd) != 0, + 'Check invalid input file when executed from command line') + + # Return + return + + # Test cssrcdetect from Python + def _test_python(self): + """ + Test cssrcdetect from Python + """ + # Allocate empty cssrcdetect script + srcdetect = cscripts.cssrcdetect() + + # Check that saving does not nothing + srcdetect['outmodel'] = 'cssrcdetect_py0.xml' + srcdetect['outds9file'] = 'cssrcdetect_py0.reg' + srcdetect['logfile'] = 'cssrcdetect_py0.log' + srcdetect.logFileOpen() + srcdetect.save() + + # Check than an empty model definition file has been created + self._check_model_file('cssrcdetect_py0.xml', 0) + + # Check that clearing does not lead to an exception or segfault + srcdetect.clear() + + # Now set cssrcdetect parameters + srcdetect['inmap'] = self._map + srcdetect['outmodel'] = 'cssrcdetect_py1.xml' + srcdetect['outds9file'] = 'cssrcdetect_py1.reg' + srcdetect['srcmodel'] = 'POINT' + srcdetect['bkgmodel'] = 'NONE' + srcdetect['threshold'] = 10.0 + srcdetect['logfile'] = 'cssrcdetect_py1.log' + srcdetect['chatter'] = 2 + + # Run cssrcdetect script and save models + srcdetect.logFileOpen() # Make sure we get a log file + srcdetect.run() + srcdetect.save() + + # Check model file + self._check_model_file('cssrcdetect_py1.xml', 2) + + # Execute cssrcdetect script again, now with a higher chatter level + # than before + srcdetect['outmodel'] = 'cssrcdetect_py2.xml' + srcdetect['outds9file'] = 'cssrcdetect_py2.reg' + srcdetect['bkgmodel'] = 'IRF' + srcdetect['logfile'] = 'cssrcdetect_py2.log' + srcdetect['chatter'] = 3 + srcdetect.logFileOpen() # Needed to get a new log file + srcdetect.execute() + + # Check model file + self._check_model_file('cssrcdetect_py2.xml', 3) + + # Return + return + + # Check model file + def _check_model_file(self, filename, number): + """ + Check model file + """ + # Open models + models = gammalib.GModels(filename) + + # Check number of models + self.test_value(models.size(), number, 'Check for '+str(number)+ + ' models in XML file') + + # Return + return diff --git a/test/test_python_cscripts.py b/test/test_python_cscripts.py index 1dfebed..177dd47 100755 --- a/test/test_python_cscripts.py +++ b/test/test_python_cscripts.py @@ -40,6 +40,7 @@ import test_csresmap import test_csroot2caldb import test_cssens import test_csspec +import test_cssrcdetect import test_cstsdist import test_cstsmapmerge import test_cstsmapsplit @@ -126,6 +127,7 @@ def test(installed=False, debug=False): test_csroot2caldb.Test(), test_cssens.Test(), test_csspec.Test(), + test_cssrcdetect.Test(), test_cstsdist.Test(), test_cstsmapmerge.Test(), test_cstsmapsplit.Test(),