dva_survey_daily_files_rasterscans.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "6cc67b8d",
   "metadata": {},
   "source": [
    "## Read in separate raw data files and create files of individual raster scans\n",
    "### A. Ordog, June 2022\n",
    "### NOTE: please do not modify\n",
    "#### July 07, 2022: added in/out directories\n",
    "#### July 12, 2022: added option of more than one raster scan for evening and morning\n",
    "#### Aug 15, 2022: removed old code no longer needed, included check of elevation range during scans"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "39d1f524",
   "metadata": {},
   "source": [
    "## Import packages and define survey day and file directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1aa737b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dva_sdhdf_combine_v3\n",
    "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 astropy.modeling import models, fitting\n",
    "\n",
    "day ='38'\n",
    "Nm = 0 # number of morning raster scans\n",
    "Ne = 1 # number of evening raster scans\n",
    "\n",
    "#### Change the directory to where the files are located\" ####\n",
    "#dir_in = '/home2/DATA_AO/DVA_DATA/survey_phase3_day'+day+'/'\n",
    "#dir_out = '/media/ordoga/DVA_data/survey_phase3_day'+day+'/'\n",
    "#dir_out = '/media/ordoga/15m_band1_survey/dva_phase3/survey_phase3_day'+day+'/'\n",
    "##############################################################\n",
    "\n",
    "dir_in = '/srv/data/dva/survey_raster_noise_corr/'\n",
    "#dir_out = '/media/ordoga/DVA_data/survey_phase3_day'+day+'/'\n",
    "dir_out = '/srv/data/dva/survey_raster_noise_corr/'\n",
    "##############################################################"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "56d940fb",
   "metadata": {},
   "source": [
    "## Read in the file listing raster scan start and stop times:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f39809e",
   "metadata": {},
   "outputs": [],
   "source": [
    "raster1_start = []\n",
    "raster2_start = []\n",
    "raster1_stop = []\n",
    "raster2_stop = []\n",
    "RA1 = []\n",
    "RA2 = []\n",
    "DEC1 = []\n",
    "DEC2 = []\n",
    "\n",
    "rast_list = ['a','b','c']\n",
    "\n",
    "for i in range(0,Ne):\n",
    "    \n",
    "    if Ne == 1:\n",
    "        rast1name = 'dva_survey_phase3_day0'+day+'_raster1.txt'\n",
    "    else:\n",
    "        rast1name = 'dva_survey_phase3_day0'+day+'_raster1'+rast_list[i]+'.txt'\n",
    "        \n",
    "    with open(dir_out+rast1name) as fp:\n",
    "        for line in fp:  \n",
    "            raster1_start.append(line.split()[3])\n",
    "            raster1_stop.append(line.split()[4])\n",
    "            RA1.append(float(line.split()[1])*360/24)\n",
    "            DEC1.append(float(line.split()[2]))\n",
    "            \n",
    "for i in range(0,Nm):\n",
    "    \n",
    "    if Nm == 1:\n",
    "        rast2name = 'dva_survey_phase1_day0'+day+'_raster2.txt'\n",
    "    else:\n",
    "        rast2name = 'dva_survey_phase1_day0'+day+'_raster2'+rast_list[i]+'.txt'\n",
    "        \n",
    "    with open(dir_out+rast2name) as fp:\n",
    "        for line in fp:  \n",
    "            raster2_start.append(line.split()[3])\n",
    "            raster2_stop.append(line.split()[4])\n",
    "            RA2.append(float(line.split()[1])*360/24)\n",
    "            DEC2.append(float(line.split()[2]))\n",
    "\n",
    "print('Evening raster starts:')\n",
    "print(raster1_start)\n",
    "print('Evening raster stops:')\n",
    "print(raster1_stop)\n",
    "print('')\n",
    "print('Morning raster starts:')\n",
    "print(raster2_start)\n",
    "print('Morning raster stops:')\n",
    "print(raster2_stop)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "62b2e2cc",
   "metadata": {},
   "source": [
    "## Make files for each raster scan\n",
    "\n",
    "### ***Note: skip this if files already made"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5e30df0",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "imp.reload(dva_sdhdf_combine_v3)\n",
    "fsamp = 2\n",
    "rast_list = ['a','b','c']\n",
    "\n",
    "for i in range(0,Ne):\n",
    "#for i in range(0,1):\n",
    "    \n",
    "    if Ne == 1:\n",
    "        rast1name = 'dva_survey_phase1_day0'+day+'_raster1_median'\n",
    "    else:\n",
    "        rast1name = 'dva_survey_phase1_day0'+day+'_raster1'+rast_list[i]\n",
    "        \n",
    "    dva_sdhdf_combine_v3.combine(dir_in,dir_out,raster1_start[i],raster1_stop[i],\n",
    "                                 rast1name,freq_s=fsamp,freq_avg=True)\n",
    "\n",
    "for i in range(0,Nm):\n",
    "#for i in range(1,Nm):   \n",
    "    if Nm == 1:\n",
    "        rast2name = 'dva_survey_phase1_day0'+day+'_raster2'\n",
    "    else:\n",
    "        rast2name = 'dva_survey_phase1_day0'+day+'_raster2'+rast_list[i]\n",
    "        \n",
    "    dva_sdhdf_combine_v3.combine(dir_in,dir_out,raster2_start[i],raster2_stop[i],\n",
    "                                 rast2name,freq_s=fsamp,freq_avg=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "5c1ab615",
   "metadata": {},
   "source": [
    "## Read in raster scan files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f2d22c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "rast_list = ['a','b','c']\n",
    "\n",
    "dec_rast1   = []; dec_rast2   = []\n",
    "ra_rast1    = []; ra_rast2    = []\n",
    "el_rast1    = []; el_rast2    = []\n",
    "az_rast1    = []; az_rast2    = []\n",
    "t_rast1     = []; t_rast2     = []\n",
    "noise_rast1 = []; noise_rast2 = []\n",
    "RR_rast1    = []; RR_rast2    = []\n",
    "LL_rast1    = []; LL_rast2    = []\n",
    "reRL_rast1  = []; reRL_rast2  = []\n",
    "imRL_rast1  = []; imRL_rast2  = []\n",
    "\n",
    "for i in range(0,Ne):\n",
    "    \n",
    "    print(i)\n",
    "    \n",
    "    if Ne == 1:\n",
    "        rast1name = 'dva_survey_phase1_day0'+day+'_raster1.h5'\n",
    "    else:\n",
    "        rast1name = 'dva_survey_phase1_day0'+day+'_raster1'+rast_list[i]+'.h5'\n",
    "        \n",
    "    file = h5py.File(dir_out+rast1name,'r')\n",
    "    dataset1 = file['data']['beam_0']['band_SB0']['scan_0'] \n",
    "    dec_rast1.append(dataset1['metadata']['declination'])\n",
    "    ra_rast1.append(dataset1['metadata']['right_ascension'])\n",
    "    el_rast1.append(dataset1['metadata']['elevation'])\n",
    "    az_rast1.append(dataset1['metadata']['azimuth'])\n",
    "    t_rast1.append(dataset1['metadata']['utc'])\n",
    "    noise_rast1.append(dataset1['metadata']['noise_state'])  \n",
    "    RR_rast1.append(dataset1['data'][:,0,:])\n",
    "    LL_rast1.append(dataset1['data'][:,1,:])\n",
    "    reRL_rast1.append(dataset1['data'][:,2,:])\n",
    "    imRL_rast1.append(dataset1['data'][:,3,:])\n",
    "        \n",
    "\n",
    "for i in range(0,Nm):\n",
    "    \n",
    "    print(i)\n",
    "    \n",
    "    if Nm == 1:\n",
    "        rast2name = 'dva_survey_phase1_day0'+day+'_raster2.h5'\n",
    "    else:\n",
    "        rast2name = 'dva_survey_phase1_day0'+day+'_raster2'+rast_list[i]+'.h5'\n",
    "        \n",
    "    file = h5py.File(dir_out+rast2name,'r')\n",
    "    dataset2 = file['data']['beam_0']['band_SB0']['scan_0']\n",
    "    dec_rast2.append(dataset2['metadata']['declination'])\n",
    "    ra_rast2.append(dataset2['metadata']['right_ascension'])\n",
    "    el_rast2.append(dataset2['metadata']['elevation'])\n",
    "    az_rast2.append(dataset2['metadata']['azimuth'])\n",
    "    t_rast2.append(dataset2['metadata']['utc'])\n",
    "    noise_rast2.append(dataset2['metadata']['noise_state'])  \n",
    "    RR_rast2.append(dataset2['data'][:,0,:])\n",
    "    LL_rast2.append(dataset2['data'][:,1,:])\n",
    "    reRL_rast2.append(dataset2['data'][:,2,:])\n",
    "    imRL_rast2.append(dataset2['data'][:,3,:])\n",
    "    \n",
    "freq_rast = file['data']['beam_0']['band_SB0']['frequency'][:]/1e6\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97c1fca2",
   "metadata": {},
   "outputs": [],
   "source": [
    "ra_rast1_deg = []\n",
    "t_plt_rast1 = []\n",
    "\n",
    "ra_rast2_deg = []\n",
    "t_plt_rast2 = []\n",
    "\n",
    "for i in range(0,Ne):\n",
    "    \n",
    "    ra_rast1_deg.append(ra_rast1[i]*360./24.)\n",
    "    t_plt_rast1.append(Time(t_rast1[i], format='isot',scale='utc').mjd)\n",
    "    \n",
    "for i in range(0,Nm):\n",
    "    \n",
    "    ra_rast2_deg.append(ra_rast2[i]*360./24.)\n",
    "    t_plt_rast2.append(Time(t_rast2[i], format='isot',scale='utc').mjd)\n",
    "\n",
    "print(t_plt_rast1)\n",
    "print(t_plt_rast2)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "34e6a256",
   "metadata": {},
   "source": [
    "## Plot raster scans:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "546cd6fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 12\n",
    "freq_plt = 650. # in MHz\n",
    "df = freq_rast[1]-freq_rast[0]\n",
    "wf = np.where(abs(freq_rast-freq_plt)<df)[0][0]\n",
    "\n",
    "if (Ne > 1 or Nm > 1):\n",
    "    fig1, axs = plt.subplots(2,2,figsize=(12,12))\n",
    "    for i in range(0,2):\n",
    "        \n",
    "        t1_plt = Time(raster1_start[i], format='isot',scale='utc').mjd\n",
    "        t2_plt = Time(raster1_stop[i], format='isot',scale='utc').mjd\n",
    "        wplt1 = np.where((t_plt_rast1[i] >= t1_plt) & (t_plt_rast1[i] <= t2_plt))[0]\n",
    "        print(wplt1)\n",
    "        \n",
    "        plot1 = axs[i,0].scatter( (ra_rast1_deg[i][wplt1]-RA1[i])*np.cos(dec_rast1[i][wplt1]*np.pi/180.), dec_rast1[i][wplt1],\n",
    "                                s=20,c=10.*np.log10(LL_rast1[i][wplt1,wf]),cmap='viridis',vmin=72, vmax=74)\n",
    "        axs[i,0].scatter(0,DEC1[i],s=30,color='black')\n",
    "        \n",
    "        \n",
    "        t1_plt = Time(raster2_start[i], format='isot',scale='utc').mjd\n",
    "        t2_plt = Time(raster2_stop[i], format='isot',scale='utc').mjd\n",
    "        wplt2 = np.where((t_plt_rast2[i] >= t1_plt) & (t_plt_rast2[i] <= t2_plt))[0]\n",
    "\n",
    "        plot2 = axs[i,1].scatter( (ra_rast2_deg[i][wplt2]-RA2[i])*np.cos(dec_rast2[i][wplt2]*np.pi/180.), dec_rast2[i][wplt2],\n",
    "                                s=20,c=10.*np.log10(LL_rast2[i][wplt2,wf]),cmap='viridis',vmin=72, vmax=74)\n",
    "        axs[i,1].scatter(0,DEC2[i],s=30,color='black')\n",
    "        \n",
    "        for j in range(0,2):\n",
    "            axs[i,j].set_xlabel('$\\Delta$ RAcos($\\delta$) (deg.)',fontsize=fs)\n",
    "            axs[i,j].set_ylabel('Dec. (deg.)',fontsize=fs)\n",
    "            axs[i,j].tick_params(axis=\"x\", labelsize=fs)\n",
    "            axs[i,j].tick_params(axis=\"y\", labelsize=fs)\n",
    "            axs[i,j].set_xlim(5,-5)\n",
    "            axs[i,j].set(aspect='equal')\n",
    "            axs[i,j].grid()\n",
    "    \n",
    "        \n",
    "else:\n",
    "    fig1, axs = plt.subplots(1,2,figsize=(12,6))    \n",
    "    t1_plt = Time(raster1_start[0], format='isot',scale='utc').mjd\n",
    "    t2_plt = Time(raster1_stop[0], format='isot',scale='utc').mjd\n",
    "    wplt1 = np.where((t_plt_rast1[0] >= t1_plt) & (t_plt_rast1[0] <= t2_plt))\n",
    "    t1_plt = Time(raster2_start[0], format='isot',scale='utc').mjd\n",
    "    t2_plt = Time(raster2_stop[0], format='isot',scale='utc').mjd\n",
    "    wplt2 = np.where((t_plt_rast2[0] >= t1_plt) & (t_plt_rast2[0] <= t2_plt))\n",
    "    dB1 = 68\n",
    "    dB2 = 74\n",
    "    plot1 = axs[0].scatter( (ra_rast1_deg[0][wplt1]-RA1[0])*np.cos(dec_rast1[0][wplt1]*np.pi/180.), dec_rast1[0][wplt1],\n",
    "                            s=20,c=10.*np.log10(LL_rast1[0][wplt1,wf]),cmap='viridis',vmin=72, vmax=74)\n",
    "    axs[0].scatter(0,DEC1[0],s=30,color='black')\n",
    "\n",
    "    plot2 = axs[1].scatter( (ra_rast2_deg[0][wplt2]-RA2[0])*np.cos(dec_rast2[0][wplt2]*np.pi/180.), dec_rast2[0][wplt2],\n",
    "            s=20,c=10.*np.log10(LL_rast2[0][wplt2,wf]),cmap='viridis',vmin=72, vmax=74)\n",
    "    axs[1].scatter(0,DEC2,s=30,color='black')\n",
    "\n",
    "    for i in range(0,2):\n",
    "        axs[i].set_xlabel('$\\Delta$ RAcos($\\delta$) (deg.)',fontsize=fs)\n",
    "        axs[i].set_ylabel('Dec. (deg.)',fontsize=fs)\n",
    "        axs[i].tick_params(axis=\"x\", labelsize=fs)\n",
    "        axs[i].tick_params(axis=\"y\", labelsize=fs)\n",
    "        axs[i].set_xlim(5,-5)\n",
    "        axs[i].set(aspect='equal')\n",
    "        axs[i].grid()\n",
    "    \n",
    "#plt.savefig('../DVA2/PLOTS/rasters_day45.jpg')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "18e7e675",
   "metadata": {},
   "source": [
    "## Check elevation range of each raster scan:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "725693ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "ii = 2900\n",
    "print(freq_rast[ii])\n",
    "plt.plot(t_plt_rast2[0],el_rast2[0])\n",
    "#plt.plot(t_plt_rast1[0],10*np.log10(LL_rast1[0][:,ii]))\n",
    "#plt.ylim(70,90)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "499bccde",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57ee473e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d4cf66f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "3a52ba69",
   "metadata": {},
   "source": [
    "# Further binning frequency channels for select scan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "892447f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def freq_avg(step,i0):\n",
    "\n",
    "    RR_avg = np.empty([len(t_set),int(len(freq[i0:-1])/step)])\n",
    "    LL_avg = np.empty([len(t_set),int(len(freq[i0:-1])/step)])\n",
    "    reRL_avg = np.empty([len(t_set),int(len(freq[i0:-1])/step)])\n",
    "    imRL_avg = np.empty([len(t_set),int(len(freq[i0:-1])/step)])\n",
    "    freq_avg = np.empty([int(len(freq[i0:-1])/step)])\n",
    "\n",
    "    for i in range(0,int(len(freq[i0:-1])/step)-1):\n",
    "    \n",
    "        i1 = i0+i*step\n",
    "        i2 = i0+(i+1)*step\n",
    "    \n",
    "        freq_avg[i] = np.nanmean(freq[i1:i2])\n",
    "        #print(freq[i1:i2],freq_avg[i])\n",
    "    \n",
    "        RR_avg[:,i] = np.nanmedian(RR_set[:,i1:i2],axis=1)\n",
    "        LL_avg[:,i] = np.nanmedian(LL_set[:,i1:i2],axis=1)\n",
    "        reRL_avg[:,i] = np.nanmedian(reRL_set[:,i1:i2],axis=1)\n",
    "        imRL_avg[:,i] = np.nanmedian(imRL_set[:,i1:i2],axis=1)\n",
    "    \n",
    "    return RR_avg, LL_avg, reRL_avg, imRL_avg, freq_avg\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d9cb300",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "file = h5py.File('/srv/data/dva/survey_raster_noise_corr/dva_survey_phase3_day040_raster2_gain_corr.h5','r')\n",
    "freq = file['data']['beam_0']['band_SB0']['frequency'][:]/1e6\n",
    "\n",
    "dataset = file['data']['beam_0']['band_SB0']['scan_0']\n",
    "\n",
    "dec_set = dataset['metadata']['declination']\n",
    "ra_set = dataset['metadata']['right_ascension']\n",
    "el_set =dataset['metadata']['elevation']\n",
    "az_set = dataset['metadata']['azimuth']\n",
    "t_set = dataset['metadata']['utc']\n",
    "noise_set = dataset['metadata']['noise_state']\n",
    "int_time = dataset['metadata']['integration_time']\n",
    "corrupt = dataset['metadata']['corrupted']\n",
    "trim_flag = dataset['metadata']['trim_scan_flag']\n",
    "\n",
    "# Add the spectrometer power data to the corresponding arrays:\n",
    "RR_set = dataset['data'][:,0,:]\n",
    "LL_set = dataset['data'][:,1,:]\n",
    "reRL_set = dataset['data'][:,2,:]\n",
    "imRL_set = dataset['data'][:,3,::]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bee0e12",
   "metadata": {},
   "outputs": [],
   "source": [
    "step = 12\n",
    "i0 = 81\n",
    "RR_avg, LL_avg, reRL_avg, imRL_avg, freq_bin = freq_avg(step,i0)\n",
    "#print(len(freq_bin))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c2fefb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#outfiles = '/media/ordoga/15m_band1_survey/dva_phase3/survey_phase3_day54/'\n",
    "file_ex = 'dva_survey_phase3_day040_raster2_gain_corr.h5'\n",
    "outname = 'dva_survey_phase3_day040_raster2_gain_corr_1MHzMed.h5'\n",
    "\n",
    "cmd2 = 'cp '+dir_out+file_ex+' '+dir_out+outname\n",
    "\n",
    "os.system(cmd2)\n",
    "file2 = h5py.File(dir_out+outname,'r+')\n",
    "\n",
    "del file2['data']['beam_0']['band_SB0']\n",
    "\n",
    "# Create band and scan groups:\n",
    "file2['data']['beam_0'].create_group(\"band_SB0\")\n",
    "file2['data']['beam_0']['band_SB0'].create_group(f\"scan_0\")\n",
    "print('0')\n",
    "\n",
    "# Create power dataset:\n",
    "dat = np.empty((len(t_set), 4, len(freq_bin[0:-1])), dtype=float)\n",
    "file2['data']['beam_0']['band_SB0']['scan_0'].create_dataset(\"data\", data=dat)\n",
    "file2['data']['beam_0']['band_SB0']['scan_0']['data'][:,0,:] = RR_avg[:,0:-1]\n",
    "file2['data']['beam_0']['band_SB0']['scan_0']['data'][:,1,:] = LL_avg[:,0:-1]\n",
    "file2['data']['beam_0']['band_SB0']['scan_0']['data'][:,2,:] = reRL_avg[:,0:-1]\n",
    "file2['data']['beam_0']['band_SB0']['scan_0']['data'][:,3,:] = imRL_avg[:,0:-1]\n",
    "print('1')\n",
    "\n",
    "# Create metadata with timestamps and coordinates:\n",
    "metadata_content = [t_set,az_set,el_set,ra_set,dec_set,int_time,corrupt,noise_set,trim_flag]  \n",
    "#print(noise)\n",
    "col_names = [\"utc\", \"azimuth\", \"elevation\", \"right_ascension\", \"declination\", \"integration_time\",\n",
    "             \"corrupted\",\"noise_state\",\"trim_scan_flag\"]\n",
    "col_types = np.dtype({'names':col_names,'formats':[\"S32\", \"f8\", \"f8\", \"f8\", \"f8\", \"f8\", \"?\", \"i8\",\"i8\"] } )\n",
    "rec_arr = np.rec.array(metadata_content,dtype=col_types)        \n",
    "file2['data']['beam_0']['band_SB0']['scan_0'].create_dataset(\"metadata\",data=rec_arr)\n",
    "print('2')\n",
    "\n",
    "# Create frequency dataset:\n",
    "file2['data']['beam_0']['band_SB0'].create_dataset(\"frequency\",dtype=\"f8\",data=freq_bin[0:-1]*1e6)\n",
    "\n",
    "# Create polarizations dataset:\n",
    "pol_labels = [b\"ReRR\",b\"ReLL\",b\"ReRL\",b\"ImRL\"]\n",
    "file2['data']['beam_0']['band_SB0'].create_dataset(\"polarization\",dtype=\"S32\",data=pol_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6c26c95",
   "metadata": {},
   "outputs": [],
   "source": [
    "filetest = file2 = h5py.File(dir_out+outname,'r+')\n",
    "freq_test = filetest['data']['beam_0']['band_SB0']['frequency'][:]/1e6\n",
    "\n",
    "dataset_test = filetest['data']['beam_0']['band_SB0']['scan_0']\n",
    "\n",
    "dec_test = dataset_test['metadata']['declination']\n",
    "ra_test = dataset_test['metadata']['right_ascension']\n",
    "el_test =dataset_test['metadata']['elevation']\n",
    "az_test = dataset_test['metadata']['azimuth']\n",
    "t_test = dataset_test['metadata']['utc']\n",
    "t_test_mjd = Time(t_test, format='isot',scale='utc').mjd\n",
    "\n",
    "RR_test = dataset_test['data'][:,0,:]\n",
    "LL_test = dataset_test['data'][:,1,:]\n",
    "reRL_test = dataset_test['data'][:,2,:]\n",
    "imRL_test = dataset_test['data'][:,3,::]\n",
    "\n",
    "fs = 12\n",
    "freq_plt = 650. # in MHz\n",
    "df = freq_test[1]-freq_test[0]\n",
    "wf = np.where(abs(freq_test-freq_plt)<df)[0][0]\n",
    "\n",
    "fig1, axs = plt.subplots(1,1,figsize=(12,12))\n",
    "\n",
    "plt.scatter(t_test_mjd,LL_test[:,wf])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2da83c4c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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
}