pylal.plotdata

1 #!/usr/bin/env python 2 3 # Copyright (C) 2011 Duncan Macleod 4 # 5 # This program is free software; you can redistribute it and/or modify it 6 # under the terms of the GNU General Public License as published by the 7 # Free Software Foundation; either version 3 of the License, or (at your 8 # option) any later version. 9 # 10 # This program is distributed in the hope that it will be useful, but 11 # WITHOUT ANY WARRANTY; without even the implied warranty of 12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 13 # Public License for more details. 14 # 15 # You should have received a copy of the GNU General Public License along 16 # with this program; if not, write to the Free Software Foundation, Inc., 17 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 18 19 """ 20 This module provides some wrapper functions for plotting TimeSeries and FrequencySeries using classes defined in pylal.plotutils. 21 """ 22 23 # ============================================================================= 24 # Preamble 25 # ============================================================================= 26 27 from __future__ import division 28 29 import datetime 30 import re 31 import itertools 32 import numpy 33 import numbers 34 35 from scipy import interpolate 36 37 import lal 38 39 from pylal import plotutils 40 from pylal import seriesutils 41 from pylal import git_version 42 from lal import rate 43 44 __author__ = "Duncan M. Macleod <duncan.macleod@ligo.org>" 45 __version__ = git_version.id 46 __date__ = git_version.date 47 48 # ============================================================================= 49 # Plot TimeSeries 50 # ============================================================================= 51

52 -def plottimeseries(series, outfile, t0=0, zeroline=False, **kwargs):

53 """ 54 Plot a REALXTimeSeries. 55 """ 56 # construct list of series 57 if hasattr(series, "__contains__"): 58 serieslist = series 59 else: 60 serieslist = [series] 61 62 # get limits 63 xlim = kwargs.pop("xlim", None) 64 ylim = kwargs.pop("ylim", None) 65 if xlim: 66 start,end = xlim 67 else: 68 start = min(float(s.epoch) for s in serieslist) 69 end = max(float(s.epoch + s.data.length*s.deltaT) for s in serieslist) 70 71 # get axis scales 72 logx = kwargs.pop("logx", False) 73 logy = kwargs.pop("logy", False) 74 75 # get legend loc 76 loc = kwargs.pop("loc", 0) 77 alpha = kwargs.pop("alpha", 0.8) 78 79 # get colorbar options 80 hidden_colorbar = kwargs.pop("hidden_colorbar", False) 81 82 # get savefig option 83 bbox_inches = kwargs.pop("bbox_inches", None) 84 85 # 86 # get labels 87 # 88 89 xlabel = kwargs.pop("xlabel", None) 90 if xlabel: 91 unit = 1 92 if not xlabel: 93 unit, timestr = plotutils.time_axis_unit(end-start) 94 if not t0: 95 t0 = start 96 t0 = lal.LIGOTimeGPS(t0) 97 if int(t0.gpsNanoSeconds) == 0: 98 xlabel = datetime.datetime(*lal.GPSToUTC(int(t0))[:6])\ 99 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 100 xlabel = "Time (%s) since %s (%s)" % (timestr, xlabel, int(t0)) 101 else: 102 xlabel = datetime.datetime(*lal.GPSToUTC(t0.gpsSeconds)[:6])\ 103 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 104 xlabel = "Time (%s) since %s (%s)"\ 105 % (timestr, 106 xlabel.replace(" UTC",".%.3s UTC" % t0.gpsNanoSeconds),\ 107 t0) 108 t0 = float(t0) 109 ylabel = kwargs.pop("ylabel", "Amplitude (counts)") 110 title = kwargs.pop("title", "") 111 subtitle = kwargs.pop("subtitle", "") 112 113 # 114 # set default line params 115 # 116 117 color = kwargs.pop("color", None) 118 if isinstance(color, str): 119 color = color.split(',') 120 if len(color) == 1: 121 color = [color[0]]*len(serieslist) 122 kwargs2 = dict() 123 kwargs2.update(kwargs) 124 if kwargs.has_key("marker"): 125 kwargs.setdefault("markersize", 5) 126 kwargs2.setdefault("markersize", 2) 127 else: 128 kwargs.setdefault("linestyle", "-") 129 kwargs.setdefault("linewidth", "1") 130 kwargs2.setdefault("linewidth", "0.1") 131 132 # 133 # make plot 134 # 135 136 allnames = [s.name for s in serieslist] 137 namedseries = [s for s in serieslist if not re.search("(min|max)\Z",s.name)] 138 139 plot = plotutils.SimplePlot(xlabel, ylabel, title, subtitle) 140 for i,(series,c) in enumerate(itertools.izip(namedseries,\ 141 plotutils.default_colors())): 142 if color: 143 c = color[serieslist.index(series)] 144 x = numpy.arange(series.data.length) * series.deltaT +\ 145 float(series.epoch) - float(t0) 146 x = x.astype(float) 147 x /= unit 148 d = series.data.data 149 if logy and ylim: 150 numpy.putmask(d, d==0, ylim[0]-abs(ylim[0])*0.01) 151 plot.add_content(x, d, color=c,\ 152 label=plotutils.display_name(series.name), **kwargs) 153 # find min/max and plot 154 for i,name in enumerate(allnames): 155 for ext in ["min", "max"]: 156 if re.match("%s[- _]%s" % (re.escape(series.name), ext), name): 157 if color: 158 c = color[i] 159 series2 = serieslist[i] 160 x2 = numpy.arange(series2.data.length) * series2.deltaT\ 161 + float(series2.epoch) - float(t0) 162 x2 /= unit 163 d2 = series2.data.data 164 if logy and ylim: 165 numpy.putmask(d2, d2==0, ylim[0]-abs(ylim[0])*0.01) 166 plot.ax.plot(x2, d2, color=c, **kwargs2) 167 plot.ax.fill_between(x2, d, d2, alpha=0.1,\ 168 color=c) 169 170 # finalize 171 plot.finalize(loc=loc, alpha=alpha) 172 if hidden_colorbar: 173 plotutils.add_colorbar(plot.ax, visible=False) 174 175 # add zero line 176 axis_lims = plot.ax.get_ylim() 177 if zeroline: 178 plot.ax.plot([0, 0], [axis_lims[0], axis_lims[1]], 'r--', linewidth=2) 179 plot.ax.set_ylim([ axis_lims[0], axis_lims[1] ]) 180 181 # set logscale 182 if logx: 183 plot.ax.xaxis.set_scale("log") 184 if logy: 185 plot.ax.yaxis.set_scale("log") 186 plot.ax._update_transScale() 187 188 # format axes 189 if xlim: 190 xlim = (numpy.asarray(xlim).astype(float)-t0)/unit 191 plot.ax.set_xlim(xlim) 192 if ylim: 193 plot.ax.set_ylim(ylim) 194 195 plotutils.set_time_ticks(plot.ax) 196 plotutils.set_minor_ticks(plot.ax, x=False) 197 198 # save and close 199 plot.savefig(outfile, bbox_inches=bbox_inches,\ 200 bbox_extra_artists=plot.ax.texts) 201 plot.close()

