HI Survey

This tutorial goes through the data reduction of position switched observations of the 21 cm HI line. The data used for this tutorial is part of the New Reference Catalog of Extragalactic HI Observations by Karen O’Neil. For more details about how the observations were set up, please refer to the GBTdocs HI Position Switched (psw) Spectrum tutorial or the survey article.

Some basic information about the observations. The observations use position switching, and in some cases, the data is recorded without firing a noise diode, so that it is not possible to derive a system temperature from those records alone. In those cases, we will use observations close in time when the noise diodes where fired.

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
import numpy as np

# This module is used for custom plotting.
import matplotlib.pyplot as plt

# We will use this module to compare with published results.
import pandas as pd

# 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/hi-survey/data/AGBT04A_008_02.raw.acs/AGBT04A_008_02.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
220	3C286	0.0	OffOn	1	1.400000	1.400000	1	2	6	1	185.2806	82.0246
221	3C286	0.0	OffOn	2	1.400000	1.400000	1	2	6	1	187.2136	81.9980
222	3C286	0.0	OffOn	1	1.400000	1.400000	1	2	6	1	193.8331	81.8413
223	3C286	0.0	OffOn	2	1.400000	1.400000	1	2	6	1	195.6766	81.7788
224	3C286	0.0	OffOn	1	1.400000	1.400000	1	2	6	1	195.5182	80.2910
225	3C286	0.0	OffOn	2	1.400000	1.400000	1	2	5	1	199.9358	81.6005
226	3C286	0.0	OffOn	1	1.400000	1.400000	1	2	6	1	200.8333	80.0265
227	3C286	0.0	OffOn	2	1.400000	1.400000	1	2	6	1	205.9471	81.2609
228	B1328+254	0.0	OffOn	1	1.400000	1.400000	1	2	6	1	207.5257	73.9844
229	B1328+254	0.0	OffOn	2	1.400000	1.400000	1	2	6	1	210.9600	75.1584
230	B1345+125	0.0	OffOn	1	1.400000	1.400000	1	2	18	1	193.2738	62.0703
231	B1345+125	0.0	OffOn	2	1.400000	1.400000	1	2	18	1	195.6794	63.3244
244	B1345+125	0.0	OffOn	1	1.400000	1.400000	1	2	18	1	200.9543	60.9284
245	B1345+125	0.0	OffOn	2	1.400000	1.400000	1	2	18	1	203.5311	62.0587
246	B1345+125	0.0	OffOn	1	1.400000	1.400000	1	2	18	1	204.3408	60.3930
247	B1345+125	0.0	OffOn	2	1.400000	1.400000	1	2	18	1	206.9763	61.4655
248	B1345+125	0.0	OffOn	1	1.370000	1.370000	1	2	18	1	208.0408	59.6983
249	B1345+125	0.0	OffOn	2	1.370000	1.370000	1	2	18	1	210.7215	60.7060
250	B1345+125	0.0	Track	1	1.370000	1.370000	1	2	3	1	215.2430	59.6119
251	B1345+125	0.0	Track	1	1.370000	1.370000	1	2	3	1	215.9692	59.4175
263	U8091	213.0	OffOn	1	1.420405	1.420405	1	2	30	1	241.1339	50.6393
264	U8091	213.0	OffOn	2	1.420405	1.420405	1	2	30	1	240.9795	50.7395
265	U8091	213.0	Track	1	1.420405	1.420405	1	2	3	1	243.4807	49.8234
266	U8249	2541.0	OffOn	1	1.420405	1.420405	1	2	10	1	241.7065	50.2801
267	U8249	2541.0	OffOn	1	1.420405	1.420405	1	2	30	1	242.6892	49.6343
268	U8249	2541.0	OffOn	2	1.420405	1.420405	1	2	30	1	242.5451	49.7333
269	U8249	2541.0	Track	1	1.420405	1.420405	1	2	3	1	244.9435	48.0787
270	U8503	4676.0	OffOn	1	1.420405	1.420405	1	2	30	1	270.1741	60.2485
271	U8503	4676.0	OffOn	2	1.420405	1.420405	1	2	30	1	269.9396	60.5449
272	U8091	213.0	OffOn	1	1.420405	1.420405	1	2	30	1	253.5780	41.4420
273	U8091	213.0	OffOn	2	1.420405	1.420405	1	2	30	1	253.4586	41.5517
274	U8091	213.0	Track	1	1.420405	1.420405	1	2	3	1	255.5793	39.5541
275	U8249	2541.0	OffOn	1	1.420405	1.420405	1	2	30	1	266.5469	46.8377
276	U8249	2541.0	OffOn	2	1.420405	1.420405	1	2	30	1	266.3738	47.0388
277	U8249	2541.0	Track	1	1.420405	1.420405	1	2	3	1	267.9426	45.1863
278	U9965	4524.0	OffOn	1	1.420405	1.420405	1	2	30	1	221.5582	67.7182
279	U9965	4524.0	OffOn	2	1.420405	1.420405	1	2	30	1	221.2070	67.8070
280	U9965	4524.0	Track	1	1.420405	1.420405	1	2	3	1	222.9262	67.3657
281	U10351	891.0	OffOn	1	1.420405	1.420405	1	2	30	1	221.2955	77.4668
282	U10351	891.0	OffOn	2	1.420405	1.420405	1	2	30	1	220.4924	77.5922
283	U10351	891.0	Track	1	1.420405	1.420405	1	2	3	1	227.9032	76.2644
284	U9007	4618.0	OffOn	1	1.420405	1.420405	1	2	30	1	248.6632	38.3497
285	U9007	4618.0	OffOn	2	1.420405	1.420405	1	2	30	1	248.5549	38.4433
286	U9007	4618.0	Track	1	1.420405	1.420405	1	2	3	1	250.5574	36.6787
287	U9007	5257.0	OffOn	1	1.420405	1.420405	1	2	30	1	234.2956	48.0209
288	U9007	5257.0	OffOn	2	1.420405	1.420405	1	2	30	1	234.1669	48.0937
289	U9803	5257.0	OffOn	1	1.420405	1.420405	1	2	30	1	264.9610	57.9997
290	U9803	5257.0	OffOn	2	1.420405	1.420405	1	2	30	1	264.7276	58.2483
291	U9803	5257.0	Track	1	1.420405	1.420405	1	2	3	1	266.5221	56.2884
292	U10351	891.0	OffOn	1	1.420405	1.420405	1	2	30	1	252.1644	67.0847
293	U10351	891.0	OffOn	2	1.420405	1.420405	1	2	30	1	251.8134	67.3038
294	U10351	891.0	Track	1	1.420405	1.420405	1	2	3	1	255.3764	64.9184
295	U10629	2980.0	OffOn	1	1.420405	1.420405	1	2	30	1	233.9472	66.8551
296	U10629	2980.0	OffOn	2	1.420405	1.420405	1	2	30	1	233.5991	66.9846
297	U10629	2980.0	OffOn	1	1.420405	1.420405	1	2	30	1	238.3584	65.0614
298	U10629	2980.0	OffOn	2	1.420405	1.420405	1	2	30	1	238.0441	65.2006
299	U10629	2980.0	Track	1	1.420405	1.420405	1	2	3	1	241.4405	63.6131
300	U11017	4644.0	OffOn	1	1.420405	1.420405	1	2	30	1	235.4605	76.2074
301	U11017	4644.0	OffOn	2	1.420405	1.420405	1	2	30	1	234.7726	76.3879
302	U11017	4644.0	OffOn	1	1.420405	1.420405	1	2	30	1	241.5197	74.3514
303	U11017	4644.0	OffOn	2	1.420405	1.420405	1	2	30	1	240.9590	74.5471
304	U11017	4644.0	Track	1	1.420405	1.420405	1	2	3	1	242.6641	73.9425
305	U11017	4644.0	Track	1	1.420405	1.420405	1	2	1	1	242.9359	73.8408
306	U11017	4644.0	Track	1	1.420405	1.420405	1	2	3	1	246.0511	72.5738
307	U11461	3122.0	OffOn	1	1.420405	1.420405	1	2	30	1	168.0531	59.8975
308	U11461	3122.0	OffOn	2	1.420405	1.420405	1	2	30	1	167.8894	59.8843
309	U11461	3122.0	OffOn	1	1.420405	1.420405	1	2	30	1	173.4698	60.2443
310	U11461	3122.0	OffOn	2	1.420405	1.420405	1	2	30	1	173.3111	60.2376
311	U11461	3122.0	Track	1	1.420405	1.420405	1	2	3	1	178.1207	60.3802
312	U11578	4601.0	OffOn	1	1.420405	1.420405	1	2	30	1	156.3573	58.8036
313	U11578	4601.0	OffOn	2	1.420405	1.420405	1	2	30	1	156.1796	58.7731
314	U11578	4601.0	Track	1	1.420405	1.420405	1	2	3	1	160.4469	59.4464
315	U11578	4601.0	Track	1	1.420405	1.420405	1	2	3	1	161.0040	59.5231
316	U11627	4864.0	OffOn	1	1.420405	1.420405	1	2	30	1	157.8065	55.3372
317	U11627	4864.0	OffOn	2	1.420405	1.420405	1	2	30	1	157.6549	55.3107
318	U11627	4864.0	OffOn	1	1.420405	1.420405	1	2	30	1	162.4625	56.0655
319	U11627	4864.0	OffOn	2	1.420405	1.420405	1	2	30	1	162.3199	56.0465
320	U11627	4864.0	Track	1	1.420405	1.420405	1	2	3	1	166.7983	56.5802
321	U11992	3592.0	Track	1	1.420405	1.420405	1	2	3	1	124.4762	54.3242
322	U11992	3592.0	OffOn	1	1.420405	1.420405	1	2	30	1	125.6320	54.9128
323	U11992	3592.0	OffOn	2	1.420405	1.420405	1	2	30	1	125.4565	54.8251

