Position-Switched Data Reduction

Position-Switched Data Reduction#

This notebook shows how to use dysh to calibrate an OnOff observation. It retrieves and calibrates position-switch scans using GBTFITSLoad.getps(), which returns a ScanBlock object. OffOn observations can be reduced the same way.

You can find a copy of this tutorial as a Jupyter notebook here or download it by right clicking here and selecting “Save Link As”.

Loading Modules#

We start by loading the modules we will use for the data reduction.

For display purposes, we use the static (non-interactive) matplotlib backend in this tutorial. However, you can tell matplotlib to use the ipympl backend to enable interactive plots. This is only needed if working on jupyter lab or notebook.

# Set interactive plots in jupyter.
#%matplotlib ipympl

# These modules are required for the data reduction.
from dysh.fits.gbtfitsload import GBTFITSLoad
from astropy import units as u

# These modules are only used to download the data.
from pathlib import Path
from dysh.util.download import from_url

Data Retrieval#

Download the example SDFITS data, if necessary.

url = "http://www.gb.nrao.edu/dysh/example_data/positionswitch/data/AGBT05B_047_01/AGBT05B_047_01.raw.acs/AGBT05B_047_01.raw.acs.fits"
savepath = Path.cwd() / "data"
savepath.mkdir(exist_ok=True) # Create the data directory if it does not exist.
filename = from_url(url, savepath)

Data Loading#

Next, we use GBTFITSLoad to load the data, and then its summary method to inspect its contents.

sdfits = GBTFITSLoad(filename)

sdfits.summary()

SCAN	OBJECT	VELOCITY	PROC	PROCSEQN	RESTFREQ	DOPFREQ	# IF	# POL	# INT	# FEED	AZIMUTH	ELEVATION
51	NGC5291	4386.0	OnOff	1	1.420405	1.420405	1	2	11	1	198.3431	18.6427
52	NGC5291	4386.0	OnOff	2	1.420405	1.420405	1	2	11	1	198.9306	18.7872
53	NGC5291	4386.0	OnOff	1	1.420405	1.420405	1	2	11	1	199.3305	18.3561
54	NGC5291	4386.0	OnOff	2	1.420405	1.420405	1	2	11	1	199.9157	18.4927
55	NGC5291	4386.0	OnOff	1	1.420405	1.420405	1	2	11	1	200.3042	18.0575
56	NGC5291	4386.0	OnOff	2	1.420405	1.420405	1	2	11	1	200.8906	18.1860
57	NGC5291	4386.0	OnOff	1	1.420405	1.420405	1	2	11	1	202.3275	17.3853
58	NGC5291	4386.0	OnOff	2	1.420405	1.420405	1	2	11	1	202.9192	17.4949

Data Reduction#

Single Scan#

Next we calibrate one scan of the position switched observations. We will start with scan 51, a single spectral window and polarization.

If you don’t want to calibrate, add, calibrate=False.

ps_scan_block = sdfits.getps(scan=51, ifnum=0, plnum=0, fdnum=0)

print(f"T_sys = {ps_scan_block[0].tsys.mean():.2f}")

T_sys = 19.36

Time Averaging#

To time average the contents of a ScanBlock use its timeaverage method. Be aware that time averging will not check if the source is the same.

By default time averaging uses the following weights:

\[ \frac{T^{2}_{sys}}{\Delta\nu\Delta t} \]

with \(T_{sys}\) the system temperature, \(\Delta\nu\) the channel width and \(\Delta t\) the integration time. In dysh these are set using weights='tsys' (the default).

ta = ps_scan_block.timeaverage(weights='tsys')

Plotting#

Plot the data and use different units for the spectral axis.

ta.plot()

<dysh.plot.specplot.SpectrumPlot at 0x729f49ff0f10>

../../_images/95fd59a855ba74ab92912dcbbf0fc8f5152c6ce6a9abea9ea9df2f77e24597b5.png

Change the spectral axis units to km/s and the y-axis to mK, while also showing the spectra between 3600 and 5300 km/s, with the y-axis range between -100 and 1000 mK.

ta.plot(xaxis_unit="km/s", yaxis_unit="mK", ymin=-100, ymax=1000, xmin=3600, xmax=5300)

<dysh.plot.specplot.SpectrumPlot at 0x729f49ff0f10>

../../_images/f5cef38eba8dd1e1509f622edef8daea71a38c42f20d1cdd0e42cf7fa2573a87.png

Switch back to GHz as the spectral axis unit.

ta.plot(xaxis_unit="GHz", ymin=-100, ymax=1000, yaxis_unit="mK")

<dysh.plot.specplot.SpectrumPlot at 0x729f49ff0f10>

