dysh for GBTIDL Users#
Some key differences#
No longer tied to a single input and output file, can load as many SDFITS files as RAM allows
No longer tied to a finite number of data containers, you can have as many ScanBlock/Spectrum objects as RAM allows
Python allows for well-trackable local and global variables
Rough equivalents#
In the tables below, it is assumed you have executed the following commands in Python:
from dysh.fits import GBTFITSLoad, GBTOnline, GBTOffline
If you are running from the dysh shell, then the above commands are already available.
File I/O and Metadata Operations#
GBTIDL
dysh
Notes
online
sdf = GBTOnline(); sdf = GBTOnline([projID])
Monitor the currently online project or a specific project
projID. Load the SDFITS intosdf.offline
sdf = GBTOffline([projID])
Load data from project/session
projIDintosdf.filein, [SDFITS file]
sdfits = GBTFITSLoad( [SDFITSfile] )
This will directly load a SDFITS file into python.
dirin, [SDFITS dir]
sdfits = GBTFITSLoad( [SDFITS dir] )
This will load each SDFITS file in a directory into python, in the case it is a project directory where each file is a different VEGAS bank. (A, B, C, etc.). The implementation is identical to the single SDFITS file load.
fileout
There is no direct translation; in dysh, simply specify a filename with the
Spectrum.writeorScanBlock.writefunctions.keep
spec.write( [SDFITSfile], format=[fmt] )
Write the spectrum
specto the outputSDFITSfilespecified, using the formatfmt.summary
sdfits.summary()
In dysh, you can also input a list of scan numbers you would like to see, or modify which SDFITS columns are printed.
list
sdfits.summary( verbose=True )
Using the
verbose=Trueoption returns similar results to GBTIDL’slist. It is possible to provide a scan number asscan=x.header, scan_info
spec.meta
Return meta information about the spectrum
spec. In GBTIDL, this returns a text printout. In dysh, this returns a dictionary with all SDFITS keywords and values.
Calibration and Data Retrieval#
GBTIDL
dysh
explanation
select
sdf.select( [col1]=[val1], [col2]=[val2] )
Select subsets of the
GBTFITSLoadobjectsdfbased on the constraints oncol1andcol2given byval1andval2flag
sdf.flag( [col1]=[val1], [col2]=[val2] )
Create flagging rules of the
GBTFITSLoadobjectsdfbased on the constraints oncol1andcol2given byval1andval2. See flagging examples given in the Recipes.gettp, [scan_num]
sdfits.gettp( scan=[scan_num] )
Get a total power scan. Returns a
ScanBlock.getps, [scan_num]
sdfits.getps( scan=[scan_num] )
Retrieve and calibrate position-switched data. Returns a
ScanBlock.getfs, [scan_num]
sdfits.getfs( scan=[scan_num] )
Retrieve and calibrate frequency-switched data. Returns a
ScanBlock.getnod, [scan_num]
sdfits.getnod( scan=[scan_num] )
Retrieve and calibrate nodding data. Returns a
ScanBlock.show
spec.plot()
Plot the spectrum
spec. Calibration routines in GBTIDL do this automatically.
Spectrum Operations#
GBTIDL
dysh
Notes
baseline, bshape, bsubtract, nfit
spec.baseline( [deg], exclude=[ex_reg],remove=True/False )
Compute and optionally remove a baseline with degree
degand exclusion regionsex_reg.gsmooth, [int]
spec.smooth( “gauss”, [int] )
Convolve spec with a Gaussian kernel. Default behavior is to decimate, which isn’t true in GBTIDL.
hanning
spec.smooth( “hann” )
Convolve
specwith a Hanning window. Default behavior is to decimate, which isn’t true in GBTIDL.boxcar, [int]
spec.smooth( “box”, [int] )
Convolve
specwith a boxcar window. This matches the GBTIDL results only when integerintis an odd number. Default behavior is to decimate, which isn’t true in GBTIDLadd/subtract/divide/scale
spec + other, spec - other, spec/other, spec*other
Perform element-wise math between spectrum
specand other vector/scalarother.bias, [bias]
spec + bias
Add scalar
biasto spectrumspec.accum
average_spectra(list of spectra)
average_spectrahas to be imported first (from dysh.spectra import average_spectrum). Alternatively,spec.average(list of other spectra)will averagespecwith the spectra in the list.ave
scanblock.timeaverage(); scan.timeaverage()
Average a set of scans inside
ScanBlockscanblockorScanBasescan.stats
spec.stats()
Print statistics of spectrum
spec.fshift, vshift, xshift
spec.find_shift(other)
Find the shift between
specandother.gshift
spec.shift(shift); spec.align_to(other)
Shift
specbyshift;Find the shift betweenspecandotherand apply. Equivalent tospec.shift(spec.find_shift(other)).write_ascii, filename
spec.write(filename, format=”ascii.basic”)
Write the contents of
specto a text file in ascii format with minimum header information.
Side-by-side Examples#
The following are examples in GBTIDL and dysh that produce equivalent results.
The dysh examples assume they are being run from the dysh shell, so that the modules imported on startup are available (e.g., from astropy import units as u).
OTF Mapping#
This example is based on the observations under TGBT17A_506_11, the same observations used for the on-the-fly data reduction section of the users guide.
sdfits = GBTFITSLoad("TGBT17A_506_11.raw.vegas")
scan_block = sdfits.getsigref(scan=list(range(14,27)), ref=27, fdnum=0, ifnum=0, plnum=0)
scan_block.write("dysh.fits")
filein,"TGBT17A_506_11.raw.vegas"
fileout,"gbtidl.fits"
for s=14,26,1 do begin &$
getsigref,s,27,/avgref,/keepints &$
keep &$
endfor
Extragalactic 21 cm Observations using Position Switching#
This example was borrowed from GBTdocs position switched tutorial. A more thorough dysh example using the same data can be found in example reduction of an HI survey. The GBTIDL version of this example can only be run with a display.
sdfits = GBTFITSLoad("/home/astro-util/HIsurvey/Session02")
tp0 = sdfits.gettp(scan=299, ifnum=0, plnum=0, fdnum=0).timeaverage()
tsys0 = tp0.meta["TSYS"]
tp1 = sdfits.gettp(scan=299, ifnum=0, plnum=1, fdnum=0).timeaverage()
tsys1 = tp1.meta["TSYS"]
sigref0a = sdfits.getsigref(scan=[296], ref=295, ifnum=0, fdnum=0, plnum=0, t_sys=tsys0, units="flux", zenith_opacity=0.08).timeaverage()
sigref0b = sdfits.getsigref(scan=[298], ref=297, ifnum=0, fdnum=0, plnum=0, t_sys=tsys0, units="flux", zenith_opacity=0.08).timeaverage()
sigref0 = sigref0a.average(sigref0b)
sigref1a = sdfits.getsigref(scan=[296], ref=295, ifnum=0, fdnum=0, plnum=1, t_sys=tsys1, units="flux", zenith_opacity=0.08).timeaverage()
sigref1b = sdfits.getsigref(scan=[298], ref=297, ifnum=0, fdnum=0, plnum=1, t_sys=tsys1, units="flux", zenith_opacity=0.08).timeaverage()
sigref1 = sigref1a.average(sigref1b)
sigref = sigref0.average(sigref1)
sigref_smo = sigref.smooth(kernel="gauss", width=100, decimate=0)
region = [[1.402*u.GHz, 1.4045*u.GHz],
[1.40506*u.GHz, 1.4054*u.GHz],
[1.4072*u.GHz, 1.4115*u.GHz]]
sigref_smo.baseline(1, model="poly", include=region, remove=True)
kms = u.km/u.s
stats_b = sigref_smo[2000*kms:2500*kms].stats()
stats_r = sigref_smo[3500*kms:4000*kms].stats()
rms = (stats_b["rms"] + stats_r["rms"])/2.
cog = sigref_smo.cog(bchan=100, echan=200)
dirin,'/home/astro-util/HIsurvey/Session02'
freeze ; keep from plotting on the screen, to speed up the processing
gettp, 299, plnum=0, /quiet
tsys0 = !g.s.tsys
gettp, 299, plnum=1, /quiet
tsys1 = !g.s.tsys
for i=295,297,2 do begin &$
for p=0,1 do begin &$
if (p eq 0) then tsys=tsys0[0] else tsys=tsys1[0] &$
getsigref, i+1, i, plnum=p, tsys=tsys, unit='Jy', /quiet &$
accum &$
endfor &$
endfor
ave
gsmooth, 100, /decimate
setxunit, 'GHz'
setx, 1.401, 1.412
show
region = [1.402, 1.4045, 1.40506, 1.4054, 1.4072, 1.4115]
reg1 = xtochan(region)
Nregion, reg1
nfit, 1
bshape
baseline
velo
show
stats, 2000, 2500, ret=mystats
rms1 = mystats.rms
stats, 3500, 4000, ret=mystats
rms2 = mystats.rms
rms = (rms1 + rms2) / 2
gmeasure, 1, 0.5, brange=2815, erange=3142, rms=rms
gmeasure, 1, 0.2, brange=2815, erange=3142, rms=rms