There is a total of 81 scans in this data set, all used a single spectral window (IF), two polarizations (PLNUM), and a single feed (FDNUM).

Data Reduction

Noise Diode Temperature

The first step in calibrating the data is to determine the temperature of the noise diode(s). It is a known issue with GBT observations that the values provided are only accurate to \(\sim20\%\) (see e.g., Goddy et al 2020). The data we are working with contains eight scans of 3C286, a known calibrator source (see e.g., Perley & Butler 2017). We use the last two (scans 226 and 227), as the other seem to have problems, to derive the temperature of the noise diode(s).

Computing the Noise Diode Temperature

To derive the temperature of the noise diode(s) we use the gettcal method. This method is described in more detail in another tutorial. gettcal requires a zenith_opacity as input. Since this an L-band observation, we use a value of 0.08. The exact value does not have a significant impact on the results.

tcal = sdfits.gettcal(scan=226, ifnum=0, plnum=0, fdnum=0, zenith_opacity=0.08)

The return from gettcal is a TCal object, which is a child of Spectrum, so it has the same methods and properties, plus a name and snu (the flux density of the calibrator) properties. To get the mean value of the noise diode temperature over the inner \(80\%\) of the spectral window we use the TCal.get_tcal method.

tcal_0 = tcal.get_tcal()
print(tcal_0)

20.488287

Thus the temperature of the noise diode for polarization 0 is \(20.49\) K. We do the same for the second polarization.

tcal = sdfits.gettcal(scan=226, ifnum=0, plnum=1, fdnum=0, zenith_opacity=0.08)
tcal_1 = tcal.get_tcal()
print(tcal_1)

23.851044

Updating the Noise Diode Temperature

It is possible to provide the value of the noise diode temperature using the t_cal argument to the calibration methods. An alternative is to update the metadata of the GBTFITSLoad object with the new values. Here we use the second approach. First we need to determine where the derived noise temperatures are applicable. This is important as there are two choices for the noise diode temperature, high or low. We start by determining what noise diode was use for scan 226, we one we used to derive the noise diode temperature. We leverage the summary, and its add_columns and scan arguments. The type of noise diode used is stored in the “CALTYPE” column of an SDFITS. We also show the noise diode temperature stored in the metadata, in the “TCAL” column.

sdfits.summary(scan=226, add_columns="CALTYPE, TCAL")

SCAN	OBJECT	VELOCITY	PROC	PROCSEQN	RESTFREQ	DOPFREQ	# IF	# POL	# INT	# FEED	AZIMUTH	ELEVATION	CALTYPE	TCAL
226	3C286	0.0	OffOn	1	1.4	1.4	1	2	6	1	200.8333	80.0265	HIGH	21.317758

So, for scan 226 the “high” noise diode was used, and its value is not that far from what we found earlier. Where else was it used?

sdfits.summary(add_columns="CALTYPE")