../../_images/d9c55e3400d223c026b010b6810f27377d71dfa35f1bb24483b1bde506c13cd4.png

Baseline Subtraction#

The following code cells show how to subtract a polynomial baseline from the data. This example uses an order 2 polynomial, and excludes the regions between 3800 and 5000 km/s, where a line is detected. The use of remove=True will remove the best fit baseline model from the spectrum.

For a polynomial model one may need to normalize the frequency axis in channel space using normalize=True, to prevent poorly conditioned fits, but this will not allow you to undo the fit.

kms = u.km/u.s
ta.baseline(model="poly", degree=2, exclude=[3800*kms,5000*kms], remove=True)

When we plot the spectrum again, it has been baseline subtracted.

ta_plt = ta.plot(xaxis_unit="km/s", yaxis_unit="mK", ymin=-500, ymax=500)

../../_images/7f55788b3aff66e61f5664dd69c071e24816680d0a0b16a9f48f8a9a284616c4.png

We can inspect the best fit baseline coefficients.

print(ta.baseline_model)

<QuantityModel Polynomial1D(2, c0=0.27098987, c1=0.02448574, c2=0.01186811, domain=(np.float64(1424816838.1210938), np.float64(1374818364.0))), input_units=Hz, return_units=K>

And save the figure.

output_dir = Path.cwd() / "output"
ta_plt.savefig(output_dir / "baselined_removed.png")

Using Selection#

The following code shows how to calibrate scan 51 using selection. At this time selection does not know about signal and reference scan pairs, so the selection must include both scans, otherwise the calibration will fail.

sdfits.select(scan=[51,52])
sdfits.selection.show()

 ID    TAG      SCAN  # SELECTED
--- --------- ------- ----------
  0 c73f66aa1 [51,52]         88

sb = sdfits.getps(plnum=0, ifnum=0, fdnum=0)

ta2 = sb.timeaverage(weights='tsys')
ta2.plot(xaxis_unit="GHz", ymin=-100, ymax=800, yaxis_unit="mK", title="PLNUM=0")

<dysh.plot.specplot.SpectrumPlot at 0x729f3f56d2d0>

../../_images/5c462710d156e542be7c0e550afe074d9e43dffe15563ba7bec5ad6f731a6c0b.png

We can calibrate the other polarization, with the scan numbers already selected.

sb = sdfits.getps(plnum=1, fdnum=0, ifnum=0)
ta3 = sb.timeaverage(weights='tsys')
ta3.plot(xaxis_unit="GHz", ymin=-100, ymax=800, yaxis_unit="mK", title="PLNUM=1")

<dysh.plot.specplot.SpectrumPlot at 0x729f3da83d30>

../../_images/8e5130cec6a686d93f37875bfbb48ba0154a23007f4903bc28a4af9dc36244ed.png

Polarization Average#

Average the polarizations and plot the result.

avg = ta2.average(ta3)
avg.plot(ymin=-100, ymax=800, yaxis_unit="mK", xaxis_unit="GHz")

<dysh.plot.specplot.SpectrumPlot at 0x729f3fb17940>

../../_images/35af3aad8010f8bf7a23ab984eb358b39fa8f8fadd3d0a2d84d9d80f95b0545d.png

All Scans#

Now we leverage the power of dysh to calibrate and time average all of the scans in the data.

We start by clearing the selection, so all of the scans are available.

sdfits.selection.clear()

Then, make a ScanBlock for spectral window 0 and polarization 0. Time average the scans in the ScanBlock into a single Spectrum and then remove a baseline.

ps_scan_block_0 = sdfits.getps(ifnum=0, plnum=0, fdnum=0)
ps_ta_0 = ps_scan_block_0.timeaverage(weights='tsys')
ps_ta_0.baseline(model="poly", degree=2, exclude=[3800*kms,5000*kms], remove=True)

Now plot and compare with the result for a single scan.

ta.plot(ymin=-400, ymax=400, yaxis_unit="mK", xaxis_unit="GHz")
ps_ta_0.plot(ymin=-400, ymax=400, yaxis_unit="mK", xaxis_unit="GHz")

<dysh.plot.specplot.SpectrumPlot at 0x729f496a3520>

../../_images/c11f07735156fffcca5e8938b80310732c8b3ac31b88f1e3a60e53b5238850ec.png

../../_images/1caaad7479cabdc8f888be6e7f7a6d487b5eea927919418338950ebe22fd59cb.png

The rms in the second Figure is almost half that of the first Figure. That is because there are four pairs of position switched scans in the data, so that results in a factor of \(\sqrt{4}\) reduced noise when we average all the data.