202 203 # ============================================================================= 204 # Plot FrequencySeries 205 # ============================================================================= 206

207 -def plotfrequencyseries(series, outfile, **kwargs):

208 """ 209 Plot a (swig)LAL FrequencySeries. 210 """ 211 # construct list of series 212 if hasattr(series, "__contains__"): 213 serieslist = series 214 else: 215 serieslist = [series] 216 217 # get limits 218 xlim = kwargs.pop("xlim", None) 219 ylim = kwargs.pop("ylim", None) 220 if not xlim: 221 fmin = min(float(s.f0) for s in serieslist) 222 fmax = max(float(s.f0 + s.data.length*s.deltaF) for s in serieslist) 223 xlim = [fmin, fmax] 224 225 # get axis scales 226 logx = kwargs.pop("logx", False) 227 logy = kwargs.pop("logy", False) 228 229 # get legend loc 230 loc = kwargs.pop("loc", 0) 231 alpha = kwargs.pop("alpha", 0.8) 232 233 # get colorbar options 234 hidden_colorbar = kwargs.pop("hidden_colorbar", False) 235 236 # get savefig option 237 bbox_inches = kwargs.pop("bbox_inches", None) 238 239 # 240 # get labels 241 # 242 243 xlabel = kwargs.pop("xlabel", "Frequency (Hz)") 244 ylabel = kwargs.pop("ylabel", "Amplitude") 245 title = kwargs.pop("title", "") 246 subtitle = kwargs.pop("subtitle", "") 247 248 # 249 # set default line params 250 # 251 252 color = kwargs.pop("color", None) 253 if isinstance(color, str): 254 color = [color]*len(serieslist) 255 kwargs2 = dict() 256 kwargs2.update(kwargs) 257 if kwargs.has_key("marker"): 258 kwargs.setdefault("markersize", 5) 259 kwargs2.setdefault("markersize", 2) 260 else: 261 kwargs.setdefault("linestyle", "-") 262 kwargs.setdefault("linewidth", "1") 263 kwargs2.setdefault("linewidth", "0.1") 264 265 # 266 # make plot 267 # 268 269 allnames = [s.name for s in serieslist] 270 namedseries = [s for s in serieslist if not re.search("(min|max)\Z",s.name)] 271 272 plot = plotutils.SimplePlot(xlabel, ylabel, title, subtitle) 273 for i,(series,c) in enumerate(itertools.izip(namedseries,\ 274 plotutils.default_colors())): 275 if color: 276 c = color[serieslist.index(series)] 277 if series.f_array is not None: 278 x = series.f_array 279 else: 280 x = numpy.arange(series.data.length) * series.deltaF + series.f0 281 x = x.astype(float) 282 plot.add_content(x, series.data.data, color=c, 283 label=plotutils.display_name(series.name), **kwargs) 284 # find min/max and plot 285 for i,name in enumerate(allnames): 286 for ext in ["min", "max"]: 287 if re.match("%s[- _]%s" % (re.escape(series.name), ext), name): 288 if color: 289 c = color[i] 290 series2 = serieslist[i] 291 if series2.f_array is not None: 292 x2 = series2.f_array 293 else: 294 x2 = numpy.arange(series2.data.length) * series2.deltaF\ 295 + series2.f0 296 plot.ax.plot(x2, series2.data.data, color=c, **kwargs2) 297 # sanity check for malformed inputs 298 if series.data.data.shape == series2.data.data.shape: 299 plot.ax.fill_between(x2, series.data.data,\ 300 series2.data.data, alpha=0.1,\ 301 color=c) 302 303 # finalize 304 plot.finalize(loc=loc, alpha=alpha) 305 if hidden_colorbar: 306 plotutils.add_colorbar(plot.ax, visible=False) 307 308 # set logscale 309 if logx: 310 plot.ax.xaxis.set_scale("log") 311 if logy: 312 plot.ax.yaxis.set_scale("log") 313 plot.ax._update_transScale() 314 315 # format axes 316 if xlim: 317 plot.ax.set_xlim(xlim) 318 if ylim: 319 plot.ax.set_ylim(ylim) 320 321 plot.ax.grid(True, which="both") 322 plotutils.set_minor_ticks(plot.ax) 323 324 # save and close 325 plot.savefig(outfile, bbox_inches=bbox_inches,\ 326 bbox_extra_artists=plot.ax.texts) 327 plot.close()