SCAN	OBJECT	VELOCITY	PROC	PROCSEQN	RESTFREQ	DOPFREQ	# IF	# POL	# INT	# FEED	AZIMUTH	ELEVATION	CALTYPE
220	3C286	0.0	OffOn	1	1.400000	1.400000	1	2	6	1	185.2806	82.0246	HIGH
221	3C286	0.0	OffOn	2	1.400000	1.400000	1	2	6	1	187.2136	81.9980	HIGH
222	3C286	0.0	OffOn	1	1.400000	1.400000	1	2	6	1	193.8331	81.8413	HIGH
223	3C286	0.0	OffOn	2	1.400000	1.400000	1	2	6	1	195.6766	81.7788	HIGH
224	3C286	0.0	OffOn	1	1.400000	1.400000	1	2	6	1	195.5182	80.2910	HIGH
225	3C286	0.0	OffOn	2	1.400000	1.400000	1	2	5	1	199.9358	81.6005	HIGH
226	3C286	0.0	OffOn	1	1.400000	1.400000	1	2	6	1	200.8333	80.0265	HIGH
227	3C286	0.0	OffOn	2	1.400000	1.400000	1	2	6	1	205.9471	81.2609	HIGH
228	B1328+254	0.0	OffOn	1	1.400000	1.400000	1	2	6	1	207.5257	73.9844	HIGH
229	B1328+254	0.0	OffOn	2	1.400000	1.400000	1	2	6	1	210.9600	75.1584	HIGH
230	B1345+125	0.0	OffOn	1	1.400000	1.400000	1	2	18	1	193.2738	62.0703	HIGH
231	B1345+125	0.0	OffOn	2	1.400000	1.400000	1	2	18	1	195.6794	63.3244	HIGH
244	B1345+125	0.0	OffOn	1	1.400000	1.400000	1	2	18	1	200.9543	60.9284	HIGH
245	B1345+125	0.0	OffOn	2	1.400000	1.400000	1	2	18	1	203.5311	62.0587	HIGH
246	B1345+125	0.0	OffOn	1	1.400000	1.400000	1	2	18	1	204.3408	60.3930	HIGH
247	B1345+125	0.0	OffOn	2	1.400000	1.400000	1	2	18	1	206.9763	61.4655	HIGH
248	B1345+125	0.0	OffOn	1	1.370000	1.370000	1	2	18	1	208.0408	59.6983	HIGH
249	B1345+125	0.0	OffOn	2	1.370000	1.370000	1	2	18	1	210.7215	60.7060	HIGH
250	B1345+125	0.0	Track	1	1.370000	1.370000	1	2	3	1	215.2430	59.6119	HIGH
251	B1345+125	0.0	Track	1	1.370000	1.370000	1	2	3	1	215.9692	59.4175	LOW
263	U8091	213.0	OffOn	1	1.420405	1.420405	1	2	30	1	241.1339	50.6393	HIGH
264	U8091	213.0	OffOn	2	1.420405	1.420405	1	2	30	1	240.9795	50.7395	HIGH
265	U8091	213.0	Track	1	1.420405	1.420405	1	2	3	1	243.4807	49.8234	HIGH
266	U8249	2541.0	OffOn	1	1.420405	1.420405	1	2	10	1	241.7065	50.2801	HIGH
267	U8249	2541.0	OffOn	1	1.420405	1.420405	1	2	30	1	242.6892	49.6343	HIGH
268	U8249	2541.0	OffOn	2	1.420405	1.420405	1	2	30	1	242.5451	49.7333	HIGH
269	U8249	2541.0	Track	1	1.420405	1.420405	1	2	3	1	244.9435	48.0787	HIGH
270	U8503	4676.0	OffOn	1	1.420405	1.420405	1	2	30	1	270.1741	60.2485	HIGH
271	U8503	4676.0	OffOn	2	1.420405	1.420405	1	2	30	1	269.9396	60.5449	HIGH
272	U8091	213.0	OffOn	1	1.420405	1.420405	1	2	30	1	253.5780	41.4420	HIGH
273	U8091	213.0	OffOn	2	1.420405	1.420405	1	2	30	1	253.4586	41.5517	HIGH
274	U8091	213.0	Track	1	1.420405	1.420405	1	2	3	1	255.5793	39.5541	HIGH
275	U8249	2541.0	OffOn	1	1.420405	1.420405	1	2	30	1	266.5469	46.8377	HIGH
276	U8249	2541.0	OffOn	2	1.420405	1.420405	1	2	30	1	266.3738	47.0388	HIGH
277	U8249	2541.0	Track	1	1.420405	1.420405	1	2	3	1	267.9426	45.1863	HIGH
278	U9965	4524.0	OffOn	1	1.420405	1.420405	1	2	30	1	221.5582	67.7182	HIGH
279	U9965	4524.0	OffOn	2	1.420405	1.420405	1	2	30	1	221.2070	67.8070	HIGH
280	U9965	4524.0	Track	1	1.420405	1.420405	1	2	3	1	222.9262	67.3657	HIGH
281	U10351	891.0	OffOn	1	1.420405	1.420405	1	2	30	1	221.2955	77.4668	HIGH
282	U10351	891.0	OffOn	2	1.420405	1.420405	1	2	30	1	220.4924	77.5922	HIGH
283	U10351	891.0	Track	1	1.420405	1.420405	1	2	3	1	227.9032	76.2644	HIGH
284	U9007	4618.0	OffOn	1	1.420405	1.420405	1	2	30	1	248.6632	38.3497	HIGH
285	U9007	4618.0	OffOn	2	1.420405	1.420405	1	2	30	1	248.5549	38.4433	HIGH
286	U9007	4618.0	Track	1	1.420405	1.420405	1	2	3	1	250.5574	36.6787	HIGH
287	U9007	5257.0	OffOn	1	1.420405	1.420405	1	2	30	1	234.2956	48.0209	HIGH
288	U9007	5257.0	OffOn	2	1.420405	1.420405	1	2	30	1	234.1669	48.0937	HIGH
289	U9803	5257.0	OffOn	1	1.420405	1.420405	1	2	30	1	264.9610	57.9997	HIGH
290	U9803	5257.0	OffOn	2	1.420405	1.420405	1	2	30	1	264.7276	58.2483	HIGH
291	U9803	5257.0	Track	1	1.420405	1.420405	1	2	3	1	266.5221	56.2884	HIGH
292	U10351	891.0	OffOn	1	1.420405	1.420405	1	2	30	1	252.1644	67.0847	HIGH
293	U10351	891.0	OffOn	2	1.420405	1.420405	1	2	30	1	251.8134	67.3038	HIGH
294	U10351	891.0	Track	1	1.420405	1.420405	1	2	3	1	255.3764	64.9184	HIGH
295	U10629	2980.0	OffOn	1	1.420405	1.420405	1	2	30	1	233.9472	66.8551	HIGH
296	U10629	2980.0	OffOn	2	1.420405	1.420405	1	2	30	1	233.5991	66.9846	HIGH
297	U10629	2980.0	OffOn	1	1.420405	1.420405	1	2	30	1	238.3584	65.0614	HIGH
298	U10629	2980.0	OffOn	2	1.420405	1.420405	1	2	30	1	238.0441	65.2006	HIGH
299	U10629	2980.0	Track	1	1.420405	1.420405	1	2	3	1	241.4405	63.6131	HIGH
300	U11017	4644.0	OffOn	1	1.420405	1.420405	1	2	30	1	235.4605	76.2074	HIGH
301	U11017	4644.0	OffOn	2	1.420405	1.420405	1	2	30	1	234.7726	76.3879	HIGH
302	U11017	4644.0	OffOn	1	1.420405	1.420405	1	2	30	1	241.5197	74.3514	HIGH
303	U11017	4644.0	OffOn	2	1.420405	1.420405	1	2	30	1	240.9590	74.5471	HIGH
304	U11017	4644.0	Track	1	1.420405	1.420405	1	2	3	1	242.6641	73.9425	HIGH
305	U11017	4644.0	Track	1	1.420405	1.420405	1	2	1	1	242.9359	73.8408	HIGH
306	U11017	4644.0	Track	1	1.420405	1.420405	1	2	3	1	246.0511	72.5738	HIGH
307	U11461	3122.0	OffOn	1	1.420405	1.420405	1	2	30	1	168.0531	59.8975	HIGH
308	U11461	3122.0	OffOn	2	1.420405	1.420405	1	2	30	1	167.8894	59.8843	HIGH
309	U11461	3122.0	OffOn	1	1.420405	1.420405	1	2	30	1	173.4698	60.2443	HIGH
310	U11461	3122.0	OffOn	2	1.420405	1.420405	1	2	30	1	173.3111	60.2376	HIGH
311	U11461	3122.0	Track	1	1.420405	1.420405	1	2	3	1	178.1207	60.3802	HIGH
312	U11578	4601.0	OffOn	1	1.420405	1.420405	1	2	30	1	156.3573	58.8036	HIGH
313	U11578	4601.0	OffOn	2	1.420405	1.420405	1	2	30	1	156.1796	58.7731	HIGH
314	U11578	4601.0	Track	1	1.420405	1.420405	1	2	3	1	160.4469	59.4464	HIGH
315	U11578	4601.0	Track	1	1.420405	1.420405	1	2	3	1	161.0040	59.5231	HIGH
316	U11627	4864.0	OffOn	1	1.420405	1.420405	1	2	30	1	157.8065	55.3372	HIGH
317	U11627	4864.0	OffOn	2	1.420405	1.420405	1	2	30	1	157.6549	55.3107	HIGH
318	U11627	4864.0	OffOn	1	1.420405	1.420405	1	2	30	1	162.4625	56.0655	HIGH
319	U11627	4864.0	OffOn	2	1.420405	1.420405	1	2	30	1	162.3199	56.0465	HIGH
320	U11627	4864.0	Track	1	1.420405	1.420405	1	2	3	1	166.7983	56.5802	HIGH
321	U11992	3592.0	Track	1	1.420405	1.420405	1	2	3	1	124.4762	54.3242	HIGH
322	U11992	3592.0	OffOn	1	1.420405	1.420405	1	2	30	1	125.6320	54.9128	HIGH
323	U11992	3592.0	OffOn	2	1.420405	1.420405	1	2	30	1	125.4565	54.8251	HIGH

