diff --git a/reduction/obs_report/AR1_ptuse_obs_report_template.ipynb b/reduction/obs_report/AR1_ptuse_obs_report_template.ipynb deleted file mode 100644 index ab62a7d4..00000000 --- a/reduction/obs_report/AR1_ptuse_obs_report_template.ipynb +++ /dev/null @@ -1,784 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "# AR1 PTUSE observation report" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "## Files to be processed" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "directory = " - ] - }, - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "## Set up environment" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "%pylab inline" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "from matplotlib import cm\n", - "from matplotlib import colors\n", - "import psrchive as psr\n", - "from coast_guard import cleaners\n", - "from coast_guard import clean_utils\n", - "from coast_guard import utils\n", - "import os\n", - "import sys" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "## Load data" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "files = []\n", - "for file in os.listdir(directory):\n", - " if file.endswith(\".ar\"):\n", - " files.append(file)\n", - "files.sort()\n", - "archives = []\n", - "archives = [psr.Archive_load(directory + file) for file in files]\n", - "print archives[0]\n", - "for i in range(1, len(archives)):\n", - " print archives[i]\n", - " archives[0].append(archives[i])" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "## Display metadata" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "nBin = archives[0].get_nbin()\n", - "nChan = archives[0].get_nchan()\n", - "nPol = archives[0].get_npol()\n", - "nSubint = archives[0].get_nsubint()\n", - "obsType = archives[0].get_type()\n", - "telescopeName = archives[0].get_telescope()\n", - "sourceName = archives[0].get_source()\n", - "RA = archives[0].get_coordinates().ra()\n", - "Dec = archives[0].get_coordinates().dec()\n", - "centreFrequency = archives[0].get_centre_frequency()\n", - "bandwidth = archives[0].get_bandwidth()\n", - "DM = archives[0].get_dispersion_measure()\n", - "RM = archives[0].get_rotation_measure()\n", - "isDedispersed = archives[0].get_dedispersed()\n", - "isFaradayRotated = archives[0].get_faraday_corrected()\n", - "isPolCalib = archives[0].get_poln_calibrated()\n", - "dataUnits = archives[0].get_scale()\n", - "dataState = archives[0].get_state()\n", - "obsDuration = archives[0].integration_length()\n", - "receiverName = archives[0].get_receiver_name()\n", - "receptorBasis = archives[0].get_basis()\n", - "backendName = archives[0].get_backend_name()\n", - "lowFreq = archives[0].get_centre_frequency() - archives[0].get_bandwidth() / 2.0\n", - "highFreq = archives[0].get_centre_frequency() + archives[0].get_bandwidth() / 2.0" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, - "outputs": [], - "source": [ - "print \"nbin Number of pulse phase bins %s\" % nBin\n", - "print \"nchan Number of frequency channels %s\" % nChan\n", - "print \"npol Number of polarizations %s\" % nPol\n", - "print \"nsubint Number of sub-integrations %s\" % nSubint\n", - "print \"type Observation type %s\" % obsType\n", - "print \"site Telescope name %s\" % telescopeName\n", - "print \"name Source name %s\" % sourceName\n", - "print \"coord Source coordinates %s%s\" % (RA.getHMS(), Dec.getDMS())\n", - "print \"freq Centre frequency (MHz) %s\" % centreFrequency\n", - "print \"bw Bandwidth (MHz) %s\" % bandwidth\n", - "print \"dm Dispersion measure (pc/cm^3) %s\" % DM\n", - "print \"rm Rotation measure (rad/m^2) %s\" % RM\n", - "print \"dmc Dispersion corrected %s\" % isDedispersed\n", - "print \"rmc Faraday Rotation corrected %s\" % isFaradayRotated\n", - "print \"polc Polarization calibrated %s\" % isPolCalib\n", - "print \"scale Data units %s\" % dataUnits\n", - "print \"state Data state %s\" % dataState\n", - "print \"length Observation duration (s) %s\" % obsDuration\n", - "print \"rcvr:name Receiver name %s\" % receiverName\n", - "print \"rcvr:basis Basis of receptors %s\" % receptorBasis\n", - "print \"be:name Name of the backend instrument %s\" % backendName" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "## Clean archive, display RFI stats and mask" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "cleanRFI = archives[0].clone()\n", - "cleaner = cleaners.load_cleaner(\"surgical\")\n", - "cleaner.parse_config_string(\"chan_numpieces=1,subint_numpieces=1,chanthresh=3,subintthresh=3\")\n", - "cleaner.run(cleanRFI)\n", - "weights = cleanRFI.get_weights()\n", - "i = j = counter = channel = 0\n", - "emptyChan = [i for i in range(nChan)]\n", - "for i in range(nSubint):\n", - " channel = 0\n", - " del emptyChan[:]\n", - " for j in range(nChan):\n", - " if weights[i][j] == 0.0:\n", - " counter += 1\n", - " channel += 1\n", - " emptyChan.append(j)\n", - " subintProc = (float(channel)/float(nChan))*100\n", - " print \"Subint %s has %s channels (%.2f%%) set to zero.\" % (i, channel, subintProc)\n", - "percentage = (float(counter)/float(weights.size))*100\n", - "print \"%s data points out of %s with weights set to zero.\" % (counter, weights.size)\n", - "print \"%.2f%% data was set to zero.\" % (percentage)\n", - "fig, ax1 = plt.subplots(1, 1, figsize = [15, 10], tight_layout=\"false\")\n", - "ax1.set_title(sourceName)\n", - "ax1.set_title(\"RFI mask\", loc=\"left\")\n", - "ax1.set_ylabel(\"Channel number\")\n", - "ax1.yaxis.set_ticks(np.arange(0, nChan - 1, 200))\n", - "ax1.set_xlabel(\"Subint number\")\n", - "ax1Secondary = ax1.twinx()\n", - "ax1Secondary.set_ylabel(\"Frequency (MHz)\")\n", - "ax1Secondary.set_ylim(lowFreq, highFreq)\n", - "ax1Secondary.yaxis.set_ticks(np.arange(lowFreq, highFreq, 25))\n", - "ax1.imshow(weights.T, origin=\"lower\", aspect=\"auto\", cmap=colors.ListedColormap(['red', 'white']), interpolation=\"none\", extent = (0, nSubint - 1, 0, nChan - 1))" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "## Plot average profiles (cleaned Stokes data)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "cleanProfile = cleanRFI.clone()\n", - "cleanProfile.dedisperse()\n", - "cleanProfile.remove_baseline()\n", - "cleanProfile.centre_max_bin()\n", - "cleanProfile.tscrunch()\n", - "cleanProfile.fscrunch()\n", - "cleanProfile.convert_state(\"Stokes\")\n", - "cleanProfileData = cleanProfile.get_data()\n", - "cleanProfileData = cleanProfileData.squeeze()\n", - "cleanMinMaxVals = np.array([cleanProfileData[0, :].max(), cleanProfileData[1, :].max(), cleanProfileData[2, :].max(), cleanProfileData[3, :].max(), cleanProfileData[0, :].min(), cleanProfileData[1, :].min(), cleanProfileData[2, :].min(), cleanProfileData[3, :].min()])\n", - "dirtyProfile = archives[0].clone()\n", - "dirtyProfile.dedisperse()\n", - "dirtyProfile.remove_baseline()\n", - "dirtyProfile.centre_max_bin()\n", - "dirtyProfile.tscrunch()\n", - "dirtyProfile.fscrunch()\n", - "dirtyProfile.convert_state(\"Stokes\")\n", - "dirtyProfileData = dirtyProfile.get_data()\n", - "dirtyProfileData = dirtyProfileData.squeeze()\n", - "dirtyMinMaxVals = np.array([dirtyProfileData[0, :].min(), dirtyProfileData[1, :].min(), dirtyProfileData[2, :].min(), dirtyProfileData[3, :].min(), dirtyProfileData[0, :].max(), dirtyProfileData[1, :].max(), dirtyProfileData[2, :].max(), dirtyProfileData[3, :].max()])\n", - "fig, ((ax1, ax2)) = plt.subplots(2, 1, sharex=\"col\", figsize = [15, 10], tight_layout=\"false\")\n", - "ax1.set_title(sourceName)\n", - "ax1.text(20, dirtyProfileData.max(), \"Uncleaned data\", fontsize=\"large\")\n", - "ax1.set_xlim(0, nBin - 1)\n", - "ax1.set_ylabel(\"Flux (a.u.)\")\n", - "ax1.set_ylim((dirtyMinMaxVals.min() * 1.2), dirtyMinMaxVals.max() + 0.1 * dirtyMinMaxVals.max())\n", - "ax1.xaxis.set_ticks(np.arange(0, nBin - 1, 50))\n", - "ax1.plot(dirtyProfileData[0, :])\n", - "ax1.plot(dirtyProfileData[1, :])\n", - "ax1.plot(dirtyProfileData[2, :])\n", - "ax1.plot(dirtyProfileData[3, :])\n", - "ax1.plot(dirtyProfileData[0, :], \"k\")\n", - "ax1.plot(dirtyProfileData[1, :], \"r\")\n", - "ax1.plot(dirtyProfileData[2, :], \"g\")\n", - "ax1.plot(dirtyProfileData[3, :], \"b\")\n", - "ax2.text(20, cleanProfileData.max(), \"Cleaned data\", fontsize=\"large\")\n", - "ax2.set_xlabel(\"Pulse phase (bins)\")\n", - "ax2.set_xlim(0, nBin - 1)\n", - "ax2.set_ylabel(\"Flux (a.u.)\")\n", - "ax2.set_ylim((cleanMinMaxVals.min() * 1.2), cleanMinMaxVals.max() + 0.1 * cleanMinMaxVals.max())\n", - "ax2.xaxis.set_ticks(np.arange(0, nBin - 1, 50))\n", - "ax2.plot(cleanProfileData[0, :], \"k\")\n", - "ax2.plot(cleanProfileData[1, :], \"r\")\n", - "ax2.plot(cleanProfileData[2, :], \"g\")\n", - "ax2.plot(cleanProfileData[3, :], \"b\")" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "## Plot subintegration stacks (cleaned Coherency data)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "dirtyStack = cleanRFI.clone()\n", - "dirtyStack.dedisperse()\n", - "dirtyStack.remove_baseline()\n", - "dirtyStack.centre_max_bin()\n", - "dirtyStack.fscrunch()\n", - "dirtyStack.pscrunch()\n", - "dirtyStack_data = dirtyStack.get_data().squeeze()\n", - "cleanStack = cleanRFI.clone()\n", - "cleanStack.dedisperse()\n", - "cleanStack.remove_baseline()\n", - "cleanStack.centre_max_bin()\n", - "cleanStack.fscrunch()\n", - "cleanStack.pscrunch()\n", - "cleanStack_data = cleanStack.get_data().squeeze()\n", - "fig, ((ax1, ax2)) = plt.subplots(1, 2, sharey=\"row\", figsize = [15, 10], tight_layout=\"false\")\n", - "ax1.set_title(sourceName)\n", - "ax1.set_title(\"Uncleaned data\", loc=\"left\")\n", - "ax1.set_ylabel(\"Index\")\n", - "ax1.set_xlabel(\"Pulse phase (bin)\")\n", - "ax1.xaxis.set_ticks(np.arange(0, nBin - 1, 100))\n", - "ax1.imshow(dirtyStack_data, cmap=cm.afmhot, aspect=\"auto\", vmax=dirtyStack_data.max()*0.75, interpolation=\"none\", extent = (0, nBin - 1, 0, nSubint - 1))\n", - "ax2.set_title(sourceName)\n", - "ax2.set_title(\"Cleaned data\", loc=\"left\")\n", - "ax2.set_xlabel(\"Pulse phase (bin)\")\n", - "ax2.xaxis.set_ticks(np.arange(0, nBin - 1, 100))\n", - "ax2.imshow(cleanStack_data, cmap=cm.afmhot, aspect=\"auto\", vmax=cleanStack_data.max()*0.75, interpolation=\"none\", extent = (0, nBin - 1, 0, nSubint - 1))" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "## Plot phase vs. frequency image of flux (cleaned Coherency data)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "coherencyFluxPhase = cleanRFI.clone()\n", - "coherencyFluxPhase.dedisperse()\n", - "coherencyFluxPhase.remove_baseline()\n", - "coherencyFluxPhase.tscrunch()\n", - "coherencyFluxPhase.centre_max_bin()\n", - "coherencyFluxPhase_data = coherencyFluxPhase.get_data().squeeze()\n", - "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=\"col\", sharey=\"row\", figsize = [15, 10], tight_layout=\"false\")\n", - "fig.suptitle(sourceName, fontsize=\"large\", y=1.015)\n", - "fig.text(0.05, 1.0, \"Cleaned data\", fontsize=\"large\")\n", - "ax1.set_title(\"XX\")\n", - "ax1.set_ylabel(\"Frequency (MHz)\")\n", - "ax1.xaxis.set_ticks(np.arange(0, nBin - 1, 100))\n", - "ax1.yaxis.set_ticks(np.arange(0, nChan - 1, 50))\n", - "ax1.imshow(coherencyFluxPhase_data[0, :, :], cmap=cm.afmhot, aspect=\"auto\", interpolation=\"none\", extent = (0, nBin - 1.0, lowFreq, highFreq))\n", - "ax2.set_title(\"YY\")\n", - "ax2.xaxis.set_ticks(np.arange(0, nBin - 1, 100))\n", - "ax2.yaxis.set_ticks(np.arange(0, nChan - 1, 50))\n", - "ax2Secondary = ax2.twinx()\n", - "ax2Secondary.set_ylabel(\"Channel number\")\n", - "ax2Secondary.set_ylim(0, nChan - 1)\n", - "ax2Secondary.yaxis.set_ticks(np.arange(0, nChan - 1 , 200))\n", - "ax2.imshow(coherencyFluxPhase_data[1, :, :], cmap=cm.afmhot, aspect=\"auto\", interpolation=\"none\", extent = (0, nBin - 1.0, lowFreq, highFreq))\n", - "ax3.set_title(\"Re XY\")\n", - "ax3.set_xlabel(\"Pulse phase (bin)\")\n", - "ax3.set_ylabel(\"Frequency (MHz)\")\n", - "ax3.xaxis.set_ticks(np.arange(0, nBin - 1, 100))\n", - "ax3.yaxis.set_ticks(np.arange(0, nChan - 1, 50))\n", - "ax3.imshow(coherencyFluxPhase_data[2, :, :], cmap=cm.afmhot, aspect=\"auto\", interpolation=\"none\", extent = (0, nBin - 1.0, lowFreq, highFreq))\n", - "ax4.set_title(\"Im XY\")\n", - "ax4.set_xlabel(\"Pulse phase (bin)\")\n", - "ax4.xaxis.set_ticks(np.arange(0, nBin - 1, 100))\n", - "ax4.yaxis.set_ticks(np.arange(0, nChan - 1, 50))\n", - "ax4Secondary = ax4.twinx()\n", - "ax4Secondary.set_ylabel(\"Channel number\")\n", - "ax4Secondary.set_ylim(0, nChan - 1)\n", - "ax4Secondary.yaxis.set_ticks(np.arange(0, nChan - 1 , 200))\n", - "ax4.imshow(coherencyFluxPhase_data[3, :, :], cmap=cm.afmhot, aspect=\"auto\", interpolation=\"none\", extent = (0, nBin - 1.0, lowFreq, highFreq))" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "## Plot phase vs. frequency image of flux (cleaned Stokes data)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "stokesFluxPhase = cleanRFI.clone()\n", - "stokesFluxPhase.dedisperse()\n", - "stokesFluxPhase.remove_baseline()\n", - "stokesFluxPhase.tscrunch()\n", - "stokesFluxPhase.centre_max_bin()\n", - "stokesFluxPhase.convert_state(\"Stokes\")\n", - "stokesFluxPhase_data = stokesFluxPhase.get_data().squeeze()\n", - "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=\"col\", sharey=\"row\", figsize = [15, 10], tight_layout=\"false\")\n", - "fig.suptitle(sourceName, fontsize=\"large\", y=1.01)\n", - "fig.text(0.05, 1.0, \"Cleaned data\", fontsize=\"large\")\n", - "ax1.set_title(\"I\")\n", - "ax1.set_ylabel(\"Frequency (MHz)\")\n", - "ax1.xaxis.set_ticks(np.arange(0, nBin - 1, 100))\n", - "ax1.yaxis.set_ticks(np.arange(0, nChan - 1, 50))\n", - "ax1.imshow(stokesFluxPhase_data[0, :, :], cmap=cm.afmhot, aspect=\"auto\", interpolation=\"none\", extent = (0, nBin - 1.0, lowFreq, highFreq))\n", - "ax2.set_title(\"Q\")\n", - "ax2.xaxis.set_ticks(np.arange(0, nBin - 1, 100))\n", - "ax2.yaxis.set_ticks(np.arange(0, nChan - 1, 50))\n", - "ax2Secondary = ax2.twinx()\n", - "ax2Secondary.set_ylabel(\"Channel number\")\n", - "ax2Secondary.set_ylim(0, nChan - 1)\n", - "ax2Secondary.yaxis.set_ticks(np.arange(0, nChan - 1 , 200))\n", - "ax2.imshow(stokesFluxPhase_data[1, :, :], cmap=cm.afmhot, aspect=\"auto\", interpolation=\"none\", extent = (0, nBin - 1.0, lowFreq, highFreq))\n", - "ax3.set_title(\"U\")\n", - "ax3.set_xlabel(\"Pulse phase (bin)\")\n", - "ax3.set_ylabel(\"Frequency (MHz)\")\n", - "ax3.xaxis.set_ticks(np.arange(0, nBin - 1, 100))\n", - "ax3.yaxis.set_ticks(np.arange(0, nChan - 1, 50))\n", - "ax3.imshow(stokesFluxPhase_data[2, :, :], cmap=cm.afmhot, aspect=\"auto\", interpolation=\"none\", extent = (0, nBin - 1.0, lowFreq, highFreq))\n", - "ax4.set_title(\"V\")\n", - "ax4.set_xlabel(\"Pulse phase (bin)\")\n", - "ax4.xaxis.set_ticks(np.arange(0, nBin - 1, 100))\n", - "ax4.yaxis.set_ticks(np.arange(0, nChan - 1, 50))\n", - "ax4Secondary = ax4.twinx()\n", - "ax4Secondary.set_ylabel(\"Channel number\")\n", - "ax4Secondary.set_ylim(0, nChan - 1)\n", - "ax4Secondary.yaxis.set_ticks(np.arange(0, nChan - 1 , 200))\n", - "ax4.imshow(stokesFluxPhase_data[3, :, :], cmap=cm.afmhot, aspect=\"auto\", interpolation=\"none\", extent = (0, nBin - 1.0, lowFreq, highFreq))" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "## Plot bandpass (uncleaned Stokes data)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "dirtyBandpass = archives[0].clone()\n", - "dirtyBandpass.tscrunch()\n", - "(dirtyBandpassMean, dirtyBandpassVariance) = dirtyBandpass.get_Integration(0).baseline_stats()\n", - "dirtyBandpassMaxMinVals = np.array([dirtyBandpassMean[0, :].min(), dirtyBandpassMean[0, :].max(), dirtyBandpassMean[1, :].min(), dirtyBandpassMean[1, :].max()])\n", - "dirtyBandpassStokes = archives[0].clone()\n", - "dirtyBandpassStokes.convert_state(\"Stokes\")\n", - "dirtyBandpassStokes.tscrunch()\n", - "(dirtyBandpassStokesMean, dirtyBandpassStokesVariance) = dirtyBandpassStokes.get_Integration(0).baseline_stats()\n", - "dirtyBandpassStokesMinMaxVals = np.array([dirtyBandpassStokesMean[0, :].min(), dirtyBandpassStokesMean[0, :].max()])\n", - "fig, ((ax1, ax2, ax3 )) = plt.subplots(3, 1, sharex=\"col\", figsize = [15, 10], tight_layout=\"false\")\n", - "ax1.set_title(sourceName)\n", - "ax1.set_title(\"Unclean data\", loc=\"left\")\n", - "ax1.text(np.int(nChan*0.0125), dirtyBandpassMaxMinVals.max(), \"XX\", fontsize=\"large\")\n", - "ax1.set_xlim(0, nChan - 1)\n", - "ax1.set_ylabel(\"Flux (a.u.)\")\n", - "ax1.set_ylim(dirtyBandpassMaxMinVals.min() - (0.1 * np.abs(dirtyBandpassMaxMinVals.min())), dirtyBandpassMaxMinVals.max() + (0.1 * np.abs(dirtyBandpassMaxMinVals.max())))\n", - "ax1.xaxis.set_ticks(np.arange(0, nChan - 1, 100))\n", - "ax1.plot(dirtyBandpassMean[0, :])\n", - "ax2.text(np.int(nChan*0.0125), dirtyBandpassMaxMinVals.max(), \"YY\", fontsize=\"large\")\n", - "ax2.set_xlim(0, nChan - 1)\n", - "ax2.set_ylabel(\"Flux (a.u.)\")\n", - "ax2.set_ylim(dirtyBandpassMaxMinVals.min() - (0.1 * np.abs(dirtyBandpassMaxMinVals.min())), dirtyBandpassMaxMinVals.max() + (0.1 * np.abs(dirtyBandpassMaxMinVals.max())))\n", - "ax2.xaxis.set_ticks(np.arange(0, nChan - 1, 100))\n", - "ax2.plot(dirtyBandpassMean[1, :])\n", - "ax3.text(np.int(nChan*0.0125), dirtyBandpassStokesMinMaxVals.max(), \"Stokes I\", fontsize=\"large\")\n", - "ax3.set_xlabel(\"Channel number\")\n", - "ax3.set_ylabel(\"Flux (a.u.)\")\n", - "ax3.set_ylim(dirtyBandpassStokesMinMaxVals.min() - (0.1 * np.abs(dirtyBandpassStokesMinMaxVals.min())), dirtyBandpassStokesMinMaxVals.max() + (0.1 * np.abs(dirtyBandpassStokesMinMaxVals.max())))\n", - "ax3.xaxis.set_ticks(np.arange(0, nChan - 1, 100))\n", - "ax3Secondary = ax3.twiny()\n", - "ax3Secondary.set_frame_on(True)\n", - "ax3Secondary.patch.set_visible(False)\n", - "ax3Secondary.xaxis.set_ticks_position('bottom')\n", - "ax3Secondary.set_xlabel(\"Frequency (MHz)\")\n", - "ax3Secondary.xaxis.set_label_position('bottom')\n", - "ax3Secondary.spines['bottom'].set_position(('outward', 50))\n", - "ax3Secondary.set_xlim(lowFreq, highFreq)\n", - "ax3.plot(dirtyBandpassStokesMean[0, :])" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "## Plot bandpass (cleaned Stokes data)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "cleanBandpass = cleanRFI.clone()\n", - "cleanBandpass.tscrunch()\n", - "(cleanBandpassMean, cleanBandpassVariance) = cleanBandpass.get_Integration(0).baseline_stats()\n", - "cleanBandpassMaxMinVals = np.array([cleanBandpassMean[0, :].min(), cleanBandpassMean[0, :].max(), cleanBandpassMean[1, :].min(), cleanBandpassMean[1, :].max()])\n", - "cleanBandpassStokes = cleanBandpass.clone()\n", - "cleanBandpassStokes.convert_state(\"Stokes\")\n", - "cleanBandpassStokes.tscrunch()\n", - "(cleanBandpassStokesMean, cleanBandpassStokesVariance) = cleanBandpassStokes.get_Integration(0).baseline_stats()\n", - "cleanBandpassStokesMinMaxVals = np.array([cleanBandpassStokesMean[0, :].min(), cleanBandpassStokesMean[0, :].max()])\n", - "fig, ((ax1, ax2, ax3 )) = plt.subplots(3, 1, sharex=\"col\", figsize = [15, 10], tight_layout=\"false\")\n", - "ax1.set_title(sourceName)\n", - "ax1.set_title(\"Clean data\", loc=\"left\")\n", - "ax1.text(np.int(nChan*0.0125), cleanBandpassMaxMinVals.max(), \"XX\", fontsize=\"large\")\n", - "ax1.set_xlim(0, nChan - 1)\n", - "ax1.set_ylabel(\"Flux (a.u.)\")\n", - "ax1.set_ylim(cleanBandpassMaxMinVals.min() - (0.1 * np.abs(cleanBandpassMaxMinVals.min())), cleanBandpassMaxMinVals.max() + (0.1 * np.abs(cleanBandpassMaxMinVals.max())))\n", - "ax1.xaxis.set_ticks(np.arange(0, nChan - 1, 100))\n", - "ax1.plot(cleanBandpassMean[0, :])\n", - "ax2.text(np.int(nChan*0.0125), cleanBandpassMaxMinVals.max(), \"YY\", fontsize=\"large\")\n", - "ax2.set_xlim(0, nChan - 1)\n", - "ax2.set_ylabel(\"Flux (a.u.)\")\n", - "ax2.set_ylim(cleanBandpassMaxMinVals.min() - (0.1 * np.abs(cleanBandpassMaxMinVals.min())), cleanBandpassMaxMinVals.max() + (0.1 * np.abs(cleanBandpassMaxMinVals.max())))\n", - "ax2.xaxis.set_ticks(np.arange(0, nChan - 1, 100))\n", - "ax2.plot(cleanBandpassMean[1, :])\n", - "ax3.text(np.int(nChan*0.0125), cleanBandpassStokesMinMaxVals.max(), \"Stokes I\", fontsize=\"large\")\n", - "ax3.set_xlabel(\"Channel number\")\n", - "ax3.set_ylabel(\"Flux (a.u.)\")\n", - "ax3.set_ylim(cleanBandpassStokesMinMaxVals.min() - (0.1 * np.abs(cleanBandpassStokesMinMaxVals.min())), cleanBandpassStokesMinMaxVals.max() + (0.1 * np.abs(cleanBandpassStokesMinMaxVals.max())))\n", - "ax3.xaxis.set_ticks(np.arange(0, nChan - 1, 100))\n", - "ax3Secondary = ax3.twiny()\n", - "ax3Secondary.set_frame_on(True)\n", - "ax3Secondary.patch.set_visible(False)\n", - "ax3Secondary.xaxis.set_ticks_position('bottom')\n", - "ax3Secondary.set_xlabel(\"Frequency (MHz)\")\n", - "ax3Secondary.xaxis.set_label_position('bottom')\n", - "ax3Secondary.spines['bottom'].set_position(('outward', 50))\n", - "ax3Secondary.set_xlim(lowFreq, highFreq)\n", - "ax3.plot(cleanBandpassStokesMean[0, :])" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "## Plot dynamic baseline spectrum (uncleaned Coherency data)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "dynamicSpectrum = archives[0].clone()\n", - "mean = np.zeros((nSubint, nPol, nChan))\n", - "variance = np.zeros((nSubint, nPol, nChan))\n", - "for subint in range(nSubint):\n", - " m, v = dynamicSpectrum.get_Integration(subint).baseline_stats()\n", - " mean[subint] = m[:, :]\n", - "meanXXMin = mean[:, 0, :].min()\n", - "meanXXMax = mean[:, 0, :].max()\n", - "meanYYMin = mean[:, 1, :].min()\n", - "meanYYMax = mean[:, 1, :].max()\n", - "meanReXYMin = mean[:, 2, :].min()\n", - "meanReXYMax = mean[:, 2, :].max()\n", - "meanImXYMin = mean[:, 3, :].min()\n", - "meanImXYMax = mean[:, 3, :].max()\n", - "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=\"col\", sharey=\"row\", figsize = [15, 10], tight_layout=\"false\")\n", - "fig.suptitle(sourceName, fontsize=\"large\", y=1.01)\n", - "fig.text(0.05, 1.0, \"Uncleaned data\", fontsize=\"large\")\n", - "ax1.set_title(\"XX\")\n", - "ax1.set_ylabel(\"Channel number\")\n", - "ax1.yaxis.set_ticks(np.arange(0, nChan - 1, 200))\n", - "ax1.imshow(mean[:, 0, :].T, origin=\"lower\", aspect=\"auto\", cmap=cm.afmhot, norm=colors.SymLogNorm(10, linscale=1.0,vmin=meanXXMin, vmax=meanXXMax), interpolation=\"none\")\n", - "ax2.set_title(\"YY\")\n", - "ax2.yaxis.set_ticks(np.arange(0, nChan - 1, 200))\n", - "ax2Secondary = ax2.twinx()\n", - "ax2Secondary.set_ylabel(\"Frequency (MHz)\")\n", - "ax2Secondary.set_ylim(lowFreq, highFreq)\n", - "ax2Secondary.yaxis.set_ticks(np.arange(lowFreq, highFreq, 25))\n", - "ax2.imshow(mean[:, 1, :].T, origin=\"lower\", aspect=\"auto\", cmap=cm.afmhot, norm=colors.SymLogNorm(10, linscale=1.0,vmin=meanYYMin, vmax=meanYYMax), interpolation=\"none\")\n", - "ax3.set_title(\"Re XY\")\n", - "ax3.set_xlabel(\"Subint number\")\n", - "ax3.set_ylabel(\"Channel number\")\n", - "ax3.yaxis.set_ticks(np.arange(0, nChan - 1, 200))\n", - "ax3.imshow(mean[:, 2, :].T, origin=\"lower\", aspect=\"auto\", cmap=cm.afmhot, norm=colors.SymLogNorm(10, linscale=1.0,vmin=meanReXYMin, vmax=meanReXYMax), interpolation=\"none\")\n", - "ax4.set_title(\"Im XY\")\n", - "ax4.set_xlabel(\"Subint number\")\n", - "ax4.yaxis.set_ticks(np.arange(0, nChan - 1, 200))\n", - "ax4Secondary = ax4.twinx()\n", - "ax4Secondary.set_ylabel(\"Frequency (MHz)\")\n", - "ax4Secondary.set_ylim(lowFreq, highFreq)\n", - "ax4Secondary.yaxis.set_ticks(np.arange(lowFreq, highFreq, 25))\n", - "ax4.imshow(mean[:, 3, :].T, origin=\"lower\", aspect=\"auto\", cmap=cm.afmhot, norm=colors.SymLogNorm(10, linscale=1.0,vmin=meanImXYMin, vmax=meanImXYMax), interpolation=\"none\")" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "## Plot dynamic baseline spectrum (cleaned Stokes data)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "dynamicSpectrum = cleanRFI.clone()\n", - "dynamicSpectrum.convert_state(\"Stokes\")\n", - "mean = np.zeros((nSubint, nPol, nChan))\n", - "variance = np.zeros((nSubint, nPol, nChan))\n", - "for subint in range(nSubint):\n", - " m, v = dynamicSpectrum.get_Integration(subint).baseline_stats()\n", - " variance[subint] = v[:, :]\n", - "varianceXXMin = variance[:, 0, :].min()\n", - "varianceXXMax = variance[:, 0, :].max()\n", - "varianceYYMin = variance[:, 1, :].min()\n", - "varianceYYMax = variance[:, 1, :].max()\n", - "varianceReXYMin = variance[:, 2, :].min()\n", - "varianceReXYMax = variance[:, 2, :].max()\n", - "varianceImXYMin = variance[:, 3, :].min()\n", - "varianceImXYMax = variance[:, 3, :].max()\n", - "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=\"col\", sharey=\"row\", figsize = [15, 10], tight_layout=\"false\")\n", - "fig.suptitle(sourceName, fontsize=\"large\", y=1.01)\n", - "fig.text(0.05, 1.0, \"Cleaned data\", fontsize=\"large\")\n", - "ax1.set_title(\"I\")\n", - "ax1.set_ylabel(\"Channel number\")\n", - "ax1.yaxis.set_ticks(np.arange(0, nChan - 1, 200))\n", - "ax1.imshow(variance[:, 0, :].T, origin=\"lower\", aspect=\"auto\", cmap=cm.afmhot, norm=colors.Normalize(vmin=varianceXXMin, vmax=varianceXXMax), interpolation=\"none\")\n", - "ax2.set_title(\"Q\")\n", - "ax2.yaxis.set_ticks(np.arange(0, nChan - 1, 200))\n", - "ax2Secondary = ax2.twinx()\n", - "ax2Secondary.set_ylabel(\"Frequency (MHz)\")\n", - "ax2Secondary.set_ylim(lowFreq, highFreq)\n", - "ax2Secondary.yaxis.set_ticks(np.arange(lowFreq, highFreq, 25))\n", - "ax2.imshow(variance[:, 1, :].T, origin=\"lower\", aspect=\"auto\", cmap=cm.afmhot, norm=colors.Normalize(vmin=varianceYYMin, vmax=varianceYYMax), interpolation=\"none\")\n", - "ax3.set_title(\"U\")\n", - "ax3.set_xlabel(\"Subint number\")\n", - "ax3.set_ylabel(\"Channel number\")\n", - "ax3.yaxis.set_ticks(np.arange(0, nChan - 1, 200))\n", - "ax3.imshow(variance[:, 2, :].T, origin=\"lower\", aspect=\"auto\", cmap=cm.afmhot, norm=colors.Normalize(vmin=varianceReXYMin, vmax=varianceReXYMax), interpolation=\"none\")\n", - "ax4.set_title(\"V\")\n", - "ax4.set_xlabel(\"Subint number\")\n", - "ax4.yaxis.set_ticks(np.arange(0, nChan - 1, 200))\n", - "ax4Secondary = ax4.twinx()\n", - "ax4Secondary.set_ylabel(\"Frequency (MHz)\")\n", - "ax4Secondary.set_ylim(lowFreq, highFreq)\n", - "ax4Secondary.yaxis.set_ticks(np.arange(lowFreq, highFreq, 25))\n", - "ax4.imshow(variance[:, 3, :].T, origin=\"lower\", aspect=\"auto\", cmap=cm.afmhot, norm=colors.Normalize(vmin=varianceImXYMin, vmax=varianceImXYMax), interpolation=\"none\")" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, - "source": [ - "## Plot signal-to-noise ratio vs. subintegration (cleaned Coherency data)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "archiveSNR = cleanRFI.clone()\n", - "newIntegration = archiveSNR.get_Integration(0).total()\n", - "snrData = np.zeros(nSubint)\n", - "for iSubintegration in range (1, nSubint):\n", - " nextIntegration = archiveSNR.get_Integration(iSubintegration).total()\n", - " newIntegration.combine(nextIntegration)\n", - " profile = newIntegration.get_Profile(0, 0)\n", - " snrData[iSubintegration] = profile.snr()\n", - "fig, ax1 = plt.subplots(1, 1, sharex=\"col\", sharey=\"row\", figsize = [15, 5], tight_layout=\"false\")\n", - "fig.suptitle(sourceName, fontsize=\"large\", y=1.01)\n", - "fig.text(0.05, 1.0, \"Cleaned data\", fontsize=\"large\")\n", - "ax1.set_xlabel(\"Subint number\")\n", - "ax1.set_ylabel(\"Signal-to-noise ratio\")\n", - "ax1.plot(snrData)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 2", - "language": "python", - "name": "python2" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 2 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.9" - } - }, - "nbformat": 4, - "nbformat_minor": 1 -} diff --git a/reduction/obs_report/bf_fold_obs_report.py b/reduction/obs_report/bf_fold_obs_report.py new file mode 100644 index 00000000..862fc5e2 --- /dev/null +++ b/reduction/obs_report/bf_fold_obs_report.py @@ -0,0 +1,52 @@ +#!/usr/bin/env python + +import os +import argparse +import nbformat +from subprocess import call +from tempfile import mkdtemp + +def build_ipynb_bf_fold_obs_report_command(directory, template): + print directory + bf_fold_obs_report_template_filename = template + directory_base = os.path.basename(directory) + bf_fold_obs_report_output_filename = 'bf_fold_obs_report_%s.ipynb' % directory_base + file_in = open(bf_fold_obs_report_template_filename) + file_out = open(bf_fold_obs_report_output_filename, 'w') + nb = nbformat.reader.read(file_in) + # the format of notebooks changed... find out the version and treat accordingly + ipy_version = int(nb['nbformat']) - 1 + if ipy_version == 2: + for sheet in nb['worksheets']: + for cell in sheet['cells']: + if cell['input'].startswith('directory ='): + cell['input'] = "directory = '%s'" % directory + else: + for cell in nb['cells']: + if cell['source'].startswith('directory ='): + cell['source'] = "directory = '%s'" % directory + nbformat.write(nb, file_out) + file_in.close() + file_out.close() + return 'runipy --overwrite --skip-exceptions %s' % bf_fold_obs_report_output_filename + +def generate_ipynb_bf_fold_obs_report(directory, template): + # call obs report + obs_report_command = build_ipynb_bf_fold_obs_report_command(directory, template) + print 'Command executed:' + print obs_report_command + # todo: we'll want to know know why the returncode was not 0 + call(obs_report_command, shell=True) + print 'jupyter nbconvert --to html --template full %s' % obs_report_command.split()[-1] + call('jupyter nbconvert --to html --template full %s' % obs_report_command.split()[-1], shell=True) + +parser = argparse.ArgumentParser(description='Runs the docker for beamformer folded observation report.') +parser.add_argument('directory', nargs=1) +parser.add_argument('--template', dest='template', default='/home/kat/software/katsdpscripts/reduction/obs_report/bf_fold_obs_report_template.ipynb') +args, unknown = parser.parse_known_args() + +if args.directory[0] == '': + raise RuntimeError('Please specify a directory from which files will be loaded.') +if not os.path.isfile(args.template): + raise IOError('Template not found at %s.' % args.template) +generate_ipynb_bf_fold_obs_report(args.directory[0], args.template) diff --git a/reduction/obs_report/bf_fold_obs_report_template.ipynb b/reduction/obs_report/bf_fold_obs_report_template.ipynb new file mode 100644 index 00000000..f256b8d4 --- /dev/null +++ b/reduction/obs_report/bf_fold_obs_report_template.ipynb @@ -0,0 +1,550 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# PTUSE observation report" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print 'Report Started:'\n", + "!date" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Load Useful functions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "run ptuse_report_useful_functions.py" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Files to be processed" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "directory = '/kat/isabella/timing/PTUSE_1_20180417074239_J0437m4715/'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Set up environment" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import sys\n", + "import numpy as np\n", + "from matplotlib import cm\n", + "from matplotlib import colors\n", + "import psrchive as psr\n", + "# from coast_guard import cleaners\n", + "import time\n", + "from astropy.table import Table\n", + "import subprocess" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%pylab inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Load data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# load and add archive files\n", + "archives = load_archive_data(path=directory, verbose=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Display metadata" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nbin = archives[0].get_nbin()\n", + "nchan = archives[0].get_nchan()\n", + "npol = archives[0].get_npol()\n", + "nsubint = archives[0].get_nsubint()\n", + "obs_type = archives[0].get_type()\n", + "telescope_name = archives[0].get_telescope()\n", + "source_name = archives[0].get_source()\n", + "ra = archives[0].get_coordinates().ra()\n", + "dec = archives[0].get_coordinates().dec()\n", + "centre_frequency = archives[0].get_centre_frequency()\n", + "bandwidth = archives[0].get_bandwidth()\n", + "dm = archives[0].get_dispersion_measure()\n", + "rm = archives[0].get_rotation_measure()\n", + "is_dedispersed = archives[0].get_dedispersed()\n", + "is_faraday_rotated = archives[0].get_faraday_corrected()\n", + "is_pol_calib = archives[0].get_poln_calibrated()\n", + "data_units = archives[0].get_scale()\n", + "data_state = archives[0].get_state()\n", + "obs_duration = archives[0].integration_length()\n", + "receiver_name = archives[0].get_receiver_name()\n", + "receptor_basis = archives[0].get_basis()\n", + "backend_name = archives[0].get_backend_name()\n", + "low_freq = archives[0].get_centre_frequency() - archives[0].get_bandwidth() / 2.0\n", + "high_freq = archives[0].get_centre_frequency() + archives[0].get_bandwidth() / 2.0\n", + "subint_duration = np.rint(archives[0].get_Integration(1).get_duration())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "print_metadata(archive=archives[0])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Clean archive, display RFI stats and RFI mask" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Run coast_guard RFI cleaning routines\n", + "clean_rfi = archives[0].clone()\n", + "clean_archive_surgical(clean_rfi, chan_numpieces=1, subint_numpieces=1, chan_threshold=3, subint_threshold=3)\n", + "clean_archive_rcvrstd(clean_rfi)\n", + "clean_archive_bandwagon(clean_rfi)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "weights = clean_rfi.get_weights()\n", + "plot_mask_channels(weights, low_freq, high_freq, nchan, nsubint, target=source_name)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print_zero_weight_statistics(weights, nchan, nsubint)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot subintegration stacks (Coherency data) & average profiles (Stokes data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "clean_profile_data = prepare_profile_data(clean_rfi)\n", + "dirty_profile_data = prepare_profile_data(archives[0])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dirty_stack_data = prepare_subint_data(archives=archives[0])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "clean_stack_data = prepare_subint_data(archives=clean_rfi)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax, clean_stack, clean_profile, unclean_stack, unclean_profile = create_subint_profile_figure(fig_title='Sub-intergration stack')\n", + "plot_stokes_profiles(clean_profile_data, source_name, ax=clean_profile, nbin=nbin, label='Clean data')\n", + "plot_stokes_profiles(dirty_profile_data, source_name, ax=unclean_profile, nbin=nbin, label='Unclean data')\n", + "plot_sub_intergration_stack(dirty_stack_data, ax=unclean_stack, title=False, nbin=nbin, nsubint=nsubint, obs_duration=obs_duration)\n", + "plot_sub_intergration_stack(clean_stack_data, ax=clean_stack, title=False, nbin=nbin, nsubint=nsubint, obs_duration=obs_duration)\n", + "set_phase_label(ax=clean_stack)\n", + "set_phase_label(ax=unclean_stack)\n", + "set_subint_time_label(ax=clean_stack, intergration=obs_duration)\n", + "set_subint_time_label(ax=unclean_stack, intergration=obs_duration)\n", + "clean_stack.set_ylabel('Subint number')\n", + "unclean_stack.set_ylabel('Subint number')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot phase vs. frequency image of flux (clean Coherency data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#coherency_flux_phase_data = prepare_coherency_flux_data(archives[0])\n", + "coherency_flux_phase_data = prepare_coherency_flux_data(clean_rfi)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 2, sharex='col', sharey='row', figsize=[15, 10], tight_layout='true')\n", + "fig.suptitle(source_name, fontsize='large', y=1.015)\n", + "fig.text(0.05, 1.0, 'Clean coherency data', fontsize='large')\n", + "ax[0, 0].set_ylabel('Frequency (MHz)')\n", + "ax[1, 0].set_ylabel('Frequency (MHz)')\n", + "ax[1, 0].set_xlabel('Pulse phase (bin)')\n", + "ax[1, 1].set_xlabel('Pulse phase (bin)')\n", + "plot_phase_freq_flux(coherency_flux_phase_data[0, :, :], ax=ax[0, 0], title='XX', low_freq=low_freq, high_freq=high_freq, nbin=nbin, nchan=nchan)\n", + "plot_phase_freq_flux(coherency_flux_phase_data[1, :, :], ax=ax[0, 1], title='YY', low_freq=low_freq, high_freq=high_freq, nbin=nbin, nchan=nchan)\n", + "plot_phase_freq_flux(coherency_flux_phase_data[2, :, :], ax=ax[1, 0], title='Re(XY)', low_freq=low_freq, high_freq=high_freq, nbin=nbin, nchan=nchan)\n", + "plot_phase_freq_flux(coherency_flux_phase_data[3, :, :], ax=ax[1, 1], title='Im(XY)', low_freq=low_freq, high_freq=high_freq, nbin=nbin, nchan=nchan)\n", + "set_phase_label(ax[1, 0])\n", + "set_phase_label(ax[1, 1])\n", + "set_channel_label(ax[1, 1], nchan)\n", + "set_channel_label(ax[0, 1], nchan)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot phase vs. frequency image of flux (clean Stokes data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "stokes_flux_phase_data = prepare_sokes_phase_freq(clean_rfi)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 2, sharex='col', sharey='row', figsize=[15, 10], tight_layout='true')\n", + "fig.suptitle(source_name, fontsize='large', y=1.01)\n", + "fig.text(0.05, 1.0, 'Clean stokes data', fontsize='large')\n", + "ax[0, 0].set_ylabel('Frequency (MHz)')\n", + "ax[1, 0].set_ylabel('Frequency (MHz)')\n", + "ax[1, 0].set_xlabel('Pulse phase (bin)')\n", + "ax[1, 1].set_xlabel('Pulse phase (bin)')\n", + "plot_phase_freq_flux(stokes_flux_phase_data[0, :, :], ax=ax[0, 0], title='I', low_freq=low_freq, high_freq=high_freq, nbin=nbin, nchan=nchan)\n", + "plot_phase_freq_flux(stokes_flux_phase_data[1, :, :], ax=ax[0, 1], title='Q', low_freq=low_freq, high_freq=high_freq, nbin=nbin, nchan=nchan)\n", + "plot_phase_freq_flux(stokes_flux_phase_data[2, :, :], ax=ax[1, 0], title='U', low_freq=low_freq, high_freq=high_freq, nbin=nbin, nchan=nchan)\n", + "plot_phase_freq_flux(stokes_flux_phase_data[3, :, :], ax=ax[1, 1], title='V', low_freq=low_freq, high_freq=high_freq, nbin=nbin, nchan=nchan)\n", + "set_phase_label(ax[1, 0])\n", + "set_phase_label(ax[1, 1])\n", + "set_channel_label(ax[1, 1], nchan)\n", + "set_channel_label(ax[0, 1], nchan)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot bandpass (unclean and clean Stokes data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dirty_mean_bandpass, dirty_variance_bandpass = get_mean_bandpass(archive=archives[0], stokes=False)\n", + "dirty_stokes_mean_bandpass, dirty_stokes_variance_bandpass = get_mean_bandpass(archive=archives[0], stokes=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "clean_mean_bandpass, clean_variance_bandpass = get_mean_bandpass(archive=clean_rfi, stokes=False)\n", + "clean_stokes_mean_bandpass, clean_stokes_variance_bandpass = get_mean_bandpass(archive=clean_rfi, stokes=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 1, sharex='col', figsize=[15, 10], tight_layout='true')\n", + "set_bandpass_plot_label(dirty_mean_bandpass, target=source_name, title='Unclean data', ax=ax[0], nchan=nchan, logscale=True)\n", + "if np.max(dirty_mean_bandpass[1, :]) > np.max(dirty_mean_bandpass[0, :]):\n", + " ax[0].plot(dirty_mean_bandpass[1, :], label='YY')\n", + " ax[0].plot(dirty_mean_bandpass[0, :], label='XX')\n", + "else:\n", + " ax[0].plot(dirty_mean_bandpass[0, :], label='XX')\n", + " ax[0].plot(dirty_mean_bandpass[1, :], label='YY')\n", + "ax[0].plot(dirty_stokes_mean_bandpass[0, :], label='Stokes I') \n", + "ax[0].fill_between(x=np.linspace(0, np.shape(dirty_mean_bandpass[0, :]), num=len(dirty_mean_bandpass[0, :]), endpoint=True),\n", + " y1=dirty_mean_bandpass[1, :].min(), y2=dirty_mean_bandpass[1, :].max(),\n", + " where=(clean_mean_bandpass[0, :] == 0.0), alpha=0.5, color='grey')\n", + "ax[0].set_ylim(dirty_mean_bandpass[1, :].min(), dirty_mean_bandpass[1, :].max())\n", + "ax[0].legend()\n", + "\n", + "# Plot the clean bandpass:\n", + "set_freq_label(ax[1], low_freq, high_freq, pos='xaxis')\n", + "set_bandpass_plot_label(clean_mean_bandpass, source_name, 'Clean data', ax[1], nchan, logscale=True)\n", + "clean_mean_bandpass[clean_mean_bandpass == 0.0] = np.nan\n", + "if np.max(clean_mean_bandpass[1, :]) > np.max(clean_mean_bandpass[0, :]):\n", + " ax[1].plot(clean_mean_bandpass[1, :], label='YY')\n", + " ax[1].plot(clean_mean_bandpass[0, :], label='XX')\n", + "else:\n", + " ax[1].plot(clean_mean_bandpass[0, :], label='XX')\n", + " ax[1].plot(clean_mean_bandpass[1, :], label='YY')\n", + "ax[1].set_ylim(np.nanmin(dirty_mean_bandpass[1, :]), np.nanmax(dirty_stokes_mean_bandpass[0, :])) \n", + "clean_stokes_mean_bandpass[clean_stokes_mean_bandpass == 0.0] = np.nan \n", + "ax[1].plot(clean_stokes_mean_bandpass[0, :], label='Stokes I') \n", + "ax[1].legend()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot dynamic baseline spectrum (unclean Coherency data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "unclean_mean_spectrum, unclean_varience_spectrum = get_baseline_spectrum(archives[0], npol, nchan, nsubint)\n", + "clean_mean_spectrum, clean_variance_spectrum = get_baseline_spectrum(clean_rfi, npol, nchan, nsubint, stokes=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 2, sharex='col', sharey='row', figsize=[15, 10], tight_layout='true')\n", + "fig.suptitle(source_name, fontsize='large', y=1.01)\n", + "fig.text(0.05, 1.0, 'Unclean data', fontsize='large')\n", + "\n", + "ax[0, 0].set_ylabel('Channel number')\n", + "ax[1, 0].set_ylabel('Channel number')\n", + "ax[1, 0].set_xlabel('Subint number')\n", + "ax[1, 1].set_xlabel('Subint number')\n", + "\n", + "mean_xx_min = unclean_mean_spectrum[:, 0, :].min()\n", + "mean_xx_max = unclean_mean_spectrum[:, 0, :].max()\n", + "mean_yy_min = unclean_mean_spectrum[:, 1, :].min()\n", + "mean_yy_max = unclean_mean_spectrum[:, 1, :].max()\n", + "mean_re_xy_min = unclean_mean_spectrum[:, 2, :].min()\n", + "mean_re_xy_max = unclean_mean_spectrum[:, 2, :].max()\n", + "mean_im_xy_min = unclean_mean_spectrum[:, 3, :].min()\n", + "mean_im_xy_max = unclean_mean_spectrum[:, 3, :].max()\n", + "\n", + "plot_dynamic_spectrum(unclean_mean_spectrum[:, 0, :], ax[0, 0], 'XX', low_freq, high_freq, nbin, nchan, \n", + " mean_xx_min, mean_xx_max, clean_data=False)\n", + "plot_dynamic_spectrum(unclean_mean_spectrum[:, 1, :], ax[0, 1], 'YY', low_freq, high_freq, nbin, nchan,\n", + " mean_yy_min, mean_yy_max, clean_data=False)\n", + "plot_dynamic_spectrum(unclean_mean_spectrum[:, 2, :], ax[1, 0], 'Re(XY)', low_freq, high_freq, nbin, nchan,\n", + " mean_re_xy_min, mean_re_xy_max, clean_data=False)\n", + "plot_dynamic_spectrum(unclean_mean_spectrum[:, 3, :], ax[1, 1], 'Im(XY)', low_freq, high_freq, nbin, nchan, \n", + " mean_im_xy_min, mean_im_xy_max, clean_data=False)\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot dynamic baseline spectrum (clean Stokes data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 2, sharex='col', sharey='row', figsize=[15, 10], tight_layout='true')\n", + "fig.suptitle(source_name, fontsize='large', y=1.01)\n", + "fig.text(0.05, 1.0, 'Clean data', fontsize='large')\n", + "\n", + "variance_xx_min = clean_variance_spectrum[:, 0, :].min()\n", + "variance_xx_max = clean_variance_spectrum[:, 0, :].max()\n", + "variance_yy_min = clean_variance_spectrum[:, 1, :].min()\n", + "variance_yy_max = clean_variance_spectrum[:, 1, :].max()\n", + "variance_re_xy_min = clean_variance_spectrum[:, 2, :].min()\n", + "variance_re_xy_max = clean_variance_spectrum[:, 2, :].max()\n", + "variance_im_xy_min = clean_variance_spectrum[:, 3, :].min()\n", + "variance_im_xy_max = clean_variance_spectrum[:, 3, :].max()\n", + "\n", + "plot_dynamic_spectrum(clean_variance_spectrum[:, 0, :], ax[0, 0], 'I', low_freq, high_freq, nbin, nchan, \n", + " variance_xx_min, variance_xx_max, clean_data=True)\n", + "plot_dynamic_spectrum(clean_variance_spectrum[:, 1, :], ax[0, 1], 'Q', low_freq, high_freq, nbin, nchan,\n", + " variance_yy_min, variance_yy_max, clean_data=True)\n", + "plot_dynamic_spectrum(clean_variance_spectrum[:, 2, :], ax[1, 0], 'U', low_freq, high_freq, nbin, nchan,\n", + " variance_re_xy_min, variance_re_xy_max, clean_data=True)\n", + "plot_dynamic_spectrum(clean_variance_spectrum[:, 3, :], ax[1, 1], 'V', low_freq, high_freq, nbin, nchan, \n", + " variance_im_xy_min, variance_im_xy_max, clean_data=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot signal-to-noise ratio vs. subintegration (clean Coherency data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "snr_data = get_subint_snr(clean_rfi, nsubint)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax1 = plt.subplots(1, 1, sharex='col', sharey='row', figsize=[15, 5], tight_layout='true')\n", + "fig.suptitle(source_name, fontsize='large', y=1.01)\n", + "fig.text(0.05, 1.0, 'Clean data', fontsize='large')\n", + "ax1.set_xlabel('Subint number')\n", + "ax1.set_ylabel('Signal-to-noise ratio')\n", + "ax1.set_xlim(0, nsubint)\n", + "ax1_secondary = ax1.twiny()\n", + "ax1_secondary.set_frame_on(True)\n", + "ax1_secondary.patch.set_visible(False)\n", + "ax1_secondary.xaxis.set_ticks_position('bottom')\n", + "ax1_secondary.set_xlabel('Time (seconds)')\n", + "ax1_secondary.xaxis.set_label_position('bottom')\n", + "ax1_secondary.spines['bottom'].set_position(('outward', 50))\n", + "ax1_secondary.set_xlim(0, obs_duration)\n", + "ax1.plot(snr_data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print 'Report Ended:'\n", + "!date" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.12" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/reduction/obs_report/ptuse_report_useful_functions.py b/reduction/obs_report/ptuse_report_useful_functions.py new file mode 100644 index 00000000..ee2a925a --- /dev/null +++ b/reduction/obs_report/ptuse_report_useful_functions.py @@ -0,0 +1,504 @@ +#!/usr/bin/env python +import numpy as np +import matplotlib.pyplot as plt +from astropy.table import Table +from matplotlib import cm +from matplotlib import colors +from matplotlib.gridspec import GridSpec +import os +import psrchive as psr +from coast_guard import cleaners + +def print_metadata(archive): + """Function to print archive file metadata. + + Input: + archive: PSRCHIVE Archive object. + + Output: + print metadata in a nice table. + """ + # get metadata from file + nbin = archive.get_nbin() + nchan = archive.get_nchan() + npol = archive.get_npol() + nsubint = archive.get_nsubint() + obs_type = archive.get_type() + telescope_name = archive.get_telescope() + source_name = archive.get_source() + ra = archive.get_coordinates().ra() + dec = archive.get_coordinates().dec() + centre_frequency = archive.get_centre_frequency() + bandwidth = archive.get_bandwidth() + dm = archive.get_dispersion_measure() + rm = archive.get_rotation_measure() + is_dedispersed = archive.get_dedispersed() + is_faraday_rotated = archive.get_faraday_corrected() + is_pol_calib = archive.get_poln_calibrated() + data_units = archive.get_scale() + data_state = archive.get_state() + obs_duration = archive.integration_length() + receiver_name = archive.get_receiver_name() + receptor_basis = archive.get_basis() + backend_name = archive.get_backend_name() + low_freq = archive.get_centre_frequency() - archive.get_bandwidth() / 2.0 + high_freq = archive.get_centre_frequency() + archive.get_bandwidth() / 2.0 + subint_duration = np.rint(archive.get_Integration(1).get_duration()) + subint_duration = np.rint(archive.get_Integration(1).get_duration()) + # Print out metadata + print '=====================================================================================' + print 'Attribute Names Description Value ' + print '=====================================================================================' + print 'nbin Number of pulse phase bins %s' % nbin + print 'nchan Number of frequency channels %s' % nchan + print 'npol Number of polarizations %s' % npol + print 'nsubint Number of sub-integrations %s' % nsubint + print 'type Observation type %s' % obs_type + print 'site Telescope name %s' % telescope_name + print 'name Source name %s' % source_name + print 'coord Source coordinates %s%s' % (ra.getHMS(), dec.getDMS()) + print 'freq Centre frequency (MHz) %s' % centre_frequency + print 'bw Bandwidth (MHz) %s' % bandwidth + print 'dm Dispersion measure (pc/cm^3) %s' % dm + print 'rm Rotation measure (rad/m^2) %s' % rm + print 'dmc Dispersion corrected %s' % is_dedispersed + print 'rmc Faraday Rotation corrected %s' % is_faraday_rotated + print 'polc Polarization calibrated %s' % is_pol_calib + print 'scale Data units %s' % data_units + print 'state Data state %s' % data_state + print 'length Observation duration (s) %s' % obs_duration + print 'rcvr:name Receiver name %s' % receiver_name + print 'rcvr:basis Basis of receptors %s' % receptor_basis + print 'be:name Name of the backend instrument %s' % backend_name + +def load_archive_data(path, verbose=False): + """Function to load .ar files and convert to PSRCHIVE archive objects. + + Input: + path : full path to location of the .ar files. + verbose : option to run in verbose mode (default=False) + + Output: + archives : list of PSRCHIVE archive objects + """ + files = [] + for file in os.listdir(path): + if file.endswith('.ar'): + files.append(file) + files.sort() + archives = [] + archives = [psr.Archive_load(path + file) for file in files] + if verbose: + print '======================================================================================================' + print ' Files to be processed: ' + print '======================================================================================================' + for i in range(1, len(archives)): + archives[0].append(archives[i]) # add the .ar files (added file is archive[0]) + if verbose: + print archives[i] + return archives + +def clean_archive_bandwagon(archive_clone, bad_chan_tol=0.9, bad_sub_tol=1.0): + """Function to clean the archive files using coast_guard cleaner bandwagon. + + Input: + archive_clone: a clone of the PSRCHIVE archive object to clean. + bad_chan_tol: fraction of bad channels to be tolarated before mask (float between 0 - 1). + bad_sub_tol: fraction of bad sub-intergrations to be tolerated before mask (float between 0 - 1) + + """ + cleaner3 = cleaners.load_cleaner('bandwagon') + cleaner3.parse_config_string('badchantol=0.99,badsubtol=1.0') + cleaner3.run(archive_clone) + + +def clean_archive_rcvrstd(archive_clone, bad_channels='0:210;3896:4095', bad_frequencies=None, bad_subints=None, \ + trim_bw=0, trim_frac=0, trim_num=0): + """Function to clean the archive files using coast_guard cleaner rcvrstd. + + Input: + archive_clone: a clone of the PSRCHIVE archive object to clean. + bad_channels: bad channels to de-weight (default: band edges (0:210, 3896:4095). + bad_frequencies: bad frequencies to de-weight (default: None). + bad_subints: bad sub-ints to de-weight (default: None). + trim_bw: bandwidth of each band-edge (in MHz) to de-weight (default: None). + trim_frac: fraction of each band-edge to de-weight, float between 0 - 0.5 (default: None). + trim_num: number of channels to de-weight at each edge of the band (default: None). + """ + cleaner2 = cleaners.load_cleaner('rcvrstd') + + cleaner2.parse_config_string('badchans=%s,badfreqs=%s,badsubints=%s,trimbw=%s,trimfrac=%s,trimnum=%s' \ + %(str(bad_channels), str(bad_frequencies), str(bad_subints), \ + str(trim_bw), str(trim_frac), str(trim_num))) + cleaner2.run(archive_clone) + + +def clean_archive_surgical(archive_clone, chan_threshold=3, subint_threshold=3, chan_numpieces=1, subint_numpieces=1): + """Function to clean the archive files using coast_guard cleaner surgical. + + Input: + archive_clone: a clone of the PSRCHIVE archive object. + chan_threshold: threshold sigma of a profile in a channel. + subint_threshold: threshold sigma of a profile in a sub-intergration. + chan_numpieces: the number of equally sized peices in each channel (used for detranding in surgical) + subint_numpieces: the number of equally sized peices in each sub-int (used for detranding in surgical) + """ + cleaner = cleaners.load_cleaner('surgical') + cleaner.parse_config_string('chan_numpieces=%s,subint_numpieces=%s,chanthresh=%s,subintthresh=%s'\ + % (str(chan_numpieces), str(subint_numpieces), str(chan_threshold),\ + str(subint_threshold))) + cleaner.run(archive_clone) + + +def print_zero_weight_statistics(weights, nchan, nsubint): + """Function to determine a fraction of channels set to zero. + Input: + weights: rfi weights from coast_guard cleaned files. + nchans: number of channels in the band. + nsubint: number of sub-intergrations. + """ + counter = 0 + zero_weight_channel = 0 + for i in range(nsubint): + zero_weight_channel = 0 + for j in range(nchan): + if weights[i][j] == 0.0: + counter += 1 + zero_weight_channel += 1 # channels set with zero weights + subint_proc = (float(zero_weight_channel) / float(nchan)) * 100 + print 'Subint %s has %s channels (%.2f%%) set to zero.' % (i, zero_weight_channel, subint_proc) + percentage = (float(counter) / float(weights.size)) * 100 + print '%s data points out of %s with weights set to zero.' % (counter, weights.size) + print '%.2f%% of the data was set to zero.' % (percentage) + +def prepare_subint_data(archives): + """Funtion to prepare PSRCHIVE archive file for plotting. + + Input: + archives: PSRCHIVE archive file object + Returns: + data: dedispersed and freq-pol scrunched data. + """ + stack = archives.clone() + stack.dedisperse() + stack.remove_baseline() + stack.centre_max_bin() + stack.fscrunch() + stack.pscrunch() + data = stack.get_data().squeeze() + return data + +def prepare_profile_data(archive): + """Function to prepare average profile for stokes data. + + Input: + archive: PSRCHIVE archive file object. + + Returns: + profile_data: time and frequency scrunched profiles for stokes (I, Q, U, V). + """ + profile = archive.clone() + profile.dedisperse() + profile.remove_baseline() + profile.centre_max_bin() + profile.tscrunch() + profile.fscrunch() + profile.convert_state('Stokes') + profile_data = profile.get_data() + profile_data = profile_data.squeeze() + return profile_data + +def prepare_sokes_phase_freq(archive): + """Function to prepare average profile for stokes data. + + Input: + archive: PSRCHIVE archive file object. + + Returns: + stokes_data: time-freq flux for stokes (I, Q, U, V). + """ + stokes_flux = archive.clone() + stokes_flux.dedisperse() + stokes_flux.remove_baseline() + stokes_flux.centre_max_bin() + stokes_flux.tscrunch() + stokes_flux.convert_state('Stokes') + stokes_data = stokes_flux.get_data().squeeze() + return stokes_data + +def prepare_coherency_flux_data(archive): + """Function to prepare coherancy data. + + Input: + archive: PSRCHIVE archive file object. + + Returns: + coherency_flux_phase_data: time scrunched coherency flux data. + """ + coherency_flux_phase = archive.clone() + coherency_flux_phase.dedisperse() + coherency_flux_phase.remove_baseline() + coherency_flux_phase.centre_max_bin() + coherency_flux_phase.tscrunch() + coherency_flux_phase_data = coherency_flux_phase.get_data().squeeze() + return coherency_flux_phase_data + +def get_mean_bandpass(archive, stokes=False): + """Function to detetermine the mean bandpass. + + Input: + archive: PSRCHIVE archive file object. + stokes: option to get bandpass for stokes I, Q, U, V + Returns: + bandpass_mean: mean bandpass + bandpass_varience: variance bandpass + """ + bandpass = archive.clone() + if stokes: + bandpass.convert_state('Stokes') + bandpass.tscrunch() + (bandpass_mean, bandpass_variance) = bandpass.get_Integration(0).baseline_stats() + return bandpass_mean, bandpass_variance + + + +def get_baseline_spectrum(archive, npol, nchan, nsubint, stokes=False): + """Function to determine teh baseline spectrum + """ + dynamic_spectrum = archive.clone() + if stokes: + dynamic_spectrum.convert_state('Stokes') + mean = np.zeros((nsubint, npol, nchan)) + variance = np.zeros((nsubint, npol, nchan)) + for subint in range(nsubint): + m, v = dynamic_spectrum.get_Integration(subint).baseline_stats() + mean[subint] = m[:, :] + variance[subint] = v[:, :] + return mean, variance + + + +def get_subint_snr(archive, nsubint): + """Function to determine the subintergration S/N + Input: + archive: PSRCHIVE archive object + nsubint: number of sub-intergrations + + Return: + snr_data: s/n per sub-intergration + """ + archive_snr = archive.clone() + new_integration = archive_snr.get_Integration(0).total() + snr_data = np.zeros(nsubint) + for i_subintegration in range (1, nsubint): + next_integration = archive_snr.get_Integration(i_subintegration).total() + new_integration.combine(next_integration) + profile = new_integration.get_Profile(0, 0) + snr_data[i_subintegration] = profile.snr() + return snr_data + + +def create_subint_profile_figure(fig_title): + """Function to create a figure and axis. To be used to compare + clean and unclean subintergartion stack, profile, and bandpass. + + Input: + fig_tittle: title for the whole figure. + + Returns: + fig: the figure handle. + ax: figure axes handle. + clean_stack: subplot axes handle to plot clean subintergration stack. + clean_profile: subplot axes handle to plot clean stokes average profile. + unclean_stack: subplot axes handle to plot dirty subintergration stack. + unclean_profile: subplot axes handle to plot dirty stokes average profile. + """ + fig, ax = plt.subplots(1, 2, figsize=(20, 15), sharex='col', tight_layout='True') + # create 2 sub-figures + grid_spec_1 = GridSpec(3, 3) + grid_spec_1.update(left=0.05, right=0.48, wspace=0.05, hspace=0.1) + unclean_profile = plt.subplot(grid_spec_1[0, :]) + unclean_stack = plt.subplot(grid_spec_1[1:, :]) + unclean_profile.tick_params(labelbottom=False) + + + grid_spec_2 = GridSpec(3, 3) + grid_spec_2.update(left=0.6, right=0.98, wspace=0, hspace=0.1) + clean_profile = plt.subplot(grid_spec_2[0, :]) + clean_profile.tick_params(labelbottom=False, labelleft=False) + clean_stack = plt.subplot(grid_spec_2[1:, :]) + clean_stack.tick_params(labelleft=False) + + fig.suptitle(fig_title, fontsize=18) + return fig, ax, clean_stack, clean_profile, unclean_stack, unclean_profile + + +def plot_mask_channels(weights, low_freq, high_freq, nchan, nsubint, target): + """Function to plot image of the rfi masked channels in the band. + Input: + weights: rfi weights from coast_guard cleaned files. + low_freq: minimum frequency in the band (MHz). + high_freq: maximum frequency in the band (MHz). + nchans: number of channels. + nsubint: number of sub-intergrations. + target: name of the target source observed. + """ + fig, ax1 = plt.subplots(1, 1, figsize = [15, 10]) + ax1.set_title(target, fontsize=20) + ax1.set_title('RFI mask ', loc='left', fontsize=20) + ax1.set_ylabel('Channel number', fontsize=18) + ax1.yaxis.set_ticks(np.arange(0, nchan - 1, 200)) + ax1.set_xlabel('Subint number', fontsize=18) + ax1_secondary = ax1.twinx() + ax1_secondary.set_ylabel('Frequency (MHz)', fontsize=18) + ax1_secondary.set_ylim(low_freq, high_freq) + ax1_secondary.yaxis.set_ticks(np.arange(np.rint(low_freq), np.rint(high_freq), 25)) + ax1.imshow(weights.T, origin='lower', aspect='auto', cmap=colors.ListedColormap(['red', 'white']), \ + interpolation='none', extent=(0, nsubint - 1, 0, nchan - 1)) + +def plot_sub_intergration_stack(data, ax, nbin, nsubint, obs_duration, title=False): + """Function to plot the prepared sub intergration data on provided axis. + + Input: + data: data to be plotted + ax: axes to plot + title: title of the plot + nbin: number of bins for the plot + nsubint: number of sub-intergrations + obs_duration: observation intergration time + """ + if title: + ax.set_title(title, fontsize=14) + ax.set_xlabel('Pulse phase (bin)', fontsize=14) + ax.xaxis.set_ticks(np.arange(0, nbin - 1, 100)) + ax.imshow(data, cmap=cm.afmhot, aspect='auto', vmax=data.max() * 0.75, + interpolation='none', extent=(0, nbin - 1, 0, nsubint - 1)) + +def plot_stokes_profiles(profile_data, target, ax, nbin, label='Clean data'): + """ Function to plot the stokes profiles. + + Input: + profile_data: profiles to plot shape (I, Q, U, V) + target: target source name + ax: figure axis to plot on. + nbin: number of bins. + label: plot title label. + """ + min_max_vals = np.array([profile_data[0, :].max(), profile_data[1, :].max(), + profile_data[2, :].max(), profile_data[3, :].max(), + profile_data[0, :].min(), profile_data[1, :].min(), + profile_data[2, :].min(), profile_data[3, :].min()]) + ax.text(20, profile_data.max(), label, fontsize='large') + ax.set_title(target) + ax.set_ylabel('Flux (a.u.)') + ax.set_ylim((min_max_vals.min() * 1.2), min_max_vals.max() + 0.1 * min_max_vals.max()) + ax.set_xlim(0, nbin - 1) + ax.plot(profile_data[0, :], 'k', label='I') + ax.plot(profile_data[1, :], 'r', label='Q') + ax.plot(profile_data[2, :], 'g', label='U') + ax.plot(profile_data[3, :], 'b', label='V') + ax.xaxis.set_ticks(np.arange(0, nbin - 1, 50)) + ax.legend() + +def plot_phase_freq_flux(data, ax, title, low_freq, high_freq, nbin, nchan): + """Function to plot the phase-frequency image of the flux. + Input: + data: time-frequency data to plot + ax: axis to plot. + + """ + ax.set_title(title) + ax.xaxis.set_ticks(np.arange(0, nbin - 1, 100)) + ax.yaxis.set_ticks(np.arange(0, nchan - 1, 50)) + ax.imshow(data, cmap=cm.afmhot, aspect='auto', interpolation='none', extent=(0, nbin - 1.0, low_freq, high_freq)) + + + +def plot_dynamic_spectrum(data, ax, title, low_freq, high_freq, nbin, nchan, minimum, maximum, clean_data=True): + """Function to plot the phase-frequency image of the flux. + Input: + data: time-frequency data to plot + ax: axis to plot. + """ + ax.set_title(title) + ax.yaxis.set_ticks(np.arange(0, nchan - 1, 400)) + if clean_data: + ax.imshow(data.T, origin='lower', aspect='auto', cmap=cm.afmhot, + norm=colors.Normalize(vmin=minimum, vmax=maximum), + interpolation='none') + else: + ax.imshow(data.T, origin='lower', aspect='auto', cmap=cm.afmhot, + norm=colors.SymLogNorm(10, linscale=1.0, vmin=minimum, vmax=maximum), + interpolation='none') + + + +def set_subint_time_label(ax, intergration): + """Sets twin axes sharing the axis + + Input: + ax: axis to be shared + intergration: observation duration + """ + ax_tertiary = ax.twinx() + ax_tertiary.set_ylabel('Time (seconds)') + ax_tertiary.yaxis.set_label_position('right') + ax_tertiary.set_ylim(0, intergration) + +def set_phase_label(ax): + """Sets twin axes sharing the axis + + Input: + ax : axis to be shared + """ + ax_secondary = ax.twiny() + ax_secondary.set_frame_on(True) + ax_secondary.patch.set_visible(False) + ax_secondary.xaxis.set_ticks_position('bottom') + ax_secondary.set_xlabel('Pulse phase (degrees)') + ax_secondary.xaxis.set_label_position('bottom') + ax_secondary.spines['bottom'].set_position(('outward', 50)) + ax_secondary.set_xlim(0, 360) + _ = ax_secondary.xaxis.set_ticks(np.arange(0, 360, 20)) + +def set_channel_label(ax, nchan): + """Sets twin axes sharing the axis + + Input: + ax: axis to be shared + nchan: number of channels + """ + ax_secondary = ax.twinx() + ax_secondary.set_ylabel('Channel number') + ax_secondary.set_ylim(0, nchan - 1) + ax_secondary.yaxis.set_ticks(np.arange(0, nchan - 1 , 200)) + +def set_freq_label(ax, low_freq, high_freq, pos='xaxis'): + """ + """ + if pos == 'xaxis': + ax_secondary = ax.twiny() + elif pos == 'yaxis': + ax_secondary = ax.twinx() + ax_secondary.set_frame_on(True) + ax_secondary.patch.set_visible(False) + ax_secondary.xaxis.set_ticks_position('bottom') + ax_secondary.set_xlabel('Frequency (MHz)') + ax_secondary.xaxis.set_label_position('bottom') + ax_secondary.spines['bottom'].set_position(('outward', 50)) + ax_secondary.set_xlim(low_freq, high_freq) + _ = ax_secondary.xaxis.set_ticks(np.arange(np.rint(low_freq), np.rint(high_freq), 50)) + + +def set_bandpass_plot_label(bandpass, target, title, ax, nchan, logscale=True): + """Function to prepare bandpass plot. + """ + bandpass_max_min_vals = np.array([bandpass[0, :].min(), bandpass[0, :].max(), + bandpass[1, :].min(), bandpass[1, :].max()]) + ax.set_title(target) + ax.set_title(title, loc='left') + ax.set_xlim(0, nchan - 1) + ax.set_ylabel('Flux (a.u.)') + ax.xaxis.set_ticks(np.arange(0, nchan - 1, 200)) + if logscale: + ax.set_yscale('log')