#------------------------------------------------------------------------------- # Name: SpectralGraph # Purpose: # # Author: Kelvin # # Created: 02/01/2014 #------------------------------------------------------------------------------- import matplotlib.pyplot as plt import matplotlib.widgets as mwidgets import matplotlib.cm as cm from matplotlib.figure import Figure from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas from numpy import * from PyQt4 import QtCore import pywcs class SpectralGraph(FigureCanvas): def __init__(self, Results, maskFits, realFits, maskImage): # Create the figure containing our final image self.fig = Figure() self.channelText = self.fig.text(0.9, 0.15, '') self.velocityText = self.fig.text(0.9, 0.30, '') self.intensityText = self.fig.text(0.9, 0.45, '') self.descriptionText = self.fig.text(0.65, 0.95, '') self.ax = self.fig.add_subplot(111) self.ax.cla() #self.ax.plot(arange(0, self.lenZ), zeros(self.lenZ)) self.ax.set_xlabel('Frequency Channels (GHz)') self.ax.set_ylabel('Intensity (mJ)') self.spectralTitle = self.fig.suptitle("Spectral Characteristic") # Initialise the parent class FigureCanvas.__init__(self, self.fig) self.Results = Results # Those are our data cubes self.maskCube = maskFits[0].data self.realCube = realFits[0].data self.maskImage = maskImage # Those are the size of our data cubes [self.lenZ, self.lenY, self.lenX] = self.maskCube.shape self.maskHeader = maskFits[0].header hdr = bytes(self.maskHeader.tostring(), 'utf-8') imwcs = pywcs.WCS(header = hdr) self.freq = arange(0, self.lenZ) self.vel = imwcs.wcs.crval[2] + imwcs.wcs.cdelt[2]*(self.freq - imwcs.wcs.crpix[2]) self.velunit = str(imwcs.wcs.cunit[2], 'utf-8') self.zAxis = 'vel' fx = maskFits[0].header.get('BMAJ') fy = maskFits[0].header.get('BMIN') self.C = pi*fx*fy/(4*log(2)) self.zmin = 0 self.zmax = self.lenZ self.currentObjID = 0 self.spectralType = 'max' self.trackCursorLink = self.mpl_connect('motion_notify_event', self.trackCursor) self.changeZaxisLink = self.mpl_connect('button_press_event', self.changeZaxis) self.intensityUnit = 'Jy/beam' self.draw() self.show() def plotMax(self): self.spectralType = 'max' self.intensityUnit = 'Jy/beam' self.plotSource() def plotSum(self): self.spectralType = 'sum' self.intensityUnit = 'Jy' self.plotSource() def changeZaxis(self, event): if event.button != 1: return if self.zAxis == 'vel': self.zAxis = 'freq' else: self.zAxis = 'vel' self.plotSource() def setCurrentObjID(self, ObjID): self.currentObjID = ObjID def trackCursor(self, event): if not event.inaxes: return self.channelText.set_text('Freq. Channel: %d'%(int(event.xdata))) self.velocityText.set_text('Velocity: %.3f k%s'%(self.vel[int(event.xdata)]/1000, self.velunit)) self.intensityText.set_text('Intensity: %.3f %s'%(event.ydata, self.intensityUnit)) self.draw() def plotSource(self): self.ax.cla() if self.currentObjID <= 0 or (self.spectralType != 'max' and self.spectralType != 'sum'): self.draw() return ObjStats = self.Results.SkyObjects[self.currentObjID].stats if self.spectralType == 'max': x_peak = ObjStats['X_peak'] y_peak = ObjStats['Y_peak'] flux = self.realCube[:, y_peak, x_peak] self.descriptionText.set_text('Peak Intensity occurs at (x, y) = (%d, %d)'%(x_peak, y_peak)) self.intensityUnit = 'Jy/beam' self.ax.set_ylabel('Intensity (%s)'%self.intensityUnit) elif self.spectralType == 'sum': flux = zeros(self.freq.shape) x1 = ObjStats['X1'] x2 = ObjStats['X2'] y1 = ObjStats['Y1'] y2 = ObjStats['Y2'] for x in range(x1, x2 + 1): for y in range(y1, y2 + 1): if self.maskImage[y][x] == self.currentObjID: flux += self.realCube[:, y, x] flux = self.C*flux self.intensityUnit = 'Jy' self.ax.set_ylabel('Intensity (%s)'%self.intensityUnit) self.spectralTitle.set_text("Spectral Characteristic of Source " + str(self.currentObjID)) self.linef1 = self.ax.axvline(x = self.Results.zf1, color = "r") self.linef2 = self.ax.axvline(x = self.Results.zf2, color = "r") z1 = ObjStats['Z1'] z2 = ObjStats['Z2'] self.linef1 = self.ax.axvline(x = z1, color = "g", linestyle = "--") self.linef2 = self.ax.axvline(x = z2, color = "g", linestyle = "--") self.minflux = flux[0] self.maxflux = flux[0] for z in range(self.lenZ): if z >= self.Results.zf1 and z <= self.Results.zf2: pass else: if flux[z] < self.minflux: self.minflux = flux[z] if flux[z] > self.maxflux: self.maxflux = flux[z] self.fluxrange = self.maxflux - self.minflux self.fluxstep = 0.1*self.fluxrange self.fluxbottom = self.minflux - self.fluxstep self.fluxtop = self.maxflux + self.fluxstep self.fluxsteps = arange(self.fluxbottom, self.fluxtop, self.fluxstep) zf1_ratio = self.Results.zf1/self.lenZ zf2_ratio = self.Results.zf2/self.lenZ for f in self.fluxsteps: self.ax.axhline(y = f, xmin = zf1_ratio, xmax = zf2_ratio, color = "r") zspace = 10 self.zmin = max(0, z1 - zspace) self.zmax = min(z2 + zspace, self.lenZ) self.ax.plot(self.freq, flux) self.ax.set_xlim(0, self.lenZ) if self.zAxis == 'freq': self.ax.set_xlabel('Frequency Channels (Channel Number)') else: xticks = self.ax.get_xticks() velticklabels = [] for xtick in xticks: if xtick >= 0 and xtick < self.lenZ: velticklabel = '%.3f'%(self.vel[int(xtick)]/1000) velticklabels += [velticklabel] self.ax.set_xticklabels(velticklabels) self.ax.set_xlabel('Velocity (%s%s)'%('k',self.velunit)) self.ax.set_ylim(self.fluxbottom, self.fluxtop) self.draw() class ZoomGraph(SpectralGraph): trackCursorSignal = QtCore.pyqtSignal(object) def __init__(self, Results, maskFits, realFits, maskImage): SpectralGraph.__init__(self, Results, maskFits, realFits, maskImage) def trackCursor(self, event): self.trackCursorSignal.emit(event) def plotSource(self): self.ax.cla() if self.currentObjID <= 0 or (self.spectralType != 'max' and self.spectralType != 'sum'): self.draw() return ObjStats = self.Results.SkyObjects[self.currentObjID].stats z1 = ObjStats['Z1'] z2 = ObjStats['Z2'] zspace = 10 zmin = max(0, z1 - zspace) zmax = min(z2 + zspace, self.lenZ - 1) zrange = arange(zmin, zmax + 1) if self.spectralType == 'max': x_peak = ObjStats['X_peak'] y_peak = ObjStats['Y_peak'] flux = self.realCube[zrange, y_peak, x_peak] self.intensityUnit = 'Jy/beam' self.ax.set_ylabel('Intensity (%s)'%self.intensityUnit) elif self.spectralType == 'sum': flux = zeros(zrange.shape) x1 = ObjStats['X1'] x2 = ObjStats['X2'] y1 = ObjStats['Y1'] y2 = ObjStats['Y2'] for x in range(x1, x2 + 1): for y in range(y1, y2 + 1): if self.maskImage[y][x] == self.currentObjID: flux += self.realCube[zrange, y, x] flux = self.C*flux self.intensityUnit = 'Jy' self.ax.set_ylabel('Intensity (%s)'%self.intensityUnit) self.linef1 = self.ax.axvline(x = z1, color = "g", linestyle = "--") self.linef2 = self.ax.axvline(x = z2, color = "g", linestyle = "--") self.minflux = flux[0] self.maxflux = flux[0] for i in range(len(flux)): if zrange[i] >= self.Results.zf1 and zrange[i] <= self.Results.zf2: pass else: if flux[i] < self.minflux: self.minflux = flux[i] if flux[i] > self.maxflux: self.maxflux = flux[i] self.fluxrange = self.maxflux - self.minflux self.fluxstep = 0.1*self.fluxrange self.fluxbottom = self.minflux - self.fluxstep self.fluxtop = self.maxflux + self.fluxstep self.fluxsteps = arange(self.fluxbottom, self.fluxtop, self.fluxstep) self.ax.plot(zrange, flux) self.ax.set_xlim(zmin, zmax) if self.zAxis == 'freq': self.ax.set_xlabel('Frequency Channels (Channel Number)') else: xticks = self.ax.get_xticks() velticklabels = [] for xtick in xticks: if xtick >= 0 and xtick < self.lenZ: velticklabel = velticklabel = '%4.0f'%(self.vel[int(xtick)]/1000) velticklabels += [velticklabel] self.ax.set_xticklabels(velticklabels) self.ax.set_xlabel('Velocity (%s%s)'%('k',self.velunit)) self.ax.set_ylim(self.fluxbottom, self.fluxtop) self.spectralTitle.set_text("Source " + str(self.currentObjID)) self.draw()