dva_baseline_visualization.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 20,
   "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",
    "import matplotlib.patches as patches\n",
    "\n",
    "\n",
    "#### Change the directory to where the files are located\" ####\n",
    "day ='25'\n",
    "directory = '../DVA/Data_Files/DVA_Day_Surveys/'\n",
    "\n",
    "### Change the scan details to your current scan\n",
    "# df = 0.5        \n",
    "# dt = 6.9444431574083865e-06\n",
    "\n",
    "#TODO: actually do the scan properties definition in the beginning\n",
    "##############################################################\n",
    "def find_nearest_idx(array, value):\n",
    "    array = np.asarray(array)\n",
    "    idx = (np.abs(array - value)).argmin()\n",
    "    return idx\n",
    "\n",
    "def Apply_Mask(mask, input_array): #TODO: This needs to be updated to implement all frequencies\n",
    "    \n",
    "    masked_copy = input_array.copy()\n",
    "    mask_idx = np.where(mask == 1)\n",
    "    masked_copy[mask_idx, :] = np.nan\n",
    "    return(masked_copy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in the file listing azimuth scan start and stop times:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1052 2022-07-04T03:57:53.110 2022-07-04T04:15:55.184\n",
      "2565 2022-07-04T04:20:37.406 2022-07-04T04:38:39.310\n",
      "2566 2022-07-04T04:40:22.266 2022-07-04T04:58:24.003\n",
      "1271 2022-07-04T05:03:42.716 2022-07-04T05:21:44.818\n",
      "1272 2022-07-04T05:23:27.291 2022-07-04T05:41:29.292\n",
      "2857 2022-07-04T05:48:35.269 2022-07-04T06:06:37.317\n",
      "0124 2022-07-04T06:11:19.190 2022-07-04T06:29:21.101\n",
      "1637 2022-07-04T06:34:03.920 2022-07-04T06:52:06.005\n",
      "0342 2022-07-04T06:57:23.897 2022-07-04T07:15:25.930\n",
      "1927 2022-07-04T07:21:55.804 2022-07-04T07:39:58.034\n",
      "1928 2022-07-04T07:42:16.530 2022-07-04T08:00:18.319\n",
      "2001 2022-07-04T08:03:48.814 2022-07-04T08:21:51.062\n",
      "2002 2022-07-04T08:24:09.373 2022-07-04T08:42:11.545\n",
      "2075 2022-07-04T08:45:41.821 2022-07-04T09:03:44.188\n",
      "2148 2022-07-04T09:07:50.625 2022-07-04T09:25:52.832\n",
      "2149 2022-07-04T09:27:35.229 2022-07-04T09:45:37.044\n",
      "2150 2022-07-04T09:47:56.026 2022-07-04T10:05:57.828\n",
      "2151 2022-07-04T10:07:40.752 2022-07-04T10:25:42.662\n",
      "2152 2022-07-04T10:27:25.410 2022-07-04T10:45:27.414\n",
      "2153 2022-07-04T10:47:46.110 2022-07-04T11:05:48.028\n",
      "0858 2022-07-04T11:10:30.298 2022-07-04T11:28:32.281\n",
      "0859 2022-07-04T11:30:50.747 2022-07-04T11:48:52.934\n",
      "2372 2022-07-04T11:53:35.348 2022-07-04T12:11:37.677\n"
     ]
    }
   ],
   "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",
    "\n",
    "# Read in the data and store it in arrays:\n",
    "with open(directory+'DVAsurvey_phase1_day0'+day+'.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",
    "        \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])\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",
   "metadata": {},
   "source": [
    "## Read in scan files and stich them together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1052\n",
      "<HDF5 file \"dva_survey_raw_scan_1052.h5\" (mode r)>\n",
      "2565\n",
      "<HDF5 file \"dva_survey_raw_scan_2565.h5\" (mode r)>\n",
      "2566\n",
      "<HDF5 file \"dva_survey_raw_scan_2566.h5\" (mode r)>\n",
      "1271\n",
      "<HDF5 file \"dva_survey_raw_scan_1271.h5\" (mode r)>\n",
      "1272\n",
      "<HDF5 file \"dva_survey_raw_scan_1272.h5\" (mode r)>\n"
     ]
    }
   ],
   "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",
    "# Use one of the scans to get the list of frequencies:\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",
    "file = h5py.File(directory+'dva_survey_raw_scan_'+scan0+'.h5','r')\n",
    "freq = file['data']['beam_0']['band_SB0']['frequency'][::freq_channel_increment]/1e6\n",
    "df = freq[1] - freq[0]\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[0:5]:\n",
    "    print(i)\n",
    "    # select the file:\n",
    "    file = h5py.File(directory+'dva_survey_raw_scan_'+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\n",
    "\n",
    "#NOTE: Local change because Leo only has 5 scans on his computer\n",
    "scan_id = [1052, 2565, 2566, 1271, 1272]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in persistent RFI mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 0\n",
    "RFI_mask_idx = []\n",
    "with open('../DVA/Data_Files/RFIpersist_mask.txt') as fp:\n",
    "    for line in fp:\n",
    "        if i>0: \n",
    "            RFI_mask_idx.append(int(line.split()[0]))\n",
    "        i=i+1\n",
    "RFI_mask_idx = np.array(RFI_mask_idx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Polarized Intensity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "polarized = []\n",
    "for i,j in zip(reRL_set,imRL_set):\n",
    "    PI = np.sqrt((i**2)+(j**2))\n",
    "    polarized.append(PI)\n",
    "polarized_set = np.array(polarized)\n",
    "\n",
    "noise_idx = np.array(np.where(noise_set == 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# df = freq[1]-freq[0]\n",
    "# power_min = 65\n",
    "# power_max = 79\n",
    "# fs = 14\n",
    "\n",
    "# freq_chosen = 844\n",
    "# freq_idx = find_nearest_idx(freq, freq_chosen)\n",
    "# freq_measured = freq[freq_idx]\n",
    "# time_idx = (0,1802)\n",
    "# abs_threshold = 5e5\n",
    "\n",
    "def DVA_Singlescan_Time_Cross_Section(scan_chosen):\n",
    "\n",
    "    scan_id_plot = scan_chosen\n",
    "    scan_idx = np.where(np.array(scan_id) == scan_id_plot)\n",
    "    scan_duration_idx = np.where((t_plt>=scan_start_mjd[scan_idx]) & (t_plt<=scan_stop_mjd[scan_idx]))[0]\n",
    "\n",
    "    scan_baseline = np.nanmedian(polarized_set[scan_duration_idx,freq_idx])\n",
    "\n",
    "    possible_RFI_idxes = np.where(polarized_set[:, freq_idx] >= scan_baseline*3)\n",
    "    polarized_set_masked = polarized_set.copy()\n",
    "    polarized_set_masked[possible_RFI_idxes, :] = np.nan\n",
    "\n",
    "\n",
    "    fig,axs1 = plt.subplots(1,1,figsize=(16,6))  \n",
    "    axs1.plot(t_plt[scan_duration_idx], polarized_set[scan_duration_idx,freq_idx],color='red', linewidth = 0.5, label = \"Possible RFI\")\n",
    "    axs1.plot(t_plt[scan_duration_idx], polarized_set_masked[scan_duration_idx,freq_idx],color='black', linewidth = 0.65, label = \"Masked Polarized Set\")\n",
    "\n",
    "\n",
    "    axs1.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "    axs1.fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "    axs1.set_xlabel('Time (UTC)',fontsize=fs)\n",
    "    axs1.set_ylabel(\"Power\")\n",
    "    axs1.set_xlim(t_plt[scan_duration_idx][0],t_plt[scan_duration_idx][-1])\n",
    "    axs1.set_title(\"Polarized set of scan number ({:.0f}) with baseline {:.0f}\".format(scan_chosen, scan_baseline))\n",
    "\n",
    "    # axs1.axhline(abs_threshold,xmin=0,xmax=3,c=\"green\",linewidth=1,zorder=0, label = \"Absolute Threshold\")\n",
    "    axs1.axhline(scan_baseline*3,xmin=0,xmax=3,c=\"purple\",linewidth=1,zorder=0, label = \"Relative Threshold\")\n",
    "    axs1.axhline(scan_baseline,xmin=0,xmax=3,c=\"red\",linewidth=1,zorder=0, label = \"Scan Baseline\")\n",
    "    axs1.set_ylim(ymax=1e6, ymin=0)\n",
    "    axs1.tick_params(axis='both',labelsize=fs)   \n",
    "    axs1.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Plot_Baseline_Bandwitdh(scan_chosen):\n",
    "    baseline_set = []\n",
    "    scan_id_plot = scan_chosen\n",
    "    scan_idx = np.where(np.array(scan_id) == scan_id_plot)\n",
    "    scan_duration_idx = np.where((t_plt>=scan_start_mjd[scan_idx]) & (t_plt<=scan_stop_mjd[scan_idx]))[0]\n",
    "    for freq_idx in range(0, len(freq)):\n",
    "        freq_baseline = np.nanmedian(polarized_set[scan_duration_idx,freq_idx])\n",
    "        baseline_set.append(freq_baseline)\n",
    "\n",
    "    fig,axs1 = plt.subplots(1,1,figsize=(16,6))  \n",
    "\n",
    "    axs1.plot(freq,baseline_set,color='blue') \n",
    "    axs1.set_title(\"Baseline Distribution of polarized scan {:.0f}\".format(scan_chosen))\n",
    "    axs1.set_ylim(ymax=5e6, ymin=0)\n",
    "    axs1.set_xlim(freq[0],freq[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "polarized_set[noise_idx, :] = np.nan\n",
    "polarized_set[:,RFI_mask_idx] = np.nan\n",
    "\n",
    "LL_set[noise_idx, :] = np.nan\n",
    "LL_set[:,RFI_mask_idx] = np.nan\n",
    "\n",
    "RR_set[noise_idx, :] = np.nan\n",
    "RR_set[:,RFI_mask_idx] = np.nan\n",
    "\n",
    "# interact(DVA_Singlescan_Time_Cross_Section,  scan_chosen = scan_id)\n",
    "# interact(Plot_Baseline_Bandwitdh,  scan_chosen = scan_id)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Combine_Into_RFI_Events(possible_RFI):\n",
    "    Combined_RFI_events = [] \n",
    "    remaining_idxes = len(possible_RFI)\n",
    "\n",
    "    current_idx = 1\n",
    "    rfi_start = possible_RFI[current_idx]\n",
    "    rfi_duration = 1\n",
    "\n",
    "    while remaining_idxes >= 2:\n",
    "        remaining_idxes = (len(possible_RFI) - current_idx)\n",
    "        previous_idx = current_idx - 1\n",
    "        idx_gap = (possible_RFI[current_idx] - possible_RFI[previous_idx])\n",
    "        if(idx_gap <= 2):\n",
    "            rfi_duration += 1\n",
    "            current_idx += 1\n",
    "        else:\n",
    "            #Close current event\n",
    "            Combined_RFI_events.append([rfi_start, rfi_duration])\n",
    "            #Initiate next event starting on the current_idx\n",
    "            rfi_duration = 1\n",
    "            current_idx += 1\n",
    "            rfi_start = possible_RFI[current_idx]\n",
    "\n",
    "    return Combined_RFI_events #[idx, idx_duration]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spectrum Scan Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RFI_Start_Found(time_array, spectrum_power_array, idx, slope_threshold):\n",
    "    RFI_Start_Found = False\n",
    "    spectrum_gradient = np.gradient(spectrum_power_array, 3)\n",
    "    current_slope = spectrum_gradient[idx-1]\n",
    "    next_slope = spectrum_gradient[idx]\n",
    "    if((current_slope <= slope_threshold) and (next_slope >= slope_threshold)):\n",
    "        RFI_Start_Found = True\n",
    "        # print(\"RFI start information:\\n - current slope:\", current_slope,\"next slope:\", next_slope, \"Slope threshold:\", slope_threshold)\n",
    "    return RFI_Start_Found\n",
    "\n",
    "\n",
    "def RFI_End_Found(spectrum_power_array, rfi_end_idx, RFI_start_value):\n",
    "    retVal = False\n",
    "    if spectrum_power_array[rfi_end_idx] <= RFI_start_value:\n",
    "        retVal = True\n",
    "    return retVal\n",
    "\n",
    "def RFI_Spectrum_Scan(time_array, spectrum_power_array, freq_idx, slope_threshold): #NOTE: This might be an issue if I have a nan value on my array\n",
    "    scan_bandwidth = 10 #[MHz]\n",
    "    start_idx = int(freq_idx - int((scan_bandwidth/df)/2))  #This changes the bandwidth from MHz to idxes\n",
    "    end_idx = int(freq_idx + int((scan_bandwidth/df)/2))\n",
    "    RFI_confirmed = False\n",
    "    RFI_spectral_thickness = 0\n",
    "    for rfi_start_idx in range(start_idx, end_idx):                                             #Scan over the entire bandwith interval\n",
    "        if(RFI_Start_Found(time_array, spectrum_power_array, rfi_start_idx, slope_threshold)):     #If I found the starting pattern\n",
    "            RFI_start_value = spectrum_power_array[rfi_start_idx]\n",
    "            rfi_minimum_length = 10     #This exists to avoid the algorithm from fake crossing near the start due to fluctuations.\n",
    "            for rfi_end_idx in range(rfi_start_idx + rfi_minimum_length, end_idx):                                   #Finish looking at the bandwidth interval looking for the end\n",
    "                # if(Spectrum_RFI_End_Found(time_array, spectrum_power_array, rfi_end_idx, slope_threshold)):     #If I found the ending pattern          TODO: I have a problem in the end found algorithm\n",
    "                if(RFI_End_Found(spectrum_power_array, rfi_end_idx, RFI_start_value)):     #If the RFI crosses below the starting value         TODO: I have a problem in the end found algorithm\n",
    "                    RFI_confirmed = True      \n",
    "                    break\n",
    "                    # RFI_spectral_thickness = (rfi_end_idx - rfi_start_idx)\n",
    "                else:\n",
    "                    continue\n",
    "            if(RFI_confirmed):                                                                          #and stop looping throught the interval for efficiency\n",
    "                break\n",
    "        else:\n",
    "            rfi_end_idx = end_idx\n",
    "    return RFI_confirmed, rfi_start_idx, (rfi_end_idx)#RFI_spectral_thickness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RFI_Excision stuff below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "def DVA_Find_Possible_RFI_Events(scan_chosen, freq_idx, baseline_multiplier):\n",
    "    scan_id_plot = scan_chosen\n",
    "    scan_idx = np.where(np.array(scan_id) == scan_id_plot)\n",
    "    scan_duration_idx = np.where((t_plt>=scan_start_mjd[scan_idx]) & (t_plt<=scan_stop_mjd[scan_idx]))[0]\n",
    "\n",
    "    scan_baseline = np.nanmedian(polarized_set[scan_duration_idx,freq_idx])\n",
    "    scan_threshold = scan_baseline*baseline_multiplier\n",
    "\n",
    "    possible_RFI_idxes = np.where(polarized_set[:, freq_idx] >= scan_threshold)\n",
    "    possible_RFI_events = Combine_Into_RFI_Events(possible_RFI_idxes[0])\n",
    "    return possible_RFI_events  #Returns [time_idx, idx_duration]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TODO:\n",
    "- Look at other frequencies\n",
    "- Loop through all the frequencies of a single scan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of possible RFI event: 49\n",
      "Number of confirmed RFI regions: 37\n"
     ]
    }
   ],
   "source": [
    "def RFI_Verification(possible_RFI_events, freq_slope_threshold, event, freq_idx):\n",
    "    rfi_confirmed = False\n",
    "    for time_idx in range(possible_RFI_events[event][0], possible_RFI_events[event][0]+ possible_RFI_events[event][1]):\n",
    "        event_verification_result = RFI_Spectrum_Scan(t_plt, polarized_set[time_idx, :], freq_idx, freq_slope_threshold)\n",
    "        if event_verification_result[0]:\n",
    "            rfi_confirmed = True\n",
    "            break\n",
    "    return rfi_confirmed, event_verification_result\n",
    "\n",
    "\n",
    "def RFI_Detection(scan, freq_slope_threshold, freq_chosen, baseline_multiplier):\n",
    "    freq_idx = find_nearest_idx(freq, freq_chosen)\n",
    "    freq_measured = freq[freq_idx]\n",
    "\n",
    "\n",
    "    confirmed_RFI_results = []\n",
    "    possible_RFI_events = DVA_Find_Possible_RFI_Events(scan, freq_idx, baseline_multiplier)\n",
    "    print(\"Number of possible RFI event:\", len(possible_RFI_events))\n",
    "    for event in range(0, len(possible_RFI_events)-1):\n",
    "        rfi_confirmed, event_verification_result = RFI_Verification(possible_RFI_events, freq_slope_threshold, event, freq_idx)      \n",
    "        if rfi_confirmed:  \n",
    "            # DETERMINE RFI REGION --------------------------------------------------------------------------------------------------------\n",
    "            t1_plt = possible_RFI_events[event][0]                                      #Start time     [idx]\n",
    "            t2_plt = possible_RFI_events[event][0] + possible_RFI_events[event][1]      #End time       [idx]\n",
    "            freq1 = event_verification_result[1]                                        #Start freq     [idx]\n",
    "            freq2 = event_verification_result[2]                                        #End freq       [idx]\n",
    "            \n",
    "            confirmed_RFI_results.append([t1_plt, t2_plt, freq1, freq2])\n",
    "    print(\"Number of confirmed RFI regions:\", len(confirmed_RFI_results))\n",
    "    return confirmed_RFI_results\n",
    "\n",
    "confirmed_RFI_results = RFI_Detection(scan = 1052, freq_slope_threshold = 1e5, freq_chosen = 844, baseline_multiplier = 3)\n",
    "# interact(RFI_Detection(scan = 1045))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d32dbdc9c03f47d5800a0cd341168498",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=17, description='rfi_number', max=34), Output()), _dom_classes=('widget-…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.DVA_Plot_RFI(rfi_number)>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def DVA_Plot_RFI(rfi_number):\n",
    "    #Arbitraty constants for visualization   \n",
    "    time_framing = 50\n",
    "    freq_framing = 50\n",
    "    time_buffer = 10\n",
    "    freq_buffer = 10\n",
    "\n",
    "    # DETERMINE RFI REGION --------------------------------------------------------------------------------------------------------\n",
    "    t1_plt = confirmed_RFI_results[rfi_number][0]\n",
    "    t2_plt = confirmed_RFI_results[rfi_number][1]\n",
    "    freq1 = confirmed_RFI_results[rfi_number][2]\n",
    "    freq2 = confirmed_RFI_results[rfi_number][3]\n",
    "    print(\"RFI starting frequency: {:.2f}.\\nRFI ending frequency: {:.2f}\".format(freq[freq1], freq[freq2]))\n",
    "\n",
    "    # PLOTTING WATERFALL RFI ------------------------------------------------------------------------------------------------------\n",
    "\n",
    "    RFI_duration_time = t_plt[t1_plt-time_framing: t2_plt+time_framing]\n",
    "    RFI_freq_bandwith = freq[freq1 - freq_framing: freq2 + freq_framing]\n",
    "\n",
    "    fig,axs = plt.subplots(1,2,figsize=(25,10)) \n",
    "    fs = 16\n",
    "\n",
    "    # Create a Rectangle patch ---------------------------------------------------------------------------------------------------\n",
    "    RFI_detected = patches.Rectangle((t_plt[t1_plt], freq[freq1]), np.abs(t_plt[t2_plt] - t_plt[t1_plt]), np.abs(freq[freq2] - freq[freq1]), linewidth=1.5, edgecolor='r', facecolor='none')\n",
    "    RFI_detected_LL = patches.Rectangle((t_plt[t1_plt], freq[freq1]), np.abs(t_plt[t2_plt] - t_plt[t1_plt]), np.abs(freq[freq2] - freq[freq1]), linewidth=1.5, edgecolor='r', facecolor='none')\n",
    "    # Add the patch to the Axes0\n",
    "    axs[0].add_patch(RFI_detected)\n",
    "    axs[1].add_patch(RFI_detected_LL)\n",
    "\n",
    "    #Polarized Waterfall ---------------------------------------------------------------------------------------------------------------------------\n",
    "    power_min = np.nanmedian((10.*np.log10(polarized_set[t1_plt-time_buffer: t2_plt+time_buffer,freq1-freq_buffer: freq2+freq_buffer].T)))\n",
    "    power_max = np.nanmax((10.*np.log10(polarized_set[t1_plt-time_buffer: t2_plt+time_buffer,freq1-freq_buffer: freq2+freq_buffer].T)))\n",
    "\n",
    "    im_pol = axs[0].imshow(10.*np.log10(polarized_set[t1_plt-time_framing: t2_plt+time_framing,freq1-freq_framing: freq2+freq_framing].T),aspect='auto',vmin=power_min,vmax=power_max,origin='lower',\n",
    "            extent=[RFI_duration_time[0],RFI_duration_time[-1],RFI_freq_bandwith[0],RFI_freq_bandwith[-1]])\n",
    "\n",
    "    divider = make_axes_locatable(axs[0])\n",
    "    cax = divider.append_axes('right', size='2%', pad=0.05)\n",
    "    cbar = fig.colorbar(im_pol, cax=cax, orientation='vertical')\n",
    "    cbar.ax.tick_params(labelsize=fs) \n",
    "    cbar.set_label('Power (dB)', fontsize=fs)\n",
    "\n",
    "    axs[0].set_xlim(RFI_duration_time[0],RFI_duration_time[-1])\n",
    "    axs[0].tick_params(axis='both', labelsize=fs)\n",
    "    axs[0].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "    axs[0].fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "    axs[0].set_xlabel('Time (UTC)',fontsize=fs)        \n",
    "    axs[0].set_ylabel('Frequency (MHz)',fontsize=fs)\n",
    "    axs[0].set_title('Polarized Set at Frequency {:.2f}'.format(freq_measured), size = fs)\n",
    "\n",
    "    #LL Waterfall ---------------------------------------------------------------------------------------------------------------------------\n",
    "    power_min = np.nanmedian((10.*np.log10(LL_set[t1_plt-time_buffer: t2_plt+time_buffer,freq1-freq_buffer: freq2+freq_buffer].T)))\n",
    "    power_max = np.nanmax((10.*np.log10(LL_set[t1_plt-time_buffer: t2_plt+time_buffer,freq1-freq_buffer: freq2+freq_buffer].T)))\n",
    "\n",
    "    im_LL = axs[1].imshow(10.*np.log10(LL_set[t1_plt-time_framing: t2_plt+time_framing,freq1-freq_framing: freq2+freq_framing].T),aspect='auto',vmin=power_min,vmax=power_max,origin='lower',\n",
    "            extent=[RFI_duration_time[0],RFI_duration_time[-1],RFI_freq_bandwith[0],RFI_freq_bandwith[-1]])\n",
    "\n",
    "    divider2 = make_axes_locatable(axs[1])\n",
    "    cax2 = divider2.append_axes('right', size='2%', pad=0.05)\n",
    "    cbar2 = fig.colorbar(im_LL, cax=cax2, orientation='vertical')\n",
    "    cbar2.ax.tick_params(labelsize=fs) \n",
    "    cbar2.set_label('Power (dB)', fontsize=fs)\n",
    "\n",
    "    axs[1].set_xlim(RFI_duration_time[0],RFI_duration_time[-1])\n",
    "    axs[1].tick_params(axis='both', labelsize=fs)\n",
    "    axs[1].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "    axs[1].fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "    axs[1].set_xlabel('Time (UTC)',fontsize=fs)        \n",
    "    axs[1].set_ylabel('Frequency (MHz)',fontsize=fs)\n",
    "    axs[1].set_title('LL Set at Frequency {:.2f}'.format(freq_measured), size = fs)\n",
    "\n",
    "interact(DVA_Plot_RFI, rfi_number = (0,34))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.8 ('base')",
   "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.8.8"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "ad2bdc8ecc057115af97d19610ffacc2b4e99fae6737bb82f5d7fb13d2f2c186"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}