Open in Colab: https://colab.research.google.com/github/casangi/casadocs/blob/v6.5.2/docs/notebooks/image_analysis.ipynb
Image Analysis¶
The task viewer is deprecated in lieu of imview and msview, which contain the same functionality. Please invoke the imview (msview) task for running the CASA Viewer to visualize images or image cubes (visibility data).
CASA Images¶
CASA images are stored as tables and can be accessed with CASA tasks and tools. Image metadata can be listed and edited with the imhead task. Further processing includes the computation of statistics including spectral indices and polarization properties, transformation onto different spatial coordinates, spatial resolutions, and spectral frames, and many other processes (see the following section on Dealing with Images for a description of tasks that operate on CASA images).
Image Headers
Image Headers contain metadata on the observation – e.g. the observing date, pointing position, object observed, etc., and the resulting image – e.g. the restoring beam size, image intensity units, spatial coordinate system, spectral parameters, stokes parameters, etc. Header metadata tells the user what is in the image, and is used by imview and other tasks to set the data array on the correct spatial and spectral coordinates, assign the intensity values correctly, and otherwise properly handle the data cube.
Image Headers can be accessed and edited via the imhead task and the msmd tool. Header data can also be inspected with the casabrowser. See the page on Image Headers for further details.
Image Axes / Velocity Systems
CASA images typically have the following axis order (python indices are zero-based): Axis 0 = RA, 1 = DEC, 2 = Stokes, 3 = Frequency. The spatial axes can alternately contain GLON/GLAT or other coordinate systems. The spectral axis of images created in CASA is always in frequency units. In addition, one or more velocity systems can be added to relabel the spectral axis. When images are imported into CASA from FITS files rather than generated within CASA itself, the above conventions may not apply. See the page on Image Import and Export for further details on importing and exporting FITS files.
The spatial and spectral axes in CASA images can be modified using CASA tasks and tools described in the Reformat Images page.
Image Masks
Internal Image Masks are stored as Boolean True/False cubes within the images themselves. There can be multiple masks stored in each data cube and one of them is defined to be the ‘default’ mask. The default mask is the one visible when the image is displayed, e.g. in the CASA Viewer, and that is applied for operations on images. All masks have labels, such as mask0 etc. and they can be selected by specifying the image name followed by the mask name, separated by a colon. For example, ‘mask1’ in ‘image.im’ is used when specifying the image as ‘image.im:mask1’. Available masks can be listed with the task makemask which can also assign any mask as the default. The same task can also be used to export masks into separate CASA zero/non-zero cubes and to import such cubes as Boolean masks inside images. In addition, makemask enables the creation of masks from image regions. More information on masks is provided on the Image Masks and LEL Masks sections.
CASA Regions
CASA Regions can be specified through simple lists in LEL (e.g. region = 'box\[[108, 108,], [148, 148]]') or through CASA Region Text Format (CRTF) files, which are text files that contain one or more regions with specific shapes (e.g. ellipses and rectangles), sizes, and other properties. These files can be used to specify the region of an image in which to operate, and they can easily be modified by the user or converted to CASA image masks (Boolean data cubes) using the makemask task.
More information on CRTF files is available on the Region Files section.
Dealing with Images¶
Image cubes in CASA can be manipulated and analyzed in various ways mainly using tasks with an ‘im’ prefix and with the image CASA tool. Frequently, the tasks and tools handle CASA, FITS, and MIRIAD images, but we recommend using images in the CASA format.
In the following pages, useful image analysis tasks are introduced that span import/export tasks, image information, reformatting, mathematical operations, and spatial and spectral fitting. Available image analysis tasks include:
imhead — summarize and manipulate the “header” information in a CASA image
imsubimage — Create a (sub)image from a region of the image
imcontsub — perform continuum subtraction on a spectral-line image cube
imfit — image plane Gaussian component fitting
immath — perform mathematical operations on or between images
immoments — compute the moments of an image cube
impv — generate a position-velocity diagram along a slit
imstat — calculate statistics on an image or part of an image
imval — extract the data and mask values from a pixel or region of an image
imtrans — reorder the axes of an image or cube
imcollapse — collapse image along one or more axes by aggregating pixel values along that axis
imregrid — regrid an image onto the coordinate system of another image
imreframe — change the frame in which the image reports its spectral values
imrebin — rebin an image
specsmooth — 1-dimensional smooth images in the spectral and angular directions
imsmooth — 2-dimensional smooth images in the spectral and angular directions
specfit — fit 1-dimensional Gaussians, polynomial, and/or Lorentzians models to an image or image region
specflux — Report details of an image spectrum.
plotprofilemap — Plot spectra at their position
rmfit — Calculation of rotation measures
spxfit — Calculation of Spectral Indices and higher order polynomials
makemask — image mask handling
slsearch — query a subset of the Splatalogue spectral line catalog
splattotable — convert a file exported from Splatalogue to a CASA table
importfits — import a FITS image into a CASA image format table
exportfits — write out an image in FITS format
There are other tasks which are useful during image analysis. These include:
imview — there are useful region statistics and image cube slice and profile capabilities in the viewer
Common Task Parameters¶
Certain parameters are present in many image analysis tasks. These include:
imagename
The imagename parameter is used to specify the image(s) on which a task should operate. In most tasks, this will be a string containing the image name, but in some tasks, this can be a list of strings, as for example, in immath. Most image analysis tasks accept both CASA images and FITS images, although we recommend working with CASA images.
outfile
The outfile parameter specifies the name (in string format) of the file that the task should output. This parameter is only present in tasks that produce processed files (typically images) as output. It will therefore not be present for tasks that return python dictionaries, arrays, or other data types.
axes
The axes parameter is used to specify the image axes that the task should operate on, and the user should input a list of integers for this (e.g. “axes = [0,1]”). CASA images typically have the following axis order (python indices are zero-based): Axis 0 = RA, 1 = DEC, 2 = Stokes parameter, and 3 = Frequency. The imhead task can be used to confirm the axis specifications in the data cube of interest, and the axes may differ from the above sequence, particularly when using FITS data cubes or CASA images that were converted from FITS files. In the examples, we assume the above axis order.
To obtain statistics across RA and DEC for each velocity channel, the user would run the imstat task (imstat stands for “image statistics”) with “axes = [0,1]”. To obtain statistics over the spectral axis, one would run imstat with axes = [3].
box, chans, stokes
The box, chans, and stokes parameters are used to select parts of an image cube for the task to operate on. If a box is applied, the task will operate only on a specific spatial region (e.g. box = ‘100,100,200,200’ will only operate on pixels in the range (100,100) <= (x,y) <= (200,200) ). If specific channels are specified through chans, the task will select that segment of the spectral axis (e.g. chans = ‘30~45’ will operate on channels 30 through 45). In the same way, stokes selects specific Stokes parameter axes, as e.g. stokes = ‘I’. Further detail is provided in the Image Selection Parameters section.
mask
The mask parameter tells the task to operate on specific segments of the image cube, as set by a mask. The input for the mask parameter may be a conditional statement in LEL string format (e.g. mask = ‘ “ngc5921.im > 0.5’, which selects all pixels in that image that have values larger than 0.5 and zeros out all other pixels), or may be a Boolean True/False cube or an Integer zero/non-zero cube. The task will not operate on pixels that are “masked”, or zeroed out. See the Image Masks page for more detail and examples of usage.
stretch
This parameter can be True or False, with a default value of False. Set stretch = True when applying a single-plane mask to a full image cube. As an example, if you have a mask in a single spectral channel image that you wish to apply to all spectral channels in a cube, you would “stretch” the mask over all of the channels. The mask can also be stretched over all Stokes parameter planes for polarization images.
Returned Python Dictionaries
Many image analysis tasks return python dictionaries with information that is also printed to the logger. The dictionaries can be assigned to a variable and then used later for other scripting purposes. In the following the output of imstat is assigned to the python dictionary ‘test_stats’:
CASA <20>: test_stats=imstat(imagename='test.image')
CASA <21>: test
Out[21]:
{'blc': array([0, 0, 0, 0], dtype=int32),
'blcf': '17:45:40.899, -29.00.18.780, I, 1.62457e+10Hz',
'max': array([ 0.49454519]),
'maxpos': array([32, 32, 0, 0], dtype=int32),
'maxposf': '17:45:40.655, -29.00.15.580, I, 1.62457e+10Hz',
'mean': array([ 0.00033688]),
'medabsdevmed': array([ 0.]),
'median': array([ 0.]),
'min': array([-0.0174111]),
'minpos': array([15, 42, 0, 0], dtype=int32),
'minposf': '17:45:40.785, -29.00.14.580, I, 1.62457e+10Hz',
'npts': array([ 4096.]),
'q1': array([ 0.]),
'q3': array([ 0.]),
'quartile': array([ 0.]),
'rms': array([ 0.00906393]),
'sigma': array([ 0.00905878]),
'sum': array([ 1.37985568]),
'sumsq': array([ 0.3365063]),
'trc': array([63, 63, 0, 0], dtype=int32),
'trcf': '17:45:40.419, -29.00.12.480, I, 1.62457e+10Hz'}
Image Import/Export¶
The exportfits and importfits tasks enable conversion between CASA images and FITS data. The exportfits task allows you to write your CASA image to a FITS file that other packages can read, and the importfits task converts existing FITS files into CASA images. While many image analysis tasks can operate on FITS files, we recommend converting to CASA images for processing and analysis purposes.
Export CASA Image to FITS (exportfits)
The exportfits task is used to export a CASA image to FITS format. The inputs are:
#exportfits :: Convert a CASA image to a FITS file
imagename = '' #Name of input CASA image
fitsimage = '' #Name of output image FITS
#file
velocity = False #Use velocity (rather than
#frequency) as spectral axis
optical = False #Use the optical (rather than
#radio) velocity convention
bitpix = -32 #Bits per pixel
minpix = 0 #Minimum pixel value (if
#minpix > maxpix, value is
#automatically determined)
maxpix = -1 #Maximum pixel value (if
#minpix > maxpix, value is
#automatically determined)
overwrite = False #Overwrite pre-existing
#imagename
dropstokes = False #Drop the Stokes axis?
stokeslast = True #Put Stokes axis last in
#header?
history = True #Write history to the FITS
#image?
dropdeg = False #Drop all degenerate axes (e.g.
#Stokes and/or Frequency)?
Alert: The spectral axis of CASA images is nearly always in frequency rather than velocity. Velocities are computed only as a secondary mapping of the spectral channels with respect to a rest frequency. If velocity units are desired and the user sets velocity = True, exportfits will write the spectral axis in velocity units instead of in frequency units. The exportfits task will not output a FITS file with multiple spectral coordinate systems.
As a simple example of an exportfits command, the following will write the CASA image (‘ngc5921.clean.image’) as a FITS file (‘ngc5921.clean.fits’). In this case, the default parameter values will be adopted, so that the resulting FITS file will have the same axis order, number of pixels, etc. as the original CASA image.
exportfits(imagename='ngc5921.clean.image',outfile='ngc5921.clean.fits')
In some cases, the user may wish to use the dropstokes, stokeslast, and/or dropdeg parameters in order for the FITS image to be compatible with certain external applications. The dropdeg parameter will remove the frequency axis if it has a length of one channel, and/or it will drop the Stokes axis if that has a length of one (i.e. only one Stokes parameter is present). This would be useful, for example, for continuum data so that other programs will interpret it as a 2-D image rather than a cube.
See exportfits in the Global Task List for examples in which these and other parameters are specified.
FITS Image Import (importfits)
The importfits task enables the user to import a FITS image into CASA image table format. It is not essential to generate a CASA image file if you intend to simply view the image, as the CASA viewer can read FITS images, however we recommend importing to CASA image format for analyzing images with CASA. The inputs for importfits are:
#importfits :: Convert an image FITS file into a CASA image
fitsimage = '' #Name of input image FITS file
imagename = '' #Name of output CASA image
whichrep = 0 #If fits image has multiple
#coordinate reps, choose one.
whichhdu = 0 #If its file contains
#multiple images, choose one.
zeroblanks = True #Set blanked pixels to zero (not NaN)
overwrite = False #Overwrite pre-existing imagename
defaultaxes = False #Add the default 4D
#coordinate axes where they are missing
defaultaxesvalues = [] #List of values to assign to
#added degenerate axes when
#defaultaxes=True (ra,dec,freq,stokes)
As a simple example, the following command would create a CASA image named ‘ngc5921.clean.image’ from the FITS file ‘ngc5921.clean.fits’:
importfits(fitsimage='ngc5921.clean.fits',imagename='ngc5921.clean.image')
See importfits in the Global Task List for more complex examples.
Extracting data from an image (imval)
The imval task will extract the values of the data and mask from a specified region of an image and place in the task return value as a Python dictionary. The inputs are:
#imval :: Get the data value(s) and/or mask value in an image.
imagename = '' #Name of the input image
region = '' #Image Region. Use viewer
box = '' #Select one or more box regions
chans = '' #Select the channel(spectral) range
stokes = '' #Stokes params to image (I,IV,IQU,IQUV)
Area selection using box region is detailed in the Image Selection Parameters section. By default, box=’ ‘ will extract the image information at the reference pixel on the direction axes. Plane selection is controlled by chans and stokes. By default, chans=’ ‘ and stokes=’ ‘ will extract the image information in all channels and Stokes planes.For instance,
xval = imval('myimage', box='144,144', stokes='I' )
will extract the Stokes I value or spectrum at pixel 144,144, while
xval = imval('myimage', box='134,134.154,154', stokes='I' )
will extract a 21 by 21 pixel region. Extractions are returned in NumPy arrays in the return value dictionary, plus some extra elements describing the axes and selection:
CASA <2>: xval = imval('ngc5921.demo.moments.integrated')
CASA <3>: xval
Out[3]:
{'axes': [[0, 'Right Ascension'],
[1, 'Declination'],
[3, 'Frequency'],
[2, 'Stokes']],
'blc': [128, 128, 0, 0],
'data': array([ 0.89667124]),
'mask': array([ True], dtype=bool),
'trc': [128, 128, 0, 0],
'unit': 'Jy/beam.km/s'}
extracts the reference pixel value in this 1-plane image. Note that the ‘data’ and ‘mask’ elements are NumPy arrays, not Python lists. To extract a spectrum from a cube:
CASA <8>: xval = imval('ngc5921.demo.clean.image',box='125,125')
CASA <9>: xval
Out[9]:
{'axes': [[0, 'Right Ascension'],
[1, 'Declination'],
[3, 'Frequency'],
[2, 'Stokes']],
'blc': [125, 125, 0, 0],
'data': array([ 8.45717848e-04, 1.93370355e-03, 1.53750915e-03,
2.88399984e-03, 2.38683447e-03, 2.89159478e-04,
3.16268904e-03, 9.93389636e-03, 1.88773088e-02,
3.01138610e-02, 3.14478502e-02, 4.03211266e-02,
3.82498614e-02, 3.06552909e-02, 2.80734301e-02,
1.72479432e-02, 1.20884273e-02, 6.13593217e-03,
9.04005766e-03, 1.71429547e-03, 5.22095338e-03,
2.49114982e-03, 5.30831399e-04, 4.80734324e-03,
1.19265869e-05, 1.29435991e-03, 3.75700940e-04,
2.34788167e-03, 2.72604497e-03, 1.78467855e-03,
9.74952069e-04, 2.24676146e-03, 1.82263291e-04,
1.98463408e-06, 2.02975096e-03, 9.65532148e-04,
1.68218743e-03, 2.92119570e-03, 1.29359076e-03,
-5.11484570e-04, 1.54162932e-03, 4.68662125e-04,
-8.50282842e-04, -7.91683051e-05, 2.95954203e-04,
-1.30133145e-03]),
'mask': array([ True, True, True, True, True, True, True, True, True,
True, True, True, True, True, True, True, True, True,
True, True, True, True, True, True, True, True, True,
True, True, True, True, True, True, True, True, True,
True, True, True, True, True, True, True, True, True, True], dtype=bool),
'trc': [125, 125, 0, 45],
'unit': 'Jy/beam'}
To extract a region from the plane of a cube:
CASA <13>: xval = imval('ngc5921.demo.clean.image',box='126,128,130,129',chans='23')
CASA <14>: xval
Out[14]:
{'axes': [[0, 'Right Ascension'],
[1, 'Declination'],
[3, 'Frequency'],
[2, 'Stokes']],
'blc': [126, 128, 0, 23],
'data': array([[ 0.00938627, 0.01487772],
[ 0.00955847, 0.01688832],
[ 0.00696965, 0.01501907],
[ 0.00460964, 0.01220793],
[ 0.00358087, 0.00990202]]),
'mask': array([[ True, True],
[ True, True],
[ True, True],
[ True, True],
[ True, True]], dtype=bool),
'trc': [130, 129, 0, 23],
'unit': 'Jy/beam'}
CASA <15>: print xval['data'][0][1]
0.0148777160794
In this example, a rectangular box was extracted, and you can see the order in the array and how to address specific elements.
Image Headers¶
As summarized in the CASA Images page, an image header contains information on the observation – e.g. the observing date, pointing position, object observed, etc., and the resulting image – e.g. the restoring beam size, image intensity units, spatial coordinate system, spectral parameters, stokes parameters, etc.. Header metadata can also store notes on the observation and/or calibration and image processing. The header tells the user what is in the image and is used by the CASA viewer and other tasks to set the data array on the correct spatial and spectral coordinates, assign the intensity values correctly, and otherwise properly handle the data cube.
FITS image headers can be read in CASA using the listfits task, whereas CASA image headers can be read and edited using the imhead task. Additionally, the imhistory task can be used to view the history of the image, i.e. what operations or processes have been applied to it. These three tasks are described and demonstrated below.
List the Header of a FITS image
CASA can frequently read and write image FITS files directly. Nevertheless, it is advisable to convert the images to the CASA format first with importfits for some tasks and applications.
The task listfits can be used to display the Header Data Unit (HDU) of a FITS image. The input includes only the name of the of the FITS file, as follows:
#listfits :: List the HDU and typical data rows of a fits file:
fitsfile = '' #Name of input fits file
The logger will output the full FITS HDU. The example below shows the logger output for a Digital Sky Survey Image, which we have truncated somewhat due to the length of the output:
##########################################
#####Begin Task: listfits #####
listfits(fitsfile="dss.test.fits")
read fitsfile=dss.test.fits
d 29: DATE-OBS= '1998-11-24T11:83:00' /Observation: Date/Time time.
Primary Array HDU ------>>>
d 156: DATAMIN = 2701 /GetImage: Minimum returned pixel value
value has wrong data type. erted to type double.
d 157: DATAMAX = 22189 /GetImage: Maximum returned pixel value
value has wrong data type. erted to type double.
SIMPLE = T /FITS: Compliance
BITPIX = 16 /FITS: I*2 Data
NAXIS = 2 /FITS: 2-D Image Data
NAXIS1 = 891 /FITS: X Dimension
NAXIS2 = 893 /FITS: Y Dimension
EXTEND = T /FITS: File can contain extensions
DATE = '2016-11-17' /FITS: Creation Date
ORIGIN = 'STScI/MAST' /GSSS: STScI Digitized Sky Survey
SURVEY = 'POSSII-F' /GSSS: Sky Survey
REGION = 'XP061 ' /GSSS: Region Name
PLATEID = 'A2U4 ' /GSSS: Plate ID
SCANNUM = '01 ' /GSSS: Scan Number
DSCNDNUM= '00 ' /GSSS: Descendant Number
TELESCID= 3 /GSSS: Telescope ID
BANDPASS= 35 /GSSS: Bandpass Code
COPYRGHT= 'Caltech/Palomar' /GSSS: Copyright Holder
SITELAT = 33.356 /Observatory: Latitude
SITELONG= 116.863 /Observatory: Longitude
TELESCOP= 'Oschin Schmidt - D' /Observatory: Telescope
INSTRUME= 'Photographic Plate' /Detector: Photographic Plate
EMULSION= 'IIIaF ' /Detector: Emulsion
FILTER = 'RG610 ' /Detector: Filter
PLTSCALE= 67.2 /Detector: Plate Scale arcsec per mm
PLTSIZEX= 355 /Detector: Plate X Dimension mm
PLTSIZEY= 355 /Detector: Plate Y Dimension mm
PLATERA = 144.055 /Observation: Field centre RA degrees
PLATEDEC= 69.812 /Observation: Field centre Dec degrees
PLTLABEL= 'SF07740 ' /Observation: Plate Label
DATE-OBS= '1998-11-24T11:83:00' /Observation: Date/Time
EXPOSURE= 50 /Observation: Exposure Minutes
PLTGRADE= 'A ' /Observation: Plate Grade
OBSHA = 1.28333 /Observation: Hour Angle
OBSZD = 37.9539 /Observation: Zenith Distance
AIRMASS = 1.26743 /Observation: Airmass
REFBETA = 61.7761 /Observation: Refraction Coeff
REFBETAP= -0.082 /Observation: Refraction Coeff
REFK1 = -48616.4 /Observation: Refraction Coeff
REFK2 = -148442 /Observation: Refraction Coeff
CNPIX1 = 4993 /Scan: X Corner
CNPIX2 = 10823 /Scan: Y Corner
XPIXELS = 23040 /Scan: X Dimension
YPIXELS = 23040 /Scan: Y Dimension
XPIXELSZ= 15.0295 /Scan: Pixel Size microns
YPIXELSZ= 15 /Scan: Pixel Size microns
ASTRMASK= 'xp.mask ' /Astrometry: GSC2 Mask
WCSAXES = 2 /GetImage: Number WCS axes
WCSNAME = 'DSS ' /GetImage: Local WCS approximation from full plat
RADESYS = 'ICRS ' /GetImage: GSC-II calibration using ICRS system
CTYPE1 = 'RA---TAN' /GetImage: RA-Gnomic projection
CRPIX1 = 446 /GetImage: X reference pixel
CRVAL1 = 148.97 /GetImage: RA of reference pixel
CUNIT1 = 'deg ' /GetImage: degrees
CTYPE2 = 'DEC--TAN' /GetImage: Dec-Gnomic projection
CRPIX2 = 447 /GetImage: Y reference pixel
CRVAL2 = 69.6795 /GetImage: Dec of reference pixel
CUNIT2 = 'deg ' /Getimage: degrees
CD1_1 = -0.000279458 /GetImage: rotation matrix coefficient
CD1_2 = 2.15165e-05 /GetImage: rotation matrix coefficient
CD2_1 = 2.14552e-05 /GetImage: rotation matrix coefficient
CD2_2 = 0.00027889 /GetImage: rotation matrix coefficient
OBJECT = 'data ' /GetImage: Requested Object Name
DATAMIN = 2701 /GetImage: Minimum returned pixel value
DATAMAX = 22189 /GetImage: Maximum returned pixel value
OBJCTRA = '09 55 52.730' /GetImage: Requested Right Ascension (J2000)
OBJCTDEC= '+69 40 45.80' /GetImage: Requested Declination (J2000)
OBJCTX = 5438.47 /GetImage: Requested X on plate (pixels)
OBJCTY = 11269.3 /GetImage: Requested Y on plate (pixels)
END (0,0) = 4058 (0,1) = 4058
Reading and Manipulating CASA Image Headers
CASA image headers can be accessed and edited with the imhead task. The imagename and mode are the two primary parameters in the imhead task. The imhead task can be run with mode=’summary’, ‘list’, ‘get’, ‘put’, ‘add’, ‘del’, or ‘history’, and setting the mode opens up mode-specific sub-parameters. Many of these modes are described below and further documented in the imhead page of the Global Task List.
The default mode is mode=’summary’, which prints a summary of the image properties to the logger and terminal, and returns a dictionary containing header information. With mode=’summary’, imhead has the following inputs:
#imhead :: List, get and put image header parameters
imagename = '' #Name of the input image
mode = 'summary' #imhead options: add, del,
#get, history, list, put, summary
verbose = False #Give a full listing of
#beams or just a short summary?
#Only used when the image has multiple beams
#and mode='summary'.
Note that to capture the dictionary, it must be assigned as a Python variable, e.g. by running:
header_summary = imhead('ngc5921.demo.cleanimg.image',mode='summary')
Setting mode=’list’ prints all header keywords and values to the logger and terminal, and returns a dictionary containing the keywords and values. This mode does not have any sub-parameters.
The mode=’get’ setting allows the user to retrieve the value for a specified keyword hdkey:
#imhead :: List, get and put image header parameters
imagename = '' #Name of the input image
mode = 'get' #imhead options: list, summary, get, put
hdkey = '' #The FITS keyword
The mode=’put’ setting allows the user to replace the current value for a given keyword hdkey with that specified in hdvalue. There are two sub-parameters that are opened by this option:
#imhead :: List, get and put image header parameters
imagename = '' #Name of the input image
mode = 'put' #imhead options: list, summary, get, put
hdkey = '' #The FITS keyword
hdvalue = '' #Value of hdkey
Alert: Be careful when using mode=’put’. This task does not check whether the values you specify (e.g. for the axes types) are valid, and you can render your image invalid. Make sure you know what you are doing when using this option!
Examples for imhead
In the command below, we print the header summary to the logger:
CASA <51>: imhead('ngc5921.demo.cleanimg.image',mode='summary')
The logger output is the following:
#####Begin Task: imhead #####
Image name : ngc5921.demo.cleanimg.image
Object name : N5921_2
Image type : PagedImage
Image quantity : Intensity
Pixel mask(s) : None
Region(s) : None
Image units : Jy/beam
Restoring Beam : 52.3782 arcsec, 45.7319 arcsec, -165.572 deg
Direction reference : J2000
Spectral reference : LSRK
Velocity type : RADIO
Rest frequency : 1.42041e+09 Hz
Pointing center : 15:22:00.000000 +05.04.00.000000
Telescope : VLA
Observer : TEST
Date observation : 1995/04/13/00:00:00
Telescope position: [-1.60119e+06m, -5.04198e+06m, 3.55488e+06m] (ITRF)
Axis Coord Type Name Proj Shape Tile Coord value at pixel Coord incr Units
------------------------------------------------------------------------------------------------
0 0 Direction Right Ascension SIN 256 64 15:22:00.000 128.00 -1.500000e+01 arcsec
1 0 Direction Declination SIN 256 64 +05.04.00.000 128.00 1.500000e+01 arcsec
2 1 Stokes Stokes 1 1 I
3 2 Spectral Frequency 46 8 1.41279e+09 0.00 2.4414062e+04 Hz
Velocity 1607.99 0.00 -5.152860e+00 km/s
#####End Task: imhead
If the beam size per plane differs (for example, in a spectral data cube), the beam information will be displayed for the channel with the largest beam (i.e. the lowest frequency channel), the chennel with the smallest beam (i.e. the highest frequency channel), and the channel closest to the median beam size. If you set verbose=True, the beam information would be provided for each spectral channel (or each plane of the image). Running imhead with mode=’summary’ and verbose=False for a spectral data cube would print information on the restoring beams as follows:
Restoring Beams
Pol Type Chan Freq Vel
I Max 0 9.680e+08 0 39.59 arcsec x 22.77 arcsec pa=-70.57 deg
I Min 511 1.990e+09 -316516 20.36 arcsec x 12.05 arcsec pa=-65.67 deg
I Median 255 1.478e+09 -157949 27.11 arcsec x 15.54 arcsec pa=-70.36 deg
Setting mode=’list’ prints all header keywords and values to the logger and terminal, and returns a dictionary containing the keywords and values. In the following, we capture the resulting dictionary in the variable hlist, and print the variable.
CASA <52>: hlist = imhead('ngc5921.demo.cleanimg.image',mode='list')
CASA <53>: hlist
Out[53]:
{'beammajor': 52.378242492675781,
'beamminor': 45.731891632080078,
'beampa': -165.5721435546875,
'bunit': 'Jy/beam',
'cdelt1': '-7.27220521664e-05',
'cdelt2': '7.27220521664e-05',
'cdelt3': '1.0',
'cdelt4': '24414.0625',
'crpix1': 128.0,
'crpix2': 128.0,
'crpix3': 0.0,
'crpix4': 0.0,
'crval1': '4.02298392585',
'crval2': '0.0884300154344',
'crval3': 'I',
'crval4': '1412787144.08',
'ctype1': 'Right Ascension',
'ctype2': 'Declination',
'ctype3': 'Stokes',
'ctype4': 'Frequency',
'cunit1': 'rad',
'cunit2': 'rad',
'cunit3': '',
'cunit4': 'Hz',
'datamax': ' Not Known ',
'datamin': -0.010392956435680389,
'date-obs': '1995/04/13/00:00:00',
'equinox': 'J2000',
'imtype': 'Intensity',
'masks': ' Not Known ',
'maxpixpos': array([134, 134, 0, 38], dtype=int32),
'maxpos': '15:21:53.976, +05.05.29.998, I, 1.41371e+09Hz',
'minpixpos': array([117, 0, 0, 21], dtype=int32),
'minpos': '15:22:11.035, +04.31.59.966, I, 1.4133e+09Hz',
'object': 'N5921_2',
'observer': 'TEST',
'projection': 'SIN',
'reffreqtype': 'LSRK',
'restfreq': [1420405752.0],
'telescope': 'VLA'}
The values for these keywords can be queried using mode=’get’. In the following examples, we capture the return value:
CASA <53>: mybmaj = imhead('ngc5921.demo.cleanimg.image',mode='get',hdkey='beammajor')
CASA <54>: mybmaj
Out[54]: {'unit': 'arcsec', 'value': 52.378242492699997}
CASA <55>: myobserver = imhead('ngc5921.demo.cleanimg.image',mode='get',hdkey='observer')
CASA <56>: print myobserver
{'value': 'TEST', 'unit': ''}
You can set the values for keywords using mode=’put’. For example:
CASA <57>: imhead('ngc5921.demo.cleanimg.image',mode='put',hdkey='observer',hdvalue='CASA')
Out[57]: 'CASA'
CASA <58>: imhead('ngc5921.demo.cleanimg.image',mode='get',hdkey='observer')
Out[58]: {'unit': '', 'value': 'CASA'}
Image History (imhistory)
Image headers contain records of the operations applied to them, as CASA tasks append the image header with a record of what they did. This information can be retrieved via the imhistory task, and new messages can be appended using the imhistory task as well. The primary inputs are imagename and mode, with sub-parameters arising from the selected mode. To view the history of the image, the inputs are:
#imhistory :: Retrieve and modify image history
imagename = '' #Name of the input image
mode = 'list' #Mode to run in, 'list' to
#retrieve history,'append'
#to append a record to history.
verbose = True #Write history to logger if
#mode='list'?
With verbose=True (default) the image history is also reported in the CASA logger. The imhistory task returns the messages in a Python list that can be captured by a variable, e.g.
myhistory=imhistory('image.im')
It is also possible to add messages to the image headers via mode=’append’. See the `imhistory <../api/casatasks.rst>`__ page in the Global Task List for an example of appending messages to the image history.
Reformat Images¶
This section contains a description of various tasks that reformat images. These include:
imsubimage: enables the user to extract a sub-image from a larger cube,
imtrans: changes the axis order in an image,
imregrid: sets the image onto a different spatial coordinate system or spectral grid,
imreframe: changes the velocity system of an image
imrebin: rebins an image in a spatial or spectral dimension
imcollapse: collapses an image along an axis.
Extracting sub-images
The task imsubimage provides a way to extract a smaller data cube from a bigger one. The inputs are:
#imsubimage :: Create a (sub)image from a region of the image
imagename = '' #Input image name. Default is unset.
outfile = '' #Output image name. Default is unset.
box = '' #Rectangular region to select in
#direction plane. Default is to use the
#entire direction plane.
region = '' #Region selection. Default is to use the
#full image.
chans = '' #Channels to use. Default is to use all
#channels.
stokes = '' #Stokes planes to use. Default is to use
#all Stokes planes.
mask = '' #Mask to use. Default is none.
dropdeg = True #Drop degenerate axes
[ keepaxes = [] #If dropdeg=True, these are the]
#degenerate axes to keep. Nondegenerate
#axes are implicitly always kept.
verbose = True #Post additional informative messages to
#the logger
The region keyword defines the size of the smaller cube and is specified via the CASA region CRTF syntax. E.g.
region='box [ [ 100pix , 130pix] , [120pix, 150pix ] ]'
will extract the portion of the image that is between pixel coordinates (100,130) and (12,150). dropdeg=T with selection via keepaxes is useful to remove axes in the data cube that are degenerate, i.e. axes with a single plane only. A single Stokes I axis is a common example.
Reordering the Axes of an Image Cube
Sometimes data cubes can be in axis orders that are not adequate for processing. The CASA task imtrans can change the ordering of the axis:
#imtrans :: Reorder image axes
imagename = '' #Name of the input image
outfile = '' #Name of output CASA image.
order = '' #New zero-based axes order.
wantreturn = True #Return an image tool referencing the
#transposed image
The order parameter is the most important input here. It is a string of numbers that shows how axes 0, 1, 2, 3, … are mapped onto the new cube (note that the first axis has the label 0, as typical in python). E.g. order=’1032’ will reorder the input axis 0 to be axis 1 in the output, input axis 1 to be output axis 0, input axis 2 to output axis 3 (the last axis) and input axis 3 to output axis 2. Alternatively, axes can be specified by their names. E.g., to reorder an image with right ascension, declination, and frequency and reverse the first two, order=[‘declination’, ‘right ascension’, ‘frequency’] will work. The axes names can be found typing ia.coordsys.names. Minimum match is supported, so that order=[‘d’, ‘f’, ‘r’] will produce the same results.Axes can simultaneously be transposed and reversed. To reverse an axis, precede it by a ‘-‘. For example, order=’-10-32’ will reverse the direction of the first and third axis of the input image (the zeroth and second axes in the output image).Example (swap the stokes and spectral axes in an RA-Dec-Stokes-Frequency image):
imagename = 'myim.im'
outfile = 'outim.im'
order = '0132'
imtrans()
or
outfile = 'myim_2.im'
order = 132
imtrans()
or
outfile = 'myim_3.im'
order = ['r', 'd', 'f', 's']
imtrans()
or
outfile = 'myim_4.im'
order = ['rig', 'declin', 'frequ', 'stok']
imtrans()
If the outfile parameter is empty, only a temporary image is created; no output image is written to disk. The temporary image can be captured in the returned value (assuming wantreturn=True).
Regridding an Image (imregrid)
Inside the Toolkit: More complex coordinate system and image regridding operation can be carried out in the toolkit. The coordsys (cs) tool and the ia.regrid method are the relevant components.
It is occasionally necessary to regrid an image onto a new coordinate system. The imregrid task will regrid one image onto the coordinate system of another, creating an output image. In this task, the user need only specify the names of the input, template, and output images. The default inputs are:
#imregrid :: regrid an image onto a template image
imagename = '' #Name of the source image
template = 'get' #A dictionary, refcode, or name of an
#image that provides the output shape
#and coordinate system
output = '' #Name for the regridded image
asvelocity = True #Regrid spectral axis in velocity space
#rather than frequency space?
axes = [-1] #The pixel axes to regrid. -1 => all.
interpolation = 'linear' #The interpolation method. One of
#'nearest', 'linear', 'cubic'.
decimate = 10 #Decimation factor for coordinate grid
#computation
replicate = False #Replicate image rather than regrid?
overwrite = False #Overwrite (unprompted) pre-existing
#output file?
The output image will have the data in imagename regridded onto the coordinate system provided by the template parameter. template is used universally for a range of ways to define the grid of the output image:
a template image: specify an image name here and the input will be regridded to the same 3-dimensional coordinate system as the one in template. Values are filled in as blanks if they do not exist in the input. Note that the input and template images must have the same coordinate structure to begin with (like 3 or 4 axes, with the same ordering)
a coordinate system (reference code): to convert from one coordinate frame to another one, e.g. from B1950 to J2000, the template parameter can be used to specify the output coordinate system. These following values are supported: ‘J2000’, ‘B1950’, ‘B1950_VLA’, ‘GALACTIC’, ‘HADEC’, ‘AZEL’, ‘AZELSW’, ‘AZELNE’, ‘ECLIPTIC’, ‘MECLIPTIC’, ‘TECLIPTIC’, ‘SUPERGAL’
‘get’: This option returns a python dictionary in the {‘csys’: csys_record, ‘shap’: shape} format
a python dictionary: In turn, such a dictionary can be used as a template to define the final grid
Redefining the Velocity System of an Image
imreframe can be used to change the velocity system of an image. It is not applying a regridding as a change from radio to optical conventions would require, but it will change the labels of the velocity axes.
#imreframe :: Change the frame in which the image reports its spectral values
imagename = '' #Name of the input image
output = '' #Name of the output image; '' => modify input image
outframe = 'lsrk' #Spectral frame in which the frequency or velocity
#values will be reported by default
restfreq = '' #restfrequency to use for velocity values (e.g.
#'1.420GHz' for the HI line)
outframe defines the velocity frame (LSRK, BARY, etc.,) of the output image and a rest frequency should be specified to relabel the spectral axis in new velocity units.
Rebin an Image
The task imrebin allows one to rebin an image in any spatial or spectral direction:
imrebin :: Rebin an image by the specified integer factors
imagename = '' #Name of the input image
outfile = '' #Output image name.
factor = [] #Binning factors for each axis. Use
#imhead or ia.summary to determine axis
#ordering.
region = '' #Region selection. Default is to use the full
#image.
box = '' #Rectangular region to select in
#direction plane. Default is to use the entire
#direction plane.
chans = '' #Channels to use. Default is to use all
#channels.
stokes = '' #Stokes planes to use. Default is to
#use all Stokes planes. Stokes planes
#cannot be rebinned.
mask = '' #Mask to use. Default is none.
dropdeg = False #Drop degenerate axes?
crop = True #Remove pixels from the end of an axis to
#be rebinned if there are not enough to
#form an integral bin?
where factor is a list of integers that provides the numbers of pixels to be binned for each axis. The crop parameters controls how pixels at the boundaries are treated if the bin values are not multiple integers of the image dimensions.Example:
imrebin(imagename='my.im', outfile='rebinned.im', factor=[1,2,1,4], crop=T)
This leaves RA untouched, bins DEC by a factor of 2, leaves Stokes as is, and bins the spectral axis by a factor of 4. If the input image has a spectral axis with a length that is not a multiple of 4, the crop=T setting will discard the remaining 1-3 edge pixels.
Collapsing an Image Along an Axis
imcollapse allows to apply an aggregation function along one or more axes of an image. Functions supported are ‘max’, ‘mean’, ‘median’, ‘min’, ‘rms’, ‘stdev’, ‘sum’, ‘variance’ (minimum match supported). The relevant axes will then collapse to a single value or plane (i.e. they will result in a degenerate axis). The functions are specified in the function parameter of the imcollapse inputs:
#imcollapse :: Collapse image along one axis, aggregating pixel values along that axis.
imagename = '' #Name of the input image
function = '' #Function used to compute aggregation
#of pixel values.
axes = [0] #Zero-based axis number(s) or minimal
#match strings to collapse.
outfile = '' #Name of output CASA image.
box = '' #Optional direction plane box ('blcx,
#blcy, trcx trcy').
region = '' #Name of optional region file to use.
chans = '' #Optional zero-based contiguous
#frequency channel specification.
stokes = '' #Optional contiguous stokes planes
#specification.
mask = '' #Optional mask to use.
wantreturn = True #Should an image analysis tool
#referencing the collapsed image be
#returned?
wantreturn=True returns an image analysis tool containing the newly created collapsed image.Example (myimage.im is a 512x512x128x4 (ra,dec,freq,stokes; i.e. in the 0-based system, frequency is labeled as axis 2) image and we want to collapse a subimage of it along its spectral axis avoiding the 8 edge channels at each end of the band, computing the mean value of the pixels (resulting image is 256x256x1x4 in size)):
imcollapse(imagename='myimage.im', outfile='collapse_spec_mean.im',
function='mean', axis=2, box='127,127,383,383', chans='8~119')
Note that imcollapse will not smooth to a common beam for all planes if they differ. If this is desired, run imsmooth before imcollapse.
Spectral Analysis¶
Continuum Subtraction on an Image Cube (imcontsub)
One method to separate line and continuum emission in an image cube is to specify a number of line-free channels in that cube, make a linear fit to the visibilities in those channels, and subtract the fit from the whole cube. Note that the task uvcontsub serves a similar purpose but the subtraction is performed in visibility space (see UV Continuum Subtraction. The imcontsub task will perform a polynomial baseline fit to the specified channels from an image cube and subtract it from all channels. The default inputs are:
#imcontsub :: Continuum subtraction on images
imagename = '' #Name of the input image
linefile = '' #Output line image file name
contfile = '' #Output continuum image file name
fitorder = 0 #Polynomial order for the continuum estimation
region = '' #Image region or name to process see viewer
box = '' #Select one or more box regions
chans = '' #Select the channel(spectral) range
stokes = '' #Stokes params to image (I,IV,IQU,IQUV)
ALERT: imcontsub has issues when the image does not contain a spectral or stokes axis. Errors are generated when run on an image missing one or both of these axes. You will need to use the toolkit (e.g. the ia.adddegaxes method) to add degenerate missing axes to the image.
Examples for imcontsub
For example, we first make a clean image from data in which no uv-plane continuum subtraction has been performed:
#Now clean, keeping all the channels except first and last
default('clean')
vis = 'ngc5921.demo.src.split.ms'
imagename = 'ngc5921.demo.nouvcontsub'
mode = 'channel'
nchan = 61
start = 1
width = 1
imsize = [256,256]
psfmode = 'clark'
imagermode = ''
cell = [15.,15.]
niter = 6000
threshold='8.0mJy'
weighting = 'briggs'
robust = 0.5
mask = [108,108,148,148]
interactive=False
clean()
#It will have made the image:
#-----------------------------
#ngc5921.demo.nouvcontsub.image
#You can view this image
imview('ngc5921.demo.nouvcontsub.image')
Channels 0 through 4 and 50 through 60 are line-free. Continuum subtraction is then performed with:
default('imcontsub')
imagename = 'ngc5921.demo.nouvcontsub.image'
linefile = 'ngc5921.demo.nouvcontsub.lineimage'
contfile = 'ngc5921.demo.nouvcontsub.contimage'
fitorder = 1
chans = '0~4,50~60'
stokes = 'I'
imcontsub()
Computing the Moments of an Image Cube
For spectral line datasets, the output of the imaging process is an image cube, with a frequency or velocity channel axis in addition to the two sky coordinate axes. This can be most easily thought of as a series of image planes stacked along the spectral dimension. A useful product to compute is to collapse the cube into a moment image by taking a linear combination of the individual planes:
\(M_m(x_i,y_i) = \sum_k^N w_m(x_i,y_i,v_k)\,I(x_i,y_i,v_k)\)
for pixel i and channel k in the cube I. There are a number of choices to form the moment-m, usually approximating some polynomial expansion of the intensity distribution over velocity mean or sum, gradient, dispersion, skew, kurtosis, etc. There are other possibilities (other than a weighted sum) for calculating the image, such as median filtering, finding minima or maxima along the spectral axis, or absolute mean deviations. And the axis along which to do these calculations need not be the spectral axis (i.e. do moments along Dec for a RA-Velocity image). We will treat all of these as generalized instances of a “moment” map.The immoments task will compute basic moment images from a cube. The default inputs are:
#immoments :: Compute moments of an image cube:
imagename = '' #Input image name
moments = [0] #List of moments you would like to compute
axis = 'spectral' #The moment axis: ra, dec, lat, long, spectral, or stokes
region = '' #Image Region. Use viewer
box = '' #Select one or more box regions
chans = '' #Select the channel(spectral) range
stokes = '' #Stokes params to image (I,IV,IQU,IQUV)
mask = '' #mask used for selecting the area of the
#image to calculate the moments on
includepix = -1 #Range of pixel values to include
excludepix = -1 #Range of pixel values to exclude
outfile = '' #Output image file name (or root for multiple moments)
This task will operate on the input file given by imagename and produce a new image or set of images based on the name given in outfile.The moments parameter chooses which moments are calculated. The choices for the operation mode are:
moments=-1 - mean value of the spectrum
moments=0 - integrated value of the spectrum
moments=1 - intensity weighted coordinate; traditionally used to get ‘velocity fields’
moments=2 - intensity weighted dispersion of the coordinate; traditionally used to get ‘velocity dispersion’
moments=3 - median of I
moments=4 - median coordinate
moments=5 - standard deviation about the mean of the spectrum
moments=6 - root mean square of the spectrum
moments=7 - absolute mean deviation of the spectrum
moments=8 - maximum value of the spectrum
moments=9 - coordinate of the maximum value of the spectrum
moments=10 - minimum value of the spectrum
moments=11 - coordinate of the minimum value of the spectrum
The meaning of these is described in the CASA Toolkit Manual, that describes the associated image.moments tool.The axis parameter sets the axis along which the moment is “collapsed” or calculated. Choices are: ‘ra’, ‘dec’, ‘lat’, ‘long’, ‘spectral’, or ‘stokes’. A standard moment-0 or moment-1 image of a spectral cube would use the default choice ‘spectral’. One could make a position-velocity map by setting ‘ra’ or ‘dec’.The includepix and excludepix parameters are used to set ranges for the inclusion and exclusion of pixels based on values. For example, includepix=[0.05,100.0] will include pixels with values from 50 mJy to 1000 Jy, and excludepix=[100.0,1000.0] will exclude pixels with values from 100 to 1000 Jy.If a single moment is chosen, the outfile specifies the exact name of the output image. If multiple moments are chosen, then outfile will be used as the root of the output filenames, which will get different suffixes for each moment. For image cubes that contain different beam sizes for each plane, immoments will smooth all planes to the largest beam size first, then collapse to the desired moment.
Hints for using immoments
In order to make an unbiased moment-0 image, do not put in any thresholding using includepix or excludepix. This is so that the (presumably) zero-mean noise fluctuations in off-line parts of the image cube will cancel out. If you image has large biases, like a pronounced clean bowl due to missing large-scale flux, then your moment-0 image will be biased also. It will be difficult to alleviate this with a threshold, but you can try.
To make a usable moment-1 (or higher) image, on the other hand, it is critical to set a reasonable threshold to exclude noise from being added to the moment maps. Something like a few times the rms noise level in the usable planes seems to work (put into includepix or excludepix as needed). Also use chans to ignore channels with bad data.
Examples using immoments
default('immoments')
imagename = 'ngc5921.demo.cleanimg'
#Do first and second spectral moments
axis = 'spectral'
chans = ''
moments = [0,1]
#Need to mask out noisy pixels, currently done
#using hard global limits
excludepix = [-100,0.009]
outfile = 'ngc5921.demo.moments'
immoments()
#It will have made the images:
#--------------------------------------
#ngc5921.demo.moments.integrated
#ngc5921.demo.moments.weighted_coord
Other examples of NGC2403 (a moment-0 image of a VLA line dataset) and NGC4826 (a moment-1 image of a BIMA CO line dataset) are shown in the Figure below
NGC2403 VLA moment zero (left) and NGC4826 BIMA moment one (right) images as shown in the viewer.
Generating Position-Velocity Diagrams (impv)
CASA can generate position-velocity (PV) diagrams via the task impv:
#impv :: Construct a position-velocity image by choosing two points in the direction plane.
imagename = '' #Name of the input image
outfile = '' #Output image name. If empty, no image is written.
mode = 'coords' #If 'coords', use start and end values. If 'length', use
#center, length, and pa values.
width = 1 #Width of slice for averaging pixels perpendicular to the
#slice. Must be an odd positive integer or valid
#quantity. See help for details.
unit = 'arcsec' #Unit for the offset axis in the resulting image. Must be
#a unit of angular measure.
chans = '' #Channels to use.
#Channels must be contiguous. Default is to use all
#channels.
region = '' #Region selection. Default is entire image. No selection
#is permitted in the direction plane.
stokes = 'I' #Stokes planes to use. Planes must be contiguous. Default
#is to use all stokes.
mask = [] #Mask to use. Default is none.
stretch = False #Stretch the mask if necessary and possible? Default False
PV diagrams are generated by “slicing” a datacube through the RA/DEC planes. The “slit” can be defined either by start/end coordinates or by a length, center coordinate, and position angle. Averaged over the width of the ‘slit’ the image cube values are then stored in a new image with position and velocity as the two axes. The slit position is specified by a start and end pixel in the RA/DEC plane of the data cube. An angular unit can be set to define what is stored in the resulting PV image.
1-dimensional Smoothing (specsmooth)
To gain higher signal-to-noise of data cubes, one can smooth the data along one dimension. Typically this is the spectral axis. Hanning and Boxcar smoothing kernels are available in the task specsmooth:
#specsmooth :: Smooth an image region in one dimension
imagename = '' #Name of the input image
outfile = '' #Output image name.
region = '' #Region selection. Default is to use the full
#image.
box = '' #Rectangular region to select in
#direction plane. Default is to use the entire
#direction plane.
mask = '' #Mask to use. Default is none..
axis = -1 #The profile axis. Default: use the
#spectral axis if one exists, axis 0
#otherwise (<0).
function = 'hanning' #Convolution function. hanning and boxcar
#are supported functions. Minimum match
#is supported.
dmethod = 'copy' #Decimation method. '' means no
#decimation, 'copy' and 'mean' are also
#supported (minimum match).
The parameter dmethod=’copy’ allows one to only keep every nth channel, if the smoothing kernel has a width of n. Leaving this parameter empty will return the same size cube as the input and setting it to ‘mean’ will average planes using the kernel width.
Spectral Line fitting
specfit is a powerful task to perform spectral line fits in data cubes. Three types of fitting functions are currently supported, polynomials, Gaussians, and Lorentzians. specfit can fit these functions in two ways: over data that were averaged across a region (multifit=False) or on a pixel by pixel basis (multifit=True).
#specfit :: Fit 1-dimensional Gaussians and/or polynomial models to an image or image region
imagename = '' #Name of the input image
box = '' #Rectangular box in direction coordinate
#blc, trc. Default: entire image ('').
region = '' #Region of interest. Default: Do
#not use a region.
chans = '' #Channels to use. Channels must be
#contiguous. Default: all channels ('').
stokes = '' #Stokes planes to use. Planes must be
#contiguous. Default: all stokes ('').
axis = -1 #The profile axis. Default: use the
#spectral axis if one exists, axis 0
#otherwise (<0).
mask = '' #Mask to use. Default is
#none..
poly = -1 #Order of polynomial element. Default: do
#not fit a polynomial (<0).
estimates = '' #Name of file containing initial estimates.
#Default: No initial estimates ('').
ngauss = 1 #Number of Gaussian elements. Default: 1.
pampest = '' #Initial estimate of PCF profile (gaussian
#or lorentzian) amplitudes.
pcenterest = '' #Initial estimate PCF profile centers, in
#pixels.
pfwhmest = '' #Initial estimate PCF profile FWHMs, in
#pixels.
pfix = '' #PCF profile parameters to fix during fit.
pfunc = '' #PCF singlet functions to fit. 'gaussian'
#or 'lorentzian' (minimal match
#supported). Unspecified means all
#gaussians.
minpts = 0 #Minimum number of unmasked points
#necessary to attempt fit.
multifit = True #If true, fit a profile along the desired
#axis at each pixel in the specified
#region. If false, average the non-fit
#axis pixels and do a single fit to that
#average profile. Default False.
amp = '' #Name of amplitude solution image. Default:
#do not write the image ('').
amperr = '' #Name of amplitude solution error image.
#Default: do not write the image ('').
center = '' #Name of center solution image. Default: do
#not write the image ('').
centererr = '' #Name of center solution error image.
#Default: do not write the image ('').
fwhm = '' #Name of fwhm solution image. Default: do
#not write the image ('').
fwhmerr = '' #Name of fwhm solution error image.
#Default: do not write the image ('').
integral = '' #Prefix of name of integral solution image.
#Name of image will have gaussian
#component number appended. Default: do
#not write the image ('').
integralerr = '' #Prefix of name of integral error solution
#image. Name of image will have gaussian
#component number appended. Default: do
#not write the image ('').
model = '' #Name of model image. Default: do not write
#the model image ('').
residual = '' #Name of residual image. Default: do not
#write the residual image ('').
wantreturn = True #Should a record summarizing the results be
#returned?
logresults = True #Output results to logger?
gmncomps = 0 #Number of components in each gaussian
#multiplet to fit
gmampcon = '' #The amplitude ratio constraints for non-
#reference components to reference
#component in gaussian multiplets.
gmcentercon = '' #The center offset constraints (in pixels)
#for non-reference components to reference
#component in gaussian multiplets.
gmfwhmcon = '' #The FWHM ratio constraints for non-
#reference components to reference
#component in gaussian multiplets.
gmampest = [0.0] #Initial estimate of individual gaussian
#amplitudes in gaussian multiplets.
gmcenterest = [0.0] #Initial estimate of individual gaussian
#centers in gaussian multiplets, in
#pixels.
gmfwhmest = [0.0] #Initial estimate of individual gaussian
#FWHMss in gaussian multiplets, in pixels.
gmfix = '' #Parameters of individual gaussians in
#gaussian multiplets to fix during fit.
logfile = '' #File in which to log results. Default is
#not to write a logfile.
goodamprange = [0.0] #Acceptable amplitude solution range. [0.0]
#=> all amplitude solutions are
#acceptable.
goodcenterrange = [0.0] #Acceptable center solution range in pixels
#relative to region start. [0.0] => all
#center solutions are acceptable.
goodfwhmrange = [0.0] #Acceptable FWHM solution range in pixels.
#[0.0] => all FWHM solutions are
#acceptable.
sigma = '' #Standard deviation array or image name.
Polynomial Fits
Polynomials can be fit by specifying the polynomial order in poly. Negative orders will not fit any polynomials.
Lorentzian and Gaussian Fits
Gaussian and Lorentzian fits are very similar, they both require amplitude, center, and FWHM to be fully specified. All of the following discussion is thus valid for both functions. The parameter pfunc controls whether Gaussian or Lorentzian functions are to be used. Default is all Gaussians. pfunc=[‘L’, ‘G’, ‘G’, ‘L’] would use Lorentzian, Gaussian, Gaussian, and Lorentzian components in the order they appear in the estimates file (see below).
One or more single Gaussian/Lorentzian
For Gaussian and Lorentzian fits, the task will allow multiple components and specfit will try to find the best solution. The parameter space, however, is usually not uniform and to avoid local minima in the goodness-of-fit space, one can provide initial start values for the fits. This can be done either through the parameters pampest, pcenterest, and pfwhmest for the amplitudes, center, and FWHM estimates in image coordinates. pfix can take parameters that specify fixed fit values. Any combination of the characters ‘p’ (peak), ‘c’ (center), and ‘f’ (fwhm) are permitted, e.g. ‘fc’ will hold the fwhm and the center constant during the fit. Fixed parameters will have no errors associated with them in the solution. Alternatively, a file with initial values can be supplied by the estimates parameter (one Gaussian/Lorentzian parameter set per line). The file has the following format:
[peak intensity], [center], [fwhm], [optional fixed parameter string]
The first three values are required and must be numerical values. The peak intensity must be expressed in map units, while the center and fwhm must be specified in pixels. The fourth value is optional and if present, represents the parameter(s) that should be held constant during the fit (see above). An example estimates file is:
# estimates file indicating that two Gaussians should be fit
# first guassian estimate, peak=40, center at pixel number 10.5,
# fwhm = 5.8 pixels, all parameters allowed to vary during
# fit
40, 10.5, 5.8
# second Gaussian, peak = 4, center at pixel number 90.2,
# fwhm = 7.2 pixels, hold fwhm constant
4, 90.2, 7.2, f
# end file
and the output of a typical execution, e.g.
specfit(imagename='IRC10216_HC3N.cube_r0.5.image', region='specfit.crtf',
multifit=F, estimates='', ngauss=2)
(‘specfit.crtf’ is a CASA regions file, see Section D) will be
Fit :
RA : 09:47:57.49
Dec : 13.16.46.46
Stokes : I
Pixel : [146.002, 164.499, 0.000, *]
Attempted : YES
Converged : YES
Iterations : 28
Results for component 0:
Type : GAUSSIAN
Peak : 5.76 +/- 0.45 mJy/beam
Center : -15.96 +/- 0.32 km/s
40.78 +/- 0.31 pixel
FWHM : 7.70 +/- 0.77 km/s
7.48 +/- 0.74 pixel
Integral : 47.2 +/- 6.0 mJy/beam.km/s
Results for component 1:
Type : GAUSSIAN
Peak : 4.37 +/- 0.33 mJy/beam
Center : -33.51 +/- 0.58 km/s
23.73 +/- 0.57 pixel
FWHM : 15.1 +/- 1.5 km/s
14.7 +/- 1.5 pixel
Integral : 70.2 +/- 8.8 mJy/beam.km/s
If wantreturn=True (the default value), the task returns a python dictionary (here captured in a variable with the inventive name of ‘fitresults’) :
fitresults=specfit(imagename='IRC10216_HC3N.cube_r0.5.image', region='specfit.rgn', multifit=F, estimates='', ngauss=2)
The values can then be used by other python code for further processing.
Gaussian Multiplets
It is possible to fit a number of Gaussians together, as multiplets with restrictions. All restrictions are relative to a reference Gaussian (the zero’th component of each multiplet). gncomps specifies the number of Gaussians for each multiplets, and, in fact, a number of these multiplets can be fit simultaneously. gncomps=[2,4,3], e.g. fits a 2-component Gaussian, a 4-component Gaussian, and a 3-component Gaussian all at once. The initial parameter estimates can be specified with the gmampest, gmcenterest, and gmfwhmest parameters and the estimates are simply listed in the sequence of gncomps. E.g. if gncomps=[2,4,3] is specified with multiplet G0 consisting of 2 Gaussians a, b, multiplet G1 of 4 Gaussians c, d, e, f, and multiplet G2 of three Gaussians g, h, i, the parameter list in gm*est would be like gm*est=[a,b,c,d,e,f,g,h,i].Restrictions can be specified via the gmampcon parameter for the amplitude ratio (non-reference to reference), gmcentercon for the offset in pixels (to a reference), and gmfwhmcon for the FWHM ratio (non-reference to reference). A value of 0 will not constrain anything. The reference is always the zero’th component in each multiplet, in our example, Gaussians a, c, and g. They cannot be constrained. So gmncomps=[2, 4, 3], gmampcon= [ 0 , 0.2, 0 , 0.1, 4.5, 0 ], gcentercon=[24.2, 45.6, 92.7, 0 , -22.8, -33.5], and gfwhmcon=’ ‘ would constrain Gaussians b relative to a with a 24.2 pixel offset, Gaussian d to c with a amplitude ratio of 0.2 and a 45.6 pixel offset, Gaussian e to c with a offset of 92.7 pixels, etc. Restrictions will overrule any estimates.The parameters goodamprange, goodcenterrange, and goodfwhmrange can be used to limit the range of amplitude, center, and fwhm solutions for all Gaussians.
Pixel-by-pixel fits
As mentioned above, specfit can also fit spectral cubes on a pixel by pixel basis. In this case, one can choose to write none, any or all of the solution and error images for Gaussian/Lorentzian fits via the parameters amp, amperr, center, centererr, fwhm, and fwhmerr. Each Gaussian component will produce a set of images _0, _1, etc. suffixes. Writing analogous images for polynomial coefficients is not yet supported although polynomial fits when multifit=True is supported. Best fit coefficients are written to the logger. Pixels for which fits were not attempted or did not converge will be masked as bad.
Spatial Spectral Line Properties
specflux calculates the flux as a function of frequency and velocity over a selected spatial region. Flux densities of Jy/beam are being converted to Jy by properly integrating over the selected region.The input parameters of specflux are:
#specflux :: Report details of an image spectrum.
imagename = '' #Name of the input image
box = '' #Rectangular region to select in
#direction plane. Default is to use the entire
#direction plane.
region = '' #Region selection. Default is to use the full
#image.
chans = '' #Channels to use. Default is to use all
#channels.
stokes = '' #Stokes planes to use. Default is to
#use all Stokes planes.
mask = '' #Mask to use. Default
#is none.
unit = 'km/s' #Unit to use for the abscissa. Must be
#conformant with a typical spectral axis
#unit.
major = '' #Major axis of overriding restoring beam.
#If specified, must be a valid quantity.
minor = '' #Minor axis of overriding restoring beam.
#If specified, must be a valid quantity
logfile = '' #File which to write details. Default is
#to not write to a file.
The results can be written into a logfile to be plotted in other packages.
Plot Spectra on a Map (plotprofilemap)
The profilemap task enables plotting spectra according to their pointing directions (a.k.a. a profile map) in plots. The input should be CASA image,or FITS format cube. Spectra within the cube are averaged into a bin number specified with the numpanels keyword. Absent or masked data are treated according to plotmasked keyword setting.
plotprofilemap(imagename='interesting_spectralcube_casaimge.im',
figfile = 'grid_map.png',
separatepanel=F,
spectralaxis = 'velocity',
title = 'myprofilemap',
transparent = F,
showaxislabel = True,
showtick = True,
showticklabel = True,
pol=0,
numpanels='8')
Calculation of Rotation Measures
rmfit is an image domain task that accepts an input cube image containing Stokes Q and U axes and will generate the rotation measure by performing a least square fit in the image domain to obtain the best fit to the equation \(\chi = \chi_0 + RM\times \lambda^2\).
The inputs to rmfit are:
#rmfit :: Calculate rotation measure.
imagename = '' #Name(s) of the input image(s). Must be specified.
rm = '' #Output rotation measure image name. If not specified, no
#image is written.
rmerr = '' #Output rotation measure error image name. If not
#specified, no image is written.
pa0 = '' #Output position angle (degrees) at zero wavelength image
#name. If not specified, no image is written.
pa0err = '' #Output position angle (degrees) at zero wavelength error
#image name. If not specified, no image is written.
nturns = '' #Output number of turns image name. If not specified, no
#image is written.
chisq = '' #Output reduced chi squared image name. If not specified,
#no image is written.
sigma = '' #Estimate of the thermal noise. A value less than 0 means
#auto estimate.
rmfg = 0.0 #Foreground rotation measure in rad/m/m to subtract.
rmmax = 0.0 #Maximum rotation measure in rad/m/m for which to solve.
#IMPORTANT TO SPECIFY.
maxpaerr = 1e+30 #Maximum input position angle error in degrees to allow in
#solution determination.
This task generates the rotation measure image from stokes Q and U measurements at several different frequencies. You are required to specify the name of at least one image with a polarization axis containing stokes Q and U planes and with a frequency axis containing more than two pixels. The frequencies do not have to be equally spaced (i.e. the frequency coordinate can be a tabular coordinate). It will work out the position angle images for you. You may also specify multiple image names, in which case these images will first be concatenated along the spectral axis using ia.imageconcat. The requirements are that for all images, the axis order must be the same and the number of pixels along each axis must be identical, except for the spectral axis which may differ in length between images. The spectral axis need not be contiguous from one image to another. See also the imagepol.fourierrotationmeasure function for a new Fourier-based approach.Rotation measure algorithms that work robustly are few. The main problem is in trying to account for the \(n- \pi\) ambiguity (see Leahy et al.1986 - Appendix A.1) [1] and the MIRIAD manual.
But as in all these algorithms, the basic process is that for each spatial pixel, the position angle vs frequency data is fit to determine the rotation measure and the position angle at zero wavelength (and associated errors). An image containing the number of \(n- \pi\) turns that were added to the data at each spatial pixel and for which the best fit was found can be written. The reduced \(\chi^2\) image for the fits can also be written. Any combination of output images can be written.
NOTE: No assessment of curvature (i.e. deviation from the simple linear position angle - \(lambda^2\) functional form) is made.
The parameter sigma gives the thermal noise in Stokes Q and U. By default it is determined automatically using the image data. But if it proves to be inaccurate (maybe not many signal-free pixels), it may be specified. This is used for calculating the error in the position angles (via propagation of Gaussian errors).The argument maxpaerr specifies the maximum allowable error in the position angle that is acceptable. The default is an infinite value. From the standard propagation of errors, the error in the linearly polarized position angle is determined from the Stokes Q and U images (at each directional pixel for each frequency). If the position angle error for any pixel exceeds the specified value, the position angle at that pixel is omitted from the fit. The process generates an error for the fit and this is used to compute the errors in the output images.
NOTE: maxpaerr is not used to mask pixels in the output images.
The argument rmfg is used to specify a foreground RM value. For example, you may know the mean RM in some direction out of the Galaxy, then including this can improve the algorithm by reducing ambiguity. The parameter rmmax specifies the maximum absolute RM value that should be solved for. This quite an important parameter. If you leave it at the default, zero, no ambiguity handling will be used. So some a priori information should be supplied; this is the basic problem with rotation measure algorithms.
Calculation of Spectral Indices and Higher Order Polynomials
This application fits a power logarithmic polynomial or a logarithmic transformed polynomial to pixel values along a specified axis of an image or images. These functions are most often used for fitting the spectral index and higher order terms of a spectrum. A power logarithmic polynomial has the form
\(y = \frac{c_0 x}{D^{(c_1 + c_2 \ln(x/D) + c_3 \ln(x/D)^2 + c_n \ln(x/D)^{(n - 1)})}}\)
and a logarithmic transformed polynomial is simply the result of this equation after taking the natural log of both sides so that it has the form
\(\ln(y) = c_0 + c_1 \ln(x) + c_2 \ln(x/D)^2 + ... + c_n \ln(x/D)^n\)
Because the logarithm of the ordinate values must be taken before fitting a logarithmic transformed polynomial, all non-positive pixel values are effectively masked for the purposes of fitting. The coefficients of the two forms are equal to each other except that c0 in the second equation is equal to \(\ln(c_0)\) of the first. In the case of fitting a spectral index, which is traditionally represented as \(\alpha\), is equal to \(c_1\). In both cases, \(D\) is a normalization constant used so that abscissa values are closer to unity when they are sent to the fitter. This generally improves the probability that the fit will converge. This parameter may be specified via the div parameter. A value of 0 (the default) indicates that the application should determine a reasonable value for \(D\), which is determined via\(D = 10^{\int(\log10(\sqrt(\min(x)*\max(x)))}\)where min(x) and max(x) are the minimum and maximum abscissa values, respectively.The inputs are:
#spxfit :: Fit a 1-dimensional model to an image or image region
for determination of spectral index.
imagename = #Name of the input image(s)
box = '' #Rectangular box in
#direction coordinate blc, trc.
#Default: entire image ('').
region = '' #Region of interest. Default:
#Do not use a region.
chans = '' #Channels to use. Channels
#must be contiguous. Default: all channels ('').
stokes = '' #Stokes planes to
#use. Planes must be contiguous. Default:
#all stokes ('').
axis = -1 #The profile axis. Default:
#use the spectral axis if one
#exists, axis 0 otherwise (<0).
mask = '' #Mask to use. Default is none.
minpts = 1 #Minimum number of unmasked
#points necessary to attempt
#fit.
multifit = True #If true, fit a profile
#along the desired axis at each
#pixel in the specified
#region. If false, average the
#non-fit axis pixels and do
#a single fit to that average
#profile. Default False.
spxsol = '' #Name of the spectral index
#function coefficient solution
#image to write.
spxerr = '' #Name of the spectral index
#function coefficient error
#image to write.
model = '' #Name of model
#image. Default: do not write the model
#image ('').
residual = '' #Name of residual
#image. Default: do not write the
#residual image ('').
spxtype = 'plp' #Type of function to
#fit. 'plp' => power logarithmic
#polynomial, 'ltp' =>
#logarithmic transformed polynomial.
spxest = [] #Initial estimates for the
#spectral index function
#coefficients.
spxfix = [] #Fix the corresponding spectral index function
#coefficients during the fit. True=>hold fixed.
div = 0 #Divisor (numerical value or
#quantity) to use in the
#logarithmic terms of the
#plp or ltp function. 0 =>
#calculate a useful value on the fly.
wantreturn = True #Should a record summarizing
#the results be returned?
logresults = True #Output results to logger?
logfile = '' #File in which to log
#results. Default is not to write a
#logfile.
sigma = -1 #Standard deviation array or image name(s).
outsigma = '' #Name of output image used
#for standard deviation. Ignored
#if sigma is empty.
For more than a single input image or cube, all images must have the same dimensions along all axes other than the fit axis. multifit will perform a per-pixel fit, otherwise there will be a single value over the entire region.
Search for Spectral Line Rest Frequencies
The slsearch task allows the spectral line enthusiast to find their favorite spectral lines in subset of the Splatalogue spectral line catalog which is distributed with CASA. In addition, one can export custom catalogs from Splatalogue and import them to CASA using the task splattotable (next section) or tool method sl.splattotable. One can even import catalogs with lines not in Splatalogue using the same file format.The inputs to slsearch are as follows:
#slsearch :: Search a spectral line table.
tablename = '' #Input spectral line table name to
#search. If not specified, use the
#default table in the system.
outfile = '' #Results table name. Blank means do not
#write the table to disk.
freqrange = [84, 90] #Frequency range in GHz.
species = [''] #Species to search for.
reconly = False #List only NRAO recommended
#frequencies.
chemnames = [''] #Chemical names to search for.
qns = [''] #Resolved quantum numbers to search
#for.
rrlinclude = True #Include RRLs in the result set?
rrlonly = False #Include only RRLs in the result set?
intensity = -1 #CDMS/JPL intensity range. -1 -> do not
#use an intensity range.
smu2 = -1 #S*mu*mu range in Debye**2. -1 -> do
#not use an S*mu*mu range.
loga = -1 #log(A) (Einstein coefficient) range.
#-1 -> do not use a loga range.
eu = -1 #Upper energy state range in Kelvin. -1
#-> do not use an eu range.
el = -1 #Lower energy state range in Kelvin. -1
#-> do not use an el range.
verbose = True #List result set to logger (and
#optionally logfile)?
logfile = '' #List result set to this logfile (only
#used if verbose=True).
append = True #If true, append to logfile if it
#already exists, if false overwrite
#logfile if it exists. Only used if
#verbose=True and logfile not blank.
wantreturn = True #If true, return the spectralline tool
#associated with the result set.
The table is provided in the tablename parameter but if it is blank (the default), the catalog which is included with CASA will be used. Searches can be made in a parameter space with large dimensionality:
freqrange Frequency range in GHz.
species Species to search for.
reconly List only NRAO recommended frequencies.
chemnames Chemical names to search for.
qns Resolved quantum numbers to search for.
intensit CDMS/JPL intensity range.
smu2 $S\mu^{2}$ range in Debye$^2$.
loga log(A) (Einstein coefficient) range.
el Lower energy state range in Kelvin.
eu Upper energy state range in Kelvin.
rrlinclude Include RRLs in the result set?
rrlonly Include only RRLs in the result set?
Notation is as found in the Splatalogue catalog.
Example: Search for all lines of the species HOCN and HOCO\(^+\) in the 200-300GHz range:
slsearch(outfile='myresults.tbl', freqrange = [200,300],
species=['HOCN', 'HOCO+'])
The task can also return a python dictionary if assigned a variable like:
myLines = slsearch(outfile='myresults.tbl', freqrange = [200,300],
species=['HOCN', 'HOCO+'])
Convert Exported Splatalogue Catalogs to CASA Tables
In some cases, the internal spectral line catalog may not contain the lines in which one is interested. In that case, one can export a catalog from Splatalogue or even create their own ‘by hand’ (be careful to get the format exactly right though!). CASA’s task splattotable can then be used to create a CASA table that contains these lines and can be searched:
#splattotable :: Convert a downloaded Splatalogue spectral line list to a casa table.
filenames = [''] #Files containing Splatalogue lists.
table = '' #Output table name.
wantreturn = True #Do you want the task to return a spectralline tool attached to the results table?
A search in Splatalogue will return a catalog that can be saved in a file (look for the ‘Export’ section after the results on the search results page). The exported filename(s) should be entered in the filenames parameter of splattotable. The downloaded files must be in a specific format for this task to succeed. If you use the Splatalogue Export CASA fields feature, you should have no difficulties.
Bibliography¶
Leahy, J.P., Pooley, G.G., & Jagers, W.J. 1986, A&A, 156, 234 (http://adsabs.harvard.edu/abs/1986A%26A…156..234L)
Image Plane Analysis¶
Image-plane Component Fitting
Task imfit inputs are:
#imfit :: Fit one or more elliptical Gaussian components on an image region(s)
imagename = '' #Name of the input image
box = '' #Specify one or more box regions for the fit.
region = '' #Region.
chans = '' #Spectral channels on which to perform fit.
stokes = '' #Stokes parameter to fit. If blank, first stokes plane is
#used.
mask = '' #Mask to use. Default is none.
includepix = [] #Range of pixel values to include for fitting.
excludepix = [] #Range of pixel values to exclude for fitting.
residual = '' #Name of output residual image.
model = '' #Name of output model image.
estimates = '' #Name of file containing initial estimates of component
#parameters.
logfile = '' #Name of file to write fit results.
newestimates = '' #File to write fit results which can be used as initial
#estimates for next run.
complist = '' #Name of output component list table.
dooff = False #Also fit a zero level offset? Default is False
rms = -1 #RMS to use in calculation of uncertainties. Numeric or
#valid quantity (record or string). If numeric, it is
#given units of the input image. If quantity, units must
#conform to image units. If not positive, the rms of the
#residual image, in the region of the fit, is used.
noisefwhm = '' #Noise correlation beam FWHM. If numeric value,
#interpreted as pixel widths. If quantity (dictionary,
#string), it must have angular units.
imfit will return (as a Python dictionary) the results of the fit, but the results can also be written into a component list table or a logfile.
NOTE: To fit more than a single component, you must provide starting estimates for each component via the estimates file. See ‘’help imfit’’ for more details on this. A noise estimate will be calculated automatically or can be provided through the rms and noisefwhm keywords.
Examples for imfit
#First fit only a single component at a time
#This is OK since the components are well-separated and not blended
#Box around component A
xfit_A_res = imfit('b1608.demo.clean2.image',box='121,121,136,136',
newestimates='b1608.demo.clean2.newestimate')
#Now extract the fit part of the return value
xfit_A = xfit_A_res['results']['component0']
#xfit_A
#Out[7]:
#{'flux': {'error': array([ 6.73398035e-05, 0.00000000e+00, 0.00000000e+00,
#0.00000000e+00]),
#'polarisation': 'Stokes',
#'unit': 'Jy',
#'value': array([ 0.01753742, 0. , 0. , 0. ])},
#'label': '',
#'shape': {'direction': {'error': {'latitude': {'unit': 'arcsec',
#'value': 0.00041154866279462775},
#'longitude': {'unit': 'arcsec',
#'value': 0.00046695916589535109}},
#'m0': {'unit': 'rad', 'value': -2.0541102061078207}, NOTE: 'm0' and 'm1' are the coordinates of peak/controid
#'m1': {'unit': 'rad', 'value': 1.1439131060384089}, NOTE: 'm0' and 'm1' are the coordinates of peak/controid
#'refer': 'J2000',
#'type': 'direction'},
#'majoraxis': {'unit': 'arcsec', 'value': 0.29100166137741568},
#'majoraxiserror': {'unit': 'arcsec',
#'value': 0.0011186420613222663},
#'minoraxis': {'unit': 'arcsec', 'value': 0.24738110059830495},
#'minoraxiserror': {'unit': 'arcsec',
#'value': 0.0013431999725066338},
#'positionangle': {'unit': 'deg', 'value': 19.369249322401796},
#'positionangleerror': {'unit': 'rad',
#'value': 0.016663189295782171},
#'type': 'Gaussian'},
#'spectrum': {'frequency': {'m0': {'unit': 'GHz', 'value': 1.0},
#'refer': 'LSRK',
#'type': 'frequency'},
#'type': 'Constant'}}
#Now the other components
xfit_B_res = imfit('b1608.demo.clean2.image',box='108,114,120,126',
newestimates='b1608.demo.clean2.newestimate',append=True)
xfit_B = xfit_B_res['results']['component0']
xfit_C_res= imfit('b1608.demo.clean2.image',box='108,84,120,96')
xfit_C = xfit_C_res['results']['component0']
xfit_D_res = imfit('b1608.demo.clean2.image',box='144,98,157,110')
xfit_D = xfit_D_res['results']['component0']
print ""
print "Imfit Results:"
print "--------------"
print "A Flux = %6.4f Bmaj = %6.4f" % (xfit_A['flux']['value'][0],xfit_A['shape']['majoraxis']['value'])
print "B Flux = %6.4f Bmaj = %6.4f" % (xfit_B['flux']['value'][0],xfit_B['shape']['majoraxis']['value'])
print "C Flux = %6.4f Bmaj = %6.4f" % (xfit_C['flux']['value'][0],xfit_C['shape']['majoraxis']['value'])
print "D Flux = %6.4f Bmaj = %6.4f" % (xfit_D['flux']['value'][0],xfit_D['shape']['majoraxis']['value'])
print ""
Now try fitting four components together. For this we will have to provide an estimate file. We will use the clean beam for the estimate of the component sizes:
estfile=open('b1608.demo.clean2.estimate','w')
print >>estfile,'#peak, x, y, bmaj, bmin, bpa'
print >>estfile,'0.017, 128, 129, 0.293arcsec, 0.238arcsec, 21.7deg'
print >>estfile,'0.008, 113, 120, 0.293arcsec, 0.238arcsec, 21.7deg'
print >>estfile,'0.008, 113, 90, 0.293arcsec, 0.238arcsec, 21.7deg'
print >>estfile,'0.002, 151, 104, 0.293arcsec, 0.238arcsec, 21.7deg'
estfile.close()
Then, this can be used in imfit:
fit_all_res = imfit('b1608.demo.clean2.image',
estimates='b1608.demo.clean2.estimate',
logfile='b1608.demo.clean2.imfitall.log',
newestimates='b1608.demo.clean2.newestimate',
box='121,121,136,136,108,114,120,126,108,84,120,96,144,98,157,110')
#Now extract the fit part of the return values
xfit_allA = xfit_all_res['results']['component0']
xfit_allB = xfit_all_res['results']['component1']
xfit_allC = xfit_all_res['results']['component2']
xfit_allD = xfit_all_res['results']['component3']
These results are almost identical to those from the individual fits. You can see a nicer printout of the fit results in the logfile.
2-dimensional Smoothing; Image Convolution
A data cube can be smoothed across spatial dimensions with imsmooth. The inputs are:
#imsmooth :: Smooth an image or portion of an image
imagename = '' #Name of the input image. Must be
#specified.
kernel = 'gauss' #Type of kernel to use. Acceptable values
#are 'b', 'box', or 'boxcar' for a
#boxcar kernel, 'g', 'gauss', or
#'gaussian' for a gaussian kernel, 'c',
#'common', or 'commonbeam' to use the
#common beam of an image with multiple
#beams as the gaussian to which to
#convolve all the planes, 'i' or 'image'
#to use an image as the kernel.
beam = '' #Alternate way of describing a Gaussian.
#If specified, must be a dictionary with
#keys 'major', 'minor', and 'pa' (or
#'positionangle'). Do not specify beam
#if specifying major, minor, and pa.
#Example: Example: {'major': '5arcsec',
#'minor': '2arcsec', 'pa': '20deg'}.
targetres = False #If gaussian kernel, specified parameters
#are to be resolution of output image
#(True) or parameters of gaussian to
#convolve with input image (False).
major = '' #Major axis for the kernels. Standard
#quantity representation. Must be
#specified for kernel='boxcar'. Example:
#'4arcsec'.
minor = '' #Minor axis. Standard quantity
#representation. Must be specified for
#kernel='boxcar'. Example: '2arcsec'.
pa = '' #Position angle used only for gaussian
#kernel. Standard quantity
#representation. Example: '40deg'.
region = '' #Region selection. See Default is to use the full
#image.
box = '' #Rectangular region to select in
#direction plane. Default is to use the entire
#direction plane.
chans = '' #Channels to use. Default is to use all
#channels.
stokes = '' #Stokes planes to use. Default is to
#use all Stokes planes.
mask = '' #Mask to use. Default
#is none.
outfile = '' #Output image name. Must be specified.
overwrite = False #Overwrite (unprompted) pre-existing
#output file?
where the cube/image imagename will be convolved with a kernel defined in the kernel keyword. Kernels ‘gauss’ and ‘boxcar’ need the major and minor axes sizes as input, the Gaussian kernel smoothing also requires a position angle. By default, the kernel size defines the kernel itself, i.e. the data will be smoothed with this kernel. If the targetres parameter for Gaussian kernels is set to ‘True’, major and minor axes will be those from the output resolution, and the kernel will be adjusted for each plane to arrive at the final resolution. The ‘commonbeam’ kernel is to be used when the beam shape is different as a function of frequency. This option will then smooth all planes to a single beam, defined by the largest beam in the cube. With the ‘image’ kernel, one can specify an image that will serve as the convolution kernel. A scale factor can be applied, which defaults to flux conservation where units are Jy/beam or Jy/beam.km/s. For all other units, like K, the output will be scaled by the inverse of the convolution kernel. e.g., in the extreme case of a flat distribution the values before and after smoothing will be the same.
Examples:
Smoothing with a Gaussian kernel 20” by 10”
imsmooth( imagename='my.image', kernel='gauss', major='20arcsec', minor='10arcsec',targetres=T)
Smoothing using pixel coordinates and a boxcar kernel.
imsmooth( imagename='new.image', major='20pix', minor='10pix', kernel='boxcar')
Image Import/Export¶
These tasks will allow you to write your CASA image to a FITS file that other packages can read, and to import existing FITS files into CASA as an image. They can also extract data to Python data structures for further analysis.
FITS Image Export
To export your images to fits format use the exportfits task. The inputs are:
#exportfits :: Convert a CASA image to a FITS file
imagename = '' #Name of input CASA image
fitsimage = '' #Name of output image FITS file
velocity = False #Use velocity (rather than frequency) as spectral axis
optical = False #Use the optical (rather than radio) velocity convention
bitpix = -32 #Bits per pixel
minpix = 0 #Minimum pixel value
maxpix = 0 #Maximum pixel value
overwrite = False #Overwrite pre-existing imagename
dropstokes = False #Drop the Stokes axis?
stokeslast = True #Put Stokes axis last in header?
The dropstokes or stokeslast parameter may be needed to make the FITS image compatible with an external application.For example,
exportfits('ngc5921.demo.cleanimg.image','ngc5921.demo.cleanimg.image.fits')
FITS Image Import
You can also use the importfits task to import a FITS image into CASA image table format. Note, the CASA Viewer can read fits images so you don’t need to do this if you just want to look at the image. The inputs for importfits are:
#importfits :: Convert an image FITS file into a CASA image
fitsimage = '' #Name of input image FITS file
imagename = '' #Name of output CASA image
whichrep = 0 #If fits image has multiple
#coordinate reps, choose one.
whichhdu = 0 #If its file contains
#multiple images, choose one.
zeroblanks = True #Set blanked pixels to zero (not NaN)
overwrite = False #Overwrite pre-existing imagename
defaultaxes = False #Add the default 4D
#coordinate axes where they are missing
defaultaxesvalues = [] #List of values to assign to
#added degenerate axes when
#defaultaxes=True (ra,dec,freq,stokes)
For example, we can read the above image back in
importfits('ngc5921.demo.cleanimg.image.fits','ngc5921.demo.cleanimage')
Extracting data from an image
The imval task will extract the values of the data and mask from a specified region of an image and place in the task return value as a Python dictionary. The inputs are:
#imval :: Get the data value(s) and/or mask value in an image.
imagename = '' #Name of the input image
region = '' #Image Region. Use viewer
box = '' #Select one or more box regions
chans = '' #Select the channel(spectral) range
stokes = '' #Stokes params to image (I,IV,IQU,IQUV)
By default, box=’ ‘ will extract the image information at the reference pixel on the direction axes. Plane selection is controlled by chans and stokes. By default, chans=’ ‘ and stokes=’ ‘ will extract the image information in all channels and Stokes planes. For instance,
xval = imval('myimage', box='144,144', stokes='I' )
will extract the Stokes I value or spectrum at pixel 144,144, while
xval = imval('myimage', box='134,134.154,154', stokes='I' )
will extract a 21 by 21 pixel region. Extractions are returned in NumPy arrays in the return value dictionary, plus some extra elements describing the axes and selection:
CASA <2>: xval = imval('ngc5921.demo.moments.integrated')
CASA <3>: xval
Out[3]:
{'axes': [[0, 'Right Ascension'],
[1, 'Declination'],
[3, 'Frequency'],
[2, 'Stokes']],
'blc': [128, 128, 0, 0],
'data': array([ 0.89667124]),
'mask': array([ True], dtype=bool),
'trc': [128, 128, 0, 0],
'unit': 'Jy/beam.km/s'}
extracts the reference pixel value in this 1-plane image. Note that the ‘data’ and ‘mask’ elements are NumPy arrays, not Python lists. To extract a spectrum from a cube:
CASA <8>: xval = imval('ngc5921.demo.clean.image',box='125,125')
CASA <9>: xval
Out[9]:
{'axes': [[0, 'Right Ascension'],
[1, 'Declination'],
[3, 'Frequency'],
[2, 'Stokes']],
'blc': [125, 125, 0, 0],
'data': array([ 8.45717848e-04, 1.93370355e-03, 1.53750915e-03,
2.88399984e-03, 2.38683447e-03, 2.89159478e-04,
3.16268904e-03, 9.93389636e-03, 1.88773088e-02,
3.01138610e-02, 3.14478502e-02, 4.03211266e-02,
3.82498614e-02, 3.06552909e-02, 2.80734301e-02,
1.72479432e-02, 1.20884273e-02, 6.13593217e-03,
9.04005766e-03, 1.71429547e-03, 5.22095338e-03,
2.49114982e-03, 5.30831399e-04, 4.80734324e-03,
1.19265869e-05, 1.29435991e-03, 3.75700940e-04,
2.34788167e-03, 2.72604497e-03, 1.78467855e-03,
9.74952069e-04, 2.24676146e-03, 1.82263291e-04,
1.98463408e-06, 2.02975096e-03, 9.65532148e-04,
1.68218743e-03, 2.92119570e-03, 1.29359076e-03,
-5.11484570e-04, 1.54162932e-03, 4.68662125e-04,
-8.50282842e-04, -7.91683051e-05, 2.95954203e-04,
-1.30133145e-03]),
'mask': array([ True, True, True, True, True, True, True, True, True,
True, True, True, True, True, True, True, True, True,
True, True, True, True, True, True, True, True, True,
True, True, True, True, True, True, True, True, True,
True, True, True, True, True, True, True, True, True, True], dtype=bool),
'trc': [125, 125, 0, 45],
'unit': 'Jy/beam'}
To extract a region from the plane of a cube:
CASA <13>: xval = imval('ngc5921.demo.clean.image',box='126,128,130,129',chans='23')
CASA <14>: xval
Out[14]:
{'axes': [[0, 'Right Ascension'],
[1, 'Declination'],
[3, 'Frequency'],
[2, 'Stokes']],
'blc': [126, 128, 0, 23],
'data': array([[ 0.00938627, 0.01487772],
[ 0.00955847, 0.01688832],
[ 0.00696965, 0.01501907],
[ 0.00460964, 0.01220793],
[ 0.00358087, 0.00990202]]),
'mask': array([[ True, True],
[ True, True],
[ True, True],
[ True, True],
[ True, True]], dtype=bool),
'trc': [130, 129, 0, 23],
'unit': 'Jy/beam'}
CASA <15>: print xval['data'][0][1]
0.0148777160794
In this example, a rectangular box was extracted, and you can see the order in the array and how to address specific elements.
Math Operations / Statistics¶
Mathematical operations on images are typically completed using the CASA task immath, and image statistics may be derived using the CASA tasks imstat and imdev. Here, we give an overview of how these tasks are used.
Mathematical Operations on Images
The CASA task immath is useful for performing mathematical operations on images and on specific channels within images, including e.g. addition or subtraction of two cubes, squaring an image, computing a spectral index, and determining polarization angles and intensities. The inputs are:
#immath :: Perform math operations on images
imagename = '' #a list of input images
mode = 'evalexpr' #mode for math operation (evalexpr, spix, pola, poli)
expr = '' #Mathematical expression using images
varnames = '' #a list of variable names to use with the image files
outfile = 'immath_results.im' #File where the output is saved
mask = '' #Mask to use. Default is none.
region = '' #Region selection.
#Default is to use the full image.
box = '' #Rectangular region to
#select in direction plane.
#Default is to use the
#entire direction plane.
chans = '' #Channels to use.
#Default is to use all channels.
stokes = '' #Stokes planes to use.
#Default is to use all Stokes planes.
imagemd = '' #An image name from which metadata should be copied. The input
#can be either an image listed under imagename or any other
#image on disk. Leaving this parameter unset may copy header
#metadata from any of the input images, which
#one is not guaranteed.
Alert: immath does not convert any brightness units, e.g. from Jy/beam to K or vice versa. The user is responsible for making sure the images are consistent with the values in the header and the image. It is not advisable to mix input images that are in different units or have different beam sizes.
The imagename parameter must be given the name of a single image as a string (e.g. imagename=’image1.im’) or the names of multiple images in a list of strings (e.g. imagename=[‘image1.im’, ‘image2.im’] ). The immath task outputs an image file, and the name of the output file is specified using the outfile parameter.
The mode parameter selects what immath is to do. The default, mode=’evalexpr’, allows the user to specify a mathematical operation to execute on the input images through the expr sub-parameter. The mathematical expression is specified in expr as a Lattice Expression Language (LEL) string (see the page on LEL strings). The standard usage for mode=’evalexpr’ is to input a list of images into the imagename parameter, and then refer to them in the expr subparameter in LEL by the names IM0, IM1, …. For example,
immath(imagename=['image1.im','image2.im'],expr='IM0-IM1',outfile='ImageDiff.im')
would subtract the second image given from the first.
For the special modes ‘spix’, ‘pola’, and ‘poli’, the images required for the operation may need to be listed in imagename in a particular order. See examples of usage for polarization data below, paying particular attention to posted alerts.
The mathematical expression can be computed on the entire image cube, or on selected regions and image planes, which can be specified through the mask, region, box, chans, and stokes parameters. Mask specification is done using the mask parameter which can optionally contain an on-the-fly mask expression (in LEL) or point to an image with a pixel mask. In some cases, one may like to use a flat image (e.g. a moment image) mask applied to an entire cube. The stretch=True subparameter in mask allows one to expand the mask to all planes (i.e. channels or Stokes planes) of the cube. Region selection can also be carried out through the region parameter (see the pages on Region Files and Region File Format) and box parameter, while image plane selection is controlled by chans and stokes parameters.
The image metadata in the output file is adopted from another image, which can be specified through the imagemd parameter. In imagemd, input the name of the image from which the metadata should be copied and used for the output image. If left blank, the task may pick any of the input image headers, so it is better to define this parameter. In fact, the image specified in imagemd can be any image, even an image that is not part of the calculations in immath.
Detailed examples of immath usage are given below.
Examples for immath
In the following, we show a examples of immath usage. Note that the image names in the expr are assumed to refer to existing image files in the current working directory.
Simple math - Select a single plane (channel 22) of the 3-D cube:
immath(imagename='ngc5921.demo.cleanimg.image',
expr='IM0',chans='22',
outfile='ngc5921.demo.chan22.image')
Double all values in our image:
immath(imagename=['ngc5921.demo.chan22.image'],
expr='IM0*2.0',
outfile='ngc5921.demo.chan22double.image' )
Square all values in our image:
immath(imagename=['ngc5921.demo.chan22.image'],
expr='IM0^2',
outfile='ngc5921.demo.chan22squared.image' )
NOTE: The units in the output image are still claimed to be “Jy/beam”, i.e. immath will not correctly scale the units in the image for non-linear cases like this. Beware!
Subtract our image containing channel 22 from the original 3-D cube. Note that in this example, the 2-D plane (channel 22) is extended into the third dimension, so that che channel 22 image is subtracted from each plane in the 3-D cube:
immath(imagename=['ngc5921.demo.cleanimg.image','ngc5921.demo.chan22.image'],
expr='IM0-IM1',
outfile='ngc5921.demo.sub22.image')
Divide an image by another, with a threshold on one of the images:
immath(imagename=['ngc5921.demo.cleanimg.image','ngc5921.demo.chan22.image'],
expr='IM0/IM1[IM1>0.008]',
outfile='ngc5921.demo.div22.image')
You can do other mathematical operations on an image (e.g. trigonometric functions), as well as use scalar results from an image (e.g. max, min, median, mean, variance) in immath. You also have access to constants such as e() and pi(). As an example, the following expression uses the sine function, square root (sqrt), a scalar function (max), and the constant pi :
immath(imagename=['ngc5921.demo.chan22.image','ngc5921.demo.chan22squared.image'],
expr='sin(float(pi())*IM0/sqrt(max(IM1)))',
outfile='ngc5921.demo.chan22sine.image')
NOTE: Once again, the units in the output image are still claimed to be “Jy/beam”, i.e. immath will not correctly scale the units in the image for non-linear cases like this. Beware!
Region and Channel Selection - Select and save a region including the inner 1/4 of an image for channels 0 through 9 (chans=’<10’) and channels 40, 42, and 44:
default('immath')
imagename=['ngc5921.demo.cleanimg.image']
expr='IM0'
region='box[[64pix,64pix],[192pix,192pix]]'
chans='<10;40,42,44'
outfile='ngc5921.demo.inner.image'
immath()
If more than one channel is specified in the chans parameter, then the output image will contain multiple channels spanning the range from the lowest channel specified to the highest. In the example above, the output image will span channels 0 through 44, for a total of 45 channels. The channels that were not selected (in this case, channels 10 through 39 and channels 41 and 43) will be masked in the output cube. If we had set chans=’40,42,44’ then there would be 5 output channels corresponding to channels 40, 41, 42, 43, and 44 of the MS with 41 and 43 masked.
Note that the chans syntax allows the operators ‘<’, ‘<=’, ‘>’, and ‘>’. For example, the following two inputs select the same channels.
chans = '<17,>79'
chans = '<=16,>=80'
Polarization manipulation - Extract each Stokes plane from a cube into an individual image:
default('immath')
imagename = '3C129BC.clean.image'
outfile='3C129BC.I'; expr='IM0'; stokes='I'; immath();
outfile='3C129BC.Q'; expr='IM0'; stokes='Q'; immath();
outfile='3C129BC.U'; expr='IM0'; stokes='U'; immath();
outfile='3C129BC.V'; expr='IM0'; stokes='V'; immath();
Extract linearly polarized intensity and polarization position angle images:
immath(stokes='', outfile='3C129BC.P', mode='poli',
imagename=['3C129BC.Q','3C129BC.U'], sigma='0.0mJy/beam');
immath(stokes='', outfile='3C129BC.X', mode='pola',
imagename=['3C129BC.Q','3C129BC.U'], sigma='0.0mJy/beam');
ALERT: For mode=’pola’ you MUST call as a function as in this example (giving the parameters as arguments) or immath will fail.
Create a fractional linear polarization image:
default( 'immath')
imagename = ['3C129BC.I','3C129BC.Q','3C129BC.U']
outfile='3C129BC.fractional_linpol'
expr='sqrt((IM1^2 + IM2^2)/IM0^2)'
stokes=''
immath()
Create a polarized intensity image:
default( 'immath')
imagename = ['3C129BC.Q','3C129BC.U','3C129BC.V']
outfile='3C129BC.pol_intensity'
expr='sqrt(IM0^2 + IM1^2 + IM2^2)'
stokes=''
immath()
Toolkit Tricks: The following uses the toolkit. You can make a complex linear polarization (Q+iU) image using the imagepol tool:
#Make an imagepol tool and open the clean image
potool = casac.homefinder.find_home_by_name('imagepolHome')
po = potool.create()
po.open('3C129BC.clean.image')
#Use complexlinpol to make a Q+iU image
po.complexlinpol('3C129BC.cmplxlinpol')
po.close()
You can now display this in the Viewer, in particular overlay this over the intensity raster with the intensity contours. When you load the image, use the LEL:
'3C129BC.cmplxlinpol'['3C129BC.P'>0.0001]
which is entered into the LEL box at the bottom of the Load Data menu.
Using Masks in immath
The mask parameter is used inside immath to apply a mask to all the images used in expr before calculations are done (if you are curious, it uses the ia.subimage tool method to make virtual images that are then input in the LEL to the ia.imagecalc method).For example, let’s assume that we have made a single channel image using clean:
default('clean')
vis = 'ngc5921.demo.src.split.ms.contsub'
imagename = 'ngc5921.demo.chan22.cleanimg'
mode = 'channel'
nchan = 1
start = 22
step = 1
field = ''
spw = ''
imsize = [256,256]
cell = [15.,15.]
psfalg = 'clark'
gain = 0.1
niter = 6000
threshold='8.0mJy'
weighting = 'briggs'
rmode = 'norm'
robust = 0.5
mask = [108,108,148,148]
clean()
There is now a file ngc5921.demo.chan22.cleanimg.mask that is an image with values 1.0 inside the cleanbox region and 0.0 outside. We can use this to mask the clean image:
default('immath')
imagename = 'ngc5921.demo.chan22.cleanimg.image'
expr='IM0'
mask='"ngc5921.demo.chan22.cleanimg.mask">0.5'
outfile='ngc5921.demo.chan22.cleanimg.imasked'
immath()
Toolkit Tricks: Note that there are also pixel masks that can be contained in each image. These are Boolean masks, and are implicitly used in the calculation for each image in expr. If you want to use the mask in a different image not in expr, try it in mask:
#First make a pixel mask inside ngc5921.demo.chan22.cleanimg.mask
ia.open('ngc5921.demo.chan22.cleanimg.mask')
ia.calcmask('"ngc5921.demo.chan22.cleanimg.mask">0.5')
ia.summary()
ia.close()
#There is now a 'mask0' mask in this image as reported by the summary
#Now apply this pixel mask in immath
default('immath')
imagename='ngc5921.demo.chan22.cleanimg.image'
expr='IM0'
mask='mask(ngc5921.demo.chan22.cleanimg.mask)'
outfile='ngc5921.demo.chan22.cleanimg.imasked1'
immath()
Note that nominally the axes of the mask must be congruent to the axes of the images in expr. However, one exception is that the image in mask can have fewer axes (but not axes that exist but are of the wrong lengths). In this case, immath will extend the missing axes to cover the range in the images in expr. Thus, you can apply a mask made from a single channel to a whole cube.
#drop degenerate stokes and freq axes from mask image
ia.open('ngc5921.demo.chan22.cleanimg.mask')
im2 = ia.subimage(outfile='ngc5921.demo.chan22.cleanimg.mymask',dropdeg=True)
im2.summary()
im2.close()
ia.close()
#mymask has only RA and Dec axes
#Now apply this mask to the whole cube
default('immath')
imagename='ngc5921.demo.cleanimg.image'
expr='IM0'
mask='"ngc5921.demo.chan22.cleanimg.mymask">0.5'
outfile='ngc5921.demo.cleanimg.imasked'
immath()
Computing Image Statistics
The imstat task will calculate statistics on a region of an image and return the results as a value in a Python dictionary. The inputs are:
#imstat :: Displays statistical information from an image or image region
imagename = '' #Name of the input image.
axes = -1 #List of axes to evaluate statistics over. Default is
#all axes.
region = '' #Image Region or name. Use Viewer.
box = '' #Select one or more box regions.
chans = '' #Select the channel(spectral) range.
stokes = '' #Stokes params to image (I,IV,IQU,IQUV). Default '' =>
#include all
listit = True #Print stats and bounding box to logger?
verbose = False #Print additional messages to logger?
mask = '' #Mask to use. Default is none.
logfile = '' #Name of file to write fit results.
algorithm = 'classic' #Algorithm to use. Supported values are 'chauvenet',
#'classic', 'fit-half', and 'hinges-fences'. Minimum
#match is supported.
clmethod = 'auto' #Method to use for calculating classical statistics.
#Supported methods are 'auto', 'tiled', and
#'framework'. Ignored if algorithm is not 'classic'.
Area selection can be done using region and mask parameters. Plane selection is controlled by chans and stokes. The parameter axes will select the dimensions that the statistics are calculated over. Typical data cubes have axes like: RA axis 0, DEC axis 1, Velocity axis 2. So, e.g. axes=[0,1] would be the most common setting to calculate statistics per spectral channel.A typical output of imstat on a cube with axes=[0,1] and algorithm=’classic’ (default) looks like:
No region specified. Using full positional plane.
Using all spectral channels.
Using polarizations ALL
Determining stats for image IRC10216_HC3N.cube_r0.5.image
Set region from supplied region record
Statistics calculated using Classic algorithm
Regions ---
-- bottom-left corner (pixel) [blc]: [0, 0, 0, 0]
-- top-right corner (pixel) [trc]: [299, 299, 0, 63]
-- bottom-left corner (world) [blcf]: 09:48:01.492, +13.15.40.658, I, 3.63994e+10Hz
-- top-right corner (world) [trcf]: 09:47:53.299, +13.17.40.258, I, 3.63915e+10Hz
No region specified. Using full positional plane.
Using all spectral channels.
Using polarizations ALL
Selected bounding box :
[0, 0, 0, 0] to [299, 299, 0, 63] (09:48:01.492, +13.15.40.658, I, 3.63994e+10Hz to 09:47:53.299, +13.17.40.258, I, 3.63915e+10Hz)
#Frequency Frequency(Plane) Npts Sum Mean Rms Std dev Minimum Maximum
3.63993552e+10 0 9.000000e+04 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
3.63992302e+10 1 9.000000e+04 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
3.63991052e+10 2 9.000000e+04 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
3.63989802e+10 3 9.000000e+04 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
3.63988551e+10 4 9.000000e+04 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
3.63987301e+10 5 9.000000e+04 6.069948e-01 6.744386e-06 1.534640e-03 1.534634e-03 -6.355108e-03 6.166496e-03
3.63986051e+10 6 9.000000e+04 2.711720e-01 3.013023e-06 1.538071e-03 1.538077e-03 -6.165663e-03 5.862981e-03
3.63984801e+10 7 9.000000e+04 2.501259e-01 2.779177e-06 1.578049e-03 1.578056e-03 -6.771976e-03 6.272645e-03
3.63983551e+10 8 9.000000e+04 -3.706732e-01 -4.118591e-06 1.607191e-03 1.607194e-03 -8.871284e-03 6.591001e-03
where the header information provides the specifications of the data that were selected followed by the table with the frequency values of the planes, the plane numbers, Npts (the number of pixels per plane), and the Sum, Median, RMS, Standard deviations, Minimum, and Maximum of the pixel values for each plane. Similar output is provided when the data is averaged over different axes. The logger output can also be written into or appended to a log file for further processing elsewhere (logfile parameter).imstat has access to different statistics algorithms. Most of them represent different ways on how to treat distributions that are not Gaussian, in particular to eliminate outlier values from the statistics. Available algorithms are CLASSIC, where all unmasked pixels are used, FIT-HALF, where one (good) half of the distribution is being mirrored across a central value, HINGES-FENCES, where the inner quartiles plus a ‘fence’ data portion is being used, and CHAUVENET, which includes values based on the number of standard deviations from the mean. For more information, see the inline help of the imstat task.
Using the task return value
The contents of the return value of imstat are in a Python dictionary of key-value sets. For example,
xstat = imstat()
will assign this to the Python variable xstat. The keys for xstat are outlined on the imstat page.For example, an imstat call might be
default('imstat')
imagename = 'ngc5921.demo.cleanimg.image' #The NGC5921 image cube
box = '108,108,148,148' #20 pixels around the center
chans = '21' #channel 21
xstat = imstat()
In the terminal window, imstat reports:
Statistics on ngc5921.usecase.clean.image
Region ---
-- bottom-left corner (pixel) [blc]: [108, 108, 0, 21]
-- top-right corner (pixel) [trc]: [148, 148, 0, 21]
-- bottom-left corner (world) [blcf]: 15:22:20.076, +04.58.59.981, I, 1.41332e+09Hz
-- top-right corner( world) [trcf]: 15:21:39.919, +05.08.59.981, I, 1.41332e+09Hz
Values --
-- flux [flux]: 0.111799236126
-- number of points [npts]: 1681.0
-- maximum value [max]: 0.029451508075
-- minimum value [min]: -0.00612453464419
-- position of max value (pixel) [maxpos]: [124, 131, 0, 21]
-- position of min value (pixel) [minpos]: [142, 110, 0, 21]
-- position of max value (world) [maxposf]: 15:22:04.016, +05.04.44.999, I, 1.41332e+09Hz
-- position of min value (world) [minposf]: 15:21:45.947, +04.59.29.990, I, 1.41332e+09Hz
-- Sum of pixel values [sum]: 1.32267159822
-- Sum of squared pixel values [sumsq]: 0.0284534543692
Statistics ---
-- Mean of the pixel values [mean]: 0.000786836167885
-- Standard deviation of the Mean [sigma]: 0.00403944306904
-- Root mean square [rms]: 0.00411418313161
-- Median of the pixel values [median]: 0.000137259965413
-- Median of the deviations [medabsdevmed]: 0.00152346317191
-- Quartile [quartile]: 0.00305395200849
The return value in xstat is
CASA <152>: xstat
Out[152]:
{'blc': array([108, 108, 0, 21]),
'blcf': '15:22:20.076, +04.58.59.981, I, 1.41332e+09Hz',
'flux': array([ 0.11179924]),
'max': array([ 0.02945151]),
'maxpos': array([124, 131, 0, 21]),
'maxposf': '15:22:04.016, +05.04.44.999, I, 1.41332e+09Hz',
'mean': array([ 0.00078684]),
'medabsdevmed': array([ 0.00152346]),
'median': array([ 0.00013726]),
'min': array([-0.00612453]),
'minpos': array([142, 110, 0, 21]),
'minposf': '15:21:45.947, +04.59.29.990, I, 1.41332e+09Hz',
'npts': array([ 1681.]),
'quartile': array([ 0.00305395]),
'rms': array([ 0.00411418]),
'sigma': array([ 0.00403944]),
'sum': array([ 1.3226716]),
'sumsq': array([ 0.02845345]),
'trc': array([148, 148, 0, 21]),
'trcf': '15:21:39.919, +05.08.59.981, I, 1.41332e+09Hz'}
ALERT: The return dictionary currently includes NumPy array values, which have to be accessed by an array index to get the array value. To access these dictionary elements, use the standard Python dictionary syntax, e.g. xstat [][]
For example, to extract the standard deviation as a number
mystddev = xstat['sigma'][0]
print 'Sigma = '+str(xstat['sigma'][0])
Examples for imstat
To extract statistics for an image:
xstat = imstat('b1608.demo.clean2.image')
#Printing out some of these
print 'Max = '+str(xstat['max'][0])
print 'Sigma = '+str(xstat['sigma'][0])
#results:
#Max = 0.016796965152
#Sigma = 0.00033631979385
In a box around the brightest component:
xstat_A = imstat('b1608.demo.clean2.image',box='124,125,132,133')
#Printing out some of these
print 'Comp A Max Flux = '+str(xstat_A['max'][0])
print 'Comp A Max X,Y = ('+str(xstat_A['maxpos'][0])+','+str(xstat_A['maxpos'][1])+')'
#results:
#Comp A Max Flux = 0.016796965152
#Comp A Max X,Y = (128,129)
Computing a Deviation Image
The imdev task produces an output image whose value in each pixel represents the “error” or “deviation” in the input image at the corresponding pixel. The output image has the same dimensions and coordinate system as the input image, or as the selected region of the input image. The inputs are:
#imdev :: Create an image that can represent the statistical deviations of the input image.
imagename = '' #Input image name
outfile = '' #Output image file name. If left blank (the default), no
#image is written but a new image tool referencing
#the collapsed image is returned.
region = '' #Region selection. Default is to use the full image.
box = '' #Rectangular region(s) to select in direction plane.
#Default is to use the entire direction plane.
chans = '' #Channels to use. Default is to use all channels.
stokes = '' #Stokes planes to use. Default is to use all Stokes planes.
mask = '' #Mask to use. Default setting is none.
overwrite = False #Overwrite (unprompted) pre-existing output file? Ignored
#if "outfile" is left blank.
grid = [1, 1] #x,y grid spacing. Array of exactly two positive integers.
anchor = 'ref' #x,y anchor pixel location. Either "ref" to use the image
#exactly two integers.
xlength = '1pix' #Either x coordinate length of box, or diameter of circle.
#Circle is used if ylength is a reference pixel or an
#empty string.
ylength = '1pix' #y coordinate length of box. Use a circle if ylength is
#an empty string.
interp = 'cubic' #Interpolation algorithm to use. Can be "nearest", "linear",
#or "cubic". Minimum match supported.
stattype = 'sigma' #Statistic to compute. See below for the list of supported
#statistics.
statalg = 'classic' #Statistics computation algorithm to use. Supported values
#are "chauvenet" and "classic". Minimum match is supported.
Area selection can be done using the region and mask parameters. Plane selection is controlled by the chans and stokes parameters. Statistics are computed spatially: a deviation image is computed independently for each channel/Stokes plane. If the outfile parameter is left blank, the task returns an image tool referencing the resulting image; otherwise the resulting image is written to disk.
The statistic to be computed is selected using the stattype parameter. Allowed statistics are:
iqr inner quartile range (q3 - q1)
max maximum
mean mean
medabsdevmed, madm median absolute deviation from the median
median median
min minimum
npts number of points
q1 first quartile
q3 third quartile
rms rms
sigma, std standard deviation
sumsq sum of squares
sum sum
var variance
xmadm median absolute deviation from the median multipied by x, where x is the
reciprocal of Phi^-1(3/4),where Phi^-1 is the reciprocal of the quantile
function. Numerically, x = 1.482602218505602. See, eg,
https://en.wikipedia.org/wiki/Median_absolute_deviation#Relation_to_standard_deviation
The chosen statistic is calculated around a set of grid points (pixels) across the input image with grid spacing specified by the grid parameter. The size and shape of the region used to compute the statistic at each grid point is specified by the xlength and ylength parameters. If ylength is an empty string, then the region used is a circle centered on each grid point with diameter provided by xlength. Otherwise, a rectangular region with dimensions of xlength by ylength is used. These two parameters may be specified as valid quantities with recognized units (e.g., “4arcsec” or “4pix”). They may also be specified as numerical values, in which case the unit is assumed to be pixels.
The chosen statistic is calculated at every grid point in the input image, and the result is reflected at the corresponding pixel of the output image. Values at all other pixels in the output image are determined by interpolating across the grid points using the interpolation scheme given by the input parameter interp. The statalg parameter specifies the algorithm for the statistics computation. Available algorithms are CLASSIC, where all unmasked pixels are used, and CHAUVENET, which includes values based on the number of standard deviations from the mean.
Examples for imdev
Compute a “standard deviation” image using grid-spacing of 4 arcsec in the X direction and 5 arcsec in the Y direction, with linear interpolation to compute values at non-grid-pixels. Compute the standard deviation in a box of 20 x 25 arcsec.
imdev("my.image", "std.image", grid=[4,5], xlength="20arcsec", ylength="25arcsec",
stattype="sigma", interp="linear", statalg="classic")
Compute an image showing median absolute deviation (MAD) across the image, with MAD converted to an equivalent RMS value. Anchor the grid at a specific pixel [1000,1000] with grid-spacing of 10 pixels, and use circles of diameter 30 pixels for the statistical computation. Calculate the statistic using the z-score/Chauvenet algorithm by fixing the maximum z-score to determine outliers to 5. Use cubic interpolation to determine the value at non-grid-pixels. Have the task return a pointer to the output image.
myim = imdev("my.image", anchor=[1000,1000], grid=[10,10], xlength=30, ylength='',
stattype="xmadm", interp="cubic", statalg="chauvenet", zscore=5)
Image Selection Parameters¶
Many of the image analysis tasks contain a set of parameters that can be used to define a sub-section of an image that the task works on. This includes selection in the spatial coordinates, typically RA/DEC or GLON/GLAT (or image axes 0,1) which can be defined by the box parameter. Spectral selection can be achieved by the chans parameter (typically axis 3) and Stokes selection with the stokes parameter (typically axis 2). These parameters are described below and are a quick way to select sub-images. For more complicated selections, we recommend the usage of the CASA region file (CRTF file).
Region Selection
Direction (e.g. RA, Dec) areal selection in the image analysis tasks is controlled by the box parameter or through the region parameter. Note that one should either specify a region (recommended) or any of box/chans/stokes parameters. Specifying both at the same time will give priority to the command line inputs in ‘chans’ and ‘stokes’ but will keep the region file specification for the spatial region selection.
The box parameter selects spatial rectangular areas:
box = '' #Select one or more box regions
E.g. for right ascension and declination. Boxes are specified by their bottom-left corner (blc) and top-right corner (trc) as follows: blcx, blcy, trcx, trcy. At the moment, CASA only accepts pixel values. The default (empty string) is to select all pixels of an image.
Example:
box='0,0,50,50'
selects a square with 50 pixels on the side starting at 0.
box='10,20,30,40,100,100,150,150'
selects two regions at a time. The first and the last four values mark the two boxes, following (blcx1, blcy1, trcx1, trcy1, blcx2, blcy2, trcx2, trcy2), with b/t = bottom/top and l/r = left/right.
ALERT: Note that one should either specify a region (recommended) or any of box/chans/stokes. If both are specified, the following rules apply: - If the region parameter is specified as a python dictionary (e.g. such as various rg tool methods return), a binary region file, or a region-in-image, it is not permissible to specify any of the box, chans, or stokes parameters. - If the region parameter is specified to be a CRTF file name, or a CRTF region string, then the resulting region and box selection is the union of the box specification with any regions in the CRTF file/string. This is the equivalent of translating the box specification into the equivalent “box” CRTF specification and prepending that specification to the specified CRTF file/string in the region parameter.
Plane Selection (chans, stokes)
The channel, frequency, or velocity plane(s) of the image is chosen using the chans parameter:
chans = '' #Select the channel(spectral) range
Using channel numbers, it is possible to also set ranges, e.g.**
chans='0,3,4,8' #select channels 0, 3, 4, 8
chans='3~20;50,51' #channels 3 to 20 and 50 and 51
chans='<10;>=55' #channels 0 to 9 and 55 and greater (inclusively)*
*
Sometimes, as in the immoments task, the channel/plane selection is generalized to work on more than one axis type. In this case, the planes parameter is used. This behaves like chans in syntax.
chans can also be set in the CASA region format to allow settings in frequency and velocity, e.g.
chans=('range=[-50km/s,50km/s], restfreq=100GHz, frame=LSRK')
This example would even define a new velocity system independent of the one in the image itself. If the rest frequency and velocity frame within the image are being used, the latter two entries are not needed. The parentheses are needed when the call is in a single command.
A frequency selection looks like:
chans=('range=[100GHz,100.125GHz]')
The polarization plane(s) of the image is chosen with the stokes parameter:
stokes = '' #Stokes params to image (I,IV,IQU,IQUV)
stokes parameters to image. Best practice is to separate the stokes parameters by commas, but this is not essential if no ambiguity exists. Options are: ‘I’,’Q’,’U’,’V’, ‘I’,’IV’,’QU’,’IQUV’, ‘RR’,’RL’,’LR’,’LL’, ‘XX’,’YX’,’XY’,’YY’,…
Examples:
stokes='IQUV';
stokes='I,Q'
stokes='RR,RL'
ALERT: For image analysis tasks and tool methods which also accept the region parameter, the following rules apply if both the chans and region parameters are simultaneously specified: - If the region parameter is specified as a python dictionary (e.g. such as various rg tool methods return), a binary region file, or a region-in-image, it is not permissable to specify any of the box, chans, or stokes parameters. - If the region parameter is specified to be a CRTF file name, or a CRTF region string, it is only permissable to specify chans if that specification can be represented as a single contiguous channel range. In that case, the chans specification overrides any global or per-region range specification in the CRTF file/string, and is used as the global spectral range selection for all regions in the CRTF file/string.
Region tool and multi-dimensional selection
The region parameter only supports 2D representations of shapes (box, circle, etc) in general astronomical coordinate systems, such as RA and dec or pixel coordinates. Additional axes need to be specified with different parameters, such as chans or stokes.
The region tool is less constrained, as it allows multi-dimensional selection in pixel coordinates, independent of the coordinate axes. For example:
reg = rg.box([0,0,0,0],[5,6,7,8])
The return value is a big dictionary specifying everything the software needs to know about the region. The reg variable can then be specified as the value of the region parameter in ia tool methods.
Region Files¶
The region parameter points to a CASA region which can be directly specified or listed in a ImageRegion file. An ImageRegion file can be created with the CASA viewer’s region manager, or directly using the CASA region file (crtf) syntax. Typically ImageRegion files will have the suffix “.crtf” for CASA Region Text Format.
For example:
region='circle[[18h12m24s, -23d11m00s], 2.3arcsec]'
or
region='myimage.im.crtf'
to specify a region file. For the most part, the region parameter in tasks only accepts strings (e.g. file names, region shape descriptions) while the region parameter in ia tool methods only accepts python region dictionaries (e.g. produced using the rg tool).
Alert: When both the region parameter and any of box/chans/stokes are specified simultaneously, the task may perform unwanted selections. While it is safest to only specify one of these (sets of) parameters, the following rules apply:
If region is specified as a python dictionary (eg such as various rg tool methods return), a binary region file, or a region-in-image, then it is not permissable to specify any of the box, chans, or stokes parameters.
If the region parameter is specified to be a CRTF file name, or a CRTF region string, the following rules apply: - If box is specified, the resulting selection is the union of the box specification with any regions in the CRTF file/string. This is the equivalent of translating the box specification into the equivalent “box” CRTF specification and prepending that specification to the specified CRTF file/string in the region parameter.
If chans is specified, it must be able to be represented as a single contiguous range of channels. In this case, the chans specification overrides any global or per-region range specification in the CRTF file/string, and is used as the global spectral range selection for all regions in the CRTF file/string.
If stokes is specified, this specification overrides any global or per-region corr specification in the CRTF file/string, and is used as the global correlation selection for all regions in the CRTF file/string.
NOTE: The CASA image analysis tasks will determine how a region is projected on a pixel image. The current CASA definition is that when the center of a pixel is inside the region, the full pixel is considered to be included in the region. If the center of the pixel is outside the region, the full pixel will be excluded. Note that the CASA viewer behavior is not entirely consistent and for rectangles it assumes that any fractional pixel coverage will include the entire pixel. For other supported shapes (ellipses and polygons), however, the Viewer adheres to the ‘center of pixel’ definition, consistent with the image analysis tools and tasks.
For purely single-pixel work regions may not necessarily be the best choice and alternate methods may be preferable to using regions, eg. ia.topixel, ia.toworld, ia.pixelvalue.
ALERT: Some region file specifications are not recognized by the viewer, the viewer only supports rectangles (box), ellipses, and polygons.
Region File Format¶
The CASA region file format provides a flexible, easily edited set of region definitions which are accepted across CASA tasks. Region files may be written by hand or using the CASA Viewer.
ALERT: Whereas the region format is supported by all the data processing tasks, the viewer implementation is still limited to rectangles, ellipses, and some markers.
For a file to be recognized as a valid CASA region text file, the first line must contain the string:
#CRTF
“CRTF” stands for “CASA Region Text Format”. One may also include an optional version number at the end of the string, so it reads #CRTFv0; this indicates the version of the format definition.
Region files have two different kinds of definitions, “regions” and “annotations”, each of which is one line long. To indicate an annotation, a line must begin with “ann”. Lines that begin with the comment character (#) are not considered for processing or display.
The second line of a file may define global parameters that are to be used for all regions and annotations in that file, in which case the line starts with the word “global”. The parameters set here may also be overridden by keywords in a specific line, in which case the keywords pertain only to that one line.
NOTE: All regions are considered by tasks. They will be displayed by visualization tasks as well as used to create masks, etc., as appropriate. Annotations are used by display tasks, and are for visual reference only.
Some tasks, like clean, require that a region cannot be entirely outside the image.
Region Definitions
All regions lines will follow this general arrangement:
{shape} {additional parameter=value pairs}
The possible parameter/value pairs are described in more detail below. Note that most parameters beyond the shape and its coordinates can be defined globally.
Possible units for coordinates are:
sexagesimal, e.g. 18h12m24s for right ascension or -03.47.27.1 for declination
decimal degrees, e.g. 140.0342deg for both RA and Dec
radians, e.g. 2.37666rad for both RA and Dec
pixels, e.g. 204pix
Possible units of length are:
degrees, e.g. 23deg
arcminutes, e.g. 23arcmin
arcseconds, e.g. 23arcsec
radians, e.g. 0.00035rad
pixels, e.g. 23pix
Units must always be included when defining a region.
NOTE: The CASA image analysis tasks will determine how a region is projected on a pixel image. The current CASA definition is that when the center of a pixel is inside the region, the full pixel is considered to be included in the region. If the center of the pixel is outside the region, the full pixel will be excluded. Note that the CASA viewer behavior is not entirely consistent and for rectangles it assumes that any fractional pixel coverage will include the entire pixel. For other supported shapes (ellipses and polygons), however, ithe viewer adheres to the ‘center of pixel’ definition, consistent with the image analysis tools and tasks.
For purely single-pixel work regions may not necessarily be the best choice and alternate methods may be preferable to using regions, eg. ia.topixel, ia.toworld, ia.pixelvalue.
Allowed Shapes
Rectangular box; the two coordinates are two opposite corners:
box[[x1, y1], [x2, y2]]
Center box; [x, y] define the center point of the box and [x_width, y_width] the width of the sides:
centerbox[[x, y], [x_width, y_width]]
Rotated box; [x, y] define the center point of the box; [x_width, y_width] the width of the sides; rotang the rotation angle:
rotbox[[x, y], [x_width, y_width], rotang]
Polygon; there could be many [x, y] corners. If parts of the polygon overlap, then the pixels in that overlapping region are taken into account. Note that the last point will connect with the first point to close the polygon:
poly[[x1, y1], [x2, y2], [x3, y3], ...]
Circle; center of the circle [x,y], r is the radius:
circle[[x, y], r]
Annulus; center of the circle is [x, y], [r1, r2] are inner and outer radii:
annulus[[x, y], [r1, r2]]
Ellipse; center of the ellipse is [x, y]; semi-major and semi-minor axes are [bmaj, bmin]; position angle of the major axis is pa:
ellipse[[x, y], [bmaj, bmin], pa]
Annotation Definitions
In addition to the definitions for regions, above, the following are always treated as annotations:
Line; coordinates define the end points of the line:
line[[x1, y1], [x2, y2]]
Vector; coordinates define end points; second coordinate pair is location of tip of arrow:
vector[[x1, y1], [x2, y2]]
Text; coordinates define leftmost point of text string:
text[[x, y], ’my text’]
Symbol; coordinates define location of symbol:
symbol[[x, y], {symbol}]
Global Definitions
Definitions to be used throughout the region file are placed on a line beginning with ‘global’, usually at the top of the file. These definitions may also be used on any individual region or annotation line; in this case, the value defined on that line will override the predefined global (but only for that line). If a ‘global’ line occurs later in the file, subsequent lines will obey those definitions.
Coordinate reference frame:
Supported values: J2000, B1950, B1950_VLA, BMEAN, GALACTIC, ECLIPTIC, SUPERGAL, ICRS
Default: image value
coord = J2000
Frequency/velocity axis:
Possible values: REST, LSRK, LSRD, BARY, GEO, TOPO, GALACTO, LGROUP, CMB
Default: image value
frame=TOPO
Frequency/velocity range:
Possible units: GHz, MHz, kHz, km/s, Hz, channel, chan (=channel)
Default: image range
range=[min, max]
Correlation axis:
Possible values: I, Q, U, V, RR, RL, LR, LL, XX, XY, YX, YY, RX, RY, LX, LY, XR, XL, YR, YL, PP, PQ, QP, QQ, RCircular, LCircular, Linear, Ptotal, Plinear, PFtotal, PFlinear, Pangle
Default: all planes present in image
corr=[X, Y]
Velocity calculation:
Possible values: RADIO, OPTICAL, Z, BETA, GAMMA
Default: image value
veltype=RADIO
Rest frequency:
Default: image value
restfreq=1.42GHz
Line characteristics:
Possible values: any line style recognized by matplotlib: ‘-‘=solid, ‘–’=dashed, ‘:’=dotted
Default:
linewidth=1
linestyle=’-’
Symbol characteristics:
Symbol size and thickness:
symsize = 1
symthick = 1
Region, symbol, and text color:
Possible values: any color recognized by matplotlib, including hex values
Default:
color=green
color=red
Text font characteristics:
Possible values: see below
‘usetex’ is a boolean parameter that determines whether or not the text line should be interpreted as LaTeX, and would require working LaTeX, dvipng, and Ghostscript installations (equivalent to the text.usetex parameter in matplotlib).
font=Helvetica
fontsize=10pt
fontstyle=bold
usetex=True/False
Label position:
Possible values: ‘left’, ‘right’, ‘top’, ‘bottom’
Default: ‘top’
labelpos=’right’
Label color:
Default: color of associated region.
Allowed values: same as values for region colors.
labelcolor=’green’
Label offset:
Default: [0,0].
Allowed values: any positive or negative number, in units of pixels.
labeloff=[1, 1]
Allowed Additional Parameters
These must be defined per region line:
Labels: text label for a region; should be placed so text does not overlap with region boundary
label=’string’
“OR/NOT” operators: A “+” at the beginning of a line will flag it with a boolean “OR” (default), and a “-” will flag it with a boolean “NOT”. Overlapping regions will be treated according to their sequence in the file; i.e., ((((entireImage OR line1) OR line2) NOT line3) OR line4). This allows some flexibility in building “non-standard” regions. Note that a task (e.g., clean) will still consider all lines: if one wishes to remove a region from consideration, it should be commented out (“#”).
Default: OR (+)
Examples
A file with both global definitions and per-line definitions:
#CRTFv0
global coord=B1950_VLA, frame=BARY, corr=[I, Q], color=blue
# A simple circle region:
circle[[18h12m24s, -23d11m00s], 2.3arcsec]
# A box region, this one only for annotation:
ann box[[140.0342deg, -12.34243deg], [140.0360deg, -12.34320deg]]
# A rotated box region, for a particular range of velocities:
rotbox[[12h01m34.1s, 12d23m33s], [3arcmin, 1arcmin], 12deg], range=[-1240km/s, 1240km/s]
# An annular region, overriding some of the global defaults:
annulus[[17h51m03.2s, -45d17m50s], [4.12deg, 0.10deg]], corr=[I,Q,U,V], color=red, label=’My label here’
# Cuts an ellipse out of the previous regions, but only for Q and a particular frequency range:
-ellipse[[17:51:03.2, -45.17.50], [1.34deg, 0.25deg], 45rad], range=[1.420GHz, 1.421GHz], corr=[Q], color=green, label=’Removed this’
# A diamond marker, in J2000 coordinates:
symbol[[32.1423deg, 12.1412deg], D], linewidth=2, coord=J2000, symsize=2
Fonts and Symbols
Allowed Symbols
symbol |
description |
|---|---|
. |
point marker |
, |
pixel marker |
o |
circle marker |
v |
triangle_down marker |
^ |
triangle_up marker |
< |
triangle_left marker |
> |
triangle_right marker |
1 |
tri_down marker |
2 |
tri_up marker |
3 |
tri_left marker |
4 |
tri_right marker |
s |
square marker |
p |
pentagon marker |
* |
star marker |
h |
hexagon1 marker |
H |
hexagon2 marker |
plus marker |
|
x |
x marker |
D |
diamond marker |
d |
thin_diamond marker |
| |
vline marker |
_ |
hline marker |
Allowed Fonts for Linux
“Century Schoolbook L”, “Console”, “Courier”, “Courier 10 Pitch”, “Cursor”, “David CLM”, “DejaVu LGC Sans”, “DejaVu LGC Sans Condensed”, “DejaVu LGC Sans Light”, “DejaVu LGC Sans Mono”, “DejaVu LGC Serif”, “DejaVu LGC Serif Condensed”, “Dingbats”, “Drugulin CLM”, “East Syriac Adiabene”, “Ellinia CLM”, “Estrangelo Antioch”, “Estrangelo Edessa”, “Estrangelo Nisibin”, “Estrangelo Nisibin Outline”, “Estrangelo Talada”, “Fangsong ti”, “Fixed [Sony]”, “Fixed [Eten]”, “Fixed [Misc]”, “Fixed [MNKANAME]”, “Frank Ruehl CLM”, “fxd”, “Goha-Tibeb Zemen”, “goth_p”, “Gothic [Shinonome]”, “Gothic [mplus]”, “hlv”, “hlvw”, “KacstArt”, “KacstBook”, “KacstDecorative”, “KacstDigital”, “KacstFarsi”, “KacstLetter”, “KacstPoster”, “KacstQura”, “KacstQuraFixed”, “KacstQuran”, “KacstTitle”, “KacstTitleL”, “Liberation Mono”, “Liberation Sans”, “Liberation Serif”, “LKLUG”, “Lohit Bengali”, “Lohit Gujarati”, “Lohit Hindi”, “Lohit Kannada”, “Lohit Malayalam”, “Lohit Oriya”, “Lohit Punjabi”, “Lohit Tamil”, “Lohit Telugu”, “LucidaTypewriter”, “Luxi Mono”, “Luxi Sans”, “Luxi Serif”, “Marumoji”, “Miriam CLM”, “Miriam Mono CLM”, “MiscFixed”, “Monospace”, “Nachlieli CLM”, “Nimbus Mono L”, “Nimbus Roman No9 L”, “Nimbus Sans L”, “Nimbus Sans L Condensed”, “PakTypeNaqsh”, “PakTypeTehreer”, “qub”, “Sans Serif”, “Sazanami Gothic”, “Sazanami Mincho”, “Serif”, “Serto Batnan”, “Serto Jerusalem”, “Serto Jerusalem Outline”, “Serto Mardin”, “Standard Symbols L”, “sys”, “URW Bookman L”, “URW Chancery L”, “URW Gothic L”, “URW Palladio L”, “Utopia”, “Yehuda CLM”
Allowed Fonts for MacOS X
“Abadi MT Condensed Light”, “Adobe Caslon Pro”, “Adobe Garamond Pro”, “Al Bayan”, “American Typewriter”, “Andale Mono”, “Apple Braille”, “Apple Chancery”, “Apple LiGothic”, “Apple LiSung”, “Apple Symbols”, “AppleGothic”, “AppleMyungjo”, “Arial”, “Arial Black”, “Arial Hebrew”, “Arial Narrow”, “Arial Rounded MT Bold”, “Arial Unicode MS”, “Arno Pro”, “Ayuthaya”, “Baghdad”, “Baskerville”, “Baskerville Old Face”, “Batang”, “Bauhaus 93”, “Bell Gothic Std”, “Bell MT”, “Bernard MT Condensed”, “BiauKai”, “Bickham Script Pro”, “Big Caslon”, “Birch Std”, “Blackoak Std”, “Book Antiqua”, “Bookman Old Style”, “Bookshelf Symbol 7”, “Braggadocio”, “Britannic Bold”, “Brush Script MT”, “Brush Script Std”, “Calibri”, “Calisto MT”, “Cambria”, “Candara”, “Century”, “Century Gothic”, “Century Schoolbook”, “Chalkboard”, “Chalkduster”, “Chaparral Pro”, “Charcoal CY”, “Charlemagne Std”, “Cochin”, “Colonna MT”, “Comic Sans MS”, “Consolas”, “Constantia”, “Cooper Black”, “Cooper Std”, “Copperplate”, “Copperplate Gothic Bold”, “Copperplate Gothic Light”, “Corbel”, “Corsiva Hebrew”, “Courier”, “Courier New”, “Curlz MT”, “DecoType Naskh”, “Desdemona”, “Devanagari MT”, “Didot”, “Eccentric Std”, “Edwardian Script ITC”, “Engravers MT”, “Euphemia UCAS”, “Eurostile”, “Footlight MT Light”, “Franklin Gothic Book”, “Franklin Gothic Medium”, “Futura”, “Garamond”, “Garamond Premier Pro”, “GB18030 Bitmap”, “Geeza Pro”, “Geneva”, “Geneva CY”, “Georgia”, “Giddyup Std”, “Gill Sans”, “Gill Sans MT”, “Gill Sans Ultra Bold”, “Gloucester MT Extra Condensed”, “Goudy Old Style”, “Gujarati MT”, “Gulim”, “GungSeo”, “Gurmukhi MT”, “Haettenschweiler”, “Harrington”, “HeadLineA”, “Hei”, “Heiti SC”, “Heiti TC”, “Helvetica”, “Helvetica CY”, “Helvetica Neue”, “Herculanum” “Hiragino Kaku Gothic Pro”, “Hiragino Kaku Gothic ProN”, “Hiragino Kaku Gothic Std”, “Hiragino Kaku Gothic StdN”, “Hiragino Maru Gothic Pro”, “Hiragino Maru Gothic ProN”, “Hiragino Mincho Pro”, “Hiragino Mincho ProN”, “Hiragino Sans GB”, “Hobo Std”, “Hoefler Text”, “Impact”, “Imprint MT Shadow”, “InaiMathi”, “Kai”, “Kailasa”, “Kino MT”, “Kokonor”, “Kozuka Gothic Pro”, “Kozuka Mincho Pro”, “Krungthep”, “KufiStandardGK”, “Letter Gothic Std”, “LiHei Pro”, “LiSong Pro”, “Lithos Pro”, “Lucida Blackletter”, “Lucida Bright”, “Lucida Calligraphy”, “Lucida Console”, “Lucida Fax”, “Lucida Grande”, “Lucida Handwriting”, “Lucida Sans”, “Lucida Sans Typewriter”, “Lucida Sans Unicode”, “Marker Felt”, “Marlett”, “Matura MT Script Capitals”, “Meiryo”, “Menlo”, “Mesquite Std”, “Microsoft Sans Serif”, “Minion Pro”, “Mistral”, “Modern No. 20”, “Monaco”, “Monotype Corsiva”, “Monotype Sorts”, “MS Gothic”, “MS Mincho”, “MS PGothic”, “MS PMincho”, “MS Reference Sans Serif”, “MS Reference Specialty”, “Mshtakan”, “MT Extra”, “Myriad Pro”, “Nadeem”, “New Peninim MT”, “News Gothic MT”, “Nueva Std”, “OCR A Std”, “Onyx”, “Optima”, “Orator Std”, “Osaka”, “Papyrus”, “PCMyungjo”, “Perpetua”, “Perpetua Titling MT”, “PilGi”, “Plantagenet Cherokee”, “Playbill”, “PMingLiU”, “Poplar Std”, “Prestige Elite Std”, “Raanana”, “Rockwell”, “Rockwell Extra Bold”, “Rosewood Std”, “Sathu”, “Silom”, “SimSun”, “Skia”, “Stencil”, “Stencil Std”, “STFangsong”, “STHeiti”, “STKaiti”, “STSong”, “Symbol”, “Tahoma”, “Tekton Pro”, “Thonburi”, “Times”, “Times New Roman”, “Trajan Pro”, “Trebuchet MS”, “Tw Cen MT”, “Verdana”, “Webdings”, “Wide Latin”, “Wingdings”, “Wingdings 2”, “Wingdings 3”, “Zapf Dingbats”, “Zapfino”
Image Masks¶
A mask can be used to define whether part of an image is used or not. There are different options for masks:
an image cube with Boolean True (not masked) or False (masked) values: They usually live inside image cubes and are automatically applied to the data. More than one mask may exist in a cube. The task makemask can be used to access and manipulate internal Boolean masks via the image.im:mask syntax.
an image cube with zero (masked) and non-zero (not masked) values: They are their own image cubes and are applied to other image cubes when needed.
an LEL string for a condition.
Masks (mask parameter)
Using image cubes is useful to mask on a pixel by pixel basis, where False and zeros mark masked (excluded) pixels. Both versions can be converted into each other with the task makemask. Some analysis tasks show an optional stretch parameter which is useful, e.g., to expand a single plane mask to an entire cube along the spectral axis.
To use a different zero/non-zero mask (in this case ‘myimage.mask’), the parameter can be set like the following:
mask='mask(myimage.mask)'
The default boolean masks inside images will also be respected with the above syntax.
But remember that an image can have multiple Boolean masks, so to use the mask2 in an image, set the parameter as:
mask='mask(myimage.mask:mask2)'
using the syntax where the mask is specified after the colon. To see what masks are present in your image, use the makemask task.
An LEL string can be an on-the-fly (OTF) mask expression or refer to an image pixel mask.
Note that the mask file supplied in the mask parameter must have the same shape, same number of axes, and same axes length as the images supplied in the expr parameter, with one exception. The mask may be missing some of the axes, if this is the case then the mask will be expanded along these axes to become the same shape.
The following example uses the mask from file ngc5921.clean.cleanbox.mask :
mask='mask(ngc5921.clean.cleanbox.mask)'
Here, the mask is calculated to be all pixels with values larger than 0.5:
mask='"ngc5921.clean.cleanbox.mask">0.5'
Because it is an LEL expression, care must be taken to properly escape characters which LEL views as special. For details, see the LEL document. As an example, specifying
mask = "3clean_mask.im" (WILL FAIL)
will cause the image analysis application to fail, because the image name begins with a number. The solution is to escape the name properly, e.g.:
mask = "'3clean_mask.im'"
Image Mask Handling
makemask facilitates the handling of image masks in CASA. As mentioned above, there are two basic mask formats: 1) one or more Boolean masks stored internally in the image file, and 2) images with zero and non-zero image values. makemask looks like:
#makemask :: Makes and manipulates image masks
mode = 'list' #Mask method (list, copy,expand,delete,setdefaultmask)
inpimage = '' #Name of input image.
To distinguish between Boolean internal masks and zero/non-zero images, makemask uses the syntax galaxy.image:mask0 for Boolean masks within an image: in this example, the Boolean mask mask0 within the image galaxy.image. Without the colon separator, the image itself will be treated as a zero/non-zero mask.mode=’list’ lists all the internal Boolean masks that are present in an image. The default masks can be set with mode=’setdefaultmask’ and they can be deleted with the mode=’delete’. The default mask is used when an image is displayed in the viewer and is used in all analysis tasks.mode=’copy’ lets a user copy a Boolean mask from one image to another image, or to write out the Boolean mask as a zero/non-zero image. The latter format is very useful when different masks are combined or manipulated. All the image analysis tools, in particular immath are applicable for such zero/non-zero masks as they act like normal images. makemask will always attempt to regrid the input mask to the output image.In addition mode=’copy’ can be used to merge masks and also accepts regions. E.g. to create a mask from a CASA region file, the input would look like:
#makemask :: Makes and manipulates image masks
mode = 'copy' #Mask method (list, copy,expand,delete,setdefaultmask)
inpimage = 'inputimage.im' #Name of input image.
inpmask = 'region.crtf' #mask(s) to be processed: image masks,T/F internal masks
#(Need to include parent image names),regions(for copy mode)
output = 'imagemask.im' #Name of output mask (imagename or imagename:internal_maskname)
overwrite = False #overwrite output if exists?
mode=’expand’ furthermore expands a mask in the spectral domain. It regrids first then stretches the edge channels. E.g. a one plane continuum image would be stretched to all planes of a data cube.
Image Analysis Tools¶
The CASA image analysis module contains an image analysis tool with numerous methods, as well as several higher level tasks. The tasks free users from the burden of resource management, and offer what many consider to be a more user-friendly interface available via the input <taskname> CASA command. In many cases, image analysis tasks are really just simple wrappers around analogous tool methods (e.g., the imcollapse task is just a relatively simple wrapper around the ia.collapse() tool method call), although in some cases, such as with the imregrid task, the mapping is not as simple, and much more goes on “under the hood” of a task.
Overview of Image Analysis Tool Functionality
At the heart of the image analysis module is the image analysis tool. An image analysis tool provides access to CASA images. Currently only single-precision, floating-point CASA images are supported by all methods in the image analysis tool and complex-valued images are supported by many, but not all, methods.
The default, global image analysis tool is named ia. New, initially-unattached image analysis tools can be created via
my_new_ia = iatool()
Image analysis tools also provide direct (native) access to FITS and Miriad images, although such access is read-only. These foreign formats to CASA format. For optimum processing speed, it is highly recommended to convert foreign formats to CASA images.
It is important to note that many methods return new image analysis tools that are attached to an image that method has created. Even if one does not intend on using this returned tool, it is important to capture it and run ia.done() on it or it will continue to use resources unnecessarily, e.g.
new_image_tool = ia.collapse("my_collapsed.im")
#do things with new_image_tool and then run done() on it
new_image_tool.done()
Tool Manipulation
ia.close(): Detach tool from image and perform required clean up.
ia.done(): Detach tool from image and perform required clean up and optionally removed attached image.
ia.isopen(): Determines if there is an image attached to the tool.
ia.newimage(): Create a new image analysis tool using an image.
ia.newimagefromarray(): Create a new image analysis tool from a numpy array.
ia.newimagefromfile(): Create a new image analysis tool using an image.
ia.newimagefromfits(): Create a new image analysis tool using a FITS image.
ia.newimagefromimage(): Create a new image analysis tool using an image.
ia.newimagefromshape(): Create a new image analysis tool using an image shape.
ia.open(): Attach the image analysis tool to the specified image.
ia.type(): Tool type. Always returns ‘image’.
FITS Conversion
There is functionality to interconvert between CASA images and FITS files. There is also native access to FITS files.
ia.fromfits(): Convert a FITS image file to a CASA image
ia.tofits(): Convert a CASA image to a FITS file.
ImageCreation
There are various ways to create CASA images from various data structures.
ia.fromarray(): Create a CASA image from a numpy array of pixel values.
ia.fromshape(): Create a CASA image of a specified shape.
ia.maketestimage(): Create a test image from a FITS file in the CASA data repository.
Image Destruction
ia.remove(): Delete the attached image from disk.
ia.removefile(): Delete the specified image from disk.
Image Interrogation
Various metadata and pixel data can be interrogated.
ia.beamarea(): Get the image synthesized beam area.
ia.boundingbox(): Get the bounding rectangular box which circumscribes the specified region.
ia.brightnessunit(): Get the image brightness unit.
ia.commonbeam(): For an image with multiple beams, compute the size of the smallest beam that circumscribes all of the image’s beams.
ia.getchunk(): Get pixel or mask values from (a specified rectangular region of) an image.
ia.getregion(): Get pixel or mask values from a specified region of an image.
ia.haslock(): Determines if the image has a lock associated with it.
ia.history(): Get the history information from an image.
ia.miscinfo(): Retrieve “miscellaneous” metadata associated with an image.
ia.name(): Get the image name.
ia.pixelvalue(): Get the pixel and mask values at a specified location of an image.
ia.restoringbeam(): Get information about the synthesized beam(s) of an image.
ia.shape(): Get image shape.
ia.summary(): Get various metadata of an image.
Manipulation of Image Metadata
ia.lock(): Acquire a lock on the attached image.
ia.rename(): Rename the image.
ia.rotatebeam(): Rotate the synthesized beam(s) of an image through a specified angle.
ia.setbrightnessunit(): Set image brightness unit.
ia.sethistory(): Add history records to an image.
ia.setmiscinfo(): Set image miscellaneous metadata.
ia.setrestoringbeam(): Set image synthesized beam(s).
ia.unlock(): Release the image lock.
Manipulation of Image Pixel and Pixel Mask Values
ia.calc(): Replace the pixel values in the attached image with the values determined from the specified LEL expression.
ia.calcmask(): Compute a pixel mask based on an LEL expression.
ia.insert(): Insert the pixel values of another image into an image.
ia.maskhandler(): Manipulate image pixel masks.
ia.modify(): Modify an image using a model specified by a component list.
ia.putchunk(): Set pixel values (in a specified rectrangular region) of an image.
ia.putregion(): Set pixel values in a specified region of an image.
ia.replacemaskedpixels(): Set masked pixel to a specified value.
ia.set(): Set pixel or mask values.
Operations on Images
Various operations can be performed on images which result in new images.
ia.addnoise(): Add noise to an image.
ia.boxcar(): Boxcar smooth an image along a specified axis.
ia.decimate(): Remove planes of an image.
ia.collapse(): Collapse image along specified axis, computing aggregate function of pixels along that axis.
ia.convolve(): Convolve an image with an array or with another image.
ia.continuumsub(): Subtract continuum emission in a spectral line image.
ia.convolve2d(): Convolve an image with a two-dimensional kernel.
ia.crop(): Crop pixels from the edge of an image.
ia.fft(): Fast Fourier Transform (FFT) the image.
ia.hanning(): Hanning smooth an image along a specified axis.
ia.imagecalc(): Create an image from an LEL expression.
ia.imageconcat(): Concatenate multiple images along a specified axis.
ia.makecomplex(): Create a complex-valued image from two float-valued images representing the real and imaginary values.
ia.pad(): Pad the edges of an image with pixels.
ia.pv(): Create a position-velocity image.
ia.pbcor(): Construct a primary beam corrected image.
ia.rebin(): Rebin pixel values by specified factors.
ia.regrid(): Regrid an image to a specified coordinate system.
ia.rotate(): Rotate the direction coordinate of an image.
ia.sepconvolve(): Convolve an image with a separable kernel.
ia.subimage(): Create an image by specifying a region of an image.
ia.transpose(): Transpose an image.
Image Analysis
ia.convertflux(): Interconvert between peak intensity and flux density for a specified Gaussian source.
ia.decompose(): Decompose complex source into individual two dimensional models.
ia.deconvolvecomponentlist(): Deconvolve a component list from the restoring beam.
ia.findsources(): Find strong point sources in an image.
ia.fitcomponents(): Fit two-dimensional models to the direction plane(s) of an image.
ia.fitprofile(): Fit one-dimensional models along an axis image.
ia.histograms(): Compute histograms from the pixel values of an image.
ia.maxfit(): Find maximum value in the direction coordinate and do a simple parabolic fit.
ia.moments(): Compute moments of an image.
ia.statistics(): Compute image statistics using various algorithms.
ia.twopointcorrelation(): compute two point autocorrelation functions from the image
Image Coordinates
The coordinate system of an image can be manipulated. Specific coordinate system values can be directly manipulated using the CASA coordinate system tool.
ia.adddegaxes(): Add degenerate axes to an image’s coordinate system.
ia.coordmeasures(): Convert from pixel to world coordinates, and return as a measure.
ia.coordsys(): Retrieve the image coordinate system as a CASA coordinate system tool.
ia.setcoordsys(): Replace the image’s coordinate system with another.
ia.topixel(): Convert from world to pixel coordinates.
ia.toworld(): Convert from pixel to world coordinates.
Miscellaneous
ia.makearray(): Create a numpy array of specified shape and value.
FITS Conversion
exportfits: Convert a CASA image to a FITS image.
importfits: Convert a FITS image to a CASA image.
Interrogation and Manipulation of Image Metadata
imhead: Summarize, interrogate, and modify image metadata
imhistory: List and append records to image history.
Operations on Images
Various operations can be performed on images which result in new images.
imcollapse: Collapse image along specified axis, computing aggregate function of pixels along that axis.
imcontsub: Subtract continuum emission in a spectral line image.
immath: Perform mathematical operations upon images.
immoments: Compute image moments.
impbcor: Construct a primary beam corrected image.
impv: Create a position-velocity image.
imrebin: Rebin pixel values by specified factors.
imregrid: Regrid an image to a specfied coordinate system.
imsmooth: Perform various two-dimensional convolutions.
imsubimage: Create an image by specifying a region of an image.
imtrans: Transpose an image.
specsmooth: Perform various one-dimensional convolutions.
Image Analysis
imfit: Fit two-dimensional models to the direction plane(s) of an image.
imstat: Compute image statistics using various algorithms.
imval: Interrogate pixel values.
rmfit: Compute rotation measure.
specfit: Fit one-dimensional models along a specified axis of an image.
specflux: Report spectral profile and calculate spectral flux over a user-specified region.
spxfit: Fit spectral index models along a specified axis of an image.
A persistent CASA image is stored on disk. Several files and subdirectories containing the image pixel data, mask data, and metadata are stored in a directory. The name of that directory is the name of the image.To access an existing persistent image, use the ia.open() method:
ia.open("my.im")
When you are finished with the image, it is important to close the tool so it no longer uses system resources:
ia.close()
It is also possible to create temporary images, which, if small enough, are stored completely in memory and destroyed when the user is finished with them. Creating such images is usually accomplished by running one of the image creation methods, and leaving the name of the output image blank (this is usually the default). So, for example, to create an image of a specified shape, one might run:
ia.fromshape(shape=[20,20,20])
As with persistent images, it is important to close the image analysis tool when finished with temporary images. In this case, the temporary image will be destroyed.
Persistent images can, in principle, be stored in a variety of ways. For example, the image could be stored row by row; this is the way that most older generation packages store images. It makes for very fast row by row access, but very slow in other directions (e.g. extract all the profiles along the third axis of an image). A CASA image is stored with what is called tiling. This means that small multi-dimensional chunks (a tile) are stored sequentially. It means that row by row access is a little slower, but access speed is essentially the same in all directions.
Here are some simple examples using image tools.
#access the CASA "test" FITS image and write it to a CASA image named "zz"
ia.maketestimage('zz',overwrite=true)
#print a summary to the logger and capture the summary metadata in variable "summary"
summary = ia.summary()
#evaluate image statistics and save the stats info to a variable called "stats"
stats = ia.statistics()
#create a rectangular region using the rg tool
box = rg.box([10,10], [50,50])
#create a subimage of that region, and name the resulting image "zz2"
#capture the new image tool attached to "zz2" in the variable "im2"
im2 = ia.subimage('zz2', box, overwrite=true)
#get statistics for zz2 and store the results in the variable "stats2"
stats2 = im2.statistics()
print "CLEANING UP OLD zz2.amp/zz2.phase IF THEY EXIST. IGNORE WARNINGS!"
ia.removefile('zz2.amp')
ia.removefile('zz2.phase')
#FFT subimage and store amp and phase
im2.fft(amp='zz2.amp',phase='zz2.phase')
#close image tools
im2.close()
ia.close()
Foreign Images
The image analysis tool also provides native, read-only access to some foreign image formats. Presently, these are FITS (Float, Double, Short and Long pixel values are supported) and Miriad. This means that you don’t have to convert the file to native CASA format in order to access the image. For example:
#Assumes environment variable is set
pathname = os.environ.get("CASAPATH")
pathname = pathname.split()[0]
datapath1 = pathname + "/data/demo/Images/imagetestimage.fits"
#Access FITS image
ia.open(datapath1)
ia.close()
#Access Miriad image
ia.open('im.mir')
ia.close()
#create a new image tool attached to the FITS image
ims = ia.newimagefromimage(infile=datapath1)
#create a region record representing the inner quarter of an image
innerquarter=rg.box([0.25,0.25],[0.75,0.75],frac=true)
#create a subimage of the inner quarter of the FITS image
subim = ims.subimage(region=innerquarter)
#done with the tools, release resources
ia.close()
ims.close()
In general, any parameter to a task or a tool method which accepts an image name will support CASA, FITS, or Miriad images.
There are some performance penalties of which you should be aware. First, because CASA images are tiled (see above), performance is the same regardless of how the images are accessed. In contrast, FITS and Miriad images are not tiled. This means that the performance when accessing these types of images will be poorer for certain operations. e.g., extracting a profile along the third axis of an image. Second, for FITS images, masked values are indicated via a “magic value”. This means that the mask is worked out on the fly every time the image is accessed.
If you find performance is poor or if you want a writable image, then use appropriate tool methods to convert the foreign format image to a CASA image.
Virtual Images
It is possible to have an image analysis tool that is not associated with a single persistent image; these are called “virtual’’ images. For example, with ia.imagecalc(), one can create an expression which may contain many images. You can write the result of the expression to a persistent image, but if you wish, you can also just maintain the expression, evaluating it each time it is needed - nothing is ever written out to disk in this case. There are other image methods like this (the documentation for each one explains what it does). The rules are:
If you specify the outfile or equivalent parameter, then the output image is always persistent with the name specified.
If you leave the outfile or equivalent parameter unset, then if possible, a virtual image will be created. Sometimes this virtual image will be an expression as in the example above (i.e. it references other images) or a temporary image in memory, or a temporary image on disk. The ia.summary() method will list the type of image. When you ia.close() that image tool, the virtual image will be destroyed.
If you leave the outfile or equivalent parameter unset, and the called method cannot create a virtual image, it will create a persistent image with a name of its choice (sometimes input plus function name).
A virtual image can always be written to disk as a persistent image with the ia.subimage() method.
Coordinate Systems
An image contains a coordinate system. A coordinate system tool is used to manipulate a coordinate system. An image tool allows you to recover the coordinate system into a coordinate system tool via the ia.coordsys() method. You can set a new image coordinate system with the ia.setcoordsys() method.
You can do some basic world to pixel and vice versa coordinate transformations via the image tool ia.topixel(), ia.toworld(), and ia.coordmeasures() methods.
Lattice Expression Language (LEL)
LEL allows you to create mathematical expressions involving images. For example, add the corresponding pixel values of two images, or multiply the miniumum value of one image by the square root of the pixel values of another image. The LEL syntax is quite rich and is described in detail on the Lattice Expression Language section.
IMPORTANT NOTE: Image names which contain “special” characters (eg, “+”, “-“, etc) must be properly escaped. See the Lattice names subsection of the Expressions section in the aforementioned document for details.
To produce an image that is the result of an LEL computation, use the ia.calc() or ia.imagecalc() image analysis tool methods. Here are some examples.
In this example the image analysis tool is attached to the persistent image named “zz”. This image’s name is used in an LEL expression which adds the pixel values of that image to the sine of the pixel values of that image (for trigonometric LEL functions, pixel values are taken to be in radians). Note that the ia.calc() method overwrites the pixel values of the attached image with the values computed by the LEL expression. To create a new image without overwriting the pixel values of the image associated with the image tool, use the ia.imagecalc() method.
ia.maketestimage('zz', overwrite=true)
#Make the minimum value zero
ia.calc('zz + min(zz)')
ia.close()
This example demonstrates ways of dealing with image names which have special characters.
ia.maketestimage("test-im", overwrite=true)
#escape special characters using a ""
im1 = ia.imagecalc(pixels='test-im + 5')
#or surround the entire image name with quotes
im2 = ia.imagecalc(pixels='"test-im" + 5')
#or
im3 = ia.imagecalc(pixels="'test-im' + 5")
im1.close()
im2.close()
im3.close()
ia.close()
Region Selection
A region designates a subset of pixels in the image in which one is interested. The region is selected based on coordinate information. Such a selection complements on-the-fly masks in which pixels are selected based on a mathematical expression which is tested against their values (see below). Regions may be specified in several ways. The region manager tool (default rg) has several methods for generating regions. These methods generally return a dictionary representation of a region which can be used as input for the region parameter in various image analysis tool methods and tasks. A region can also be specified by the box/chans/stokes selection parameters in tasks and tool methods which accept them. Regions can also be specified in a special format known as CASA region text format. This format allows for specifying of various region shapes and spectral and polarization extents. This specification can be placed in a file, and in this case, the region parameter can be set to the name of that file and the region information will be extracted. Alternatively, the region parameter can be set directly to the CRTF specification. The complete CRTF specification can be found in the “Region File Format” section.
Pixel Masks
A pixel mask is a set of boolean values which have a one-to-one correspondence with image pixels. A value of True indicates that pixel is “good” (i.e., should be used in computations), while a value of False indicates that pixel is “bad”. For example, blanked pixels in a FITS image are treated as “bad” by CASA. When such a file is imported into a CASA image, a pixel mask is created to reflect the badness of blanked pixels in the FITS image. For persistent CASA images, pixel masks are stored in the same directory in which other image information is stored.
If an image does not have a pixel mask associated with it, all of its pixels are treated as good by CASA.
A CASA image may contain any number of pixel masks and these masks can be managed via the ia.maskhandler() image analysis tool method. If an image contains multiple pixel masks, only a maximum of one mask will be used during a run of a task or tool method. This pixel mask is known as the “default” pixel mask. The default pixel mask can be set by running ia.maskhandler(set=”pixelmaskname”). You can also indicate that none of the image pixel masks should be applied by running ia.maskhandler (set=””). In this case, all pixels are considered to be good. Pixel masks can also be viewed in the output of the ia.summary() image analysis tool method and imhead task output.
The ia.putregion() image analysis tool method run with usemask=True can be used to change the values of the default pixel mask. The image analysis tool method ia.set() can also be used to set the values of the default pixel mask. The image analysis tool method ia.calcmask() can be used to create a new pixel mask based on a boolean LEL expression.
On The Fly Pixel Masks
Most image analysis tool methods and tasks accept a parameter named mask, which represents an OTF (on-the-fly) pixel mask that is computed for use by only that tool method or task (the exception being the ia.calcmask() image analysis tool method in which case a persistent pixel mask is attached to the image; see previous section). This parameter may be specified in one of two ways:
As an LEL boolean expression, or
as a single image name, in which case, pixel values >= 0.5 are treated as True (good) values, and all others are treated as False.
If the image has a default pixel mask, the mask used in the computation is the logical AND of the OTF pixel mask and default pixel mask. For example:
ia.maketestimage('zz', overwrite=true)
#create default pixel mask for which only positive valued pixels are good
ia.calcmask("zz>0")
#compute statistics by specifying an OTF mask, which gets ANDed with
#the default pixel mask, effectively making only pixels with values between 0 and 1 "good"
#for the statistics computation
stats = ia.statistics(mask="zz < 1")
ia.close()
The mask expression must in general conform in shape and coordinates with the input image.
A useful LEL function to use in conjunction with the mask parameter is indexin(). This enables the user to specify a mask based upon selected pixel coordinates or indices rather than image values. For example:
ia.fromshape(shape=[20])
#only pixels in the specified planes along the specified axis are considered good.
#prints [False False False False True True True True True True False False False False True False False False True True]
print ia.getregion(mask='indexin(0, [4:9, 14, 18:19])',getmask=true)
ia.close()
Regions As Pixel Masks
Regions, which have previously been discussed, are just another form of an OTF pixel mask, and in fact, if one specifies the region and mask parameters simultaneously, and the associated image also has a default pixel mask, all these three types of pixel masks are just ANDed together to form the pixel mask that is used in the resulting computation. One can even convert a region specification into a persistent pixel mask by specifying the region parameter in e.g., the ia.fromimage() image analysis tool method. The created image will have a default pixel mask that is a representation of the region specified (if the initial image had a default pixel mask, then that will be ANDed with the region specification to form the default pixel mask of the resulting image).
Lattice Expression Language¶
The Lattice Expression Language (LEL) makes it possible to do arithmetic on lattices (in particular on images [which are just lattices plus coordinates]). An expression can be seen as a lattice (or image) in itself. It can be used in any operation where a normal image is used.
To summarize, the following functionality is supported:
The common mathematical, comparison, and relational operators.
An extensive list of mathematical and logical functions.
Mixed data type arithmetic and automatic data type promotion.
Support of image masks.
Masking using boolean expressions.
Handling of masks in an expression.
Support of image regions.
Interface from both Python and C++.
The first section explains the syntax. The last sections show the interface to LEL using Python or C++. At the end some examples are given. If you like, you can go straight to the examples and hopefully immediately be able to do some basic things.
LEL operates on lattices, which are a generalization of arrays. As said above a particular type of lattice is an image; they will be used most often. Because lattices can be very large and usually reside on disk, an expression is only evaluated when a chunk of its data is requested. This is similar to reading only the requested chunk of data from a disk file.
LEL is quite efficient and can therefore be used well in C++ and Python code. Note however, that it can never be as efficient as carefully constructed C++ code.
LEL Expressions
A LEL expression can be as simple or as complex as one likes using the standard arithmetic, comparison, and logical operators. Parentheses can be used to group subexpressions. The operands in an expression can be lattices, constants, functions, and condition masks.
lat1 + 10
lat1 + 2 * max(lat2,1)
amp(lat1, lat2)
lat1 + mean(img[region1])
lat1 + mean(lat2[lat2>5 && lat2<10])
The last example shows how a boolean expression can be used to form a mask on a lattice. Only the pixels fulfilling the boolean condition will be used when calculating the mean. In general the result of a LEL expression is a lattice, but it can be a scalar too. If is is a scalar, it will be handled correctly by C++ and Python functions using it as the source in, say, an assignment to another lattice. LEL fully supports masks. In most cases the mask of a subexpression is formed by AND-ing the masks of its operands. It is fully explained in a later section.
LEL supports the following data types:
Bool
Float - single precision real (which includes integers)
Double - double precision real
Complex - single precision complex
DComplex - double precision complex
All these data types can be used for scalars and lattices.
LEL will do automatic data type promotion when needed. E.g. when a Double and a Complex are used in an operation, they will be promoted to DComplex. It is also possible to promote explicitly using the conversion functions (FLOAT, DOUBLE, COMPLEX and DCOMPLEX). These functions can also be used to demote a data type (e.g. convert from Double to Float), which can sometimes be useful for better performance.
Region is a specific data type. It indicates a region of any type (in pixel or world coordinates, relative, fractional). A region can only be applied to a lattice subexpression using operator [].
Constants
Scalar constants of the various data types can be formed as follows (which is similar to Python):
A Bool constant can be given as True or False.
A Float constant can be any integer or floating-point number. For example:
3
3.14
3.14e-2
A Double constant is a floating-point number using a D for the exponent. One can also use the
DOUBLEfunction. For example:
1d2
$3.14d-2
double(2)
The imaginary part of a Complex or DComplex constant is formed by a Float or Double constant immediately followed by a lowercase i. A full complex constant is formed by adding another constant as the real part. For example:
1.5 + 2i
2i+1.5$ is identical
Note that a full complex constant has to be enclosed in parentheses when, say, a multiplication is performed on it. For example:
2 * (1.5+2i)
The functions pi() and e() should be used to specify the constants pi and e. Note that they form a Double constant, so when using for example pi with a Float lattice, it could make a lot of sense to convert pi to a Float. Otherwise the lattice is converted to a Double, which is time-consuming. However, one may have very valid reasons to convert to Double, e.g. to ensure that the calculations are accurate enough.
Operators
The following operators can be used (with their normal meaning and precedence):
Unary + and -, can not be used with Bool operands.
Unary !
Logical NOT operator, can only be used with Bool operands. For a region it forms the complement.
Binary ^, *, /, %, +, and -
% is the modulo operator. E.g.
3%1.4results in0.2and-10%3results in-1.^ is the power operator.
All operators are left-associative, except ^ which is right-associative; thus
2^1^2results in2.Operator % can only be used for real operands, while the others can be used for real and complex operands.
Operator - can also be used for regions. It forms the difference of the left and right operand.
==, ! =, >, > =, <, and < =
For Bool operands only = = and ! = can be used. A Bool operand cannot be compared with a numeric operand. The comparison operators use the norm for complex values.
&& and | | && and ||
Logical AND and OR operator.
These operators can only be used with Bool operands. When used on a region && forms the intersection, while | | forms the union.
The precedence order is:
^
unary +, -, ! *, /, %
+, -
= = ,! = , > , > = , < , < =
&&
| |
Note that ^ has a higher precedence than the unary operators.For example, -3^2 results in -9.
The operands of these operators can be 2 scalars, 2 lattices, or a lattice and a scalar. When 2 lattices are used, they should in principle conform; i.e. they should have the same shape and coordinates. However, LEL will try if it can extend one lattice to make it conformant with the other. It can do that if both lattices have coordinates and if one lattice is a true subset of the other (thus if one lattice has all the coordinate axes of the other lattice and if those axes have the same length or have length 1). If so, LEL will add missing axes and/or stretch axes with length 1.
Functions
In the following tables the function names are shown in uppercase, while the result and argument types are shown in lowercase. Note, however, that function names are case-insensitive. All functions can have scalar and/or lattice arguments.When a function can have multiple arguments (e.g. atan2), the operands are automatically promoted where needed.
Mathematical functions
Several functions can operate on real or complex arguments. The data types of such arguments are given as ‘numeric’.
Double PI()
Returns the value of pi.
Double E()
Returns the value of e.
numeric SIN(numeric)
numeric SINH(numeric)
real ASIN(real)
numeric COS(numeric)
numeric COSH(numeric)
real ACOS(real)
real TAN(real)
real TANH(real)
real ATAN(real)
real ATAN2(real y, real x)
Returns ATAN(y/x) in correct quadrant.
numeric EXP(numeric)
numeric LOG(numeric)
Natural logarithm.
numeric LOG10(numeric)
numeric POW(numeric, numeric)
The same as operator ^.
numeric SQRT(numeric)
complex COMPLEX(real, real)
Create a complex number from two reals.
complex CONJ(complex)
real REAL(numeric)
Real value itself or real part of a complex number.
real IMAG(complex)
Imaginary part of a complex number.
real NORM(numeric)
real ABS(numeric), real AMPLITUDE(numeric)
Both find the amplitude of a complex number. If the numeric argument is real, imaginary part zero is assumed.
real ARG(complex), real PHASE(complex)
Both find the phase of a complex number.
numeric MIN(numeric, numeric)
numeric MAX(numeric, numeric)
Float SIGN(real)
Returns -1 for a negative value, 0 for zero, 1 for a positive value.
real ROUND(real)
Rounds the absolute value of the number. E.g. ROUND(-1.6) = -2.
real FLOOR(real)
Works towards negative infinity. E.g. FLOOR(-1.2) = -2
real CEIL(real)
Works towards positive infinity.
real FMOD(real, real)
The same as operator %.
Note that the trigonometric functions need their arguments in radians.
Scalar result functions
The result of these functions is a scalar.
double NELEMENTS(anytype)
Return number of elements in a lattice (1 for a scalar).
double NDIM(anytype)
Return dimensionality of a lattice (0 for a scalar).
double LENGTH(anytype, real axis)
Return length of a lattice axis (returns 1 for a scalar or if axis exceeds number of axes). Axis number is 1-relative.
Bool ANY(Bool)
Is any element true?
Bool ALL(Bool)
Are all elements true?
Double NTRUE(Bool)
Number of true elements.
Double NFALSE(Bool)
Number of false elements.
numeric SUM(numeric)
Return sum of all elements.
numeric MIN(numeric)
Return minimum of all elements.
numeric MAX(numeric)
Return maximum of all elements.
real MEDIAN(real)
Return median of a lattice. For smallish lattices (max. 512*512 elements) the median can be found in 1 pass. Other lattices usually require 2 passes.
real FRACTILE(real,float)
Return the fractile of a lattice at the fraction given by the second argument. A fraction of 0.5 is the same as the median. The fraction has to be between 0 and 1. For smallish lattices (max. 512*512 elements) the fractile can be found in 1 pass. Other lattices usually require 2 passes.
real FRACTILERANGE(real,float,float)
Return the range between the fractiles at the fraction given by the second and third argument. The fractions have to be between 0 and 1 and the second fraction has to be greater than the first one. The second fraction is optional and defaults to 1-fraction1. Thus:
FRACTILERANGE(lat, 0.1)
FRACTILERANGE(lat, 0.1, 0.9)
FRACTILE(lat,0.9) - FRACTILE(lat,0.1)
are equal, be it that the last one is about twice as slow. For smallish lattices (max. 512*512 elements) the fractile range can be found in 1 pass. Other lattices usually require 2 passes.
numeric MEAN(numeric)
Return mean of all elements.
numeric VARIANCE(numeric)
Return variance.
(sum((a(i) - mean(a))**2) / (nelements(a) - 1)`)
All calculations are done in double precision.
numeric STDDEV(numeric)
Return standard deviation (the square root of the variance).
real AVDEV(numeric)
Return average deviation.
(`sum(abs(a(i) - mean(a))) / nelements(a)`). All calculations are done in double precision.
Miscellaneous functions
numeric REBIN(numeric,[f1,f2,...])
Rebins the image using the given (integer) factors. It averages the pixels in each bin with shape [f1,f2,…]. Masked-off pixels are not taken into account. If all pixels in a bin are masked off, the resulting pixel will be masked off. The length of the factor list [f1,f2,…] has to match the dimensionality of the image. The factors do not need to divide the axes lengths evenly. Each factor can be a literal value, but it can also be any expression resulting in a real scalar value. For instance, for a 3-dimensional image:
rebin(lat,[2,2,1])
will halve the size of axis 1 and 2.
real AMP(real,real)
It returns the square root of the quadrature sum of the two arguments. Thus:
amp(lat1,lat2)
gives \(\sqrt{{lat}_1^2 + {lat}_2^2}\)
This can be used to form, for example, (biased) polarized intensity images when lat1 and lat2 are Stokes Q and U images.
real PA(real,real)
It returns a `position angle’ (in degrees) from the two lattices. That is,
pa(lat1,lat2)
gives \(180/\pi*atan2(lat1, lat2)/2\)
This can be used to form, for example, linear polarization position angle images when lat1 and lat2 are Stokes Q and U images, respectively.
real SPECTRALINDEX(real,real)
It returns a the spectral index made from the two lattices. That is,
log(s1/s2) / log(f1/f2)
where s1 and s2 are the source fluxes in the lattices and f1 and f2 are the frequencies of the spectral axes of both lattices. Similar to e.g. operator + the lattices do not need to have the same shape. One can be extended/stretched as needed.
anytype VALUE(anytype)
It returns the argument without its possible mask, thus it removes the mask from the argument. The section about mask handling discusses it in more detail.
Bool MASK(anytype)
It returns the mask of the argument. The section about mask handling discusses it in more detail.
Bool ISNAN(anytype)
It tests lattice elements on a NaN value and sets the corresponding output element to T if so; otherwise to F.
anytype REPLACE(anytype, anytype)
The first argument has to be a lattice (expression). The optional second argument can be a scalar or a lattice (expression). It defaults to 0. The result of the function is a copy of the first argument, where each masked-off element in the first argument is replaced by the corresponding element in the second argument. The result’s mask is a copy of the mask of the first argument.
replace (lat1, 0)
replace (lat1, lat2)
The first example replaces each masked-off element in lat1 by 0. The second example replaces it by the corresponding element in lat2. A possible mask of lat2 is not used.
anytype IIF(Bool, anytype, anytype)
The first argument is a boolean expression. If an element in it is true, the corresponding element from the second argument is taken, otherwise from the third argument. It is similar to the ternary ?: construct in C++. E.g.
iif (lat1>0, lat1, 0) same as max(lat1,0)
iif (sum(lat1)>0, lat1, lat2)
The examples shows that scalars and lattices can be freely mixed. When all arguments are scalars, the result is a scalar. Otherwise the result is a lattice. Note that the mask of the result is formed by combining the mask of the arguments in an appropriate way as explained in the section about mask handling.
Bool INDEXIN(real axis, set indices)
The first argument is a 1-relative axis number. The second argument is a set of indices. It returns a Bool array telling for each lattice element if the index of the given axis is contained in the set of indices.
The 1-relative indices have to be given as elements with integer values enclosed in square brackets and separated by commas. Each element can be a single index, an index range as start:end, or a strided index range as start:end:stride. The elements do not need to be ordered, but in a range start must be < = end. For example:
image[indexin(2, [3,4:8,10:20:2])]
masks image such that only the pixels with an index 3, 4, 5, 6, 7, 8, 10, 12, 14, 16, 18, or 20 on the second axis are set to True.
The following special syntax exists for this function.
INDEXi IN set
where i is the axis number. So the example above can also be written as:
image[index2 in [3,4:8,10:20:2]]
Negated versions of this function exist as:
INDEXNOTIN(axis, set)
INDEXi NOT IN set
Conversion functions
Float FLOAT(real)
Convert to single precision.
Double DOUBLE(real)
Convert to double precision.
Complex COMPLEX(numeric)
Convert to single precision complex. If the argument is real, the imaginary part is set to 0.
DComplex DCOMPLEX(numeric)
Convert to double precision complex. If the argument is real, the imaginary part is set to 0.
Bool BOOLEAN(region)
Convert to boolean. This can be useful to convert a region to a boolean lattice. Only a region in pixel coordinates can be converted, so in practice only an image mask can be converted.
Note that, where necessary, up-conversions are done automatically. Usually it may only be needed to do a down-conversion (e.g. Double to Float).
Lattice names
When a lattice (e.g. an image) is used in an expression, its name has to be given. The name can be given directly if it consists of the characters -.$~and alphanumeric characters.
If the name contains other characters or if it is a reserved word (currently only T and F are reserved), it has to be escaped. Escaping can be done by preceeding the special characters with a backslash or by enclosing the string in single or double quotes. E.g.
~/myimage.data
~/myimage.data\-old
~/myimage.data-old
LEL Masks¶
A boolean mask associated with an image indicates whether a pixel is good (mask value True) or bad (mask value False). If the mask value is False, then the image pixel is not used for computation (e.g. when finding the mean of the image).
An image can have zero (all pixels are good) or more masks. One mask can be designated as the default mask. By default it will be applied to the image (from Python, designation of the default mask is handled by the ia.maskhandler method of the Image tool).
When using LEL, the basic behaviour is that the default mask is used. However, by qualifying the image name with a suffix string, it is possible to specify that no mask or another mask should be used. The suffix is a colon followed by the word nomask or the name of the alternative mask.
myimage.data
myimage.data:nomask
'myimage.data:othermask'
The first example uses the default mask (if the image has one). The second example uses no mask (thus all pixels are designated good) and the third example uses mask othermask.
Note that if the image name is enclosed in quotes, the mask name should be enclosed too. It means that a colon cannot be part of an image name.
It is also possible to use a mask from another image like
myimage.data:nomask[myotherimage::othermask]
This syntax is explained in the section describing regions.
Lattice Condition Mask
We have seen in the previous section that lattices (in this case images) can have an associated mask. These masks are stored with the image - they are persistent.
It is also possible to create transient masks when a LEL expression is executed (dawn, usually). This is done with the operator [] and a boolean expression. For example,
sum( lat1[lat1<5 && lat1>10] )
creates a mask for lat1 indicating that only its elements fulfilling the boolean condition should be taken into account in the sum function. Note that the mask is local to that part of the expression. So in the expression
sum( lat1[lat1<5 && lat1>10] ) + sum(lat1)
the second sum function takes all elements into account. Masking can also be applied to more complex expressions and it is recursive.
(lat1+lat2)[lat3<lat4]
sum( lat1[lat1<5][lat1>10] )
(lat1 + lat2[lat3<lat4]) [lat1<5]
The first example applies the mask generated by the [] operator to the expression lat1+lat2. The second example shows the recursion (which ANDs the masks). It is effectively a (slower) implementation of the first example in this subsection. In the last example, the expression inside the parentheses is only evaluated where the condition [lat1<5] is true and the resulting expression has a mask associated with it.
Please note that it is possible to select pixels on an axis by means of the function INDEXIN (or by the INDEXi IN expression) as shown in the previous section about miscellaneous functions.
Mask Handling
As explained in the previous subsections, lattices can have a mask. Examples are a mask of good pixels in an image, a mask created using a boolean condition and the operator [], or a mask defining a region within its bounding box.A pixel is bad when the image has a mask and when the mask value for that pixel is False. Functions like max ignore the bad pixels.Note that in a
MeasurementSet a False mask value indicates a good visibility. Alas this is a historically grown distinction in radio-astronomy.Image masks are combined and propagated throughout an expression. E.g. when two lattices are added, the mask of the result is formed by and-ing the masks of the two lattices. That is, the resultant mask is True where the mask of lattice one is true AND the mask of lattice 2 is True. Otherwise, the resultant mask is False.In general the mask of a subexpression is formed
by and-ing the masks of the operands. This is true for e.g. +, *, atan2, etc.. However, there are a few special cases:
The mask created by
operator[condition]is formed by and-ing the condition result, the mask of the result, and the mask of the subexpression where the condition is applied to. For example, supposelat1andlat2each have a mask. Then in
sum( lat1[lat2<5] )
the sum function will only sum those elements for which the mask of lat1 and lat2 is valid and for which the condition is true.
The logical AND operator forms the resultant mask by looking at the result and the masks of the operands.
lat1[lat1<0 && lat2>0]
Let us say both lat1 and lat2 have masks. The operand lat1<0 is true if the mask of lat1 is true and the operand evaluates to true, otherwise it is false. Apply the same rule to the operand lat2 > 0. The AND operator gives true if the left and right operands are both true. If the left operand is false, the right operand is no longer relevant. It is, in fact, 3-valued logic with the values true, false, and undefined.
Thus, the full expression generates a lattice with a mask. The mask is true when the condition in the [] operator is true, and false otherwise. The values of the output lattice are only defined where its mask is true.
The logical OR operator works the same as the AND operator. If an operand has a true value the other operand can be ignored.
The mask of the result of the
replacefunction is a copy of the mask of its first operand. The mask of the second operand is not used at all.The
iiffunction has three operands. Depending on the condition, an element from the second or third operand has to be taken. The resultant mask is formed by the mask of the condition and-ed with the appropriate elements from the masks of the second or third operand.The
valuefunction returns the value without a mask, thus it removes the mask from a value. It has the same effect as theimage:nomaskconstruction discussed above. However, thevaluefunction is more general, because it can also be applied to a subexpression.The
maskfunction returns the mask of a value. The returned value is a boolean lattice and has no mask itself. When the value has no mask, it returns a mask consisting of all True values. When applied to an image, it returns its default mask.
Consider the following more or less equivalent examples:
value(image1)[mask(image2)]
image1:nomask[mask(image2)]
image1:nomask[image2::mask0]
The first two use the default mask of image2 as the mask for image1.The latter uses mask0 of image2 as the mask for image1. It is equivalent to the first two examples if mask0 is the default mask of image2.It is possible that the entire mask of a subexpression is false. For example, if the mean of such a subexpression is taken, the result is undefined. This is fully supported by LEL, because a scalar value also has a mask associated with it. One can see a masked-off
scalar as a lattice with an all false mask. Hence an operation involving an undefined scalar results in an undefined scalar. The following functions act as described below on fully masked-off lattices:
MEDIAN, MEAN, VARIANCE, STDDEV, AVDEV, MIN, MAX result in an undefined scalar:
NELEMENTS, NTRUE, NFALSE, SUM result in a scalar with value 0.
ANY results in a scalar with value F.
ALL results in a scalar with value T.
LENGTH, NDIM ignore the mask because only the shape of the lattice matters.
You should also be aware that if you remove a mask from an image, the values of the image that were previously masked bad may have values that are meaningless.
Mask Storage
In many of the expressions we have looked at in the examples, a mask has been generated. What happens to this mask and indeed the values of the expression depends upon the implementation. If for example, the function you are invoking with LEL writes out the result, then both the mask and result will be stored. On the other hand, it is possible to just use LEL expressions but never write out the results to disk. In this case, no data or mask is written to disk. You can read more about this in the interface section.
LEL Regions¶
A region-of-interest generally specifies a portion of a lattice which you are interested in for some astronomical purpose (e.g. what is the flux density of this source). Quite a rich variety of regions are supported in CASA. There are simple regions like a box or a polygon, and compound regions like unions and intersections. Regions may contain their own ``region masks’’. For example, with a 2-d polygon, the region is defined by the vertices, the bounding box and a mask which says whether a pixel inside the bounding box is inside of the polygon or outside of the polygon.
In addition, although masks and regions are used somewhat differently by the user, a mask is really a special kind of region; they are implemented with the same underlying code.
Like masks, regions can be persistently stored in image. Within Python, regions can be created using the various methods of the rg tool. Regions can also be defined in plain text files (see Region File Format).
We saw in the previous section how the condition operator [] could be used to generate masks with logical expressions. This operator has a further talent. A region of any type can be applied to a lattice with the [] operator. You can think of the region as also effectively being a logical expression. The only difference with what we have seen before is that it results in a lattice with the shape of the region’s bounding box. If the lattice or the region (as in the polygon above) has a
mask, they are and-ed to form the result’s mask.
All types of regions supported in CASA can be used, thus:
regions in pixel or world coordinates
in absolute, relative and/or fractional units
basic regions box, ellipsoid, and polygon
compound regions union, intersection, difference, and complement.
extension of a region or group of regions to higher dimensions
masks
The documentation in the classes LCRegion, LCSlicer, and WCRegion gives you more information about the various regions.
At this moment a region can not be defined in LEL itself. It is only possible to use regions predefined in an image or another table.
A predefined region can be used by specifying its name. There are three ways to specify a region name:
tablename::regionnameThe region is looked up in the given table (which will usually be an image) in which it is stored.::regionnameThe region is looked up in the last table used in the expression.regionnameIs usually equivalent to above. However, there is no syntactical difference between the name of a region and a lattice/image. Therefore LEL will first try if the name represents a lattice or image. If not, the name is supposed to be a region name. The prefix::in the previous way tells that the name should only be looked up as a region.
Examples are
myimage.data[reg1]
(myimage.data - otherimage)[::reg1]
(myimage.data - otherimage)[myimage.data::reg1]
myimage.data:nomask[myotherimage::othermask]
In the first example region reg1 is looked up in image myimage.data. It is assumed that reg1 is not the name of an image or lattice. It results in a lattice whose shape is the shape of the bounding box of the region. The mask of the result is the and of the region mask and the lattice mask.
In the second example it is stated explicitly that reg1 is a region by using the :: syntax. The region is looked up in otherimage, because that is the last table used in the expression. The result is a lattice with the shape of the bounding box of the region.
In the third example the region is looked up in myimage.data. Note that the this and the previous example also show that a region can be applied to a subexpression.
In the fourth example we have been very cunning. We have taken advantage of the fact that masks are special sorts of regions. We have told the image myimage.data not to apply any of its own masks. We have then used the [] operator to generate a mask from the mask stored in a different image, myotherimage. This effectively applies the mask from one image to another. Apart from copying the mask, this is the only way to do this.
Unions, intersections, differences and complements of regions can be generated and stored (in C++ and Python). However, it is also possible to form a union, etc. in LEL itself. However, that can only be done if the regions have the same type (i.e. both in world or in pixel coordinates).The following operators can be used:
reg1 || reg2to form the union.reg1 && reg2to form the intersection.reg1 - reg2to form the difference.!reg1to form the complement.
The normal CASA rules are used when a region is applied:
A region in world or relative coordinates can only be applied to an image (or a subexpression resulting in an image). Otherwise there is no way to convert it to absolute pixel coordinates.
The axes of a region in world coordinates have to be axes in the image (subexpression). However, the region can have fewer axes.
If a region has fewer axes than the image or lattice the region is automatically extended to the full image by taking the full length of the missing axes.
LEL Optimization¶
When giving a LEL expression, it is important to keep an eye on performance issues.
LEL itself will do some optimization:
As said in the introduction a LEL expression is evaluated in chunks. However, a scalar subexpression is executed only once when getting the first chunk. E.g. in
lat1 + mean(lat2)
the subexpression mean(lat2) is executed only once and not over and over again when the user gets chunks.
Often the exponent 2 is used in the
powfunction (or operator^). This is optimized by using multiplication instead of using the system pow function.When LEL finds a masked-off scalar in a subexpression, it does not evaluate the other operand. Instead it sets the result immediately to a masked-off scalar. Exceptions are the operators AND and OR and function
iif, because their masks depend on the operand values.
The user can optimize by specifying the expression carefully.
It is strongly recommended to combine scalars into a subexpression to avoid unnecessary scalar-lattice operations. E.g.
2 * lat1 * pi()
should be written as
lat1 * (2 * pi())
#or
2 * pi() * lat1
because in that way the scalars form a scalar subexpression which is calculated only once. Note that the subexpression parentheses are needed in the first case, because multiplications are done from left to right. In the future LEL will be optimized to shuffle the operands when possible and needed.
It is important to be careful with the automatic data type promotion of single precision lattices. Several scalar functions (e.g. pi) produce a double precision value, so using
piwith a single precision lattice causes the lattice to be promoted to double precision. If accuracy allows it, it is much better to convertpito single precision. E.g. assumelat1andlat2are single precision lattices.
atan2(lat1,lat2) + pi()/2
The result of atan2 is single precision, because both operands are single precision. However, pi is double precision, so the result of atan2 is promoted to double precision to make the addition possible. Specifying the expression as:
atan2(lat1,lat2) + float(pi())/2
avoids that (expensive) data type promotion.
POW(LAT,2)orLAT``^``2is faster thanLAT*LAT, because it accesses latticeLATonly once.SQRT(LAT)is faster thanLAT``^``0.5orPOW(LAT,0.5)POW(U,2) + POW(V,2) < 1000``^``2is considerably faster thanSQRT(SQUARE(U) + SQUARE(V)) < 1000, because it avoids theSQRTfunction.LEL can be used with disk-based lattices and/or memory-based lattices. When used with memory-based lattices it is better to make subexpressions the first operand in another subexpression or a function. E.g.
lat1*lat2 + lat3is better thanlat3 + lat1*lat2The reason is that in the first case no copy needs to be made of the lattice data which already reside in memory. All LEL operators and functions try to reference the data of their latter operands instead of making a copy. In general this optimization does not apply to LEL expression. However, when using the true C++ interface to classes likeLatticeExprNode, one can easily use memory-based lattices. In that case it can be advantageous to pay attention to this optimization.
LEL Interface¶
There are two interfaces to LEL. One is from Python and the other from C++. It depends upon your needs which one you access. Most high level users of CASA will access LEL only via the Python interface.
Simple String Expressions
The ia.imagecalc method evaluates the expression and stores the result and mask in an output image. If you specify the output image name, it is written to a disk file of that name. If you don’t give it, the output image is never written out; it is evaluated every time an action (like ia.statistics) is requested.
im = ia.imagecalc(outfile='outimage', pixels='inimage1+inimage2');
im.statistics();
The first command creates an image file outimage filling it with the sum of the input images. The second command does statistics on that new image. Writing it as
im = ia.imagecalc(pixels='inimage1+inimage2');
im.statistics();
would do the same with the exception of creating the output image. Instead the created image is transient; it only lives as an expression and each time it is used the expression is evaluated.
We can use the method ia.calc on an already existing image. Thus
ia.open('ngc1213');
ia.calc('ngc1213^2');
would replace the pixels by the square of their value in the opened image.
Sometimes you need to double quote the file names in your expression. For example, if the images reside in a different directory as in this example.
im = ia.imagecalc ('"dir1/im1" + "/nfs/data/im2"');
C++ interface
This consists of 2 parts.
The function
commandin Images/ImageExprParse.h can be used to execute a LEL command. The result is a LatticeExprNode object. This example does the same as the Python one shown above. E.g.
LatticeExprNode seltab1 = ImageExprParse::command
("imagein1 + imagein2");
The other interface is a true C++ interface having the advantage that C++ variables can be used. Class LatticeExprNode contains functions to form an expression. The same operators and functions as in the command interface are available. For example:
Float clipValue = 10;
PagedImage<Float> image("imagein");
LatticeExpr<Float> expr(min(image,clipValue));
forms an expression to clip the image. Note that the expression is written as a normal C++ expression. The overloaded operators and functions in class LatticeExprNode ensure that the expression is formed in the correct way. Note that a LatticeExprNode object is usually automatically converted to a templated LatticeExpr object, which makes it possible to use it as a normal Lattice. So far the expression is only formed, but not evaluated. Evaluation is only done when the expression is
used in an operation, e.g. as the source of the copy operation shown below.
PagedImage<Float> imout("imageout");
imout.copyData (expr);
LEL Examples¶
The following examples show some LEL expressions (equally valid in C++ or Glish).
Note that LEL is readonly; i.e. it does not change any value in the images given. A function in the image client has to be used to do something with the result (e.g. storing in another image).
lat1+lat2– adds 2 latticesmean(myimage:nomask) --results in a scalar value giving the mean of the image. No mask is used for the image, thus all pixels are used. The scalar value can be used as a lattice. E.g. it can be used as the source in theimagefunctionreplacemaskedpixelsto set all masked-off elements of a lattice to the mean.complex(lat1,lat2) --results in a complex lattice formed bylat1as the real part andlat2as the imaginary part.min(lat1, 2*mean(lat1)) --results in a lattice wherelat1is clipped at twice its mean value.min(myimage, 2*mean(mymage[myregion])) --results in an image wheremyimageis clipped at twice the mean value of regionmyregionin the image..lat1[lat1>2*min(lat1)]– results in a lattice with a mask. Only the pixels greater than twice the minimum are valid.replace(lat1)– results in a lattice where each masked-off element inlat1is replaced by 0.iif(lat1<mean(lat1),lat1*2,lat1/2) --results in a lattice where the elements less than the mean are doubled and the elements greater or equal to the mean are divided by 2.
Here follows a sample Glish session showing some of the LEL capabilities and how Glish variables can be used in LEL.
duw01> glish -l image.g
- a := array(1:50,5,10) # make some data
- global im1 := imagefromarray('im1', a); # fill an image with it
- im1.shape()
[5 10]
- local pixels, mask
- im1.getregion(pixels, mask); # get pixels and mask
- mask[1,1] := F # set some mask elements to False
- mask[3,4] := F
- im1.putregion(mask=mask); # put new mask back
- global reg:=drm.box([1,1],[4,4]); # a box region
- im2 := imagecalc(pixels='$im1[$reg]') # read-only image applying region
- local pixels2, mask2
- im2.getregion(pixels2, mask2); # get the pixels and mask
- print pixels2
[[1:4,]
1 6 11 16
2 7 12 17
3 8 13 18
4 9 14 19]
- print mask2
[[1:4,]
F T T T
T T T T
T T T F
T T T T]
- im1.replacemaskedpixels ('mean(im2)'); # replace masked-off values
- im1.getregion (pixels2, mask2); # by mean of masked-on in im2
- print pixels2
[[1:5,]
10.0714283 6 11 16 21 26 31 36 41 46
2 7 12 17 22 27 32 37 42 47
3 8 13 10.0714283 23 28 33 38 43 48
4 9 14 19 24 29 34 39 44 49
5 10 15 20 25 30 35 40 45 50]