Package pylal :: Package dq :: Module dqPlotUtils
[hide private]
[frames] | no frames]

Source Code for Module pylal.dq.dqPlotUtils

   1  #!/usr/bin/env python 
   2   
   3  # Copyright (C) 2011 Duncan Macleod 
   4  # 
   5  # This program is free software; you can redistribute it and/or modify it 
   6  # under the terms of the GNU General Public License as published by the 
   7  # Free Software Foundation; either version 3 of the License, or (at your 
   8  # option) any later version. 
   9  # 
  10  # This program is distributed in the hope that it will be useful, but 
  11  # WITHOUT ANY WARRANTY; without even the implied warranty of 
  12  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General 
  13  # Public License for more details. 
  14  # 
  15  # You should have received a copy of the GNU General Public License along 
  16  # with this program; if not, write to the Free Software Foundation, Inc., 
  17  # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA. 
  18   
  19  # ============================================================================= 
  20  # Preamble 
  21  # ============================================================================= 
  22   
  23  from __future__ import division 
  24  import math,re,numpy,itertools,copy,matplotlib,sys,warnings 
  25   
  26  # test matplotlib backend and reset if needed 
  27  from os import getenv 
  28  _display = getenv('DISPLAY','') 
  29  _backend_warn = """No display detected, moving to 'Agg' backend in matplotlib.""" 
  30  if not _display and not matplotlib.get_backend().lower() == 'agg': 
  31    warnings.warn(_backend_warn) 
  32    matplotlib.use('Agg', warn=False) 
  33  import pylab 
  34   
  35  try: from mpl_toolkits.basemap import Basemap 
  36  except ImportError,e: warnings.warn(str(e)) 
  37  try: from mpl_toolkits import axes_grid 
  38  except ImportError,e: warnings.warn(str(e)) 
  39   
  40  from datetime import datetime 
  41  from glue import segments,git_version 
  42  from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS 
  43  from lalburst import date 
  44  from pylal import plotutils 
  45  from pylal.dq.dqTriggerUtils import def_get_time,get_column 
  46   
  47  __author__  = "Duncan Macleod <duncan.macleod@astro.cf.ac.uk>" 
  48  __version__ = "git id %s" % git_version.id 
  49  __date__    = git_version.date 
  50   
  51  """ 
  52  This module provides plotting routines for use in data quality investigations. All routines are written to work in as general a way as possible with ligolw tables and lsctables compatible columns, and to plot in a similar pythonic way to pylal.plotutils.  
  53  """ 