All, but a single scan (which we won’t use in this tutorial), used the “high” noise diode.

Now we update all of the rows, with the corresponding noise diode temperature. We have to update both polarizations separately. In this example we only have one spectral window and one beam, but if we had more, we would also need to separate those cases, as each beam can have a different noise diode and the temperature of the noise diode is a function of frequency.

To access the metadata, we use the index property of the GBTFITSLoad object, which is a pandas.DataFrame object. Then we select only those rows where “PLNUM” is equal to the value we are interested in (e.g., sdfits["PLNUM"] == 0).

sdfits.index().loc[sdfits["PLNUM"] == 0, "TCAL"] = tcal_0
sdfits.index().loc[sdfits["PLNUM"] == 1, "TCAL"] = tcal_1

Now we show how to double check that the values where updated:

sdfits["TCAL"][ (sdfits["SCAN"] == 226) & (sdfits["PLNUM"] == 0) ], tcal_0

(142    20.488287
 143    20.488287
 146    20.488287
 147    20.488287
 150    20.488287
 151    20.488287
 154    20.488287
 155    20.488287
 158    20.488287
 159    20.488287
 162    20.488287
 163    20.488287
 Name: TCAL, dtype: float32,
 np.float32(20.488287))

sdfits["TCAL"][ (sdfits["SCAN"] == 226) & (sdfits["PLNUM"] == 1) ], tcal_1

(140    23.851044
 141    23.851044
 144    23.851044
 145    23.851044
 148    23.851044
 149    23.851044
 152    23.851044
 153    23.851044
 156    23.851044
 157    23.851044
 160    23.851044
 161    23.851044
 Name: TCAL, dtype: float32,
 np.float32(23.851044))

The above shows that the values were successfully updated.

Now we can proceed calibrating the data.

Single On/Off Pair

We will start by reducing data for a single pair of position switched scans, which used the noise diodes. We will use scan 270. First, we calibrate the data for a single polarization, plnum=0. We use the getps method of GBTFITSLoad, which returns a ScanBlock. Since we are calibrating a single pair of position switched scans, the use of a ScanBlock won’t be evident, but we will see it when we calibrate multiple pairs of observations.

pssb0 = sdfits.getps(scan=270, plnum=0, ifnum=0, fdnum=0)
pssb0

ScanBlock([<dysh.spectra.scan.PSScan at 0x7b2c6ad12890>])

The return is a ScanBlock with a single PSScan in it. We can extract information from the observations by querying the different attributes of the PSScan, like the system temperature in K (tsys), or exposure time in seconds (exposure).

pssb0[0].tsys

array([25.5281268 , 25.53413116, 25.54545267, 25.51954586, 25.53610295,
       25.54238206, 25.52621493, 25.51971548, 25.50964588, 25.54143205,
       25.54019773, 25.53505961, 25.52427646, 25.53010434, 25.52924467,
       25.52770188, 25.52649021, 25.52654657, 25.5655992 , 25.54277555,
       25.53489685, 25.5230316 , 25.52542424, 25.52569588, 25.53033634,
       25.52895463, 25.52776397, 25.56483491, 25.52837916, 25.52476141])

pssb0[0].exposure

array([4.77948809, 4.77948809, 4.77948809, 4.77948809, 4.77948809,
       4.77948809, 4.77948809, 4.77948809, 4.77948809, 4.77948797,
       4.77948809, 4.77948809, 4.77948809, 4.77948809, 4.77948809,
       4.77948809, 4.77948809, 4.77948809, 4.77948785, 4.77948809,
       4.77948809, 4.77948809, 4.77948809, 4.77948809, 4.77948809,
       4.77948809, 4.77948809, 4.77948797, 4.77948797, 4.77948809])

pssb0[0].scan

271

pssb0[0].plnum

0

pssb0[0].ifnum

0

Notice that the PSScan says it has a scan number (scan) of 271. This is because dysh can tell that the on-source observation has a scan number of 271, and the off-source observation is in scan 270.

Inspecting Individual Integrations

If we want to have a look at the calibrated data, integration by integration, we can use the _calibrated attribute of the PSScan. This returns an array with rows corresponding to the integrations, and columns to the channel number.

pssb0[0]._calibrated

masked_array(
  data=[[-47.36035713079021, -16.303216080635575, 55.87725647889933, ...,
         -0.05286681449235954, 0.06121151776612735, 0.5152755476741582],
        [22.92745338310424, -30.885218208590764, -0.18564940148060624,
         ..., 0.2704048064615289, 0.6260993917462164, 1.6133043286916484],
        [1.0302216877306447, -19.202130901048342, 5.70289052057268, ...,
         0.18584532523032254, 0.9015370330507558, 0.3343565898159296],
        ...,
        [44.34781632848251, -26.495462490875568, -21.75212213155363, ...,
         -0.832781595564901, -0.2861877333504582, -2.0186235862210133],
        [-139.51144023716225, 44.00895847258403, -46.75392773723822, ...,
         0.023746246167524878, 0.5980652591112108, 0.4328228599070799],
        [-12.40044937415493, 45.724382854502814, -16.63166750837949, ...,
         -0.801084755229115, -0.5051024022592003, 0.18379687549419707]],
  mask=[[False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        ...,
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False]],
  fill_value=1e+20)

We can also retrieve the calibrated integrations as Spectrum objects using the calibrated method.

pssb0_int0 = pssb0[0].getspec(0)
pssb0_int0

<Spectrum(flux=[-47.36035713079021 ... 0.5152755476741582] K (shape=(32768,), mean=0.00836 K); spectral_axis=<SpectralAxis 
   (observer: <ITRS Coordinate (obstime=2004-04-22T07:08:15.000, location=(0.0, 0.0, 0.0) km): (x, y, z) in m
                  (882593.9465029, -4924896.36541728, 3943748.74743984)
               (v_x, v_y, v_z) in km / s
                  (0., 0., 0.)>
    target: <SkyCoord (FK5: equinox=J2000.000): (ra, dec, distance) in (deg, deg, kpc)
                (202.68562135, 32.76085866, 1000000.)
             (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
                (0., 0., 4676.)>
    observer to target (computed from above):
      radial_velocity=4686.132814909312 km / s
      redshift=0.0157553569969906
    doppler_rest=1420405400.0 Hz
    doppler_convention=optical)
  [1.39229407e+09 1.39229445e+09 1.39229483e+09 ... 1.40479293e+09
 1.40479331e+09 1.40479369e+09] Hz> (length=32768))>

Spectrum objects have a variety of methods, such as plot, smooth, and baseline. Here we use plot to look a the data.

pssb0_int0.plot()

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

The y-axis can be adjusted during the call to plot, through the ymin and ymax arguments.

pssb0_int0.plot(ymin=-5, ymax=5)

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

Since this is a single integration, there’s not much to see. Let’s work on a time average now.

Time Averaging Integrations

Time averaging can be done using the timeaverage method of a Scan or ScanBlock. 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).