328 329 # ============================================================================= 330 # Plot spectrogram 331 # ============================================================================= 332

333 -def plotspectrogram(sequencelist, outfile, epoch=0, deltaT=1, f0=0, deltaF=1,\ 334 t0=0, ydata=None, **kwargs):

335 """ 336 Plots a list of REAL?VectorSequences on a time-frequency-amplitude colour 337 map. The epochand deltaT arguments define the spacing in the x direction 338 for the given VectorSequence, and similarly f0 and deltaF in the 339 y-direction. If a list of VectorSequences is given, any of these arguments 340 can be in list form, one for each sequence. 341 342 ydata can be given as to explicitly define the frequencies at which the 343 sequences are sampled. 344 """ 345 # construct list of series 346 if not hasattr(sequencelist, "__contains__"): 347 sequencelist = [sequencelist] 348 numseq = len(sequencelist) 349 350 # format variables 351 if isinstance(epoch, numbers.Number) or isinstance(epoch, lal.LIGOTimeGPS): 352 epoch = [epoch]*numseq 353 epoch = map(float, epoch) 354 if not len(epoch) == numseq: 355 raise ValueError("Wrong number of epoch arguments given.") 356 if isinstance(deltaT, numbers.Number): 357 deltaT = [deltaT]*numseq 358 deltaT = map(float, deltaT) 359 if not len(deltaT) == numseq: 360 raise ValueError("Wrong number of deltaT arguments given.") 361 if not ydata is None: 362 if isinstance(f0, numbers.Number): 363 f0 = [f0]*numseq 364 f0 = map(float, f0) 365 if not len(f0) == numseq: 366 raise ValueError("Wrong number of f0 arguments given.") 367 if isinstance(deltaF, numbers.Number): 368 deltaF = [deltaF]*numseq 369 deltaF = map(float, deltaF) 370 if not len(deltaF) == numseq: 371 raise ValueError("Wrong number of deltaF arguments given.") 372 373 # get limits 374 xlim = kwargs.pop("xlim", None) 375 ylim = kwargs.pop("ylim", None) 376 colorlim = kwargs.pop("colorlim", None) 377 if xlim: 378 start,end = xlim 379 else: 380 start = min(epoch) 381 end = max(e + l.length * dt\ 382 for e,l,dt in zip(epoch, sequencelist, deltaT)) 383 if not ydata is None and not ylim: 384 ylim = [ydata.min(), ydata.max()] 385 386 # get axis scales 387 logx = kwargs.pop("logx", False) 388 logy = kwargs.pop("logy", False) 389 logcolor = kwargs.pop("logcolor", False) 390 391 # get legend loc 392 loc = kwargs.pop("loc", 0) 393 alpha = kwargs.pop("alpha", 0.8) 394 395 # get colorbar options 396 hidden_colorbar = kwargs.pop("hidden_colorbar", False) 397 398 # get savefig option 399 bbox_inches = kwargs.pop("bbox_inches", None) 400 401 # 402 # get labels 403 # 404 405 xlabel = kwargs.pop("xlabel", None) 406 if xlabel: 407 unit = 1 408 if not xlabel: 409 unit, timestr = plotutils.time_axis_unit(end-start) 410 if not t0: 411 t0 = start 412 t0 = lal.LIGOTimeGPS(t0) 413 if int(t0.gpsNanoSeconds) == 0: 414 xlabel = datetime.datetime(*lal.GPSToUTC(int(t0))[:6])\ 415 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 416 xlabel = "Time (%s) since %s (%s)" % (timestr, xlabel, int(t0)) 417 else: 418 xlabel = datetime.datetime(*lal.GPSToUTC(t0.gpsSeconds)[:6])\ 419 .strftime("%B %d %Y, %H:%M:%S %ZUTC") 420 xlabel = "Time (%s) since %s (%s)"\ 421 % (timestr, 422 xlabel.replace(" UTC",".%.3s UTC" % t0.gpsNanoSeconds),\ 423 t0) 424 t0 = float(t0) 425 ylabel = kwargs.pop("ylabel", "Frequency (Hz)") 426 colorlabel = kwargs.pop("colorlabel", "Amplitude") 427 title = kwargs.pop("title", "") 428 subtitle = kwargs.pop("subtitle", "") 429 430 # 431 # restrict data to the correct limits for plotting 432 # 433 434 interpolate = logy and ydata is None 435 436 for i,sequence in enumerate(sequencelist): 437 if interpolate: 438 # interpolate the data onto a log-scale 439 sequence,ydata = loginterpolate(sequence, f0[i], deltaF[i]) 440 if logy and ylim: 441 plotted = (ydata > ylim[0]) & (ydata <= ylim[1]) 442 newVectorLength = int(plotted.sum()) 443 newsequence = lal.CreateREAL8VectorSequence(sequence.length,\ 444 newVectorLength) 445 for j in range(sequence.length): 446 newsequence.data[j,:] = sequence.data[j,:][plotted] 447 del sequence 448 sequencelist[i] = newsequence 449 if len(sequencelist) and logy and ylim: 450 ydata = ydata[plotted] 451 452 # 453 # format bins 454 # 455 456 xbins = [] 457 for i in range(numseq): 458 xmin = epoch[i] 459 xmax = epoch[i] + sequencelist[i].length * deltaT[i] 460 xbins.append(rate.LinearBins(float(xmin-t0)/unit, float(xmax-t0)/unit,\ 461 2)) 462 463 ybins = [] 464 for i in range(numseq): 465 if ydata is not None: 466 ydata = numpy.asarray(ydata) 467 ymin = ydata.min() 468 ymax = ydata.max() 469 else: 470 ymin = f0[i] 471 ymax = f0[i] + sequencelist[i].vectorLength * deltaF[i] 472 if logy: 473 if ymin == 0: 474 ymin = deltaF[i] 475 ybins.append(rate.LogarithmicBins(ymin, ymax, 2)) 476 else: 477 ybins.append(rate.LinearBins(ymin, ymax, 2)) 478 479 # 480 # plot 481 # 482 483 kwargs.setdefault("interpolation", "kaiser") 484 485 plot = plotutils.ImagePlot(xlabel=xlabel, ylabel=ylabel, title=title,\ 486 subtitle=subtitle, colorlabel=colorlabel) 487 488 for sequence,x,y in zip(sequencelist, xbins, ybins): 489 data = numpy.ma.masked_where(numpy.isnan(sequence.data), sequence.data,\ 490 copy=False) 491 plot.add_content(data.T, x, y, **kwargs) 492 493 # finalize 494 plot.finalize(colorbar=True, logcolor=logcolor, minorticks=True,\ 495 clim=colorlim) 496 if hidden_colorbar: 497 plotutils.add_colorbar(plot.ax, visible=False) 498 499 # set logscale 500 if logx: 501 plot.ax.xaxis.set_scale("log") 502 if logy: 503 plot.ax.yaxis.set_scale("log") 504 plot.ax._update_transScale() 505 506 # format axes 507 if xlim: 508 xlim = (numpy.asarray(xlim).astype(float)-t0)/unit 509 plot.ax.set_xlim(xlim) 510 if ylim: 511 plot.ax.set_ylim(ylim) 512 513 # set grid and ticks 514 plot.ax.grid(True, which="both") 515 plotutils.set_time_ticks(plot.ax) 516 plotutils.set_minor_ticks(plot.ax) 517 518 # save and close 519 plot.savefig(outfile, bbox_inches=bbox_inches,\ 520 bbox_extra_artists=plot.ax.texts) 521 plot.close()

