Merging SDFITS Files

Merging SDFITS Files#

This notebook shows how to merge SDFITS files. There is no built-in function to merge SDFITS files (yet), so we will put the files we want to merge in the same directory, and then load them and write them to a single SDFITS file. For this recipe we will use a single SDFITS file to generate two files from subsets of its data. We will need to use getsigref instead of getps to pull out a spectrum.

The following dysh commands are (re)introduced (leaving out all the function arguments):

  filename = dysh_data()
  sdf = GBTFITSLoad()
  sdf.summary()
  sdf.getsigref()
  sdf.write()

Loading Modules#

We start by loading the modules we will use for this recipe.

# These modules are required for loading and writing data.
from dysh.fits import GBTFITSLoad
from dysh.log import init_logging
from astropy import units as u

# These modules are used for file I/O
from dysh.util.files import dysh_data
from pathlib import Path

Setup#

We start the dysh logging, so we get more information about what is happening. This is only needed if working on a notebook. If using the CLI through the dysh command, then logging is setup for you.

init_logging(2)

# also create a local "output" directory where temporary notebook files can be stored.
output_dir = Path.cwd() / "output"
output_dir.mkdir(exist_ok=True)

Data Retrieval#

Download the example SDFITS data, if necessary.

filename = dysh_data(test="getps")

23:06:57.940 I Resolving test=getps -> AGBT05B_047_01/AGBT05B_047_01.raw.acs/

Data Loading#

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

sdfits = GBTFITSLoad(filename)

23:06:58.027 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.

sdfits.summary()

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

Writing SDFITS Files#

To show how to merge SDFITS files, we will create two separate SDFITS files first. We use the GBTFITSLoad.write method to write the data. We will write scans 51 and 52 to two separate files.

Create an Output Directory#

The approach we will use to merge SDFITS files relies on having the files to be merged in the same directory. We create a new temporary directory ./output/merge_sdfits

tmp_dir = output_dir / "merge_sdfits"
tmp_dir.mkdir(exist_ok=True, parents=True) # Create the output directory if it does not exist.

Write the Intermediate SDFITS Files#

Now we write scans 51 and 52 to the directory we created.

sdfits.write(tmp_dir / "scan51.fits", scan=51, overwrite=True)
sdfits.write(tmp_dir / "scan52.fits", scan=52, overwrite=True)

23:06:58.094 I write: bintable 0: 44 rows, 32768 channels, 1 chunk(s) of 5000

23:06:58.126 I write: chunk [0:44] of 44 rows built (0.0s)

23:06:58.139 I write: chunk written to disk (0.0s)

23:06:58.170 I write: bintable 0: 44 rows, 32768 channels, 1 chunk(s) of 5000

23:06:58.203 I write: chunk [0:44] of 44 rows built (0.0s)

23:06:58.215 I write: chunk written to disk (0.0s)

 ID    TAG    SCAN # SELECTED
--- --------- ---- ----------
  0 26291873b   51         44
 ID    TAG    SCAN # SELECTED
--- --------- ---- ----------
  0 d381daad3   52         44

Merge SDFITS Files#

Now that we have two SDFITS files we want to merge, we can do this by loading them using GBTFITSLoad by passing the directory containing the files we want to merge. In this case, that would be ./output/merge_sdfits.

sdfits_merged = GBTFITSLoad(tmp_dir)
sdfits_merged.files

23:06:58.303 I Loaded 2 FITS files

[PosixPath('/home/docs/checkouts/readthedocs.org/user_builds/dysh/checkouts/release-1.0.0/docs/source/users_guide/output/merge_sdfits/scan51.fits'),
 PosixPath('/home/docs/checkouts/readthedocs.org/user_builds/dysh/checkouts/release-1.0.0/docs/source/users_guide/output/merge_sdfits/scan52.fits')]

sdfits_merged.summary()

SCAN	OBJECT	VELOCITY	PROC	PROCSEQN	RESTFREQ	# IF	# POL	# INT	# FEED	AZIMUTH	ELEVATION
51	NGC5291	4386.0	OnOff	1	1.420405	1	2	11	1	198.3431	18.6427
52	NGC5291	4386.0	OnOff	2	1.420405	1	2	11	1	198.9306	18.7872

Write the Merged SDFITS Files to a New File#

Now we can write the merged SDFITS files to a new file using GBTFITSLoad.write We tell it to write to a single file with the multifile=False argument. We write to a new directory to show that this is only one file now.

new_dir = output_dir / "merged_sdfits"
new_dir.mkdir(exist_ok=True)

sdfits_merged.write(new_dir / "merged_sdfits.fits", multifile=False, overwrite=True)

23:06:58.341 I write: bintable 0: 44 rows, 32768 channels, 1 chunk(s) of 5000

23:06:58.365 I write: chunk [0:44] of 44 rows built (0.0s)

23:06:58.377 I write: chunk written to disk (0.0s)

23:06:58.380 I write: bintable 0: 44 rows, 32768 channels, 1 chunk(s) of 5000

23:06:58.404 I write: chunk [0:44] of 44 rows built (0.0s)

23:06:58.418 I write: chunk written to disk (0.0s)

Load the merged SDFITS file.

sdfits_merged = GBTFITSLoad(new_dir)
sdfits_merged.summary()

SCAN	OBJECT	VELOCITY	PROC	PROCSEQN	RESTFREQ	# IF	# POL	# INT	# FEED	AZIMUTH	ELEVATION
51	NGC5291	4386.0	OnOff	1	1.420405	1	2	11	1	198.3431	18.6427
52	NGC5291	4386.0	OnOff	2	1.420405	1	2	11	1	198.9306	18.7872

You can check that it is indeed a single file by looking at the files attribute of the GBTFITSLoad object.

sdfits_merged.files

[PosixPath('/home/docs/checkouts/readthedocs.org/user_builds/dysh/checkouts/release-1.0.0/docs/source/users_guide/output/merged_sdfits/merged_sdfits.fits')]

The same principles can be applied to more complex examples, like selecting an IF for certain scans, and another for others, and then merging them all in a single file.

Final Stats#

We then use getsigref to obtain two spectra, which using the Spectrum.stats should give the same answer

A note on this, the getps procedure does not work here, because after the merge, the ON and OFF are in different BINTABLES, and getps does not support this mode.

sdfits.getsigref(51,52,ifnum=0,plnum=0,fdnum=0)[0].timeaverage().check_stats(0.09095205 * u.K)

23:06:59.907 I rms is OK

sdfits_merged.getsigref(51,52,ifnum=0,plnum=0,fdnum=0).timeaverage().check_stats(0.09095205 * u.K)

23:07:00.434 I rms is OK