timeaverage will return a Spectrum object.

ps0_spec = pssb0.timeaverage()
ps0_spec

<Spectrum(flux=[-19.314770735049287 ... 0.05942941726178433] K (shape=(32768,), mean=0.14476 K); spectral_axis=<SpectralAxis 
   (observer: <ITRS Coordinate (obstime=2004-04-22T07:08:15.000, location=(0.0, 0.0, 0.0) km): (x, y, z) in m
                  (882593.9465029, -4924896.36541728, 3943748.74743984)
               (v_x, v_y, v_z) in km / s
                  (0., 0., 0.)>
    target: <SkyCoord (FK5: equinox=J2000.000): (ra, dec, distance) in (deg, deg, kpc)
                (202.68562135, 32.76085866, 1000000.)
             (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
                (0., 0., 4676.)>
    observer to target (computed from above):
      radial_velocity=4686.132814909312 km / s
      redshift=0.0157553569969906
    doppler_rest=1420405400.0 Hz
    doppler_convention=optical)
  [1.39229407e+09 1.39229445e+09 1.39229483e+09 ... 1.40479293e+09
 1.40479331e+09 1.40479369e+09] Hz> (length=32768))>

ps0_spec.plot(ymin=-5, ymax=5)

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

The noise is lower, and there are hints of a signal. Let’s smooth the data to further reduce the noise.

Smoothing

Smoothing is done with the smooth method. By default it decimates the spectrum, so it only retains independent samples. In this case we smooth using a Gaussian kernel with a width of 16 channels.

ps0_spec_smo = ps0_spec.smooth("gauss", 16)

ps0_spec_smo.plot(ymin=-5, ymax=5)

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

We have to zoom in further to see the signal. We also limit the x-axis range.

ps0_spec_smo.plot(ymin=-0.2, ymax=0.5, xmin=1.397e9, xmax=1.400e9)

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

Polarization Averaging

While inspecting the data we saw that there are two polarizations. We can average them together to further reduce the noise by a factor \(\sqrt{2}\). The second polarization can be calibrated following the above steps, but setting plnum=1. Here we also demonstrate the use of chaining to do the data reduction. This refers to using multiple commands in a chain, like

ps1_spec_smo = sdfits.getps(scan=270, plnum=1, ifnum=0, fdnum=0).timeaverage().smooth("gauss", 16)

ps1_spec_smo.plot(ymin=-0.2, ymax=0.5, xmin=1.397e9, xmax=1.400e9)

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

Now we average the smoothed spectra for both polarizations using the average method.

ps_spec_smo = ps0_spec_smo.average([ps1_spec_smo])

ps_spec_smo.plot(ymin=-0.2, ymax=0.5, xmin=1.397e9, xmax=1.400e9)

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

Statistics

Now we will compare the noise properties of the spectra. For this we leverage the ability to slice spectra. First, we replot the spectra over the whole x-range to find a good frequency range where to compute statistics.

ps_spec_smo.plotter.reset()
ps_spec_smo.plot(ymin=-0.2, ymax=0.5)

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

We use the range between 1394 MHz and 1398 MHz, and then compute the statistics over this range using the stats method of a Spectrum.

s = slice(1394*u.MHz, 1398*u.MHz)
ps_spec_smo[s].stats()

{'mean': <Quantity 0.01366711 K>,
 'median': <Quantity 0.0138482 K>,
 'rms': <Quantity 0.02056197 K>,
 'min': <Quantity -0.04982653 K>,
 'max': <Quantity 0.0861573 K>,
 'npt': 656,
 'nan': np.int64(0)}

Now for the individual polarizations.

ps0_spec_smo[s].stats(), ps1_spec_smo[s].stats()

({'mean': <Quantity 0.00894887 K>,
  'median': <Quantity 0.00802389 K>,
  'rms': <Quantity 0.0236396 K>,
  'min': <Quantity -0.07328762 K>,
  'max': <Quantity 0.07712214 K>,
  'npt': 656,
  'nan': np.int64(0)},
 {'mean': <Quantity 0.02042203 K>,
  'median': <Quantity 0.01898732 K>,
  'rms': <Quantity 0.02968233 K>,
  'min': <Quantity -0.07752056 K>,
  'max': <Quantity 0.13824821 K>,
  'npt': 656,
  'nan': np.int64(0)})

The individual polarizations had an rms of \(\approx0.0255\) K, and the average an rms of \(0.02\) K. Thus, the average has a noise a factor of \(0.9\sqrt{2}\) lower than the individual polarizations. That is \(10\%\) higher than expected.

Baseline Subtraction

Now we will subtract a baseline from the averaged spectrum. For this we use the baseline method. It is important to use a range of frequencies that will not bias the baseline fit. We exclude the range at the low-frequency end of the spectral window, up to 1394 MHz, and the range that contains the spectral line, between 1398 and 1400 MHz. In this case we use a polynomial of order 1 as our baseline model.

exclude = [(1*u.GHz,1.394*u.GHz),(1.398*u.GHz,1.4*u.GHz)]
ps_spec_smo.baseline(1, model="poly", exclude=exclude, remove=True)
ps_spec_smo.plot(ymin=-0.2, ymax=0.5, xmin=1.397e9, xmax=1.400e9)

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

ps_spec_smo[s].stats()

{'mean': <Quantity 0.00027589 K>,
 'median': <Quantity 0.00033945 K>,
 'rms': <Quantity 0.02056173 K>,
 'min': <Quantity -0.0631958 K>,
 'max': <Quantity 0.07325353 K>,
 'npt': 656,
 'nan': np.int64(0)}

ps_spec_smo.plot(ymax=0.25, ymin=-0.1, xaxis_unit="km/s")

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

The mean and median are closer to zero now.

Multiple On/Off Pairs

Now that we understand how to process a single pair of On/Off scans, we proceed to calibrate a bunch of them. In dysh this can be accomplished by either, giving a list of scans to the calibration routines, or by selecting the scans based on another property of the data.

Using a List of Scans

First we need to figure out all of the scans for a particular source. For U8503, there is a single pair of On/Off scans, so we need to use a different source. We use 3C286, for which we have scans 220, 221, 222, 223, 224, 225, 226, 227 using OffOn with the noise diodes.

We could have figured out which scans using:

scan_list = list(set(sdfits["SCAN"][(sdfits["OBJECT"] == "3C286") & \
                     (sdfits["PROC"] == "OffOn") & (sdfits["CAL"] == "T")]))
sorted(scan_list)

[220, 221, 222, 223, 224, 225, 226, 227]

ps0_obj = sdfits.getps(scan=scan_list, plnum=0, ifnum=0, fdnum=0).timeaverage()
ps1_obj = sdfits.getps(scan=scan_list, plnum=1, ifnum=0, fdnum=0).timeaverage()
ps_obj = ps0_obj.average(ps1_obj)

ps_obj.plot()

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

Since 3C286 is a continuum source, we can see Galactic HI absorption against the continuum.

ps_obj.plot(xmin=1.420e9, xmax=1.421e9, ymin=15, ymax=25)

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

Using Selection

We can do the same by using the selection method before calling getps.

sdfits.select(object="3C286", proc="OffOn")
ps0_obj_b = sdfits.getps(plnum=0, ifnum=0, fdnum=0).timeaverage()
ps1_obj_b = sdfits.getps(plnum=1, ifnum=0, fdnum=0).timeaverage()
ps_obj_b = ps0_obj_b.average(ps1_obj_b)

ps_obj_b.plot(xmin=1.420e9, xmax=1.421e9, ymin=15, ymax=25)

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

Once you are done calibrating, you should clear the selection to have access to all the data again.

sdfits.selection.clear()

Using Arguments

We can accomplish the same by specifying the object and procedure type (proc) during the call to getps. Any of the calibration methods (getps, getfs, getnod, getsigref, subbeamnod, ) can take as argument any of the columns of an SDFITS file. When used this way, the calibration routine will only use the data that satisfies the conditions, so for example, if we use sdfits.getps(object="3C286", plnum=0, ifnum=0, fdnum=0) the calibration routine will only use data that has the column “OBJECT” equal to “3C286”.

ps0_obj_c = sdfits.getps(plnum=0, ifnum=0, fdnum=0, object="3C286", proc="OffOn").timeaverage()
ps1_obj_c = sdfits.getps(plnum=1, ifnum=0, fdnum=0, object="3C286", proc="OffOn").timeaverage()
ps_obj_c = ps0_obj_c.average(ps1_obj_c)
ps_obj_c.plot(xmin=1.420e9, xmax=1.421e9, ymin=15, ymax=25)

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

Calibration Without Noise Diodes

Most of the scans in this example did not use the noise diode. In this case we need to provide a value for the system temperature so that the data can be calibrated. For this particular observation, there is one Track observing procedure associated with the OffOn pairs. For these Track observations, the noise diode was fired, so we can use them to figure out the system temperature.

We will work on observations of U11627. For this source the Track scan is 320, and the OffOn pairs are in scans 316, 317, 318 and 319.

First we use gettp to figure out the system temperature from the Track scan.

tp0 = sdfits.gettp(scan=320, plnum=0, ifnum=0, fdnum=0).timeaverage()
tp1 = sdfits.gettp(scan=320, plnum=1, ifnum=0, fdnum=0).timeaverage()

The system temperature is stored in the meta dictionary of each Spectrum.

print(f"System temperature for plnum={tp0.meta['PLNUM']}: {tp0.meta['TSYS']:.2f} K")
print(f"System temperature for plnum={tp1.meta['PLNUM']}: {tp1.meta['TSYS']:.2f} K")

System temperature for plnum=0: 26.33 K
System temperature for plnum=1: 31.30 K

Now we use these values to calibrate the data. The system temperature is provided for the calibration methods through the t_sys argument. It is assumed to be in K.

sdfits.select(object="U11627", proc="OffOn")
ps0_wtsys = sdfits.getps(plnum=0, ifnum=0, fdnum=0, t_sys=tp0.meta['TSYS']).timeaverage()
ps1_wtsys = sdfits.getps(plnum=1, ifnum=0, fdnum=0, t_sys=tp1.meta['TSYS']).timeaverage()
ps_wtsys = ps0_wtsys.average(ps1_wtsys)

ps_wtsys.plot(ymin=-0.5, ymax=0.1)

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

Now we smooth and remove a baseline.

ps_wtsys_smo = ps_wtsys.smooth("gauss", 16)
ps_wtsys_smo.baseline(1, model="poly", exclude=[(1*u.GHz,1.393*u.GHz),(1.397*u.GHz,1.399*u.GHz)], remove=True)

ps_wtsys_smo.plot(ymin=-0.2, ymax=0.2)

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

sdfits.selection.clear()

Combining Off Spectra

In some situations we would like to have more flexibility during calibration to specify what will be the Off or reference spectrum used for calibration. In these cases we can use the GBTFITSLoad.getsigref function. This function takes as inputs a scan number, or list of them, to be used as On, and a reference spectrum, or scan number (only one in this case). Here we will combine the two Off source spectra from the previous calibration before calibrating the data.

We use gettp to produce the reference spectrum.

tp_ref0 = sdfits.gettp(scan=[316,318], plnum=0, ifnum=0, fdnum=0, t_sys=tp0.meta['TSYS']).timeaverage()
tp_ref1 = sdfits.gettp(scan=[316,318], plnum=1, ifnum=0, fdnum=0, t_sys=tp1.meta['TSYS']).timeaverage()

Now use getsigref to do the calibration. Since we specified the system temperature in the previous call to gettp, we do not need to provide it again.

sdfits.select(object="U11627", proc="OffOn")
ps0_wtsys_tpr = sdfits.getsigref(scan=[317,319], ref=tp_ref0, plnum=0, ifnum=0, fdnum=0).timeaverage()
ps1_wtsys_tpr = sdfits.getsigref(scan=[317,319], ref=tp_ref1, plnum=1, ifnum=0, fdnum=0).timeaverage()

Average both polarizations, smooth and remove a baseline like before so we can compare the results.

ps_wtsys_tpr = ps0_wtsys_tpr.average(ps1_wtsys_tpr)
ps_wtsys_tpr_smo = ps_wtsys_tpr.smooth("gauss", 16)
ps_wtsys_tpr_smo.baseline(1, model="poly", 
                          exclude=[(1*u.GHz,1.393*u.GHz),(1.397*u.GHz,1.399*u.GHz)], 
                          remove=True)

ps_wtsys_tpr_smo.plot(ymin=-0.2, ymax=0.2)

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

Now compare the noise in the end products.

s = slice(1.393*u.GHz, 1.396*u.GHz)
rms_tpr = ps_wtsys_tpr_smo[s].stats()["rms"]
rms = ps_wtsys_smo[s].stats()["rms"]
print(f"Ratio of rms: {rms_tpr/rms}")

Ratio of rms: 0.9999915517623001

In this case, using a combined reference spectrum improved the noise by an insignificant amount, \(\approx0.02\%\).

Working in Flux Units

So far we have calibrated the data to antenna temperature units, but it is also possible to work on flux units. To do so, we must provide a zenith opacity value and specify that we want the data in Jy, using the zenith_opacity and units parameters. Since this data is observed at 1.4 GHz, we use a small value for the zenith opacity, 0.08. For more details on how to find the zenith opacity for your observations, please see this guide. We also repeat the previous calibration steps.

ps0_wtsys_tpr = sdfits.getsigref(scan=[317,319], ref=tp_ref0, plnum=0, ifnum=0, fdnum=0, 
                                 units="flux", zenith_opacity=0.08).timeaverage()
ps1_wtsys_tpr = sdfits.getsigref(scan=[317,319], ref=tp_ref1, plnum=1, ifnum=0, fdnum=0, 
                                 units="flux", zenith_opacity=0.08).timeaverage()
ps_wtsys_tpr = ps0_wtsys_tpr.average(ps1_wtsys_tpr)
ps_wtsys_tpr_smo = ps_wtsys_tpr.smooth("gauss", 16)
ps_wtsys_tpr_smo.baseline(1, model="poly", 
                          exclude=[(1*u.GHz,1.393*u.GHz),(1.397*u.GHz,1.399*u.GHz)], 
                          remove=True)
ps_wtsys_tpr_smo.plot(ymin=-30, ymax=80, yaxis_unit="mJy")

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

Measuring Line Properties

Spectrum objects provide a convenience function for analysis of HI profiles based on the Curve of Growth (CoG) method by Yu et al. (2020), through the cog method. The implementation of this method can retrieve line parameters in a fully automated way, however, this is likely to work only in high signal-to-noise cases (>10). In low signal-to-noise cases, additional inputs are required to aid the method. In particular, it helps to provide a good estimate of the central velocity of the line using the vc parameter and/or restricting the range over which the method is applied, either by slicing the Spectrum or using the bchan and echan parameters.

First we apply the method blindly, ignoring the edge channels.

line_props = ps_wtsys_tpr_smo[60:-60].cog()
line_props

{'flux': <Quantity 5.07823274 Jy km / s>,
 'flux_std': <Quantity 0.34899029 Jy km / s>,
 'flux_r': <Quantity 0.5966902 Jy km / s>,
 'flux_r_std': <Quantity 0.19097394 Jy km / s>,
 'flux_b': <Quantity 4.40853276 Jy km / s>,
 'flux_b_std': <Quantity 0.23415411 Jy km / s>,
 'width': {0.25: <Quantity 49.23181793 km / s>,
  0.65: <Quantity 155.67900098 km / s>,
  0.75: <Quantity 182.29080043 km / s>,
  0.85: <Quantity 208.90260197 km / s>,
  0.95: <Quantity 246.15912821 km / s>},
 'width_std': {0.25: <Quantity 5.34507999 km / s>,
  0.65: <Quantity 10.75795749 km / s>,
  0.75: <Quantity 13.43018918 km / s>,
  0.85: <Quantity 34.65836073 km / s>,
  0.95: <Quantity 18.7902007 km / s>},
 'A_F': np.float64(7.388310972733408),
 'A_C': np.float64(3.56302916592011),
 'C_V': np.float64(4.2432437142336825),
 'rms': <Quantity 0.00695291 Jy>,
 'bchan': np.int64(738),
 'echan': np.int64(1108),
 'vel': <Quantity 4891.75339314 km / s>,
 'vel_std': <Quantity 600.65656517 km / s>,
 'vframe': 'itrs',
 'doppler_convention': 'optical'}

Plot the results along with the central velocity and the line width that encompasses 95% of the total flux.

plt.figure()
plt.plot(ps_wtsys_tpr_smo[60:].spectral_axis.to("km/s"), ps_wtsys_tpr_smo[60:].flux)
plt.axvline(line_props["vel"].value, c="r", alpha=0.5, ls=":")
plt.axvline((line_props["vel"] - line_props["width"][0.95]/2).value, c="g")
plt.axvline((line_props["vel"] + line_props["width"][0.95]/2).value, c="g")
plt.xlim(4500, 5200)
plt.xlabel("Velocity (km/s)")
plt.ylabel("Flux (Jy)");

The results are not great. Use CoG again, but providing a range of channels through the bchan and echan parameters (these can only be channel numbers).

line_props_wrange = ps_wtsys_tpr_smo[60:-60].cog(bchan=line_props["bchan"], echan=line_props["echan"])
line_props_wrange

{'flux': <Quantity 4.72722808 Jy km / s>,
 'flux_std': <Quantity 0.05107974 Jy km / s>,
 'flux_r': <Quantity 2.60136065 Jy km / s>,
 'flux_r_std': <Quantity 0.07879069 Jy km / s>,
 'flux_b': <Quantity 2.12677119 Jy km / s>,
 'flux_b_std': <Quantity 0.06786515 Jy km / s>,
 'width': {0.25: <Quantity 33.25196154 km / s>,
  0.65: <Quantity 81.13478735 km / s>,
  0.75: <Quantity 91.77541565 km / s>,
  0.85: <Quantity 107.73635843 km / s>,
  0.95: <Quantity 158.27934713 km / s>},
 'width_std': {0.25: <Quantity 1.36866178 km / s>,
  0.65: <Quantity 1.55611853 km / s>,
  0.75: <Quantity 1.61419184 km / s>,
  0.85: <Quantity 2.87004329 km / s>,
  0.95: <Quantity 8.13591869 km / s>},
 'A_F': np.float64(1.2231502225104227),
 'A_C': np.float64(1.0622926871036238),
 'C_V': np.float64(3.2400000916649785),
 'rms': <Quantity 0.00695291 Jy>,
 'bchan': np.int64(738),
 'echan': np.int64(1108),
 'vel': <Quantity 4833.02405876 km / s>,
 'vel_std': <Quantity 248.51383112 km / s>,
 'vframe': 'itrs',
 'doppler_convention': 'optical'}

line_props_wrange = ps_wtsys_tpr_smo[60:-60].cog(bchan=750, echan=1100)
line_props_wrange

{'flux': <Quantity 4.72159151 Jy km / s>,
 'flux_std': <Quantity 0.05255238 Jy km / s>,
 'flux_r': <Quantity 2.54269603 Jy km / s>,
 'flux_r_std': <Quantity 0.09197529 Jy km / s>,
 'flux_b': <Quantity 2.17711324 Jy km / s>,
 'flux_b_std': <Quantity 0.07543114 Jy km / s>,
 'width': {0.25: <Quantity 33.2525423 km / s>,
  0.65: <Quantity 81.1362044 km / s>,
  0.75: <Quantity 91.77701855 km / s>,
  0.85: <Quantity 107.73824009 km / s>,
  0.95: <Quantity 155.62190763 km / s>},
 'width_std': {0.25: <Quantity 2.68090569 km / s>,
  0.65: <Quantity 1.55624457 km / s>,
  0.75: <Quantity 2.81406913 km / s>,
  0.85: <Quantity 2.87009341 km / s>,
  0.95: <Quantity 10.75401215 km / s>},
 'A_F': np.float64(1.1679208853011962),
 'A_C': np.float64(1.08315990940368),
 'C_V': np.float64(3.2400000916665244),
 'rms': <Quantity 0.00693929 Jy>,
 'bchan': np.int64(750),
 'echan': np.int64(1100),
 'vel': <Quantity 4835.71843257 km / s>,
 'vel_std': <Quantity 244.60438983 km / s>,
 'vframe': 'itrs',
 'doppler_convention': 'optical'}

plt.figure()
plt.plot(ps_wtsys_tpr_smo[60:].spectral_axis.to("km/s"), ps_wtsys_tpr_smo[60:].flux)
plt.axvline(line_props_wrange["vel"].value, c="r", alpha=0.5, ls=":")
plt.axvline((line_props_wrange["vel"] - line_props_wrange["width"][0.95]/2).value, c="g")
plt.axvline((line_props_wrange["vel"] + line_props_wrange["width"][0.95]/2).value, c="g")
plt.xlim(4500, 5200)
plt.xlabel("Velocity (km/s)")
plt.ylabel("Flux (Jy)");

Much better!

Putting it all Together

Now we can use what we have learned to process all of the observations. We will loop over all the objects, calibrating the data. For the system temperature, we use the last Track observation associated with a given object. Since there may have been issues with the previous Track observations if they had to be repeated (for example, for U11578 the first Track observation did not fire the noise diode). We process “U8503” individually, because there is no “Track” scan for it but the OffOn observations did use the noise diode.

We put all the steps in a function, process. This makes it easier to reuse the code, and modify it.

def process(sdfits, object, results, track=True):
    """
    Function to calibrate the AGBT04A_008_02 observations.
    This function was heavily tailored to work with this
    observations and there is no guarantee it would work with
    other data.
    """
    
    o = object

    if track:
        # Use only the last Track observation for every object.
        tp0 = sdfits.gettp(ifnum=0, plnum=0, fdnum=0, object=o, proc="Track")[-1].timeaverage()
        tp1 = sdfits.gettp(ifnum=0, plnum=1, fdnum=0, object=o, proc="Track")[-1].timeaverage()
    
        # Calibrate using the system temperature of the Track scan.
        ps0 = sdfits.getps(ifnum=0, plnum=0, fdnum=0, object=o, proc="OffOn", 
                           t_sys=tp0.meta["TSYS"], units="flux", zenith_opacity=0.08).timeaverage()
        ps1 = sdfits.getps(ifnum=0, plnum=1, fdnum=0, object=o, proc="OffOn", 
                           t_sys=tp1.meta["TSYS"], units="flux", zenith_opacity=0.08).timeaverage()

    else:
        # Calibrate computing the system temperature from the Off scan.
        ps0 = sdfits.getps(ifnum=0, plnum=0, fdnum=0, object=o, proc="OffOn", 
                           units="flux", zenith_opacity=0.08).timeaverage()
        ps1 = sdfits.getps(ifnum=0, plnum=1, fdnum=0, object=o, proc="OffOn", 
                           units="flux", zenith_opacity=0.08).timeaverage()

    # Average polarizations.
    ps = ps0.average(ps1)

    # Smooth.
    ps_smo = ps.smooth("gauss", 16)

    # Determine if the Galactic HI line is present, and if so,
    # ignore it during the baseline fit.
    idx0 = np.argmin(abs(ps_smo.spectral_axis - 1.420*u.GHz))
    idxf = np.argmin(abs(ps_smo.spectral_axis - 1.421*u.GHz))
    idx0,idxf = np.sort([idx0,idxf])
    if idx0 == 0 or idx0 == len(ps_smo.data) - 1 or idxf == 0 or idxf == len(ps_smo.data) - 1:
        exclude=[(0,100),(500,1250),(2047-100,2047)]
    else:
        exclude=[(0,100),(500,1250),(idx0,idxf),(2047-100,2047)]

    # Baseline subtraction.
    ps_smo.baseline(degree=1, model="poly", exclude=exclude, remove=True)

    # Measure line properties using Curve of Growth.
    # Ignore 100 channels in each edge, and compute the
    # CoG over the inner (750,1250) channels.
    cog = ps_smo[100:-100].cog(bchan=750, echan=1250, width_frac=[0.95])

    # Save the measured line properties.
    # In particular, the object name, flux and its width.
    results["name"].append(o.replace("U", "UGC "))
    for k in ["flux", "flux_error"]:
        results[k].append(cog[k.replace("error", "std")].to("Jy km/s").value)
    for k in ["width", "width_error"]:
        results[k].append(cog[k.replace("error", "std")][0.95].to("km/s").value)

measured = {"name": [],
            "flux": [],
            "flux_error": [],
            "width": [],
            "width_error": [],
           }

sdfits.selection.clear()
# Start at object 4 since the previous ones were calibration observations.
for o in sdfits.udata("OBJECT")[4:]:

    if o ==  "U8503":
        track = False
    else:
        track = True
    
    process(sdfits, o, measured, track=track)

Compare with the Literature

We can compare the results obtained here with those listed by the author in this page. The results are also provided as a text table. We download this table and parse its contents. We save the flux and \(W_{20}\) (the width of the line at 20% of the peak flux) values.

table_file = from_url("https://greenbankobservatory.org/~koneil/HIsurvey/results_all.dat", savepath)

data = {"name": [],
        "flux": [],
        "flux_error": [],
        "width": [],
        "width_error": []
       }

# Open the file.
with open(table_file) as f:
    lines = f.readlines()
    # Loop over lines extracting the data we are interested in.
    # The column numbers are provided in the header of the text file.
    for line in lines:
        if line.lstrip().startswith("\\"):
            continue
        data["name"].append(" ".join(line[1:14].strip().split()))
        try:
            data["flux"].append(float(line[49:53].strip()))
        except ValueError:
            data["flux"].append(np.nan)
        data["flux_error"].append(float(line[54:57].strip()))
        try:
            data["width"].append(float(line[58:61].strip()))
        except ValueError:
            data["width"].append(np.nan)
        data["width_error"].append(float(line[62:64].strip()))

Now create pandas.DataFrame objects to manipulate the measured and literature results. We use DataFrame as it provides a convenient way of handling the data. While loading the data, we sort it by “name”. We will also remove sources that do not appear in both data sets.

# Convert to DataFrame and sort.
df_obs = pd.DataFrame.from_dict(measured).sort_values(by="name")
df_lit = pd.DataFrame.from_dict(data).sort_values(by="name")

# Find sources present in both data sets.
shared_names = set(df_obs["name"]) & set(df_lit["name"])

# Keep only the sources found above.
df_obs_s = df_obs[df_obs["name"].isin(shared_names)]
df_lit_s = df_lit[df_lit["name"].isin(shared_names)]

For example, for object 11627 the flux is \(2.3\pm0.0\) Jy km s\(^{-1}\).

# With Pandas.
print(df_lit_s.loc[df_lit_s["name"] == "UGC 11627"])

print("\n") # Blank line.

# Without Pandas.
idx = data["name"].index("UGC 11627")
print(f"Flux for UGC 11627: {data['flux'][idx]} +- {data['flux_error'][idx]} Jy km/s")

         name  flux  flux_error  width  width_error
94  UGC 11627   2.3         0.0  126.0          9.0


Flux for UGC 11627: 2.3 +- 0.0 Jy km/s

Now that we have the results, plot them.

plt.figure(dpi=150)
# Plot a 1-to-1 line using the observed values.
plt.plot(sorted(df_obs_s["flux"]), sorted(df_obs_s["flux"]), "k--")
plt.errorbar(df_obs_s["flux"], df_lit_s["flux"], 
             xerr=df_obs_s["flux_error"], 
             yerr=df_lit_s["flux_error"], 
             marker=".", ls="", ms=5)
for n,o,l in zip(df_obs_s["name"], df_obs_s["flux"], df_lit_s["flux"]):
    plt.text(o, l, n, ha="center", va="bottom")
plt.tight_layout()
plt.xlabel("Measured Flux (Jy km s$^{-1}$)")
plt.ylabel("Literature Flux (Jy km s$^{-1}$)");

There are significant differences between the tabulated values and those measured from this data set. At the time of writing we do not understand the origin of the difference. As far as we can tell, processing the data in GBTIDL produces equivalent results as those presented here (the differences in flux measurements are <1 Jy km s\(^{-1}\)). Stay tuned for an update.

plt.figure(dpi=150)
# Plot a 1-to-1 line using the observed values.
plt.plot(sorted(df_obs_s["width"]), sorted(df_obs_s["width"]), "k--")
plt.errorbar(df_obs_s["width"], df_lit_s["width"], 
             xerr=df_obs_s["width_error"], 
             yerr=df_lit_s["width_error"], 
             marker=".", ls="", ms=5)
for n,o,l in zip(df_obs_s["name"], df_obs_s["width"], df_lit_s["width"]):
    plt.text(o, l, n, ha="center", va="bottom")
plt.tight_layout()
plt.xlabel("Measured Width (km s$^{-1}$)")
plt.ylabel("Literature Width (km s$^{-1}$)");

The widths agree much better, since they do not depend on the flux calibration.