1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 from __future__ import division
24 import math,re,numpy,itertools,copy,matplotlib,sys,warnings
25
26
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 """
60
61 """
62 Customise the figure parameters.
63 """
64
65
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
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
97
98
99
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
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
116
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
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
147 for line in ax.get_xticklines() + ax.get_yticklines():
148 line.set_markersize(10)
149 line.set_markeredgewidth(1)
150
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
184 for i,w in enumerate(words):
185
186 if w in acro:
187 words[i] = w.upper()
188
189 elif w in unit:
190 words[i] = '$(%s)$' % w
191
192 elif w in sub:
193 words[i] = '%s$_{\mbox{\\small %s}}$' % (w[0], w[1:])
194
195 elif w in greek:
196 words[i] = '$\%s$' % w
197
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
204 else:
205 if w.isupper():
206 words[i] = w
207 else:
208 words[i] = w.title()
209
210 words[i] = re.sub('(?<!\\\\)_', '\_', words[i])
211
212 return ' '.join(words)
213
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
222 zeroGPS = 722820.0
223
224
225
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
237 datenum = gpstime/86400 + zeroGPS
238
239 return datenum
240
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
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
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
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
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
318 if col in trig.__slots__:
319 return getattr(trig, col)
320
321
322 try:
323 return eval('trig.get_%s()', col)
324
325
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
333
334
335
336
337
338 cmin, cmax = clim
339
340
341
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
378
379 self.colorbar.set_label(label)
380 self.colorbar.draw_all()
381
388
389
390 div = axes_grid.make_axes_locatable(self.ax)
391 cax = div.new_horizontal(size="2%", pad=0.05, visible=False)
392
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
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
477 if hidden_colorbar:
478 self.add_hidden_colorbar()
479
480 PlotSegmentsPlot.add_hidden_colorbar = add_hidden_colorbar
481
482
483
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
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
518 if hidden_colorbar:
519 self.add_hidden_colorbar()
520
521
522 self.add_legend_if_labels_exist(loc=loc)
523
524 leg = self.ax.get_legend()
525
526 if leg:
527 legfr = leg.get_frame()
528 legfr.set_alpha(0.5)
529
530
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
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
574 p = []
575
576
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
584 except ValueError:
585 cmin = 1
586 cmax = 10
587
588
589 if logz:
590 cmin = numpy.math.log(cmin, base)
591 cmax = numpy.math.log(cmax, base)
592
593
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
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
617
618 self.add_colorbar(p[-1], clim=[cmin, cmax], log=logz,\
619 label=self.color_label, base=10)
620
621
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
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
640
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
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
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
694 except ValueError:
695 cmin = 1
696 cmax = 10
697
698
699 if logz:
700 cmin = numpy.math.log(cmin, base)
701 cmax = numpy.math.log(cmax, base)
702
703
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
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
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
760 self.add_colorbar(p[-1], clim=[cmin, cmax], log=logz,\
761 label=self.color_label, base=10)
762
763
764
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
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
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
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
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
877
878
879
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
891 if cumulative:
892 y = y[::-1].cumsum()[::-1]
893
894
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
905 if logy:
906 y = y.astype(float)
907 numpy.putmask(y, y==0, ymin)
908
909
910 if colorbar:
911 plot_kwargs['color'] = col
912
913
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
923 if logx:
924 self.ax.set_xscale('log')
925 if logy:
926 self.ax.set_yscale('log')
927
928
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
937 self.add_legend_if_labels_exist(loc=loc)
938
939
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
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
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
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
983 v = map(int, numpy.version.version.split('.'))
984 if v[1] >= 1: width = width[:-1]
985
986
987 if logy:
988 ymin = (base**-1)*5
989 else:
990 ymin = 0
991
992
993 legends = []
994 plot_list = []
995 for i, (data_set, plot_kwargs) in \
996 enumerate(itertools.izip(self.data_sets, self.kwarg_sets)):
997
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
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
1020
1021
1022
1023 self.ax.bar(x, y, **plot_kwargs)
1024
1025
1026 if logx:
1027 self.ax.set_xscale('log')
1028 if logy:
1029 self.ax.set_yscale('log')
1030
1031
1032 self.ax.set_ybound(lower=ymin)
1033
1034
1035 self.add_legend_if_labels_exist(loc=loc)
1036
1037 leg = self.ax.get_legend()
1038 if leg: leg.get_frame().set_alpha(0.5)
1039
1040
1041
1042 if hidden_colorbar:
1043 self.add_hidden_colorbar()
1044
1045 VerticalBarHistogram.add_hidden_colorbar = add_hidden_colorbar
1046
1047
1048
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
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
1082 plot_kwargs.setdefault('markersize', 5)
1083 markersizes.append(plot_kwargs['markersize'])
1084 plot_kwargs['markersize'] = min(20, plot_kwargs['markersize']*markerscale)
1085
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
1096 self.add_legend_if_labels_exist(loc=loc)
1097 leg = self.ax.get_legend()
1098
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
1108 for i,plot in enumerate(plots):
1109 l = plot[0]
1110 l.set_markersize(markersizes[i])
1111
1112
1113 if leg:
1114 legfr = leg.get_frame()
1115 legfr.set_alpha(0.5)
1116
1117
1118 if hidden_colorbar: self.add_hidden_colorbar()
1119
1120
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
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
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
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
1191 self.add_colorbar(p[-1], clim=[cmin, cmax], log=logz,\
1192 label=self.color_label, base=10)
1193
1194
1195 if logx: self.ax.set_xscale('log')
1196 if logy: self.ax.set_yscale('log')
1197
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
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
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
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
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
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
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
1278 m.drawmapboundary()
1279
1280
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
1289 parallels = numpy.arange(-90., 120., 30.)
1290 m.drawparallels(parallels, labels=plabels,\
1291 labelstyle='+/-', latmax=90)
1292
1293
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
1302 if projection in ['ortho']:
1303 for lon,lat in zip([0.]*len(parallels), parallels):
1304
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
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
1322 self.add_legend_if_labels_exist(loc=loc, scatterpoints=1)
1323
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
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
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
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
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
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
1392 except ValueError:
1393 cmin = 1
1394 cmax = 10
1395
1396
1397 if logz:
1398 cmin = numpy.math.log(cmin, base)
1399 cmax = numpy.math.log(cmax, base)
1400
1401
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
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
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
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
1433 m.drawmapboundary()
1434
1435
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
1444 parallels = numpy.arange(-90., 120., 30.)
1445 m.drawparallels(parallels, labels=plabels,\
1446 labelstyle='+/-', latmax=90)
1447
1448
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
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
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
1474 self.set_colorbar(p[-1], [cmin, cmax], logz, base)
1475
1476
1477 self.add_legend_if_labels_exist(loc=loc, scatterpoints=1)
1478
1479
1480
1481
1482
1483
1484
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
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
1542 if x_format=='time':
1543 unit, timestr = time_unit(kwargs['xlim'][1]-kwargs['xlim'][0])
1544
1545
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
1566 set_rcParams()
1567
1568
1569 xlim = kwargs.pop('xlim', None)
1570 ylim = kwargs.pop('ylim', None)
1571 calendar_time = kwargs.pop('calendar_time', False)
1572
1573
1574 logx = kwargs.pop('logx', False)
1575 logy = kwargs.pop('logy', False)
1576
1577
1578 loc = kwargs.pop('loc', 'best')
1579
1580
1581 hidden_colorbar = kwargs.pop('hidden_colorbar', False)
1582
1583
1584 bbox_inches = kwargs.pop('bbox_inches', 'tight')
1585
1586
1587 plot = DataPlot(xlabel, ylabel, title, subtitle)
1588
1589
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
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
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
1622 plot.finalize(logx=logx, logy=logy, loc=loc, hidden_colorbar=hidden_colorbar)
1623
1624
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
1643 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True)
1644 if ylim:
1645 plot.ax.set_ylim(tuple(ylim))
1646
1647
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
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
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
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
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
1761 if seglist==None:
1762 seglist = segments.segmentlist()
1763 else:
1764 seglist = segments.segmentlist(seglist)
1765
1766
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
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
1805 vetoLivetime = livetime-float(abs(seglist))
1806
1807
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
1822 if flag:
1823 flag = flag.replace('_','\_')
1824 else:
1825 flag = 'Unknown'
1826
1827
1828 set_rcParams()
1829
1830
1831 xlim = kwargs.pop('xlim', None)
1832 ylim = kwargs.pop('ylim', None)
1833
1834
1835 logx = kwargs.pop('logx', False)
1836 logy = kwargs.pop('logy', False)
1837
1838
1839 hidden_colorbar = kwargs.pop('hidden_colorbar', False)
1840
1841
1842 bbox_inches = kwargs.pop('bbox_inches', 'tight')
1843
1844
1845 fill = kwargs.pop('fill', False)
1846
1847
1848 cumulative = kwargs.pop('cumulative', False)
1849 rate = kwargs.pop('rate', False)
1850
1851
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
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
1886 plot = LineHistogram(xlabel, ylabel, title, subtitle, zlabel)
1887
1888
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
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
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
1923 if xlim:
1924 plot.ax.set_xlim(xlim)
1925 if ylim:
1926 plot.ax.set_ylim(ylim)
1927
1928
1929 set_ticks(plot.ax)
1930
1931
1932 plot.savefig(outfile, bbox_inches=bbox_inches,\
1933 bbox_extra_artists=plot.ax.texts)
1934
1935
1936
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
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
2038 get_time = []
2039 for t in tables:
2040 get_time.append(def_get_time(t.tableName))
2041
2042
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
2055 if not zero:
2056 zero = start
2057
2058
2059 unit,timestr = time_unit(end-start)
2060
2061
2062 if seglist:
2063 segs = segments.segmentlist(seglist)
2064 if not seglist:
2065 segs = segments.segmentlist()
2066
2067
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
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
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
2089 if seglist:
2090 tdata = get_column(triggers, 'time')
2091
2092
2093 vetoData = []
2094 nvetoData = []
2095 for j,tab in enumerate(tables):
2096 vetoData.append({})
2097 nvetoData.append({})
2098
2099 if seglist:
2100 tdata = get_column(tab, 'time')
2101
2102 for i,col in enumerate(columns):
2103 nvetoData[j][col] = get_column(tab, col).astype(float)
2104
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
2124 whitenedFlag = kwargs.pop('whitened', False)
2125 if zcolumn and whitenedFlag:
2126
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
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
2144 flatenedFlag = kwargs.pop('filter', False)
2145 if zcolumn and flatenedFlag:
2146
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
2158
2159 flatenedFlag = kwargs.pop('flaten', False)
2160 if zcolumn and flatenedFlag:
2161
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
2171 minmaxmedianFlag = kwargs.pop('minmaxmedian', False)
2172 if minmaxmedianFlag:
2173 uniqXvalues = numpy.unique1d(nvetoData[j][xcolumn])
2174
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
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
2199 repackTrig = {}
2200 for yVal in uniqYvalues:
2201 repackTrig[yVal] = {}
2202 for iBin in range(len(tedges)-1):
2203 repackTrig[yVal][iBin] = []
2204
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
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
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
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
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
2273 loudest = {}
2274 if len(columns)==4 and\
2275 len(nvetoData[0][columns[0]])+len(vetoData[0][columns[0]])>=1:
2276
2277 vetomax = 0
2278 if len(vetoData[0][columns[3]])>=1:
2279 vetomax = vetoData[0][columns[3]].max()
2280 nvetomax = 0
2281
2282 if len(nvetoData[0][columns[3]])>=1:
2283 nvetomax = nvetoData[0][columns[3]].max()
2284 if vetomax == nvetomax == 0:
2285 pass
2286
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
2297 if flag:
2298 flag = flag.replace('_', '\_')
2299 else:
2300 flag = 'Unknown'
2301
2302
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
2316 set_rcParams()
2317
2318
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
2340 logx = kwargs.pop('logx', False)
2341 logy = kwargs.pop('logy', False)
2342 logz = kwargs.pop('logz', False)
2343
2344
2345 hidden_colorbar = kwargs.pop('hidden_colorbar', False)
2346
2347
2348 bbox_inches = kwargs.pop('bbox_inches', 'tight')
2349
2350
2351 detchar = kwargs.pop('detchar', False)
2352 dcthresh = float(kwargs.pop('dcthreshold', 10))
2353
2354
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
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
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
2373 plot.finalize(logx=logx, logy=logy, hidden_colorbar=hidden_colorbar)
2374
2375 elif len(columns)==4:
2376
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
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
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
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
2402 if loudest:
2403 plot.ax.plot([loudest[columns[0]]], [loudest[columns[1]]],\
2404 marker='*', color='gold', markersize=15)
2405
2406
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
2415 set_ticks(plot.ax, calendar_time=calendar_time)
2416
2417
2418 plot.savefig(outfile, bbox_inches=bbox_inches,\
2419 bbox_extra_artists=plot.ax.texts)
2420 pylab.close(plot.fig)
2421
2422
2423
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
2451 set_rcParams()
2452
2453
2454 xlim = kwargs.pop('xlim', None)
2455 ylim = kwargs.pop('ylim', None)
2456
2457
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
2464 logx = kwargs.pop('logx', False)
2465 logy = kwargs.pop('logy', False)
2466
2467
2468 bbox_inches = kwargs.pop('bbox_inches', 'tight')
2469
2470
2471 hidden_colorbar = kwargs.pop('hidden_colorbar', False)
2472
2473
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
2498 plot = VerticalBarHistogram(xlabel, ylabel, title, subtitle)
2499
2500
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
2506 plot.finalize(num_bins=num_bins, logx=logx, logy=logy,\
2507 hidden_colorbar=hidden_colorbar)
2508
2509
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
2519 plot.savefig(outfile, bbox_inches=bbox_inches,\
2520 bbox_extra_artists=plot.ax.texts)
2521 pylab.close(plot.fig)
2522
2523
2524
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
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
2578 unit, timestr = time_unit(end-start)
2579
2580
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
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
2601 logx = kwargs.pop('logx', False)
2602 logy = kwargs.pop('logy', False)
2603
2604
2605 average = kwargs.pop('average',average)
2606
2607
2608 hidden_colorbar = kwargs.pop('hidden_colorbar', False)
2609
2610
2611 bbox_inches = kwargs.pop('bbox_inches', 'tight')
2612
2613
2614 if not bins:
2615 bins = [[0,float('inf')]]
2616 ybins = [map(float, bin) for bin in bins]
2617
2618
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
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
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
2673 set_rcParams()
2674
2675
2676 plot = ScatterPlot(xlabel, ylabel, title, subtitle)
2677
2678
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
2692 plot.finalize(logx=logx, logy=logy, hidden_colorbar=hidden_colorbar)
2693
2694
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
2701 set_ticks(plot.ax, calendar_time)
2702
2703
2704 plot.savefig(outfile, bbox_inches=bbox_inches,\
2705 bbox_extra_artists=plot.ax.texts)
2706 pylab.close(plot.fig)
2707
2708
2709
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
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
2765 unit, timestr = time_unit(end-start)
2766
2767
2768 if not etg:
2769 if re.search('burst', tableName):
2770 etg = 'Burst'
2771 elif re.search('inspiral', tableName):
2772 etg = 'Inspiral'
2773 elif re.search('ringdown', tableName):
2774 etg = 'Ringdown'
2775 else:
2776 etg = 'Unknown'
2777
2778
2779 calendar_time = kwargs.pop('calendar_time', False)
2780 if calendar_time:
2781 xlim = kwargs.pop('xlim', [gps2datenum(float(start)),\
2782 gps2datenum(float(end))])
2783 else:
2784 xlim = kwargs.pop('xlim', [float(start-zero)/unit, float(end-zero)/unit])
2785 ylim = kwargs.pop('ylim', None)
2786
2787
2788 logx = kwargs.pop('logx', False)
2789 logy = kwargs.pop('logy', False)
2790
2791
2792 average = kwargs.pop('average',average)
2793
2794
2795 hidden_colorbar = kwargs.pop('hidden_colorbar', False)
2796
2797
2798 bbox_inches = kwargs.pop('bbox_inches', 'tight')
2799
2800
2801 if not bins:
2802 bins = [[0,float('inf')]]
2803 ybins = [map(float, bin) for bin in bins]
2804
2805
2806 tbins = {}
2807 rate = {}
2808 rms = {}
2809 for bin in ybins:
2810 if calendar_time:
2811 tbins[bin[0]] = list(gps2datenum(numpy.arange(float(start), float(end),\
2812 average)))
2813 else:
2814 tbins[bin[0]] = list(numpy.arange(0,float(end-start), average)/unit)
2815 rate[bin[0]] = list(numpy.zeros(len(tbins[bin[0]])))
2816 rms[bin[0]] = list(numpy.zeros(len(tbins[bin[0]])))
2817
2818 for trig in triggers:
2819 x = int(float(getTrigAttribute(trig, 'time')-start)//average)
2820 y = getTrigAttribute(trig, bincolumn)
2821 z = getTrigAttribute(trig, rmscolumn)
2822 for bin in ybins:
2823 if bin[0] <= y < bin[1]:
2824 rms[bin[0]][x] += z*z
2825 rate[bin[0]][x] += 1
2826 break
2827
2828
2829 for bin in ybins:
2830 for x in range(len(tbins[bin[0]])):
2831 if rate[bin[0]][x] :
2832 rms[bin[0]][x] = math.sqrt(rms[bin[0]][x]/rate[bin[0]][x])
2833
2834
2835 if logy:
2836 for bin in ybins:
2837 removes = 0
2838 numtbins = len(tbins[bin[0]])
2839 for rbin in xrange(0,numtbins):
2840 if rms[bin[0]][rbin-removes]==0:
2841 rms[bin[0]].pop(rbin-removes)
2842 tbins[bin[0]].pop(rbin-removes)
2843 removes+=1
2844
2845
2846 etg = etg.replace('_', '\_')
2847 zero = LIGOTimeGPS('%.3f' % zero)
2848 if zero.nanoseconds==0:
2849 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero))[:6])\
2850 .strftime("%B %d %Y, %H:%M:%S %ZUTC")
2851 else:
2852 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero.seconds))[:6])\
2853 .strftime("%B %d %Y, %H:%M:%S %ZUTC")
2854 tlabel = tlabel.replace(' UTC', '.%.3s UTC' % zero.nanoseconds)
2855 xlabel = kwargs.pop('xlabel',\
2856 'Time (%s) since %s (%s)' % (timestr, tlabel, zero))
2857 ylabel = kwargs.pop('ylabel', 'RMS')
2858 title = kwargs.pop('title', '%s triggers binned by %s'\
2859 % (etg, display_name(bincolumn)))
2860 if start and end:
2861 subtitle = '%s-%s' % (start, end)
2862 else:
2863 subtitle = " "
2864 subtitle = kwargs.pop('subtitle', subtitle)
2865
2866
2867 set_rcParams()
2868
2869
2870 plot = ScatterPlot(xlabel, ylabel, title, subtitle)
2871
2872
2873 for bin in ybins:
2874 if logy:
2875 if len(rms[bin[0]])>0:
2876 plot.add_content(tbins[bin[0]], rms[bin[0]],\
2877 label='-'.join(map(str, bin)), **kwargs)
2878 else:
2879 plot.add_content([1],[0.1], label='-'.join(map(str, bin)),\
2880 visible=False)
2881 else:
2882 plot.add_content(tbins[bin[0]], rms[bin[0]], label='-'.join(map(str, bin)),\
2883 **kwargs)
2884
2885
2886 plot.finalize(logx=logx, logy=logy, hidden_colorbar=hidden_colorbar)
2887
2888
2889 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True)
2890 plot.ax.set_xlim(xlim)
2891 if ylim:
2892 plot.ax.set_ylim(ylim)
2893
2894
2895 set_ticks(plot.ax, calendar_time=calendar_time)
2896
2897
2898 plot.savefig(outfile, bbox_inches=bbox_inches,\
2899 bbox_extra_artists=plot.ax.texts)
2900 pylab.close(plot.fig)
2901
2902
2903
2904
2905 -def plot_segments(segdict, outfile, start=None, end=None, zero=None,
2906 keys=None, highlight_segments=None, **kwargs):
2907
2908 """
2909 Plot the segments contained within the glue.segments.segmentlistdict
2910 segdict to the given path string outfile. The list keys can be given to
2911 guarantee the order of the segments on the y-axis. x-axis limits can be
2912 controlled using start, end and zero. The glue.segments.segmentlist object
2913 highlight_segments can be given to highlight a number of segments.
2914
2915 Arguments:
2916
2917 segdict : glue.segments.segmentlistdict
2918
2919 """
2920
2921 if not start:
2922 minstart = lambda seglist: min(segdict[key][0][0] for key in segdict.keys()\
2923 if segdict[key])
2924 start = min(minstart(segdict[key]) for key in segdict.keys())
2925 if not end:
2926 maxend = lambda seglist: max(segdict[key][0][1] for key in segdict.keys()\
2927 if segdict[key])
2928 end = max(maxend(segdict[key]) for key in segdict.keys())
2929 if not zero:
2930 zero = start
2931
2932
2933 unit, timestr = time_unit(end-start)
2934
2935
2936 zero = LIGOTimeGPS('%.3f' % zero)
2937 if zero.nanoseconds==0:
2938 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero))[:6])\
2939 .strftime("%B %d %Y, %H:%M:%S %ZUTC")
2940 else:
2941 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero.seconds))[:6])\
2942 .strftime("%B %d %Y, %H:%M:%S %ZUTC")
2943 tlabel = tlabel.replace(' UTC', '.%.3s UTC' % zero.nanoseconds)
2944 xlabel = kwargs.pop('xlabel',\
2945 'Time (%s) since %s (%s)' % (timestr, tlabel, zero))
2946 ylabel = kwargs.pop('ylabel', "")
2947 title = kwargs.pop('title', '')
2948 subtitle = kwargs.pop('subtitle', '')
2949
2950
2951 calendar_time = kwargs.pop('calendar_time', False)
2952 if calendar_time:
2953 xlim = kwargs.pop('xlim', [gps2datenum(float(start)),\
2954 gps2datenum(float(end))])
2955 else:
2956 xlim = kwargs.pop('xlim', [float(start-zero)/unit, float(end-zero)/unit])
2957
2958
2959 labels_inset = kwargs.pop('labels_inset', False)
2960
2961
2962 bbox_inches = kwargs.pop('bbox_inches', 'tight')
2963
2964
2965 hidden_colorbar = kwargs.pop('hidden_colorbar', False)
2966
2967 if keys:
2968
2969 keys = [key.replace('_','\_').replace('\\_','\_') for key in keys]
2970 segkeys = segdict.keys()
2971 newdict = segments.segmentlistdict()
2972 for key in segkeys:
2973 newkey = key.replace('_','\_').replace('\\\_','\_')
2974 newdict[newkey] = segdict[key]
2975 segdict = newdict
2976
2977
2978 set_rcParams()
2979
2980 plot = PlotSegmentsPlot(xlabel, ylabel, title, subtitle, t0=zero, unit=unit,\
2981 calendar_time=calendar_time)
2982 plot.add_content(segdict, keys, **kwargs)
2983 plot.finalize(labels_inset=labels_inset, hidden_colorbar=hidden_colorbar)
2984
2985
2986 if highlight_segments:
2987 for seg in highlight_segments:
2988 plot.highlight_segment(seg)
2989
2990
2991 plot.ax.set_xlim(xlim)
2992
2993 plot.ax.grid(True,which='major')
2994 plot.ax.grid(True,which='majorminor')
2995
2996 set_ticks(plot.ax, calendar_time)
2997
2998 plot.savefig(outfile, bbox_inches=bbox_inches,\
2999 bbox_extra_artists=plot.ax.texts)
3000 pylab.close(plot.fig)
3001
3007
3008 """
3009 Parse ConfigParser.ConfigParser section for plot parameters. Sections should
3010 be name '[plot xcolumn-ycolumn-zcolumn]' e.g.
3011 '[plot time-peak_frequency-snr]'. Returns a pair of dicts with the
3012 following keys:
3013
3014 columns:
3015
3016 xcolumn : [ string | None ]
3017 column string to plot on x-axis
3018 ycolumn : [ string | None ]
3019 column string to plot on y-axis
3020 zcolumn : [ string | None ]
3021 column string to plot on z-axis
3022
3023 params:
3024
3025 xlim : list
3026 [xmin, xmax] pair for x-axis limits
3027 ylim : list
3028 [ymin, ymax] pair for y-axis limits
3029 zlim : list
3030 [zmin, zmax] pair for z-axis limits
3031 clim : list
3032 [cmin, cmax] pair for colorbar limits
3033 logx : bool
3034 True / False to plot log scale on x-axis
3035 logy : bool
3036 True / False to plot log scale on y-axis
3037 logz : bool
3038 True / False to plot log scale on z-axis
3039 """
3040
3041 columns = {'xcolumn':None, 'ycolumn':None, 'zcolumn':None, 'rankcolumn':None}
3042 params = {}
3043
3044 plot = re.split('[\s-]', section)[1:]
3045 if len(plot)>=1:
3046 columns['xcolumn'] = plot[0]
3047 if len(plot)>=2:
3048 columns['ycolumn'] = plot[1]
3049 if len(plot)>=3:
3050 columns['zcolumn'] = plot[2]
3051 if len(plot)>=4:
3052 columns['rankcolumn'] = plot[3]
3053
3054 limits = ['xlim', 'ylim', 'zlim', 'clim', 'exponents', 'constants',\
3055 'color_bins']
3056 filters = ['poles', 'zeros']
3057 bins = ['bins']
3058 booleans = ['logx', 'logy', 'logz', 'cumulative', 'rate', 'detchar',\
3059 'greyscale', 'zeroindicator', 'normalized', 'include_downtime',\
3060 'calendar_time', 'fill', 'hidden_colorbar']
3061 values = ['dcthresh','amplitude','num_bins','linewidth','average','s']
3062
3063
3064 params = {}
3065 for key,val in cp.items(section):
3066 val = val.rstrip('"').strip('"')
3067 if key in limits:
3068 params[key] = map(float, val.split(','))
3069 elif key in filters:
3070 params[key] = map(complex, val.split(','))
3071 elif key in bins:
3072 params[key] = map(lambda p: map(float, p.split(',')), val.split(';'))
3073 elif key in booleans:
3074 params[key] = cp.getboolean(section, key)
3075 elif key in values:
3076 params[key] = float(val)
3077 else:
3078 params[key] = val
3079
3080 return columns, params
3081
3082
3083
3084
3085 -def plot_color_map(data, outfile, data_limits=None, x_format='time',\
3086 y_format='frequency', z_format='amplitude', zero=None,\
3087 x_range=None, y_range=None, **kwargs):
3088
3089 """
3090 Plots data in a 2d color map.
3091 """
3092
3093 if not (isinstance(data, list) or isinstance(data, tuple)):
3094 data_sets = [data]
3095 data_limit_sets = [data_limits]
3096 else:
3097 data_sets = data
3098 data_limit_sets = data_limits
3099 if not len(data_sets) == len(data_limit_sets):
3100 raise AttributeError("You have given %d data sets and %d limit sets! Eh?"\
3101 % (len(data_sets), len(data_limit_sets)))
3102
3103
3104 for i,data_limits in enumerate(data_limit_sets):
3105 if not data_limits:
3106 numrows, numcols = numpy.shape(data_sets[i])
3107 if kwargs.has_key('xlim'):
3108 xlim = kwargs['xlim']
3109 else:
3110 xlim = [0, numrows]
3111 if kwargs.has_key('ylim'):
3112 ylim = kwargs['ylim']
3113 elif pylab.rcParams['image.origin'] == 'upper':
3114 ylim = [numrows, 0]
3115 else:
3116 ylim = [0, numrows]
3117
3118 data_limit_sets[i] = [xlim[0], xlim[1], ylim[0], ylim[1]]
3119
3120
3121 xlim = kwargs.pop('xlim', None)
3122 ylim = kwargs.pop('ylim', None)
3123 zlim = kwargs.pop('zlim', None)
3124 clim = kwargs.pop('clim', None)
3125 calendar_time = kwargs.pop('calendar_time', False)
3126
3127
3128 logx = kwargs.pop('logx', False)
3129 logy = kwargs.pop('logy', False)
3130 logz = kwargs.pop('logz', False)
3131
3132
3133 bbox_inches = kwargs.pop('bbox_inches', 'tight')
3134
3135
3136 for i,(data, data_limits) in enumerate(zip(data_sets, data_limit_sets)):
3137 if logx and xlim:
3138 shape = numpy.shape(data)
3139 xmin, xmax = data_limits[0], data_limits[1]
3140 if x_range==None:
3141 x_range = numpy.logspace(numpy.log10(xmin), numpy.log10(xmax),\
3142 num.shape[-1])
3143 condition = (x_range>xlim[0]) & (x_range<=xlim[1])
3144 newx = numpy.where(condition)[0]
3145 data2 = numpy.resize(data, (len(newx), shape[-2]))
3146 for j in xrange(newx):
3147 data2[:,j] = data[newx[j],:]
3148 data = data2.transpose()
3149
3150 if logy and ylim:
3151 shape = numpy.shape(data)
3152 ymin, ymax = data_limits[2], data_limits[3]
3153 if y_range==None:
3154 y_range = numpy.logspace(numpy.log10(ymin), numpy.log10(ymax),\
3155 num=shape[-2])
3156 condition = (y_range>ylim[0]) & (y_range<=ylim[1])
3157 newy = numpy.where((condition))[0]
3158 data2 = numpy.resize(data, (len(newy), shape[-1]))
3159 for j in xrange(shape[-1]):
3160 data2[:,j] = data[:,j][condition]
3161 data = data2
3162 data_limits[2:] = ylim
3163 data_sets[i] = data
3164
3165
3166 columns = list(map(str.lower, [x_format, y_format, z_format]))
3167 limits = [xlim, ylim, zlim]
3168 for i,lim in enumerate(limits):
3169 if lim:
3170 limits[i] = list(map(float, lim))
3171 labels = ["", "", "", "", ""]
3172
3173
3174 if 'time' in columns and not zero:
3175 i = columns.index('time')
3176 if limits[i]:
3177 zero = limits[i][0]
3178 else:
3179 zero = min(data_limits[0] for data_limits in data_limit_sets)
3180
3181
3182 for i,(col,c) in enumerate(zip(columns, ['x', 'y', 'z'])):
3183 if col.lower() == 'time' and limits[i]:
3184 unit, timestr = time_unit(float(limits[i][1]-limits[i][0]))
3185 zero = LIGOTimeGPS('%.3f' % zero)
3186 if zero.nanoseconds==0:
3187 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero))[:6])\
3188 .strftime("%B %d %Y, %H:%M:%S %ZUTC")
3189 else:
3190 tlabel = datetime(*date.XLALGPSToUTC(LIGOTimeGPS(zero.seconds))[:6])\
3191 .strftime("%B %d %Y, %H:%M:%S %ZUTC")
3192 tlabel = tlabel.replace(' UTC', '.%.3s UTC' % zero.nanoseconds)
3193 labels[i] = kwargs.pop('%slabel' % c, 'Time (%s) since %s (%s)'\
3194 % (timestr, tlabel, zero))
3195 if calendar_time:
3196 limits[i] = gps2datenum(numpy.asarray(limits[i]))
3197 else:
3198 limits[i] = (numpy.asarray(limits[i])-float(zero))/unit
3199 for i,data_limits in enumerate(data_limit_sets):
3200 if calendar_time:
3201 data_limit_sets[i][0] = gps2datenum(float(data_limits[0]))
3202 data_limit_sets[i][1] = gps2datenum(float(data_limits[1]))
3203 else:
3204 data_limit_sets[i][0] = float(data_limits[0]-zero)/unit
3205 data_limit_sets[i][1] = float(data_limits[1]-zero)/unit
3206 else:
3207 labels[i] = kwargs.pop('%slabel' % c, display_name(col))
3208
3209 labels[3] = kwargs.pop('title', "")
3210 labels[4] = kwargs.pop('subtitle', "")
3211
3212
3213 set_rcParams()
3214
3215
3216 plot = ColorMap(*labels)
3217
3218
3219 for i, (data, data_limits) in enumerate(zip(data_sets, data_limit_sets)):
3220 plot.add_content(data, data_limits, aspect='auto')
3221
3222
3223 plot.finalize(logx=logx, logy=logy, logz=logz, clim=clim, origin='lower')
3224
3225 if limits[0] is not None and len(limits[0])==2:
3226 plot.ax.set_xlim(limits[0])
3227 if limits[1] is not None and len(limits[1])==2:
3228 plot.ax.set_ylim(limits[1])
3229
3230
3231 set_ticks(plot.ax, calendar_time=calendar_time)
3232
3233
3234 plot.savefig(outfile, bbox_inches=bbox_inches,\
3235 bbox_extra_artists=plot.ax.texts)
3236 pylab.close(plot.fig)
3237
3242
3243 """
3244 Plot significance drop for each channel relative to the application of
3245 HVeto round veto segments.
3246 """
3247
3248
3249 channels = startsig.keys()
3250 for c in channels:
3251 if c not in endsig.keys():
3252 raise AttributeError("Significance lists do not match.")
3253 channels.sort()
3254
3255
3256 wch,wsig = max(startsig.items(), key=lambda x: x[1])
3257
3258
3259 kwargs.pop('xlim', None)
3260 ylim = kwargs.pop('ylim', (0, wsig+1))
3261 xlabel = kwargs.pop('xlabel', "")
3262 ylabel = kwargs.pop('ylabel', "Significance")
3263 title = kwargs.pop('title', "Coincidence significance drop plot")
3264 subtitle = kwargs.pop('subtitle',\
3265 "Winner: %s, significance: %s" % (wch, wsig))
3266
3267 kwargs.setdefault('linestyle', '-')
3268 kwargs.setdefault('marker', 'o')
3269 color = kwargs.pop('color', None)
3270
3271
3272 set_rcParams()
3273 pylab.rcParams.update({"figure.figsize":[24,6], "xtick.labelsize": 8})
3274
3275
3276 hidden_colorbar = kwargs.pop('hidden_colorbar', False)
3277
3278
3279 bbox_inches = kwargs.pop('bbox_inches', 'tight')
3280
3281
3282 plot = DataPlot(xlabel, ylabel, title, subtitle)
3283
3284 kwargs.setdefault('markersize', 10)
3285
3286
3287 for i,c in enumerate(channels):
3288 s = startsig[c]
3289 e = endsig[c]
3290 col = color and color or s>e and 'b' or 'r'
3291 plot.add_content([i,i], [s,e], color=col, **kwargs)
3292
3293
3294 plot.finalize(logx=False, logy=False, hidden_colorbar=hidden_colorbar)
3295
3296
3297 plot.ax.set_xlim(-1, len(channels))
3298 plot.ax.set_xticks(numpy.arange(0,len(channels)))
3299 plot.ax.set_xticklabels([c.replace('_','\_') for c in channels])
3300 for i,t in enumerate(plot.ax.get_xticklabels()):
3301 t.set_rotation(315)
3302
3303 t.set_verticalalignment('top')
3304 t.set_horizontalalignment('left')
3305
3306
3307
3308 plot.ax.set_ylim(ylim)
3309
3310
3311 plot.ax.xaxis.grid(False)
3312
3313 plot.fig.subplots_adjust(bottom=0.3)
3314
3315
3316 plot.savefig(outfile, bbox_inches=bbox_inches,\
3317 bbox_extra_artists=plot.ax.texts)
3318 pylab.close(plot.fig)
3319
3320
3321
3322
3323 -def plot_sky_positions(skyTable, outfile, format='radians', zcolumn=None,\
3324 projection=None, centre=None, detectors=[], lines={},\
3325 range=[(0,0), (1,1)], **kwargs):
3326
3327 """
3328 Plot latitude against longitude for the given ligolw table skyTable into the
3329 given outfile. Uses the mpl_toolkits basemap module to plot the sky sphere
3330 in a variety of projections, or simply a scatter plot if projection=None.
3331 Can plot lines and detector positions on top if given.
3332
3333 Arguments:
3334
3335 skyTable : glue.ligolw.table.Table
3336 ligolw table containing triggers or SkyPosition objects
3337 outfile : string
3338 string path for output plot
3339
3340 Keyword arguments:
3341
3342 zcolumn : string
3343 valid column of ligolw table to use for colorbar (optional).
3344 format : [ 'radians' | 'degrees' ]
3345 str identifying format of longtiude/ra, latitude/dec colums in table.
3346 Plot is always drawn in degrees.
3347 projection : str
3348 type of spherical projection to use, if any. See matplotlib Basemap
3349 documentation for details, recommended: 'ortho', 'hammer'.
3350 centre : tuple
3351 (longitude, latitude) pair on which to centre plot. Latitude centring
3352 only works for certain projections.
3353 detectors : list
3354 list of detector prefixes to plot, e.g. ['H1', 'L1'].
3355 lines : dict
3356 dict of name:table pairs from which to plot lines (+ marker) on top of
3357 points
3358 range : tuple
3359 ((xmin, ymin), (xmax, ymax)) tuples for lower left and upper right
3360 corners, in range 0-1, e.g. (0,0) is full lower left, (1,1) full upper
3361 right corners for full spherical projection. Only works with certain
3362 projections.
3363
3364 Unnamed keyword arguments:
3365
3366 logx : [ True | False ]
3367 boolean option to display x-axis in log scale. Not applicable when using
3368 projection.
3369 logy : [ True | False ]
3370 boolean option to display y-axis in log scale. Not applicable when using
3371 projection.
3372 logz : [ True | False ]
3373 boolean option to display z-axis in log scale.
3374 xlim : tuple
3375 (xmin, xmax) limits for x-axis (longitude). Triggers outside range are
3376 removed.
3377 ylim : tuple
3378 (ymin, ymax) limits for y-axis (latitude). Triggers outside range are
3379 removed.
3380 zlim : tuple
3381 (zmin, zmax) limits for z-axis. Triggers outside range are removed.
3382 clim : tuple
3383 (cmin, cmax) limits for color scale. Triggers outside range are moved
3384 onto boundary.
3385 xlabel : string
3386 label for x-axis
3387 ylabel : string
3388 label for y-axis
3389 zlabel : string
3390 label for z-axis
3391 title : string
3392 title for plot
3393 subtitle : string
3394 subtitle for plot
3395
3396 All other given arguments will be passed to matplotlib.axes.Axes.scatter.
3397 """
3398
3399 format = format.lower()
3400
3401
3402 if projection:
3403 xlabel = kwargs.pop("xlabel", "")
3404 ylabel = kwargs.pop("ylabel", "")
3405 else:
3406 xlabel = kwargs.pop("xlabel", "Longitude (degrees of arc)")
3407 ylabel = kwargs.pop("ylabel", "Latitude (degrees of arc)")
3408 if zcolumn: tmp = zcolumn
3409 else: tmp = ""
3410 zlabel = kwargs.pop("zlabel", display_name(tmp))
3411
3412
3413 if re.search('burst', skyTable.tableName.lower()):
3414 etg = 'Burst'
3415 elif re.search('inspiral', skyTable.tableName.lower()):
3416 etg = 'Inspiral'
3417 elif re.search('ringdown', skyTable.tableName.lower()):
3418 etg = 'Ringdown'
3419 else:
3420 etg = 'Unknown'
3421 etg = etg.replace('_', '\_')
3422
3423
3424 title = kwargs.pop("title", "")
3425 subtitle = kwargs.pop("subtitle", "")
3426
3427
3428 xlim = kwargs.pop('xlim', None)
3429 ylim = kwargs.pop('ylim', None)
3430 zlim = kwargs.pop('zlim', None)
3431 clim = kwargs.pop('clim', None)
3432
3433
3434 logx = kwargs.pop('logx', False)
3435 logy = kwargs.pop('logy', False)
3436 logz = kwargs.pop('logz', False)
3437
3438
3439 bbox_inches = kwargs.pop('bbox_inches', 'tight')
3440
3441
3442 detData = []
3443 detCol = []
3444 for ifo in detectors:
3445 d = XLALTools.cached_detector[inject.prefix_to_name[ifo]]
3446 if format=='radians':
3447 detData.append((numpy.degrees(d.vertexLongitudeRadians),\
3448 numpy.degrees(d.vertexLatitudeRadians)))
3449 else:
3450 detData.append((d.vertexLongitudeRadians, d.vertexLatitudeRadians))
3451 if ifo in plotsegments.PlotSegmentsPlot.color_code.keys():
3452 detCol.append(plotsegments.PlotSegmentsPlot.color_code[ifo])
3453 else:
3454 detCol.append('cyan')
3455
3456
3457
3458 lineData = []
3459 lineCol = []
3460 for line in lines:
3461 if isinstance(lines, dict):
3462 label = line
3463 line = lines[label]
3464 else:
3465 label = '_'
3466
3467 if re.search('multi_inspiral', line.tableName):
3468 loncol = 'ra'
3469 latcol = 'dec'
3470 else:
3471 loncol = 'longitude'
3472 latcol = 'latitude'
3473
3474 lon = []
3475 lat = []
3476 for p in line:
3477 if format=='radians':
3478 lon.append(numpy.degrees(getattr(p, loncol)))
3479 lat.append(numpy.degrees(getattr(p, latcol)))
3480 else:
3481 lon.append(getattr(p, loncol))
3482 lat.append(getattr(p, latcol))
3483 lineData.append((label, [lon, lat]))
3484 if len(lineData):
3485 for c in plotutils.default_colors():
3486 if c=='b': continue
3487 lineCol.append(c)
3488 if len(lineCol)==len(lineData): break
3489
3490
3491 columns = ['longitude', 'latitude']
3492 if re.search('multi_inspiral', skyTable.tableName):
3493 columns[0] = 'ra'
3494 columns[1] = 'dec'
3495 if zcolumn:
3496 columns.append(zcolumn)
3497
3498 data = {}
3499 for col in columns:
3500 data[col] = []
3501 for i,p in enumerate(skyTable):
3502
3503 if xlim and not xlim[0] <= getTrigAttribute(p, columns[0]) <= xlim[1]:
3504 continue
3505 if ylim and not ylim[0] <= getTrigAttribute(p, columns[1]) <= ylim[1]:
3506 continue
3507 if zlim and not zlim[0] <= getTrigAttribute(p, columns[2]) <= zlim[1]:
3508 continue
3509
3510 for col in columns[0:2]:
3511 if format=='radians':
3512 data[col].append(numpy.degrees(getTrigAttribute(p, col)))
3513 else:
3514 data[col].append(getTrigAttribute(p, col))
3515
3516 if zcolumn:
3517 data[columns[2]].append(getTrigAttribute(p, columns[2]))
3518
3519
3520 set_rcParams()
3521 if projection not in [None, 'hammer', 'moll', 'robin', 'sinu', 'cyl', 'merc',\
3522 'mill', 'gall']:
3523 pylab.rcParams.update({"figure.figsize":[8,6]})
3524
3525
3526 if projection:
3527
3528 if zcolumn:
3529 plot = ColorbarSkyPositionsPlot(xlabel, ylabel, zlabel,\
3530 title, subtitle)
3531 plot.add_content(data[columns[0]], data[columns[1]], data[zcolumn],\
3532 **kwargs)
3533 for p,ifo,c in zip(detData, detectors, detCol):
3534 lon, lat = p
3535 plot.add_content([lon], [lat], [1], label=ifo, facecolor=c,\
3536 edgecolor=c, **kwargs)
3537 for line,c in zip(lineData, lineCol):
3538 label, data = line
3539 plot.add_content(data[0], data[1], label=label, marker='+',\
3540 edgecolor=c, s=4)
3541 plot.finalize(projection=projection, centre=centre, range=range,\
3542 logz=logz, clim=clim)
3543
3544 else:
3545 plot = SkyPositionsPlot(xlabel, ylabel, title, subtitle)
3546 plot.add_content(data[columns[0]], data[columns[1]], label='_', **kwargs)
3547 for p,ifo,c in zip(detData, detectors, detCol):
3548 lon, lat = p
3549 plot.add_content([lon], [lat], label=ifo, facecolor=c, edgecolor=c,\
3550 **kwargs)
3551 for line,c in zip(lineData, lineCol):
3552 label, data = line
3553 plot.add_content(data[0], data[1], label=label, marker='+',\
3554 edgecolor=c, s=4)
3555 plot.finalize(projection=projection, centre=centre, range=range)
3556
3557
3558 else:
3559 if zcolumn:
3560 plot = ColorbarScatterPlot(xlabel, ylabel, zlabel, title,\
3561 subtitle)
3562 plot.add_content(data[columns[0]], data[columns[1]], data[columns[2]],\
3563 **kwargs)
3564 plot.finalize(logx=logx, logy=logy, logz=logz, clim=clim)
3565 else:
3566 plot = ScatterPlot(xlabel, ylabel, title, subtitle)
3567 plot.add_content(data[columns[0]], data[columns[1]], **kwargs)
3568 plot.finalize(logx=logx, logy=logy)
3569
3570
3571 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True)
3572 plot.ax.set_xlim(xlim)
3573 plot.ax.set_ylim(ylim)
3574
3575
3576 plot.savefig(outfile, bbox_inches=bbox_inches,\
3577 bbox_extra_artists=plot.ax.texts)
3578 pylab.close(plot.fig)
3579