day_survey_visualization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bc08fb2b",
   "metadata": {},
   "source": [
    "# Do daily visual checks of data quality\n",
    "## L. Caffarello, July 2022\n",
    "### Updates by A. Ordog, August-September 2022\n",
    "### August 25, 2022 (AO): \n",
    "#### - changed filenames to account for survey phases\n",
    "#### - changed frequency channels read in to all channels (instead of every 12th)\n",
    "#### - included 5 MHz to each side of the central frequency in the single-scan, 1D timeseries plots\n",
    "#### - included persistent RFI mask indication (grey background) on timeseries plots\n",
    "### September 12, 2022 (AO):\n",
    "#### - option to mask out persistent RFI\n",
    "#### - indicate elevation of scans\n",
    "#### - scan IDs indicated in waterfall plot columns\n",
    "### January 5, 2023 (AO):\n",
    "#### - option to plot scan in 2D in Galactic coordinates\n",
    "### January 11, 2023 (A):\n",
    "#### - fixed formatting for reading in phase and day (only need to select it in first cell)\n",
    "#### - included azimuth axis on 1D single-scan time-domain plots\n",
    "#### - added option to change frequency (default 520 MHz) and colorscale in 2D Galactic plot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1c606e0",
   "metadata": {},
   "source": [
    "## Import packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1cb01fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import imp\n",
    "import os\n",
    "import subprocess\n",
    "import h5py\n",
    "import numpy as np\n",
    "from astropy.time import Time\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.interpolate import griddata\n",
    "import datetime\n",
    "import matplotlib.dates as mdates\n",
    "from matplotlib.dates import HourLocator as HourLocator\n",
    "from matplotlib.dates import MinuteLocator as MinuteLocator\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from astropy import units as u\n",
    "from astropy.time import TimeDelta\n",
    "from ipywidgets import interact\n",
    "from ipywidgets import fixed\n",
    "from ipywidgets import interactive_output\n",
    "from astropy.coordinates import SkyCoord,EarthLocation, AltAz, ICRS, Galactic, FK4, FK5\n",
    "\n",
    "#### Change the directory to where the files are located\" ####\n",
    "day = 10\n",
    "phase = 3\n",
    "directory = '/srv/data/dva/survey_azimuth_scans/'\n",
    "#directory = '/srv/data/dva/survey_azimuth_scans/day_45_lightning/'\n",
    "#directory = '../DVA/Data_Files/DVA_Day_Surveys/'\n",
    "#directory = '/media/ordoga/15m_band1_survey/dva_phase1/survey_phase1_day'+day+'/'\n",
    "##############################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07126407",
   "metadata": {},
   "source": [
    "## Read in the file listing azimuth scan start and stop times:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba279f8a",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "scan_id = []    # The scan id number\n",
    "scan_start = []  # Start time of the scan (UTC)\n",
    "scan_stop = []   # Stop time of the scan (UTC)\n",
    "scan_el = []\n",
    "\n",
    "# Read in the data and store it in arrays:\n",
    "with open(directory+'DVAsurvey_phase'+str(phase)+'_day'+f\"{day:03}\"+'.txt') as fp:\n",
    "    for line in fp:       \n",
    "        scan_id.append(int(line.split()[0]))\n",
    "        scan_start.append(line.split()[1]+'T'+line.split()[2][0:12])\n",
    "        scan_stop.append(line.split()[3]+'T'+line.split()[4][0:12])\n",
    "        scan_el.append(line.split()[5][0:2])\n",
    "        \n",
    "# Print out the scan numbers with their start and stop times:\n",
    "for i in range(0,len(scan_id)):\n",
    "    print(f\"{scan_id[i]:04}\",scan_start[i],scan_stop[i],scan_el[i])\n",
    "\n",
    "# Convert start and stop times to Modified Julian Day (MJD).\n",
    "# This is needed for plotting and for selecting out data collected\n",
    "# between particular times:\n",
    "scan_start_mjd = Time(scan_start, format='isot',scale='utc').mjd\n",
    "scan_stop_mjd  = Time(scan_stop,  format='isot',scale='utc').mjd\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47786e9c",
   "metadata": {},
   "source": [
    "## Read in scan files and stitch them together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b36dc5f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "t_set = []\n",
    "az_set = []\n",
    "dec_set = []\n",
    "ra_set = []\n",
    "el_set = []\n",
    "noise_set = []\n",
    "trim_flag = []\n",
    "\n",
    "scan0 = f\"{scan_id[0]:04}\"\n",
    "\n",
    "freq_channel_increment = 1 #TODO: I'll have to change this to 1 once I am sure I'm reading the data correctly\n",
    "\n",
    "# Use one of the scans to get the list of frequencies:\n",
    "file = h5py.File(directory+'dva_survey_phase'+str(phase)+'_raw_'+scan0+'.h5','r')\n",
    "freq = file['data']['beam_0']['band_SB0']['frequency'][::freq_channel_increment]/1e6\n",
    "\n",
    "# Create empty arrays for the power data:\n",
    "RR_set = np.empty([0,len(freq)])\n",
    "LL_set = np.empty([0,len(freq)])\n",
    "reRL_set = np.empty([0,len(freq)])\n",
    "imRL_set = np.empty([0,len(freq)])\n",
    "\n",
    "# Loop through all the scans in the \"scan_num\" list:\n",
    "for i in scan_id:\n",
    "#for i in scan_id[0:5]:\n",
    "    print(i)\n",
    "    # select the file:\n",
    "    file = h5py.File(directory+'dva_survey_phase'+str(phase)+'_raw_'+f\"{i:04}\"+'.h5','r')\n",
    "    print(file)\n",
    "    \n",
    "    # access the correct location in the file structure:\n",
    "    dataset = file['data']['beam_0']['band_SB0']['scan_0']\n",
    "    \n",
    "    # Add the position and time data to the corresponding arrays:\n",
    "    dec_set = np.concatenate([dec_set,dataset['metadata']['declination']])\n",
    "    ra_set = np.concatenate([ra_set,dataset['metadata']['right_ascension']])\n",
    "    el_set = np.concatenate([el_set,dataset['metadata']['elevation']])\n",
    "    az_set = np.concatenate([az_set,dataset['metadata']['azimuth']])\n",
    "    t_set = np.concatenate([t_set,dataset['metadata']['utc']])\n",
    "    noise_set = np.concatenate([noise_set,dataset['metadata']['noise_state']]) #This is a \"mask\" for noise regions 1 = noise 0=all good\n",
    "    trim_flag = np.concatenate([trim_flag,dataset['metadata']['trim_scan_flag']])\n",
    "    \n",
    "    # Add the spectrometer power data to the corresponding arrays:\n",
    "    RR_set = np.concatenate([RR_set,dataset['data'][:,0,::freq_channel_increment]],axis=0)\n",
    "    LL_set = np.concatenate([LL_set,dataset['data'][:,1,::freq_channel_increment]],axis=0)\n",
    "    reRL_set = np.concatenate([reRL_set,dataset['data'][:,2,::freq_channel_increment]],axis=0)\n",
    "    imRL_set = np.concatenate([imRL_set,dataset['data'][:,3,::freq_channel_increment]],axis=0)\n",
    "    \n",
    "t_plt = Time(t_set, format='isot',scale='utc').mjd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41f5a3c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(t_set)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba424964",
   "metadata": {},
   "source": [
    "## Read in persistent RFI mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3562669",
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 0\n",
    "RFI_mask_idx = []\n",
    "with open('/srv/data/dva/RFIpersist_mask/RFIpersist_mask.txt') as fp:\n",
    "#with open('../DVA2/DATA/PersistRFI_v01_Jul13.txt') as fp:\n",
    "    for line in fp:\n",
    "        if i>0: \n",
    "            #print(line)\n",
    "            RFI_mask_idx.append(int(line.split()[0]))\n",
    "        i=i+1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fa8f724",
   "metadata": {},
   "source": [
    "## Updated visualization, including single-scan option\n",
    "### Changes made by AO:\n",
    "#### - power bounds in dB to see lower power signal\n",
    "#### - changed to log scale and added RR to time series plots and spectra\n",
    "#### - added grids to 1D plots\n",
    "#### - added option for single-scan plots (both waterfall and 1D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4648ac28",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = freq[1]-freq[0]\n",
    "power_min = 65\n",
    "power_max = 79\n",
    "fs = 14\n",
    "\n",
    "def DVA_Waterfall_View(RFI_masked):  \n",
    "    \n",
    "    #TODO: implement toggling between RR, LL,reRL, imRL\n",
    "       \n",
    "    fig,axs = plt.subplots(1,1,figsize=(15,10)) \n",
    "    \n",
    "    for i in range(0,len(scan_id)):\n",
    "        w = np.where((t_plt>=scan_start_mjd[i]) & (t_plt<=scan_stop_mjd[i]))[0]\n",
    "        extent = [scan_start_mjd[i],scan_stop_mjd[i],freq[0],freq[-1]]\n",
    "        \n",
    "        if scan_el[i] == '49':\n",
    "            textclr = 'blue'\n",
    "        else:\n",
    "            textclr = 'red'\n",
    "                \n",
    "        data_plot = 10*np.log10(LL_set[w,:].T)\n",
    "        if RFI_masked:\n",
    "            data_plot[RFI_mask_idx,:] = np.nan\n",
    "    \n",
    "        im = axs.imshow(data_plot,aspect='auto',vmin=power_min,vmax=power_max,\n",
    "                        origin='lower',extent=extent,cmap='viridis')\n",
    "\n",
    "        axs.text(scan_start_mjd[i]+3e-3,1038,f\"{scan_id[i]:04}\",rotation=45,fontsize=fs,color=textclr)\n",
    "    \n",
    "    divider = make_axes_locatable(axs)\n",
    "    cax = divider.append_axes('right', size='2%', pad=0.05)\n",
    "    cbar = fig.colorbar(im, cax=cax, orientation='vertical')\n",
    "    cbar.ax.tick_params(labelsize=fs) \n",
    "    cbar.set_label('Power (dB)', fontsize=fs)\n",
    "\n",
    "    axs.set_xlim(t_plt[0],t_plt[-1])\n",
    "    axs.set_ylim(freq[0],freq[-1])\n",
    "    axs.tick_params(axis='both', labelsize=fs)\n",
    "    axs.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "    axs.fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "    axs.set_xlabel('Time (UTC)',fontsize=fs)        \n",
    "    axs.set_ylabel('Frequency (MHz)',fontsize=fs)\n",
    "    \n",
    "    return\n",
    "\n",
    "def DVA_Waterfall_View_singlescan(scan_chosen, RFI_masked):\n",
    "    \n",
    "    #print(axmap)\n",
    "    #TODO: implement toggling between RR, LL,reRL, imRL\n",
    "    scan_id_plot = scan_chosen\n",
    "    scan_idx = np.where(np.array(scan_id) == scan_id_plot)    \n",
    "    w = np.where((t_plt>=scan_start_mjd[scan_idx]) & (t_plt<=scan_stop_mjd[scan_idx]))[0]\n",
    "        \n",
    "    fig,axs = plt.subplots(1,1,figsize=(15,10))     \n",
    "    extent = [scan_start_mjd[scan_idx][0],scan_stop_mjd[scan_idx][0],freq[0],freq[-1]]\n",
    "    \n",
    "    data_plot = 10*np.log10(LL_set[w,:].T)\n",
    "    if RFI_masked:\n",
    "        data_plot[RFI_mask_idx,:] = np.nan\n",
    "    \n",
    "    im = axs.imshow(data_plot,aspect='auto',vmin=power_min,vmax=power_max,\n",
    "                    origin='lower',extent=extent,cmap='viridis')\n",
    "\n",
    "    divider = make_axes_locatable(axs)\n",
    "    cax = divider.append_axes('right', size='2%', pad=0.05)\n",
    "    cbar = fig.colorbar(im, cax=cax, orientation='vertical')\n",
    "    cbar.ax.tick_params(labelsize=fs) \n",
    "    cbar.set_label('Power (dB)', fontsize=fs)\n",
    "\n",
    "    axs.set_xlim(extent[0],extent[1])\n",
    "    axs.set_ylim(freq[0],freq[-1])\n",
    "    axs.tick_params(axis='both', labelsize=fs)\n",
    "    axs.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "    axs.fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "    axs.set_xlabel('Time (UTC)',fontsize=fs)        \n",
    "    axs.set_ylabel('Frequency (MHz)',fontsize=fs)\n",
    "    axs.set_title('Scan '+str(scan_id_plot)+' , el = '+str(scan_el[scan_idx[0][0]]),fontsize=fs)\n",
    "    \n",
    "    return\n",
    "\n",
    "\n",
    "def DVA_Cross_Sections(freq_chosen, time_idx, freq_cross_section,RFI_masked):\n",
    "    \n",
    "    freq_measured = np.where(abs(freq-freq_chosen)<df)[0][0]\n",
    "    freq_below = np.where((freq>freq_chosen-5) & (freq<freq_chosen))[0]\n",
    "    freq_above = np.where((freq<freq_chosen+5) & (freq>freq_chosen))[0]\n",
    "   \n",
    "    fig,axs1 = plt.subplots(1,1,figsize=(16,6))  \n",
    "\n",
    "    if(freq_cross_section):\n",
    "        \n",
    "        data_plot_L = 10*np.log10(LL_set[time_idx,:])\n",
    "        data_plot_R = 10*np.log10(RR_set[time_idx,:])\n",
    "        \n",
    "        if RFI_masked:\n",
    "            data_plot_L[RFI_mask_idx] = np.nan\n",
    "            data_plot_R[RFI_mask_idx] = np.nan\n",
    "\n",
    "        axs1.plot(freq,data_plot_L,label='LL',color='blue')\n",
    "        axs1.plot(freq,data_plot_R,label='RR',color='red')        \n",
    "        axs1.vlines(freq_chosen, 0 , 100e9, color = 'purple')\n",
    "        axs1.set_ylim(power_min,power_max) \n",
    "        axs1.set_xlim(350,1050)\n",
    "        axs1.set_xlabel('Frequency',fontsize=fs)\n",
    "        axs1.tick_params(axis='both',labelsize=fs)\n",
    "    else:\n",
    "        axs1.scatter(t_plt, 10*np.log10(LL_set[:,freq_measured]),label='LL',s=0.8,zorder=1,color='blue')\n",
    "        axs1.scatter(t_plt, 10*np.log10(RR_set[:,freq_measured]),label='RR',s=0.8,zorder=1,color='red')\n",
    "          \n",
    "        axs1.vlines(t_plt[time_idx], 0 , 100e9, color = 'purple')\n",
    "        \n",
    "        for i in range(0,len(scan_id)):\n",
    "            if scan_el[i] == '49':\n",
    "                elclr = 'C0'\n",
    "            else:\n",
    "                elclr = 'C1'        \n",
    "            axs1.axvspan(scan_start_mjd[i],scan_stop_mjd[i],color=elclr,alpha=0.2,zorder=0)\n",
    "            axs1.text(scan_start_mjd[i]+3e-3,power_max+0.5,f\"{scan_id[i]:04}\",rotation=45,fontsize=fs)\n",
    "\n",
    "        axs1.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "        axs1.set_ylim(power_min,power_max)\n",
    "        axs1.fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "        axs1.set_xlabel('Time (UTC)',fontsize=fs)\n",
    "        axs1.set_xlim(t_plt[0],t_plt[-1])\n",
    "        axs1.tick_params(axis='both',labelsize=fs)\n",
    "        axs1.tick_params(axis='both',labelsize=fs)\n",
    "        \n",
    "        if freq_measured in RFI_mask_idx:\n",
    "            axs1.set_facecolor('lightgray')        \n",
    "        \n",
    "    axs1.set_ylabel('Power',fontsize=fs)\n",
    "    axs1.legend(fontsize=fs,markerscale=5)\n",
    "    axs1.grid()\n",
    "    \n",
    "    return\n",
    "\n",
    "    \n",
    "def DVA_Cross_Sections_singlescan(freq_chosen,scan_chosen,RFI_masked,freq_cross_section,time_idx):\n",
    "    \n",
    "    freq_measured = np.where(abs(freq-freq_chosen)<df)[0][0]\n",
    "    freq_below = np.where((freq>freq_chosen-5) & (freq<freq_chosen))[0]\n",
    "    freq_above = np.where((freq<freq_chosen+5) & (freq>freq_chosen))[0]\n",
    "          \n",
    "    scan_id_plot = scan_chosen\n",
    "    scan_idx = np.where(np.array(scan_id) == scan_id_plot)\n",
    "    \n",
    "    w = np.where((t_plt>=scan_start_mjd[scan_idx]) & (t_plt<=scan_stop_mjd[scan_idx]))[0]\n",
    "         \n",
    "    fig,axs1 = plt.subplots(1,1,figsize=(16,6)) \n",
    "            \n",
    "    if(freq_cross_section):\n",
    "            \n",
    "        data_plot_L = 10*np.log10(LL_set[w[time_idx],:])\n",
    "        data_plot_R = 10*np.log10(RR_set[w[time_idx],:])\n",
    "        \n",
    "        if RFI_masked:\n",
    "            data_plot_L[RFI_mask_idx] = np.nan\n",
    "            data_plot_R[RFI_mask_idx] = np.nan\n",
    "            \n",
    "        axs1.plot(freq,data_plot_L, label='LL',color='blue')\n",
    "        axs1.plot(freq,data_plot_R, label='RR',color='red')\n",
    "        axs1.vlines(freq_chosen, 0 , 100e9, color = 'purple')\n",
    "        axs1.set_ylim(power_min,power_max)\n",
    "        axs1.set_xlim(350,1030)\n",
    "        axs1.set_xlabel('Frequency',fontsize=fs)\n",
    "        axs1.tick_params(axis='both',labelsize=fs)\n",
    "    else:\n",
    "        axs1.plot(t_plt[w], 10*np.log10(LL_set[w,freq_measured]), label='LL',color='blue',zorder=1)\n",
    "        axs1.plot(t_plt[w], 10*np.log10(RR_set[w,freq_measured]), label='RR',color='red',zorder=1)\n",
    "        for ifreq in freq_above:\n",
    "            axs1.plot(t_plt[w], 10*np.log10(LL_set[w,ifreq]), alpha=0.5,color='C0',zorder=0)\n",
    "            axs1.plot(t_plt[w], 10*np.log10(RR_set[w,ifreq]), alpha=0.5,color='C1',zorder=0)\n",
    "        for ifreq in freq_below:\n",
    "            axs1.plot(t_plt[w], 10*np.log10(LL_set[w,ifreq]), alpha=0.5,color='C9',zorder=0)\n",
    "            axs1.plot(t_plt[w], 10*np.log10(RR_set[w,ifreq]), alpha=0.5,color='salmon',zorder=0)\n",
    "        axs1.vlines(t_plt[w][time_idx], 0 , 100e9, color = 'purple')\n",
    "        axaz = axs1.twiny()\n",
    "        axaz.set_xlim(az_set[w][0],az_set[w][-1])\n",
    "        axaz.tick_params(axis='x',labelsize=fs)\n",
    "\n",
    "        axs1.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "        axs1.set_ylim(power_min,power_max) #AO changed this to log scale limits\n",
    "        axs1.fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "        axs1.set_xlabel('Time (UTC)',fontsize=fs)\n",
    "        axs1.set_xlim(t_plt[w][0],t_plt[w][-1])\n",
    "        axs1.tick_params(axis='both',labelsize=fs)\n",
    "        axs1.set_title('Scan '+str(scan_id_plot)+' , el = '+str(scan_el[scan_idx[0][0]]),fontsize=fs)\n",
    "            \n",
    "        if freq_measured in RFI_mask_idx:\n",
    "             axs1.set_facecolor('lightgray')\n",
    "        \n",
    "    axs1.set_ylabel('Power',fontsize=fs)\n",
    "    axs1.legend(fontsize=fs,markerscale=5)\n",
    "    axs1.grid()\n",
    "    \n",
    "    return\n",
    "    \n",
    "def make_map(scan_pick,freq_pick=520.,dB1=74.,dB2=78.,*args,**kwargs):\n",
    "\n",
    "    fs=12\n",
    "    spts = 5\n",
    "    frq_idx = np.where(abs(freq-freq_pick)<df)[0][0]\n",
    "    \n",
    "    scan_idx = np.where(np.array(scan_id) == scan_pick)[0][0]\n",
    "    print('scan: ', scan_id[scan_idx])\n",
    "    print('frequency: ',freq_pick)\n",
    "    print('power range (dB): ',dB1,dB2)\n",
    "    \n",
    "    start_plt = Time(scan_start[scan_idx], format='isot',scale='utc').mjd\n",
    "    stop_plt = Time(scan_stop[scan_idx], format='isot',scale='utc').mjd        \n",
    "    t_idx = np.where((t_plt>=start_plt) & (t_plt<=stop_plt))[0]\n",
    "        \n",
    "    loc = EarthLocation(lat = 49.32*u.deg, lon = -119.62*u.deg, height = 546*u.m)\n",
    "    AltAzcoord = SkyCoord(alt=el_set[t_idx]*u.deg, az=az_set[t_idx]*u.deg, \n",
    "                      obstime=Time(t_set[t_idx],format='isot',scale='utc',location=loc),\n",
    "                      frame = 'altaz', location = loc)\n",
    "    l_scan = AltAzcoord.galactic.l.deg\n",
    "    b_scan = AltAzcoord.galactic.b.deg    \n",
    "    gal = SkyCoord(l_scan, b_scan, frame='galactic', unit=u.deg)\n",
    "    \n",
    "    wnoise = np.where(noise_set[t_idx] == 1)[0]\n",
    "    data_plot = 10*np.log10(LL_set[t_idx,frq_idx])\n",
    "    data_plot[wnoise] = np.nan\n",
    "    \n",
    "    fig = plt.figure(figsize=(14,7))\n",
    "    plt.subplot(111, projection='aitoff')\n",
    "    plt.scatter(-gal.l.wrap_at('180d').radian, gal.b.radian,s=spts,c=data_plot,\n",
    "                vmin=dB1,vmax=dB2,alpha=1)\n",
    "    positions = (-150*np.pi/180.,-120*np.pi/180.,-90*np.pi/180.,-60*np.pi/180.,-30*np.pi/180.,0,\n",
    "                    30*np.pi/180.,60*np.pi/180.,90*np.pi/180.,120*np.pi/180.,150*np.pi/180.)\n",
    "    labels = (\"150\",\"120\",\"90\",\"60\",\"30\",\"0\",\"330\",\"300\",\"270\",\"240\",\"210\")\n",
    "    plt.xticks(positions, labels,fontsize=fs)\n",
    "    plt.yticks(fontsize=fs)\n",
    "    plt.xlabel('Galactic Longitude (degrees)',fontsize=fs)\n",
    "    plt.ylabel('Galactic Latitude (degrees)',fontsize=fs)\n",
    "    plt.grid(True)\n",
    "\n",
    "    return \n",
    "    \n",
    "        \n",
    "def DVA_Visualization(waterfall_enabled,singlescan):\n",
    "\n",
    "    if waterfall_enabled:\n",
    "        if singlescan:\n",
    "            interact(DVA_Waterfall_View_singlescan,scan_chosen = scan_id, RFI_masked = False)\n",
    "        else:\n",
    "            interact(DVA_Waterfall_View, RFI_masked = False)\n",
    "    else:\n",
    "        if singlescan:\n",
    "            interact(DVA_Cross_Sections_singlescan, freq_chosen = (350, 1030, df), \n",
    "                     scan_chosen = scan_id, RFI_masked = False, freq_cross_section = False,\n",
    "                    time_idx = (0,1802))\n",
    "        else:\n",
    "            interact(DVA_Cross_Sections, freq_chosen = (350, 1030, df), \n",
    "                     freq_cross_section = False, time_idx = (0,len(t_plt)-1), RFI_masked = False)\n",
    "            \n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60a153a0",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "interact(DVA_Visualization,waterfall_enabled = True, singlescan = False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "676a6045",
   "metadata": {},
   "outputs": [],
   "source": [
    "make_map(2412,freq_pick=350,dB1=72,dB2=76)\n",
    "#make_map(2412)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4864d76",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6d39164",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5915290",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6b9c7b9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "72c266d8",
   "metadata": {},
   "source": [
    "## Leo's original code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a53806df",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = freq[1]-freq[0]\n",
    "\n",
    "def DVA_Waterfall_View():\n",
    "    #TODO: Add another DVA_Waterfall interactive_function such that I can change between LL_Set, RR_set, etc...\n",
    "    power_min = 66 #AO changed from 70\n",
    "    power_max = 78\n",
    "\n",
    "    fig,axs = plt.subplots(1,1,figsize=(15,10)) \n",
    "    fs = 16\n",
    "    \n",
    "    for i in range(0,len(scan_id)):\n",
    "        w = np.where((t_plt>=scan_start_mjd[i]) & (t_plt<=scan_stop_mjd[i]))[0]\n",
    "        extent = [scan_start_mjd[i],scan_stop_mjd[i],freq[0],freq[-1]]\n",
    "    \n",
    "        im = axs.imshow(10*np.log10(LL_set[w,:].T),aspect='auto',vmin=power_min,vmax=power_max,\n",
    "                        origin='lower',extent=extent,cmap='viridis')\n",
    "    \n",
    "    #im = axs.imshow(10.*np.log10(LL_set.T),aspect='auto',vmin=power_min,vmax=power_max,origin='lower',\n",
    "    #            extent=[t_plt[0],t_plt[-1],freq[0],freq[-1]])\n",
    "\n",
    "    divider = make_axes_locatable(axs)\n",
    "    cax = divider.append_axes('right', size='2%', pad=0.05)\n",
    "    cbar = fig.colorbar(im, cax=cax, orientation='vertical')\n",
    "    cbar.ax.tick_params(labelsize=fs) \n",
    "    cbar.set_label('Power (dB)', fontsize=fs)\n",
    "\n",
    "    axs.set_xlim(t_plt[0],t_plt[-1])\n",
    "    axs.set_ylim(freq[0],freq[-1])\n",
    "    axs.tick_params(axis='both', labelsize=fs)\n",
    "    axs.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "    axs.fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "    axs.set_xlabel('Time (UTC)',fontsize=fs)        \n",
    "    axs.set_ylabel('Frequency (MHz)',fontsize=fs)\n",
    "\n",
    "\n",
    "def DVA_Cross_Sections(freq_chosen, time_idx, freq_cross_section):\n",
    "    freq_measured = np.where(abs(freq-freq_chosen)<df)[0][0]\n",
    "\n",
    "    fs = 12    \n",
    "    fig,axs1 = plt.subplots(1,1,figsize=(16,6))  \n",
    "    # TODO: use the log of LL_set_clean\n",
    "    if(freq_cross_section):\n",
    "        power_min = 66 #AO added limits\n",
    "        power_max = 78\n",
    "        # AO changed to log scale and added RR:\n",
    "        axs1.plot(freq,10*np.log10(LL_set[time_idx,:]), label='LL')\n",
    "        axs1.plot(freq,10*np.log10(RR_set[time_idx,:]), label='RR')\n",
    "        axs1.vlines(freq_chosen, 0 , 100e9, color = 'red')\n",
    "        axs1.set_ylim(power_min,power_max) #AO changed this to log scale limits\n",
    "        axs1.set_xlim(350,1050)\n",
    "        #axs1.set_ylim(np.min(LL_set[:,freq_measured]), np.max(LL_set[time_idx,:]))\n",
    "        axs1.set_xlabel('Frequency',fontsize=fs)\n",
    "    else:\n",
    "        power_min = 66 #AO added limits\n",
    "        power_max = 78\n",
    "        # AO changed to log scale and added RR:\n",
    "        axs1.scatter(t_plt, 10*np.log10(LL_set[:,freq_measured]), label='LL',s=0.5)\n",
    "        axs1.scatter(t_plt, 10*np.log10(RR_set[:,freq_measured]), label='RR',s=0.5)\n",
    "        axs1.vlines(t_plt[time_idx], 0 , 100e9, color = 'red')\n",
    "\n",
    "        axs1.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "        axs1.set_ylim(power_min,power_max) #AO changed this to log scale limits\n",
    "        axs1.fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "        axs1.set_xlabel('Time (UTC)',fontsize=fs)\n",
    "        axs1.set_xlim(t_plt[0],t_plt[-1])\n",
    "    axs1.set_ylabel('Power',fontsize=fs)\n",
    "    axs1.legend()\n",
    "    axs1.grid() # AO added grid\n",
    "\n",
    "def DVA_Visualization(waterfall_enabled):\n",
    "    if waterfall_enabled:\n",
    "        interact(DVA_Waterfall_View)\n",
    "    else:\n",
    "        interact(DVA_Cross_Sections, freq_chosen = (350, 1000, df), freq_cross_section = False, time_idx = (0,len(t_plt)-1))\n",
    "\n",
    "\n",
    "interact(DVA_Visualization, waterfall_enabled = True)\n",
    "\n",
    "#TODO: plot power in DB"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "ad2bdc8ecc057115af97d19610ffacc2b4e99fae6737bb82f5d7fb13d2f2c186"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}