Frequencly Asked Questions

What is the netCDF format.

netCDF is a portable binary dataformat, a close cousin of HDF, originally designed by NOOA? In fact, the storage model in netCDF-4 is now HDF5, so you can now use the h5dump program as well as the native (and more readable) ncdump program to look at an ascii representation of the data.

A simple way to get at the raw LMT data is the following snippet of python code:

import netCDF4
nc = netCDF4.Dataset(filename)
rawdata = nc.variables['Data.Integrate.Data'][:]
nc.close()

where rawdata is now a 2D array, shaped (ntime,nchan) in the python sense (nchan runs fastest). The LMT software will typically calibrate and convert these data to the more common FITS format.

What is SDFITS format.

SDFITS is a legal FITS file, but it has a BINTABLE (a spreadsheet) extension. This BINTABLE stores data in rows and columns. Typically each row is a spectrum with meta-data (e.g. ra,dec,beam,pol,band,etc.) and the number of rows are the different observations (different ra,dec,beam,pol,band, etc.).

What is a typical workflow for OTF data reduction

First find out which obsnums belong to a given observation, where a common sky and spectral are. Perhaps the lmtinfo.py can be useful too:

lmtinfo.py data

Lets say you find there are only three: 85776 85778 85824. You start with one:

SLpipeline.sh obsnum=85776

this will attempt to make a cube of the line specified by the VLSR in the header, and make some decisions on where to place the baseline fitting. Inspect the waterfall cube (the .wt.fits file), or the plots from view_spec_file. Witht this edit the settings in lmtoy_85776.rc. Simply edit and re-run:

lmtoy_reduce obsnum=85776

Once satisfied with the fits cube, copy this for the other obsum, and reduce them

for o in 85778 85824; do
    cp lmtoy_85776.rc lmtoy_${o}.rc
    lmtoy_reduce.sh obsnum=${o}
done

Inspection is needed again, pixels may have to be removed etc. Once satisfied, combine them:

lmtoy_combine.sh obsnum=85776,85778,85824

How to make a beam as used by the gridder

There is no option in grid_data.py yet, but the -a flag in spec_driver_fits will make a beam. Pick any SpecFile, select only 1 pixel (the -u flag), although this is not essential, and inspect the -w beam.fits as the beam. The -o file can be discarded. Example:

spec_driver_fits -i IRC_79448.nc -o beam.1.fits -w beam.fits \
   -u 0 -z 4 -s 1 -x 4 -y 4 -f 1 -r 1 -n 256 -0 1.1 -1 1 -2 2 -b -1 -a -l 2 -c 1

CAVEAT: this option is still under development.

Cubes in RA-DEC or GLON-GLAT

The default mode of OTF observing is a scan in Alt-Az with with a known parallactic angle gets converted into RA-DEC space. You will find the WCS of the FITS to be in RA---SFL and DEC--SFL.

In MIRIAD the task cotra will transform coordinate values between coordinate systems, whereas the regrid task can, using a template, convert between equatorial and galactic coordinates.

In CASA use the imregrid task. For example converting from equatorial to galactic

imregrid(imagename="input.im", output="output.im", template="GALACTIC")

though if the input image has the SFL projection, CASA cannot reliably handle this. The tool cs.setprojection() can be used temporarely to allow rotation of the image, viz.:

ia.open('input.im')
csys = ia.coordsys()
csys.setprojection('TAN')
im2=ia.regrid('input2.im',csys=csys.torecord(),overwrite=True)
im2.close()
ia.close()
imregrid('input2.im','GALACTIC','output.im',overwrite=True)

A better approach is to obtain the correct Glon-Glat value for each sample, and let the regridder convolve them correctly using the SFL projection method. This option is being added to spec_driver_fits yet.

Conversion of intensity units from K to Jy/beam

By default the intensity units at the LMT are Kelvin, but another common unit system is Jy/beam.

Cubes in VLSR or FREQ

Both MIRIAD and CASA have conversion routines between them, but they will depend on other keywords to be correct. WCS can be a tricky thing, especially if you need good accuracy. Currently LMT cubes have the following keywords that influence the WCS:

BUNIT   = 'K       '           /
CTYPE1  = 'RA---SFL'           /
CTYPE2  = 'DEC--SFL'           /
CUNIT2  = 'deg     '           /
CTYPE3  = 'VELO-LSR'           /
CUNIT3  = 'm/s     '           /
EQUINOX =                2000. /
RADESYS = 'FK5     '           /
RESTFRQ =        115271204000. / Header.LineData.LineRestFrequency
SPECSYS = 'LSRK    '           / could be wrong (check ? Header.Source.VelSys)

In MIRIAD:

MIRIAD supports the VOBS concept (Velocity of the observatory w.r.t. restframe), but although this value is in the RAW header, it is not passed on to the SpecFile

In the RAW headers we can find (e.g.)

Header.Source.Velocity = 463 ;
Header.Source.VelSys = 0 ;
Header.Sky.ObsVel = -19.9564373466461 ;
Header.Sky.BaryVel = -11.6443975487706 ;

shell examples

rm -rf irc.mir
fits in=IRC_79448.fits out=irc.mir op=xyin
imhead in=irc.mir
velsw in=irc.mir axis=freq
fits in=irc.mir out=irc.mir.fits op=xyout

Importing this FREQ based file, makes it work in CASA.

importfits('irc.mir.fits','irc.mir.fits.im',overwrite=True)
imhead ('irc.mir.fits.im')
exportfits('irc.mir.fits.im','irc.mir.fits.im.fits',velocity=True,overwrite=True)

but the following one does not work

exportfits('irc.mir.fits.im','irc.mir.fits.im.fits',velocity=True,overwrite=True,optical=True)

Miriad also differentiates between CELLSCAL=’CONSTANT’ and CELLSCAL=’1/F’

in CASA:

Our LMT fits file is linear in frequency and velocity (VELO-LSR, not FELO-LSR) yes, exportfits complains about non-linear axis unless we say optical=True. issue?

importfits('IRC_79448.fits','irc.im',overwrite=True)
imhead('irc.im')
exportfits('irc.im','irc.fits',velocity=True,optical=True,overwrite=True)

The reference pixel is 350.303, but in exportfits I see ALTRVAL=349.894 It works fine of CTYPE3 is VRAD, and not the current VELO-LSR, but VELO_LSR is not a recognized axis name, so it sticks to the (correct) m/s, but doesn’t know about FREQ.

Baseline subtraction

Normally spectral data will need to be baseline subtracted before conversion to a SpecFile or SDFITS file. In particular for an SDFITS file this is not strictly needed, if you take the data to another package that can also do this more flexibly.

Collecting pipeline results

For variables that have been placed in the lmtoy_OBSNUM.rc file, there is a fast way to collect a table summarizing the variables using the nemopars program. Here is an example:

cd $WORK_LMT/2021-S1-MX-3

nemopars  src,obsnum,skytime,inttime,tint_on,nrows */lmtoy_*.rc | sort  > log1
awk '{print $3/($5+$6*1.6)}' log1 | tabhist - 1 1 3 bins=20 resid=f qac=t


nemopars  skytime */lmtoy_*.rc | tabhist -

In the last example the total skytime will be histogrammed, and the screen output should contain the line with the Sum (in seconds)