|
|
|
____________________
|
|
|
|
Advanced Methods
Advanced
Methods
|
Changing coordinate attributes with wcsattrib
There may be times when you will want to change the coordinate systems
displayed by your ACSIS data. For instance, changing between frequency
and velocity coordinates for spectra, or to offsets (given in
arcseconds) schemes. This is best done with the wcsattrib
command in KAPPA.
To change
to Offset coordinates:
> wcsattrib <filename> set skyrefis origin
> wcsattrib file set 'format(1)' 's.*'
> wcsattrib file set 'format(2)' 's.*'
The last two lines are necessary to have the offset units in
arcseconds, otherwise offsets will be in radians. To go back to your
original settings:
> wcsattrib <filename> set skyrefis ignored
> wcsattrib file clear 'format(1)'
> wcsattrib file clear 'format(2)'
To change to frequency coordinates:
> wcsattrib <filename> set 'system(3)' freq
Here, system(3) refers to the third axis which is now set to frequency space. To go back to radio velocity:
> wcsattrib <filename> set 'system(3)' vrad
To change to Galactic coordinates:
> wcsattrib <filename> set 'system(1)' galactic
Note that this does not resample the image. However, all applications which are now used on <filename> will now report positions in Galactic coordinates rather than J2000. To change back you need:
> wcsattrib <filename> set 'system(1)' 'equatorial(j2000.0)'
To change to a Geocentric rest frame:
By default, the rest frame of the spectral axis is displayed in the topocentric rest frame (at the Earth's surface at the time the observations were made) so that the velocity from the Earth's rotation is taken into account. There may be instances, such as when observing planets, when you need to use a Geocentric (centre of the Earth) rest frame. To switch to a Geocentric rest frame use:
> wcsattrib <filename> set stdofrest geoc
|
Despiking
data cubes
Interference spikes are sometimes found in ACSIS/HARP data. If
you suspect you have a spike in your data cube the following
procedure describes a method which allows you to identify the position
of the spike and to mask it out. Spikes should be quite conspicuous in
the data cube as they occur at the same frequency throughout the
observed cube and at (almost) uniform intensity.
It is probably best to eliminate spikes in your data cube from the time
series, i.e. before running makecube
or any other procedures on your data. The first thing to do is to
collapse your time series cube into a single spectrum. You'll need to
do this twice, once for each axis. To accentuate the presence of a
spike over data the collapse should use the rms estimator,
> collapse in=a20070113_00022_00_0001 estimator=rms axis=3 out=temp1
> collapse in=temp1 estimator=rms axis=2 out=spike
The time series cubes are constructed such that the spectrum is in axis=1
so we need to collapse the other axes to get an averaged spectrum. The
figure below shows the result of the two collapse
commands above but with estimator=mean
to produce a coadd of the whole cube.
The figure below shows the effect of using estimator=rms
on collapsing the cubes. Spikes (2 are seen in this cube) are accentuated
over the real signal
in the data (note the logarithmic scale). The science signal is seen at
around 0 km/s.
For some regions where the line strength might be quite
bright in places, collapsing using the mean may not produce a
spike signal which can be differentiated from data (although the real
line would have to be bright and narrow to trouble the eye). In this
example the spike has a signal of approximately +/- 3-4 K whereas the
science signal is, on average, of the order of 1.5 K. Interestingly,
the mean spectrum seems to miss the presence of a possible spike at
-175 km/s. As the interference
spike is
localised to particular and narrow range of frequencies (and
may
be
present
across the whole cube at more or less uniform intensity), it
will
always have a large rms deviation from the mean measured across the
whole cube,
as demonstrated above.
The next step involves locating the background level lying underneath
the spikes. The challenge here is to identify the spikes but preserve
the real science signal. There is a command within the CUPID
package, called findback,
which does that. The trick to using findback is to
set a size scale, with the parameter box,
which is going to be smaller than any real features but larger than a
spike. We find that values of order 5 or 7 suffice but it is best that
you experiment with this parameter with your own data. In the
following example a value of 5 is used (remember to type CUPID
into the command line first if you haven't done so already),
> thresh in=spike out=spike_thr thrlo=0 newlo=bad thrhi=1e4 newhi=bad
> findback in=spike_thr box=5 out=spike_back
The first command is there to mask out any large values that
may have been generated at the edges of the spectrum (where
the
noise blows the rms to very high values) by the collapse
routine. The upper threshold value is set to 10,000 (check that your
spike doesn't peak above this value).
The figure below shows a close up of the region where
the spikes are found in the time series. The original spike spectrum
from above is in red, with the background calculated from findback
overlayed in blue.
The red spikes remain
prominent
above the blue, interpolated, background while the science signal at 0
km/s is fitted to as part of the background. This background can now be
removed from the original leaving behind just the spikes.
> sub in1=spike_thr in2=spike_back out=spike_only
This produces a relatively flat spectrum whose only deviations are the
spikes. To make the mask we use the KAPPA
ffclean
command which attempts to fit a flat background from which any large
deviations are set to be BAD,
> ffclean in=spike_only out=spike_gone box=5 clip='[3,5]'
The clip
parameter determines the sigma clipping used to decide which values are
set to BAD.
Two values means that two passes are made. The first 3-sigma clip will
get rid of the largest deviations, while a more stringent second pass
at 5-sigma should remove any deviant pixels which remain.
A before/after comparison is made in the close up image above. The
animated gif shows the spike_only
spectrum, i.e. before ffclean,
and the blue spectrum overlaid is the same after ffclean. The
spikes have been set to BAD
(they are still red) and the region where the science signal is remains
preserved (blue).
The next step is to grow the 1D mask up into a 3D mask. This can be
done using the wcsalign
command, using the initial time series cube as a reference for the
dimensionality of the resultant cube.
> wcsalign in=spike_gone out=spike_gone3d ref=a20070113_00022_00_0001
method=near accept
The final step is to copy the locators of BAD pixels to
the original time series cube,
> copybad in=a20070113_00022_00_0001 ref=spike_gone3d out=a20070113_00022_00_0001_clean
Now run this cube through the makecube
command as described
in the main part of the document and proceed as normal. Your new cube
should be free from spikes.
|
Converting data to CLASS format
If you prefer to reduce your data cubes using the CLASS software then you need to convert the sdf files into appropriate fits files. At the time of writing, CLASS is not able to read in the standard fits files that are generated by JAC instrumentation. However, the CONVERT routines are able to produce the correct encodings. Conversion is done with the ndf2fits command.
However, given that single spectra need to be read into CLASS you need to extract each spectrum from the cube. This can be done with the KAPPA command ndfcopy. To extract the (i,j)-th spectrum from a cube use the following commands:
> ndfcopy 'cube(i,j,)' out=spec_ij
> ndf2fits spec_copy encoding=fits-class out=spec_ij_class.fits
Finally, you must open CLASS and convert each of these fits files to a CLASS spectrum which can then be saved into a CLASS data file using the commands:
> file spec_ij_class class.dat new
> set angle sec
> LAS\FITS READ spec_ij_class.fits
> WRITE 1
Clearly, with large data cubes this can become quite tiresome so a script is the most economic method of performing this conversion. A separate guide on how to do this has been provided by Jan Wouterloot in this web page.
|
Combining DAS and ACSIS data
With several years of data from the DAS available in the archive (and in your drawers...?) it is conceivable that at some point it will become necessary to combine DAS data with the current format delivered by ACSIS. To enable this, the DAS data needs to be converted to NDF. There is a straightforward method of doing this which requires reading the DAS data into SPECX and then converting it using the CONVERT routine specx2ndf.
To combine DAS and ACSIS data:
- Read the DAS spectrum into SPECX;
- Write the spectrum out to a specx file or specx map file (exporting to a FITS cube will give erroneous results);
- Use specx2ndf on the specx file just created;
You now have a NDF file on which you can run all the starlink routines, including wcsmosaic in KAPPA which will allow you to combine several NDFs together.
NDF files from a normal specx file (filename.sdf) can be displayed with GAIA or SPLAT by typing
> splat filename.SPECTRUM1
or
> gaia filename.SPECTRUM1
(and SPECTRUM2 etc if the file contains more spectra).
If filename.sdf is a specx map file (also for a single position), the
resulting NDF cube can be displayed in the normal
way with GAIA and SPLAT.
By the middle of 2008 it is planned to have a SMURF command available which
can be used to directly read DAS spectra, bypassing
SPECX.
|
How to identify bad receptors
There are two other ways to identify receptors belonging to a pixel in a map (and
subsequently mask and remove the bad receptors from the map).
When making the map with the SMURF command makecube, use the parameter
outcat as in the example below (a grid position switch map):
> makecube in=a20071123_00045_01_0001.sdf out=a20071123_00045.sdf outcat=a20071123_00045.FIT autogrid
The file a20071123_00045.FIT is a catalogue file which can be used to plot receptor names on the map using the KAPPA commands
display and listshow:
> display 'a20071123_00045(,,1)'
> listshow a20071123_00045.FIT label
The other way to identify bad receptors is to load the map in
GAIA,
> gaia a20071123_00045.sdf
and then click on Data-Servers --> Local Catalogs and load the catalogue a20071123_00045.FIT
This will plot symbols (default large circles) at the positions of the
receptors and open a new window with search results (see below).
In the latter window one can change the size of the plotted symbol by clicking
on Options --> Set Plot Symbols (e.g. from 4 to 0.2), which allows clicking on
the pixels to show the spectrum. One can also change the symbol from circle to
cross. Clicking on a cross will show the receptor in the window with search
results. This window shows all 900 individual measurements made to obtain this
map and coadded by makecube
|
Hybrid (440 and 1860 MHz) data
Observations with a bandwidth of 440 or 1860 MHz consist of two ACSIS subbands
(see the ACSIS guide).
Until August 2006 these subbands were merged by the internal ACSIS software. However this was not done correctly (spectra often showed unrealistic features or high noise in the overlap region), so this automatic merging was turned
off.
Therefore the subbands should now be merged by using the
KAPPA command wcsmosaic, after truncating the noisy edges of the spectra either
by using the command ndfcopy or
by specifying a section (velocity range) of the spectrum in
wcsmosaic. The user has to decide
whether to subtract a baseline from the individual subbands before merging, or not, depending on e.g. the width of the expected spectral lines.
|
ORAC-DR
ORAC-DR is an automated astronomical data reduction pipeline which can also be used to reduce ACSIS data.
Output from ORAC-DR
are cubes and various further processed data which
currently (December 2007) are not yet quite finalized:
Baseline-fitted cubes, truncated cubes (bad outer parts of spectra removed),
integrated intensity maps, channel maps, group files (coadded cubes of data
at the same position and frequency).
Also pointing observations are reduced by ORAC-DR.
ORAC-DR can be started with
> oracdr_acsis
> setenv ORAC_DATA_IN
> setenv ORAC_DATA_OUT
> oracdr options
Where options can be chosen from the help file:
> oracdr -h
|
CSH scripts
Sometimes it is easier to reduce data by using CSH scripts with
KAPPA
commands, rather than repeating the same
KAPPA
commands for every source or doing an
excessive amount of clicking in GAIA or
SPLAT.
The script should start with
#!/bin/csh -f
source /star/etc/login >& /dev/null
source /star/etc/cshrc >& /dev/null
kappa >& /dev/null
Some scripts are available at the JCMT and in Hilo (in /local/progs/bin).
For help type
> sp_help
Examples:
sp_get (copy acsis file from observation directory).
sp_axis (set an old-style spectral axis (pixel, freq, vradetc).
sp_trim (trim noisy channels off the edges).
sp_coadd (coadd spectra).
sp_merge (frequency-merge two cubes).
sp_plot (Plot spectrum of file 1, with optionally file 2 overlayed in red).
sp_bsl (remove a baseline).
sp_bin (bin over a number of channels).
sp_stats (combined stats over one or more ranges).
sp_col (collapse cube over velocity/frequency range).
|
Position-Velocity diagrams
Position-Velocity diagrams can be made using the
KAPPA
display
command.
> kappa
> display cube'(0,,)' fill
will show a plot in Declination direction direction through the centre of the
cube. With e.g. cube'(10,,)' or cube'(,-10,)' plots at an offset in Right
Ascencion or Declination are shown.
The "fill"
parameter causes the image to be expanded to fill the screen on both
axes - otherwise one often gets a very long thin image.
Note that the default LUT of display results in a greyscale plot. Enter a lut*
(e.g. lutcol)
command before display
to use other LUT's (see kaphelp).
If one wants to see a 2D image with velocity on one axis, and position
along some arbitrary line on the sky on the other axis, one can do e.g.:
> kappa
> rotate cube0 axes=\[1,2\] angle=44 out=cuber
> display cuber'(0,,)' fill
This rotates the cube spatially by 44 degrees (in this example) and then
displays an image
of a slice perpendicular to the rotated first pixel axis.
|
|
|
