\n",
" \n",
" 
\n",
"
\n",
" "
],
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#Get raw spectra and standard deviations\n",
"specs = sdfits.rawspectra(0,0)\n",
"stdevs = np.std(specs,axis=1)\n",
"\n",
"\n",
"#Organize into scan and switching state.\n",
"#There are 2 scans for the target and reference pointings, 2 calibration diode states, and 2 polarizations.\n",
"stdevs = np.reshape(stdevs, (2,-1,4))\n",
"\n",
"nrows = stdevs.shape[1]\n",
"\n",
"#Inspect the data (last integration is removed, it has issues)\n",
"for scan in range(2):\n",
" for sw_state in range(4):\n",
" plt.plot(stdevs[scan,:-1,sw_state],label=f'State {(4*scan)+sw_state}')\n",
" \n",
"plt.xlabel('Integration #')\n",
"plt.ylabel('sigma')\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "1cc1b7bf",
"metadata": {},
"source": [
"We can see that the 4 states corresponding to the OFF scan have a significant jump corresponding to the GPS L3 RFI. It does not appear to start until the 40th integration, so we will use that as our cutoff to calculate the statistics of the good data, and the thresholds to flag by."
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "ff6985a4",
"metadata": {},
"outputs": [],
"source": [
"flag_mask = np.zeros(stdevs.shape)\n",
"cutoff = 40\n",
"\n",
"mean = np.nanmean(stdevs[:,:cutoff,:],axis=1)\n",
"spread = 3 * np.nanstd(stdevs[:,:cutoff,:],axis=1)"
]
},
{
"cell_type": "markdown",
"id": "02813c6f",
"metadata": {},
"source": [
"Now we create our flagging mask of zeros and ones, where a one corresponds to a flag to be applied."
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "2c150694",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[176, 178, 185, 186, 187, 190, 191, 194, 195, 198, 199, 201, 202, 203, 206, 207, 208, 209, 210, 211, 214, 215, 216, 217, 218, 219, 222, 223, 226, 227, 230, 231, 234, 235, 236, 237, 238, 239, 240, 242, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 484, 486]\n"
]
}
],
"source": [
"flag_mask = np.zeros(stdevs.shape)\n",
"\n",
"mean = np.expand_dims(mean,axis=1)\n",
"spread = np.expand_dims(spread,axis=1)\n",
"\n",
"flag_mask[stdevs > mean+spread] = 1\n",
"flag_mask = flag_mask.flatten()\n",
"\n",
"flag_rows = np.where(flag_mask==1)[0].tolist()\n",
"print(flag_rows)"
]
},
{
"cell_type": "markdown",
"id": "f591df7a",
"metadata": {},
"source": [
"We apply the flags using the \"row\" keyword, and see that the RFI is removed, along with a drop in the exposure time to 112 seconds instead of the original 150."
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "efe82013",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exposure: 112.15101951854204 s\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "f21e85f6b8734964af0070ea08edb854",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"VBox(children=(HBox(children=(Button(description='Clear All Regions', style=ButtonStyle(), tooltip='Clear all …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"sdfits.flag(row=flag_rows)\n",
"\n",
"\n",
"ps = sdfits.getps(plnum=0, ifnum=0, fdnum=0).timeaverage()\n",
"print(f\"Exposure: {ps.meta['EXPOSURE']} s\")\n",
"\n",
"ps.plot();"
]
},
{
"cell_type": "markdown",
"id": "abc66330-e709-464e-b48d-c3c646e58c56",
"metadata": {},
"source": [
"## Flags vs. Masks\n",
"\n",
"A few words on flagging vs. masking. Flagging sets the rule.\n",
"``apply_flags`` or calibrating with ``apply_flags=True`` (the default) sets the data mask.\n",
"By default ``GBTFITSLoad`` will not read flag files. The keywords ``skipflags=`` and ``flag_vegas=`` control which flags are applied. ``flag_vegas`` can also be given as a parameter to any calibration routine, as below with ``gettp``.\n",
"\n",
"\n",
"For a given spectrum the `stats` function\n",
"reports the number of \"nan\", which should be the number of channels that were masked as a result of flagging.\n",
"\n",
"**Note:** If you change the ``Spectrum.mask`` attribute, it will be reflected in the plot, but *not* in the flags of the ``Spectrum`` if you write it out. \n",
"Such a mask change should be considered temporary; to correctly change both the mask and the flag, use a flagging operation."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "a246c300-707b-469d-ba2e-090a84c75cf7",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"12:14:06.506 I Loaded 2 FITS files\n",
"12:14:06.508 I Index loaded from .index file (44/93 columns). Missing columns (TCAL, WCS, calibration metadata, etc.) will be automatically loaded from FITS file when first accessed.\n",
"12:14:06.838 I Using TSYS column\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "3bb4b8f3a359444d9de4bf230a7ecce0",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"VBox(children=(HBox(children=(Button(description='Clear All Regions', style=ButtonStyle(), tooltip='Clear all …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"sdf1 = GBTFITSLoad(dysh_data(test=\"AGBT22A_325_15\"))\n",
"sp1 = sdf1.gettp(scan=281, ifnum=0, plnum=0, fdnum=8, flag_vegas=True).timeaverage()\n",
"sp1.plot();"
]
},
{
"cell_type": "markdown",
"id": "a12a6eaa-c8af-4e53-ac2d-d73d8167b505",
"metadata": {},
"source": [
"Looking carefully in this plot, you will see numerous \"gaps\" in the spectrum, where the so-called [VEGAS spurs](https://dysh.readthedocs.io/en/latest/reference/glossary.html#term-VEGAS-spurs) were automatically flagged. \n",
"Using the ``stats`` function you will see there are 31 in this spectrum of 1024 channels. \n",
"These are \"nan\" values in the spectrum, and the ``sp1.mask`` array has been set appropriately.\n",
"\n",
"In the corresponding AGBT22A_325_15.raw.vegas.A.flag file you will see the following:\n",
"\n",
"```\n",
"[header]\n",
"created = Tue Aug 20 09:38:32 2024\n",
"version = 1.0\n",
"created_by = sdfits\n",
"[flags]\n",
"#RECNUM,SCAN,INTNUM,PLNUM,IFNUM,FDNUM,BCHAN,ECHAN,IDSTRING\n",
"*|281|*|*|0|8|0,32,64,96,128,160,192,224,256,288,320,352,384,416,448,480,544,576,608,640,672,704,736,768,800,832,864,896,928,960,992|0,32,64,96,128,160,192,224,256,288,320,352,384,416,448,480,544,576,608,640,672,704,736,768,800,832,864,896,928,960,992|VEGAS_SPUR\n",
"*|281|*|*|0|10|0,32,64,96,128,160,192,224,256,288,320,352,384,416,448,480,544,576,608,640,672,704,736,768,800,832,864,896,928,960,992|0,32,64,96,128,160,192,224,256,288,320,352,384,416,448,480,544,576,608,640,672,704,736,768,800,832,864,896,928,960,992|VEGAS_SPUR\n",
"\n",
"\n",
"```\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "c667850a-6c48-464c-8a7d-b344694d9b8e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'mean': \n",
"