522

523 -def loginterpolate(sequence, y0, deltaY, N=None):

524 """ 525 Interpolate a REAL?VectorSequence into logarithmic spacing in the second 526 dimension (y-axis, or frequency). 527 """ 528 # work out y-axis limits and new y-range 529 if not N: 530 N = sequence.vectorLength 531 ylin = numpy.arange(sequence.vectorLength)*deltaY + y0 532 ymin = y0 and y0 or deltaY 533 ymax = y0 + sequence.vectorLength * deltaY 534 ylog = numpy.logspace(numpy.log10(ymin*1.01), numpy.log10(ymax*0.99), num=N,\ 535 endpoint=False) 536 537 # make new sequence 538 datatype = seriesutils.typecode(type(sequence)) 539 TYPESTR = seriesutils._typestr[datatype] 540 func = getattr(lal, "Create%sVectorSequence" % TYPESTR) 541 out = func(sequence.length, N) 542 543 # interpolate columnwise 544 for i in range(sequence.length): 545 intplt = interpolate.interp1d(ylin, sequence.data[i,:]) 546 out.data[i,:] = intplt(ylog) 547 548 return out,ylog

549 550 # ============================================================================= 551 # Plot histogram of series 552 # ============================================================================= 553

554 -def plothistogram(serieslist, outfile, nonzero=False, num_bins=100,\ 555 cumulative=False, bar=False, **kwargs):