54 55 # ============================================================================= 56 # Set plot parameters aux helper functions 57 # ============================================================================= 58 59 -def set_rcParams():
60 61 """ 62 Customise the figure parameters. 63 """ 64 65 # customise plot appearance 66 pylab.rcParams.update({"text.usetex": True, 67 "text.verticalalignment": "center", 68 "lines.linewidth": 2, 69 "xtick.labelsize": 16, 70 "ytick.labelsize": 16, 71 "axes.titlesize": 22, 72 "axes.labelsize": 16, 73 "axes.linewidth": 1, 74 "grid.linewidth": 1, 75 "legend.fontsize": 16, 76 "legend.loc": "best", 77 "figure.figsize": [12,6], 78 "figure.dpi": 80, 79 "image.origin": 'lower', 80 "axes.grid": True, 81 "axes.axisbelow": False })
82
83 -def set_ticks(ax, calendar_time=False):
84 85 """ 86 Format the x- and y-axis ticks to ensure minor ticks appear when needed 87 and the x-axis is set for spaces of 4 rather than 5. 88 89 Arguments: 90 91 ax : matplotlib.axes.AxesSubplot 92 Axes object to format 93 """ 94 95 # 96 # make sure we get minor ticks if there are no major ticks in the range 97 # 98 99 # xticks 100 ticks = list(ax.get_xticks()) 101 xlim = ax.get_xlim() 102 for i,tick in enumerate(ticks[::-1]): 103 if not xlim[0] <= tick <= xlim[1]: ticks.pop(-1) 104 if len(ticks)<=1: 105 ax.xaxis.set_minor_formatter(pylab.matplotlib.ticker.ScalarFormatter()) 106 107 # yticks 108 ticks = list(ax.get_yticks()) 109 ylim = ax.get_ylim() 110 for i,tick in enumerate(ticks[::-1]): 111 if not ylim[0] <= tick <= ylim[1]: ticks.pop(-1) 112 if len(ticks)<=1: 113 ax.yaxis.set_minor_formatter(pylab.matplotlib.ticker.ScalarFormatter()) 114 115 # set xticks in time format, python2.5 is not new enough for 116 # flexibility, recoding part of AutoDateFormatter to get it 117 if calendar_time: 118 dateLocator = pylab.matplotlib.dates.AutoDateLocator() 119 dateLocator.set_axis(ax.xaxis) 120 dateLocator.refresh() 121 scale = float( dateLocator._get_unit() ) 122 if ( scale == 365.0 ): 123 dateFormatter = pylab.matplotlib.dates.DateFormatter("$%Y$") 124 elif ( scale == 30.0 ): 125 dateFormatter = pylab.matplotlib.dates.DateFormatter("$%y$/$%m$") 126 elif ( (scale == 1.0) or (scale == 7.0) ): 127 dateFormatter = pylab.matplotlib.dates.DateFormatter("$%m$/$%d$") 128 elif ( scale == (1.0/24.0) ): 129 dateFormatter = pylab.matplotlib.dates.DateFormatter("$%d$-$%H$") 130 elif ( scale == (1.0/(24*60)) ): 131 dateFormatter = pylab.matplotlib.dates.DateFormatter("$%H$:$%M$") 132 elif ( scale == (1.0/(24*3600)) ): 133 dateFormatter = pylab.matplotlib.dates.DateFormatter("$%H$:$%M$") 134 135 ax.xaxis.set_major_locator(dateLocator) 136 ax.xaxis.set_major_formatter(dateFormatter) 137 138 else: 139 # set xticks for 4 hours rather than 5 140 xticks = ax.get_xticks() 141 if len(xticks)>1 and xticks[1]-xticks[0]==5: 142 ax.xaxis.set_major_locator(\ 143 pylab.matplotlib.ticker.MultipleLocator(base=2)) 144 145 146 # set tick linewidth 147 for line in ax.get_xticklines() + ax.get_yticklines(): 148 line.set_markersize(10) 149 line.set_markeredgewidth(1)
150
151 -def display_name(columnName):
152 153 """ 154 Format the string columnName (e.g. xml table column) into latex format for 155 an axis label. 156 157 Examples: 158 159 >>> display_name('snr') 160 'SNR' 161 162 >>> display_name('bank_chisq_dof') 163 'Bank $\\chi^2$ DOF' 164 165 Arguments: 166 167 columnName : string 168 string to format 169 """ 170 171 acro = ['snr', 'ra','dof', 'id', 'ms', 'far'] 172 greek = ['alpha', 'beta', 'gamma', 'delta', 'epsilon', 'zeta', 'eta',\ 173 'theta', 'iota', 'kappa', 'lamda', 'mu', 'nu', 'xi', 'omicron',\ 174 'pi', 'rho', 'sigma', 'tau', 'upsilon', 'phi', 'chi', 'psi', 'omega'] 175 unit = ['ns'] 176 sub = ['flow', 'fhigh', 'hrss', 'mtotal', 'mchirp'] 177 178 words = [] 179 for w in re.split('\s', columnName): 180 if w.isupper(): words.append(w) 181 else: words.extend(re.split('_', w)) 182 183 # parse words 184 for i,w in enumerate(words): 185 # get acronym in lower case 186 if w in acro: 187 words[i] = w.upper() 188 # get numerical unit 189 elif w in unit: 190 words[i] = '$(%s)$' % w 191 # get character with subscript text 192 elif w in sub: 193 words[i] = '%s$_{\mbox{\\small %s}}$' % (w[0], w[1:]) 194 # get greek word 195 elif w in greek: 196 words[i] = '$\%s$' % w 197 # get starting with greek word 198 elif re.match('(%s)' % '|'.join(greek), w): 199 if w[-1].isdigit(): 200 words[i] = '$\%s_{%s}$''' % tuple(re.findall(r"[a-zA-Z]+|\d+",w)) 201 elif w.endswith('sq'): 202 words[i] = '$\%s^2$' % w.rstrip('sq') 203 # get everything else 204 else: 205 if w.isupper(): 206 words[i] = w 207 else: 208 words[i] = w.title() 209 # escape underscore 210 words[i] = re.sub('(?<!\\\\)_', '\_', words[i]) 211 212 return ' '.join(words)
213
214 -def gps2datenum(gpstime):
215 216 """ 217 Convert GPS time into floating in standard python time format 218 (days since Jan 01 0000), don't seem to use properly leap seconds 219 """ 220 221 # set time of 0 GPS in datenum units 222 zeroGPS = 722820.0 223 ## correct for leap seconds assuming that all time stamps are within 224 # a range not including a leap 225 # select the first time stamp 226 if isinstance(gpstime,float) or isinstance(gpstime,int): 227 repTime = gpstime 228 else: 229 if len(gpstime)>0: 230 repTime = gpstime[0] 231 else: 232 return gpstime 233 234 zeroGPS = zeroGPS + float(date.XLALLeapSeconds(LIGOTimeGPS(0)) -\ 235 date.XLALLeapSeconds(LIGOTimeGPS(repTime)))/86400 236 # convert to datenum (days) 237 datenum = gpstime/86400 + zeroGPS 238 239 return datenum
240
241 -def calendar_time_unit(duration):
242 243 if duration > 63072000: 244 return "year" 245 elif duration > 5184000: 246 return "year/month" 247 elif duration > 604800: 248 return "month/day" 249 elif duration > 7200: 250 return "day-hour" 251 elif duration > 600: 252 return "hour:minute" 253 else: 254 return "hour:minute:second"
255
256 -def time_unit(duration):
257 258 """ 259 Work out renormalisation for the time axis, makes the label more 260 appropriate. Returns unit (in seconds) and string descriptor 261 262 Example: 263 264 >>> time_unit(100) 265 (1, 'seconds') 266 267 >>> time_unit(604800) 268 (86400, 'days') 269 270 Arguments: 271 272 duration : float 273 plot duration to normalise 274 """ 275 276 # set plot time unit whether it's used or not 277 if (duration) < 1000: 278 unit = 1 279 elif (duration) < 20000: 280 unit = 60 281 elif (duration) >= 20000 and (duration) < 604800: 282 unit = 3600 283 else: 284 unit = 86400 285 timestring = {1:'seconds', 60:'minutes',3600:'hours',86400:'days'} 286 287 return unit, timestring[unit]
288
289 -def getTrigAttribute(trig, col):
290 291 """ 292 Returns the value of trig.col or trig.get_col() for the given string col, 293 and the object trig. If col='time' is given, trig.get_peak() is returned for 294 *Burst* objects, trig.get_end() for *Inspiral* objects and trig.get_start() 295 for *Ringdown* objects. Raises KeyError if cannot execute. 296 297 Arguments: 298 299 trig : [ lsctables.SnglBurst | lscatbles.SnglInspiral | 300 lsctables.SnglRingdown ] 301 xml table entry from which to extract parameter 302 col : string 303 glue.ligolw.table column name to extract 304 """ 305 306 # return time 307 if col=='time': 308 if re.search('Burst', str(trig)): 309 return trig.get_peak() 310 elif re.search('Inspiral', str(trig)): 311 return trig.get_end() 312 elif re.search('Ringdown', str(trig)): 313 return trig.get_start() 314 else: 315 return trig.time 316 317 # return simple column 318 if col in trig.__slots__: 319 return getattr(trig, col) 320 321 # return get_XXX() parameter 322 try: 323 return eval('trig.get_%s()', col) 324 325 # if we get here, panic 326 except: 327 raise KeyError, "Column '%s' not found in %s." % (col, type(trig))
328
329 -def add_colorbar(self, mappable, cax=None, clim=None, log=False, base=10,\ 330 label=None, **kwargs):
331 332 # set axes for colorbar 333 #if not cax: 334 # div = axes_grid.make_axes_locatable(self.ax) 335 # cax = div.new_horizontal(size="2%", pad=0.05, visible=True) 336 #cax = div.append_axes("right", size="2%", pad=0.05) 337 338 cmin, cmax = clim 339 340 # construct colorbar tick formatter, using logic not supported before python 341 # 2.5 342 colorticksize = None 343 if log and numpy.power(base,cmax)-numpy.power(base,cmin) > 4: 344 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos:\ 345 numpy.power(base,x)>=1 and\ 346 "$%d$" % round(numpy.power(base, x)) or\ 347 numpy.power(base,x)>=0.01 and\ 348 "$%.2f$" % numpy.power(base, x) or\ 349 numpy.power(base,x)<0.01 and "$%f$") 350 elif log and numpy.power(base,cmax)-numpy.power(base,cmin) > 0.4: 351 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos: "$%.2f$"\ 352 % numpy.power(base, x)) 353 elif log: 354 colorticksize = 'small' 355 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos: '$%.2e$'\ 356 % numpy.power(base, x)) 357 elif not log and cmax-cmin > 4: 358 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos:\ 359 x>=1 and "$%d$" % round(x) or\ 360 x>=0.01 and "$%.2f$" % x or\ 361 x==0.0 and "$0$" or\ 362 x and "$%f$") 363 else: 364 formatter = None 365 366 if clim: 367 colorticks = numpy.linspace(cmin, cmax, 4) 368 else: 369 colorticks = None 370 371 self.colorbar = self.ax.figure.colorbar(mappable, cax=cax, format=formatter,\ 372 ticks=colorticks, pad=0.005,\ 373 fraction=0.02, aspect=40, **kwargs) 374 if colorticksize is not None: 375 for l in self.colorbar.ax.yaxis.get_ticklabels(): 376 l.set_size(colorticksize) 377 # colorbar = pylab.colorbar(mappable, cax=cax, format=formatter,\ 378 #ticks=colorticks, **kwargs) 379 self.colorbar.set_label(label) 380 self.colorbar.draw_all()
381
382 #try: 383 # self.colorbar = colorbar 384 #except: 385 # return colorbar 386 387 -def add_hidden_colorbar(self):
388 389 # set axes for colorbar 390 div = axes_grid.make_axes_locatable(self.ax) 391 cax = div.new_horizontal(size="2%", pad=0.05, visible=False)
392
393 # ============================================================================= 394 # Abstract classes for plots 395 # ============================================================================= 396 397 # ============================================================================= 398 # Class for segment plot 399 400 -class PlotSegmentsPlot(plotutils.BasicPlot):
401 402 """ 403 Horizontal bar segment plot. Based originally on PlotSegmentsPlot class in 404 pylal/bin/plotsegments.py 405 """ 406 407 color_code = {'H1':'r', 'H2':'b', 'L1':'g', 'V1':'m', 'G1':'k'} 408
409 - def __init__(self, xlabel="", ylabel="", title="", subtitle="", t0=0, unit=1,\ 410 calendar_time=False):
411 """ 412 Create a fresh plot. Provide t0 to provide a reference time to use as 413 zero. 414 """ 415 plotutils.BasicPlot.__init__(self, xlabel, ylabel, title) 416 self.ax.set_title(title, x=0.5, y=1.025) 417 self.ax.text(0.5, 1.035, subtitle, horizontalalignment='center', 418 transform=self.ax.transAxes, verticalalignment='top') 419 self.segdict = segments.segmentlistdict() 420 self.keys = [] 421 if calendar_time: 422 self._time_transform =\ 423 lambda seg: segments.segment(gps2datenum(float(seg[0])),\ 424 gps2datenum(float(seg[1]))) 425 else: 426 self._time_transform =\ 427 lambda seg: segments.segment(float(seg[0]-t0)/unit,\ 428 float(seg[1]-t0)/unit)
429 430 431
432 - def add_content(self, segdict, keys=None, t0=0, unit=1):
433 if not keys: 434 keys = sorted(segdict.keys()) 435 for key in keys: 436 self.segdict[key] = segdict[key] 437 self.keys.extend(keys)
438
439 - def highlight_segment(self, seg, **plot_args):
440 """ 441 Highlight a particular segment with dashed lines. 442 """ 443 a,b = self._time_transform(seg) 444 plot_args.setdefault('linestyle', '--') 445 plot_args.setdefault('color','r') 446 self.ax.axvline(a, **plot_args) 447 self.ax.axvline(b, **plot_args)
448 449 @plotutils.method_callable_once
450 - def finalize(self, labels_inset=False, hidden_colorbar=False):
451 452 for row,key in enumerate(self.keys): 453 if self.color_code.has_key(key): 454 col = self.color_code[key] 455 else: 456 col = 'b' 457 for seg in self.segdict[key]: 458 a,b = self._time_transform(seg) 459 self.ax.fill([a, b, b, a, a],\ 460 [row-0.4, row-0.4, row+0.4, row+0.4, row-0.4], 'b') 461 if labels_inset: 462 self.ax.text(0.01,(row+1)/(len(self.keys)+1), re.sub('\\+_+','\_',key),\ 463 horizontalalignment='left', verticalalignment='center',\ 464 transform=self.ax.transAxes, backgroundcolor='white',\ 465 bbox=dict(facecolor='white', alpha=0.5, edgecolor='none')) 466 467 ticks = pylab.arange(len(self.keys)) 468 self.ax.set_yticks(ticks) 469 if labels_inset: 470 self.ax.set_yticklabels(ticks, color='white') 471 else: 472 self.ax.set_yticklabels([re.sub(r'\\+_+', '\_', k)\ 473 for k in self.keys], size='small') 474 self.ax.set_ylim(-1, len(self.keys)) 475 476 # add hidden colorbar 477 if hidden_colorbar: 478 self.add_hidden_colorbar()
479 480 PlotSegmentsPlot.add_hidden_colorbar = add_hidden_colorbar
481 482 # ============================================================================= 483 # Class for standard scatter plot 484 485 -class ScatterPlot(plotutils.BasicPlot):
486 487 """ 488 A simple scatter plot, taking x- and y-axis data. 489 """ 490
491 - def __init__(self, xlabel="", ylabel="", title="", subtitle=""):
492 plotutils.BasicPlot.__init__(self, xlabel, ylabel, title) 493 self.ax.set_title(title, x=0.5, y=1.025) 494 self.ax.text(0.5, 1.035, subtitle, horizontalalignment='center', 495 transform=self.ax.transAxes, verticalalignment='top') 496 self.x_data_sets = [] 497 self.y_data_sets = [] 498 self.kwarg_sets = []
499
500 - def add_content(self, x_data, y_data, **kwargs):
501 self.x_data_sets.append(x_data) 502 self.y_data_sets.append(y_data) 503 self.kwarg_sets.append(kwargs)
504 505 @plotutils.method_callable_once
506 - def finalize(self, loc='best', logx=False, logy=False, hidden_colorbar=False):
507 # make plot 508 for x_vals, y_vals, plot_kwargs, c in \ 509 itertools.izip(self.x_data_sets, self.y_data_sets, self.kwarg_sets,\ 510 plotutils.default_colors()): 511 plot_kwargs.setdefault("marker", "o") 512 plot_kwargs.setdefault("s", 20) 513 plot_kwargs.setdefault("linewidth", 1) 514 515 self.ax.scatter(x_vals, y_vals, c=c, **plot_kwargs) 516 517 # add hidden colorbar to make plots the right aspect 518 if hidden_colorbar: 519 self.add_hidden_colorbar() 520 521 # add legend if there are any non-trivial labels 522 self.add_legend_if_labels_exist(loc=loc) 523 524 leg = self.ax.get_legend() 525 # set transparent legend 526 if leg: 527 legfr = leg.get_frame() 528 legfr.set_alpha(0.5) 529 530 # set axes 531 if logx: 532 self.ax.set_xscale('log') 533 if logy: 534 self.ax.set_yscale('log')
535 536 ScatterPlot.add_hidden_colorbar = add_hidden_colorbar
537 538 # ============================================================================= 539 # Class for scatter plot with colorbar 540 541 -class ColorbarScatterPlot(plotutils.BasicPlot):
542 543 """ 544 A scatter plot of x- versus y-data, coloured by z-data. A single colorbar 545 is used, from the union of the ranges of z-data for each content set, 546 unless strict limits are passed as clim=(cmin, cmax) to finalize(). 547 """ 548
549 - def __init__(self, xlabel="", ylabel="", zlabel="", title="", subtitle=""):
550 plotutils.BasicPlot.__init__(self, xlabel, ylabel, title) 551 self.ax.set_title(title, x=0.5, y=1.025) 552 self.ax.text(0.5, 1.035, subtitle, horizontalalignment='center', 553 transform=self.ax.transAxes, verticalalignment='top') 554 555 self.x_data_sets = [] 556 self.y_data_sets = [] 557 self.z_data_sets = [] 558 self.rank_data_sets = [] 559 self.kwarg_sets = [] 560 self.color_label = zlabel
561
562 - def add_content(self, x_data, y_data, z_data, rank_data, **kwargs):
563 self.x_data_sets.append(x_data) 564 self.y_data_sets.append(y_data) 565 self.z_data_sets.append(z_data) 566 self.rank_data_sets.append(rank_data) 567 self.kwarg_sets.append(kwargs)
568 569 @plotutils.method_callable_once
570 - def finalize(self, loc='best', logx=False, logy=False, logz=False, clim=None,\ 571 base=10):
572 573 # set up things we'll need later 574 p = [] 575 576 # set colorbar limits 577 if clim: 578 cmin,cmax = clim 579 else: 580 try: 581 cmin = min([min(z_vals) for z_vals in self.z_data_sets])*0.99 582 cmax = max([max(z_vals) for z_vals in self.z_data_sets])*1.01 583 # catch no triggers 584 except ValueError: 585 cmin = 1 586 cmax = 10 587 588 # reset logs 589 if logz: 590 cmin = numpy.math.log(cmin, base) 591 cmax = numpy.math.log(cmax, base) 592 593 # make plot 594 for x_vals, y_vals, z_vals, rank_vals, plot_kwargs in\ 595 itertools.izip(self.x_data_sets, self.y_data_sets, self.z_data_sets,\ 596 self.rank_data_sets, self.kwarg_sets): 597 598 if logz: z_vals = [numpy.math.log(z, base) for z in z_vals] 599 600 plot_kwargs.setdefault("marker", "o") 601 plot_kwargs.setdefault("vmin", cmin) 602 plot_kwargs.setdefault("vmax", cmax) 603 604 # sort data by z-value 605 zipped = zip(x_vals, y_vals, z_vals, rank_vals) 606 zipped.sort(key=lambda (x,y,z,r): r) 607 x_vals, y_vals, z_vals, rank_vals = map(list, zip(*zipped)) 608 609 p.append(self.ax.scatter(x_vals, y_vals, c=z_vals, **plot_kwargs)) 610 611 if len(p)<1: 612 p.append(self.ax.scatter([1], [1], c=[cmin], vmin=cmin,\ 613 vmax=cmax, visible=False)) 614 615 616 # write colorbar 617 #self.set_colorbar(p[-1], [cmin, cmax], logz, base) 618 self.add_colorbar(p[-1], clim=[cmin, cmax], log=logz,\ 619 label=self.color_label, base=10) 620 621 # set axes 622 if logx: 623 self.ax.set_xscale('log') 624 if logy: 625 self.ax.set_yscale('log') 626 627 self.add_legend_if_labels_exist(loc=loc) 628 629 # set transparent legend 630 leg = self.ax.get_legend() 631 if leg: 632 legfr = leg.get_frame() 633 legfr.set_alpha(0.5)
634
635 - def set_colorbar(self, mappable, clim=None, log=True, base=10):
636 637 cmin, cmax = clim 638 639 # construct colorbar tick formatter, using logic not supported before python 640 # 2.5 641 if log and numpy.power(base,cmax)-numpy.power(base,cmin) > 4: 642 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos:\ 643 numpy.power(base,x)>=1 and\ 644 "$%d$" % round(numpy.power(base, x)) or\ 645 numpy.power(base,x)>=0.01 and\ 646 "$%.2f$" % numpy.power(base, x) or\ 647 numpy.power(base,x)<0.01 and "$%f$") 648 elif log and numpy.power(base,cmax)-numpy.power(base,cmin) > 0.4: 649 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos: "$%.2f$"\ 650 % numpy.power(base, x)) 651 elif log: 652 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos: "$%.4e$"\ 653 % numpy.power(base, x)) 654 else: 655 formatter = pylab.matplotlib.ticker.FuncFormatter(lambda x,pos: "$%.4e$"\ 656 % x) 657 658 if clim: 659 colorticks = numpy.linspace(cmin, cmax, 4) 660 self.colorbar = self.ax.figure.colorbar(mappable, format=formatter,\ 661 ticks=colorticks) 662 else: 663 self.colorbar = self.ax.figure.colorbar(mappable, format=formatter) 664 665 self.colorbar.set_label(self.color_label) 666 self.colorbar.draw_all()
667 668 ColorbarScatterPlot.add_colorbar = add_colorbar
669 670 # ============================================================================= 671 # Extension of ColorbarScatterPlot to plot DetChar-style scatter plot 672 673 -class DetCharScatterPlot(ColorbarScatterPlot):
674 675 """ 676 A 'DetChar' style scatter plot, whereby those triggers under a threshold 677 on the colour column are plotted much smaller than others, allowing line 678 features to be shown easily. 679 """ 680 681 @plotutils.method_callable_once
682 - def finalize(self, loc='best', logx=False, logy=False, logz=True, base=10,\ 683 clim=None, rankthreshold=None):
684 685 p = [] 686 # set colorbar limits 687 if clim: 688 cmin,cmax = clim 689 else: 690 try: 691 cmin = min([min(z_vals) for z_vals in self.z_data_sets])*0.99 692 cmax = max([max(z_vals) for z_vals in self.z_data_sets])*1.01 693 # catch no triggers 694 except ValueError: 695 cmin = 1 696 cmax = 10 697 698 # reset logs 699 if logz: 700 cmin = numpy.math.log(cmin, base) 701 cmax = numpy.math.log(cmax, base) 702 703 # make plot 704 for x_vals, y_vals, z_vals, rank_vals, plot_kwargs in\ 705 itertools.izip(self.x_data_sets, self.y_data_sets, self.z_data_sets,\ 706 self.rank_data_sets, self.kwarg_sets): 707 plot_kwargs.setdefault("vmin", cmin) 708 plot_kwargs.setdefault("vmax", cmax) 709 plot_kwargs.setdefault("s", 15) 710 plot_kwargs.setdefault("marker", "o") 711 plot_kwargs.setdefault("edgecolor", "k") 712 713 if logz: 714 z_vals = [numpy.math.log(z, base) for z in z_vals] 715 716 # sort data by z-value 717 zipped = zip(x_vals, y_vals, z_vals, rank_vals) 718 zipped.sort(key=lambda (x,y,z,r): r) 719 x_vals, y_vals, z_vals, rank_vals = map(list, zip(*zipped)) 720 721 bins = ['low', 'high'] 722 x_bins = {} 723 y_bins = {} 724 z_bins = {} 725 for bin in bins: 726 x_bins[str(bin)] = [] 727 y_bins[str(bin)] = [] 728 z_bins[str(bin)] = [] 729 for i in xrange(len(x_vals)): 730 if rankthreshold and rank_vals[i] < rankthreshold: 731 x_bins[bins[0]].append(float(x_vals[i])) 732 y_bins[bins[0]].append(float(y_vals[i])) 733 z_bins[bins[0]].append(float(z_vals[i])) 734 else: 735 x_bins[bins[1]].append(float(x_vals[i])) 736 y_bins[bins[1]].append(float(y_vals[i])) 737 z_bins[bins[1]].append(float(z_vals[i])) 738 739 # plot bins 740 for i,bin in enumerate(bins): 741 if bin == bins[0]: 742 args = copy.deepcopy(plot_kwargs) 743 args['s']/=4 744 if not (args.has_key('marker') and args['marker']=='x'): 745 args['edgecolor'] = 'none' 746 if len(x_bins[bin])>=1: 747 p.append(self.ax.scatter(x_bins[bin], y_bins[bin], c=z_bins[bin],\ 748 **args)) 749 else: 750 if len(x_bins[bin])>=1: 751 p.append(self.ax.scatter(x_bins[bin], y_bins[bin], c=z_bins[bin],\ 752 **plot_kwargs)) 753 args = plot_kwargs 754 755 if len(p)<1: 756 p.append(self.ax.scatter([1], [1], c=[cmin], vmin=cmin,\ 757 vmax=cmax, visible=False)) 758 759 # write colorbar 760 self.add_colorbar(p[-1], clim=[cmin, cmax], log=logz,\ 761 label=self.color_label, base=10) 762 #self.set_colorbar(p[-1], [cmin, cmax], logz, base) 763 764 # set axes 765 if logx: 766 self.ax.set_xscale('log') 767 if logy: 768 self.ax.set_yscale('log') 769 770 self.add_legend_if_labels_exist(loc=loc)
771
772 # ============================================================================= 773 # Class for line histogram 774 775 -class LineHistogram(ColorbarScatterPlot, plotutils.BasicPlot):
776 777 """ 778 A simple line histogram plot. The values of each histogram bin are plotted 779 using pylab.plot(), with points centred on the x values and height equal 780 to the y values. 781 782 Cumulative, and rate options can be passeed to the finalize() method to 783 format each trace individually. 784 785 """ 786
787 - def __init__(self, xlabel="", ylabel="", title="", subtitle="",\ 788 colorlabel=""):
789 plotutils.BasicPlot.__init__(self, xlabel, ylabel, title) 790 self.ax.set_title(title, x=0.5, y=1.025) 791 self.ax.text(0.5, 1.035, subtitle, horizontalalignment='center', 792 transform=self.ax.transAxes, verticalalignment='top') 793 self.data_sets = [] 794 self.livetimes = [] 795 self.kwarg_sets = [] 796 self.color_label = colorlabel 797 self.color_values = []
798
799 - def add_content(self, data, livetime=1, cval=None, **kwargs):
800 self.data_sets.append(data) 801 self.livetimes.append(livetime) 802 self.kwarg_sets.append(kwargs) 803 self.color_values.append(cval)
804 805 @plotutils.method_callable_once
806 - def finalize(self, loc='best', num_bins=100, cumulative=False, rate=False,\ 807 logx=False, logy=False, fill=False, base=10, xlim=None,\ 808 colorbar=None, clim=None, hidden_colorbar=False):
809 810 # determine binning 811 if xlim: 812 min_stat, max_stat = xlim 813 else: 814 min_stat, max_stat = plotutils.determine_common_bin_limits(self.data_sets) 815 816 if min_stat!=max_stat: 817 max_stat *= 1.001 818 if logx: 819 bins = numpy.logspace(numpy.math.log(min_stat, base),\ 820 numpy.math.log(max_stat, base),\ 821 num_bins+1, endpoint=True) 822 else: 823 bins = numpy.linspace(min_stat, max_stat, num_bins + 1, endpoint=True) 824 else: 825 bins = [min_stat] 826 827 if logy: 828 ymin = 1E-100 829 else: 830 ymin = 0 831 832 # get colors 833 if colorbar: 834 if isinstance(colorbar, str): 835 ctype = colorbar 836 cmap = pylab.matplotlib.colors.LinearSegmentedColormap('clrs',\ 837 pylab.matplotlib.cm.jet._segmentdata) 838 else: 839 ctype, cmap = colorbar 840 assert re.match('(lin|log)', ctype, re.I),\ 841 "colorbar must have type 'linear', or 'log'" 842 try: 843 colors = pylab.matplotlib.colors.makeMappingArray(100000,\ 844 cmap) 845 except IndexError: 846 xind = numpy.linspace(0, 1, 100000) 847 colors = numpy.clip(numpy.array(cmap(xind), dtype=numpy.float), 0, 1) 848 849 if clim: 850 cmin,cmax = clim 851 else: 852 try: 853 cmin = min(self.color_values) 854 cmax = max(self.color_values) 855 except ValueError: 856 cmin = 1 857 cmax = 10 858 859 if re.search('log', ctype, re.I): 860 cmin = numpy.math.log(cmin, base) 861 cmax = numpy.math.log(cmax, base) 862 self.color_values = [numpy.math.log(x, base) for x in self.color_values] 863 color_idx = lambda x: int((x-cmin)/(cmax-cmin)*(len(colors)-1)) 864 self.color_values = [colors[color_idx(x)] for x in self.color_values] 865 866 # add hidden colorbar to make plots the right aspect 867 elif hidden_colorbar: 868 self.add_hidden_colorbar() 869 870 p = [] 871 for i, (data_set, livetime, plot_kwargs, col) in\ 872 enumerate(itertools.izip(self.data_sets, self.livetimes,\ 873 self.kwarg_sets, self.color_values)): 874 875 # 876 # make histogram 877 # 878 879 # get version 880 v = [int(i) for i in numpy.version.version.split('.')] 881 if v[1] < 1: 882 y, x = numpy.histogram(data_set, bins=bins) 883 elif v[1] < 3: 884 y, x = numpy.histogram(data_set, bins=bins, new=True) 885 x = x[:-1] 886 else: 887 y, x = numpy.histogram(data_set, bins=bins) 888 x = x[:-1] 889 890 # get cumulative sum 891 if cumulative: 892 y = y[::-1].cumsum()[::-1] 893 894 # convert to rate 895 if rate: 896 if livetime > 0: 897 y = y/livetime 898 ymin /= livetime 899 else: 900 y = numpy.empty(y.shape) 901 y.fill(numpy.nan) 902 ymin = numpy.nan 903 904 # reset zeros on logscale, tried with numpy, unreliable 905 if logy: 906 y = y.astype(float) 907 numpy.putmask(y, y==0, ymin) 908 909 # set color 910 if colorbar: 911 plot_kwargs['color'] = col 912 913 # plot 914 if fill: 915 plot_kwargs.setdefault('linewidth', 1) 916 plot_kwargs.setdefault('alpha', 0.8) 917 p.append(self.ax.plot(x, y, **plot_kwargs)) 918 self.ax.fill_between(x, 1E-100, y, **plot_kwargs) 919 else: 920 p.append(self.ax.plot(x, y, **plot_kwargs)) 921 922 # set axes 923 if logx: 924 self.ax.set_xscale('log') 925 if logy: 926 self.ax.set_yscale('log') 927 928 # set colorbar 929 if colorbar: 930 self.add_colorbar(self.ax.scatter([1], [1], c=1, cmap=cmap,\ 931 vmin=cmin, vmax=cmax, visible=False),\ 932 clim=[cmin, cmax], label=self.color_label) 933 elif hidden_colorbar: 934 self.add_hidden_colorbar() 935 936 # add legend if there are any non-trivial labels 937 self.add_legend_if_labels_exist(loc=loc) 938 939 # fix legend 940 leg = self.ax.get_legend() 941 if leg: 942 for l in leg.get_lines(): 943 l.set_linewidth(4)
944 945 LineHistogram.add_hidden_colorbar = add_hidden_colorbar
946 947 # ============================================================================= 948 # Extension of VerticalBarHistogram to include log axes 949 950 -class VerticalBarHistogram(plotutils.VerticalBarHistogram):
951
952 - def __init__(self, xlabel="", ylabel="", title="", subtitle=""):
953 plotutils.BasicPlot.__init__(self, xlabel, ylabel, title) 954 self.ax.set_title(title, x=0.5, y=1.025) 955 self.ax.text(0.5, 1.035, subtitle, horizontalalignment='center', 956 transform=self.ax.transAxes, verticalalignment='top') 957 self.data_sets = [] 958 self.kwarg_sets = []
959 960 @plotutils.method_callable_once
961 - def finalize(self, num_bins=20, normed=False, logx=False, logy=False,\ 962 base=10, hidden_colorbar=False, loc='best'):
963 964 # determine binning 965 min_stat, max_stat = plotutils.determine_common_bin_limits(self.data_sets) 966 if logx and min_stat!=0 and max_stat!=0: 967 bins = numpy.logspace(numpy.math.log(min_stat, base),\ 968 numpy.math.log(max_stat, base),\ 969 num_bins+1, endpoint=True) 970 else: 971 bins = numpy.linspace(min_stat, max_stat, num_bins + 1, endpoint=True) 972 973 # determine bar width; gets silly for more than a few data sets 974 if logx: 975 width = [bins[i+1]-bins[i] for i in xrange(len(bins)-1)] 976 width.append(width[-1]) 977 else: 978 width = (1 - 0.1 * len(self.data_sets)) * (bins[1] - bins[0]) 979 980 width = numpy.asarray(width)/2 981 982 # get version 983 v = map(int, numpy.version.version.split('.')) 984 if v[1] >= 1: width = width[:-1] 985 986 # set base of plot in log scale 987 if logy: 988 ymin = (base**-1)*5 989 else: 990 ymin = 0 991 992 # make plot 993 legends = [] 994 plot_list = [] 995 for i, (data_set, plot_kwargs) in \ 996 enumerate(itertools.izip(self.data_sets, self.kwarg_sets)): 997 # set default values 998 plot_kwargs.setdefault("alpha", 0.6) 999 plot_kwargs.setdefault("align", "center") 1000 plot_kwargs.setdefault("width", width) 1001 if logy: 1002 plot_kwargs.setdefault("bottom", ymin) 1003 else: 1004 plot_kwargs.setdefault("bottom", ymin) 1005 1006 # make histogram 1007 if v[1] < 1: 1008 y, x = numpy.histogram(data_set, bins=bins, normed=normed) 1009 elif v[1] < 3: 1010 y, x = numpy.histogram(data_set, bins=bins, new=True, normed=normed) 1011 x = x[:-1] 1012 else: 1013 y, x = numpy.histogram(data_set, bins=bins, normed=normed) 1014 x = x[:-1] 1015 1016 if logy: 1017 y = y-ymin 1018 1019 # stagger bins for pure aesthetics 1020 #x += 0.1 * i * max_stat / num_bins 1021 1022 # plot 1023 self.ax.bar(x, y, **plot_kwargs) 1024 1025 # set axes 1026 if logx: 1027 self.ax.set_xscale('log') 1028 if logy: 1029 self.ax.set_yscale('log') 1030 1031 # set logy minimum 1032 self.ax.set_ybound(lower=ymin) 1033 1034 # add legend if there are any non-trivial labels 1035 self.add_legend_if_labels_exist(loc=loc) 1036 # set transparent legend 1037 leg = self.ax.get_legend() 1038 if leg: leg.get_frame().set_alpha(0.5) 1039 1040 1041 # add hidden colorbar 1042 if hidden_colorbar: 1043 self.add_hidden_colorbar()
1044 1045 VerticalBarHistogram.add_hidden_colorbar = add_hidden_colorbar
1046 1047 # ============================================================================= 1048 # Class for time series plot 1049 1050 -class DataPlot(plotutils.BasicPlot):
1051 1052 """ 1053 Time-series data plot. Just a nice looking line plot. 1054 """ 1055
1056 - def __init__(self, xlabel="", ylabel="", title="", subtitle=""):
1057 plotutils.BasicPlot.__init__(self, xlabel, ylabel, title) 1058 self.ax.set_title(title, x=0.5, y=1.025) 1059 self.ax.text(0.5, 1.035, subtitle, horizontalalignment='center', 1060 transform=self.ax.transAxes, verticalalignment='top') 1061 self.x_data_sets = [] 1062 self.y_data_sets = [] 1063 self.kwarg_sets = []
1064
1065 - def add_content(self, x_data, y_data, **kwargs):
1066 self.x_data_sets.append(numpy.ma.asarray(x_data)) 1067 self.y_data_sets.append(numpy.ma.asarray(y_data)) 1068 self.kwarg_sets.append(kwargs)
1069 1070 @plotutils.method_callable_once
1071 - def finalize(self, loc='best', logx=False, logy=False, hidden_colorbar=False):
1072 # make plot 1073 plots = [] 1074 markersizes = [] 1075 markerscale = 4 1076 1077 for x_vals, y_vals, plot_kwargs, c in \ 1078 itertools.izip(self.x_data_sets, self.y_data_sets,\ 1079 self.kwarg_sets, plotutils.default_colors()): 1080 1081 # magnify the markers on the legend 1082 plot_kwargs.setdefault('markersize', 5) 1083 markersizes.append(plot_kwargs['markersize']) 1084 plot_kwargs['markersize'] = min(20, plot_kwargs['markersize']*markerscale) 1085 # plot 1086 if logx: 1087 x_vals = x_vals.astype(float) 1088 numpy.putmask(x_vals, x_vals==0, 1E-100) 1089 if logy: 1090 y_vals = y_vals.astype(float) 1091 numpy.putmask(y_vals, y_vals==0, 1E-100) 1092 1093 plots.append(self.ax.plot(x_vals, y_vals, **plot_kwargs)) 1094 1095 # add legend if there are any non-trivial labels 1096 self.add_legend_if_labels_exist(loc=loc) 1097 leg = self.ax.get_legend() 1098 # magnify the lines on the legend 1099 if leg: 1100 try: 1101 for l in leg.get_lines(): 1102 if l.get_linewidth(): 1103 l.set_linewidth(4) 1104 except AttributeError: 1105 pass 1106 1107 # reset markersizes on plot 1108 for i,plot in enumerate(plots): 1109 l = plot[0] 1110 l.set_markersize(markersizes[i]) 1111 1112 # set transparent legend 1113 if leg: 1114 legfr = leg.get_frame() 1115 legfr.set_alpha(0.5) 1116 1117 # add hidden colorbar 1118 if hidden_colorbar: self.add_hidden_colorbar() 1119 1120 # set axes 1121 if logx: self.ax.set_xscale('log') 1122 if logy: self.ax.set_yscale('log')
1123 1124 DataPlot.add_hidden_colorbar = add_hidden_colorbar
1125 1126 # ============================================================================= 1127 # Class for color map plot 1128 1129 -class ColorMap(ColorbarScatterPlot):
1130
1131 - def __init__(self, xlabel="", ylabel="", zlabel="", title="", subtitle=""):
1132 plotutils.BasicPlot.__init__(self, xlabel, ylabel, title) 1133 self.color_label = zlabel 1134 self.ax.set_title(title, x=0.5, y=1.025) 1135 self.ax.text(0.5, 1.035, subtitle, horizontalalignment='center', 1136 transform=self.ax.transAxes, verticalalignment='top') 1137 self.data_sets = [] 1138 self.extent_sets = [] 1139 self.kwarg_sets = []
1140
1141 - def add_content(self, data, extent, **kwargs):
1142 self.data_sets.append(numpy.asarray(data)) 1143 self.extent_sets.append(extent) 1144 self.kwarg_sets.append(kwargs)
1145 1146 @plotutils.method_callable_once
1147 - def finalize(self, loc='best', logx=False, logy=False, logz=False, clim=None,\ 1148 origin=pylab.rcParams['image.origin'], base=10):
1149 1150 # set colorbar limits 1151 if clim: 1152 cmin,cmax = clim 1153 else: 1154 cmin = min([z_vals.min() for z_vals in self.data_sets])*0.99 1155 cmax = max([z_vals.max() for z_vals in self.data_sets])*1.01 1156 1157 # reset logs 1158 if logz: 1159 cmin = numpy.math.log(cmin, base) 1160 cmax = numpy.math.log(cmax, base) 1161 1162 p = [] 1163 1164 for i, (data_set, extent, plot_kwargs) in \ 1165 enumerate(itertools.izip(self.data_sets, self.extent_sets,\ 1166 self.kwarg_sets)): 1167 1168 plot_kwargs.setdefault("vmin", cmin) 1169 plot_kwargs.setdefault("vmax", cmax) 1170 plot_kwargs.setdefault("norm", None) 1171 plot_kwargs.setdefault("interpolation", "kaiser") 1172 1173 if logz: 1174 if base==10: 1175 data_set = numpy.log10(data_set) 1176 elif base==numpy.e: 1177 data_set = numpy.log(data_set) 1178 elif base==2: 1179 data_set = numpy.log2(data_set) 1180 else: 1181 raise AttributeError("Can only use base = 2, e, or 10.") 1182 1183 p.append(self.ax.imshow(data_set, extent=extent, origin=origin,\ 1184 **plot_kwargs)) 1185 1186 if len(p)==0: 1187 p.append(self.ax.imshow([[1]], vmin=cmin,\ 1188 vmax=cmax, visible=False)) 1189 1190 # write colorbar 1191 self.add_colorbar(p[-1], clim=[cmin, cmax], log=logz,\ 1192 label=self.color_label, base=10) 1193 1194 # set axes 1195 if logx: self.ax.set_xscale('log') 1196 if logy: self.ax.set_yscale('log')
1197
1198 # ============================================================================= 1199 # Projection for sky position plotting 1200 1201 -class SkyPositionsPlot(ScatterPlot):
1202 1203 """ 1204 Plot of sky positions plotted onto mpl_toolkits basemap projection. 1205 """ 1206 1207 @plotutils.method_callable_once
1208 - def finalize(self, loc='lower right', projection='moll', centre=None,\ 1209 range=[(0, 0), (1, 1)]):
1210 1211 """ 1212 Finalize and draw plot. Can only be called once. 1213 1214 Keyword arguments: 1215 1216 toc : [ str | int ] 1217 compatible value for legend loc, default: 'lower right' 1218 projection : str 1219 valid projection name for mpl_toolkits.basemap.Basemap 1220 centre : list 1221 (ra, dec) in degrees for point to centre on plot, only works for 1222 some projections 1223 range : list 1224 [(xmin, ymin), (xmax, ymax)] pair of tuples for lower left and upper 1225 right corners of plot area. Full areai (entire sphere) is 1226 [(0,0), (1,1)]. Narrower range is zoomed in, wider range zoomed out. 1227 """ 1228 1229 # set up projection 1230 projection = projection.lower() 1231 if not centre: 1232 centre = (0,0) 1233 centre = [(-(centre[0]-180)+180) % 360, centre[1]] 1234 m1 = Basemap(projection=projection, lon_0=centre[0], lat_0=centre[1],\ 1235 resolution=None, ax=self.ax) 1236 1237 # set up range 1238 width = m1.urcrnrx 1239 height = m1.urcrnry 1240 range = [list(range[0]), list(range[1])] 1241 maprange = [[-width/2, -height/2], [width/2, height/2]] 1242 if range[0][0] != None: maprange[0][0] = width/2*(range[0][0]-1) 1243 if range[0][1] != None: maprange[0][1] = height/2*(range[0][1]-1) 1244 if range[1][0] != None: maprange[1][0] = width/2*(range[1][0]) 1245 if range[1][1] != None: maprange[1][1] = height/2*(range[1][1]) 1246 1247 # set projection 1248 m = Basemap(projection=projection, lon_0=centre[0], lat_0=centre[1],\ 1249 llcrnrx=maprange[0][0], llcrnry=maprange[0][1],\ 1250 urcrnrx=maprange[1][0], urcrnry=maprange[1][1],\ 1251 resolution=None, ax=self.ax) 1252 1253 # turn on 'fulldisk' property when using full disk 1254 if maprange == [[-width/2, -height/2], [width/2, height/2]]: 1255 m._fulldisk = True 1256 1257 xrange = (m.llcrnrx, m.urcrnrx) 1258 yrange = (m.llcrnry, m.urcrnry) 1259 1260 # make plot 1261 for x_vals, y_vals, plot_kwargs, c in\ 1262 itertools.izip(self.x_data_sets, self.y_data_sets, self.kwarg_sets,\ 1263 plotutils.default_colors()): 1264 1265 plot_kwargs.setdefault("marker", "o") 1266 plot_kwargs.setdefault("edgecolor", "k") 1267 plot_kwargs.setdefault("facecolor", c) 1268 1269 # project data 1270 convert = lambda x: (x>=180 and x-360) or (x<180 and x) 1271 x_vals = [-convert(x) for x in x_vals] 1272 if projection in ['moll', 'hammer', 'orth']: 1273 convert = lambda x: (x>=0 and x) or (x<0 and x+360) 1274 x_vals, y_vals = m(x_vals, y_vals) 1275 m.scatter(x_vals, y_vals, **plot_kwargs) 1276 1277 # finish projection 1278 m.drawmapboundary() 1279 1280 # set labels 1281 if projection in ['ortho']: 1282 plabels = [0, 0, 0, 0] 1283 mlabels = [0, 0, 0, 0] 1284 else: 1285 plabels = [1, 0, 0, 0] 1286 mlabels = [0, 0, 0, 0] 1287 1288 # draw parallels 1289 parallels = numpy.arange(-90., 120., 30.) 1290 m.drawparallels(parallels, labels=plabels,\ 1291 labelstyle='+/-', latmax=90) 1292 1293 # draw meridians 1294 if projection in ['moll', 'hammer', 'ortho']: 1295 meridians = numpy.arange(0., 360., 45.) 1296 else: 1297 meridians = numpy.arange(-180, 181, 45) 1298 m.drawmeridians(meridians, labels=mlabels,\ 1299 latmax=90, labelstyle='+/-') 1300 1301 # label parallels for certain projections 1302 if projection in ['ortho']: 1303 for lon,lat in zip([0.]*len(parallels), parallels): 1304 #if lonrange[0]<=lon<=lonrange[1] and latrange[0]<=lat<=latrange[1]: 1305 x, y = m(lon, lat) 1306 lon, lat = m1(x, y, inverse=True) 1307 if x<=10**20 and y<=10**20\ 1308 and xrange[0]<x<xrange[1] and yrange[0]<=y<=yrange[1]: 1309 m.ax.text(x, y, r"$%0.0f^\circ$" % lat) 1310 1311 # label meridians for all projections 1312 for lon,lat in zip(meridians, [0.]*len(meridians)): 1313 tlon = (-(lon-180)+180) % 360 1314 x, y = m(lon, lat) 1315 lon, lat = m1(x, y, inverse=True) 1316 1317 if x<=10**20 and y<=10**20\ 1318 and xrange[0]<x<xrange[1] and yrange[0]<=y<=yrange[1]: 1319 m.ax.text(x, y, r"$%0.0f^\circ$" % tlon) 1320 1321 # set legend 1322 self.add_legend_if_labels_exist(loc=loc, scatterpoints=1)
1323
1324 -class ColorbarSkyPositionsPlot(ColorbarScatterPlot):
1325 1326 @plotutils.method_callable_once
1327 - def finalize(self, loc='lower right', projection='moll', centre=None,\ 1328 range=[(0, 0), (1, 1)], logz=False, clim=None,\ 1329 base=10):
1330 1331 """ 1332 Finalize and draw plot. Can only be called once. 1333 1334 Keyword arguments: 1335 1336 toc : [ str | int ] 1337 compatible value for legend loc, default: 'lower right' 1338 projection : str 1339 valid projection name for mpl_toolkits.basemap.Basemap 1340 centre : list 1341 (ra, dec) in degrees for point to centre on plot, only works for 1342 some projections 1343 range : list 1344 [(xmin, ymin), (xmax, ymax)] pair of tuples for lower left and upper 1345 right corners of plot area. Full area (entire sphere) is 1346 [(0,0), (1,1)]. Narrower range is zoomed in, wider range zoomed out. 1347 logz : [ True | False ] 1348 plot z-axis in log scale 1349 clim : 1350 """ 1351 1352 p = [] 1353 1354 # set up projection 1355 projection = projection.lower() 1356 if not centre: 1357 centre = (0,0) 1358 m1 = Basemap(projection=projection, lon_0=centre[0], lat_0=centre[1],\ 1359 resolution=None, ax=self.ax) 1360 1361 # set up range 1362 width = m1.urcrnrx 1363 height = m1.urcrnry 1364 range = [list(range[0]), list(range[1])] 1365 maprange = [[-width/2, -height/2], [width/2, height/2]] 1366 if range[0][0] != None: maprange[0][0] = width/2*(range[0][0]-1) 1367 if range[0][1] != None: maprange[0][1] = height/2*(range[0][1]-1) 1368 if range[1][0] != None: maprange[1][0] = width/2*(range[1][0]) 1369 if range[1][1] != None: maprange[1][1] = height/2*(range[1][1]) 1370 1371 # set projection 1372 m = Basemap(projection=projection, lon_0=centre[0], lat_0=centre[1],\ 1373 llcrnrx=maprange[0][0], llcrnry=maprange[0][1],\ 1374 urcrnrx=maprange[1][0], urcrnry=maprange[1][1],\ 1375 resolution=None, ax=self.ax) 1376 1377 # turn on 'fulldisk' property when using full disk 1378 if maprange == [[-width/2, -height/2], [width/2, height/2]]: 1379 m._fulldisk = True 1380 1381 xrange = (m.llcrnrx, m.urcrnrx) 1382 yrange = (m.llcrnry, m.urcrnry) 1383 1384 # set colorbar limits 1385 if clim: 1386 cmin,cmax = clim 1387 else: 1388 try: 1389 cmin = min([min(z_vals) for z_vals in self.z_data_sets])*0.99 1390 cmax = max([max(z_vals) for z_vals in self.z_data_sets])*1.01 1391 # catch no triggers 1392 except ValueError: 1393 cmin = 1 1394 cmax = 10 1395 1396 # reset logs 1397 if logz: 1398 cmin = numpy.math.log(cmin, base) 1399 cmax = numpy.math.log(cmax, base) 1400 1401 # make plot 1402 for x_vals, y_vals, z_vals, plot_kwargs in\ 1403 itertools.izip(self.x_data_sets, self.y_data_sets, self.z_data_sets,\ 1404 self.kwarg_sets): 1405 1406 if logz: z_vals = [numpy.math.log(z, base) for z in z_vals] 1407 1408 plot_kwargs.setdefault("marker", "o") 1409 plot_kwargs.setdefault("edgecolor", "k") 1410 plot_kwargs.setdefault("vmin", cmin) 1411 plot_kwargs.setdefault("vmax", cmax) 1412 1413 # sort data by z-value 1414 zipped = zip(x_vals, y_vals, z_vals) 1415 zipped.sort(key=lambda (x,y,z): z) 1416 x_vals, y_vals, z_vals = map(list, zip(*zipped)) 1417 1418 # project data 1419 if projection in ['moll', 'hammer', 'orth']: 1420 convert = lambda x: (x>=0 and x) or (x<0 and x+360) 1421 else: 1422 convert = lambda x: (x>=180 and x-360) or (x<180 and x) 1423 x_vals = [convert(x) for x in x_vals] 1424 x_vals, y_vals = m(x_vals, y_vals) 1425 1426 # pull out z value if given 1427 if plot_kwargs.has_key('facecolor'): 1428 self.ax.scatter(x_vals, y_vals, **plot_kwargs) 1429 else: 1430 p.append(self.ax.scatter(x_vals, y_vals, c=z_vals, **plot_kwargs)) 1431 1432 # finish projection 1433 m.drawmapboundary() 1434 1435 # set labels 1436 if projection in ['ortho']: 1437 plabels = [0, 0, 0, 0] 1438 mlabels = [0, 0, 0, 0] 1439 else: 1440 plabels = [1, 0, 0, 0] 1441 mlabels = [0, 0, 0, 1] 1442 1443 # draw parallels 1444 parallels = numpy.arange(-90., 120., 30.) 1445 m.drawparallels(parallels, labels=plabels,\ 1446 labelstyle='+/-', latmax=90) 1447 1448 # draw meridians 1449 if projection in ['moll', 'hammer', 'ortho']: 1450 meridians = numpy.arange(0., 360., 45.) 1451 else: 1452 meridians = numpy.arange(-180, 181, 45) 1453 m.drawmeridians(meridians, labels=mlabels,\ 1454 latmax=90, labelstyle='+/-') 1455 1456 # label parallels for certain projections 1457 if projection in ['ortho']: 1458 for lon,lat in zip([0.]*len(parallels), parallels): 1459 x, y = m(lon, lat) 1460 if x<=10**20 and y<=10**20\ 1461 and xrange[0]<x<xrange[1] and yrange[0]<=y<=yrange[1]: 1462 self.ax.text(x, y, r"$%0.0f^\circ$" % lat) 1463 # label meridians for certain projections 1464 if projection in ['moll', 'hammer', 'ortho']: 1465 for lon,lat in zip(meridians, [0.]*len(meridians)): 1466 tlon = lon 1467 while tlon < 0: tlon += 360 1468 x, y = m(lon, lat) 1469 if x<=10**20 and y<=10**20\ 1470 and xrange[0]<x<xrange[1] and yrange[0]<=y<=yrange[1]: 1471 self.ax.text(x, y, r"$%0.0f^\circ$" % tlon) 1472 1473 # write colorbar 1474 self.set_colorbar(p[-1], [cmin, cmax], logz, base) 1475 1476 # add legend if there are any non-trivial labels 1477 self.add_legend_if_labels_exist(loc=loc, scatterpoints=1)
1478
1479 # ============================================================================= 1480 # Wrappers to translate ligolw table into plot 1481 # ============================================================================= 1482 1483 # ============================================================================= 1484 # Plot time series of data set 1485 1486 -def plot_data_series(data, outfile, x_format='time', zero=None, \ 1487 zeroindicator=False, **kwargs):
1488 1489 """ 1490 Plot the time series / spectrum of a given set (or given sets) of data. 1491 1492 Arguments: 1493 1494 data : list 1495 list of (ChannelName,x_data,y_data) tuples with channel name (or data 1496 source) and time/freq, amplitude arrays for each channel. Channels are 1497 plotted in the order given. 1498 outfile : str 1499 output plot path 1500 1501 Keyword Arguments: 1502 1503 x_format : [ 'time' | 'frequency' ] 1504 type of data for x_axis, allows formatting of axes 1505 zero : [ float | int | LIGOTimeGPS ] 1506 time around which to centre time series plot 1507 zeroindicator : [ False | True ] 1508 indicate zero time with veritcal dashed line, default: False 1509 1510 Unnamed keyword arguments: 1511 1512 logx : [ True | False ] 1513 boolean option to display x-axis in log scale. 1514 logy : [ True | False ] 1515 boolean option to display y-axis in log scale. 1516 xlim : tuple 1517 (xmin, xmax) limits for x-axis 1518 ylim : tuple 1519 (ymin, ymax) limits for y-axis 1520 xlabel : string 1521 label for x-axis 1522 ylabel : string 1523 label for y-axis 1524 title : string 1525 title for plot 1526 subtitle : string 1527 subtitle for plot 1528 1529 All other given arguments will be passed to matplotlib.axes.Axes.plot. 1530 """ 1531 1532 # format times 1533 if not kwargs.has_key('xlim') or not kwargs['xlim']: 1534 start = min([min(d[1]) for d in data]) 1535 end = max([max(d[1]) for d in data]) 1536 kwargs['xlim'] = [start,end] 1537 1538 if not zero: 1539 zero = kwargs['xlim'][0] 1540 1541 # set plot time unit whether it's used or not 1542 if x_format=='time': 1543 unit, timestr = time_unit(kwargs['xlim'][1]-kwargs['xlim'][0]) 1544 1545 # set labels 1546 if x_format=='time': 1547 zero = LIGOTimeGPS('%.3f' % zero) 1548 if zero.nanoseconds==0: 1549 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero))[:6])\ 1550 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 1551 else: 1552 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero.seconds))[:6])\ 1553 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 1554 tlabel = tlabel.replace(' UTC', '.%.3s UTC' % zero.nanoseconds) 1555 xlabel = kwargs.pop('xlabel',\ 1556 'Time (%s) since %s (%s)' % (timestr, tlabel, zero)) 1557 title = kwargs.pop('title', 'Time series') 1558 else: 1559 xlabel = kwargs.pop('xlabel', 'Frequency (Hz)') 1560 title = kwargs.pop('title', 'Frequency spectrum') 1561 1562 ylabel = kwargs.pop('ylabel', 'Amplitude') 1563 subtitle = kwargs.pop('subtitle', '') 1564 1565 # customise plot appearance 1566 set_rcParams() 1567 1568 # get limits 1569 xlim = kwargs.pop('xlim', None) 1570 ylim = kwargs.pop('ylim', None) 1571 calendar_time = kwargs.pop('calendar_time', False) 1572 1573 # get axis scales 1574 logx = kwargs.pop('logx', False) 1575 logy = kwargs.pop('logy', False) 1576 1577 # get legend loc 1578 loc = kwargs.pop('loc', 'best') 1579 1580 # get colorbar options 1581 hidden_colorbar = kwargs.pop('hidden_colorbar', False) 1582 1583 # get savefig option 1584 bbox_inches = kwargs.pop('bbox_inches', 'tight') 1585 1586 # generate plot object 1587 plot = DataPlot(xlabel, ylabel, title, subtitle) 1588 1589 # set plot params 1590 style = kwargs.pop('style', '-') 1591 if style in ['-', '--', '-.', ':']: 1592 kwargs.setdefault('linestyle', style) 1593 kwargs.setdefault('linewidth', 2) 1594 kwargs.pop('marker', None) 1595 else: 1596 kwargs.setdefault('marker', style) 1597 kwargs.setdefault('markersize', 5) 1598 kwargs.setdefault('linewidth', 0) 1599 kwargs.pop('linestyle', ' ') 1600 1601 # get uniq data sets that aren't errors 1602 allchannels = [] 1603 channels = [] 1604 for i,(c,_,_) in enumerate(data): 1605 allchannels.append(c) 1606 if not re.search('(min|max)\Z', str(c)): 1607 channels.append((i,c)) 1608 1609 # add data 1610 for i,c in channels: 1611 x_data,y_data = data[i][1:] 1612 if x_format=='time': 1613 if calendar_time: 1614 x_data = gps2datenum(numpy.array(map(float, x_data))) 1615 else: 1616 x_data = (x_data.astype(float)-float(zero))/unit 1617 lab = str(c) 1618 if lab != '_': lab = lab.replace('_', '\_') 1619 plot.add_content(x_data, y_data, label=lab,**kwargs) 1620 1621 # finalize plot 1622 plot.finalize(logx=logx, logy=logy, loc=loc, hidden_colorbar=hidden_colorbar) 1623 1624 # plot errors 1625 l = None 1626 for (i,channel),c in itertools.izip(channels, plotutils.default_colors()): 1627 try: 1628 minidx = allchannels.index('%s_min' % str(channel)) 1629 maxidx = allchannels.index('%s_max' % str(channel)) 1630 except ValueError: 1631 continue 1632 y = [] 1633 for idx in [minidx,maxidx]: 1634 x_data,y_data = data[idx][1:] 1635 y.append(y_data) 1636 if x_format=='time': 1637 x_data = (numpy.array(map(float, x_data))-float(zero))/unit 1638 l = l!=None and l or float(kwargs.pop('linewidth', 1))/4 1639 plot.ax.plot(x_data, y_data, color=c, linewidth=l, **kwargs) 1640 plot.ax.fill_between(x_data, y[1], y[0], color=c, alpha=0.25) 1641 1642 # set axes 1643 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True) 1644 if ylim: 1645 plot.ax.set_ylim(tuple(ylim)) 1646 1647 # FIXME add zero indicator 1648 if x_format=='time': 1649 axis_lims = plot.ax.get_ylim() 1650 if zeroindicator: 1651 plot.ax.plot([0, 0], [axis_lims[0], axis_lims[1]], 'r--', linewidth=2) 1652 plot.ax.set_ylim([ axis_lims[0], axis_lims[1] ]) 1653 1654 # set x axis 1655 if xlim and calendar_time: 1656 plot.ax.set_xlim([gps2datenum(float(xlim[0])),\ 1657 gps2datenum(float(xlim[1]))]) 1658 elif xlim: 1659 plot.ax.set_xlim([ float(xlim[0]-zero)/unit, float(xlim[1]-zero)/unit ]) 1660 else: 1661 # set global axis limits 1662 if xlim: 1663 plot.ax.set_xlim(xlim) 1664 if ylim: 1665 plot.ax.set_ylim(ylim) 1666 1667 1668 set_ticks(plot.ax, calendar_time=calendar_time) 1669 1670 plot.savefig(outfile, bbox_inches=bbox_inches,\ 1671 bbox_extra_artists=plot.ax.texts) 1672 pylab.close(plot.fig)
1673
1674 # ============================================================================= 1675 # Plot a histogram of any column 1676 1677 -def plot_trigger_hist(triggers, outfile, column='snr', num_bins=1000,\ 1678 color_column=None, color_bins=10,\ 1679 seglist=None, flag='unknown', start=None, end=None,\ 1680 livetime=None, etg=None, **kwargs):
1681 1682 """ 1683 Wrapper for dqPlotUtils.LineHistogram to plot a histogram of the value in 1684 any column of the ligolw table triggers. If a glue.segments.segmentlist 1685 seglist is given, the histogram is presented before and after removal of 1686 triggers falling inside any segment in the list. 1687 1688 Arguments: 1689 1690 triggers : glue.ligolw.table.Table 1691 ligolw table containing triggers 1692 outfile : string 1693 string path for output plot 1694 1695 Keyword arguments: 1696 1697 column : string 1698 valid column of triggers table to plot as histrogram 1699 num_bins : int 1700 number of histogram bins to use 1701 seglist : glue.segments.segmentlist 1702 list of segments with which to veto triggers 1703 flag : string 1704 display name of segmentlist, normally the name of the DQ flag 1705 start : [ float | int | LIGOTimeGPS] 1706 GPS start time (exclude triggers and segments before this time) 1707 end : [ float | int | LIGOTimeGPS] 1708 GPS end time (exclude triggers and segments after this time) 1709 livetime : [ float | int | LIGOTimeGPS ] 1710 span of time from which triggers and segments are valid, used to 1711 display histogram counts in terms of rate (Hz) for easy comparisons 1712 etg : string 1713 display name of trigger generator, defaults based on triggers tableName 1714 1715 Unnamed keyword arguments: 1716 1717 cumulative : [ True | False ] 1718 plot cumulative histogram 1719 rate : [ True | False ] 1720 plot rate histogram (normalises with given or calculated livetime) 1721 fill : [ True | False ] 1722 fill below the histogram curves, default colors: 1723 red (vetoed), green (not vetoed). 1724 logx : [ True | False ] 1725 boolean option to display x-axis in log scale. 1726 logy : [ True | False ] 1727 boolean option to display y-axis in log scale. 1728 xlim : tuple 1729 (xmin, xmax) limits for x-axis 1730 ylim : tuple 1731 (ymin, ymax) limits for y-axis 1732 xlabel : string 1733 label for x-axis 1734 ylabel : string 1735 label for y-axis 1736 title : string 1737 title for plot 1738 subtitle : string 1739 subtitle for plot 1740 greyscale : [ True | False ] 1741 use (non-greyscale) colour scheme suitable for greyscale plots 1742 1743 All other given arguments will be passed to matplotlib.axes.Axes.plot and 1744 matplotlib.axes.Axes.fill_between. 1745 """ 1746 1747 get_time = def_get_time(triggers.tableName) 1748 1749 # calculate livetime 1750 if not start or not end: 1751 times = [get_time(t) for t in triggers] 1752 if not start: 1753 start = min(times) 1754 if not end: 1755 end = max(times) 1756 if not livetime: 1757 livetime = end-start 1758 livetime = float(livetime) 1759 1760 # format seglist 1761 if seglist==None: 1762 seglist = segments.segmentlist() 1763 else: 1764 seglist = segments.segmentlist(seglist) 1765 1766 # format colours 1767 logz = kwargs.pop('logz', False) 1768 clim = kwargs.pop('clim', None) 1769 1770 if color_column: 1771 colData = get_column(triggers, color_column).astype(float) 1772 if not clim: 1773 if len(colData)>0: 1774 clim = [colData.min(), colData.max()] 1775 else: 1776 clim = [1,10] 1777 if isinstance(color_bins, int): 1778 num_color_bins = color_bins 1779 if logz: 1780 color_bins = numpy.logspace(numpy.math.log10(clim[0]),\ 1781 numpy.math.log10(clim[1]),\ 1782 num=num_color_bins) 1783 else: 1784 color_bins = numpy.linspace(0, clim[1],\ 1785 num=num_color_bins) 1786 else: 1787 num_color_bins = len(color_bins) 1788 else: 1789 num_color_bins = 1 1790 1791 # get data 1792 tdata = get_column(triggers, 'time') 1793 preData = [] 1794 postData = [] 1795 for i in range(num_color_bins): 1796 if color_column: 1797 preData.append(get_column(triggers, column).astype(float)\ 1798 [colData>=color_bins[i]]) 1799 else: 1800 preData.append(get_column(triggers, column).astype(float)) 1801 postData.append([p for j,p in enumerate(preData[i])\ 1802 if tdata[j] not in seglist]) 1803 1804 # get veto livetime 1805 vetoLivetime = livetime-float(abs(seglist)) 1806 1807 # set some random plot parameters 1808 greyscale = kwargs.pop('greyscale', False) 1809 if seglist: 1810 color = ['r','g'] 1811 label = ['Before vetoes', 'After vetoes'] 1812 else: 1813 color = ['b'] 1814 label = ['_'] 1815 if greyscale: 1816 color = ['k','k'] 1817 linestyle = ['-','--'] 1818 else: 1819 linestyle = ['-','-'] 1820 1821 # fix names for latex 1822 if flag: 1823 flag = flag.replace('_','\_') 1824 else: 1825 flag = 'Unknown' 1826 1827 # customise plot appearance 1828 set_rcParams() 1829 1830 # get limits 1831 xlim = kwargs.pop('xlim', None) 1832 ylim = kwargs.pop('ylim', None) 1833 1834 # get axis scales 1835 logx = kwargs.pop('logx', False) 1836 logy = kwargs.pop('logy', False) 1837 1838 # get colorbar options 1839 hidden_colorbar = kwargs.pop('hidden_colorbar', False) 1840 1841 # get savefig option 1842 bbox_inches = kwargs.pop('bbox_inches', 'tight') 1843 1844 # get fill 1845 fill = kwargs.pop('fill', False) 1846 1847 # get extras 1848 cumulative = kwargs.pop('cumulative', False) 1849 rate = kwargs.pop('rate', False) 1850 1851 # set labels 1852 xlabel = kwargs.pop('xlabel', display_name(column)) 1853 if rate and cumulative: 1854 ylabel = kwargs.pop('ylabel', 'Cumulative rate (Hz)') 1855 elif rate: 1856 ylabel = kwargs.pop('ylabel', 'Rate (Hz)') 1857 elif not rate and cumulative: 1858 ylabel = kwargs.pop('ylabel', 'Cumulative number') 1859 elif not rate and not cumulative: 1860 ylabel = kwargs.pop('ylabel', 'Number') 1861 1862 # get ETG 1863 if not etg: 1864 if re.search('burst', triggers.tableName.lower()): 1865 etg = 'Burst' 1866 elif re.search('inspiral', triggers.tableName.lower()): 1867 etg = 'Inspiral' 1868 elif re.search('ringdown', triggers.tableName.lower()): 1869 etg = 'Ringdown' 1870 else: 1871 etg = 'Unknown' 1872 1873 title = '%s triggers' % (etg.replace('_','\_')) 1874 if seglist: 1875 title += ' and %s segments' % (flag) 1876 title = kwargs.pop('title', title) 1877 if start and end: 1878 subtitle = '%d-%d' % (start, math.ceil(end)) 1879 else: 1880 subtitle = "" 1881 subtitle = kwargs.pop('subtitle', subtitle) 1882 zlabel = kwargs.pop('zlabel',\ 1883 color_column and display_name(color_column) or "") 1884 1885 # generate plot object 1886 plot = LineHistogram(xlabel, ylabel, title, subtitle, zlabel) 1887 1888 # add each data set 1889 if color_column: 1890 for i in range(len(preData)): 1891 plot.add_content(preData[i], livetime=livetime, cval=color_bins[i],\ 1892 linestyle=linestyle[0], label=label[0], **kwargs) 1893 else: 1894 plot.add_content(preData[0], livetime=livetime, color=color[0],\ 1895 linestyle=linestyle[0], label=label[0], **kwargs) 1896 if seglist: 1897 plot.add_content(postData[0], livetime=vetoLivetime, color=color[1],\ 1898 linestyle=linestyle[1], label=label[1], **kwargs) 1899 1900 1901 1902 # finalize plot with histograms 1903 if not num_bins: num_bins=100 1904 if color_column: 1905 colorbar = logz and 'log' or 'linear' 1906 else: 1907 colorbar = None 1908 plot.finalize(num_bins=num_bins, logx=logx, logy=logy, cumulative=cumulative,\ 1909 rate=rate, fill=fill, colorbar=colorbar,\ 1910 hidden_colorbar=hidden_colorbar) 1911 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True) 1912 1913 # set lower y axis limit 1914 if rate: 1915 ymin = 1/livetime 1916 elif logy: 1917 ymin = 0.5 1918 else: 1919 ymin = plot.ax.get_ylim()[0] 1920 plot.ax.set_ybound(lower=ymin) 1921 1922 # set global axis limits 1923 if xlim: 1924 plot.ax.set_xlim(xlim) 1925 if ylim: 1926 plot.ax.set_ylim(ylim) 1927 1928 # set global ticks 1929 set_ticks(plot.ax) 1930 1931 # save figure 1932 plot.savefig(outfile, bbox_inches=bbox_inches,\ 1933 bbox_extra_artists=plot.ax.texts)
1934
1935 # ============================================================================= 1936 # Plot one column against another column coloured by any third column 1937 1938 -def plot_triggers(triggers, outfile, reftriggers=None, xcolumn='time',\ 1939 ycolumn='snr', zcolumn=None, rankcolumn=None, etg=None, start=None, end=None,\ 1940 zero=None, seglist=None, flag=None, **kwargs):
1941 1942 """ 1943 Plots ycolumn against xcolumn for columns in given 1944 Sngl{Burst,Inspiral}Table object triggers, coloured by the zcolumn 1945 highlighting those entries falling inside one of the entries in the 1946 glue.segments.segmentlist object segments, if given. 1947 1948 'time' given as a column name is a special case, since s and ns times are 1949 stored separately in the SnglTable structures. In this case the 1950 trigger.get_xxxx() function is called. 1951 1952 Arguments: 1953 1954 triggers : glue.ligolw.table.Table 1955 ligolw table containing triggers 1956 outfile : string 1957 string path for output plot 1958 1959 Keyword arguments: 1960 1961 xcolumn : string 1962 valid column of triggers table to plot on x-axis 1963 ycolumn : string 1964 valid column of triggers table to plot on y-axis 1965 zcolumn : string 1966 valid column of triggers table to use for colorbar (optional). 1967 rankcolumn : string 1968 valid column of triggers table to use for ranking events (optional). 1969 etg : string 1970 display name of trigger generator, defaults based on triggers tableName 1971 start : [ float | int | LIGOTimeGPS ] 1972 GPS start time of plot 1973 end : [ float | int | LIGOTimeGPS ] 1974 GPS end time of plot 1975 zero : [ float | int | LIGOTimeGPS ] 1976 time around which to centre plot 1977 seglist : glue.segments.segmentlist 1978 list of segments with which to veto triggers 1979 flag : string 1980 display name of segmentlist, normally the name of the DQ flag 1981 1982 Unnamed keyword arguments: 1983 1984 detchar : [ True | False ] 1985 use 'DetChar' style for scatter plot with colorbar, triggers below given 1986 dcthreshold are small with no edges, whilst other triggers are normal 1987 dcthreshold : float 1988 threshold below which scatter points are small with no edges when using 1989 DetChar plotting style 1990 logx : [ True | False ] 1991 boolean option to display x-axis in log scale. 1992 logy : [ True | False ] 1993 boolean option to display y-axis in log scale. 1994 logz : [ True | False ] 1995 boolean option to display z-axis in log scale. 1996 xlim : tuple 1997 (xmin, xmax) limits for x-axis. Triggers outside range are removed. 1998 ylim : tuple 1999 (ymin, ymax) limits for y-axis. Triggers outside range are removed. 2000 zlim : tuple 2001 (zmin, zmax) limits for z-axis. Triggers outside range are removed. 2002 clim : tuple 2003 (cmin, cmax) limits for color scale. Triggers outside range are moved 2004 onto boundary. 2005 xlabel : string 2006 label for x-axis 2007 ylabel : string 2008 label for y-axis 2009 zlabel : string 2010 label for z-axis 2011 title : string 2012 title for plot 2013 subtitle : string 2014 subtitle for plot 2015 greyscale : [ True | False ] 2016 use (non-greyscale) colour scheme suitable for greyscale plots 2017 2018 All other given arguments will be passed to matplotlib.axes.Axes.scatter. 2019 """ 2020 2021 from pylal import plotutils 2022 2023 # test multiple tables 2024 if not len(triggers)==0 and \ 2025 (isinstance(triggers[0], tuple) or isinstance(triggers[0], list)): 2026 assert not zcolumn,\ 2027 "Can only plot single table when using colorbar plot" 2028 tables = [t[1] for t in triggers] 2029 tablelabel = [t[0] for t in triggers] 2030 for i,t in enumerate(tablelabel): 2031 if t!='_': 2032 tablelabel[i] = t.replace('_','\_') 2033 else: 2034 tables = [triggers] 2035 tablelabel = '_' 2036 2037 # get time column 2038 get_time = [] 2039 for t in tables: 2040 get_time.append(def_get_time(t.tableName)) 2041 2042 # set start and end time if needed 2043 if not start or not end: 2044 times = [get_time[i](t) for i in xrange(len(tables)) for t in tables[i]] 2045 if not start and len(times)>=1: 2046 start = int(math.floor(min(times))) 2047 elif not start: 2048 start = 0 2049 if not end and len(times)>=1: 2050 end = int(math.ceil(max(times))) 2051 elif not end: 2052 end = 1 2053 2054 # set zero time 2055 if not zero: 2056 zero = start 2057 2058 # get time params 2059 unit,timestr = time_unit(end-start) 2060 2061 # set up segmentlist 2062 if seglist: 2063 segs = segments.segmentlist(seglist) 2064 if not seglist: 2065 segs = segments.segmentlist() 2066 2067 # get axis limits 2068 xlim = kwargs.pop('xlim', None) 2069 ylim = kwargs.pop('ylim', None) 2070 zlim = kwargs.pop('zlim', None) 2071 clim = kwargs.pop('clim', zlim) 2072 calendar_time = kwargs.pop('calendar_time', None) 2073 2074 # set up columns 2075 columns = list(map(str.lower, [xcolumn, ycolumn])) 2076 if zcolumn: 2077 columns.append(zcolumn.lower()) 2078 if rankcolumn: 2079 columns.append(rankcolumn.lower()) 2080 else: columns.append(zcolumn.lower()) 2081 2082 # set up limits 2083 limits = [xlim, ylim, zlim, None] 2084 for i,col in enumerate(columns): 2085 if re.search('time\Z', col) and not limits[i]: 2086 limits[i] = [start,end] 2087 2088 # get veto info 2089 if seglist: 2090 tdata = get_column(triggers, 'time') 2091 2092 # get all data 2093 vetoData = [] 2094 nvetoData = [] 2095 for j,tab in enumerate(tables): 2096 vetoData.append({}) 2097 nvetoData.append({}) 2098 # get veto info 2099 if seglist: 2100 tdata = get_column(tab, 'time') 2101 # get data 2102 for i,col in enumerate(columns): 2103 nvetoData[j][col] = get_column(tab, col).astype(float) 2104 # apply limits and vetoes 2105 condition = True 2106 for i,col in enumerate(columns): 2107 if limits[i]: 2108 condition = condition & (limits[i][0] <= nvetoData[j][col])\ 2109 & (nvetoData[j][col] <= limits[i][1]) 2110 for col in nvetoData[j].keys(): 2111 nvetoData[j][col] = nvetoData[j][col][condition] 2112 if seglist: 2113 vetoData[j][col] = numpy.asarray([d for i,d in\ 2114 enumerate(nvetoData[j][col]) if\ 2115 tdata[i] in seglist]) 2116 else: 2117 vetoData[j][col] = numpy.array([]) 2118 2119 data = {} 2120 2121 2122 2123 # normalize zcolumn by time-averaged value 2124 whitenedFlag = kwargs.pop('whitened', False) 2125 if zcolumn and whitenedFlag: 2126 # get ref data if provided 2127 refData = {} 2128 if reftriggers: 2129 for i,col in enumerate(columns): 2130 refData[col] = get_column(reftriggers, col).astype(float) 2131 else: 2132 for i,col in enumerate(columns): 2133 refData[col] = nvetoData[0][col] 2134 2135 uniqYvalues = numpy.unique1d(refData[ycolumn]) 2136 # building look back table by hand, is included in unique1d for numpy >= v1.3 2137 for yVal in uniqYvalues: 2138 backTable = numpy.where(yVal == nvetoData[0][ycolumn]) 2139 zMedian = numpy.median(refData[zcolumn][yVal == refData[ycolumn]]) 2140 for iTrig in backTable[0]: 2141 nvetoData[0][zcolumn][iTrig] /= zMedian 2142 2143 # filter zcolumn by provided poles/zeros filter as a function of ycolumn 2144 flatenedFlag = kwargs.pop('filter', False) 2145 if zcolumn and flatenedFlag: 2146 # get filter params 2147 polesList = kwargs.pop('poles', None) 2148 zerosList = kwargs.pop('zeros', None) 2149 amplitude = kwargs.pop('amplitude', 1) 2150 nvetoData[0][zcolumn] *= amplitude 2151 for filtPole in polesList: 2152 nvetoData[0][zcolumn] /= abs(nvetoData[0][ycolumn] - filtPole) 2153 for filtZero in zerosList: 2154 nvetoData[0][zcolumn] *= abs(nvetoData[0][ycolumn] - filtZero) 2155 nvetoData[0][zcolumn].astype(float) 2156 2157 # flaten zcolumn by 1/sqrt of sum given rational fraction mononomes as a 2158 # function of ycolumn 2159 flatenedFlag = kwargs.pop('flaten', False) 2160 if zcolumn and flatenedFlag: 2161 # get filter params 2162 expList = kwargs.pop('exponents', None) 2163 constList = kwargs.pop('constants', None) 2164 filter = numpy.zeros(len(nvetoData[0][zcolumn])) 2165 for iTerm, exponent in enumerate(expList): 2166 filter += pow(constList[iTerm]*numpy.power(nvetoData[0][ycolumn],expList[iTerm]),2) 2167 filter = numpy.sqrt(filter) 2168 nvetoData[0][zcolumn] /= filter 2169 2170 # median/min/max of ycolumn binned by exact xcolumn values 2171 minmaxmedianFlag = kwargs.pop('minmaxmedian', False) 2172 if minmaxmedianFlag: 2173 uniqXvalues = numpy.unique1d(nvetoData[j][xcolumn]) 2174 # building look back table by hand, is included in unique1d for numpy >= v1.3 2175 for xVal in uniqXvalues: 2176 backTable = numpy.where(xVal == nvetoData[j][xcolumn]) 2177 if len(backTable[0]) > 3: 2178 nvetoData[j][ycolumn][backTable[0][0]] =\ 2179 numpy.median(nvetoData[j][ycolumn][xVal == nvetoData[j][xcolumn]]) 2180 nvetoData[j][ycolumn][backTable[0][1]] =\ 2181 numpy.min(nvetoData[j][ycolumn][xVal == nvetoData[j][xcolumn]]) 2182 nvetoData[j][ycolumn][0][backTable[0][2]] =\ 2183 numpy.max(nvetoData[j][ycolumn][xVal == nvetoData[j][xcolumn]]) 2184 for iTrig in backTable[0][3:]: 2185 nvetoData[j][ycolumn][iTrig] = numpy.nan 2186 2187 # down-sample (xaxis) the triggers by plotting only median z-value over averaging window 2188 medianDuration = float(kwargs.pop('medianduration', 0)) 2189 stdDuration = float(kwargs.pop('stdduration', 0)) 2190 if medianDuration or stdDuration: 2191 uniqYvalues = numpy.unique1d(nvetoData[0][ycolumn]) 2192 if medianDuration: 2193 avDuration = medianDuration 2194 elif stdDuration: 2195 avDuration = stdDuration 2196 tedges = numpy.arange(float(start), float(end), avDuration) 2197 tedges = numpy.append(tedges, float(end)) 2198 # array for repacking triggers into time-frequency bins 2199 repackTrig = {} 2200 for yVal in uniqYvalues: 2201 repackTrig[yVal] = {} 2202 for iBin in range(len(tedges)-1): 2203 repackTrig[yVal][iBin] = [] 2204 # building new set of triggers 2205 newTrigs = {} 2206 newTrigs[xcolumn] = [] 2207 newTrigs[ycolumn] = [] 2208 newTrigs[zcolumn] = [] 2209 for iTrig in range(len(nvetoData[0][zcolumn])): 2210 trigVal = nvetoData[0][ycolumn][iTrig] 2211 trigBin = math.floor((nvetoData[0][xcolumn][iTrig]-float(start))/avDuration) 2212 repackTrig[trigVal][trigBin].append(nvetoData[0][zcolumn][iTrig]) 2213 for yVal in uniqYvalues: 2214 for iBin in range(len(tedges)-1): 2215 # keep only if at least a few elements, std doesn't make sens otherwise 2216 if medianDuration and len(repackTrig[yVal][iBin]) > 0: 2217 newTrigs[xcolumn].append( (tedges[iBin]+tedges[iBin+1])/2 ) 2218 newTrigs[ycolumn].append(yVal) 2219 newTrigs[zcolumn].append(numpy.median(numpy.array(repackTrig[yVal][iBin]))) 2220 elif stdDuration and len(repackTrig[yVal][iBin]) > 5: 2221 newTrigs[xcolumn].append( (tedges[iBin]+tedges[iBin+1])/2 ) 2222 newTrigs[ycolumn].append(yVal) 2223 newTrigs[zcolumn].append(numpy.std(numpy.array(repackTrig[yVal][iBin]))/ \ 2224 numpy.median(numpy.array(repackTrig[yVal][iBin]))) 2225 if newTrigs[zcolumn][-1] == 0: 2226 newTrigs[zcolumn][-1] = numpy.nan 2227 for i,col in enumerate(columns): 2228 nvetoData[0][col] = numpy.array(newTrigs[col]) 2229 2230 # get limits 2231 for i,col in enumerate(columns): 2232 if not limits[i]: 2233 limits[i] = [0,0] 2234 for j in xrange(len(tables)): 2235 data[col] = numpy.concatenate((nvetoData[j][col], vetoData[j][col])) 2236 if len(data[col])>=1: 2237 limits[i][0] = min(data[col].min()*0.99, limits[i][0]) 2238 limits[i][1] = max(data[col].max()*1.01, limits[i][1]) 2239 2240 # renormalise time and set time axis label unless given 2241 if re.search('time\Z', col): 2242 renormalise = True 2243 if kwargs.has_key('xlabel'): renormalise = False 2244 if renormalise: 2245 for j in xrange(len(tables)): 2246 if calendar_time: 2247 vetoData[j][col] = gps2datenum(vetoData[j][col]) 2248 nvetoData[j][col] = gps2datenum(nvetoData[j][col]) 2249 else: 2250 vetoData[j][col] = (vetoData[j][col]-float(zero))/unit 2251 nvetoData[j][col] = (nvetoData[j][col]-float(zero))/unit 2252 if calendar_time: 2253 limits[i] = [gps2datenum(float(limits[i][0])),\ 2254 gps2datenum(float(limits[i][1]))] 2255 else: 2256 limits[i] = [float(limits[i][0]-zero)/unit,\ 2257 float(limits[i][1]-zero)/unit] 2258 2259 # set labels 2260 label = {} 2261 for i,col in enumerate(['xcolumn', 'ycolumn', 'zcolumn']): 2262 if i >= len(columns): continue 2263 if re.search('time\Z', columns[i]): 2264 zerostr = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero))[:6])\ 2265 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 2266 lab = 'Time (%s) since %s (%s)' % (timestr, zerostr, zero) 2267 label[columns[i]] = kwargs.pop('%slabel' % col[0], lab) 2268 else: 2269 label[columns[i]] = kwargs.pop('%slabel' % col[0],\ 2270 display_name(columns[i])) 2271 2272 # find loudest event 2273 loudest = {} 2274 if len(columns)==4 and\ 2275 len(nvetoData[0][columns[0]])+len(vetoData[0][columns[0]])>=1: 2276 # find loudest vetoed event 2277 vetomax = 0 2278 if len(vetoData[0][columns[3]])>=1: 2279 vetomax = vetoData[0][columns[3]].max() 2280 nvetomax = 0 2281 # find loudest unvetoed event 2282 if len(nvetoData[0][columns[3]])>=1: 2283 nvetomax = nvetoData[0][columns[3]].max() 2284 if vetomax == nvetomax == 0: 2285 pass 2286 # depending on which one is loudest, find loudest overall event 2287 elif vetomax > nvetomax: 2288 index = vetoData[0][columns[3]].argmax() 2289 for col in columns: 2290 loudest[col] = vetoData[0][col][index] 2291 else: 2292 index = nvetoData[0][columns[3]].argmax() 2293 for col in columns: 2294 loudest[col] = nvetoData[0][col][index] 2295 2296 # fix flag for use with latex 2297 if flag: 2298 flag = flag.replace('_', '\_') 2299 else: 2300 flag = 'Unknown' 2301 2302 # get ETG 2303 if not etg: 2304 if re.search('burst', tables[0].tableName.lower()): 2305 etg = 'Burst' 2306 elif re.search('inspiral', tables[0].tableName.lower()): 2307 etg = 'Inspiral' 2308 elif re.search('ringdown', tables[0].tableName.lower()): 2309 etg = 'Ringdown' 2310 else: 2311 etg = 'Unknown' 2312 else: 2313 etg = etg.replace('_', '\_') 2314 2315 # customise plot appearance 2316 set_rcParams() 2317 2318 # set title 2319 title = '%s triggers' % etg 2320 if seglist: 2321 title += ' and %s segments' % (flag) 2322 title = kwargs.pop('title', title) 2323 2324 if len(columns)==4 and loudest: 2325 subtitle = "Loudest event by %s:" % display_name(columns[-1]) 2326 for col in columns: 2327 maxcol = loudest[col] 2328 if re.search('time\Z', col) and renormalise: 2329 maxcol = maxcol*unit+zero 2330 loudstr = "%s=%.2f" % (display_name(col), maxcol) 2331 if not re.search(loudstr, subtitle): 2332 subtitle += ' %s' % loudstr 2333 elif start and end: 2334 subtitle = '%s-%s' % (start, end) 2335 else: 2336 subtitle = '' 2337 subtitle = kwargs.pop('subtitle', subtitle) 2338 2339 # get axis scales 2340 logx = kwargs.pop('logx', False) 2341 logy = kwargs.pop('logy', False) 2342 logz = kwargs.pop('logz', False) 2343 2344 # get colorbar options 2345 hidden_colorbar = kwargs.pop('hidden_colorbar', False) 2346 2347 # get savefig option 2348 bbox_inches = kwargs.pop('bbox_inches', 'tight') 2349 2350 # get detchar plot params 2351 detchar = kwargs.pop('detchar', False) 2352 dcthresh = float(kwargs.pop('dcthreshold', 10)) 2353 2354 # get greyscale param 2355 greyscale = kwargs.pop('greyscale', False) 2356 if greyscale and not kwargs.has_key('cmap'): 2357 kwargs['cmap'] = pylab.matplotlib.colors.LinearSegmentedColormap('clrs',\ 2358 pylab.matplotlib.cm.hot._segmentdata) 2359 2360 # initialise standard scatter plot 2361 if len(columns)==2: 2362 plotutils.default_colors = lambda: itertools.cycle(('b', 'r', 'g', 'c', 'm', 'y', 'k')) 2363 plot = ScatterPlot(label[columns[0]], label[columns[1]], title, subtitle) 2364 for j in xrange(len(tables)): 2365 if len(nvetoData[j][columns[0]])>=1: 2366 plot.add_content(nvetoData[j][columns[0]], nvetoData[j][columns[1]],\ 2367 label=tablelabel[j], **kwargs) 2368 # add veto triggers 2369 if len(vetoData[j][columns[0]])>=1: 2370 plot.add_content(vetoData[j][columns[0]], vetoData[j][columns[1]],\ 2371 label=tablelabel[j], marker='x', color='r') 2372 # finalise 2373 plot.finalize(logx=logx, logy=logy, hidden_colorbar=hidden_colorbar) 2374 # initialise scatter plot with colorbar 2375 elif len(columns)==4: 2376 # initialize color bar plot 2377 if detchar: 2378 plot = DetCharScatterPlot(label[columns[0]], label[columns[1]],\ 2379 label[columns[2]], title, subtitle) 2380 else: 2381 plot = ColorbarScatterPlot(label[columns[0]], label[columns[1]],\ 2382 label[columns[2]], title, subtitle) 2383 2384 # add non veto triggers 2385 if len(nvetoData[0][columns[0]])>=1: 2386 plot.add_content(nvetoData[0][columns[0]], nvetoData[0][columns[1]],\ 2387 nvetoData[0][columns[2]], nvetoData[0][columns[3]],\ 2388 **kwargs) 2389 # add veto triggers 2390 if len(vetoData[0][columns[0]])>=1: 2391 plot.add_content(vetoData[0][columns[0]], vetoData[0][columns[1]],\ 2392 vetoData[0][columns[2]], vetoData[0][columns[2]],\ 2393 marker='x', edgecolor='r',\ 2394 **kwargs) 2395 # finalise 2396 if detchar: 2397 plot.finalize(logx=logx, logy=logy, logz=logz, clim=clim,\ 2398 rankthreshold=dcthresh) 2399 else: 2400 plot.finalize(logx=logx, logy=logy, logz=logz, clim=clim) 2401 # add loudest event to plot 2402 if loudest: 2403 plot.ax.plot([loudest[columns[0]]], [loudest[columns[1]]],\ 2404 marker='*', color='gold', markersize=15) 2405 2406 # set axes 2407 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True) 2408 2409 if limits[0]: 2410 plot.ax.set_xlim(limits[0]) 2411 if limits[1]: 2412 plot.ax.set_ylim(limits[1]) 2413 2414 # reset ticks 2415 set_ticks(plot.ax, calendar_time=calendar_time) 2416 2417 # get both major and minor grid lines 2418 plot.savefig(outfile, bbox_inches=bbox_inches,\ 2419 bbox_extra_artists=plot.ax.texts) 2420 pylab.close(plot.fig)
2421
2422 # ============================================================================= 2423 # Plot a histogram of segment duration 2424 2425 -def plot_segment_hist(segs, outfile, keys=None, num_bins=100, coltype=int,\ 2426 **kwargs):
2427 2428 """ 2429 segments. 2430 Plots a histogram of segment duration for the glue.segments.segmentlist 2431 2432 Arguments: 2433 2434 segs : [ glue.segments.segmentlist | glue.segments.segmentlistdict ] 2435 list of segments with which to veto triggers, use dict for multiple 2436 datasets 2437 outfile : string 2438 string path for output plot 2439 2440 Keyword arguments: 2441 2442 flag : string 2443 display name for segments, normally the name of the DQ flag 2444 logx : [ True | False ] 2445 boolean option to display x-axis in log scale. 2446 logy : [ True | False ] 2447 boolean option to display y-axis in log scale. 2448 """ 2449 2450 # customise plot appearance 2451 set_rcParams() 2452 2453 # get limits 2454 xlim = kwargs.pop('xlim', None) 2455 ylim = kwargs.pop('ylim', None) 2456 2457 # get labels 2458 xlabel = kwargs.pop('xlabel', 'Length of segment (seconds)') 2459 ylabel = kwargs.pop('ylabel', 'Number of segments') 2460 title = kwargs.pop('title', 'Segment Duration Histogram') 2461 subtitle = kwargs.pop('subtitle', "") 2462 2463 # get axis scale 2464 logx = kwargs.pop('logx', False) 2465 logy = kwargs.pop('logy', False) 2466 2467 # get savefig option 2468 bbox_inches = kwargs.pop('bbox_inches', 'tight') 2469 2470 # get colorbar option 2471 hidden_colorbar = kwargs.pop('hidden_colorbar', False) 2472 2473 # format mutltiple segments 2474 if isinstance(segs,list): 2475 if keys and isinstance(keys, 'str'): 2476 segs = segments.segmentlistdict({keys:segs}) 2477 keys = [keys] 2478 elif keys: 2479 segs = segments.segmentlistdict({keys[0]:segs}) 2480 else: 2481 segs = segments.segmentlistdict({'_':segs}) 2482 keys = ['_'] 2483 else: 2484 if not keys: 2485 keys = segs.keys() 2486 for i,flag in enumerate(keys): 2487 keys[i] = display_name(flag) 2488 if keys[i]!=flag: 2489 segs[keys[i]] = segs[flag] 2490 del segs[flag] 2491 2492 if kwargs.has_key('color'): 2493 colors = kwargs.pop('color').split(',') 2494 else: 2495 colors = zip(*zip(range(len(keys)), plotutils.default_colors()))[-1] 2496 2497 # generate plot object 2498 plot = VerticalBarHistogram(xlabel, ylabel, title, subtitle) 2499 2500 # add each segmentlist 2501 for flag,c in zip(keys, colors): 2502 plot.add_content([float(abs(seg)) for seg in segs[flag]],\ 2503 label=flag, color=c, **kwargs) 2504 2505 # finalize plot with histograms 2506 plot.finalize(num_bins=num_bins, logx=logx, logy=logy,\ 2507 hidden_colorbar=hidden_colorbar) 2508 2509 # set limits 2510 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True) 2511 if ylim: 2512 ylim = map(float, ylim) 2513 plot.ax.set_ylim(ylim) 2514 if xlim: 2515 xlim = map(float, xlim) 2516 plot.ax.set_xlim(xlim) 2517 2518 # save figure 2519 plot.savefig(outfile, bbox_inches=bbox_inches,\ 2520 bbox_extra_artists=plot.ax.texts) 2521 pylab.close(plot.fig)
2522
2523 # ============================================================================= 2524 # Plot rate versus time in bins 2525 2526 -def plot_trigger_rate(triggers, outfile, average=600, start=None, end=None,\ 2527 zero=None, bincolumn='peak_frequency', bins=[],\ 2528 etg='Unknown', **kwargs):
2529 2530 """ 2531 Plot rate versus time for the given ligolw table triggers, binned by the 2532 given bincolumn using the bins list. 2533 2534 Arguments: 2535 2536 triggers : glue.ligolw.table 2537 LIGOLW table containing a list of triggers 2538 outfile : string 2539 string path for output plot 2540 2541 Keyword arguments: 2542 2543 average : float 2544 Length (seconds) of rate segment 2545 start : [ float | int | LIGOTimeGPS ] 2546 GPS start time 2547 end : [ float | int | LIGOTimeGPS ] 2548 GPS end time 2549 zero : [ float | int | LIGOTimeGPS ] 2550 GPS time to use for 0 on time axis 2551 bincolumn : string 2552 valid column of the trigger table to use for binning 2553 bins : list 2554 list of tuples defining the rate bins 2555 etg : string 2556 display name of trigger generator 2557 logy : [ True | False ] 2558 boolean option to display y-axis in log scale 2559 ylim : tuple 2560 (ymin, ymax) limits for rate axis 2561 """ 2562 2563 tableName = triggers.tableName.lower() 2564 get_time = def_get_time(tableName) 2565 2566 # set start and end times 2567 if not start and not end: 2568 times = [get_time(t) for t in triggers] 2569 if not start: 2570 start = min(times) 2571 if not end: 2572 end = max(times) 2573 2574 if not zero: 2575 zero = start 2576 2577 # set plot time unit whether it's used or not 2578 unit, timestr = time_unit(end-start) 2579 2580 # set ETG 2581 if not etg: 2582 if re.search('burst', tableName): 2583 etg = 'Burst' 2584 elif re.search('inspiral', tableName): 2585 etg = 'Inspiral' 2586 elif re.search('ringdown', tableName): 2587 etg = 'Ringdown' 2588 else: 2589 etg = 'Unknown' 2590 2591 # get limits 2592 calendar_time = kwargs.pop('calendar_time', False) 2593 if calendar_time: 2594 xlim = kwargs.pop('xlim', [gps2datenum(float(start)),\ 2595 gps2datenum(float(end))]) 2596 else: 2597 xlim = kwargs.pop('xlim', [float(start-zero)/unit, float(end-zero)/unit]) 2598 ylim = kwargs.pop('ylim', None) 2599 2600 # get axis scales 2601 logx = kwargs.pop('logx', False) 2602 logy = kwargs.pop('logy', False) 2603 2604 # get averaging time 2605 average = kwargs.pop('average',average) 2606 2607 # get colorbar options 2608 hidden_colorbar = kwargs.pop('hidden_colorbar', False) 2609 2610 # get savefig option 2611 bbox_inches = kwargs.pop('bbox_inches', 'tight') 2612 2613 # format ybins 2614 if not bins: 2615 bins = [[0,float('inf')]] 2616 ybins = [map(float, bin) for bin in bins] 2617 2618 # bin data 2619 tbins = {} 2620 rate = {} 2621 for bin in ybins: 2622 if calendar_time: 2623 tbins[bin[0]] = list(gps2datenum(numpy.arange(float(start), float(end),\ 2624 average))) 2625 else: 2626 tbins[bin[0]] = list(numpy.arange(0,float(end-start), average)/unit) 2627 rate[bin[0]] = list(numpy.zeros(len(tbins[bin[0]]))) 2628 2629 for trig in triggers: 2630 x = int(float(getTrigAttribute(trig, 'time')-start)//average) 2631 y = getTrigAttribute(trig, bincolumn) 2632 for bin in ybins: 2633 if bin[0] <= y < bin[1]: 2634 if x>=len(rate[bin[0]]): 2635 print "trigger after end time, something is wrong", x, len(rate[bin[0]]) 2636 continue 2637 rate[bin[0]][x] += 1/average 2638 break 2639 2640 # if logscale includes zeros, pylab.scatter will break, so remove zeros 2641 if logy: 2642 for bin in ybins: 2643 removes = 0 2644 numtbins = len(tbins[bin[0]]) 2645 for rbin in xrange(0,numtbins): 2646 if rate[bin[0]][rbin-removes]==0: 2647 rate[bin[0]].pop(rbin-removes) 2648 tbins[bin[0]].pop(rbin-removes) 2649 removes+=1 2650 2651 # set labels 2652 etg = etg.replace('_', '\_') 2653 zero = LIGOTimeGPS('%.3f' % zero) 2654 if zero.nanoseconds==0: 2655 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero))[:6])\ 2656 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 2657 else: 2658 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero.seconds))[:6])\ 2659 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 2660 tlabel = tlabel.replace(' UTC', '.%.3s UTC' % zero.nanoseconds) 2661 xlabel = kwargs.pop('xlabel',\ 2662 'Time (%s) since %s (%s)' % (timestr, tlabel, zero)) 2663 ylabel = kwargs.pop('ylabel', 'Rate (Hz)') 2664 title = kwargs.pop('title', '%s triggers binned by %s'\ 2665 % (etg, display_name(bincolumn))) 2666 if start and end: 2667 subtitle = '%s-%s' % (start, end) 2668 else: 2669 subtitle = " " 2670 subtitle = kwargs.pop('subtitle', subtitle) 2671 2672 # customise plot appearance 2673 set_rcParams() 2674 2675 # generate plot object 2676 plot = ScatterPlot(xlabel, ylabel, title, subtitle) 2677 2678 # plot rates 2679 for bin in ybins: 2680 if logy: 2681 if len(rate[bin[0]])>0: 2682 plot.add_content(tbins[bin[0]], rate[bin[0]],\ 2683 label='-'.join(map(str, bin)), **kwargs) 2684 else: 2685 plot.add_content([1],[0.1], label='-'.join(map(str, bin)),\ 2686 visible=False) 2687 else: 2688 plot.add_content(tbins[bin[0]], rate[bin[0]], label='-'.join(map(str, bin)),\ 2689 **kwargs) 2690 2691 # finalise plot 2692 plot.finalize(logx=logx, logy=logy, hidden_colorbar=hidden_colorbar) 2693 2694 # set limits 2695 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True) 2696 plot.ax.set_xlim(xlim) 2697 if ylim: 2698 plot.ax.set_ylim(ylim) 2699 2700 # normalize ticks 2701 set_ticks(plot.ax, calendar_time) 2702 2703 # save 2704 plot.savefig(outfile, bbox_inches=bbox_inches,\ 2705 bbox_extra_artists=plot.ax.texts) 2706 pylab.close(plot.fig)
2707
2708 # ============================================================================= 2709 # Plot RMS versus time in bins 2710 2711 -def plot_trigger_rms(triggers, outfile, average=600, start=None, end=None,\ 2712 zero=None, rmscolumn='snr', bincolumn='peak_frequency',\ 2713 bins=[], etg='Unknown', **kwargs):
2714 2715 """ 2716 Plot RMS versus time for the given ligolw table triggers, binned by the 2717 given bincolumn using the bins list. 2718 2719 Arguments: 2720 2721 triggers : glue.ligolw.table 2722 LIGOLW table containing a list of triggers 2723 outfile : string 2724 string path for output plot 2725 2726 Keyword arguments: 2727 2728 average : float 2729 Length (seconds) of RMS segment 2730 start : [ float | int | LIGOTimeGPS ] 2731 GPS start time 2732 end : [ float | int | LIGOTimeGPS ] 2733 GPS end time 2734 zero : [ float | int | LIGOTimeGPS ] 2735 GPS time to use for 0 on time axis 2736 rmscolumn : string 2737 valid column of the trigger table to RMS over 2738 bincolumn : string 2739 valid column of the trigger table to use for binning 2740 bins : list 2741 list of tuples defining the rate bins 2742 etg : string 2743 display name of trigger generator 2744 logy : [ True | False ] 2745 boolean option to display y-axis in log scale 2746 ylim : tuple 2747 (ymin, ymax) limits for rate axis 2748 """ 2749 2750 tableName = triggers.tableName.lower() 2751 get_time = def_get_time(tableName) 2752 2753 # set start and end times 2754 if not start and not end: 2755 times = [get_time(t) for t in triggers] 2756 if not start: 2757 start = min(times) 2758 if not end: 2759 end = max(times) 2760 2761 if not zero: 2762 zero = start 2763 2764 # set plot time unit whether it's used or not 2765 unit, timestr = time_unit(end-start) 2766 2767 # set ETG 2768 if not etg: 2769 if re.search('burst',