 News |
FAQ
|
| Search
|
Home
|
| Getting Started | Documentation | Glish | Learn More | Programming | Contact Us |
|
| Version 1.6 Build 363 |
|
T.J. Cornwell
This chapter describes how to make single dish images in AIPS++. It should be read in conjunction with several other chapters in Getting Results, including the chapters which cover: i) use of the viewer, and; ii) image analysis.
The examples here are worked from GBT data but the principles apply to other telescopes. Single dish calibration and imaging is still being improved and so some areas may change over time. In particular, the calibration approach used here will be replaced shortly.
GBT data in are read into AIPS++ using the gbtmsfiller tool. At the shell command line:
gbtmsfiller project=/home/gbtdata/pnt_prime_13 backend=DCR fillrawpointing=True minscan=2360 maxscan=2415
This will create a MeasurementSet call pnt_prime_13_DCR in the current directory.
Note that at the time of writing, the data filled by the filler has some omissions. To fix these, use the following script:
include 'gbtcalutils.g'
gtbcalutils.fixpnt('pnt_prime_13_DCR');
Data in SDFITS format may be read in AIPS++ using the general SDFITS filler ms, which forms part of the ms tool. This is achieved as:
- m:=sdfitstoms(fitsfile='/aips++/data/demo/CALS00FEB.CUV',
msfile='cals00feb.ms');
# Create an AIPS++ uv-data file from a UVFITS file
- m.done();
Interactive editing is possible on almost all graphical views of the data which are possible in msplot. To enable the editing option, it is necessary to set parameter edit=T when the msplot tool is created.
include 'msplot.g'
msp:=msplot('pnt_prime_13_DCR', edit=T)
The GUI panel controlling flagging in msplot allows global flag masks to be set, such as that for all polarizations if any one is bad, as well as access to the list of accumulated flags, reversal of flags already applied, and other related utility operations. Raster displays of data are possible also, and are recommended for spectral line data editing.
You may also view and edit the data at a lower level: that of the tables themselves. To browse the MeasurementSet table using a graphical browser, do:
include 'table.g'
t:=table('pnt_prime_13_DCR')
t.browse();
Note that editing of individual values using the browser is possible (but unlikely to be needed unless there is something structurally wrong with the tables).
As noted above, AIPS++ data are stored in a MeasurementSet (MS). The observed data are stored in the DATA (complex) or FLOAT_DATA (real) column in the MAIN table of the MS; these are the raw data as loaded from the VLA archive or from a UVFITS file. Associated with the DATA column are related columns to hold the most recent version of the calibrated data (CORRECTED_DATA), and the latest data plane representation of the source or field model (MODEL_DATA). The latter two columns are filled in by the calibrater and imager tools. The observed DATA do not change during reduction, but the related columns do, and should be viewed in this context. When plotting the data using a tool such as msplot, the data sub-types (observed, corrected and model) can be selected individually.
Calibration of single dish data is still under development in
AIPS++
he plan is that calibration will be unified for single
dish and synthesis data. Meanwhile, the powers of glish and the table
system can be used to calibrate the data. For example, the following
script provides a simple means of performing on-off calibration.
sdcalibrate := function(msname, width=100) {
#
# Use imager to add the optional columns
#
include 'imager.g';
im:=imager(msname);im.done();
#
# Open the table for read/write
#
tab := table(msname, readonly=F);
#
# Get data for all cal phases
#
phases := [=];
phases[1] := tab.query('NRAO_GBT_RECEIVER_ID==0&&NRAO_GBT_STATE_ID==0');
phases[2] := tab.query('NRAO_GBT_RECEIVER_ID==0&&NRAO_GBT_STATE_ID==1');
phases[3] := tab.query('NRAO_GBT_RECEIVER_ID==1&&NRAO_GBT_STATE_ID==0');
phases[4] := tab.query('NRAO_GBT_RECEIVER_ID==1&&NRAO_GBT_STATE_ID==1');
times := [=];
data := [=];
for (i in 1:4) {
times[i] := phases[i].getcol('TIME');
times[i] -:= min(times[i]);
data[i] := phases[i].getcol('FLOAT_DATA');
}
#
# This seems to be the one that we want
#
cal := data[4] - data[3];
#
# Smooth
#
scal := smooth(cal, width);
#
# Plot the smoothed cal
#
include 'pgplotter.g';
p:=pgplotter(background='white', foreground='black');
p.plotxy(times[1], scal);
p.lab('Time (s)', 'Smoothed cal', 'Smoothed calibration signal vs time');
p.postscript(spaste(msname, '.cal.ps'));
#
# Recalibrate
#
data[2] := data[2]/scal;
#
p.clear();
p.plotxy(times[1], data[2]);
p.lab('Time (s)', 'Flux', 'Final calibrated flux vs time');
p.postscript(spaste(msname, '.final.ps'));
#
# Put back into the corrected data for second phase
#
note('Overwriting phase [0,1] in the CORRECTED_DATA column', priority='WARN');
print phases[2].putcol('CORRECTED_DATA', complex(data[2]));
#
tab.flush();
tab.done();
note(msname,' has been calibrated');
}
where smooth is a running-median smoothing function:
#
# Smooth a function using a running median filter
# (expensive but robust)
#
smooth := function(x, width=5) {
include 'statistics.g';
newx := array(0.0, length(x));
for (i in 1:length(x)) {
newx[i] := median(x[max(1,(i-width)):min(length(x), (i+width))]);
if(i%100==1) {
print i, 'of ', length(x), ':', x[i], '->', newx[i];
}
}
return newx;
}
Imaging is performed using the imager tool, which supports a broad range of imaging modes and deconvolution algorithms, as described in the imaging chapter in Getting Results and in the User Reference Manual.
The imager tool, which can be used to make images for each field or source, requires initial selection of the data and initialization of the imaging coordinates. Imaging, deconvolution and image restoration operations are then possible, as described in more detail in the imager documentation.
A single-field, continuum GBT unnormalized image can be made as follows:
include 'imager.g';
#
myimager:=imager('pnt_prime_13_DCR');
#
# Use the data from one cal phase only
#
myimager.setdata(fieldid=1,spwid=1,
msselect='NRAO_GBT_RECEIVER_ID==0&&NRAO_GBT_STATE_ID==1');
#
# Set the direction for the center of the image
#
dir:=dm.direction('J2000', '20:20:00.00', '+040.00.00');
#
myimager.setimage(nx=400, ny=400,
cellx='2arcmin',celly='2arcmin',
stokes='I',
doshift=T, phasecenter=dir, spwid=1)
myimager.setoptions(gridfunction='PB')
#
# Fill in the weights
#
myimager.weight('natural');
#
# Make the image
#
myimager.makeimage(image='scanimage', type='singledish');
myimager.close();
This image is not normalized by the sampling density. The following additional steps rectify that:
# Make the coverage image
#
myimager.makeimage(image='scanweight', type='coverage');
#
# Display the coverage image
#
include 'image.g';
imcov := image('scanweight');
s:=0;
imcov.statistics(s);
print s;
threshold := s.max / 10.0;
print "Thresholding at ", threshold;
#
# Normalize the image, threshold coverage image to avoid nonsampled points
#
im:=imagecalc('pnt_prime_13.bigimage',
pixels=spaste('scanimage[scanweight>', threshold,
']/scanweight[scanweight>', threshold, ']'))
#
# View the image
im.view(raster=T, axislabels=T);
The coverage image is simply the gridded weights. Note the use of the image calculator imagecalc to normalize one image by another.
At this point, you may have many files in the local directory. The catalog may be used to browse these files:
include 'catalog.g'; dc.gui()
The process of making a single dish image is very straightforward. The data are converted into the coordinate system of the image (set via imager.setimage), and then added to the image using a gridding function. A variety of convolution functions can be used: BOX, SF, or PB. The first is a simple nearest neighbor gridding, the second uses a prolate spheroidal wavefunction, and the last uses a primary beam model appropriate to the telescope being used. The primary beam is optimal in the least squares sense but degrades the resolution, whereas the prolate spheroidal function can avoid the resolution degradation but the noise level will be higher.
The imager documentation gives many more details of the options available in imaging.
The imaging process for a single dish telescope may be simulated using the simulator. It is easiest to do this from an existing MeasurementSet. Make a copy and then do the following on the copy:
include 'simulator.g';
#
mysimulator:=simulatorfromms('m31sim.ms');
#
# Use the SD gridding machine with the primary beam
# as the convolution function
#
mysimulator.setoptions(gridfunction='PB', ftmachine='SD');
#
# Predict the model fluxes
#
mysimulator.predict(modelimage='gbtm31.model');
#
# Close
#
mysimulator.close();
Noise and calibration effects may be added. The suite of available error formats will be expanded over time.
There are some nice images to use in the data repository (usually at e.g. /aips++/data). For example, to make the image used above, do the following:
include 'image.g'
#
# Open a pre-existing image with the correct coordinates
#
im:=image('scanweight');
cs:=im.coordsys();
#
# Make an image from the FITS file of the model
#
model:=imagefromfits('gbtm31temp.model', '/aips++/data/demo/M31.model.fits');
#
# Now construct a new image and insert the pixels from the model
#
imm31:=imagefromshape('gbtm31.model', shape=[256,256,1,1], csys=cs);
imm31.putchunk(model.getchunk());
imm31.view();
At this point you should have an interesting image to look at. To understand some of the steps involved in generating this image, you may wish to review the top level documentation, particularly Getting Started in AIPS++.