556 """ 557 Plots a line histogram of the data contained in the series, or list of 558 series'. 559 """ 560 # construct list of series 561 if not hasattr(serieslist, "__contains__"): 562 serieslist = [serieslist] 563 564 # get limits 565 xlim = kwargs.pop("xlim", None) 566 ylim = kwargs.pop("ylim", None) 567 568 # set labels 569 xlabel = kwargs.pop("xlabel", "Amplitude") 570 ylabel = kwargs.pop("ylabel", "Fraction of data") 571 title = kwargs.pop("title", "") 572 subtitle = kwargs.pop("subtitle","") 573 574 # get axis scales 575 logx = kwargs.pop("logx", False) 576 logy = kwargs.pop("logy", False) 577 578 # get colorbar options 579 hidden_colorbar = kwargs.pop("hidden_colorbar", False) 580 581 # get savefig option 582 bbox_inches = kwargs.pop("bbox_inches", "tight") 583 584 # get fill 585 fill = kwargs.pop("fill", False) 586 587 # get legend loc 588 loc = kwargs.pop("loc", 0) 589 alpha = kwargs.pop("alpha", 0.8) 590 591 # get colors to use 592 color = kwargs.pop("color", None) 593 if isinstance(color, str): 594 color = color.split(',') 595 if len(color) == 1: 596 color = [color[0]]*len(serieslist) 597 598 599 # 600 # plot 601 # 602 603 plot = plotutils.LineHistogram(xlabel, ylabel, title, subtitle) 604 605 for i, series in enumerate(serieslist): 606 data = series.data.data 607 if nonzero: 608 data = data[data!=0] 609 if color: 610 kwargs["color"] = color[i] 611 plot.add_content(data, normalize=data.size, label=series.name, **kwargs) 612 plot.finalize(loc=loc, alpha=alpha, logx=logx, logy=logy, bar=bar,\ 613 xlim=xlim, num_bins=num_bins, cumulative=cumulative) 614 if hidden_colorbar: 615 plotutils.add_colorbar(plot.ax, visible=False) 616 617 plot.ax.autoscale_view(tight=True, scalex=True, scaley=True) 618 619 if xlim: 620 plot.ax.set_xlim(xlim) 621 if ylim: 622 plot.ax.set_ylim(ylim) 623 624 plotutils.set_minor_ticks(plot.ax) 625 626 # save and close 627 plot.savefig(outfile, bbox_inches=bbox_inches,\ 628 bbox_extra_artists=plot.ax.texts) 629 plot.close()

630

Source Code for Module pylal.plotdata