image

class image[source]

Operations on images

An Image tool provides access to  images. Currently only single precision floating point  images are supported by all methods in the Image and complex-valued images are supported by many, but not all, methods.

Image tools also provide direct (native) access to  and Miriad images. You can also convert these foreign formats to  format (for optimum processing speed).

It is important to note that many methods return new image 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 done() on it or it will continue to use resources unnecessarily, eg

Methods Summary

adddegaxes

This method adds degenerate axes (i.e.

addnoise

This function adds noise to the image.

beamarea

Get the area of the image’s restoring beam.

beamforconvolvedsize

Determine the size of the beam necessary to convolve with the given source to reach the given convolved (source+beam) size.

boundingbox

This function finds the bounding box of a region of interest when it is applied to a particular image.

boxcar

This application performs boxcar convolution of one axis of an image defined by

brightnessunit

This function gets the image brightness unit.

calc

This function is used to evaluate a mathematical expression involving  images, assigning the result to the current (already existing) image.

calcmask

This method is used to create a new pixel mask via a boolean LEL expression.

close

This function closes the .

collapse

This method collapses an image along a specified axis or set of axes of length N pixels to a single pixel on each specified axis.

commonbeam

Determine a beam to which all beams in an image can be convolved.

continuumsub

This function packages the relevant image tool functionality for simple specification and application of image plane continuum subtraction.

convertflux

This function interconverts between peak intensity and flux density for a Gaussian component.

convolve

This function performs Fourier-based convolution of an image by the specified kernel.

convolve2d

This function performs Fourier-based convolution of an  using the provided 2D kernel.

coordmeasures

You can use this function to get the world coordinates for a specified absolute pixel coordinate in the image.

coordsys

This function returns the Coordinate System of an image in a Coordsys tool.

crop

This method crops masked slices from the perimeter of an image.

decimate

This application removes planes along the specified axis of an image.

decompose

This function is an image decomposition tool that performs several tasks, with the end result being that a strongly blended image is separated into components - both in the sense that it determines the parameters for each component (assuming a Gaussian model) and that it physically assigns each pixel in the image to an individual object.

deconvolvecomponentlist

This method deconvolves (a record representation of) a Componentlist tool from the restoring beam, returning (a record representation of) a new Componentlist tool.

deconvolvefrombeam

This is a helper function.

deviation

This application creates an image that reflects the statistics of the input image. The output image has

dohistory

This allows control over if tool methods record history of what parameters they were called with to the input and/or output image.

done

When the user no longer needs to use an , calling this function will free up its resources.

fft

This method fast Fourier Transforms the supplied image to the Fourier plane.

findsources

This function finds strong point sources in the image.

fitcomponents

OVERVIEW

fitprofile

This application simultaneously fits any number of gaussian singlets, any number of lorentzian singlets, and any number of gaussian multiplets, and/or a polynomial to one dimensional profiles using the non-linear, least squares Levenberg-Marquardt algorithm.

fitsheader

This method constructs a FITS header dictionary which describes the current image.

fromarray

This application converts a numerical numpy array of any size and dimensionality into a CASA image.

fromcomplist

This method allows one to create an image based on a component list.

fromfits

This function is used to convert a FITS disk image file (Float, Double, Short, Long are supported) to an  .

fromimage

This function applies a  to an , creates a new  containing the (sub)image, and associates the  with it.

fromrecord

You can convert an associated image to a record (torecord) or imagepol tool functions will sometimes give you a record.

fromshape

This function creates a CASA image with the specified shape.

getchunk

This function returns the pixels (or optionally the pixel mask) from the attached image between blc and trc, inclusive.

getprofile

This application returns information on a one-dimensional profile taken along a specified image axis.

getregion

This function recovers the image pixel or pixel mask values in the given region of interest.

getslice

This function returns a 1-D slice (the pixels and opionally the pixel mask) from the attached image.

hanning

This application performs Hanning convolution of one axis of an image defined by

haslock

This function can be used to find out whether the image has a read or a write lock set.

histograms

This method computes histograms of the pixel values in the image.

history

This method interogates the history of an image.

image

image method

imagecalc

This method is used to evaluate a mathematical expression involving existing images.

imageconcat

This application is used to concatenate two or more input CASA images into one output image.

insert

This function inserts the specified image (or part of it) into the image referenced by this tool.

isconform

Returns True if the shape, coordinate system, and axes order of the specified image matches the current image.

isopen

This method returns True if the image tool has an attached image.

ispersistent

This function can be used to find out whether the image is persistent on disk or not.

lock

This function can be used to acquire a Read or a Read/Write lock on the .

makearray

This function takes two arguments.

makecomplex

This function combines the current image with another image to make a complex image.

maketestimage

This function converts a FITS file resident in the  system into a  image.

maskhandler

This method is used to manage or handle pixel masks .

maxfit

This function finds the pixel with the maximum value in the region, and then uses function findsources to generate a Componentlist with one component.

miscinfo

A   can accumulate miscellaneous information during its lifetime.

modify

This function applies a model of the sky to the image.

moments

The primary goal of this function is to enable you to analyze a multi-dimensional image by generating moments of a specified axis.

name

This function returns the name of the  By default, this function returns the full absolute path of the .

newimage

This method is identical to ia.newimagefromfile().

newimagefromarray

This application converts a numpy array of any size into a CASA image.

newimagefromfile

This method returns an image analysis tool associated with the specified image.

newimagefromfits

This function is used to convert a FITS disk image file (Float, Double, Short, Long are supported) to an  .

newimagefromimage

This function applies a  to a disk , creates a new  containing the (sub)image, and associates a new  with it.

newimagefromshape

This function creates a CASA image with the specified shape.

open

This method detaches from the current image (if an image is attached to the tool), and reattaches it (opens) to the new image.

pad

This method pads the directional plane of an image with a specified number of pixels on each side.

pbcor

Correct an image for primary beam attenuation using an image of the primary beam pattern.

pixeltype

This application returns the data type of the pixels of the attached image as a string.

pixelvalue

This function gets the value of the image and the mask at the specified pixel coordinate.

putchunk

This method puts an array into the image (which must be writable, eg it will fail on FITS images).

putregion

This function replaces data and/or pixel mask values in the image in the specified region.

pv

Create a position-velocity image by specifying either two points between which a slice is taken in the direction coordinate or a center, position angle, and length describing the slice.

rebin

This application rebins the current image by the specified integer binning factors for each axis.

regrid

This function regrids the current image onto a grid specified by the given coordinate system.

remove

This function first closes the  which detaches it from its underlying .

removefile

This function deletes the specified image file.

rename

This function renames the  associated with the .

replacemaskedpixels

This application replaces the values of all pixels whose total input mask (default input  and OTF mask) is bad (F) with the specified value.

restoringbeam

This function returns the restoring beam(s), if any, of the attached image.

rotate

This function rotates two axes of an image.

rotatebeam

This method rotates the attached image’s beam(s) counterclockwise through the specified angle.

sepconvolve

This function does Fourier-based convolution of an  by a specified separable kernel.

set

This function replaces data and/or mask values within the image in the specified .

setbrightnessunit

This function sets the image brightness unit.

setcoordsys

This method replaces the coordinate system in the image.

sethistory

A   can accumulate history information from an input  file or by you writing something into it explicitly with this function.

setmiscinfo

A CASA image can include user-specified or miscellaneous metadata.

setrestoringbeam

This method sets the restoring beam(s) for an image or removes an existing beam(s).

shape

The shape of an image is a vector holding the length of each axis of the image.

statistics

This method computes statistics from the pixel values in the image.

subimage

This function copies all or part of the image to another on-the-fly Image tool.

summary

This function summarizes various metadata such as shape, Coordinate System, restoring beams, and masks.

tofits

This function converts the image into a  file.

topixel

This method converts from world to pixel coordinates.

torecord

You can convert an associated image to a record for manipulation or passing it to inputs of other methods of other tools.

toworld

This method converts between pixel and world coordinates.

transpose

This method transposes the axes in the input image to the specified order.

twopointcorrelation

This function computes two-point auto-correlation functions from an image.

type

This function returns the string ’image’.

unlock

This function releases any lock set on the  (and also flushes any outstanding I/O to disk).

adddegaxes(outfile='', direction=False, spectral=False, stokes='', linear=False, tabular=False, overwrite=False, silent=False)[source]

This method adds degenerate axes (i.e. axes of length 1) of the specified type. Sometimes this can be useful although you will generally need to modify the coordinate system of the added axis to give it the coordinate you want (do this with the Coordsys ). This method supports both float and complex valued images.

You specify which type of axes you want to add. You can’t add an axis type that already exists in the image. For the Stokes axis, the allowed value (a string such as I, Q, XX, RR) can be found in the Coordsys newcoordsys function documentation.

If outfile is given, the image is written to the specified disk file. If outfile is unset, the on-the-fly Image  returned by the function is associated with a temporary image. This temporary image may be in memory or on disk, depending on its size. When you destroy the generated Image  (with the done function) this temporary image is deleted.

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • direction (bool=False) - Add direction axes?

  • spectral (bool=False) - Add spectral axis?

  • stokes (string='') - Add Stokes axis? Default is empty string.

  • linear (bool=False) - Add linear axis?

  • tabular (bool=False) - Add tabular axis?

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • silent (bool=False) - Skip silently existing axes?

Returns

image

Examples

#
print "\t----\t adddegaxes Ex 1 \t----"
ia.maketestimage() 
print ia.shape()
#[113L, 76L]
mycs=ia.coordsys()
print mycs.axiscoordinatetypes()
#['Direction', 'Direction']
mycs.done()
im2 = ia.adddegaxes(spectral=True)
print im2.shape()
#[113L, 76L, 1L]
mycs=im2.coordsys()
print mycs.axiscoordinatetypes()
['Direction', 'Direction', 'Spectral']
mycs.done()
im3 = im2.adddegaxes(stokes='Q')
print im3.shape()
#[113L, 76L, 1L, 1L]
mycs = im3.coordsys()
print mycs.axiscoordinatetypes()
#['Direction', 'Direction', 'Spectral', 'Stokes']
mycs.done()
im2.done()
im3.done()
ia.close()
#


In this example, all the images are virtual (temporary images).
addnoise(type='normal', pars=[0.0, 1.0], region='', zero=False, seeds='')[source]

This function adds noise to the image. You may zero the image first before the noise is added if you wish.

The noise can be drawn from one of many distributions.

For each distribution, you must supply the type via the type argument (minimum match is active) and parameters via the pars argument. Briefly:

  • binomial – the binomial distribution models successfully drawing items from a pool. Specify two parameters, \(n\) and \(p\), respectively. \(n\) is the number of items in the pool, and \(p\), is the probability of each item being successfully drawn. It is required that \(n \> 0\) and \(0 \le p \le 1\).

  • discreteuniform – models a uniform random variable over the closed interval. Specify two parameters, the low and high values, respectively. The low parameter is the lowest possible return value and the high parameter is the highest. It is required that \(low \< high\).

  • erlang – Specify two parameters, the mean and variance, respectively. It is required that the mean is non-zero and the variance is positive.

  • geometric – Specify one parameter, the probability. It is required that \(0 \le probability \< 1\).

  • hypergeometric – Specify two parameters, the mean and the variance. It is required that the variance is positive and that the mean is non-zero and not bigger than the square-root of the variance.

  • normal – Specify two parameters, the mean and the variance. It is required that the variance is positive.

  • lognormal – Specify two parameters, the mean and the variance. It is required that the supplied variance is positive and that the mean is non-zero.

  • negativeexponential – Supply one parameter, the mean.

  • poisson – Specify one parameter, the mean. It is required that the mean is non-negative.

  • uniform – Model a uniform random variable over a closed interval. Specify two parameters, the low and high values. The low parameter is the lowest possible return value and the high parameter can never be returned. It is required that \(low \< high\).

  • weibull – Specify two parameters, alpha and beta. It is required that the alpha parameter is not zero.

    The random number generator seeds may be specified as an array of integers. Only the first two values are used. If none or a single value is provided, the necessary remaining value(s) are generated based on the current time, using the algorithm

    seedBase = 1e7\*MJD
    seed[1] = (Int)seedBase;
    # and if seed[0] is also not supplied
    seed[0] = (Int)((1e7\*(seedBase - seed[1])))
    

    where MJD is the Modidfied Julian Day.

Parameters

  • type (string='normal') - Type of distribution, normal

  • pars (doubleVec=[0.0, 1.0]) - Parameters of distribution

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • zero (bool=False) - Zero image first?

  • seeds (intVec='') - Seeds to use for the random number generator. If not specified, seeds are created based on the current time.

Returns

bool

Examples

ia.maketestimage()
ia.addnoise(type='normal', pars=[0.5, 1], zero=True, seeds=[1,2]) 
ia.statistics()
ia.close()


A test image is created, zeroed, and noise of mean 0.5 and variance 1
from a normal distribution added. Specifying the same combination of seeds
ensures the same random number sequence is created each time addnoise is called.
To have different sequences created during the same casapy session, use the default
value (which is an empty list).
beamarea(channel=-1, polarization=-1)[source]

Get the area of the image’s restoring beam. If a non-negative channel and non-negative polarization are specified, the area for the beam associated with that channel and polarization will be returned. The return value will be a dictionary containing the keys ’arcsec2’ and ’pixels’, and the associated values will be the beam area in arcsec2 and in pixels, respectively. If both channel and polarization are set to negative values, then a dictionary with the same keys will be returned, and the values will be either scalars (if the image has a single traditional beam) or arrays if the image has multiple beams. In the latter case, the arrays will have shapes indicative of the number of channels and number of polarizations. If the image has a spectral axis but not a polarization axis, the returned arrays will be one dimensional and have a length equal to the number of channels. Similarly, if the image has a polarization axis but not a spectral axis, the arrays will be one dimensional and have a lenghts equal to the number of polarizations. If the image has both a spectral and polarization axis, the returned arrays will be two dimensional with shape (m, n) where m is the length of the first channel or polarization axis, and n is the length of the second channel or polarization axis. So, if an image has shape [200, 200, 10, 4] with 10 channels and 4 stokes, the returned arrays will have shapes of (10, 4) representing the spectral axis as the first axis and the polarization axis as the second. If, instead, the image shape is [200, 200, 4, 10] again with 10 channels and 4 stokes, the shape of the returned arrays will be (4, 10) since the polarization axis precedes the spectral axis.

Parameters

  • channel (int=-1) - The zero-based spectral channel number for a per-plane beam. Default -1

  • polarization (int=-1) - The zero-based polarization plane number for a per-plane beam. Default -1

Returns

record

beamforconvolvedsize(source='', convolved='')[source]

Determine the size of the beam necessary to convolve with the given source to reach the given convolved (source+beam) size. Because the problem is completely specified by the input parameters, no image needs to be attached to the associated tool; eg ia.open() need not be called prior to calling this method.

Parameters

  • source (variant='') - Three quantities that define the deconvolved source major axis, minor axis and position angle

  • convolved (variant='') - Three quantities that define the convolved source (source+beam) major axis, minor axis and position angle. Do not specify if beam is specified.

Returns

record

Examples

# get the beam necessary to convolve the specified source with to achieve the target convolved source size (source convolved with beam).
beam = ia.beamforconvolvedsize(source=["1arcsec", "1arcsec", "0deg"], convolved="3arcsec", "2arcsec", "45deg"])
boundingbox(region='')[source]

This function finds the bounding box of a region of interest when it is applied to a particular image. Both float and complex valued images are supported. It is returned in a record which has fields ’blc’, ’trc’, ’inc’, ’bbShape’, ’regionShape’, ’imageShape’, ’blcf’ and ’trcf’ containing the bottom-left corner, the top-right corner (in absolute image pixel coordinates), the increment (stride) of the region, the shape of the boundingbox, the shape of the region, the shape of the image, the blc in formatted absolute world coordinates and the trc in formatted absolute world coordinates, respectively.

Note that the shape of the bounding box will be different from the shape of the region if a non-unit stride (increment) is involved (see the example below).

Note that the integer size of the elements in blc, trc, inc, regionShape, bbShape, and imageShape are 32 bits, even on a 64 bit machine. This means that, on 64 bit machines, you may have to convert them to 64 bit ints using eg numpy.int64, before being able to use them as direct input to other methods such as ia.getchunk().

Parameters

  • region ({record, string}='') - Region selection. Default is to use the full image.

Returns

record

Examples

#
print "\t----\t boundingbox Ex 1 \t----"
ia.maketestimage()                           # Create image tool
x=['3pix','6pix','9pix','6pix','5pix','5pix','3pix'] # X vector in abs pixels
y=['3pix','4pix','7pix','9pix','7pix','5pix','5pix'] # Y vector in abs pixels
mycs = ia.coordsys()
r1=rg.wpolygon(x=x,y=y,csys=mycs.torecord()) # Create polygonal world region
mycs.done()
bb = ia.boundingbox(r1)                      # Find bounding box
print bb
#{'regionShape': array([7, 7]), 'trc': array([9, 9]),
# 'imageShape': array([113, 76]),
# 'blcf': '00:00:27.733, -00.06.48.000',
# 'trcf': '00:00:24.533, -00.05.36.000', 'bbShape': array([7, 7]),
# 'blc': array([3, 3]), 'inc': array([1, 1])}
ia.close()
#
boxcar(outfile='', region='', mask='', axis=-1, width=2, drop=True, dmethod='copy', overwrite=False, stretch=False)[source]

This application performs boxcar convolution of one axis of an image defined by

z[i] = (y[i] + y[i+i] + … + y[i+w-1])/w

where z[i] is the value at pixel i in the box car smoothed image, y[k] is the pixel value of the input image at pixel k, and w is a postivie integer representing the width of the boxcar in pixels. Both float and complex valued images are supported. The length of the axis along which the convolution is to occur must be at least w pixels in the selected region, unless decimation using the mean function is chosen in which case the axis length must be at least 2w (see below). Masked pixel values are set to zero prior to convolution. All nondefault pixel masks are ignored during the calculation. The convolution is done in the image domain (i.e., not with an FFT).

If drop=False (no decimation), the length of the output axis will be equal to the length of the input axis - w + 1. The pixel mask, ORed with the OTF mask if specified, is copied from the selected region of the input image to the output image. Thus for example, if the selected region in the input image has six planes along the convolution axis, if the specified boxcar width is 2, and if the pixel values, which are all unmasked, on a slice along this axis are [1, 2, 5, 10, 17, 26], then the corresponding output slice will be of length five and the output pixel values will be [1.5, 3.5, 7.5, 13.5, 21.5].

If drop=True and dmethod=”copy”, the output image is the image calculated if drop=True, except that only every wth plane is kept. Both the pixel and mask values of these planes are copied directly to the output image, without further processing. Thus for example, if the selected region in the input image has six planes along the convolution axis, the boxcar width is chosen to be 2, and if the pixel values, which are all unmasked, on a slice along this axis are [1, 2, 5, 10, 17, 26], the corresponding output pixel values will be [1.5, 7.5, 21.5].

If drop=True and dmethod=”mean”, first the image described in the drop=False case is calculated. Then, the ith plane of the output image is calculated by averaging the iw to the (i+1)*w-1 planes of this intermediate image. Thus, for example, if the selected region in the input image has six planes along the convolution axis, the boxcar width is chosen to be 2, and if the pixel values, which are all unmasked, on a slice along this axis are [1, 2, 5, 10, 17, 26], then the corresponding output pixel values will be [2.5, 10.5]. Any pixels at the end of the plane of the intermediate image that do not fall into a complete bin of width w are ignored. Masked values are taken into consideration when forming this average, so if one of the values is masked, it is not used in the average. If at least one of the values in the intermediate image bin is not masked, the corresponding output pixel will not be masked.

The smoothed image is written to disk with name outfile, if specified. If not, no image is written but the image is still accessible via the returned image analysis tool (see below).

This method always returns an image analysis tool which is attached to the smoothed image. This tool should always be captured and closed after any desired manipulations have been done. Closing the tool frees up system resources (eg memory), eg,

smoothedim = ia.boxcar(...)
# do things (or not) with smoothedim
...
# close the returned tool promptly upon finishing with it.
smoothedim.done()

Parameters

  • outfile (string='') - Output image file name. Default is none.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • axis (int=-1) - Zero-based axis along which to convolve. ia.coordsys().names() gives the order of the axes in the image. Less than 0 means use the spectral axis if there is one, if not an exception is thrown.

  • width (int=2) - Width of the boxcar in pixels.

  • drop (bool=True) - Drop every nth pixel on output, where n is the width of the boxcar?

  • dmethod (string='copy') - If drop=True, method to use in plane decimation. “copy”: direct copy of every second plane, “m(ean)”: average planes n*i through n*(i+1) - 1 (inclusive) in the smoothed, non-decimated image to form plane i in the output image.

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • stretch (bool=False) - Stretch the mask if necessary and possible? See help par.stretch. Default False

Returns

image

Examples

ia.open("mynonsmoothed.im")
# smooth the spectral axis by 3 pixels, say it's axis 2 and only
# write every other pixel
boxcar = ia.boxcar(outfile="myboxcarsmoothed.im", axis=2, drop=True,
width=3, dmethod="c" overwrite=True)
# done with input
ia.done()
# do something with the output image, get statistics say
stats = boxcar.statistics()
# close the result image
boxcar.done()
brightnessunit()[source]

This function gets the image brightness unit. Both float and complex valued images are supported.

calc(pixels='', verbose=True)[source]

This function is used to evaluate a mathematical expression involving  images, assigning the result to the current (already existing) image. Both float and complex valued images are supported, although the image which results from the calculation must have the same type of pixel values as the image already attached to the tool. That is, one cannot create a complex valued image using this method if the associated ia tool is currently attached to a float valued image. It complements the imagecalc function which returns a newly constructed on-the-fly image tool. See which describes the the syntax and functionality in detail.

If the expression, supplied via the pixels argument, is not a scalar, the shapes and coordinates of the image and expression must conform.

If the image (that associated with the tool) has a , then only pixels for which that mask is good will be changed. See the function maskhandler for managing image .

Note that when multiple image are used in the expression, there is no garauntee about which of those images will be used to create the header of the output image. Therefore, one may have to modify the output header as needed if the input headers differ.

See the related functions set and putregion.

Parameters

  • pixels (string='') - LEL expression

  • verbose (bool=True) - Emit possibly useful messages.

Returns

bool

Examples

#
print "\t----\t calc Ex 1 \t----"
ia.maketestimage('aF', overwrite=True)
ia.calc('min(aF, (min(aF)+max(aF))/2)')
ia.calc('1.0')
ia.close()
#



The first example shows that there are 2 {\cf min} functions.  One with a
single argument returning the minimum value of that image.  The other
with 2 arguments returning an image containing ``aF'' data clipped at
the value of the 2nd argument.   The second example sets all good
pixels to unity.



#
print "\t----\t calc Ex 2 \t----"
ia.maketestimage('aD', overwrite=True)       # create some
ia.close()
ia.maketestimage('aF', overwrite=True)       # image files
ia.close()
ia.maketestimage('bF', overwrite=True)       # for use
ia.close()
ia.maketestimage('aC', overwrite=True)       # in
ia.close()
ia.maketestimage()
ia.calc('sin(aD)+(aF\*2)+min(bF)+real(aC)')   # the example
ia.close()
#


This shows a mixed type expression.  The real part of the complex image
``aC''  is used in an expression that otherwise uses Float type.
calcmask(mask='', name='', asdefault=True)[source]

This method is used to create a new pixel mask via a boolean LEL expression.

See http://casa.nrao.edu/aips2_docs/notes/223/index.shtml which describes the the syntax and functionality of LEL in detail. Also in this document is a description of ways to escape image names that contain certain non-alphanumeric characters so they are compatible with LEL syntax.

If the expression is not a scalar, the shapes and coordinates of the image and expression must conform. If the expression is a scalar then the entire pixel mask will be set to that value.

By specifying the name parameter to be an empty string, the method automatically determines the name of the output mask. Otherwise, the output mask is named the value specified by the name parameter. If a mask by that name already exists, it is overwritten. Use method ia.summary() to view existing mask names.

The asdefault parameter specifies if the new mask should be set as the default pixel mask of the image. This value is set to True by default.

Parameters

  • mask (string='') - Mask to use. Default is none.

  • name (string='') - Mask name. Default is auto new name.

  • asdefault (bool=True) - Make specified mask the default mask?

Returns

bool

Examples

#
print "\t----\t calcmask Ex 1 \t----"
ia.maketestimage('zz', overwrite=True)
subim = ia.subimage()                # Make "another" image
ia.calcmask('T')                     # Specify 'True' mask as a string
ia.calcmask('zz\>0')                  # Mask of zz ignored
ia.calcmask('mask(zz) && zz\>0')      # Mask of zz included
ia.calcmask(subim.name(True)+'\>min('+subim.name(True)+')') # Use tool names
ia.calcmask('zz\>min(zz:nomask)')  # Mask of zz not used in scalar function
subim.done()
ia.close()
#



The first calcmask example is the equivalent of {\cf
ia.set(pixelmask=1)}.  It sets the entire mask to True.

The second example creates a new \pixelmask\ which is True when
the pixel values in image {\sff zz} are greater than 0.  

Now for some subtlety.  Read carefully !  Any LEL expression can be
thought of as having a value and a mask.  Usually the value is Float and
the mask Boolean.  In this case, because the expression is Boolean
itself, the value is also Boolean.  The expression mask would just be
the mask of {\sff zz}.  Now what {\stfaf calcmask} does is create a mask
from the expression value (which is Boolean) and discards the expression
mask.  Therefore, the resulting mask is independent of any mask
that {\sff zz} might have.

If you wish the mask of the expression be honoured as well,
then you can do as in the third example.   It says the output \pixelmask\ 
will be True if the current \pixelmask\ of {\sff zz} is True and the expression
value is True.

The fourth example is like the second, except that we use the pixel
values associated with the on-the-fly {\stf subim} Image tool  disk file.  Note one further
subtlety here.  When the scalar function {\cf min} evaluates a value
from {\cf subim.name()}, which in this case is just {\cf zz}, the default
mask of {\cf subim.name()} {\it will} be used.  All the scalar
functions look at the mask.  If you didn't want the mask to be used
you can use the special {\cf :nomask} syntax shown in the final
example.
close()[source]

This function closes the . This means that it detaches the tool from its  (flushing all the changes first). The  is “null” after this change (it is not destroyed) and calling any  other than open will result in an error.

collapse(function='', axes='0', outfile='', region='', box='', chans='', stokes='', mask='', overwrite=False, stretch=False)[source]

This method collapses an image along a specified axis or set of axes of length N pixels to a single pixel on each specified axis. Both float valued and complex valued images are supported. It computes a user-specified aggregate function for pixel values along the specified axes, and places those values in the single remaining plane of those axes in the output image. The method returns an image analysis tool containing the newly-created collapsed image. Valid choices of aggregate functions are: ‘flux’ (see below for constraitns), ‘madm’ (median absolute deviation from the median), ‘max’, ‘mean’, ‘median’, ‘min’, ‘npts’, ‘rms’, ‘stddev’, ‘sum’, ‘variance’ and ‘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). Minimal unique matching is supported for the function parameter (e.g. function = ‘r’ will compute the rms of the pixel values, ‘med’ will compute the median, etc.).

If one specifies function=’flux’, the following constraints must be True:

  1. The image must have a direction coordinate,

  2. The image must have at least one beam,

  3. The specified axes must be exactly the direction coordinate axes,

  4. Only one of the non-directional axes may be non-degenerate,

  5. The iamge brightness unit must be conformant with x*yJy/beam, where x is an optional unit (such as km/s for moments images) and y is an optional SI prefix.

Axes may be specified as a single integer or an array of integers indicating the zero-based axes along which to collapse the image. Axes may also be specified as a single or array of strings which minimally and uniquely match (ignoring case) world axis names in the image (e.g. ‘dec’ for collapsing along the declination axis or [‘ri’, ‘d’] for collapsing along both the right ascension and declination axes).

If outfile is not specified (or contains only whitespace characters), no image is written but the collapsed image is still accessible via the image analysis tool this method always returns (which references the collapsed image). If the returned object is not wanted, it should still be captured and destroyed via its done() method. If this is not done, there is no guarantee as to when the Python garbage collector will delete it. If the returned object is wanted, it should still be deleted as soon as possible for the same reasons, e.g.

collapsed_image = ia.collapse(…) begin{verbatim} # do things (or not) with the collapsed_image and when finished working with the object, do end{verbatim} collapsed_image.done()

The reference pixel of the collapsed axis is set to 0 and its reference value is set to the mean of the the first and last values of that axis in the specified region of the input image. The reference value is the world coordinate value of the reference pixel. For instance, if an axis to be collapsed were to be the frequency axis, in the collapsed image, the reference value would be the mean value of the frequency range spanned, and would be stored in pixel 0.

If the input image has per plane beams, the beam at the origin of the subimage determined by the selected region is arbitrarily made the global beam of the output image. In general, the user should understand the pitfalls of collapsing images with multiple beams (i.e. that employing an aggregate function on pixels with varying beam sizes more often than not leads to ill-defined results). Convolution to a common beam is not performed automatically as part of the preprocessing before the actual rebinning occurs. In such cases, therefore, the user should probably first convolve the input image with a common restoring beam so that each plane has the same resolution, and/or use imsmooth to smooth the data to have the same beam.

Parameters

  • function (string='') - Aggregate function to apply. This can be set one of flux, madm, max, mean, median, min, npts, rms, stddev, sum, variance, xmadm. Must be specified.

  • axes (variant='0') - Zero-based axis number (specified as a list or integer) along which to collapse the specified image. Default value is 0.

  • outfile (string='') - 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 ({string, record}='') - Region selection. Default is to use the full image.

  • box (string='') - Rectangular region to select in direction plane. Default is to use the entire direction plane.

  • chans (string='') - Channels to use. Channels must be contiguous. Default is to use all channels.

  • stokes (string='') - Stokes planes to use. Planes specified must be contiguous. Default is to use all Stokes planes.

  • mask (string='') - Mask to use. Default setting is none.

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file? Ignored if “outfile” is left blank.

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default value is False.

Returns

image

Examples

# myimage.im is a 512x512x128x4 (ra,dec,freq,stokes) image
ia.open("myimage.im")
# 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.
collapsed = ia.collapse(outfile="collapse_spec_mean.im", function="mean", axes=2, box="127,127,383,383", chans="8~119")
# manipulate collapsed
collapsed.done()
commonbeam()[source]

Determine a beam to which all beams in an image can be convolved. If the image does not have a beam, an exception will be thrown. If the image has a single beam, that beam will be returned. If the image has multiple beams, this will be the beam with the largest area in the image beam set if all the other beams can be convolved to that beam. If not, this is guaranteed to be the minimum area beam to which all beams in the set can be convolved if all but one of the beams in the set can be convolved to the beam in the set with the largest area. Otherwise, the returned beam may or may not be the smallest possible beam to which all the beams in the set can be convolved.

continuumsub(outline='', outcont='continuumsub.im', region='', channels=[-1], pol='', fitorder=0, overwrite=False)[source]

This function packages the relevant image tool functionality for simple specification and application of image plane continuum subtraction. All that is required of the input image is that it have a non-degenerate spectral axis.

The user specifies region, the region of the input image over which continuum subtraction is desired (otherwise the whole image will be treated); channels, the subset of channels on the spectral axis to use in the continuum estimation, specified as a vector; fitorder, the polynomial order to use in the estimation. Optionally, output line and continuum images may be written by specifying outline and outcont, respectively. If outline is not specified, a virtual image tool is all that is produced. If outcont is not specified, the output continuum image will be written in ’continuumsub.im’. Note that the pol parameter is no longer supported; one should use the region parameter if polarization selection is desired, in conformance with other ia tool methods.

Parameters

  • outline (string='') - Output line image filename. Default is unset.

  • outcont (string='continuumsub.im') - Output continuum image filename

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • channels (intVec=[-1]) - Channels to use for continuum estimation. Default is all.

  • pol (string='') - THIS PARAMETER IS NO LONGER SUPPORTED. USE THE region PARAMETER TO CHOOSE WHICH POLARIZATIONS YOU WOULD LIKE TO PROCESS

  • fitorder (int=0) - Polynomial order for continuum estimation

  • overwrite (bool=False) - Auto-overwrite output files if they exist?

Returns

image

Examples

Fit a second order polynomial (fitorder=2) to channels 3-8 and 54-60 (python's
range function includes the lower limit and excludes the upper limit)
to an RA x Dec x Frequency x Stokes cube.

ia.open("myimage.im")

# select stokes plane 1 on which to perform the fit, as well
# as a box of pixels with blc=25,25 and trc=75,75 in the direction
# plane, and channels 0 to 100. This will be the portion of the cube
# from which the fit is subtracted
reg = rb.box(blc=[25, 25, 0, 1], trc=[75, 75, 100, 1])
csub = ia.continuumsub(region=reg, channels=range(3,9)+range(54,61), fitorder=2)

# do stuff with original image (ia) and csub image tools as necessary and
# finally close them
ia.done()
csub.done()
convertflux(value='0Jy/beam', major='1arcsec', minor='1arcsec', type='Gaussian', topeak=True, channel=-1, polarization=-1)[source]

This function interconverts between peak intensity and flux density for a Gaussian component. The image must hold a restoring beam.

Parameters

  • value (variant='0Jy/beam') - Flux density to convert. Must be specified.

  • major (variant='1arcsec') - Major axis of component. Must be specified.

  • minor (variant='1arcsec') - Minor axis of component. Must be specified.

  • type (string='Gaussian') - Type of component. String from Gaussian, Disk.

  • topeak (bool=True) - Convert to peak or integral flux desnity

  • channel (int=-1) - Channel to use if and only if image has per plane beams.

  • polarization (int=-1) - Zero-based polarization number to use for beam if and only if image has per plane beams.

Returns

record

Examples

#
print "\t----\t convertflux Ex 1 \t----"
ia.maketestimage('in.im', overwrite=True);
p1 = qa.quantity('1mJy/beam')
i1 = ia.convertflux(p1, major='30arcsec', minor='10arcsec', topeak=False);
p2 = ia.convertflux(i1, major='30arcsec', minor='10arcsec', topeak=True)
print 'peak, integral, peak = ', p1, i1, p2
#peak, integral, peak =  {'value': 1.0, 'unit': 'mJy/beam'}
#                        {'value': 0.00016396129551656742, 'unit': 'Jy'}
#                        {'value': 0.0010000000000000002, 'unit': 'Jy/beam'}

ia.close()
#
convolve(outfile='', kernel='', scale=-1.0, region='', mask='', overwrite=False, stretch=False)[source]

This function performs Fourier-based convolution of an image by the specified kernel. The input image must have real-valued pixels.

If outfile is specified, a persistent image is written to the specified disk file. If outfile is left set to the empty string, the on-the-fly image analysis tool generated by this function is associated with a temporary image. This temporary image may be stored in memory or on disk, depending on its size. When the user destroys the generated image analysis tool with the close() or done() method, the temporary image is deleted.

The kernel is provided as a multi-dimensional array or as the filename of a persistent image. The provided kernel can have fewer dimensions than the image being convolved. In this case, it will be padded with degenerate axes. An exception will be thrown if the kernel has more dimensions than the image. No additional scaling of the kernel is provided.

The scaling of the output image is determined by the value of the scale parameter. If this is left unset, then the kernel is normalized to unit sum. If scale is not left unset, then the convolution kernel will be scaled (multiplied) by this value.

Masked pixels will be assigned the value 0.0 before convolution.

The output mask is the logical AND of the input image’s default pixel mask (if any) and the OTF mask. Any other pixel masks associated with the input image will not be copied. The maskhandler method should be used if there is a need to copy other masks too.

See also the other convolution methods convolve2d(), sepconvolve(), and hanning() for more specialized convolution algorithms.

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • kernel (variant='') - Convolution kernel - An array or an image filename. Must be specified by the user.

  • scale (double=-1.0) - Scale factor. The default behavior is to autoscale (specified with -1.0).

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. The default setting is none.

  • overwrite (bool=False) - Overwrite (unprompted) the pre-existing output file?

  • stretch (bool=False) - Stretch the mask if necessary and possible?

Returns

image

Examples

#
print "\t----\t convolve Ex 1 \t----"
# This example presupposes the existence of an input image file, testdata, and a kernel image file, kerneldata.
# Open the input image file:
ia.open(infile='testdata')
# Set up a region to be operated upon (in this case, the whole image):
r1 = rg.box()
# Perform the convolution:
im2 = ia.convolve (outfile = 'convout', overwrite = True, region = r1, kernel = 'kerneldata')
ia.close()
im2.done()

#
print "\t----\t convolve Ex 2 \t----"
# This example uses an array as the convolution kernel, and presupposes the existence of an input image file, testdata, which we first open:
ia.open(infile='testdata')
# Next, create a Python array of some kind to use as a convolution kernel. For example:
from numpy import arange
kernelarray = arange(10)\*\*3
# Set up a region to be operated upon (in this case, the whole image):
r1 = rg.box()
# Perform the convolution:
im3 = ia.convolve (outfile = 'convout2', overwrite = True, region = r1, kernel = kernelarray)
ia.close()
im3.done()
#
convolve2d(outfile='', axes=[0, 1], type='gaussian', major='0deg', minor='0deg', pa='0deg', scale=-1, region='', mask='', overwrite=False, stretch=False, targetres=False, beam='')[source]

This function performs Fourier-based convolution of an  using the provided 2D kernel.

If outfile is left unset, the image is written to the specified disk file. If outfile is not given, the newly constructed on-the-fly Image  is associated with a temporary image. This temporary image may be stored in memory or on disk, depending on its size. When the user destroys the on-the-fly Image  (with the done function) this temporary image is deleted.

The user specifies which 2 pixel axes of the image are to be convolved via the axes argument. The pixels must be square or an error will result.

The user specifies the type of convolution kernel with type (minimum match is supported); currently only ’gaussian’ is available.

The user specifies the parameters of the convolution kernel via the arguments major, minor, and pa. These arguments can be specified in one of three ways:

  • Quantity - for example major=qa.quantity(1, ’arcsec’) Note that you pixel units can be used, viz. major=qa.quantity(1, ’pix’), see below.

  • String - for example minor=’1km’ (i.e. one that the Quanta quantity function accepts).

  • Numeric - for example major=10. In this case, the units of major and minor are assumed to be in pixels. Using pixel units allows the user to convolve unlike axes (see one of the provided example for this use case). For the position angle, units of degrees are assumed.

The interpretation of major and minor depends upon the kernel type.

  • Gaussian - major and minor are the Full Width at Half Maximum (FWHM) of the major and minor axes of the Gaussian.

The position angle is measured North through East when a plane holding a celestial coordinate (the usual astronomical convention) is convolved. For other axis/coordinate combinations, a positive position angle is measured from +x to +y in the absolute pixel coordinate frame (x is the first axis that is specified, with argument axes).

In the case of a Gaussian, the beam parameter offers an alternate way of describing the convolving Gaussian. If used, neither major, stfaf minor, nor pa can be specified. The beam parameter must have exactly three fields: “major”, “minor”, and “pa” (or “positionangle”). This is, not coincidentally, the record format for the output of ia.restoringbeam().

The scaling of the output image is determined by the argument scale. If this is left unset then autoscaling will be invoked.

If the user is not convolving the sky, then autoscaling means that the convolution kernel will be normalized to have unit volume so as to conserve flux.

If the user is convolving the sky, then there are two cases for which autoscaling is useful:

Firstly, if the input image units are Jy/pixel, then the output image will have units of Jy/beam and be appropriately scaled. In addition, the restoring beam of the output image will be the same as the convolution kernel.

Secondly,if the input image units are Jy/beam, then the output image will also have units of Jy/beam and be appropriately scaled. In addition, the restoring beam of the output image will be the convolution of the input image restoring beam and the convolution kernel. In the case of an image with per-plane beams, for each plane, the kernel is convolved with the appropriate beam and the result is associated with that plane in the output image.

If the user sets a value for scale, then the convolution kernel will be scaled by this value. Note that it has peak of unity before the application of this scale factor.

If the units on the original image include Jy/beam, the units on the output image will be rescaled by the ratio of the input and output beams as well as rescaling by the area of convolution kernel.

If the units on the original image include K, then only the image convolution kernel rescaling is done.

If targetres=True and type=”gaussian” and the input image has a restoring beam, this method will interpret the values of major, minor, and pa as the resolution of the final image and will calculate the parameters of the Gaussian to use in the convolution so that this target resolution is achieved.

Masked pixels will be assigned the value 0.0 before convolution. The output mask is the intersection (logical AND) of the default input  (if any) and the OTF mask. Any other input  will not be copied. The function maskhandler can be used if there is a need to copy other masks too.

See also the other convolution functions:

convolve, hanning, and sepconvolve.

Parameters

  • outfile (string='') - Output image file name. The default value is unset.

  • axes (intVec=[0, 1]) - Axes to convolve. The default setting is [0,1].

  • type (string='gaussian') - Type of convolution kernel to be used.

  • major (variant='0deg') - Major axis, Quantity, string, numeric (e.g. 10arcsec, 20pix, 3km, etc.). This must be specified by the user.

  • minor (variant='0deg') - Minor axis, Quantity, string, numeric (e.g. 10arcsec, 20pix, 3km, etc.). This must be specified by the user.

  • pa (variant='0deg') - Position Angle, Quantity, string, numeric (e.g. 10deg). The default value is 0deg.

  • scale (double=-1) - Scale factor. The default setting (-1) is to autoscale.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. The default option is none.

  • overwrite (bool=False) - Overwrite (unprompted) the pre-existing output file?

  • stretch (bool=False) - Stretch the mask if necessary and possible?

  • targetres (bool=False) - If True and type=”gaussian”, major, minor, and pa are interpreted as the image resolution that the user wants to achieve.

  • beam (record='') - Alternate way of describing a Gaussian. Must have fields “major”, “minor”, and “pa” (or “positionangle”)

Returns

image

Examples

#
print "\t----\t convolve2d Ex 1 \t----"
ia.maketestimage('xy',overwrite=True)         # Create a simple RA/DEC test image
# Convolve axes 0 and 1 of the test image with a 20x10-arcsec, 45-degree Gaussian:
im2 = ia.convolve2d(outfile='xy.con', axes=[0,1], type='gauss',
                    major='20arcsec', minor='10arcsec', pa='45deg',
                    overwrite=True);
# Clean up, by destroying the im2 tool and close the image tool:
im2.done()
ia.close()
#

ia.fromarray(outfile='xypf', pixels=ia.makearray(0, [64, 64, 4, 64]),
             overwrite=True)         # Create a simple RA/DEC/Pol/Freq test dataset
print "!!!EXPECT WARNING REGARDING INVALID SPATIAL RESTORING BEAM!!!"
# Convolve axes 0 and 3 of the test dataset with a 20x10-pixel, 45-degree Gaussian:
im2 = ia.convolve2d(outfile='xypf.con', axes=[0,3], type='gauss',
                    major='20pix', minor='10pix', pa='45deg',
                    overwrite=True);
# Note that pixel units must be used in the above because axes 0 and 3 are unlike.
# Clean up, by destroying the im2 tool and close the image tool:
im2.done()
ia.close()
#
coordmeasures(pixel=[-1], dframe='cl', sframe='cl')[source]

You can use this function to get the world coordinates for a specified absolute pixel coordinate in the image. You specify a pixel coordinate (0-rel) for each axis in the image.

If you supply fewer pixel values then there are axes in the image, your value will be padded out with the reference pixel for the missing axes. Excess values will be ignored.

The parameters dframe and sframe allow one to specify to which reference frame the direction and spectral measures, respectively, should be converted. These values are case-insensitive. “native” means use the native reference frame of the coordinate in question. “cl” means use the conversion layer frame if one exists (if not, the native frame will be used).

The world coordinate is returned as a record of measures. This function is just a wrapper for the Coordsys tool toworld function (invoked with argument format=’m’). Please see its documentation for discussion about the formatting and meaning of the measures.

This Image  function adds two additional fields to the return record.

The mask field contains the value of the image  at the specified position. It is either T (pixel is good) or F (pixel is masked as bad or the specified position was off the image).

The intensity field contains the value of the image (at the nearest pixel to that given) and its units. This is actually stored as a Quantity. This field does not exist if the specified pixel coordinate is off the image.

Parameters

  • pixel (doubleVec=[-1]) - Absolute pixel coordinate. Default is reference pixel.

  • dframe (string='cl') - Direction reference frame to which to convert the direction data. Case insensitive. “cl” means use the conversion layer, if present, of the image direction coordinate. “native” means use the native native direction frame of the image. Other examples are “J2000”, “B1950”, “GALACTIC”, etc.

  • sframe (string='cl') - Spectral reference frame to which to convert the spectral data. Case insensitive. “cl” means use the conversion layer, if present, of the image spectral coordinate. “native” means use the native spectral reference frame of the image. Other examples are “LSRK”, “CMB”, “LGROUP”, etc.

Returns

record

Examples

#
print "\t----\t coordmeasures Ex 1 \t----"
ia.maketestimage('myimage',overwrite=True)
s = ia.shape()
for i in range(len(s)):
  s[i] = 0.5\*s[i]
meas = ia.coordmeasures(s)
print meas.keys()                   # Get names of fields in record
#['intensity', 'mask', 'measure']
print meas['intensity']
#{'value': 1.39924156665802, 'unit': 'Jy/beam'}
print meas['measure']['direction']
#{'type': 'direction',
# 'm1': {'value': 5.817764003289323e-05, 'unit': 'rad'},
# 'm0': {'value': -5.8177644130875234e-05, 'unit': 'rad'}, 'refer': 'J2000'}
dir = meas['measure']['direction']  # Get direction coordinate
me.doframe(me.observatory('ATCA'))  # Set location on earth
me.doframe(me.epoch('utc','16jun1999/12:30:20'))  # Set epoch
azel = me.measure(dir,'azel')       # Convert to azimuth/elevation
print 'az,el=', qa.angle(azel['m0']), qa.angle(azel['m1'])  # Format nicely
#az,el= +105.15.47 -024.22.57
meas2=ia.coordmeasures()            # defaults to reference pixel
print meas2['intensity']
#{'value': 2.5064315795898438, 'unit': 'Jy/beam'}
print meas2['measure']['direction']
#{'type': 'direction',
# 'm1': {'value': 0.0, 'unit': 'rad'},
# 'm0': {'value': 0.0, 'unit': 'rad'}, 'refer': 'J2000'}
dir = meas2['measure']['direction'] # Get direction coordinate
me.doframe(me.observatory('ATCA'))  # Set location on earth
me.doframe(me.epoch('utc','16jun1999/12:30:20'))   # Set epoch
azel = me.measure(dir,'azel')       # Convert to azimuth/elevation
print 'az,el=', qa.angle(azel['m0']), qa.angle(azel['m1'])
#az,el= +105.16.05 -024.23.00
#



In this example we first find the world coordinates of the centre of the
image.  Then we use the Measures \tool\ {\stf me} to convert the
{\cf direction coordinate} field from J2000 to an azimuth and elevation
at a particular location at a particular time. 
coordsys(axes=[-1])[source]

This function returns the Coordinate System of an image in a Coordsys tool. Both float and complex valued images are supported.

By default, the Coordinate System describes all of the axes in the image. If you desire, you can select a subset of the axes, thus reducing the dimensionality of the Coordinate System. This may be useful if you are supplying a Coordinate System to the functions fromarray or fromshape.

Parameters

  • axes (intVec=[-1]) - Axes to which the Coordinate System pertains. Default is all axes.

Returns

coordsys

Examples

#
print "\t----\t coordsys Ex 1 \t----"
ia.maketestimage('hcn',overwrite=True)
ia.summary()
mycs = ia.coordsys([0,1])
imshape = ia.shape()
ia.fromshape(outfile='test', shape=imshape, csys=mycs.torecord(), overwrite=True)
ia.summary()
mycs.done()
ia.close()
#



In this example, we create a Coordinate System pertaining to the first
two axes of the image and then we create a new (empty) 2D image with
this Coordinate System using the {\cf fromshape} function.
crop(outfile='', axes='', overwrite=False, region='', box='', chans='', stokes='', mask='', stretch=False, wantreturn=True)[source]

This method crops masked slices from the perimeter of an image. The axes parameter specifies which axes to consider. Axes not specified will not be cropped. An empty array implies that all axes should be considered. If wantreturn is True, an image analysis tool attached to the output image is returned. If False, none is returned.

Parameters

  • outfile (string='') - Output image name. If not specified, no persistent image is created.

  • axes (intVec='') - Axes to crop. Empty array means consider all axes.

  • overwrite (bool=False) - Overwrite the output if it exists? Default False

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • box (string='') - Rectangular region to select in direction plane. Default is to use the entire direction plane.

  • chans (string='') - Channels to use. Default is to use all channels.

  • stokes (string='') - Polarization selection. Default is all.

  • mask (string='') - Mask to use. Default is none.

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

  • wantreturn (bool=True) - Return an image analysis tool attached to the created subimage?

Returns

image

Examples

# myimage is of shape 20, 20, 20 with only the inner 16 x 14 x 12 pixels unmasked
ia.open("myimage")
# crop masked slices on all axes
cropped = ia.crop()
# returns [16, 14, 12]
cropped.shape()
cropped.done()
# crop only the masked slices at the edges of the image along axis 1
cropped2 = ia.crop(outfile="", axes=[1])
ia.done()
# returns [20, 14, 20]
cropped2.shape()
cropped2.done()
decimate(outfile='', axis=0, factor=1, method='copy', region='', mask='', overwrite=False, stretch=False)[source]

This application removes planes along the specified axis of an image. It supports both float valued and complex valued images. The factor parameter represents the factor by which to reduce the number of planes.

The method parameter represents how to calculate the pixel values of the output image. A value of method=”copy” means that every factorth plane of the selected region in the input image will be directly copied to the corresponding plane in the output image. So, if one wanted to copy every third spectral plane in the input image to the output image, one would specify factor=3 and method=”copy”. If the selected region along the specified axis had 11 planes, then there would be 4 output planes which would map to planes 0, 3, 6, and 9 of the specified region of input image. A value of method=”mean” indicates that each of factor number of planes in the range starting at each factorth plane should be averaged to produce the corresponding output plane. So, if one specified factor=3 and method=”mean” along an axis of the selected region of the input image which had 11 pixels, the corresponding axis in the output image would have three pixels and the pixel values for each of those output planes would corresponding to averaging along that axis planes 0-2, 3-5, and 6-8 of the selected region of the input image. Note that the remaining planes, 9 and 10, in the selected region of the input image would be ignored because the last interval must have exactly factor number of planes in order to be included in the output image.

The coordinate system of the output image takes into account the decimation; that is, along the decimated axis, the increment of the output image is factor times that of the input image, and the reference pixel of the output image is located at pixel 1/factor times the reference pixel in the input image.

This method returns an image analysis tool which references the output image. If this tool is not desired, one should capture it anyway and then close() it immediately to free up resources.

Images with multiple beams are not supported; please convolve a multi-beam image to a single resolution before running this application.

Parameters

  • outfile (string='') - Output image file name. If empty, a persistent image is not created.

  • axis (int=0) - Axis along which to remove planes.

  • factor (int=1) - Reduce number of planes by this factor.

  • method (string='copy') - Method to use for calculating pixel values of output. Supported values are “copy” or “mean”.

  • region ({string, record}='') - Region selection. Default is to use the full image.

  • mask (string='') - Mask to use. Default setting is none.

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file? Ignored if “outfile” is left blank.

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default value is False.

Returns

image

Examples

# Copy verbatim every 5th plane of axis 2 of the input image
ia.open("myim.im")
decimated = ia.decimate("dec1.im", axis=2, factor=5, method="copy")
# do stuff with decimated and then close it
decimated.close()

# Decimate by averaging every 7 planes of the input image along axis 2
decimated = ia.decimate("dec2.im", axis=2, factor=7, method="mean")
# do stuff with decimated and then close it
decimated.close()
# close input ia tool
ia.close()
decompose(region='', mask='', simple=False, threshold=-1, ncontour=11, minrange=1, naxis=2, fit=True, maxrms=-1, maxretry=-1, maxiter=256, convcriteria=0.0001, stretch=False)[source]

This function is an image decomposition tool that performs several tasks, with the end result being that a strongly blended image is separated into components - both in the sense that it determines the parameters for each component (assuming a Gaussian model) and that it physically assigns each pixel in the image to an individual object. The products of these two operations are called the component list and the component map, respectively. The fitting process (which determines the component list) and the pixel-decomposition process (which determines the component map) are designed to work cooperatively to increase the efficiency and accuracy of both.

The algorithm behind the decomposition is based on the function clfind, described in Williams et al 1994, which uses a contouring procedure whereby a closed contour designates a separate component. The program first separates the image into clearly distint ’regions’ of blended emission, then contours each region to determine the areas constituting each component and passes this information on to the fitter, which determines the component list.

The contour deblending can optionally be replaced with a simpler local maximum scan, and the fitting can be replaced with a moment-based estimation method to speed up calculations on very large images or if either primary method causes trouble, but in general this will impede the accuracy of the fit.

The function works with both two and three dimensional images.

The return value is a record (or dictionary) that has 3 keys: ’components’, ’blc’, ’trc’. The ’components’ element is a matrix each row of which contains the gaussian parameters of the component fitted. The ’blc’ element is a matrix of the bottom left corners (blc) of the regions found. Each row correspond to a region blc. The ’trc’ element is a matrix of the top right corners (trc) of the regions found. Each row correspond to a region trc. Please Note that the returned blc’s and trc’s are relative to region defined by the user. A blc of [0,0] implies the bottom left of the region selected and not the bottom left of the image. Obviously if no region is defined then it is the bottom left of the image.

Parameters

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • simple (bool=False) - Skip contour deblending and scan for local maxima

  • threshold (double=-1) - Value of minimum positive contour. Must be set and nonnegative.

  • ncontour (int=11) - Number of contours to use in deblending (>= 2)

  • minrange (int=1) - Minimum number of closed contours in a component (> 0)

  • naxis (int=2) - Max number of perpendicular steps between contiguous pixels. Values of 1, 2 or 3 are allowed.

  • fit (bool=True) - Fit to the components after deblending?

  • maxrms (double=-1) - Maximum RMS of fit residuals to not retry fit (> 0). Default is unset.

  • maxretry (int=-1) - Maximum number of times to retry the fit (>= 0). Default is unset.

  • maxiter (int=256) - Maximum number of iterations allowed in a single fit (> 0)

  • convcriteria (double=0.0001) - Criterion to establish convergence (>=0)

  • stretch (bool=False) - Stretch the mask if necessary and possible?

Returns

record

Examples

#
print "\t----\t decompose Ex 1 \t----"
ia.maketestimage()
out=ia.decompose(threshold=2.5, maxrms=1.0)
#Attempt 1: Converged after 21 iterations
#Attempt 1: Converged after 15 iterations
#1: Peak: 17.955  Mu: [0.000327928, 8.62573e-05]
#               Axes: [0.00175981, 0.00142841]  Rotation: 1.29539
#2: Peak: 19.8093  Mu: [1.67927e-06, -0.000374393]
#                Axes: [0.00179054, 0.00132541]  Rotation: 1.78404
#3: Peak: 10.1155  Mu: [6.28252, -7.09688e-05]
#                Axes: [0.00180877, 0.00104523]  Rotation: 1.78847
print out['components']
#[[  1.79549522e+01   3.27928370e-04   8.62573434e-05   1.75980886e-03
#    8.11686337e-01   1.29538655e+00]
# [  1.98093319e+01   1.67927124e-06  -3.74393392e-04   1.79054437e-03
#    7.40229547e-01   1.78403902e+00]
# [  1.01155214e+01   6.28252172e+00  -7.09688029e-05   1.80877140e-03
#    5.77867746e-01   1.78847444e+00]]
print out['blc']
#[[37 31]
# [47 25]
# [67 33]]
print out['trc']
#[[54 47]
# [66 38]
# [78 40]]
ia.close()
#
deconvolvecomponentlist(complist='', channel=-1, polarization=-1)[source]

This method deconvolves (a record representation of) a Componentlist tool from the restoring beam, returning (a record representation of) a new Componentlist tool. If there is no restoring beam, a fail is generated.

Currently, only deconvolution of Gaussian components is supported.

For images with per-plane beam, the user must choose which beam is used for the deconvolution by setting channel and/or polarization. Only a single beam is used to deconvolve all components.

See also functions setrestoringbeam and restoringbeam.

Parameters

  • complist (record='') - Componentlist to deconvolve

  • channel (int=-1) - Zero-based channel number to use for beam for per plane images. Not used if the image has a single beam.

  • polarization (int=-1) - Zero-based polarization number to use for beam for per plane images. Not used if the image has a single beam.

Returns

record

Examples

#
print "\t----\t deconvolvecomponentlist Ex 1 \t----"
ia.maketestimage()
r = ia.fitcomponents()
cl1 = r['results']                     # cl1 and cl2 are record representations
r = ia.fitcomponents()
cl1 = r['results']                      # cl1 and cl2 are record representations
cl2 = ia.deconvolvecomponentlist(cl1)  #   of componentlists
print cl1, cl2
cl.fromrecord(cl2)                     # set componentlist tool with record
ia.close()
cl.close()
#
deconvolvefrombeam(source='', beam='')[source]

This is a helper function. It is to provide a way to deconvolve gaussians from other gaussians if that is what is needed for example removing a beam Gaussian from a Gaussian source. To run this function the tool need not be attached to an image.

The return value is a record that contains the fit param and the return value is a boolean which is set to True if fit model is a point source

Parameters

  • source (variant='') - Three quantities that define the source majoraxis, minoraxis and Position angle

  • beam (variant='') - Three quantities that define the beam majoraxis, minoraxis and Position angle

Returns

record

Examples

#
print "\t----\t deconvolvefrombeam Ex 1 \t----"
ia.maketestimage()
recout=ia.deconvolvefrombeam(source=['5arcmin', '3arcmin', '20.0deg'], beam=['50arcsec','30arcsec', '15deg'])
ia.close()
print 'Is pointsource ', recout['return']
print 'major=',recout['fit']['major']
print 'minor=',recout['fit']['minor']
print 'pa=',recout['fit']['pa']
deviation(outfile='', region='', mask='', overwrite=False, stretch=False, grid=[1, 1], anchor='ref', xlength='1pix', ylength='1pix', interp='cubic', stattype='sigma', statalg='classic', zscore=-1, maxiter=-1)[source]

This application creates an image that reflects the statistics of the input image. The output image has the same dimensions and coordinate system as the (selected region in) input image. The grid parameter describes how many pixels apart “grid” pixels are. Statistics are computed around each grid pixel. Grid pixels are limited to the direction plane only; independent statistics are computed for each direction plane (ie at each frequency/stokes pixel should the input image happen to have such additional axes). Using the xlength and ylength parameters, one may specify either a rectangular or circular region around each grid point that defines which surrounding pixels are used in the statistic computation for individual grid points. If the ylength parameter is the empty string, then a circle of diameter provided by xlength centered on the grid point is used. If ylength is not empty, then a rectangular box of dimensions xlength x ylength centered on the grid pixel is used. These two parameters may be specified in pixels, using either numerical values or valid quantities with “pix” as the unit (eg “4pix”). Otherwise, they must be specified as valid angular quantities, with recognized units (eg “4arcsec”). As with other region selections in CASA, full pixels are included in the computation even if the specified region includes only a fraction of that pixel. BEWARE OF MACHINE PRECISION ISSUES, because you may get a smaller number of pixels included in a region than you expect if you specify, eg, an integer number of pixels. In such cases, you probably want to specify that number plus a small epsilon value (eg “2.0001pix” rather than “2pix”) to mitigate machine precision issues when computing region extents.

The output image is formed by putting the statistics calculated at each grid point at the corresponding grid point in the output image. Interpolation of these output values is then used to compute values at non-grid-point pixels. The user may specify which interpolation algorithm to use for this computation using the interp parameter.

The input image pixel mask is copied to the output image. If interpolation is performed, output pixels are masked where the interpolation fails.

ANCHORING THE GRID

The user may choose at which pixel to “anchor” the grid. For example, if one specifies grid=[4,4] and anchor=[0,0], grid points will be located at pixels [0,0], [0,4], [0,8] … [4,0], [4,4], etc. This is exactly the same grid that would be produced if the user specified anchor=[4,4] or anchor=[20,44]. If the user specifies anchor=[1, 2] and grid=[4,4], then the grid points will be at pixels [1,2], [5,2], [9,2]… [5,2], [5,6], etc. and the resulting grid is the same as it would be if the user specified eg anchor=[9,10] or anchor=[21, 18]. The value “ref”, which is the default, indicates that the reference pixel of the input image should be used to anchor the grid. The x and y values of this pixel will be rounded to the nearest integer if necessary.

SUPPORTED STATISTICS AND STATISTICS ALGORITHMS

One may specify which statistic should be represented using the stattype parameter. The following values are recognized (minimum match supported):

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

Using the statalg parameter, one may also select whether to use the Classical or Chauvenet/ZScore statistics algorithm to compute the desired statistic (see the help for ia.statistics() or imstat for a full description of these algorithms).

Parameters

  • outfile (string='') - 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 ({string, record}='') - Region selection. Default is to use the full image.

  • mask (string='') - Mask to use. Default setting is none.

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file? Ignored if “outfile” is left blank.

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default value is False.

  • grid (intVec=[1, 1]) - x,y grid spacing. Array of exactly two positive integers.

  • anchor ({string, intVec}='ref') - x,y anchor pixel location. Either “ref” to use the image reference pixel, or an array of exactly two integers.

  • xlength ({string, int}='1pix') - Either x coordinate length of box, or diameter of circle. Circle is used if ylength is empty string.

  • ylength ({string, int}='1pix') - y coordinate length of box. Use a circle if ylength is empty string.

  • interp (string='cubic') - Interpolation algorithm to use. One of “nearest”, “linear”, or “cubic”. Minimum match supported.

  • stattype (string='sigma') - Statistic to compute. See full description for supported statistics.

  • statalg (string='classic') - Statistics computation algorithm to use. Supported values are “chauvenet” and “classic”, Minimum match is supported.

  • zscore (double=-1) - For chauvenet, this is the target maximum number of standard deviations data may have to be included. If negative, use Chauvenet’s criterion. Ignored if algorithm is not “chauvenet”.

  • maxiter (int=-1) - For chauvenet, this is the maximum number of iterations to attempt. Iterating will stop when either this limit is reached, or the zscore criterion is met. If negative, iterate until the zscore criterion is met. Ignored if algortihm is not “chauvenet”.

Returns

image

Examples

# compute standard deviations in circles of diameter 10arcsec around
# grid pixels spaced every 4 x 5 pixels and anchored at pixel [30, 40],
# and use linear interpolation to compute values at non-grid-pixels
ia.open("my.im")
zz = ia.deviation("sigma.im", grid=[4, 5], anchor=[30, 40], xlength="10arcsec", stattype="sigma", interp="lin", statalg="cl")

# compute median of the absolute deviations from the median values using
# the z-score/Chauvenet algorithm, by fixing the maximum z-score to determine outliers to 5.
# Use cubic interpolation to compute values for non-grid-point pixels. Use a rectangular region
# of dimensions 5arcsec x 20arcsec centered on each grid point as the region in which to include
# pixels for the computation of stats for that grid point.
zz = ia.deviation("madm.im", grid=[4, 5], anchor=[30, 40], xlength="5arcsec", ylength="20arcsec, stattype="madm", interp="cub", statalg="ch", zscore=5)
dohistory(enable=True)[source]

This allows control over if tool methods record history of what parameters they were called with to the input and/or output image. By default, tool methods will write to the image history. By explicitly disabling history writing, tool methods will not write to the history. When created, an ia tool will have history writing enabled. Note that the setting is specific to the individual tool, and that methods such as open(), close(), done(), fromshape(), etc do not implicitly change the internal state of whether or not history writing by methods is enabled. One can explicitly enable/disable history writing even if the tool is not yet attached to an image. In the case where a method returns a new ia tool, the method will not write history to the output image, but the returned tool will have history writing enabled, so that running methods on the returned tool will cause history to be written to the image attached to that tool, or any new image created by running methods on that tool, unless dohistory(False) is explicitly run on that tool prior to running other methods.

IMPORTANT NOTE: This setting will not affect the behavior of ia.sethistory(); this tool method will always write to the history, no matter if ia.dohistory(False) was run prior.

Parameters

  • enable (bool=True) - Enable history writing of tool methods? True means yes, False means no.

Returns

bool

done(remove=False, verbose=True)[source]

When the user no longer needs to use an , calling this function will free up its resources. That is, it destroys the . This means that the user can no longer call any functions on the  after it has been done.

If the Image  is associated with a disk file, then (unlike the close function, the user can also choose to delete that by setting remove=True. By default, any associated disk file is not deleted.

Note that this function is different from the close function because the latter does not destroy the . For example, the user can use the open function straight after the close function on the same .

Parameters

  • remove (bool=False) - Delete the associated disk file as well?

  • verbose (bool=True) - Send a progress report to the logger?

Returns

bool

Examples

#
print "\t----\t done Ex 1 \t----"
# Make a test image and create tool subim:
ia.maketestimage('myfile',overwrite=True)
subim = ia.subimage('myfile2',overwrite=True)
# Check that subim exists as intended by attempting to display its summary:
subim.summary()      # This displays a summary of the dataset.
# Use done to destroy the subim tool:
subim.done()
# Check that the subim tool has been detached as intended, by attempting to display its summary:
subim.summary()      # This should now throw an error.
ia.summary()     # This still works, though, as the ia tool is still open, and the dataset is still available.
ia.close()
#
fft(real='', imag='', amp='', phase='', axes=[-1], region='', mask='', stretch=False, complex='')[source]

This method fast Fourier Transforms the supplied image to the Fourier plane. If the axes parameter is left unset, then the direction plane of the image (if there is one) is transformed. Otherwise, the user can specify which axes are to be transformed. Note that if the direction plane is to be transformed, both axes associated with it must be specified.

The user specifies which form is desired in the result by specifying the desired output image file name(s).

Before the FFT is performed, any masked pixels are set to values of zero. The output mask is the result of ANDing the default input pixel mask (if any) and the OTF mask. Any other input pixel masks will not be copied. The method maskhandler() can be used if there is a need to copy other masks.

The following rules are used to set the brightness units of the output images. 1. The phase image always has units of radians. For the other output images, 2. if the input image has units of Jy/beam or Jy/pixel (ie it is an image-plane image), the output (uv-plane) images will have units of Jy. In the case of the input image having a synthesized beam, the beam will be copied to the output images (which is important for transforming back). 3. If the input image has units of Jy (ie is a uv-plane image), the output images will have either units of Jy/beam or Jy/pixel, depending on if the input image has a beam or not.

For some transformations (e.g., UV domain to image domain transforms), it is not possible to automatically generate an expected coordinate system for the output image(s); only the FFT numerics are performed and the coordinate system is generated using generic conventions.

Parameters

  • real (string='') - Output real image file name.

  • imag (string='') - Output imaginary image file name.

  • amp (string='') - Output amplitude image file name.

  • phase (string='') - Output phase image file name.

  • axes (intVec=[-1]) - Specify the pixel axes that are to undergo the FFT. The default option (-1) is to transform the sky plane(s).

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - The mask to be used. The default option is none.

  • stretch (bool=False) - Stretch the mask if it is necessary and possible.

  • complex (string='') - Output complex valued image file name.

Returns

bool

Examples

#
print "\t----\t fft Ex 1 \t----"
# Create a test image:
ia.maketestimage('gc.small', overwrite=True)
# Perform an FFT on the sky plane of the test image,
# writing out just the resulting real and amplitude images:
ia.fft(real='r.im', amp='a.im')
# Close the image tool when done:
ia.close()
# Lastly, clean up the example output files:
ia.removefile('r.im')
ia.removefile('a.im')
#
findsources(nmax=20, cutoff=0.1, region='', mask='', point=True, width=5, negfind=False)[source]

This function finds strong point sources in the image. The sources are returned in a record that can be used by a Componentlist .

An efficient method is used to locate sources under the assumption that they are point-like and not too close to the noise. Only sources with a peak greater than the cutoff fraction of the strongest source will be found. Only positive sources will be found, unless the negfind=T whereupon positive and negative sources will be found.

After the list of point sources has been made, you may choose to make a Gaussian fit for each one (point=F) so that shape information can be recovered as well. You can specify the half-width of the fitting grid with argument width which defaults to 5 (fitting grid would then be [11,11] pixels). If you set width=0, this is a signal that you would still like Gaussian components returned, but a default width should be used for the Gaussian shapes. The default is such that the component is circular with a FWHM of width pixels.

Thus, if point=T, the components in the returned Componentlist are Point components. If point=F then Gaussian components are returned.

The  must be 2-dimensional and it must hold a region of the sky. Any degenerate trailing dimensions in the region are discarded.

See also the function fitcomponents (for which findsources can provide an initial estimate).

Parameters

  • nmax (int=20) - Maximum number of sources to find, > 0

  • cutoff (double=0.1) - Fractional cutoff level

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • point (bool=True) - Find only point sources?

  • width (int=5) - Half-width of fit grid when point=F

  • negfind (bool=False) - Find negative sources as well as positive?

Returns

record

Examples

#
print "\t----\t findsources Ex 1 \t----"
ia.maketestimage()
clrec = ia.findsources(nmax=5, cutoff=0.5)
print clrec
#



All sources stronger than 0.5 of the strongest will be found.
We use the Componentlist GUI to look at the strongest component.
fitcomponents(box='', region='', chans='', stokes='', mask='', includepix=[-1], excludepix=[-1], residual='', model='', estimates='', logfile='', append=True, newestimates='', complist='', overwrite=False, dooff=False, offset=0.0, fixoffset=False, stretch=False, rms='', noisefwhm='', summary='')[source]

OVERVIEW

This application is used to fit one or more two dimensional gaussians to sources in an image as well as an optional zero-level offset. Fitting is limited to a single polarization but can be performed over several contiguous spectral channels. If the image has a clean beam, the report and returned dictionary will contain both the convolved and the deconvolved fit results.

When dooff is False, the method returns a dictionary with keys named ’converged’, ’pixelsperarcsec’, ’results’, and ’deconvolved’. The value of ’converged’ is a boolean array which indicates if the fit converged on a channel by channel basis. The value of ’pixelsperarcsec’ is a two element double array with the absolute values of the direction coordinate pixel increments (longitude-like and latitude-like coordinate, respectively) in arcsec. The value of ’results’ is a dictionary representing a component list reflecting the fit results. In the case of an image containing beam information, the sizes and position angles in the ’results’ dictionary are those of the source(s) convolved with the restoring beam, while the same parameters in the ’deconvolved’ dictionary represent the source sizes deconvolved from the beam. In the case where the image does not contain a beam, ’deconvolved’ will be absent. Both the ’results’ and ’deconvolved’ dictionaries can be read into a component list tool (default tool is named cl) using the fromrecord() method for easier inspection using tool methods, eg

cl.fromrecord(res[’results’])

although this only works if the flux density units are conformant with Jy.

There are also values in each component subdictionary not used by cl.fromrecord() but meant to supply additional information. There is a ’peak’ subdictionary for each component that provides the peak intensity of the component. It is present for both ’results’ and ’deconvolved’ components. There is also a ’sum’ subdictionary for each component indicated the simple sum of pixel values in the the original image enclosed by the fitted ellipse. There is a ’channel’ entry in the ’spectrum’ subdictionary which provides the zero-based channel number in the input image for which the solution applies. In addtion, if the image has a beam(s), then there will be a ’beam’ subdictionary associated with each component in both the ’results’ and ’deconvolved’ dictionaries. This subdictionary will have three keys: ’beamarcsec’ will be a subdictionary giving the beam dimensions in arcsec, ’beampixels’ will have the value of the beam area expressed in pixels, and ’beamster’ will have the value of the beam area epressed in steradians. Also, if the image has a beam(s), in the component level dictionaries will be an ’ispoint’ entry with an associated boolean value describing if the component is consistent with a point source. Each component level dictionary will have a ’pixelcoords’ entry which has the value of a two element numeric array which provides the direction pixel coordinates of the fitted position.

If dooff is True, in addtion to the specified number of gaussians, a zero-level offset will also be fit. The initial estimate for this offset is specified using the offset parameter. Units are assumed to be the same as the image brightness units. The zero level offset can be held constant during the fit by specifying fixoffset=True. In the case of dooff=True, the returned dictionary contains two additional keys, ’zerooff’ and ’zeroofferr’, which are both dictionaries containing ’unit’ and ’value’ keys. The values associated with the ’value’ keys are arrays containing the the fitted zero level offset value and its error, respectively, for each channel. In cases where the fit did not converge, these values are set to NaN. The value associated with ’unit’ is just the image brightness unit.

The region can either be specified by a box(es) or a region. Ranges of pixel values can be included or excluded from the fit. If specified using the box parameter, multiple boxes can be given using the format box=”blcx1, blcy1, trcx1, trcy1, blcx2, blcy2, trcx2, trcy2, … , blcxN, blcyN, trcxN, trcyN” where N is the number of boxes. In this case, the union of the specified boxes will be used.

If specified, the residual and/or model images for successful fits will be written.

If an estimates file is not specified, an attempt is made to estimate initial parameters and fit a single Gaussian. If a multiple Gaussian fit is desired, the user must specify initial estimates via a text file (see below for details).

The user has the option of writing the result of the fit to a log file, and has the option of either appending to or overwriting an existing file.

The user has the option of writing the (convolved) parameters of a successful fit to a file which can be fed back to fitcomponents() as the estimates file for a subsequent run.

The user has the option of writing the fit results in tabular format to a file whose name is specified using the summary parameter.

If specified and positive, the value of rms is used to calculate the parameter uncertainties, otherwise, the rms in the selected region in the relevant channel is used for these calculations.

The noisefwhm parameter represents the noise-correlation beam FWHM. If specified as a quantity, it should have angular units. If specified as a numerical value, it is set equal to that number of pixels. If specified and greater than or equal to the pixel size, it is used to calculate parameter uncertainties using the correlated noise equations (see below). If it is specified but less than a pixel width, the the uncorrelated noise equations (see below) are used to compute the parameter uncertainties. If it is not specified and the image has a restoring beam(s), the the correlated noise equations are used to compute parameter uncertainties using the geometric mean of the relevant beam major and minor axes as the noise-correlation beam FWHM. If noisefwhm is not specified and the image does not have a restoring beam, then the uncorrelated noise equations are used to compute the parameter uncertainties.

SUPPORTED UNITS

Currently only images with brightness units conformant with Jy/beam, Jy.km/s/beam, and K are fully supported for fitting. If your image has some other base brightness unit, that unit will be assumed to be equivalent to Jy/pixel and results will be calculated accordingly. In particular, the flux density (reported as Integrated Flux in the logger and associated with the “flux” key in the returned component subdictionary(ies)) for such a case represents the sum of pixel values.

Note also that converting the returned results subdictionary to a component list via cl.fromrecord() currently only works properly if the flux density units in the results dictionary are conformant with Jy. If you need to be able to run cl.fromrecord() on the resulting dictionary you can first modify the flux density units by hand to be (some prefix)Jy and then run cl.fromrecord() on that dictionary, bearing in mind your unit conversion.

If the input image has units of K, the flux density of components will be reported in units of [prefix]Kradrad, where prefix is an SI prefix used so that the numerical value is between 1 and 1000. To convert to units of Kbeam, determine the area of the appropriate beam, which is given by pi/(4ln(2))*bmajbmin, where bmaj and bmin are the major and minor axes of the beam, and convert to steradians (=radrad). This value is included in the beam portion of the component subdictionary (key ’beamster’). Then divide the numerical value of the logged flux density by the beam area in steradians. So, for example

# run on an image with K brightness units
res = imfit(...)
# get the I flux density in K\*beam of component 0
comp = res['results']['component0']
flux\_density_kbeam = comp['flux']['value'][0]/comp['beam']['beamster']

FITTING OVER MULTIPLE CHANNELS

For fitting over multiple channels, the result of the previous successful fit is used as the estimate for the next channel. The number of gaussians fit cannot be varied on a channel by channel basis. Thus the variation of source structure should be reasonably smooth in frequency to produce reliable fit results.

MASK SPECIFICATION

Mask specification can be done using an LEL expression. For example

mask = ’”myimage”’ will use only pixels with values greater than 5.

INCLUDING AND EXCLUDING PIXELS

Pixels can be included or excluded from the fit based on their values using these parameters. Note that specifying both is not permitted and will cause an error. If specified, both take an array of two numeric values.

ESTIMATES

Initial estimates of fit parameters may be specified via an estimates text file. Each line of this file should contain a set of parameters for a single gaussian. Optionally, some of these parameters can be fixed during the fit. The format of each line is

peak intensity, peak x-pixel value, peak y-pixel value, major axis, minor axis, position angle, fixed

The fixed parameter is optional. The peak intensity is assumed to be in the same units as the image pixel values (eg Jy/beam). The peak coordinates are specified in pixel coordinates. The major and minor axes and the position angle are the convolved parameters if the image has been convolved with a clean beam and are specified as quantities. The fixed parameter is optional and is a string. It may contain any combination of the following characters ’f’ (peak intensity), ’x’ (peak x position), ’y’ (peak y position), ’a’ (major axis), ’b’ (axial ratio, R = (major axis FWHM)/(minor axis FWHM)), ’p’ (position angle). NOTE: One cannot hold the minor axis fixed without holding the major axis fixed. If the major axis is not fixed, specifying “b” in the fixed string will hold the axial ratio fixed during the fit.

In addition, lines in the file starting with a # are considered comments.

An example of such a file is:

# peak intensity must be in map units
120, 150, 110, 23.5arcsec, 18.9arcsec, 120deg  
90, 60, 200, 46arcsec, 23arcsec, 140deg, fxp

This is a file which specifies that two gaussians are to be simultaneously fit, and for the second gaussian the specified peak intensity, x position, and position angle are to be held fixed during the fit.

ERROR ESTIMATES

Error estimates are based on the work of Condon 1997, PASP, 109, 166. Key assumptions made are:

  • The given model (elliptical Gaussian, or elliptical Gaussian plus constant offset) is an adequate representation of the data - An accurate estimate of the pixel noise is provided or can be derived (see above). For the case of correlated noise (e.g., a CLEAN map), the fit region should contain many “beams” or an independent value of rms should be provided. - The signal-to-noise ratio (SNR) or the Gaussian component is large. This is necessary because a Taylor series is used to linearize the problem. Condon (1997) states that the fractional bias in the fitted amplitude due to this assumption is of order 1/(SS), where S is the overall SNR of the Gaussian with respect to the given data set (defined more precisely below). For a 5 sigma “detection” of the Gaussian, this is a 4 - All (or practically all) of the flux in the component being fit falls within the selected region. If a constant offset term is simultaneously fit and not fixed, the region of interest should be even larger. The derivations of the expressions summarized in this note assume an effectively infinite region.

Two sets of equations are used to calculate the parameter uncertainties, based on if the noise is correlated or uncorrelated. The rules governing which set of equations are used have been described above in the description of the noisefwhm parameter.

In the case of uncorrelated noise, the equations used are

f(A) = f(I) = f(M) = f(m) = ks(x)/M = ks(y)/m = (s(p)/sqrt(2))*((MM - mm)/(Mm)) = sqrt(2)/S

where s(z) is the uncertainty associated with parameter z, f(z) = s(z)/abs(z) is the fractional uncertainty associated with parameter z, A is the peak intensity, I is the flux density, M and m are the FWHM major and minor axes, p is the position angle of the component, and k = sqrt(8ln(2)). s(x) and s(y) are the direction uncertainties of the component measured along the major and minor axes; the resulting uncertainties measured along the principle axes of the image direction coordinate are calculated by propagation of errors using the 2D rotation matrix which enacts the rotation through the position angle plus 90 degrees. S is the overall signal to noise ratio of the component, which, for the uncorrelated noise case is given by

S = (A/(khr))*sqrt(piMm)

where h is the pixel width of the direction coordinate and r is the rms noise (see the discussion above for the rules governing how the value of r is determined).

For the correlated noise case, the same equations are used to determine the uncertainties as in the uncorrelated noise case, except for the uncertainty in I (see below). However, S is given by

S = (A/(2rN)) sqrt(Mm) (1 + ((NN/(MM)))**(a/2)) (1 + ((NN/(mm)))**(b/2))

where N is the noise-correlation beam FWHM (see discussion of the noisefwhm parameter for rules governing how this value is determined). “” indicates exponentiation and a and b depend on which uncertainty is being calculated. For sigma(A), a = b = 3/2. For M and x, a = 5/2 and b = 1/2. For m, y, and p, a = 1/2 and b = 5/2. f(I) is calculated in the correlated noise case according to

f(I) = sqrt( f(A)*f(A) + (NN/(Mm))*(f(Mf(M) + f(m)*f(m))) )

Note well the following caveats: - Fixing Gaussian component parameters will tend to cause the parameter uncertainties reported for free parameters to be overestimated. - Fitting a zero level offset that is not fixed will tend to cause the reported parameter uncertainties to be slightly underestimated. - The parameter uncertainties will be inaccurate at low SNR (a  10 - If the fitted region is not considerably larger than the largest component that is fit, parameter uncertainties may be mis-estimated. - An accurate rms noise measurement, r, for the region in question must be supplied. Alternatively, a sufficiently large signal-free region must be present in the selected region (at least about 25 noise beams in area) to auto-derive such an estimate. - If the image noise is not statistically independent from pixel to pixel, a reasonably accurate noise correlation scale, N, must be provided. If the noise correlation function is not approximately Gaussian, the correlation length can be estimated using

N = sqrt(2ln(2)/pi)* double-integral(dx dy C(x,y))/sqrt(double-integral(dx dy C(x, y) C(x,y)))

where C(x,y) is the associated noise-smoothing function - If fitted model components have significan spatial overlap, the parameter uncertainties are likely to be mis-estimated (i.e., correlations between the parameters of separate components are not accounted for). - If the image being analyzed is an interferometric image with poor uv sampling, the parameter uncertainties may be significantly underestimated.

The deconvolved size and position angle errors are computed by taking the maximum of the absolute values of the differences of the best fit deconvolved value of the given parameter and the deconvolved size of the eight possible combinations of (FWHM major axis +/- major axis error), (FWHM minor axis +/- minor axis error), and (position andle +/- position angle error). If the source cannot be deconvolved from the beam (if the best fit convolved source size cannot be deconvolved from the beam), upper limits on the deconvolved source size are sometimes reported. These limits simply come from the maximum major and minor axes of the deconvolved gaussians taken from trying all eight of the aforementioned combinations. In the case none of these combinations produces a deconvolved size, no upper limit is reported.

EXAMPLE:

Here is how one might fit two gaussians to multiple channels of a cube using the fit from the previous channel as the initial estimate for the next. It also illustrates how one can specify a region in the associated continuum image as the region to use as the fit for the channel.

imagename = "co_cube.im"
# specify region using region from continuum
region = "continuum.im:source.rgn"
chans = "2~20"
# only use pixels with positive values in the fit
excludepix = [-1e10,0]
# estimates file contains initial parameters for two Gaussians in channel 2
estimates = "initial_estimates.txt"
logfile = "co_fit.log"
# append results to the log file for all the channels
append = "True"
ia.open(imagename)
ia.fitcomponents(region=region, chans=chans, excludepix=excludepix, estimates=estimates, logfile=logfile, append=append)

Parameters

  • box (string='') - Rectangular region(s) to select in direction plane. Default is to use the entire direction plane.

  • region (variant='') - Region selection. Default is to use the full image.

  • chans (variant='') - Channels to use. Default is 0 (first plane).

  • stokes (string='') - The stokes planes to use. Default is to use the first stokes plane.

  • mask (variant='') - Mask to use. Default is none.

  • includepix (doubleVec=[-1]) - Range of pixel values to include. Default is to include all pixels.

  • excludepix (doubleVec=[-1]) - Range of pixel values to exclude. Default is to exclude no pixels.

  • residual (string='') - Name of the residual image to write. Default is not to write the residual.

  • model (string='') - Name of the model image to write. Default is not to write the model.

  • estimates (string='') - Name of the input estimates file. Default is to auto-estimate in which case a single gaussian will be fit.

  • logfile (string='') - File in which to log results. Default is not to write a logfile.

  • append (bool=True) - Append results to logfile? Logfile must be specified. Default is to append. False means overwrite existing file if it exists.

  • newestimates (string='') - File to which to write results in “estimates” format suitable as estimates input for another run. Default is do not write an estimates file.

  • complist (string='') - Output component list table name. Default is do not write a component list table.

  • overwrite (bool=False) - Overwrite component list if it already exists. Default is False.

  • dooff (bool=False) - Also fit a zero level offset? Default is False

  • offset (double=0.0) - Initial estimate of zero-level offset. Only used if doff is True. Default is 0.0

  • fixoffset (bool=False) - Keep the zero level offset fixed during fit? Default is False

  • stretch (bool=False) - Stretch the mask if necessary and possible?

  • rms ({int, double, record, string}='') - 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 ({int, double, record, string}='') - Noise correlation beam FWHM. If numeric value, interpreted as pixel widths. If quantity (dictionary, string), it must have angular units.

  • summary (string='') - File name to which to write table of fit parameters.

Returns

record

fitprofile(box='', region='', chans='', stokes='', axis=-1, mask='', ngauss=1, poly=-1, estimates='', minpts=1, multifit=False, model='', residual='', amp='', amperr='', center='', centererr='', fwhm='', fwhmerr='', integral='', integralerr='', stretch=False, logresults=True, pampest='', pcenterest='', pfwhmest='', pfix='', gmncomps=0, gmampcon='', gmcentercon='', gmfwhmcon='', gmampest=[0.0], gmcenterest=[0.0], gmfwhmest=[0.0], gmfix='', spxtype='', spxest='', spxfix='', div='0', spxsol='', spxerr='', logfile='', append=True, pfunc='', goodamprange=[0.0], goodcenterrange=[0.0], goodfwhmrange=[0.0], sigma='', outsigma='', planes='')[source]

This application simultaneously fits any number of gaussian singlets, any number of lorentzian singlets, and any number of gaussian multiplets, and/or a polynomial to one dimensional profiles using the non-linear, least squares Levenberg-Marquardt algorithm. A description of the fitting algorithm may be found in AIPS++ Note 224 (http://www.astron.nl/casacore/trunk/casacore/doc/notes/224.html) and in Numerical Recipes by W.H. Press et al., Cambridge University Press. A gaussian/lorentzian singlet is a gaussian/lorentzian whose parameters (amplitude, center position, and width) are all independent from any other feature that may be simultaneously fit. A gaussian multiplet is a set of two or more gaussian lines in which at least one (and possibly two or three) parameter of each line is dependent on the parameter of another, single (reference) profile in the multiplet. For example, one can specify a doublet in which the amplitude of the first line is 0.6 times the amplitude of the zeroth line and/or the center of the first line is 20 pixels from the center of the zeroth line, and/or the fwhm of the first line is identical (in pixels) to that of the zeroth line. There is no limit to the number of components one can specify in a multiplet (except of course that the number of parameters to be fit should be significantly less than the number of data points), but there can be only a single reference profile in a multiplet to which to tie constraints of parameters of the other profiles in the set.

Additionally, a power logarithmic polynomial (plp) or a logarithmic tranformed polynomial (ltp) can be fit. In this case, each of these functions cannot be fit simultaneously with any other supported function. 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 = c0*x/div**(c1 + c2*ln(x/div) + c3*ln(x/div)**2 + … + cn*ln(x/div)**(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) = c0 + c1*ln(x/div) + c2*ln(x/div)**2 + … + cn*ln(x/div)**n

The coefficients of the two forms correspond with each other except that c0 in the second equation is equal to ln(c0) of the first. In the case of fitting a spectral index, the spectral index, traditionally represented as alpha, is equal to c1.

In both cases, div is a numerical value used to scale abscissa values so they 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 div, which is determined via

div = 10**int(log10(sqrt(min(x)*max(x))))

where min(x) and max(x) are the minimum and maximum abscissa values, respectively.

So, for example, if S(nu) is proportional to nu**alpha and you expect alpha to be near -0.8 and the value of S(nu) is 1.5 at 1e9 Hz and your image(s) have spectral units of Hz, you would specify spxest=[1.5, -0.8] and div=1e9 when fitting a plp function, or spxest=[0.405, -0.8] and div=1e9 if fitting an ltp function.

More details of fitting all of these functions are described in following sections.

A CAUTIONARY NOTE Note that the likelihood of getting a reliable solution increases with the number of good data points as well as the goodness of the initial estimate. It is possible that the first solution found might not be the best one, and so, if a solution is found, it is recommended that the fit be repeated using the solution of the previous fit as the initial estimatE for the new fit. This process should be repeated until the solutions from one fit to the next differ only insignificantly. The convergent solution is very likely the best solution.

AXIS The axis parameter indicates on which axis profiles should be fit; a value <0 indicates the spectral axis should be used, or if one does not exist, that the zeroth axis should be used.

MINIMUM NUMBER OF PIXELS The minpts parameter indicates the minimum number of unmasked pixels that must be present in order for a fit to be attempted. When multifit=True, positions with too few good points will be masked in any output images.

ONE FIT OF REGION AVERAGE OR PIXEL BY PIXEL FIT The multifit parameter indicates if profiles should be fit at each pixel in the selected region (True), or if the profiles in that region should be averaged and the fit done to that average profile (False).

POLYNOMIAL FITTING The order of the polynomial to fit is specified only via the poly parameter. If poly<0, no polynomial will be fit. No initial estimates of coefficients can be specified; these are determined automatically.

GAUSSIAN SINGLET FITTING In the absence of an estimates file and no estimates being specified by the p*est parameters, and gmncomps=0 or is empty, the ngauss parameter indicates the maximum number of gaussian singlets that should be fit. The initial estimates of the parameters for these gaussians will be attempted automatically in this case. If it deems appropriate, the fitter will fit fewer than this number. In the case where an estimates file is supplied, ngauss is ignored (see below). ngauss is also ignored if the p*est parameters are specified or if gmncomps is not an empty array or, if an integer, is greater than zero. If estimates is not specified or the p*est parameters are not specified and ngauss=0, gmncomps is empty or 0, and poly<0, an error will occur as this indicates there is nothing to fit.

One can specify initial estimates of gaussian singlet parameters via an estimates file or the pampest, pcenterest, pfwhmest, and optionally, the pfix parameters. The latter is the recommended way to specify these estimates as support for estimates files may be deprecated in the future. No matter which option is used, an amplitude initial estimate must always be nonzero. A negative fwhm estimate will be silently changed to positve.

SPECIFYING INITIAL ESTIMATES FOR GAUSSIAN AND LORENTZIAN SINGLETS (RECOMMENDED METHOD) One may specify initial estimates via the pampest, pcenterest, and pfwhmest parameters. In the case of a single gaussian or lorentzian singlet, these parameters can be numbers. pampest must be specified in image brightness units, pcenterest must be given in the number of pixels from the zeroth pixel, and pfwhmest must be given in pixels. Optionally pfix can be specified and in the case of a single gaussian or lorentzian singlet can be a string. In it is coded which parameters should be held constant during the fix. Any combination of “p” (amplitude), “c” (center), or “f” (fwhm) is allowed; eg pfix=”pc” means fix both the amplitude and center during the fit. In the case of more than one gaussian and/or lorentzian singlets, these parameters must be specified as arrays of numbers. The length of the arrays indicates the number of singlets to fit and must be the same for all the p*est parameters.

If no parameters are to be fixed for any of the singlets, pfix can be set to the empty string. However, if at least one parameter of one singlet is to be fixed, pfix must be an array of strings and have a length equal to the p*est arrays. Singlets which are not to have any parameters fixed should be represented as an empty string in the pfix array. So, for example, if one desires to fit three singlets and fix the fwhm of the middle one, one must specify pfix=[“”, “f”, “”], the empty strings indicating no parameters of the zeroth and second singlet should be held constant.

In the case of multifit=True, the initial estimates, whether from the p*est parameters or from a file (see below), will be applied to the location of the first fit. This is normally the bottom left corner of the region selected. If masked, not enough good points to perform a fit, or the attempted fit fails, the fitting proceeds to the next pixel with the pixel value of the lowest numbered axis changing the fastest. Once a successful fit has been performed, subsequent fits will use the results of a fit for a nearest pixel for which a previous fit was successful as the initial estimate for the parameters at the current location. The fixed parameter string will be honored for every fit performed when multifit=True.

One specifies what type of PCF profile to fit via the pfunc parameter. A PCF function is one that can be parameterized by a peak, center, and FWHM, as both gaussian and lorentzian singlets can. If all singlets to be fit are gaussians, one can set pfunc equal to the empty string and all snglets will be assumed to be gaussians. If at least one lorentzian is to be fit, pfunc must be specified as a string (in the case of a single singlet) or an array of strings (in the case of multiple singlets). The position of each string corresponds to the positions of the initial estimates in the p*est and pfix arrays. Minimal match (“g”, “G”, “l”, or “L”) is supported. So, if one wanted to simultaneously fit two gaussian and two lorentzian singlets, the zeroth and last of which were lorentzians, one would specify pfunc=[“L”, “G”, “G”, “L”].

ESTIMATES FILE FOR GAUSSIAN SINGLETS (NONRECOMMENDED METHOD) Initial estimates for gaussian singlets can be specified in an estimates file. Estimates files may be deprecated in the future in favor of the p*est parameters, so it is recommended users use those parameters instead. If an estimates file is desired to be used, the p*est parameters must be 0 or empty and mgncomps must be 0 or empty. Only gaussian singlets can be specified in an estimates file. If one desires to fit one or more gaussian multiplets and/or one or more lorentzian singlets simultaneously, the p*est parameters must be used to specify the initial parameters of all gaussian singlets to fit; one cannot use an estimates file in this case. If an estimates file is specified, a polynomial can be fit simultaneously by specifying the poly parameter. The estimates file must contain initial estimates of parameters for all gaussian singlets to be fit. The number of gaussian singlets to fit is gotten from the number of estimates in the file. The file can contain comments which are indicated by a “#” at the beginning of a line. All non-comment lines will be interpreted as initial estimates. The format of such a line is

[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 image brightness units, while the center must be specified in pixels offset from the zeroth pixel, 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. Any combination of the characters ‘p’ (peak), ‘c’ (center), and ‘f’ (fwhm) are permitted, eg “fc” means hold the fwhm and the center constant during the fit. Fixed parameters will have no error associated with them. Here is an example file:

begin{verbatim} # 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 end{verbatim}

GAUSSIAN MULTIPLET FITTING Any number of gaussian multiplets, each containing any number of two or more components, can be simultaneously fit, optionally with a polynomial and/or any number of gaussian and/or lorentzian singlets, the only caveat being that the number of parameters to be fit should be significantly less than the number of data points. The gmncomps parameter indicates the number of multiplets to fit and the number of components in each multiplet. In the case of a single multiplet, an integer (>1) can be specified. For example, mgncomps=4 means fit a single quadruplet of gaussians. In the case of 2 or more multiplets, and array of integers (all >1) must be specified. For example, gmncomps=[2, 4, 3] means 3 seperate multiples are to be fit, the zeroth being a doublet, the first being a quadruplet, and the second being a triplet.

Initial estimates of all gaussians in all multiplets are specified via the gm*est parameters which must be arrays of numbers. The order starts with the zeroth component of the zeroth multiplet to the last component of the zeroth multiplet, then the zeroth component of the first multiplet to the last compoenent of the first multiplet, etc to the zeroth component of the last multiplet to the last element of the last multiplet. The zeroth element of a multiplet is defined as the reference component of that multiplet and has the special significance that it is the profile to which all constraints of all other profiles in that multiplet are referenced (see below). So, in our example of gmncomps=[2, 4, 3], gmampest, gmcenterest, and gmfwhmest must each be nine (the total number of individual gaussian profiles summed over all multiplets) element arrays. The zeroth, second, and sixth elements represent parameters of the reference profiles in the zeroth, first, and second multiplet, respectively.

The fixed relationships between the non-reference profile(s) and the reference profile of a multiplet are specified via the gmampcon, gmcentercon, and gmfwhmcon parameters. At least one, and any combination, of constraints can be specified for any non-reference component of a multiplet. The amplitude ratio of a non-reference line to that of the reference line is set in gmampcon. The ratio of the fwhm of a non-reference line to that of the reference line is set in gmfwhmcon. The offset in pixels of the center position of a non-reference line to that of the reference line is set in gmcentercon. In the case where a parameter is not constrained for any non-reference line of any multiplet, the value of the associated parameter must be 0. In the case of a single doublet, a constraint may be specified as a number or an array of a single number. For example, mgncomps=2 and gmampcon=0.65 and gmcentercon=[32.4] means there is a single doublet to fit where the amplitude ratio of the first to the zeroth line is constained to be 0.65 and the center of the first line is constrained to be offset by 32.4 pixels from the center of the zeroth line. In cases of a total of three or more gaussians, the constraints parameters must be specified as arrays with lengths equal to the total number of gaussians summed over all multiplets minus the number of reference lines (one per multiplet, or just number of multiplets, since reference lines cannot be constrained by themselves). In the cases where an array must be specified but a component in that array does not have that constraint, 0 should be specified. Here’s an example

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] gfwhmcon=””

In this case we have our previous example of one doublet, one quadruplet, and one triplet. The first component of the doublet has the constraint that its center is offset by 24.2 pixels from the zeroth (reference) component. The first component of the quadruplet is constrained to have an amplitude of 0.2 times that of the quadruplet’s zeroth component and its center is constrained to be offset by 45.6 pixels from the reference component. The second component of the quadruplet is constained to have its center offset by 92.7 pixels from the associated reference component and the third component is constrained to have an amplitude of 0.1 times that of the associated reference component. The first component of the triplet is constrained to have an amplitude of 4.5 times that of its associated reference component and its center is constrained to be offset by -22.8 pixels from the reference component’s center. The second component of the triplet is constrained to have its center offset by -33.5 pixels from the center of the reference component. No lines have FWHM constraints, so the empty string can be given for that parameter. Note that using 0 to indicate no constraint for line center means that one cannot specify a line centered at the same position as the reference component but having a different FWHM from the reference component. If you must specify this very unusual case, try using a very small positive (or even negative) value for the center constraint.

Note that when a parameter for a line is constrained, the corresponding value for that component in the corresponding gm*est array is ignored and the value of the constrained parameter is automatically used instead. So let’s say, for our example above, we had specified the following estimates:

gmampest = [ 1, .2, 2, .1, .1, .5, 3, 2, 5] gmcenterest = [20, 10 , 30, 45.2, 609 , -233, 30, -859, 1]

Before any fitting is done, the constraints would be taken into account and these arrays would be implicitly rewritten as:

gmampest = [ 1, .2, 2, .4, .1, .2, 3, 13.5, 5 ] gmcenterest = [20, 44.2, 30, 75.6, 127.7, -233, 30, 7.2, -3.5]

The value of gmfwhmest would be unchanged since there are no FWHM constraints in this example.

In addition to be constrained by values of the reference component, parameters of individual components can be fixed. Fixed parameters are specified via the gmfix parameter. If no parameters are to be fixed, gmfix can be specified as the empty string or a zero element array. In the case where any parameter is to be fixed, gmfix must be specified as an array of strings with length equal to the total number of components summed over all multiplets. These strings encode which parameters to be fixed for the corresponding components. If a component is to have no parameters fixed, an empty string is used. In other cases one or more of any combination of parameters can be fixed using “p”, “c”, and/or “f” described above for fixing singlet parameters. There are a couople of special cases to be aware of. In the case where a non-reference component parameter is constrained and the corresponding reference component parameter is set as fixed, that parameter in the non-reference parameter will automatically be fixed even if it was specified not to be fixed in the gmfix array. This is the only way the constraint can be honored afterall. In the converse case of when a constrained parameter of a non-reference component is specified as fixed, but the corresponding parameter in the reference component is not specified to be fixed, an error will occur. Fixing an unconstrained parameter in a non-reference component is always legal as is fixing any combination of parameters in a reference component (with the above caveat that corresponding constrained parameters in non-reference components will be silently held fixed as well).

The same rules that apply to singlets when multifit=True apply to multiplets.

LIMITING RANGES FOR SOLUTION PARAMETERS In cases of low (or no) signal to noise spectra, it is still possible for the fit to converge, but often to a nonsensical solution. The astronomer can use her knowledge of the source to filter out obviously bogus solutions. Any solution which contains a NaN value as a value or error in any one of its parameters is automatically marked as invalid.

One can also limit the ranges of solution parameters to known “good” values via the goodamprange, goodcenterrange, and goodfwhmrange parameters. Any combination can be specified and the limit constraints will be ANDed together. The ranges apply to all PCF components that might be fit; choosing ranges on a component by component basis is not supported. If specified, an array of exactly two numerical values must be given to indicate the range of acceptable solution values for that parameter. goodamprange is expressed in terms of image brightness units. goodcenterrange is expressed in terms of pixels from the zeroth pixel in the specified region. goodfwhmrange is expressed in terms of pixels (only non-negative values should be given for FWHM range endpoints). In the case of a multiple-PCF fit, if any of the corresponding solutions are outside the specified ranges, the entire solution is considered to be invalid.

In addition, solutions for which the absolute value of the ratio of the amplitude error to the amplitude exceeds 100 or the ratio of the FWHM error to the FWHM exceeds 100 are automatically marked as invalid.

POWER LOGARITHMIC POLYNOMIAL AND LOGARITHMIC TRANSFORMED POLYNOMIAL FITTING Fitting of a sngle power logarithmic polynomial or a single logarithmic transformed polynomial function is supported. No other functions may be fit simultaneously with either of these; if parameters relating to other functions are supplied simultaneously with parameters relating to these functions, an exception will occur. For details of the functional forms, see the introduction of this document.

The set of c0 … cn coefficients (as defined previously) can be solved for. Initial estimates for the c values should be supplied via the plpest or ltpest parameters, depending on which form is being fit. The number of values given in this array will be the number of coeffecients that are solved for. One may specify which coefficients should be held fixed during the fit in the plpfix or ltpfix array. If supplied, this array should have the same number of elements as its respective initial estimates array. A value of True means the corresponding coefficient will be held fixed during the fit. An empty array indicates that no parameters will be held fixed. This is the default.

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.

INCLUDING STANDARD DEVIATIONS OF PIXEL VALUES If the standard deviations of the pixel values in the input image are known and they vary in the image (eg they are higher for pixels near the edge of the band), they can be included in the sigma parameter. This parameter takes either an array or an image name. The array or image must have one of three shapes: 1. the shape of the input image, 2. the same dimensions as the input image with the lengths of all axes being one except for the fit axis which must have length corresponding to its length in the input image, or 3. be one dimensional with lenght equal the the length of the fit axis in the input image. In cases 2 and 3, the array or pixels in sigma will be replicated such that the image that is ultimately used is the same shape as the input image. The values of sigma must be non-negative. It is only the relative values that are important. A value of 0 means that pixel should not be used in the fit. Other than that, if pixel A has a higher standard deviation than pixel B, then pixel A is noisier than pixel B and will receive a lower weight when the fit is done. The weight of a pixel is the usual

weight = 1/(sigma*sigma)

In the case of multifit=False, the sigma values at each pixel along the fit axis in the hyperplane perpendicular to the fit axis which includes that pixel are averaged and the resultant averaged standard deviation spectrum is the one used in the fit. Internally, sigma values are normalized such that the maximum value is 1. This mitigates a known overflow issue.

One can write the normalized standard deviation image used in the fit but specifying its name in outsigma. This image can then be used as sigma for subsequent runs.

RETURNED DICTIONARY STRUCTURE The returned dictionary has a (necessarily) complex structure. First, there are keys “xUnit” and “yUnit” whose values are the abscissa unit and the ordinate unit described by simple strings. Next there are arrays giving a broad overview of the fit quality. These arrays have the shape of the specified region collapsed along the fit axis with the axis corresponding to the fit axis having length of 1:

attempted: a boolean array indicating which fits were attempted (eg if too few unmasked points, a fit will not be attempted). converged: a boolean array indicating which fits converged. False if the fit was not attempted. valid: a boolean array indicating which solutions fall within the specified valid ranges of parameter space (see . section LIMITING RANGES FOR SOLUTION PARAMETERS for details). niter: an int array indicating the number of iterations for each profile, <0 if the fit did not converge ncomps: the number of components (gaussian singlets + lorentzian singlets + gaussian multiplets + polynomial) fit for the profile, . <0 if the fit did not converge direction: a string array containing the world direction coordinate for each profile

There is a “type” array having number of dimensions equal to the number of dimensions in the above arrays plus one. The shape of the first n-1 dimensions is the same as the shape of the above arrays. The length of the last dimension is equal to the number of components fit. The values of this array are strings describing the components that were fit at each possition (“POLYNOMIAL”, “GAUSSIAN” in the case of gaussian singlets, “LORENTZIAN” in the case of lorentzian singlets, and “”GAUSSIAN MULTPLET”).

If any gaussian singlets were fit, there will be a subdictionary accessible via the “gs” key which will have subkeys “amp”, “ampErr”, “center”, “centerErr”, “fwhm”, “fwhmErr, “integral”, and “integralErr”. Each of these arrays will have one more dimension than the overview arrays described above. The shape of the first n-1 dimensions will be the same as the shape of the arrays described above, while the final dimension will have length equal to the maximum number of gaussian singlets that were fit. Along this axis will be the corresponding fit result or associated error (depending on the array’s associated key) of the fit for that singlet component number. In cases where the fit did not converge, or that particular component was excluded from the fit, a value of NAN will be present.

If any lorentzian singlets were fit, their solutions will be accessible via the “ls” key. These arrays follow the same rules as the “gs” arrays described above.

If any gaussian multiplets were fit, there will be subdictionaries accessible by keys “gm0”, “gm1”, …, “gm{n-1}” where n is the number of gaussian muliplets that were fit. Each of these dictionaries will have the same arrays described above for gaussian singlets. The last dimension will have length equal to the number of components in that particular multiplet. Each pixel along the last axis will be the parameter solution value or error for that component number in the multiplet, eg the zeroth pixel along that axis contains the parameter solution or error for the reference component of the multiplet.

The polynomial coefficient solutions and errors are not returned, although they are logged.

If a power logarithmic polynomial was fit, there will be a subdictionary accessible via the “plp” key which will have subkeys “soltuion” and “error” which will each have an array value. Each of these arrays will have one more dimension than the overview arrays described above. The shape of the first n-1 dimensions will be the same as the shape of the overview arrays described above, while the final dimension will have length equal to the number of parameters that were fit. Along this axis will be the corresponding fit result or associated error (depending on the array’s associated key) of the fit. In cases where the fit was not attempted or did not converge, a value of NAN will be present.

OUTPUT IMAGES In addition to the returned dictionary, optionally one or more of any combination of output images can be written. The model and residual parameters indicate the names of the model and residual images to be written; blank values inidcate that these images should not be written.

One can also write none, any or all of the solution and error images for gaussian singlet, lorentzian singlet, and gaussian multiplet fits via the parameters amp, amperr, center, centererr, fwhm, fwhmerr, integral, integralerr when doing multi-pixel fits. For a power logarithmic polynomial or a logarithmic transformed polynomial fit, plpsol or ltpsol and plperr or ltpsol are the names of the solution and error images to write, respectively.

These images contain the arrays described for the associated parameter solutions or errors described in previous sections. Each component is written to a different image, and each image is distiguished by the component it represents by its name ending in an uderscore and the relevant component number (“_0”, “_1”, etc). In the case of Gaussian multiplets, the image name ends with the number of the mulitplet group followed by the number of the component in that group (eg “_3_4” represents component 4 in multiplet group 3). In the case of lorentzian singlets, “_ls” is appended to the image names (but before the identifying component number), in the case of gaussian multiplets. Similarly “_gm” is included in the name of Gaussian multiplet images. Pixels for which fits were not attempted, did not converge, or converged but have values of NaN (not a number) or INF (infinity) will be masked as bad.

Writing analogous images for polynomial coefficients is not supported.

Parameters

  • box (string='') - Rectangular region to select in direction plane. Default is to use the entire direction plane.

  • region ({string, string, record}='') - Region selection. Default is to use the full image.

  • chans (string='') - Channels to use. Channels must be contiguous. Default is to use all channels..

  • stokes (string='') - Stokes planes to use. Planes must be contiguous. Default is to use all stokes planes.

  • axis (int=-1) - The profile axis. Default: use the spectral axis if one exists, axis 0 otherwise (<0).

  • mask ({string, stringVec}='') - Mask to use. Default is none.

  • ngauss (int=1) - Number of Gaussian elements. Default: 1.

  • poly (int=-1) - Order of polynomial element. Default: do not fit a polynomial (<0).

  • estimates (string='') - Name of file containing initial estimates. Default: No initial estimates (“”).

  • minpts (int=1) - Minimum number of unmasked points necessary to attempt fit.

  • multifit (bool=False) - 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.

  • model (string='') - Name of model image. Default: do not write the model image (“”).

  • residual (string='') - Name of residual image. Default: do not write the residual image (“”).

  • amp (string='') - Prefix of name of amplitude solution image. Name of image will have gaussian component number appended. Default: do not write the image (“”).

  • amperr (string='') - Prefix of name of amplitude error solution image. Name of image will have gaussian component number appended. Default: do not write the image (“”).

  • center (string='') - Prefix of name of center solution image. Name of image will have gaussian component number appended. Default: do not write the image (“”).

  • centererr (string='') - Prefix of name of center error solution image. Name of image will have gaussian component number appended. Default: do not write the image (“”).

  • fwhm (string='') - Prefix of name of FWHM solution image. Name of image will have gaussian component number appended. Default: do not write the image (“”).

  • fwhmerr (string='') - Prefix of name of FWHM error solution image. Name of image will have gaussian component number appended. Default: do not write the image (“”).

  • integral (string='') - Prefix of name of integral solution image. Name of image will have gaussian component number appended. Default: do not write the image (“”).

  • integralerr (string='') - Prefix of name of integral error solution image. Name of image will have gaussian component number appended. Default: do not write the image (“”).

  • stretch (bool=False) - Stretch the mask if necessary and possible?

  • logresults (bool=True) - Output results to logger?

  • pampest ({double, doubleVec}='') - Initial estimate PCF profile amplitudes.

  • pcenterest ({double, doubleVec}='') - Initial estimate PCF profile centers, in pixels.

  • pfwhmest ({double, doubleVec}='') - Initial estimate PCF profile FWHMs, in pixels.

  • pfix ({string, stringVec}='') - PCF parameters to fix during fit. Any combination of “p”, “c”, or “f”.

  • gmncomps ({int, intVec}=0) - Number of components in each gaussian multiplet to fit

  • gmampcon ({double, doubleVec}='') - The amplitude ratio constraints for non-reference components to reference component in gaussian multiplets.

  • gmcentercon ({double, doubleVec}='') - The center offset constraints (in pixels) for non-reference components to reference component in gaussian multiplets.

  • gmfwhmcon ({double, doubleVec}='') - The FWHM ratio constraints for non-reference components to reference component in gaussian multiplets.

  • gmampest (doubleVec=[0.0]) - Initial estimate of individual gaussian amplitudes in gaussian multiplets.

  • gmcenterest (doubleVec=[0.0]) - Initial estimate of individual gaussian centers in gaussian multiplets, in pixels.

  • gmfwhmest (doubleVec=[0.0]) - Initial estimate of individual gaussian FWHMss in gaussian multiplets, in pixels.

  • gmfix ({string, stringVec}='') - Parameters of individual gaussians in gaussian multiplets to fix during fit.

  • spxtype (string='') - Type of function to fit. “plp” $=>$ power logarithmic polynomial, “ltp” $=>$ logarithmic transformed polynomial.

  • spxest (doubleVec='') - REQUIRED. Initial estimates as array of numerical values for the spectral index function coefficients. eg [1.5, -0.8] if fitting a plp function thought to be close to 1.5*(x/div)**(-0.8) or [0.4055, -0.8] if fitting an lpt function thought to be close to ln(1.5) - 0.8*ln(x/div).

  • spxfix (boolVec='') - Fix the corresponding spectral index function coefficients during the fit. True$=>$hold fixed.

  • div (variant='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.

  • spxsol (string='') - Name of the spectral index function coefficient solution image to write.

  • spxerr (string='') - Name of the spectral index function coefficient error image to write.

  • logfile (string='') - File in which to log results. Default is not to write a logfile.

  • append (bool=True) - Append results to logfile? Logfile must be specified. Default is to append. False means overwrite existing file if it exists.

  • pfunc ({string, stringVec}='') - PCF singlet functions to fit. “gaussian” or “lorentzian” (minimal match supported). Unspecified means all gaussians.

  • goodamprange (doubleVec=[0.0]) - Acceptable amplitude solution range. [0.0] $=>$ all amplitude solutions are acceptable.

  • goodcenterrange (doubleVec=[0.0]) - Acceptable center solution range in pixels relative to region start. [0.0] $=>$ all center solutions are acceptable.

  • goodfwhmrange (doubleVec=[0.0]) - Acceptable FWHM solution range in pixels. [0.0] $=>$ all FWHM solutions are acceptable.

  • sigma ({string, intVec, intArray, doubleVec, doubleArray}='') - Standard deviation array or image name.

  • outsigma (string='') - Name of output image used for standard deviations. Ignored if sigma is empty.

  • planes (intVec='') - Planes along fit axis to use in the fit. Empty means use all planes. All values must be non-negative.

Returns

record

Examples

ia.open("myspectrum.im")
res = ia.fitprofile(ngauss=2, box="3,3,4,5", poly=2, multifit=True)
fitsheader(retstr=False, exclude='')[source]

This method constructs a FITS header dictionary which describes the current image. The exclude parameter can either be a string or a list of strings. It indicates FITS header keywords to exclude. This may be useful for the HISTORY keyword because the HISTORY can be very long. If the retstr parameter is set to True, then a string will be returned representing the FITS header. The exculde parameter has no effect when the header information is returned as a single string. An new IMAGENME keyword is added to the header. The value for IMAGENME is the string returned by name method, but the name is truncated to the maximum size for a FITS header value, 68 characters. The image name is truncated from the beginning of the string to reduce it to 68 characters, so the value will contain either the full name or the final 68 characters of the name.

Parameters

  • retstr (bool=False) - return a string representing the FITS header instead of a record

  • exclude (stringVec='') - Header fields to exclude from the return record. Fields cannot be excluded when retstr=True.

Returns

variant

fromarray(outfile='', pixels='', csys='', linear=False, overwrite=False, log=True, type='f')[source]

This application converts a numerical numpy array of any size and dimensionality into a CASA image. It will create float, double, complex-float, and complex-double valued images.

The image analysis tool on which this method is called will reference the created image; if this tool referenced another image before this call, that image will no longer be referenced by the tool after the creation of the new image. If you would rather have a new image analysis tool returned, keeping the one on which this method is called unaltered, use newimagefromarray() instead. If outfile is specified, a persistent image is written to disk, if not, the image tool on which this method was called will reference a temporary image (either in memory or on disk, depending on its size) that will be deleted when the tool is closed.

The type parameter controls the data type/precision of the pixel values of the created image. ’f’ indicates that float precision point (32 bit precision) pixel values should be writted. ’d’ indicates that double precision (64 bit precision) pixel values should be written. If the input array has complex (as opposed to real) values, then complex pixel values, with each of the real and imaginary parts having the specified precision, will be written. Array values will be cast automatically to the specified precision, so that the precision of the input array values may be increased, decreased, or unchanged depending on the input array type.

The coordinate system, provided as a a dictionary (use eg, cs.torecord() to do that), is optional. If specified, it must have the same number of dimensions as the pixels array. Call the naxes() method on the coordinate system tool to see how many dimensions the coordinate system has. A coordinate system can be created from scratch using the coordinate system (cs) tool and methods therein, but often users prefer to use a coordinate system from an already existing image. This can be gotten using ia.coordsys() which returns a coordinate system tool. A torecord() call on that tool will result in a python dictionary describing the coordinate system which is the necessary format for the csys input parameter of ia.fromarray().

If csys is not specified, a default coordinate system will be created. If linear=False (the default), the created coordinate system will have standard RA/DEC/Stokes/Spectral Coordinate axes depending upon the shape of the pixels array (Stokes axis must be no longer than 4 pixels and the spectral axis may precede the Stokes axis if eg, shape=[64,64,32,4]. Extra dimensions are given linear coordinates. If linear=True, then all the resulting coordinates are linear with the axes represent lengths. In this case each axis will have a value of 0.0 at its center pixel. The increment of each axis will be 1.0 km.

The method returns True if creation of the image was successful, False otherwise, so a check can be made programmatically if the image creation was successful.

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • pixels (variant='') - Numeric array

  • csys (record='') - Coordinate System. Default is to construct an appropriate coordinate system given the array shape.

  • linear (bool=False) - Make a linear Coordinate System if csys not given

  • overwrite (bool=False) - Overwrite the output image if it already exists?

  • log (bool=True) - Write image creation messages to logger

  • type (string='f') - Pixel data type to write. ‘f’ (float precision) or ‘d’ (double precision)

Returns

bool

Examples

# make an image with a default RA/Dec/Stokes/Frequency coordinate system
# having all pixels set to 2.5.
ary = ia.makearray(v=2.5, shape=[64, 64, 4, 128])
# the ia tool does not need to reference an image in this case (ie open()
# need not have been called), if it does reference another image, that reference
# will be lost and replaced with a reference to the newly created image.
res = ia.fromarray(outfile='test.data', pixels=ary, overwrite=True)
if res:
    # perform operations on the newly created image if desired and make sure
    # to close it when done to free up system resources (eg memory)
    ia.shape()
ia.done()
fromcomplist(outfile='', shape='', cl='', csys='', overwrite=False, log=True, cache=True)[source]

This method allows one to create an image based on a component list. A component list is a list of simple models (point sources, Gaussians, disks, etc) that describe the sky brightness (cf the component list (cl) tool). Images that can be described in this way normally require significantly less space to store than traditional images in which all the pixel values are stored. For a component list image, pixel values are computed “on the fly”. Pixel values can be cached by specifying cache=True (the default value) while the image is attached to an image tool, which permits faster access to them after they are computed the first time. The trade off to caching is that resources such as memory and disk space must be used to cache the pixel values.

The image is constrained to have two, three, or four dimensions. One must specify an image shape (the dimensionality of which must adhere to this constraint). One may also supply a coordinate system specification using the csys parameter. If a coordinate system is not specified, a default coordinate system is used. If specified, the coordinate system must have a direction coordinate which has two pixel axes. It can also have a spectral and/or polarization coordinate. The maximum length of the polarization coordinate is four pixels, and the world coordinate values of the polarization coordinate are constrained to be in the set of stokes parameters I, Q, U, and V.

As is common with image creation methods, specifying an empty string for the outfile parameter results in a tempoary image being created that will be deleted when either the done() or close() method is run on the tool. By specifying a non-empty string, the image is saved to disk and can be opened with the open() method later. A persistent component list image is composed of a component list table that has metadata describing the image-related information, such as the coordinate system and the shape, as well as a history (log) table.

Because pixel values are computed from the component models, altering pixel values is not supported. So methods such as putchunk(), putregion(), and addnoise() will fail on component list images when trying to modify pixel values. However, persistent image masks and on the fly masks are fully supported.

The brightness unit of a component list image is constrained to be “Jy/pixel”. Attempts to modify this value using setbrightnessunit() will fail.

Component list images do not support synthesized beams; attempting to run setrestoringbeam() on a component list image to add a beam(s) will fail.

One can easily create an image in which the pixel values are persistently stored from a component list image by running methods such as fromimage(), subimage(), tofits(), etc. In general, any method run on a component list image that creates a new image will create a non-component list image (eg, a traditional CASA Paged image or Temporary image) in which the pixel values are explicitly stored.

DISK MODELS

Pixels with centers inside the disk will have the same values, even if a pixel straddles the edge of the disk. Pixels with straddle the edge of the disk which have centers outside the disk are given values of zero. Thus, one should not expect the flux density of the disk to be exactly the provided value to the component list; for a given size disk, the computed flux density will be closer to the expected value for images with smaller pixels.

Parameters

  • outfile (string='') - Name of output image file. Default is unset.

  • shape (intVec='') - Shape of image. Must be specified and have either two, three, or four elements.

  • cl ({string, record}='') - Component list as a dictionary or name of component list table. Must be specified

  • csys (record='') - Coordinate System. Default uses a four dimensional “standard” RA/Dec/Stokes/Freq coordinate system.

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • log (bool=True) - Write image creation messages to logger

  • cache (bool=True) - Cache pixel values for use while image is in use?

Returns

bool

Examples

# use an existing component list and the coordinate system from an existing image
# to create a component list image. To prevent the size of the files associated with
# the component list from increasing, open the list read only
cl.open("mymodel.cl", nomodify=True)
clrec = cl.torecord()
cl.done()
# say this image has axes RA x Dec x Stokes x Freq
ia.open("myimage.im")
csys = ia.coordsys().torecord()
ia.done()
ia.fromcomplist(outfile="mycomplist.im", shape=[512, 512, 4, 20], cl=clrec,
    csys=csys, cache=True)
# get statistics
stats = ia.statistics()
# do other operations, like convolving the image
im2 = ia.convolve2d(outfile="convolved.im", axes=[0,1], type="gauss", "
    major='20arcsec', minor='10arces', pa='45deg')
# be sure to close the image tools when done with them
ia.done()
im2.done()
fromfits(outfile='', infile='', whichrep=0, whichhdu=0, zeroblanks=False, overwrite=False)[source]

This function is used to convert a FITS disk image file (Float, Double, Short, Long are supported) to an  . If outfile is given, the image is written to the specified disk file. If outfile is unset, the Image  is associated with a temporary image. This temporary image may be in memory or on disk, depending on its size. When you close the Image  (with the close function) this temporary image is deleted.

This function reads from the FITS primary array (when the image is at the beginning of the FITS file; whichhdu=0), or an image extension (when the image is elsewhere in the FITS file, whichhdu \(\>\) 0).

By default, any blanked pixels will be converted to a mask value which is False, and a pixel value that is NaN. If you set zeroblanks=T then the pixel value will be zero rather than NaN. The mask will still be set to False. See the function replacemaskedpixels if you need to replace masked pixel values after you have created the image.

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • infile (string='') - Input FITS disk file name. Must be specified.

  • whichrep (int=0) - If this FITS file contains multiple coordinate representations, which one should we read (0-based)

  • whichhdu (int=0) - If this FITS file contains multiple images, which one should we read (0-based).

  • zeroblanks (bool=False) - If there are blanked pixels, set them to zero instead of NaN

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

Returns

bool

Examples

#
print "\t----\t fromfits Ex 1 \t----"
datapath=pathname+'/data/demo/Images/imagetestimage.fits'
ia.fromfits('./myimage', datapath, overwrite=True)
print ia.summary()
s = ia.miscinfo()
print s.keys()  #  prints any unrecognized field names
ia.close()
#



The FITS image is converted to a \casa\ \imagefile\ and access is
provided via the default \imagetool\ called {\stf ia}.  Any FITS
header keywords which were not recognized or used are put in the
miscellaneous information bucket accessible with the miscinfo function.  In
the example we list the names of the fields in this record.
fromimage(outfile='', infile='', region='', mask='', dropdeg=False, overwrite=False)[source]

This function applies a  to an , creates a new  containing the (sub)image, and associates the  with it.

The input image file may be in native , , or Miriad format. Look for more information on foreign images.

If outfile is given, the (sub)image is written to the specified disk file.

If outfile is unset, the Image  actually references the input image file. So if you deleted the input image disk file, it would render this  useless. When you close this  (with the close function) the reference connection is broken.

Sometimes it is useful to drop axes of length one (degenerate axes). Use the dropdeg argument if you want to do this.

The output mask is the combination (logical OR) of the default input  (if any) and the OTF mask. Any other input  will not be copied. Use function maskhandler if you need to copy other masks too.

See also the subimage function.

Parameters

  • outfile (string='') - Output sub-image file name. Default is unset.

  • infile (string='') - Input image file name. Must be specified.

  • region (variant='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • dropdeg (bool=False) - Drop degenerate axes

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

Returns

bool

Examples

#
print "\t----\t fromimage Ex 1 \t----"
innerquarter = rg.box([0.25,0.25],[0.75,0.75],frac=True)
ia.close()
ia.fromimage(outfile='image.small', infile='test.data', region=innerquarter, overwrite=True)
ia.close()
#


The specified \region\ takes a quarter by area of the first two axes
of the image, and all pixels of other axes.
fromrecord(record='', outfile='')[source]

You can convert an associated image to a record (torecord) or imagepol tool functions will sometimes give you a record. This function (fromrecord) allows you to set the contents of an image tool to the content of the record. This and torecord are used for deserialization and serialization.

Parameters

  • record (record='') - Record containing the image

  • outfile (string='') - The name of the diskfile to be created for image from record

Returns

bool

Examples

#
print "\t----\t fromrecord Ex 1 \t----"
ia.maketestimage('image.large', overwrite=True)
rec=ia.torecord()
ia.close()
ia.fromrecord(rec, "testimage")
fromshape(outfile='', shape=[0], csys='', linear=False, overwrite=False, log=True, type='f')[source]

This function creates a CASA image with the specified shape. All the pixel values in the image are set to 0. One may create an image with float precision pixels (type=’f’), complex float precision pixels (type=’c’), double precision pixels (type=’d’), or complex double precision pixels (’cd’). To use a numpy array of values to create an image, use ia.fromarray(). To make a 2-D image from a packaged FITS file, use ia.maketestimage().

If outfile is given, the image is written to the specified disk file. If outfile is unset, the image analysis tool is associated with a temporary image. This temporary image may be in memory or on disk, depending on its size. When you close the image analysis tool (with the ia.close() method, the temporary image is deleted.

The coordinate system, provided as a coordinate system tool record, is optional. If provided, it must be dimensionally consistent with the specified shape.

If the coordinate system is not provided, a default coordinate system will be created. If linear=False (the default), then it is a standard RA/DEC/Stokes/Spectral coordinate system depending exactly upon the shape (the Stokes axis must be no longer than 4 pixels and spectral axis may occur prior to the Stokes axis if eg, shape=[64,64,32,4]. Extra dimensions are given linear coordinates. If linear=True, then the coordinate system will have linear coordinates.

The method returns True if successful, False otherwise.

Parameters

  • outfile (string='') - Name of output image file. Default is unset.

  • shape (intVec=[0]) - Shape of image. Must be specified.

  • csys (record='') - Coordinate System. Default is unset.

  • linear (bool=False) - Make a linear Coordinate System if csys not given?

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • log (bool=True) - Write image creation messages to logger

  • type (string='f') - Data type of image pixels. ‘f’ (float), ‘c’ (complex float), ‘d’ (double), ‘cd’ (complex double).

Returns

bool

Examples

# Create a 64 x 64 x 128 RA-Dec-Frequency image
ia.fromshape('test2.data', [64,64,128], overwrite=True)
mycs = ia.coordsys(axes=[0,1])
ia.close()
# create temporary image using first and third axes of
# previously created image
ia.fromshape(shape=[10, 20], csys=mycs.torecord())
mycs.done()
ia.close()


The first example creates a zero-filled \imagefile\ named {\sff
test2.data} of shape [64,64,128].  If you examine the header with {\stff
ia.summary()} you will see the RA/DEC/Spectral coordinate information. 
In the second example, a Coordinate System describing the first and second axes
of the image {\sff test.data} is created and used to create a 2D temporary image.

# create a 4-D image populated with noise with some interesting metadata
ia.fromshape("noisy.im", [100, 100, 4, 20])
# add the noise
ia.addnoise()
# add a beam
ia.setrestoringbeam(major="4arcmin", minor="3arcmin", pa="20deg")
# add brightness unit
ia.setbrightnessunit("Jy/beam")
ia.done()
# set an interesting object name
imd.open("noisy.im")
imd.set("object","thor")
imd.done()
getchunk(blc=[-1], trc=[-1], inc=[1], axes=[-1], list=False, dropdeg=False, getmask=False)[source]

This function returns the pixels (or optionally the pixel mask) from the attached image between blc and trc, inclusive. Images with float, complex float, double, and complex double precision pixel values are supported. An increment may be specified with inc. Note that if too many pixel values are retrieved, swapping may occur, result in a decrease in performance, since the pixel values are stored in memory.

Any illegal blc values are set to zero. Any illegal trc values are set to the end of the image. If any trc values are less than corresponding blc values, all the pixel values for that axis are returned. Any illegal inc values are set to unity.

The axes parameter can be used to reduce the dimensionality of the output array. It specifies which pixel axes of the image over which to average the data. For example, consider a 3-D image, with axes=[0,1] and all other parameters set to their defaults. The result would be a 1-D vector, a profile along the third axis, with the data averaged over the first two axes.

A related function is getregion(), which retrieves the pixel values or pixel mask from a potentially more complex region. Method getchunk() is retained because it is faster and therefore preferable for repeated operations in loops if the pixel mask is not required and the region is a simple box.

If getmask=True, the return value is the pixel mask values, rather than the pixel values.

Parameters

  • blc (intVec=[-1]) - Bottom-Left-Corner (beginning) of pixel section. Default is start of image.

  • trc (intVec=[-1]) - Top-Right-Corner (end) of pixel section. Default is end of image.

  • inc (intVec=[1]) - increment (stride) along axes

  • axes (intVec=[-1]) - Axes to average over. Default is none.

  • list (bool=False) - List bounding box to logger?

  • dropdeg (bool=False) - Drop degenerate axes?

  • getmask (bool=False) - Get the pixel mask rather than the pixel values

Returns

any

Examples

Suppose that we have a 3-dimensional image called {\sff im}. Then:


#
print "\t----\t getchunk Ex 1 \t----"
ia.fromshape(shape=[64,64,128])
pix = ia.getchunk()                      # all pixels
ia.calcmask('T')                         # give image a mask
pix = ia.getchunk([1,1,1], [10,10,1])    # 10 by 10 section of plane # 1
pix = ia.getchunk([1,1], [1,1])          # first spectrum
pix = ia.getchunk(inc=[1,5])             # all planes, decimated by 5 in y
mask = ia.getchunk(getmask=True)            # Get pixelmask
ia.close()
#
getprofile(axis=-1, function='mean', region='', mask='', unit='', stretch=False, spectype='default', restfreq='', frame='', logfile='')[source]

This application returns information on a one-dimensional profile taken along a specified image axis. The region of interest is collapsed (a’la ia.collapse() along all axes orthogonal to the one specified, and) the specified aggregate function is applied to these pixels to generate the returned values.

The aggregate function must be one of the functions supported by ia.collapse; ie, ’flux’, ’madm’, ’max’, ’mean’, ’median’, ’min’, ’rms’, ’stdev’, ’sum’, ’variance’, and ’xmadm’. See the help for ia.collapse() for details regarding these functions. Minimum match and case insenstivity is supported. In addition, single binary (addition, subtraction, multiplication, and division) operations of these functions are supported, eg function=’maxmin’ will return data that is the product of the maximum and the mininum for each plane along the specified axis.

One may specify the unit of the returned coordinate values. Unless axis is the spectral axis, unit must be conformant with the corresponding axis unit in the image coordinate system or it must be ’pixel’ which signifies, pixel, rather than world, coordinate values should be calculated. If axis is the spectral axis, unit may be a velocity unit (assuming the coordinate system has a rest frequency or restfreq is specified) or a length unit. In these cases, the returned coordinate values will be converted to velocity or wavelength, respectively.

The parameter spectype may be used to specify the velocity or wavelength type for the returned coordinate values if profile is taken along spectral axis. Supported (minimum match, case insensitive) values) are “relativistic velocity”, “beta”, “radio velocity”, “optical velocity”, “wavelength”, “air wavelength”, “default”. The “default” value is equivalent to “relativistic” if unit is a velocity unit or “wavelength” if unit is a length unit.

The restfreq parameter allows one to set the rest frequency for the coordinates to be returned if axis is the spectral axis and unit is a velocity unit. If blank, the rest frequency associated with the image coordinate system is used.

The frame allows one to specify which kinematic reference frame that the returned coordinate values should be calculated in. It is only used if axis is the spectral axis and unit is unspecified or is specified and a frequency unit. If blank, the reference frame associated with the image coordinate system is used.

The returned dictionary contains the keys:

values: one-dimensional array along the specified axis containing values resulting from applying the specified aggregate function to corresponding pixels at the same location along that axis. mask: one-dimensional array of booleans of the resulting mask after applying the aggregate function, formed in the same way as that formed by ia.collapse. coords One-dimensional array of corresponding coordinate values along the specified axis in the specified unit (or the unit associated with the axis in the image coordinate system if unspecified). xUnit The unit used for calculating the values the coords array.

Parameters

  • axis (int=-1) - Axis along which to determine profile. Must be specified

  • function (string='mean') - Aggregate function to apply for collapse along axes orthogonal to specified axis.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (string='') - Mask to use. Default is none.

  • unit (string='') - Unit of the returned abscissa values. Must either be ‘pixel’ or be conformant with image axis unit unless axis is the spectral axis. Default is the unit associated with axis in the image coordinate system.

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

  • spectype (string='default') - Velocity or wavelength type if profile taken along spectral axis. Supported (minimum match, case insensitive) values are “relativistic velocity”, “beta”, “radio velocity”, “optical velocity”, “wavelength”, “air wavelength”, “default”.

  • restfreq (variant='') - Rest frequency to use when calculating coordinate values. Used only if axis is spectral axis and unit is not the unit associated with the axis in the coordinate system. Empty string means use the rest frequency associated with the image coordinate system

  • frame (string='') - Reference frame to use when calculating coordinate values. Used only if axis is spectral axis and unit is not the unit associated with the axis in the coordinate system. Empty string means use the reference frame associated with the image coordinate system

  • logfile (string='') - File to which to write profile.

Returns

record

Examples

ia.open('myimage')
# get the max pixel values along axis 2
res = ia.getprofile(axis=2, function='max')

# axis 2 is the spectral axis. Get the minimum pixel values along this axis
# and specify that the returned coordinate values should be optical velocities
# in km/s
res2 = ia.getprofile(axis=2, function='min', unit='km/s', spectype='optical')

# get the result of the maximum and the median of the absolute deviation from the median
# along this axis
res3 = ia.getprofile(axis=2, function='max/madm')

ia.done()
getregion(region='', axes=[-1], mask='', list=False, dropdeg=False, getmask=False, stretch=False)[source]

This function recovers the image pixel or pixel mask values in the given region of interest. Regardless of the shape of the specified, the shape of the pixels and pixelmask arrays must necessarily be the bounding box of the specified region. If the region extends beyond the image, it is truncated.

Recall that the recovered pixelmask will reflect both the pixelmask stored in the image, and the region (their masks are ’anded’ together).

The argument axes can be used to reduce the dimensionality of the output array. It specifies which pixel axes of the image to average the data over. For example, consider a 3-D image. With axes=[0,1] and all other arguments left at their defaults, the result would be a 1-D vector, a profile along the third axis, with the data averaged over the first two axes.

This method differs in a couple of ways from the getchunk() method. For example, the specified region can be much more complex (eg, a union of polygons) than the limited, simple regions that can be specified in getchunk(), which must be rectangular. On the other hand, getregion() is less effective than getchunk() for the same region specification. So if one is interested in say, iterating through an image, getting a regular hyper-cube of pixels and doing something with them, then getchunk() will be faster. This would be especially noticeable if you iterated line by line through a large image (and of course, in both cases, retrieving very large regions will become very resource intensive, as these returned arrays are completely stored in memory).

Parameters

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • axes (intVec=[-1]) - Axes to average over. Default is none.

  • mask (variant='') - Mask to use. Default is none.

  • list (bool=False) - List the bounding box to the logger

  • dropdeg (bool=False) - Drop degenerate axes

  • getmask (bool=False) - Get the pixel mask rather than pixel values

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

Returns

variant

Examples

Suppose that we have a 3-dimensional image called {\sff cube} and wish
to recover the pixel from a simple regular region.


#
print "\t----\t getregion Ex 1 \t----"
ia.fromshape('cube', [64,64,64], overwrite=True)
#r1=rg.box(blc=[10,10,10],trc=[30,40]) # Create region
r1=rg.box([10,10,10],[30,40,40]) # Create region
pixels=ia.getregion(r1)
ia.close()
#
getslice(x='', y='', axes=[0, 1], coord=[-1], npts=0, method='linear')[source]

This function returns a 1-D slice (the pixels and opionally the pixel mask) from the attached image. The slice is constrained to lie in a plane of two cardinal axes (e.g. XY or YZ). Interpolation is permitted between pixels, and a set of interpolation schemes is available.

The slice is specified as a polyline giving the x and y coordinates and the axes of the plane holding that slice. The absolute pixel coordinates of the other axes may be specified using the coord parameter. If not specified, these values default to pixel 0 on the relevant axes.

The npts parameter allows the number of values to be returned.

The method parameter allows specification of the interpolation method to be used. Allowed values are ’nearest’, ’linear’, and ’cubic’. In the case of an image with complex valued pixels, the interpolation is done independently on the real and imaginary values. For example, the linearly interpolated value midway between pixels with values of 1 + 5j and 2 + 7j would be 1.5 + 6j.

The return value is a dictionary with keys ’pixels’ (interpolated pixel values), ’mask’ (interpolated mask), ’xpos’ (x-location in absolute pixel coordinates), ’ypos’ (y-location in absolute pixel coordinates), ’distance’ (distance along slice in pixels), and ’axes’ (the x and y axes of slice).

Parameters

  • x (doubleVec='') - Polyline x vertices in absolute pixel coordinates

  • y (doubleVec='') - Polyline y vertices in absolute pixel coordinates

  • axes (intVec=[0, 1]) - Pixel axes of plane holding slice. Default is first two axes.

  • coord (intVec=[-1]) - Specify pixel coordinate for other axes. Default is first pixel.

  • npts (int=0) - Number of points in slice. Default is auto determination.

  • method (string='linear') - The interpolation method, String from ‘nearest’, ‘linear’, ‘cubic’

Returns

record

Examples

Suppose that we have a 2-dimensional image. Then:


#
print "\t----\t getslice Ex 1 \t----"
ia.maketestimage();
rec = ia.getslice (x=[1,20], y=[2,30])     # SLice from [1,2] -\> [20,30]
print rec.keys()
#['distance', 'xpos', 'axes', 'mask', 'ypos', 'pixel']
rec = ia.getslice (x=[1,20,25,11], y=[2,30,32,40]) # Polyline slice
ia.close()
#
hanning(outfile='', region='', mask='', axis=-10, drop=True, overwrite=False, stretch=False, dmethod='copy')[source]

This application performs Hanning convolution of one axis of an image defined by

z[i] = 0.25y[i-1] + 0.5y[i] + 0.25y[i+1] (equation 1)

where z[i] is the value at pixel i in the hanning smoothed image, and y[i-1], y[i], and y[i+1] are the values of the input image at pixels i-1, i, and i+1 respectively. It supports both float and complex valued images. The length of the axis along which the convolution is to occur must be at least three pixels in the selected region. Masked pixel values are set to zero prior to convolution. All nondefault pixel masks are ignored during the calculation.

The convolution is done in the image domain (i.e., not with an FFT).

If drop=False, the length of the output axis will be the same as that of the input axis. The output pixel values along the convolution axis will be related to those of the input values according to equation 1, except the first and last pixels. In that case,

z[0] = 0.5(y[0] + y[1])

and,

z[N-1] = 0.5(y[N-2] + y[N-1])

where N is the number of pixels along the convolution aixs. The pixel mask, ORed with the OTF mask if specified, is copied from the selected region of the input image to the output image. Thus for example, if the selected region in the input image has six planes along the convolution axis, and if the pixel values, which are all unmasked, on a slice along this axis are [1, 2, 5, 10, 17, 26], the corresponding output pixel values will be [1.5, 2.5, 5.5, 10.5, 17.5, 21.5].

If drop=True and dmethod=”copy”, the output image is the image calculated if drop=True, except that only the odd-numbered planes are kept. Furthermore, if the number of planes along the convolution axis in the selected region of the input image is even, the last odd number plane is also discarded. Thus, if the selected region has N pixels along the convolution axis in the input image, along the convolution axis the output image will have (N-1)/2 planes if N is odd, or (N-2)/2 planes if N is even. In this case, the pixel and mask values are copied directly, without further processing. Thus for example, if the selected region in the input image has six planes along the convolution axis, and if the pixel values, which are all unmasked, on a slice along this axis are [1, 2, 5, 10, 17, 26], the corresponding output pixel values will be [2.5, 10.5].

If drop=True and dmethod=”mean”, first the image described in the drop=False case is calculated. The first plane and last plane(s) of that image are then discarded as described in the drop=True, dmethod=”copy” case. Then, the ith plane of the output image is calculated by averaging the (2i)th and (2i + 1)th planes of the intermediate image.Thus for example, if the selected region in the input image has six planes along the convolution axis, and if the pixel values, which are all unmasked, on a slice along this axis are [1, 2, 5, 10, 17, 26], the corresponding output pixel values will be [4.0, 14.0]. Masked values are taken into consideration when forming this average, so if one of the values is masked, it is not used in the average. If at least one of the values in the input pair is not masked, the corresponding output pixel will not be masked.

The hanning smoothed image is written to disk with name outfile, if specified. If not, no image is written but the image is still accessible via the returned image analysis tool (see below).

This method always returns an image analysis tool which is attached to the hanning smoothed image. This tool should always be captured and closed after any desired manipulations have been done. Closing the tool frees up system resources (eg memory), eg,

hanning_image = ia.hanning(…)

# do things (or not) with hanning\_image 
... 
# close the returned tool promptly upon finishing with it.

hanning_image.done()

See also the other convolution functions convolve2d, sepconvolve and convolve.

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • axis (int=-10) - Zero based axis to convolve. ia.coordsys().names() gives the order of the axes in the image. Less than 0 means use the spectral axis if there is one, if not an exception is thrown.

  • drop (bool=True) - Drop every other pixel on output?

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

  • dmethod (string='copy') - If drop=True, method to use in plane decimation. “c(opy)”: direct copy of every second plane, “m(ean)”: average planes 2*i and 2*i+1 in the smoothed, non-decimated image to form plane i in the output image.

Returns

image

Examples

ia.open("mynonsmoothed.im")
# smooth the spectral axis, say it's axis 2 and only write every other pixel
hanning = ia.hanning(outfile="myhanningsmoothed.im", axis=2, drop=True, overwrite=True)
# done with input
ia.done()
# do something with the output image, get statistics say
stats = hanning.statistics()
# close the result image
hanning.done()
haslock()[source]

This function can be used to find out whether the image has a read or a write lock set. It is not of general user interest. It returns a vector of Booleans of length 2. Position 1 says whether a read lock is set, position 2 says whether a write lock is set.

In general locking is handled automatically, with a built in lock release cycle. However, this function can be useful in scripts when a file is being shared between more than one process. See also functions unlock and lock.

histograms(axes=[-1], region='', mask='', nbins=25, includepix=[-1], cumu=False, log=False, stretch=False)[source]

This method computes histograms of the pixel values in the image. The values are returned in a dictionary.

The chunk of the image over which you compute the histograms is specified by a vector of axis numbers (argument axes). For example, consider a 3-dimensional image for which you specify axes=[0,2]. The histograms would be computed for each XZ (axes 0 and 2) plane in the image. You could then examine those histograms as a function of the Y (axis 1) axis. Or perhaps you set axes=[2], whereupon you could examine the histogram for each Z (axis 2) profile as a function of X and Y location in the image.

You have control over the number of bins for each histogram (nbins). The bin width is worked out automatically for each histogram and may vary from histogram to histogram (the range of pixel values is worked out for each chunk being histogrammed).

You have control over which pixels are included in the histograms via the includepix argument. This vector specifies a range of pixel values to be included in the histograms. If you only give one value for this, say includepix=[b], then this is interpreted as includepix=[-abs(b),abs(b)]. If you specify an inclusion range, then the range of pixel intensities over which the histograms are binned is given by this range too. This is a way to make the bin width the same for each histogram.

You can control if the histogram is cumulative or non-cumulative via the cumu parameter.

You have countrol over how the bin counts are returned. If log = False, the actual counts are returned. If True, the values returned are the log10 values of the actual counts.

The results are returned as a dictionary. The counts (field “counts”) and the abscissa values (field “values”) for all bins in each histogram are returned. The shape of the first dimension of those arrays contained in those fields is nbins. The number and shape of the remaining dimensions are those of the display axes(the axes in the image for which you did not compute the histograms). For example, if one has a three dimensional image and sets axes=[2], the display axes are 0 and 1, so the shape of each counts and values array is then [nbins,nx,ny], where nx and ny are the length of the zeroth and first axes, respectively.

In addition, the mean (field “mean”) and standard deviation (field “sigma”) computed using the data in each histogram is returned. The shape of these arrays is equal to the shape of the display axes. So,

Parameters

  • axes (intVec=[-1]) - List of axes to compute histograms over. Default is all axes.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • nbins (int=25) - Number of bins in histograms, > 0

  • includepix (doubleVec=[-1]) - Range of pixel values to include. Default is to include all pixels.

  • cumu (bool=False) - If T the bin values are cumulative.

  • log (bool=False) - If True, the returned counts values will be the log10 values of the actual counts, if False, the actual counts will be returned.

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

Returns

record

Examples

# obtain a histogram using the entire image
ia.maketestimage()
res = ia.histograms()
ia.close()

# obtain histograms for each plane along axis 1 with each
# histogram having 30 bins. Only pixel values in the range
# -0.001 to 0.001 are used in computing the histograms and the
# statistics. The counts in the returned data structure represent
# the cumulative number of data points in the current bin and in
# bins less than the current bin.
ia.open("myimage.im")
r = ia.histograms(axes=[0,2],nbins=30,includepix=1e-3,cumu=True)
ia.close()
history(list=True)[source]

This method interogates the history of an image.

The history is returned as an array of strings, where each element represents an individual history entry. If True, the list parameter will also cause the history to be emitted by the logger.

Note that entries can be permanently added to the image history by using the ia.sethistory() method.

Parameters

  • list (bool=True) - List history to the logger?

Returns

stringVec

Examples

#
print "\t----\t history Ex 1 \t----"
ia.maketestimage()
ia.history()                       # List history to logger
h = ia.history(list=False)             # Recover history in variable h
ia.history(list=True, browse=False)       # List history to logger
#
image()[source]

image method

imagecalc(outfile='', pixels='', overwrite=False, imagemd='', prec='float')[source]

This method is used to evaluate a mathematical expression involving existing images. It fully supports float, double, and complex float, and complex double valued images. The syntax of the expression supplied via the pixels parameter (in what is called the Lattice Expression Language, or LEL) is explained in detail in . This is a rich mathematical language with allows all manner of mathematical operations to be applied to images.

Any image files embedded in the expression may be native  or  (but not yet Miriad) images.

If successful, this method always returns an image analysis tool that references the image resulting from the calculation. This returned tool should always be captured and closed as soon as the user is done with it to free up system resources (eg, memory). The image analysis tool on which the method is called (eg the ia tool when one runs ia.imagecalc()) remains unaltered, eg it still refers to the same image it did prior to the imagecalc() call.

Values of the returned tool are evaluated “on demand”. That is, only when a method is run on the returned tool are the necessary values computed. And in fact, the values have to be reevaluated for each operation (method call). This means that there is a small performance hit for using the returned tool rather than the image written to disk and that none of the images which were used in the expression should be deleted while the returned tool is in use because they must be accessed for calculating the expression each time an operation of the returned tool is performed. These limitations do not apply to the ouput image if one is specified with the outfile parameter; it is a genuine CASA image with numerical values. If outfile is blank, no ouput image is written (although the resulting image can still be accessed via the returned image analysis tool as described below).

Normally you should just write the image, close the returned tool, and open the results image with the default ia tool and operate on it. If you are interested in conserving disk space, you don’t need to keep the result of the calculation around for very long, and/or you are only going to do a small number of operations on the result image, should you set outfile=””.

Note that when multiple image are used in the expression, there is no garauntee about which of those images will be used to create the metadata of the output image, unless imagemd is specified. If imagemd is specified, the following rules of metadata copying will be followed:

  1. The pixel data type of the image specified by imagemd and the output image must be the same. 2. The metadata copied include the coordinate system (and so of course the dimensionality of the output image must correspond to the coordinate system to be copied), the image_info record (which contains things like the beam(s)), the misc_info record (should one exist in the image specified by imagemd), and the units. 3. If the output image is a spectral image, the brightness units are set to the empty string. 4. If the ouptut image is a polarization angle image, the brightness unit is set to “deg” and the stokes coordinate is set to have a single plane of type of Pangle.

The precision (float or double) of the output image pixels is determined by the value specified by the prec parameter. The domain (real or complex) of the output pixel values is determined from the expression.

Parameters

  • outfile (string='') - Output image file name. If blank the resulting image is not written, but it can still be accessed via the returned image analysis tool.

  • pixels (string='') - LEL expression. Must be specified. For example “myimage1.im + myimage2.im”.

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • imagemd (string='') - The image from which metadata should be copied. Default means no gaurantee from which image is used.

  • prec (string='float') - Precision for the output image pixels. ‘float’ or ‘double’ (minimum match supported)

Returns

image

Examples

# Suppose aF and bF are images with single precision and we want
# to determine the result of the following expression:
# aF + min(float($\pi$, mean(bF))
#  
# In this case, the images aF and bF do not need to have the same shapes and
# coordinates, because only the mean(bF) results in the mean of all pixel
# values in bF.  If aF has single precision pixel values, the resulting image
# will as well. This expression first computes the scalar value of the minimum
# of $\pi$ and the mean of the pixel values of bF. That scalar is then
# added to the value of each pixel in aF. In the code below, the
# result is written to image cF which can be accessed immediately
# via the returned image analysis tool captured in the variable myim.
# If the expression is masked, that mask will be copied to the new image.

# create images of different sizes to operate on
ia.fromshape('aF',[10,10],overwrite=True)
ia.fromshape('bF',[10,20,30],overwrite=True)
# close the ia tool to free up resources
ia.done()
# at each pixel in bF, take the minimum of that pixel value and pi and add
# the resulting value to the corresponding pixel in af
# note that only the subset of pixels in bF that correspond to those in aF
# are used; the resulting image has the same size as the smaller image, aF,
# used in the input
myim = ia.imagecalc(outfile='cF', pixels='aF + min(float(pi()), mean(bF))',
                 overwrite=True)
# confirm the resulting image has the same size as aF, should be [10, 10]
myim.shape()
# close the myim tool to free up system resources
myim.done()
imageconcat(outfile='', infiles='', axis=-1, relax=False, tempclose=True, overwrite=False, reorder=False, mode='paged')[source]

This application is used to concatenate two or more input CASA images into one output image. For example, if you have two image cubes which are contiguous along one axis (say a spectral axis) and you would like to glue them together along this axis, then this function is the appropriate thing to use.

The axis parameter is used to specify which zero-based axis the images should be concatenated along. A negative value indicates that the spectral axis should be used. If a negative value is given but there is no spectral axis, an exception will be thrown. The zero-based order of the axes of an image can be determined from ia.coordsys().names().

If successful, this method will return an image analysis tool referencing the concatenated image. Even if it is not wanted, the returned tool should be captured and closed as soon as the user is finished with it to free up system resources (eg memory).

If outfile is given, the image is written to the specified disk file. If outfile is unset, the on-the-fly image tool created by the method actually references all of the input files. So if you deleted any of the input image disk files, it would render this tool useless. Note that a blank outfile name may only be given in the case where mode=’p’ (see below). When you destroy this tool done() method, the reference connections are broken.

The input and output images must be of the same dimensionality. Therefore, if you wish to concatenate 2-D images into a 3-D image, the 2-D images must have a third axis (of length unity) so that the output image coordinates are known along the concatenation axis.

The input images are concatenated in the order in which they are listed unless the reorder parameter is set to True. If True, the images are reordered if necessary so that the world coordinate values along the selected axis monotonically increase or decrease. The direction of the increment is determined by the first listed image. If reorder=True, the world coordinate ranges of the images along the selected axis are not permitted to overlap, and the signs of the increments for this axis in all images must be the same. If reorder=False, the coordinate system of the first listed image is used as the coordinate system for the output image. If reorder=True, the coordinate system of the first image in the list of the reordered images is used as the coordinate system of the output image. Setting reorder=True can be especially useful if the infiles are specified using a wildcard character(s).

If relax=False, the input images are checked to see that they are contiguous along the concatenation axis and an exception is thrown if they are not. In addition, the coordinate descriptors (e.g. reference pixel, reference value etc) for the non-concatenation axes must be the same or an error will result.

The input disk image files may be in native (CASA or FITS) formats.

The contiguous criterion and coordinate descriptor equality criteria can be relaxed by setting relax=True whereupon only warnings will be issued. Dimension and shape must still be the same though. When the concatenation axis is not contiguous (but still monotonically increasing or decreasing) and relax=True, a tabular coordinate will be used to correctly describe the axis. But be aware that it means adjacent pixels are not regularly spaced. However, methods like toworld() and topixel() will correctly interconvert world and pixel coordinates.

In giving the input image names, the infiles argument can be a single string if you wild card it with standard shell symbols. For example, infiles=’cena_???.*’, where the ‘?’ represents one character and ‘*’ any number of characters.

Otherwise, you must input a vector of strings such as infiles=”cena1 cena2 cena3”. An input such as infiles=’files1,file2’ will be interpreted as one string naming one file and an exception will be thrown. The reason for this is that although the latter could be parsed to extract two file names by recognizing comma delimiters, it is not possible because an expression such as infiles=’cena.{a,b}’ (meaning files of name ‘cena.a’ and ‘cena.b’) would confuse such parsing (you would get two files of name cena.{a} and {sff b}.

You can look at the coordinate system of the output image using the ia.summary() method to ensure it’s correct.

The parameter tempclose is, by default, True. This means that all internal reference copies of the input images are kept closed until they are needed. Then they are opened temporarily and then closed again. This enables you to effectively concatenate as many images as you like without encountering any operating system open file number limits. However, it comes at some performance loss, because opening and closing all those files takes time. If you are concatenating a smallish number of files, you might use tempclose=False. This will leave all internal reference copies permanently open, but performance, if you don’t hit the file limit, will be better.

The mode parameter controls the format of the output image. If mode=’p’, then a “paged” (normal persistent CASA format) image in which all the pixel values and metadata are copied to a single image is created. All the other modes will create a “virtual concat” image. In this format, a directory, named outfile, is created to store metadata about each image. The metadata are placed in a file in this directory named imageconcat.json, which is a plain text file in JSON format. In the json file are pointers to the actual disk images that make up the virtual concat image. In the case of mode=’c’ (copy), alll the disk images are copied to the outfile directory and the json file references the images in the outfile directory. In the case of mode=’m’ (move), the disk images are instead moved to the outfile directory, and again the json file references the images in the outfile directory. These two modes make it simple to deal with the virtual concat image as a single unit, for example, if the user needs to move the image or tar it up, the resulting moved image or untarred file will be accessible by the image tool as a valid image. The final option, mode=’n’ (no copy) means the output directory only contains the json file and references the input image paths as given. In this case, moving the outfile directory without moving the reference images will normally cause problems, especially if the input images were given using relative path names. So, one must think carefully if mode=’n’ will be the correct option for their use case in the long run.

This method requires multiple images which are specified with the infiles parameter. Therefore calling ia.open() is not necessary, although calling imageconcat() using an already open image analysis tool will work and the state of that tool (eg the image it references) will not be changed.

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • infiles (variant='') - List of input casaimage files to concatenate; wild cards accepted. Default is empty string.

  • axis (int=-1) - Concatenation pixel axis. Use ia.coordsys().names() to get a list of axes. A negative value means use the spectral axis if there is one, if not an exception is thrown.

  • relax (bool=False) - Relax constraints that axis coordinate descriptors match

  • tempclose (bool=True) - Keep all lattices closed until needed

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • reorder (bool=False) - Automatically reorder the images if necessary.

  • mode (string='paged') - Type of output image. “paged”, “movevirtual”, “nomoveevirtual”, “copyvirtual”. Minimal match supported.

Returns

image

Examples

# Create three images to concatenate together.
ia.fromshape('im.1',[10,10,10],overwrite=True)
ia.fromshape('im.2',[10,10,10],overwrite=True)
ia.fromshape('im.3',[10,10,10],overwrite=True)
ia.done()
# now concatenate.
# The three images have the same shape along the axes not to be
# concatenated as they must. relax=True means that the contiguity
# constraint along the concatenated axis is not imposed (if it were
# the call would fail because the spectral axes of the input images
# are not contiguous).
bigim = ia.imageconcat(outfile='bigimage', infiles='im.1 im.2 im.3',
                       axis=2, relax=True, tempclose=False, overwrite=True)
# be sure to call done() on the return tool to free up system resources.
bigim.done()
insert(infile='', region='', locate=[-1], verbose=False)[source]

This function inserts the specified image (or part of it) into the image referenced by this tool. The specified image may be given via argument infile as a disk file name (it may be in native , , or Miriad format; Look for more information on foreign images).

If the locate vector is not given, then the images are aligned (to an integer pixel shift) by their reference pixels.

If locate vector is given, then those values that are given, give the absolute pixel in the output (this) image of the bottom left corner of the input (sub)image. For those values that are not given, the input image is symmetrically placed in the output image.

The image referenced by this tool is modified in place; no new image is created. The method returns True if successful.

Parameters

  • infile (string='') - Name of image to be inserted.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • locate (doubleVec=[-1]) - Location of input image in output image. Default is centrally located.

  • verbose (bool=False) - Emit informational messages to logger?

Returns

bool

Examples

#
print "\t----\t insert Ex 1 \t----"
ia.maketestimage('myfile.insert',overwrite=True)
ia.close()
ia.fromshape(shape=[200,200])
ia.insert(infile='myfile.insert')       # Align by reference pixel
ia.newimagefromfile('myfile.insert')
ia.insert(infile=im2.name(), locate=[]) # Align centrally
# This time align axis 0 as given and axis 1 centrally
ia.insert(infile='myfile.insert', locate=[20])
ia.close()                                    # close default tool and
isconform(other='')[source]

Returns True if the shape, coordinate system, and axes order of the specified image matches the current image.

Parameters

  • other (string='') - name of image to test

Returns

bool

Examples

ia.isconform("my_mystery.image")
isopen()[source]

This method returns True if the image tool has an attached image.

ispersistent()[source]

This function can be used to find out whether the image is persistent on disk or not. There is a subtle difference from the image being virtual. For example, a virtual image which references another which is on disk is termed persistent.

lock(writelock=False, nattempts=0)[source]

This function can be used to acquire a Read or a Read/Write lock on the . It is not of general user interest.

In general locking is handled automatically, with a built in lock release cycle. However, this function can be useful in scripts when a file is being shared between more than one process. See also functions unlock and haslock.

Parameters

  • writelock (bool=False) - Acquire a read/write (T) or a readonly (F) lock

  • nattempts (int=0) - Number of attempts, > 0. Default is unlimiited.

Returns

bool

Examples

#
print "\t----\t lock Ex 1 \t----"
ia.maketestimage('xx', overwrite=True)
ia.lock(writelock=True)
ia.unlock()
ia.lock(writelock=False)
ia.close(remove=True)
#


This acquires a read/write lock on the file. Then we unlock it
and acquire a readonly lock.
makearray(v=0.0, shape=[0])[source]

This function takes two arguments. The first argument is the initial value for the new array. The second is a vector giving the lengths of the dimensions of the array.

Parameters

  • v (double=0.0) - Value with which to initial array elements

  • shape (intVec=[0]) - Vector containing array dimensions.

Returns

any

Examples

A three dimensional array that is initialized to all zeros. Each of
the three dimensions of the cube has a length of four.


#
print "\t----\t makearray Ex 1 \t----"
cube = ia.makearray(0,[4,4,4])
#
makecomplex(outfile='', imag='', region='', overwrite=False)[source]

This function combines the current image with another image to make a complex image. The current image (i.e. that associated with this Image must have real valued pixels). The image used for generating the imaginary part of the pixel values is specified using the imag parameter, and it must persistent. The image attached to this tool and the image specified using the imag parameter must have the same precision, or else an exception will be thrown. If both are float precision, the resulting image will have float precision pixel values. If both are double precision, the resulting image will have double precision pixel values. The coordinate systems of the two input images must be conformant. The metadata written to the resulting image is copied from the image attached to this tool.

Parameters

  • outfile (string='') - Output Complex (disk) image file name

  • imag (string='') - Imaginary image file name

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

Returns

bool

Examples

#
print "\t----\t makecomplex Ex 1 \t----"
ia.maketestimage('imag.im',overwrite=True)  #imaginary image
ia.close()
ia.maketestimage('real.im',overwrite=True)  #assoc. real image
ia.makecomplex('complex.im', 'imag.im', overwrite=True)
ia.close()
#
maketestimage(outfile='', overwrite=False)[source]

This function converts a FITS file resident in the  system into a  image.

If outfile is given, the image is written to the specified disk file. If outfile is unset, the Image tool is associated with a temporary image. This temporary image may be in memory or on disk, depending on its size. When you close the Image tool (with the close function) this temporary image is deleted.

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

Returns

bool

Examples

#
print "\t----\t maketestimage Ex 1 \t----"
ia.maketestimage()     # make virtual image
ia.close()
ia.maketestimage('tmp', overwrite=True)
ia.close()             # close to unlock disk image
#
maskhandler(op='default', name='')[source]

This method is used to manage or handle pixel masks . A CASA image may contain zero, one, or more pixel masks. Any of these masks can be designated the default pixel mask assoicated with the given image. The default mask is acted upon and/or used by CASA applications. For example, if ia.statistics() will exclude pixels which are masked as bad (False) from the calculations.

This method does not modify the individual boolean values of any masks.

The op parameter is used to specify the behaviour. In all cases, the specified operation can be specified by a three character string. Supported values of the op parameter are:

’default’: Retrieve the name of the default pixel mask associated with the image. A one element array containing the empty string is returned if the image has no default mask.

’get’: Retrieve the name(s) of all the existing pixel masks. Note that the ia.summary() method may also be used to view the pixel masks associated with an image.

’set’: Set the default pixel mask to the value specified by the name parameter. If the value of the name parameter is the empty string, then the default mask is unset (ie, all the pixels will be treated as being unmasked).

’delete’: Delete the pixel mask(s) specified by the name parameter. To delete more than one mask, the name parameter can be an array of strings. Any supplied pixel mask name that does not exist is silently ignored.

’rename’: Rename the mask specified as the first element of the name array parameter (name[0]) to the value specified in the second element of the name array parameter (name[1]]. In this case, the name array parameter must have exactly two elements.

’copy’: Copy the mask specified in the first element of the name array parameter (name[0]) to the a mask name specified in the second element of the name array parameter (name[1]]. In this case, the name array parameter must have exactly two elements. A mask from another image can be copied by using the imagename:maskname syntax for the first element in the name array, eg, ’myimage:mask2’.

Parameters

  • op (string='default') - The operation. One of ‘set’, ‘delete’, ‘rename’, ‘get’, ‘copy’ or ‘default’

  • name (stringVec='') - Name of mask or masks.

Returns

stringVec

Examples

#
print "\t----\t maskhandler Ex 1 \t----"
ia.maketestimage('g1.app', overwrite=True)
ia.calcmask('T', name='mask1')
ia.close()
ia.maketestimage('myimage', overwrite=True)
ia.calcmask('T')                         # Create some masks
ia.calcmask('T', name='mask1')
ia.calcmask('T', name='mask2')
names = ia.maskhandler('get')            # Get the mask names
print names
#['mask0', 'mask1', 'mask2']
name = ia.maskhandler('default')         # Get the default mask name
print name
#mask2
ia.maskhandler('set', ['mask1'])         # Make 'mask1' the default mask
ia.maskhandler('set', [''])              # Unset the default mask
ia.maskhandler('delete', ['mask1'])      # Delete 'mask1'
ia.calcmask('T', name='mask1')           # Make another 'mask1'
ia.maskhandler('delete', ['mask0', 'mask1'])# Delete 'mask0' and 'mask1'
ia.calcmask('T', name='mask1')
ia.maskhandler('rename', ['mask1', 'mask0'])# Rename 'mask1' to 'mask0'

# Copy 'mask1' from image 'g1.app'  to 'mask10' in image 'myimage'
ia.maskhandler('copy', ['g1.app:mask1', 'mask10'])
ia.removefile('g1.app')                  # Cleanup
ia.close()
#
maxfit(region='', point=True, width=5, negfind=False, list=True)[source]

This function finds the pixel with the maximum value in the region, and then uses function findsources to generate a Componentlist with one component. The component will be of type Point (point=T) or Gaussian (point=F).

If negfind=F the maximum pixel value is found in the region and fit. If negfind=T the absolute maximum pixel value is found in the region and fit.

See function findsources for a description of arguments point and width.

See also the function fitcomponents.

Parameters

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • point (bool=True) - Find only point sources?

  • width (int=5) - Half-width of fit grid when point=F

  • negfind (bool=False) - Find negative sources as well as positive?

  • list (bool=True) - List the fitted parameters to the logger?

Returns

record

Examples

#
print "\t----\t maxfit Ex 1 \t----"
ia.maketestimage()
clrec = ia.maxfit()
print clrec          # There is only one component
ia.close()
#
miscinfo()[source]

A   can accumulate miscellaneous information during its lifetime. This information is stored in a record called the miscinfo record. For example, the  filler puts header keywords it doesn’t otherwise use into the miscinfo record. This miscinfo record is not guaranteed to have any entries, so it’s up to you to check for any fields that you require.

You can also put things into this record (see setmiscinfo) yourself, to keep information that the system might not otherwise store for you.

When the image is written out to , the items in the miscinfo record are written to the  file as keywords with the corresponding record field name.

modify(model='', region='', mask='', subtract=True, list=True, stretch=False)[source]

This function applies a model of the sky to the image. You can add or subtract the model which is contained in a Componentlist tool.

The pixel values are only changed where the total mask (combination of the default pixel mask [if any] and the OTF mask) is good (True). If the computation fails for a particular pixel (e.g. coordinate undefined) that pixel will be masked bad.

DISK MODELS

Pixels with centers inside the disk will have the same values, even if a pixel straddles the edge of the disk. Pixels with straddle the edge of the disk which have centers outside the disk are given values of zero. Thus, one should not expect the flux density of the disk to be exactly the provided value to the component list; for a given size disk, the computed flux density will be closer to the expected value for images with smaller pixels.

Parameters

  • model (record='') - Record representation of a ComponentList model

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • subtract (bool=True) - Subtract or add the model

  • list (bool=True) - List informative messages to the logger

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

Returns

bool

Examples

#
print "\t----\t modify Ex 1 \t----"
ia.maketestimage()
clrec = ia.fitcomponents()
ia.modify(clrec['results'])
ia.close()
#





In this example we subtract the model returned by the fitcomponents function.
moments(moments=[0], axis=-10, region='', mask='', method='', smoothaxes=[-1], smoothtypes='', smoothwidths=[0.0], includepix=[-1], excludepix=[-1], peaksnr=3.0, stddev=0.0, doppler='RADIO', outfile='', smoothout='', overwrite=False, drop=True, stretch=False)[source]

The primary goal of this function is to enable you to analyze a multi-dimensional image by generating moments of a specified axis. This is a time-honoured spectral-line analysis technique used for extracting information about spectral lines.

You can generate one or more output moment images. The return value of this function is an on-the-fly Image  holding the first of the output moment images.

The word ’moment’ is used loosely here. It refers to collapsing an axis (the moment axis) to one pixel and setting the value of that pixel (for all of the other non-collapsed axes) to something computed from the data values along the moment axis. For example, take an RA-DEC-Velocity cube, collapse the velocity axis by computing the mean intensity at each RA-DEC pixel. This function offers many different moments and a variety of automatic methods to compute them.

We try to make a distinction between a ’moment’ and a ’method’. This boundary is a little blurred, but it claims to refer to the distinction between what you are computing, and how the pixels that were included in that computation were selected. For example, a ’moment’ would be the average value of some pixel values in a spectrum. A ’method’ for selecting those pixels would be a simple pixel value range specifying which pixels should be included.

There are many available moments, and you specify each one with an integer code as it would get rather cumbersome to refer to them via strings. In the list below, the value of the \(i\)th pixel of the spectrum is \(I_i\), the coordinate of this pixel is \(v_i\) (of course it may not be velocity), and there are \(n\) pixels in the spectrum. The available moments are:

  • \(-1\) – the mean value of the spectrum

    \[{ {1\over n} {\sum {I_i}}}\]
  • 0 – the integrated value of the spectrum

    \[M\_0 = \Delta v \sum I\_i\]

    where \(\Delta v\) is the width (in world coordinate units) of a pixel along the moment axis

  • 1 – the intensity weighted coordinate (this is traditionally used to get ’velocity fields’)

    \[M\_1 = { {\sum {I_i v\_i}} \over {M_0}}\]
  • 2 – the intensity weighted dispersion of the coordinate (this is traditionally used to get ’velocity dispersion fields’)

    \[\sqrt{ { {\sum {I_i \left(v_i - M\_1\right)^2}} \over {M_0}}}\]
  • 3 – the median of \(I\)

  • 4 – the median coordinate. Here we treat the spectrum as a probability distribution, generate the cumulative distribution, and then find the coordinate corresponding to the 50% value. This moment is not very robust, but it is useful for quickly generating a velocity field in a way that is not sensitive to noise. However, it will only give sensible results under certain conditions. The generation of the cumulative distribution and the finding of the 50% level really only makes sense if the cumulative distribution is monotonic. This essentially means only selecting pixels which are positive or negative. For this reason, this moment type is only supported with the basic method (see below – i.e. no smoothing, no windowing, no fitting) with a pixel selection range that is either all positive, or all negative

  • 5 – the standard deviation about the mean of the spectrum

    \[\sqrt{ {1\over {\left(n-1\right)}} \sum{\left(I_i - \bar{I}\right)^2 }}\]
  • 6 – the root mean square of the spectrum

    \[\sqrt{ {1 \over n} \sum{I_i^2}}\]
  • 7 – the absolute mean deviation of the spectrum

    \[{1 \over n} \sum {|(I_i - \bar{I})|}\]
  • 8 – the maximum value of the spectrum

  • 9 – the coordinate of the maximum value of the spectrum

  • 10 – the minimum value of the spectrum

  • 11 – the coordinate of the minimum value of the spectrum

The purpose of the smoothing functionality is purely to provide a mask. Thus, you can smooth the input image, apply a pixel include or exclude range, and generate a smoothed mask which is then applied before the moments are generated. The smoothed data are not used to compute the actual moments; that is always done from the original data.

The basic method is to just compute moments directly from the pixel values. This can be modified by applying pixel value inclusion or exclusion ranges (arguments includepix and excludepix).

You can then also convolve the image (arguments smoothaxes, smoothtypes, and smoothwidths) and find a mask based on the inclusion or exclusion ranges applied to the convolved image. This mask is then applied to the unsmoothed data for moment computation.

The window method (invoked with argument method=’window’) does no pixel-value-based selection. Instead a window is found (hopefully surrounding the spectral line feature) and only the pixels in that window are used for computation. This window can be found from the convolved or unconvolved image (arguments smoothaxes, smoothtypes, and smoothwidths).

The moments are always computed from the unconvolved data. The window can be found (for each spectrum) automatically. The automatic methods are via Bosma’s converging mean algorithm (method=’window’) or by fitting Gaussians and taking \(\pm 3\sigma\) as the window (method=’window,fit’). In Bosma’s algorithm, an initial guess for a range of pixels surrounding a spectral feature is refined by widening until the mean of the pixels outside of the range converges (to the noise).

The fit method (method=’fit’) fits Gaussians to spectral features automatically. The moments are then computed from the Gaussian fits (not the data themselves).

  • outfile - If you are creating just one moment image, and you specify outfile, then the image is created on disk with this name. If you leave outfile empty then a temporary image is created. In both cases, you can access this image with the returned Image . If you are making more than one moment image, then theses images are always created on disk. If you specify outfile then this is the root for the output file names. If you don’t specify it, then the input image name is used as the root.

  • smoothing - If you smooth the image to generate a mask, you specify the kernel widths via the smoothwidths argument in the same way as in the sepconvolve function. See it for details.

  • stddev - Some of the automatic methods also require an estimate of the noise level in the image. This is used to assess whether a spectrum is purely noise or not, and whether there is any signal worth digging out. If you don’t give it via the stddev argument, it will be worked out automatically from a Gaussian fit to the bins above 25% from a histogram of the entire image.

  • includepix, excludepix - The vectors given by arguments includepix and excludepix specify a range of pixel values for which pixels are either included or excluded. They are mutually exclusive; you can specify one or the other, but not both. If you only give one value for either of these, say includepix=b, then this is interpreted as includepix=[-abs(b),abs(b)].

    The convolving point-spread function is normalized to have a volume of unity. This means that point sources are depressed in value, but extended sources that are large with respect to the PSF remain essentially on the same intensity scale; these are the structures you are trying to find with the convolution so this is what you want. If you convolve the image, then arguments like includepix select based upon the convolved image pixel values. If you are having trouble getting these right, you can output the convolved image (smoothout) and assess the validity of your pixel ranges. Note also that if you are Hanning convolving (usually used on a velocity axis), then the width for this kernel must be 3 pixels (triangular smoothing kernels of other widths have no valid theoretical basis).

  • doppler - If you compute the moments along a spectral axis, it is conventional to compute the world coordinate (needed for moments 0, 1 and 2) along that axis in “km/s”. The argument doppler lets you specify what doppler convention the velocity will be calculated in. You can choose from doppler=radio, optical, True. See function summary for the definitions of these codes. For other moment-axis types, the world coordinate is computed in the native units.

  • mask - The total input mask is the combination of the default  (if any) and the OTF mask. Once this mask has been established, then the moment method may make additional pixel selections.

  • drop - If this is True (the default) then the moment axis is dropped from the output image. Otherwise, the output images have a moment axis of unit length and coordinate information that is the same as for the input image. This coordinate information may be totally meaningless for the moment images.

Finally, if you ask for a moment which requires the coordinate to be computed for each profile pixel (these are the intensity weighted mean coordinate [moment 1] and the intensity weighted dispersion of the coordinate [moment 2]), and the profile axis is not separable then there will be a performance loss. Examples of non-separable axes are RA and Dec. If the axis is separable (e.g. a spectral axis) there is no penalty. In the latter case, the vector of coordinates for one profile is the same as the vector for another profile, and it can be precomputed (once).

Note that this function has no “virtual” output file capability. All output files are written to disk. The output mask for these images is good (T) unless the moment method fails to generate a value (e.g. the total input pixel mask was all bad for the profile) in which case it will be bad (F).

If an image has multiple (per-channel beams) and the moment axis is equal to the spectral axis, each channel will be convolved with a beam that is equal to the beam having the largest area in the beamset prior to moment determination.

Parameters

  • moments (intVec=[0]) - List of moments that you would like to compute. Default is integrated spectrum.

  • axis (int=-10) - The moment axis. Default is the spectral axis if there is one.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • method (stringVec='') - List of windowing and/or fitting functions you would like to invoke. Vector of strings from ‘window’ and ‘fit’. The default is to not invoke the window or fit functions.

  • smoothaxes (intVec=[-1]) - List of axes to smooth. Default is no smoothing.

  • smoothtypes (variant='') - List of smoothing kernel types, one for each axis to smooth. Vector of strings from ‘gauss’, ‘boxcar’, ‘hanning’. Default is no smoothing.

  • smoothwidths (doubleVec=[0.0]) - List of widths (full width for boxcar, full width at half maximum for gaussian, 3 for Hanning) in pixels for the smoothing kernels. Vector of numeric. Default is no smoothing.

  • includepix (doubleVec=[-1]) - Range of pixel values to include. Vector of 1 or 2 doubles. Default is include all pixel.

  • excludepix (doubleVec=[-1]) - Range of pixel values to exclude. Default is exclude no pixels.

  • peaksnr (double=3.0) - The SNR ratio below which the spectrum will be rejected as noise (used by the window and fit functions only)

  • stddev (double=0.0) - Standard deviation of the noise signal in the image (used by the window and fit functions only)

  • doppler (string='RADIO') - Velocity doppler definition for velocity computations along spectral axes

  • outfile (string='') - Output image file name (or root for multiple moments). Default is input + an auto-determined suffix.

  • smoothout (string='') - Output file name for convolved image. Default is don’t save the convolved image.

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • drop (bool=True) - Drop moments axis from output images?

  • stretch (bool=False) - Stretch the mask if necessary and possible?

Returns

image

Examples

#
print "\t----\t moments Ex 1 \t----"
ia.fromshape(shape=[32,32,32,32]) # replace with your own cube
im2 = ia.moments(moments=[-1,1,2], axis=2, smoothaxes=[0,1,2],
                 smoothtypes=["gauss","gauss","hann"],
                 smoothwidths=[5.0,5.0,3], excludepix=[1e-3],
                 smoothout='smooth', overwrite=True)
im2.done()
ia.close()
#



In this example, standard moments (average intensity, weighted velocity
and weighted velocity dispersion) are computed via the convolve (spatially
convolved by gaussians and spectrally by a Hanning kernel) and clip
method (we exclude any pixels with absolute value less than $0.001$).
The output file names are automatically created for us and
the convolved image is saved.   The returned image tool holds the first
moment image.
name(strippath=False)[source]

This function returns the name of the  By default, this function returns the full absolute path of the . You can strip this path off if you wish with the strippath argument and just recover the  name itself.

Parameters

  • strippath (bool=False) - Strip off the path before the actual file name?

Returns

string

Examples

#
print "\t----\t name Ex 1 \t----"
ia.maketestimage('g1.app', overwrite=True)
print ia.name(strippath=False)
#/casa/code/xmlcasa/implement/images/scripts/g1.app
print ia.name(strippath=True)
#g1.app
ia.close()
#
newimage(infile='')[source]

This method is identical to ia.newimagefromfile(). The description of how it works is in the help for that method.

Parameters

  • infile (string='') - Input image file name

Returns

image

newimagefromarray(outfile='', pixels='', csys='', linear=False, overwrite=False, log=True, type='f')[source]

This application converts a numpy array of any size into a CASA image.

If outfile is specified, the image is written to the specified (persistent) disk file. If outfile is unset, the returned image tool is associated with a temporary image. This temporary image may be in memory or on disk, depending on its size. In this case, when the close() or done() method is called on the returned image tool, the associated temporary image is deleted.

The type parameter controls the data type/precision of the pixel values of the created image. ’f’ indicates that float precision point (32 bit precision) pixel values should be writted. ’d’ indicates that double precision (64 bit precision) pixel values should be written. If the input array has complex (as opposed to real) values, then complex pixel values, with each of the real and imaginary parts having the specified precision, will be written. Array values will be cast automatically to the specified precision, so that the precision of the input array values may be increased, decreased, or unchanged depending on the input array type.

The coordinate system, provided as a a dictionary (use eg, cs.torecord() to do that), is optional. If specified, it must have the same number of dimensions as the pixels array. Call the naxes() method on the coordinate system tool to see how many dimensions the coordinate system has. A coordinate system can be created from scratch using the coordinate system (cs) tool and methods therein, but often users prefer to use a coordinate system from an already existing image. This can be gotten using ia.coordsys() which returns a coordinate system tool. A torecord() call on that tool will result in a python dictionary describing the coordinate system which is the necessary format for the csys input parameter of ia.fromarray().

If csys is not specified, a default coordinate system will be created. If linear=False (the default), the created coordinate system will have standard RA/DEC/Stokes/Spectral Coordinate axes depending upon the shape of the pixels array (Stokes axis must be no longer than 4 pixels and the spectral axis may precede the Stokes axis if eg, shape=[64,64,32,4]. Extra dimensions are given linear coordinates. If linear=True, then all the resulting coordinates are linear with the axes represent lengths. In this case each axis will have a value of 0.0 at its center pixel. The increment of each axis will be 1.0 km.

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • pixels (variant='') - A numeric array is required.

  • csys (record='') - Coordinate System. Default is unset.

  • linear (bool=False) - Make a linear Coordinate System if csys not given

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • log (bool=True) - Write image creation messages to logger

  • type (string='f') - Pixel data type to write. ‘f’ (float precision) or ‘d’ (double precision)

Returns

image

Examples

#
print "\t----\t newimagefromarray Ex 1 \t----"
im1=ia.newimagefromarray(outfile='test.data',
                         pixels=ia.makearray(0, [64, 64, 4, 128]),
                         overwrite=True)
cs1 = im1.coordsys(axes=[0,1])
im1.done()
im2 = ia.newimagefromarray(pixels=ia.makearray(1.0, [32, 64]),
                           csys=cs1.torecord())
cs1.done()
im2.done()
#



The first example creates a zero-filled \imagefile\ named {\sff
test.data} which is of shape [64,64,4,128].  If you examine the header
with {\stff ia.summary()} you will see the default
RA/DEC/Stokes/Frequency coordinate information.  In the second
example, a Coordinate System describing the first two axes of the
image {\sff test.data} is created and used to create a 2D image
temporary image.
newimagefromfile(infile='')[source]

This method returns an image analysis tool associated with the specified image. Constructing a image analysis tool in addition to the default ia tool allows the user to operate on multiple images simultaneously. All ia.newimagefrom() methods share this behavior.

The parameter infile may refer to a CASA image, a Miriad image, or a FITS image. FITS images of types Float are supported.

When finished with the newly created tool, the user should close it to free up system resources (eg memory).

ia.newimage() is an alias for this method.

Parameters

  • infile (string='') - Input image file name

Returns

image

Examples

# This is one way to copy a FITS image into an already extant CASA image
# of the same shape (ia.subimage() is more effecient, but this example is
# meant to demonstrate ia.newimagefromfile()

# note that the ia tool is not attached to an image after the first command,
# the fitsimage tool is
fitsimage = ia.newimagefromfile("myimage.fits")
# now attach the target CASA image to the ia tool
ia.open("myimage.im")
# copy pixel values
ia.putchunk(fitsimage.getchunk())
# copy the coordinate system
ia.setcoordsys(fitsimage.coordsys().torecord())
# copy other miscellaneous things
ia.setbrightnessunit(fitsimage.getbrightnessunit())
ia.setmiscinfo(fitsimage.miscinfo())
# be sure to call done() on both tools to free up memory
ia.done()
fitsimage.done()
newimagefromfits(outfile='', infile='', whichrep=0, whichhdu=0, zeroblanks=False, overwrite=False)[source]

This function is used to convert a FITS disk image file (Float, Double, Short, Long are supported) to an  . If outfile is given, the image is written to the specified disk file. If outfile is unset, the on-the-fly Image  returned by this function is associated with a temporary image. This temporary image may be in memory or on disk, depending on its size. When you destroy the on-the-fly Image  (with the done function) this temporary image is deleted.

This function reads from the FITS primary array (when the image is at the beginning of the FITS file; whichhdu=0), or an image extension (when the image is elsewhere in the FITS file, whichhdu \(\>\) 0).

By default, any blanked pixels will be converted to a mask value which is False, and a pixel value that is NaN. If you set zeroblanks=T then the pixel value will be zero rather than NaN. The mask will still be set to False. See the function replacemaskedpixels if you need to replace masked pixel values after you have created the image.

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • infile (string='') - Input FITS disk file name. Required.

  • whichrep (int=0) - If this FITS file contains multiple coordinate representations, which one should we read

  • whichhdu (int=0) - If this FITS file contains multiple images, which one should we read (0-based).

  • zeroblanks (bool=False) - If there are blanked pixels, set them to zero instead of NaN

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

Returns

image

Examples

#
print "\t----\t newimagefromfits Ex 1 \t----"
# Assume we can find test fits file using
# CASAPATH environment variable
pathname=os.environ.get("CASAPATH")
pathname=pathname.split()[0]
datapath=pathname+'/data/demo/Images/imagetestimage.fits'
im1=ia.newimagefromfits('./myimage', datapath, overwrite=True)
print im1.summary()
print im1.miscinfo()
print 'fields=', im1.miscinfo().keys()
im1.done()
#



The FITS image is converted to a \casa\ \imagefile\ and access is
provided via the \imagetool\ called {\stf im1}.  Any FITS header
keywords which were not recognized or used are put in the
miscellaneous information bucket accessible with the miscinfo function.  In
the example we list the names of the fields in this record.
newimagefromimage(infile='', outfile='', region='', mask='', dropdeg=False, overwrite=False)[source]

This function applies a  to a disk , creates a new  containing the (sub)image, and associates a new  with it.

The input disk image file may be in native , (Float, Double, Short, Long are supported), or Miriad format. Look for more information on foreign images.

If outfile is given, the (sub)image is written to the specified disk file.

If outfile is unset, the Image  actually references the input image file. So if you deleted the input image disk file, it would render this  useless. When you destroy this on-the-fly  (with the done function) the reference connection is broken.

Sometimes it is useful to drop axes of length one (degenerate axes). Use the dropdeg argument if you want to do this.

The output mask is the combination (logical OR) of the default input  (if any) and the OTF mask. Any other input  will not be copied. Use function maskhandler if you need to copy other masks too.

See also the subimage function.

Parameters

  • infile (string='') - Input image file name. Required.

  • outfile (string='') - Output sub-image file name. Default is unset.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • dropdeg (bool=False) - Drop degenerate axes

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

Returns

image

Examples

#
print "\t----\t newimagefromimage Ex 1 \t----"
ia.maketestimage('test1',overwrite=True)
ia.maketestimage('test2',overwrite=True)
print ia.name()
#[...]test2
im1=ia.newimagefromimage('test1')
print im1.name()
#[...]test1
print im1.summary()
im2=ia.newimagefromimage('test2')
print im2.name()
#[...]test2
print im1.name()
#[...]test1
ia.close()
im1.done()
im2.done()
#
newimagefromshape(outfile='', shape=[0], csys='', linear=False, overwrite=False, log=True, type='f')[source]

This function creates a CASA image with the specified shape. It is similar to ia.fromshape(), but instead returns a new image analysis tool attached to the new image, rather than attaching the new image to the current tool. All the pixel values in the image are set to 0. One may create an image with float precision pixels (type=’f’), complex float precision pixels (type=’c’), double precision pixels (type=’d’), or complex double precision pixels (’cd’). To use a numpy array of values to create an image, use ia.(newimage)fromarray(). To make a 2-D image from a packaged FITS file, use ia.maketestimage().

If outfile is given, the image is written to the specified disk file. If outfile is unset, the image analysis tool is associated with a temporary image. This temporary image may be in memory or on disk, depending on its size. When you close the image analysis tool (with the ia.close() method, the temporary image is deleted.

The coordinate system, provided as a coordinate system tool record, is optional. If provided, it must be dimensionally consistent with the specified shape.

If the coordinate system is not provided, a default coordinate system will be created. If linear=False (the default), then it is a standard RA/DEC/Stokes/Spectral coordinate system depending exactly upon the shape (the Stokes axis must be no longer than 4 pixels and spectral axis may occur prior to the Stokes axis if eg, shape=[64,64,32,4]. Extra dimensions are given linear coordinates. If linear=True, then the coordinate system will have linear coordinates.

Parameters

  • outfile (string='') - Name of output image file. Default is unset.

  • shape (intVec=[0]) - Shape of image. Required.

  • csys (record='') - Record describing Coordinate System. Default is unset.

  • linear (bool=False) - Make a linear Coordinate System if csys not given?

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • log (bool=True) - Write image creation messages to logger

  • type (string='f') - Type of image. ‘f’ means Float, ‘c’ means complex.

Returns

image

Examples

#
print "\t----\t newimagefromshape Ex 1 \t----"
im1=ia.newimagefromshape('test2.data', [64,64,128], overwrite=True)
cs1 = im1.coordsys(axes=[0,2])
im1.done()
im2 = ia.newimagefromshape(shape=[10, 20], csys=cs1.torecord())
cs1.done()
im2.done()
#



The first example creates a zero-filled \imagefile\ named {\sff
test.data} of shape [64,64,128].  If you examine the header with
{\stff ia.summary()} you will see the RA/DEC/Spectral coordinate
information.  In the second example, a Coordinate System describing
the first and third axes of the image {\sff test2.data} is created and
used to create a 2D temporary image.
open(cache=True)[source]

This method detaches from the current image (if an image is attached to the tool), and reattaches it (opens) to the new image.

The input image file may be in native CASA, FITS, or Miriad format. In the case of CASA images, images with float, complex float, double, and complex double valued pixels are supported. Note that only FITS images with float valued pixels are supported.

The cache parameter applies only to component list images and indicates if pixel values should be cached after they are computed for faster retrieval. It is not used for other image types.

Parameters

  • infile (path='') - image file name

  • cache (bool=True) - Cache pixel values for use while image is in use? Ignored if image being opened is not a ComponentListImage.

Returns

bool

Examples

#
print "\t----\t open Ex 1 \t----"
ia.maketestimage('anotherimage',overwrite=True) #first make 2nd image
ia.close()
ia.maketestimage('myimage',overwrite=True)      #open image myimage
ia.open('anotherimage')               # attach tool to 'anotherimage'
ia.close()
#


The {\stff open} function first closes the old \imagefile.
pad(outfile='', npixels=1, value=0, padmask=False, overwrite=False, region='', box='', chans='', stokes='', mask='', stretch=False, wantreturn=True)[source]

This method pads the directional plane of an image with a specified number of pixels on each side. The numerical and mask values of the padding pixels may also be specified. If a region is selected, a subimage of that region is created and then padded with the specified pixel parameters. Thus, padding an image of shape (ra, dec, freq) = (512, 512, 10) specifying npixels = 3 results in an image of size (518, 518, 10), with the blc of the directional plane of the original pixel set corresponding to the directional pixel of (3, 3) in the output. If wantreturn is True, an image analysis tool attached to the output image is returned. If False, none is returned.

Parameters

  • outfile (string='') - Output image name. If not specified, no persistent image is created.

  • npixels (int=1) - Number of pixels with which to pad each side of the direction plane.

  • value (double=0) - Value given to the padding pixels.

  • padmask (bool=False) - Value of the mask for the padding pixels. True$=>$good (unmasked), False$=>$bad (masked).

  • overwrite (bool=False) - Overwrite the output if it exists? Default False

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • box (string='') - Rectangular region to select in direction plane. Default is to use the entire direction plane.

  • chans (string='') - Channels to use. Default is to use all channels.

  • stokes (string='') - Stokes planes to use. Default is to use all stokes planes.

  • mask (string='') - Mask to use. Default is none.

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

  • wantreturn (bool=True) - Return an image analysis tool attached to the created subimage?

Returns

image

Examples

ia.fromshape("",[50, 50, 10])
# pad it with 5 pixels of value 2.5 all unmasked
padded = ia.pad(npixels=5, value=2.5, padmask=True)
ia.done()
# returns [60, 60, 10]
paddedshape = padded.shape()
padded.done()
pbcor(pbimage='', outfile='', overwrite=False, box='', region='', chans='', stokes='', mask='', mode='divide', cutoff=-1.0, stretch=False)[source]

Correct an image for primary beam attenuation using an image of the primary beam pattern. The primary beam pattern can be provided as an image, in which case 1. it must have the same shape as the input image and its coordinate system must be the same, or 2. it must be a 2-D image in which case its coordinate system must consist of a (2-D) direction coordinate which is the same as the direction coordinate in the input image and its direction plane must be the same shape as that of the input image. Alternatively, pbimage can be an array of pixel values in which case the same dimensionality and shape constraints apply. An image tool referencing the corrected image is returned. The corrected image will also be written to disk if outfile is not empty (and overwrite=True if outfile already exists). One can choose between dividing the image by the primary beam pattern (mode=”divide”) or multiplying the image by the primary beam pattern (mode=”multiply”). One can also choose to specify a cutoff limit for the primary beam pattern. For mode=”divide”, for all pixels below this cutoff in the primary beam pattern, the output image will be masked. In the case of mode=”multiply”, all pixels in the output will be masked corresponding to pixels with values greater than the cutoff in the primary beam pattern. A negative value for cutoff means that no cutoff will be applied, which is the default.

Parameters

  • pbimage (variant='') - Name of the primary beam image which must exist or array of values for the pb response. Default “”

  • outfile (string='') - Output image name. If empty, no image is written. Default “”

  • overwrite (bool=False) - Overwrite the output if it exists? Default False

  • box (string='') - Rectangular region(s) to select in direction plane. Default is to use the entire direction plane.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • chans (string='') - Channels to use. Default is to use all channels.

  • stokes (string='') - Stokes planes to use. Default is to use all stokes planes.

  • mask (variant='') - Mask to use. Default is none.

  • mode (string='divide') - Divide or multiply the image by the primary beam image. Minimal match supported. Default “divide”

  • cutoff (float=-1.0) - PB cutoff. If mode is “d”, all values less than this will be masked. If “m”, all values greater will be masked. Less than 0, no cutoff. Default no cutoff

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

Returns

image

Examples

ia.open("attenuated.im")
pbia = ia.pbcor(pbimage="mypb.im", outname="pbcorred.im", mode="divide", cutoff=0.1)
ia.done()
# do stuff with ia tool attached to pb image and close it
pbia.done()
pixeltype()[source]

This application returns the data type of the pixels of the attached image as a string. The possible values are: “float” which indicates real valued, floating point, 32 bit pixel values, “complex” which indicates complex valued, floating point, 32 bit (for each of the real and imaginary parts) pixel values, “double” which indicates real valued, floating point, 64 bit pixel values, and “dcomplex” which indicates complex valued, floating point, 64 bit (for each of the real and imaginary parts) pixel values.

pixelvalue(pixel=[-1])[source]

This function gets the value of the image and the mask at the specified pixel coordinate. The values are returned in a record with fields ’value’, ’mask’ and ’pixel’. The value is returned as a quantity, the mask as a Bool (T is good). The ’pixel’ field holds the actual pixel coordinate used.

If the specified pixel coordinate is off the image, “{}” is returned.

Excessive elements in pixel are silently discarded. Missing elements are given the (nearest integer) value of the reference pixel. This is reflected in the output record ’pixel’ field.

Parameters

  • pixel (intVec=[-1]) - Pixel coordinate

Returns

record

Examples

#
print "\t----\t pixelvalue Ex 1 \t----"
ia.maketestimage();
ia.pixelvalue()
#{'mask': True,
# 'pixel': array([55, 37]),
# 'value': {'unit': 'Jy/beam', 'value': 2.5064315795898438}}
print ia.pixelvalue([-1,-1])
# {}
print ia.pixelvalue([9])
#{'mask': True,
# 'pixel': array([ 9, 37]),
# 'value': {'unit': 'Jy/beam', 'value': 0.14012207090854645}}
print ia.pixelvalue([9,9,9])
#{'mask': True,
# 'pixel': array([9, 9]),
# 'value': {'unit': 'Jy/beam', 'value': -0.45252728462219238}}
ia.close()
#
putchunk(pixels='', blc=[-1], inc=[1], list=False, locking=True, replicate=False)[source]

This method puts an array into the image (which must be writable, eg it will fail on FITS images). If there is a default pixel mask, it is ignored. It is the complement of the getchunk() method. The blc, trc, and increment (inc) may be specified. If they are unspecified, they default to the beginning of the image and an increment of one.

Any illegal blc values are set to zero. Any illegal inc values are set to unity.

An error will result if an attempt is made to put an array the extends beyond the image edge (i.e., it is not truncated or decimated).

If there are fewer axes in the array than in the image, the array is assumed to have trailing axes of length unity. Thus, if you have a 2D array and want to put it in as the YZ plane rather than the XY plane, you must ensure that the shape of the array is [1,nx,ny].

However, the replicate parameter can be used to replicate the array throughout the image (from the blc to the trc). For example, if a 2D array is provided for a 3D image, it can be replicated along the third axis by setting replicate=True. The replication is done from the specified blc to the end of the image. Method putregion() can be used to terminate the replication at a trc value.

The locking parameter controls two aspects. If True, then after the method is called, the image is unlocked (so some other process can acquire a lock) and it is indicated that the image has changed. The reason for having this argument is that the unlocking and updating processes are quite expensive. If putchunk() is called repeatedly in eg, a loop, it is advisable to set this parameter to True.

A related function is putregion(), which supports putting the pixel and mask values into a more complex region. Method putchunk() is faster and therefore preferable for repeated operation in loops if the pixel mask is not required.

See also the methods set() and calc() which can also be used to change pixel values.

Parameters

  • pixels (variant='') - Numeric array. Required input.

  • blc (intVec=[-1]) - Bottom-Left-Corner (start) of location in image. Default is start of image.

  • inc (intVec=[1]) - increment (stride) along axes

  • list (bool=False) - List bounding box to logger?

  • locking (bool=True) - Unlock image after use?

  • replicate (bool=False) - Replicate array through image

Returns

bool

Examples

We can clip all pixels to be {\tt <= } 5 as follows.


#
print "\t----\t putchunk Ex 1 \t----"
ia.fromshape(shape=[10,10])   # create an example image
pix = ia.getchunk()           # get pixels to modify from example image
for i in range(len(pix)):
  pix[i] = list(pix[i])       # convert tuple to list so it can be modified
  for j in range(len(pix[i])):
    pix[i][j] = i\*10 + j
  pix[i] = tuple(pix[i])      # convert list back to tuple
ia.putchunk(pix)              # put pixels back into example image
print pix                     # pixels have values 0-99
pix2 = ia.getchunk()          # get all pixels into an array (again)
for i in range(len(pix2)):
  pix2[i] = list(pix2[i])     # convert tuple to list so it can be modified
  for j in range(len(pix2[i])):
    if pix2[i][j] \> 5:
      pix2[i][j] = 5          # clip values to 5
  pix2[i] = tuple(pix2[i])    # convert list back to tuple
ia.putchunk(pix2)             # put array back into image
print ia.getchunk()
ia.close()
#



The above example shows how you could clip an image to a value.  If
all the pixels didn't easily fit in memory, you would iterate through
the image chunk by chunk to avoid exhausting virtual memory.  Better
would be to do this via LEL through function calc.

Suppose we wanted to set the fifth XY plane to 1.

We could do so as follows:


#
print "\t----\t putchunk Ex 2 \t----"
ia.fromshape(shape=[10,10,10])
imshape = ia.shape()
pix = ia.makearray(1, [imshape[0],imshape[1]])
ia.putchunk(pix, blc=[0,0,4])
print ia.getchunk()[0:3]
ia.close()
#



Suppose we wanted to set the first YZ plane to 2.



#
print "\t----\t putchunk Ex 3 \t----"
ia.fromshape(shape=[10,10,10])
imshape = ia.shape()
pix = ia.makearray(2, [1,imshape[1],imshape[2]])
ia.putchunk(pix)
print ia.getchunk()[0:3]
ia.close()
#
putregion(pixels='', pixelmask='', region='', list=False, usemask=True, locking=True, replicate=False)[source]

This function replaces data and/or pixel mask values in the image in the specified region. The pixels and/or pixelmask arrays must be the shape of the bounding box, and the whole bounding box is replaced in the image. The region is only used to specify the bounding box. If the region extends beyond the image, it is truncated. If the pixels or pixelmask array shapes do not match the bounding box, an error will result. Values in the pixels array must share the same domain as the pixel values in the image. If the pixels array contains real values and the image pixels contain complex values (or vice versa), an exception will be thrown.

When you put a pixel mask, it either replaces the current default pixel mask, or is created.

The usemask parameter is only relevant when you are putting pixel values and there is a pixel mask (meaning also the one you might have just put in place). If usemask=True, then only pixels for which the mask is good (True) are altered. If usemask=False, then all the pixels in the region are altered.

The replicate parameter can be used to replicate the array throughout the image (from the blc to the trc). For example, if a two dimensional array is provided for a three dimensional image, it can be replicated along the third axis by setting replicate=True. The replication is done in the specified region.

The locking parameter controls two things. If True, then after the method is called, the image is unlocked (so some other process can acquire a lock) and it is indicated that the image has changed. The reason for this parameter is that the unlocking and updating processes are quite expensive. If putregion() is being called multiple times, in a for loop, for example, it is recommended to set locking=True (and to consider using putchunk() instead).

See the related functions putchunk, set and calc.

Parameters

  • pixels (variant='') - The pixel values. Default is none.

  • pixelmask (variant='') - The pixel mask values. Default is none.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • list (bool=False) - List the bounding box and any mask creation to the logger

  • usemask (bool=True) - Honour the mask when putting pixels

  • locking (bool=True) - Unlock image after use?

  • replicate (bool=False) - Replicate array through image

Returns

bool

Examples

Suppose that we have a 2-dimensional image.  First we recover the pixel
and \pixelmask\ values from a polygonal region.  Then we change the values in
the array that are within the region to zero and replace the data. 



#
print "\t----\t putregion Ex 1 \t----"
ia.maketestimage()                         # Attach an image to image tool
x = ['3pix','6pix','9pix','6pix','5pix','5pix','3pix'] # X vector abs pixels
y = ['3pix','4pix','7pix','9pix','7pix','5pix','3pix'] # Y vector abs pixels
mycs = ia.coordsys()
r1 = rg.wpolygon(x,y,csys=mycs.torecord()) # Create polygonal world region
mycs.done()
pixels = ia.getregion(r1)                  # Recover pixels
pixelmask = ia.getregion(r1, getmask=True)    # and mask
for i in range(len(pixels)):
  pixels[i] = list(pixels[i])              # convert tuple to list for mods
  for j in range(len(pixels[i])):
    if pixelmask[i][j]:
      pixels[i][j] = 0                     # Set pixels where mask is T to zero
  pixels[i] = tuple(pixels[i])             # convert list back to tuple
ia.putregion(pixels=pixels, pixelmask=pixelmask,
             region=r1)                    # Replace pixels only
ia.close()
#
pv(outfile='', start='', end='', center='', length='', pa='', width='1', unit='arcsec', overwrite=False, region='', chans='', stokes='', mask='', stretch=False, wantreturn=True)[source]

Create a position-velocity image by specifying either two points between which a slice is taken in the direction coordinate or a center, position angle, and length describing the slice. The spectral extent of the resulting image will be that provided by the region specification or the entire spectral range of the input image if no region is specified. One may not specify a region in direction space; that is accomplished by specifying the start and end points or the center, length, and position angle of the slice. The parameters start and end may be specified as two element arrays of numerical values, in which case these values will be interpreted as pixel locations in the input image. Alternatively, they may be expressed as arrays of two strings each representing the direction. These strings can either represent quantities (eg [“40.5deg”, “0.5rad”) or be sexigesimal format (eg [“14h20m20.5s”,”-30d45m25.4s”], [“14:20:20.5s”,”-30.45.25.4”]). In addition, they may be expressed as a single string containing the longitude and latitude values and optionally a reference frame value, eg “J2000 14:20:20.5s -30.45.25.4”. The center parameter is specified in the same way. The length parameter may be specified as a single numerical value, in which case it is interpreted as the length in pixels, or a valid quantity, in which case it must have units conformant with the direction axes units. The pa (position angle) parameter must be specified as a valid quantity with angular units. The position angle is interpreted in the usual astronomical sense; ie measured from north through east. Either start/end or center/pa/length must be specified; if a parameter from one of these sets is specified, a parameter from the other set may not be specified. In either case, the end points of the segment must fail within the input image, and they both must be at least 2 pixels from the edge of the input image to facilite rotation (see below).

One may specify a width, which is the number of pixels centered along and perpendicular to the direction slice that are used for averaging along the slice. The width may be specified as an integer, in which case it must be positive and odd. Alternatively, it may be specified as a valid quantity string (eg, “4arcsec”) or quantity record (eg qa.quantity(“4arcsec”). In this case, units must be conformant to the direction axes units (usually angular units) and the specified quantity will be rounded up, if necessary, to the next highest equivalent odd integer number of pixels. The default value of 1 represents no averaging. A value of 3 means average one pixel on each side of the slice and the pixel on the slice. Note that this width is applied to pixels in the image after it has been rotated (see below for a description of the algorithm used). The end points of the specified segment must fail within the input image, and they both must be at least 2 pixels from the edge of the input image to facilite rotation (see below).

One may specify the unit for the angular offset axis.

A True value for the wantreturn parameter indicates that an image analysis tool attached to the created image should be returned. Nothing is returned if wantreturn is False, but then outfile should be specified (unless perhaps you are debugging).

Internally, the image is first rotated, padding first if necessary to include relevant pixels that would otherwise be excluded by the rotation operation, so that the slice is horizontal, with the starting pixel left of the ending pixel. Then, the pixels within the specified width of the slice are averaged and the resulting image is written and/or returned. The output image has a linear coordinate in place of the direction coordinate of the input image, and the corresponding axis represents angular offset with the center pixel having a value of 0.

The equivalent coordinate system, with a (usually) rotated direction coordinate (eg, RA and Dec) is written to the output image as a table record. It can be retrieved using the table tool as shown in the example below.

Parameters

  • outfile (string='') - Output image name. If empty, no image is written. Default “”

  • start (variant='') - The starting point in the direction plane (array of two values). If specified, end must also be specified and none of center, pa, nor length may be specified.

  • end (variant='') - The ending point in the direction plane (array of two values). If specified, start must also be specified and none of center, pa, nor length may be specified.

  • center (variant='') - The center point in the direction plane (array of two values). If specified, length and pa must also be specified and neither of start nor end may be specified.

  • length (variant='') - The length of the segment in the direction plane. If specified, center and pa must also be specified and neither of start nor end may be specified.

  • pa (variant='') - The position angle of the segment in the direction plane, measured from north through east. If specified, center and length must also be specified and neither of start nor end may be specified.

  • width (variant='1') - Width in pixels for averaging pixels perpendicular to the slice. Must be an odd integer >= 1 (1 means only use the pixels along the slice), or a quantity which will be rounded up if necessary so it corresponds to the next highest odd number of pixels.

  • unit (string='arcsec') - Unit for the offset axis in the resulting image. Must be a unit of angular measure.

  • overwrite (bool=False) - Overwrite the output if it exists?

  • region ({record, string}='') - Region selection. Default is to use the full image. No selection is permitted in the direction plane.

  • chans (string='') - Channels to use. Channels must be contiguous. Default is to use all channels.

  • stokes (string='') - Stokes planes to use. Planes must be contiguous. Default is to use all stokes planes.

  • mask (variant='') - Mask to use. Default is none.

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

  • wantreturn (bool=True) - Return an image analysis tool attached to the created image?

Returns

image

Examples

ia.open("my_spectral_cube.im")
# create a pv image with the position axis running from ra, dec pixel positions of [45, 50] to [100, 120]
# in the input image
mypv = ia.pv(outfile="pv.im", start=[45,50], end=[100,120], wantreturn=True)
ia.done()
# analyze the pv image, such as get statistics
pvstats = mypv.statistics()
# when done, close the tool to release system resources
mypv.done()

# get the alternate coordinate system information
tb.open("pv.im")
alternate_csys_record = tb.getkeyword("misc")["secondary_coordinates"]
tb.done()
rebin(outfile='', bin='', region='', mask='', dropdeg=False, overwrite=False, stretch=False, crop=False)[source]

This application rebins the current image by the specified integer binning factors for each axis. It supports both float valued and complex valued images. The corresponding output pixel value is the average of the input pixel values. The output pixel will be masked bad if there were no good input pixels. A polarization axis cannot be rebinned.

The binning factors array must contain at least one element and no more elements than the number of input image axes. If the number of elements specified is less than the number of image axes, then the remaining axes not specified are not rebinned. All specified values must be positive. A value of one indicates that no rebinning of the associated axis will occur.

Binning starts from the origin pixel of the bounding box of the selected region or the origin pixel of the input image if no region is specified. The value of crop is used to determine how to handle cases where there are pixels at the end of the axis that do not form a complete bin. If crop=True, extra pixels at the end of the axis are discarded. If crop=False, the remaining pixels are averaged into the final bin along that axis. Should the length of the axis to be rebinned be an integral multiple of the associated binning factor, the value of crop is irrelevant.

A value of dropdeg=True will result in the output image not containing axes that are degenerate in the specified region or in the input image if no region is specified. Note that, however, the binning factors array must still account for degenerate axes, and the binning factor associated with a degenerate axis must always be 1.

If outfile is given, the image is written to the specified disk file. If outfile is unset, the Image  is associated with a temporary image. This temporary image may be in memory or on disk, depending on its size. When you destroy the on-the-fly Image  returned by this function (with the done function) this temporary image is deleted.

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • bin (intVec='') - Binning factors for each axis

  • region ({string, record}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • dropdeg (bool=False) - Drop degenerate axes

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

  • crop (bool=False) - Remove pixels from the end of an axis to be rebinned if there are not enough to form an integral bin?

Returns

image

Examples

#
print "\t----\t rebin Ex 1 \t----"
ia.maketestimage();
im2 = ia.rebin(bin=[2,3]);
im2.done()
ia.close()
#
regrid(outfile='', shape=[-1], csys='', axes=[-1], region='', mask='', method='linear', decimate=10, replicate=False, doref=True, dropdeg=False, overwrite=False, force=False, asvelocity=False, stretch=False)[source]

This function regrids the current image onto a grid specified by the given coordinate system. The shape of the output image may also be specified.

The coordinate system must be specified via a cs tool (using cs.torecord()). It is optional; if not specified, the coordinate system from the input image (ie, the one to which you are applying the regrid function) is used. The order of the coordinates and axes in the output image is always the same as the input image. It simply finds the relevant coordinate in the supplied coordinate system in order to determine the regridding parameters. The supplied coordinate system must have at least as many coordinates as are required to accomodate the axes that are to be regridded (eg, if the first two axes are to be regridded, and these belong to a direction coordinate, one direction coordinate in the supplied coordinate system is required). Coordinates pertaining to axes that are not being regridded are supplied from the input image, not the specified coordinate system.

Reference changes are handled (eg, J2000 to B1950, LSR to TOPO). In general, the conversion machinery attempts to work out how sophisticated it needs to be (eg, is the regridding being done from LSR to LSR or from LSR to TOPO). However, it errs on the side of caution if the conversion machine requires more information than it actually needs. For full frame conversions, one needs to know things like the position on the earth’s surface where the observation occurred, direction (celestial coordinates) of observation, and time of observation.

If you get such errors and you are doing a frame conversion, then that means you must insert some extra information into the coordinate system of your image. Most likely it’s the time (in which case you can use cs.setepoch()) and/or position (in which case you can use cs.settelescope()) that are missing. If you get these errors and you are certain that you are not specifying a frame change (eg, regrid LSR to LSR) then try setting doref=False. This will (silently) bypass all possible frame conversions. Note that if you are requesting a frame conversion and you set doref=False, no warnings will be emitted and the output image will likely be nonsensical.

If you regrid a plane holding a direction coordinate and the units are Jy/pixel, then the output is scaled to conserve flux (roughly; just one scale factor at the reference pixel is computed).

Regridding of complex-valued images is supported. The real and imaginary parts are regridded independently and the resulting regridded pixel values are combined to form the regridded, complex-valued image.

A variety of interpolation schemes are provided (you need only specify the first three characters to the method parameter). The cubic interpolation is substantially slower than linear, and often the improvement is modest. By default, linear interpolation is used.

You specify the shape of the output image using the shape parameter and which output axes you want to regrid. Note that a Stokes axis cannot be regridded (you will get a warning if you try).

The axes parameter cannot be used to discard axes from the output image; it can only be used to specify which output axes are going to be regridded and which are not. Any axis that you are not regridding must have the same output shape as the input image shape for that axis.

The axes parameter can also be used to specify the order in which the output axes are regridded. This may give you significant performance benefits. For example, imagine we are going to regrid a spectral-line cube of shape [512,512,1204] to shape [256,256,32]. If you specified axes=[0,1,2] then first, the direction axes would be regridded for each of the 1024 pixels (and stored in a temporary image). Then each spectral profile at each spatial location in the temporary image would be regridded to 32 pixels. You could speed this process up significantly by setting axes=[2,0,1]. In this case, first each spectral profile would be regridded to 32 pixels, and then each plane of the 32 pixels would be regridded. Note that the order of axes does not affect the order of the shape parameter; ie, it should be given in the natural pixel axis order of the image ()[256,256,32] in both cases of this example).

You can also specify a region to be applied to the input image. If you do this, you need to be careful with the output shape for non-regridded axes (must match that of the region - use function ia.boundingbox() to determine that).

If the outfile parameter is specified, the image is written to the specified disk file. If this parameter is unset, the on-the-fly image analysis tool returned by this method is associated with a temporary image. This temporary image may be in memory or on disk, depending on its size. When you destroy the on-the-fly image analysis tool (with either the ia.close() or ia.done() methods), this temporary image is deleted.

The replicate parameter can be used to simply replicate pixels rather than regridding them. Normally replicate=False, for every output pixel, its world coordinate is computed and the corresponding input pixel found (then a little interpolation grid is generated). If replicate=True, then for every output axis, a vector of regularly sampled input pixels is generated (based on the ratio of the output and input axis shapes). So this just means the pixels get replicated (by whatever interpolation scheme you use) rather than regridded in world coordinate space. This process is much faster, but its not a True world coordinate based regrid.

As decribed above, when replicate=False, a coordinate is computed for each output pixel; this is an expensive operation. The decimate parameter allows you to decimate the computation of that coordinate grid to a sparse grid, which is then filled in via fast interpolation. The default is decimate=10. The number of pixels per axis in the sparse grid is the number of output pixels for that axis divided by the decimation factor. A factor of 10 does pretty well. You may find that for very non-linear coordinate systems (e.g. very close to the pole) that you have to reduce the decimation factor. You may also have to reduce the decimation factor if the number of pixels in the output image along an axis to be regridded is less than about 50, or the output image may be completely masked.

If one of the axes to be regridded is a spectral axis and asvelocity=True, the axis will be regridded to match the velocity, not the frequency, description of the template spectral coordinate. Thus the output pixel values will correspond only to the velocity, not the frequency, of the output axis.

Sometimes it is useful to drop axes of length one (degenerate axes). Setting the dropdeg parameter to True will do that. It will discard the axes from the input image. Therefore the output shape and coordinate system that you supply must be consistent with the input image after the degenerate axes are dropped.

The force parameter can be used to force all specified axes to be regridded, even if the algorithm determines that they don’t need to be (because the input and output coordinate information is identical).

The cs tool has a useful method, cs.setreferencelocation(), that can be used to keep a specific world coordinate in the center of an image when regridding (see example below).

The output pixel mask will be True (good) unless the regridding failed to find a value for that output pixel in which case it will be False. For example, if the total input mask (default input pixel mask plus OTF mask) for all of the relevant input pixels were masked bad then the output pixel would be masked (False).

MULTIPLE AXIS COORDINATES LIMITATION Some cooordinates contain multiple axes. For example, a direction coordinate contains both longitude-like and latitude-like axes. A linear coordinate can also contain multiple axes. When you regrid *any* axis from a coordinate which contains multiple axes, you must fully specify the coordinate information for all axes in that coordinate in the coordinate system that you provide. For example, if you have a linear coordinate with two axes and you want to regrid axis one only. In the coordinate system you provide, the coordinate information for axis two (not being regridded) must correctly be a copy from the input coordinate system (it won’t be filled in for you).

If an image has per-plane beams and one attempts to regrid the spectral axis, an exception is thrown.

IMPORTANT NOTE ABOUT FLUX CONSERVATION in general regridding is inaccurate for images that the angular resolution is poorly sampled. A check is done for such cases and a warning message is emitted if a beam present. However, no such check is done if there is no beam present. To add a restoring beam to an image, use ia.setrestoringbeam().

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • shape (intVec=[-1]) - Shape of output image. Default is input shape.

  • csys (record='') - Coordinate System for output image. Default is input image coordinate system.

  • axes (intVec=[-1]) - The output pixel axes to regrid. Default is all.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • method (string='linear') - The interpolation method. String from ‘nearest’, ‘linear’, ‘cubic’.

  • decimate (int=10) - Decimation factor for coordinate grid computation

  • replicate (bool=False) - Replicate image rather than regrid?

  • doref (bool=True) - Turn on reference frame changes

  • dropdeg (bool=False) - Drop degenerate axes

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • force (bool=False) - Force specified axes to be regridded

  • asvelocity (bool=False) - Regrid spectral axis in velocity space rather than frequency space?

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

Returns

image

Examples

#
print "\t----\t regrid Ex 1 \t----"
ia.maketestimage('radio.image', overwrite=True)  
ia.maketestimage('optical.image', overwrite=True)
mycs = ia.coordsys();     # get optical image co-ordinate system
ia.open('radio.image')
imrr = ia.regrid(outfile='radio.regridded', csys=mycs.torecord(),
                  shape=ia.shape(), overwrite=True)
#viewer()
mycs.done()
imrr.done()
ia.close()
#



In this example, we regrid a radio image onto the grid of an optical
image - this probably (if the optical FITS image was correctly labelled
!!) will involve a projection change (optical images are usually TAN
projection, radio usually SIN).  
remove(done=False, verbose=True)[source]

This function first closes the  which detaches it from its underlying . It then deletes that . If done=False, the  is still viable, and can be used with function open to open a new . Otherwise the  is destroyed. If verbose=True, the logger will receive a progress report.

Parameters

  • done (bool=False) - Destroy this tool after deletion

  • verbose (bool=True) - Send a progress report to the logger.

Returns

bool

Examples

#
print "\t----\t remove Ex 1 \t----"
ia.maketestimage('myimage',overwrite=True)
ia.close()
ia.maketestimage('myotherimage',overwrite=True)
ia.close()
ia.open('myimage')          # Attach to 'myimage'
ia.remove(F)                # Close imagetool and delete 'myimage'
ia.open('myotherimage')     # Open new imagefile 'myotherimage'
ia.remove()
print "!!!EXPECT THE FOLLOWING TO GENERATE AN ERROR MESSAGE!!!"
ia.open('myimage')          # 'myimage' was deleted above
ia.close()
#
removefile(file='')[source]

This function deletes the specified image file.

Parameters

  • file (string='') - Name of image file/directory to be removed. Must be specified.

Returns

bool

Examples

#
print "\t----\t removefile Ex 1 \t----"
ia.maketestimage('myimage',overwrite=True)
ia.close()
ia.removefile('myimage')                    # remove image 'myimage'
ia.maketestimage('myimage',overwrite=False) # error here if 'myimage' exists
ia.close()
ia.removefile('myimage')
#
rename(name='', overwrite=False)[source]

This function renames the  associated with the . If a file with name name already exists, you can overwrite it with the argument overwrite; otherwise a fail will result.

Parameters

  • name (string='') - The new image file name

  • overwrite (bool=False) - Overwrite target file if it already exists

Returns

bool

Examples

Example:
#
print ' ---- rename Ex 1, rename a new image ----'
ia.maketestimage('myimage',overwrite=True)
print ia.name(strippath=True)
#myimage
ia.rename('newimage', overwrite=True)
print ia.name(strippath=True)
#newimage
#
print ' ---- rename Ex 2, rename an existing image ----'
ia.open('originalimage',overwrite=True)
ia.rename('newimage', overwrite=True)
ia.close()
#
replacemaskedpixels(pixels='', region='', mask='', update=False, list=False, stretch=False)[source]

This application replaces the values of all pixels whose total input mask (default input  and OTF mask) is bad (F) with the specified value. It supports both float valued and compplex valued images.

If the argument update is F (the default), the actual  is left unchanged. That is, masked pixels remain masked. However, if you set update=T then the  will be updated so that the  will now be T (good) where the total input mask was F (bad).

See maskhandler for information on how to set the default .

There are a few ways in which you can specify what to replace the masked pixel values by.

  • First, you can give the pixels argument a simple numeric scalar (e.g. pixels=1.0). Then, all masked values will be replaced by that value.

  • Second, you can give a scalar expression string (e.g. pixels=’min(myimage)’). Then, all masked values will be replaced by the scalar that results from the expression. If the scalar expression is illegal (e.g. in the expression pixels=’min(myimage)’ there were no good pixels in myimage) then the value 0 is used for replacement.

  • Third, you can give a expression string which has the same shape as the  you are applying the function to. For example, putting pixels=’myotherimage’ means replace all masked pixels in this  with the equivalent pixel in the  called myotherimage.

    Your expression might be quite complex, and you can think of it as producing another masked lattice. However, in the replace process, the mask of that expression lattice is ignored. Thus, only the mask of the  you are replacing and the pixel values of the expression lattice are relevant.

    The expression must conform with the subimage formed by applying the  to the image (i.e. that associated with this Image ). If you use the mask argument as well, the  is applied to it as well (see examples).

Parameters

  • pixels (variant='') - The new value(s), Numeric scalar or LEL expression

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • update (bool=False) - Update mask as well?

  • list (bool=False) - List the bounding box to the logger

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

Returns

bool

Examples

#
print "\t----\t replacemaskedpixels Ex 1 \t----"
ia.maketestimage('zz1',overwrite=True)
ia.calcmask('zz1<0')
ia.replacemaskedpixels(0.0)
ia.replacemaskedpixels('min(zz1)')
ia.close()
#



These examples replace all masked pixels by the specified scalar.  In
the second case, the scalar comes from a LEL expression operating on
{\sff zz1} (or it could be from an LEL expression operating on some
other image).
restoringbeam(channel=-1, polarization=-1, mbret='list')[source]

This function returns the restoring beam(s), if any, of the attached image.

If mbret=”matrix” (or “m”), an exception will be thrown if the attached image does not have per-plane beams, or if channel or polarization is specified to be non-negative. See below for more information on the mbret parameter.

If mbret=”list” (or “l”) and the attached image has no restoring beam(s), an empty dictionary is returned.

If the image has a single traditional restoring beam and mbret=”list”, the beam is returned as a dictionary no matter what channel and polarization values are specified. This dictionary has keys “major”, “minor”, andi “positionangle”, and each of these fields contains a dictionary with “value” and “unit” keys which provide the quantity associated with each beam parameter.

If mbret=”list”, and if the image has per-plane beams, and the image has both spectral and polarization axes, and if both channel and polarization are set to non-negative values that are less than the number of planes along each parameter’s representative axis, the beam for that particular channel/polarization pair is returned. If at least one of these values is set to greater or equal to the number of planes along that parameter’s representative axis, and that axis exists and is not degenerate, then an exception will be thrown. In the case of a non-extant axis or a degenerate axis, the parameter associated with that axis can be set to any value and a beam will be returned, because in that case there is exactly one beam for each plane along the other parameter’s representative axis, and so there is no ambiguity regarding which beam to return. In the case where both spectral and polarization axes exist and are not degenerate, an exception will be thrown if one of channel or polarization is set to a non-negative value and the other is set to a negative value. In the above cases in which a beam is returned, the returned dictionary will have the same structure as that previously described in the single beam case.

If the image contains multiple beams and both channel and polarization are specified to be negative, the structure of the returned dictionary depends on the specified value of mbret. Supported values of this parameter are “list” and “matrix” (case insensitive, mimimum match supported). In both cases, the returned dictionary will contain the keys “nChannels”, which contains an integer value equal to the number of channels, and “nStokes”, which contains an integer value equal to the number of polarizations. In the case where one axis doesn’t exist, the associated value will be 1.

In the case of mbret=”list”, the returned dictionary will contain the key “beams”, which contains a sub-dictionary of information for all beams. This subdictionary contains keys “*0” through “*(c-1)”, where c is the value associated with the “nChannels” key. Each of these keys has an associated value which is a subdictionary containing keys “*0” through “*(p-1)”, where p is the value associated associated with the “nStokes” key. Each of these keys has an assciated value of a beam dictionary, the structure of which is the same as described previously for a single beam image, which is associated with that particular channel/polarization pair.

In the case of mbret=”matrix”, the returned dictionary will have keys “major”, “minor”, and “pa” which hold values representing the major axes, the minor axes, and the position angle, respectively, of all beams. Each of these fields contains a dictionary which represents a quantity with keys “unit” and “value”. The “unit” field contains a string representing the unit of the associated value. The “value” field contains a matrix with values corresponding to the associated parameter for the associated beam. In the case where an image does not have both a spectral and polarization axis, these matrices will have shape (n, 1), where n is the number of planes on the extant axis. In the case where an image has both spectral and polarization axes, the matrices will have shape (m, n) where m is the number of channels if the spectral axis precedes the polarization axis, or the number of polarization planes otherwise, and n is the number of planes on the axis not represented by m. So, an image with a spectral axis that has 5 channels and precedes the polarization axis that has 4 planes, the returned matrix shape will be (5, 4). In the case where the number of spectral and polarization planes are the same as in the previous example, but the polarization axis precedes the spectral axis, the matrix shape will be (4, 5). Thus, the matrix shape pair follows the order of the image spectral and polarization axes. This mode is useful for returning to the python UI the major axes, minor axes, and position angles directly as numpy arrays.

The restoring beam(s) of an image may be set with image tool method setrestoringbeam.

Parameters

  • channel (int=-1) - The zero-based spectral channel number for a per-plane beam. Default -1

  • polarization (int=-1) - The zero-based polarization plane number for a per-plane beam. Default -1

  • mbret (string='list') - Determines the format of the dictionary returned if there are multiple beams and both channel and polarization are negative. Supported values are “list” and “matrix” (case insensitive, minimum match supported).

Returns

record

Examples

#
print "\t----\t restoringbeam Ex 1 \t----"
ia.maketestimage()      
print ia.restoringbeam()
#{'major': {'unit': 'arcsec', 'value': 53.500004857778549},
# 'minor': {'unit': 'arcsec', 'value': 34.199998900294304},
# 'positionangle': {'unit': 'deg', 'value': 6.0}}
ia.close()
#
rotate(outfile='', shape=[-1], pa='0deg', region='', mask='', method='cubic', decimate=0, replicate=False, dropdeg=False, overwrite=False, stretch=False)[source]

This function rotates two axes of an image. These axes are either those associated with a Direction coordinate or with a Linear coordinate. The Direction coordinate takes precedence. If rotating a Linear coordinate, it must hold precisely two axes.

The method is that the Coordinate is rotated and then the input image is regridded to the rotated Coordinate System.

If the image brightness units are Jy/pixel then the output is scaled to conserve flux (roughly; just one scale factor at the reference pixel is computed).

A variety of interpolation schemes are provided (you need only specify the first three characters to method). The cubic interpolation is substantially slower than linear. By default you get cubic interpolation.

You can specify the shape of the output image (shape). However, all axis that are not regrided retain the same output shape as the input image shape for that axis. Only the direction coordinate axes are regridded.

You can also specify a  to be applied to the input image. If you do this, you need to be careful with the output shape for non-regridded axes (must match that of the region - use function boundingbox to find that out).

If outfile is given, the image is written to the specified disk file. If outfile is unset, the on-the-fly Image  returned by this function is associated with a temporary image. This temporary image may be in memory or on disk, depending on its size. When you destroy the on-the-fly Image  (with the done function) this temporary image is deleted.

The argument replicate can be used to simply replicate pixels rather than regridding them. Normally (replicate=F), for every output pixel, its world coordinate is computed and the corresponding input pixel found (then a little interpolation grid is generated). If you set replicate=T, then what happens is that for every output axis, a vector of regularly sampled input pixels is generated (based on the ratio of the output and input axis shapes). So this just means the pixels get replicated (by whatever interpolation scheme you use) rather than regridded in world coordinate space. This process is much faster, but its not a True world coordinate based regrid.

As decribed above, when replicate is False, a coordinate is computed for each output pixel; this is an expensive operation. The argument decimate allows you to decimate the computation of that coordinate grid to a sparse grid, which is then filled in via fast interpolation. The default for decimate is 0 (no decimation). The number of pixels per axis in the sparse grid is the number of output pixels for that axis divided by the decimation factor. A factor of 10 does pretty well. You may find that for very non-linear coordinate systems (e.g. very close to the pole) that you have to reduce the decimation factor.

The output  will be good (T) unless the regridding failed to find a value for that output pixel in which case it will be bad (F). For example, if the total input mask (default input  plus OTF mask) for all of the relevant input pixels were masked bad then the output pixel would be masked bad (F).

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • shape (intVec=[-1]) - Shape of output image. Default is shape of input image.

  • pa (variant='0deg') - Angle by which to rotate. Default is no rotation.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • method (string='cubic') - The interpolation method. String from ‘nearest’, ‘linear’, or ‘cubic’.

  • decimate (int=0) - Decimation factor for coordinate grid computation

  • replicate (bool=False) - Replicate image rather than regrid?

  • dropdeg (bool=False) - Drop degenerate axes

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

Returns

image

Examples

ia.maketestimage()
imr=ia.rotate(outfile="rotated.im", pa='45deg')
imr.done()
ia.close()


In this example, we rotate the direction coordinate axes (RA/Dec) of a
test image by 45 degress and regrid the image onto the axes.
rotatebeam(angle='0deg')[source]

This method rotates the attached image’s beam(s) counterclockwise through the specified angle. This is the same thing as increasing the position angle(s) of the beam(s) by the specified angle. If the image does not have a beam, no changes to the image are made. If the image has multiple beams, all the beams are rotated through the same angle.

Parameters

  • angle (variant='0deg') - Angle by which to rotate image’s beam(s). Default is no rotation.

Returns

bool

Examples

# rotate any and all beams in the image (increase their position angles) by 30 degrees.
ia.open("my.im")
ia.rotatebeam("30deg")
ia.done()
sepconvolve(outfile='', axes=[-1], types=[''], widths='', scale=-1, region='', mask='', overwrite=False, stretch=False)[source]

This function does Fourier-based convolution of an  by a specified separable kernel.

If outfile is given, the image is written to the specified disk file. If outfile is unset, the on-the-fly Image  returned by this function is associated with a temporary image. This temporary image may be in memory or on disk, depending on its size. When you destroy the Image  (with the done function) this temporary image is deleted.

You specify which axes of the image you wish to convolve, by what kernel of what width. The kernel types can be shortened to ’gauss’, ’hann’ and ’box’.

You specify the widths of the convolution kernels via the argument widths. The values can be specified as a vector of three different types.

  • Quantity - for example widths=qa.quantity(“1arcsec 0.00001rad”). Note that you can use pixel units, viz. widths=qa.quantity(“10pix 0.00001rad”) see below.

  • String - for example widths=”1km 2arcsec” (i.e. a string that qa.quantity() accepts).

  • Numeric - for example widths=[10,20]. In this case, the units of the widths are assumed to be in pixels.

The interpretation of widths depends upon the kernel type.

  • Gaussian - the specified width is the full-width at half-maximum.

  • Boxcar (tophat) - the specified width is the full width.

  • Hanning - The kernel is \(z[i] = 0.25\*y[i-1] + 0.5\*y[i] + 0.25\*y[i+1]\). The width is always 3 pixels, regardless of what you give (but you still have to give it !).

The scaling of the output image is determined by the argument scale. If you leave it unset, then autoscaling will be invoked which means that the convolution kernels will all be normalized to have unit volume to as to conserve flux.

If you do not leave scale unset, then the convolution kernel will be scaled by this value (it has peak unity before application of this scale factor).

Masked pixels will be assigned the value 0.0 before convolution. The output mask is the combination (logical OR) of the default input  (if any) and the OTF mask. Any other input  will not be copied. Use function maskhandler if you need to copy other masks too.

See also the other convolution functions convolve2d, convolve and hanning.

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • axes (intVec=[-1]) - Axes to convolve. Default is [0,1,…].

  • types (stringVec=['']) - Type of convolution kernel. Vector of strings from ‘boxcar’, ‘gaussian’, and ‘hanning’. Default is appropriately sized vector of ‘gaussian’.

  • widths (variant='') - Convolution kernel widths, Vector of numeric, quantity or string

  • scale (double=-1) - Scale factor. Default is autoscale.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

Returns

image

Examples

#
print "\t----\t sepconvolve Ex 1 \t----"
ia.maketestimage('xyv',overwrite=True)
im2 = ia.sepconvolve(outfile='xyv.con', axes=[0,1], types=["gauss","box"], widths=[10,20], overwrite=True)
im2.done()
ia.close()
#
set(pixels='', pixelmask=-1, region='', list=False)[source]

This function replaces data and/or mask values within the image in the specified . You can think of it as a simplified version of the image calculator.

Unlike the calc function, you can only set a scalar value for all pixels in the specified . For example, it can be useful to set a whole image to one value, or a mask in a  to one value. Although you could do that with the related functions putregion and putchunk, you would have to make an array of the shape of the image and if that is large, it could be resource expensive.

The value for the pixels is specified with the pixels argument. It can be given as either a Lattice Expression Language (or LEL) expression string or a simple numeric scalar. See for a detailed description of the LEL expression syntax. If you give a LEL expression it must be a scalar expression.

Note that any default mask is ignored by this function when you set pixel values. This is different from calc where the extant mask is honoured.

The value for the pixel mask is specified with the pixelmask argument (T, F, unset). If it’s unset then the mask is not changed.

If you specify pixelmask= T or F, then the mask that is affected is the current default mask (see maskhandler). If there is no mask, a mask is created for you and made the default mask.

Parameters

  • pixels (variant='') - The pixel value, LEL scalar expression or numeric scalar. Default is unset.

  • pixelmask (int=-1) - The pixel mask value. Either 0 or 1 if set. Default is unset.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • list (bool=False) - List the bounding box and any mask creation to the logger

Returns

bool

Examples

#
print "\t----\t set Ex 1 \t----"
ia.maketestimage('yy',overwrite=True)
ia.fromshape('xx', [10,20], overwrite=True)
r1 = rg.box([2,2],[6,8])         # Make a box region
ia.set(pixels=1.0)               # Set all pixels to 1
ia.set(pixels='2.0', region=r1)  # Set all pixels to 2 in the region
ia.set(pixels='min(yy)')         # Set all pixels to minimum of image yy
                                 # Set pixels in region to minimum of image xx
ia.set(pixels='min('+ia.name(strippath=True)+')', region=r1)
ia.set(pixelmask=True)              # Set mask to all T
ia.set(pixels=0, pixelmask=False, region=r1)  #Set pixels and mask in region
ia.close()
#
setbrightnessunit(unit='')[source]

This function sets the image brightness unit. Both float and complex valued images are supported. You can get the brightness unit with function brightnessunit.

Parameters

  • unit (string='') - New brightness unit

Returns

bool

Examples

#
print "\t----\t setbrightnessunit Ex 1 \t----"
ia.fromshape(shape=[10,10])
ia.setbrightnessunit('km')
print ia.brightnessunit()
#km
#
setcoordsys(csys='')[source]

This method replaces the coordinate system in the image. Coordinate systems are manipulated with the cs (coordintate system) tool. The coordinate system can be recovered from an image via the coordsys() method of the image analysis (ia) tool.

Note that changing the coordinate system using the cs tool has no effect on the original image, until it is replaced with this method; the value returned by coordsys() is a copy of, not a reference to, the image’s coordinate system.

Parameters

  • csys (record='') - Record describing new Coordinate System

Returns

bool

Examples

#
print "\t----\t setcoordsys Ex 1 \t----"
ia.fromshape(shape=[10,20])          # Make image
mycs = ia.coordsys();                # Recover Coordinate System
incr = mycs.increment('n');          # Get increment as numeric vector
incrn = incr['numeric']
for i in range(len(incrn)):
  incrn[i] = 2\*incrn[i]
incr['numeric']=incrn
mycs.setincrement(value=incr);       # Set new increment in Coordinate System
ia.setcoordsys(mycs.torecord());     # Set new Coordinate System in image
mycs.done()
ia.close()
#
sethistory(origin='', history='')[source]

A   can accumulate history information from an input  file or by you writing something into it explicitly with this function. Each element of the input vector is one line of history. The new history is appended to the old.

You can recover the history information with function history.

Parameters

  • origin (string='') - Used to set message origin. Default is image::sethistory.

  • history (stringVec='') - New history

Returns

bool

Examples

#
print "\t----\t sethistory Ex 1 \t----"
ia.maketestimage('myfile',overwrite=True)  
h = ia.history()
# Adds three lines, 'I', 'like' and 'fish'
ia.sethistory(origin="sethistory", history=["I","like","fish"])
ia.close()
#
setmiscinfo(info='')[source]

A CASA image can include user-specified or miscellaneous metadata. This metadata is stored in a data structure referred to as a miscinfo record. For example, the FITS reader ia.fromfits() puts header keywords it doesn’t otherwise use into such a record. The miscinfo record is not required to be populated, though.

This method sets the miscinfo record of an image. Note that this method will overwrite, not add to, an existing miscinfo record. Thus if the user wishes to augment an existing record, they must first capture the existing record using the ia.miscinfo() method, modify the captured record, and then replace the record in the image using setmiscinfo() by passing it the modified record. The FITS writer will attempt to write all the fields in the miscinfo record to the FITS file it creates. It can do so for scalars and 1-dimensional arrays. Records will be omitted, and multi-dimensional arrays will be flattened into one dimensional arrays.

Parameters

  • info (record='') - Miscellaneous REPLACEMENT header

Returns

bool

Examples

#
print "\t----\t setmiscinfo Ex 1 \t----"
ia.maketestimage('myfile',overwrite=True)  
info = ia.miscinfo()            # capture the miscinfo record
info['extra'] = 'a test entry'  # add our information
ia.setmiscinfo(info)            # put it back into the image
ia.close()
#
setrestoringbeam(major='', minor='', pa='Not specified', beam='', remove=False, log=True, channel=-1, polarization=-1, imagename='')[source]

This method sets the restoring beam(s) for an image or removes an existing beam(s).

An image must have exactly one of the following states:

  1. An image can have a single “traditional” or global beam. In that case, the beam applies to every channel and polarization in the image.

  2. If an image has more than one spectral channel or more than one polarization, it can have a set of beams. In this case, each channel and/or polarization will have its own beam.

  3. An image can have neither a global beam nor a beam set.

It is never permissible for an image to have both a traditional (global) beam and a set of per-plane beams. Task and method behavior is undefined in that case and any resulting products are considered corrupt.

RULES FOR BEAM MODIFICATION

  1. If remove=True, any existing beam(s) are removed.

  2. Else if imagename is specified, the beam(s) from the specified image will be copied to the image being accessed by the image tool. Multiple beams may be copied, but the two images must have the same number of frequency channels and polarization planes. If not, an exception is thrown.

  3. Else if the beam parameter is specified, it will be used. It must be fully specified. It must have exactly three items with keys “major”, “minor”, and “positionangle”. Each of these keys must be a proper quantity dictionary with keys “value” and “unit”. The units for all three items should be angular. The “major” and “minor” items must have non-negative values. If any of these conditions is not met, an exception will be thrown.

  4. Else the “major”, “minor”, and “pa” parameters must be specified. Any or all of these may be a quantity string (eg, “2arcsec”, a quantity dictionary with keys “value” and “unit”. or a numerical value. “major” and “minor” must have non-negative values. If any of these conditions is not met, an exception will be thrown. In the case of a quantity being specified the “unit” should be an angular unit or no unit. In the case of no unit or a numerical value specified, if the image already has a beam then the corresponding unit for that parameter in the current beam will be used. If the image has no beam, then “arcsec” is used for “major” and “minor”, and “deg” is used for “pa”.

If an image has no beams, a traditional (global) beam can be added by setting both channel and polarization to negative values.

If an image has no beams, a set of per-plane beams can be added by setting either or both channel and/or polarization to a non-negative value. In this case, a number of per-plane beams are added consistent with the image and they are all set to be the same with parameters equal to those specified by either the beam or major/minor/pa parameters.

If an image has a traditional beam, it can be modified by setting both channel and polarization to negative values. If one or both is not set to a negative value, an exception is thrown, and nothing is modified.

If an image has a set of per plane beams, one at a time of these can be modified by setting the appropriate channel number and/or polarization number. All the per-plane beams can be modified to the same values in one go by setting both channel and polarization to negative values. Also, in the case where an image has multiple channels, the beams associated with all channels for a given polarization can be modified to the same beam by setting polarization equal to the desired polarization plane number and by setting channel to a negative value. Similarly, in the case where an image has multiple polarizations, the beams associated with all polarizations for a given spectral channel can be modified to the same beam by setting channel equal to the desired spectral channel number and by setting polarization to a negative value.

Parameters

  • major (variant='') - Major axis FWHM, Quantity or float (e.g., 1arcsec). Default is unset.

  • minor (variant='') - Minor axis FWHM, Quantity or float (e.g., 1arcsec). Default is unset.

  • pa (variant='Not specified') - Position angle, Quantity or float (e.g., ‘5deg’). Default is unset.

  • beam (record='') - The complete restoring beam (output of restoringbeam()). Default is unset.

  • remove (bool=False) - Delete the restoring beam?

  • log (bool=True) - Write new beam values to the logger?

  • channel (int=-1) - Zero-based channel number for which to set a per plane beam. If the image has a traditional beam, set to less than zero. Default -1.

  • polarization (int=-1) - Zero-based polarization number for which to set a per plane beam. If the image has a traditional beam, set to less than zero. Default -1.

  • imagename (string='') - Copy the beam(s) from the specified image to this image. If multiple beams, the current image must be able to hold a beam set of the shape in the specified image.

Returns

bool

Examples

ia.maketestimage('hcn',overwrite=True)      
rb = ia.restoringbeam()        # returns beam in record
print rb
#{'major': {'unit': 'arcsec', 'value': 53.500004857778549},
# 'minor': {'unit': 'arcsec', 'value': 34.199998900294304},
# 'positionangle': {'unit': 'deg', 'value': 6.0}}
rb['minor']['value'] = 12.5
# new beam specified in record
# NOTE This will not work for an image with multiple beams
ia.setrestoringbeam(beam=rb)   
print ia.restoringbeam()
#{'major': {'unit': 'arcsec', 'value': 53.500004857778549},
# 'minor': {'unit': 'arcsec', 'value': 12.5},
# 'positionangle': {'unit': 'deg', 'value': 6.0}}

# beam specified using parameter
# NOTE This will only work for an image with a traditional beam
ia.setrestoringbeam(major='36arcsec') 
print ia.restoringbeam()
#{'major': {'unit': 'arcsec', 'value': 36.0},
# 'minor': {'unit': 'arcsec', 'value': 12.5},
# 'positionangle': {'unit': 'deg', 'value': 6.0}}
ia.setrestoringbeam(remove=True)
print ia.restoringbeam()
#{}
ia.setrestoringbeam(major='53.5arcsec',minor='34.2arcsec',pa='6deg')
print ia.restoringbeam()
#{'major': {'unit': 'arcsec', 'value': 53.5},
# 'minor': {'unit': 'arcsec', 'value': 34.200000000000003},
# 'positionangle': {'unit': 'deg', 'value': 6.0}}
ia.close()

# Copy all beams from an image with multiple beams to another
# image with the same number of channels and polarizations

ia.open("multibeam.im")
ib = iatool()
ib.open("target.im")

# ensure target has no beam(s) at start, not always necessary
# but it doesn't hurt to do it.
ib.setrestoringbeam(remove=True)
# Now copy the beams. This only will work correctly if both images
# have the same number of channels and polarizations. nchan is set to
# the number of channels and npol is set to the number of polarizations
for c in range(nchan):
    for p in range(npol):
        beam = ia.restoringbeam(channel=c, polarization=p)
        ib.setrestoringbeam(beam=beam, channel=c, polarization=p)

ia.done()
ib.done()
shape()[source]

The shape of an image is a vector holding the length of each axis of the image. Although this information is also available in the summary function, it is so useful that it can be obtained directly. Both Float and Complex valued images are supported.

statistics(axes=[-1], region='', mask='', includepix=[-1], excludepix=[-1], list=False, force=False, disk=False, robust=False, verbose=False, stretch=False, logfile='', append=True, algorithm='classic', fence=-1, center='mean', lside=True, zscore=-1, maxiter=-1, clmethod='auto', niter=3)[source]

This method computes statistics from the pixel values in the image. You can then list them and retrieve them (into a record) for further analysis. This method supports only real valued images.

The names of the fields in the returned record are summarized below:

begin{itemize}

item {stfaf npts} - the number of unmasked points used

item {stfaf sum} - the sum of the pixel values: $sum I_i$

item {stfaf flux} - flux or flux density, see below for details

item {stfaf sumsq} - the sum of the squares of the pixel values: $sum I_i^2$

item {stfaf mean} - the mean of pixel values: $bar{I} = sum I_i / n$

item {stfaf sigma} - the standard deviation about the mean: $sigma^2 = (sum I_i - bar{I})^2 / (n-1)$

item {stfaf rms} - the root mean square: $sqrt {sum I_i^2 / n}$

item {stfaf min} - minimum pixel value

item {stfaf max} - the maximum pixel value

item {stfaf median} - the median pixel value (if {stfaf robust=T})

item {stfaf medabsdevmed} - the median of the absolute deviations from the median (if {stfaf robust=T})

item {stfaf quartile} - the inter-quartile range (if {stfaf robust=T}). Find the points which are 25% largest and 75% largest (the median is 50% largest).

item {stfaf q1} - The first quartile. Reported only if robust=T.

item {stfaf q3} - The third quartile. Reported only if robust=T.

item {stfaf blc} - the absolute pixel coordinate of the bottom left corner of the bounding box of the region of interest. If ‘region’ is unset, this will be the bottom left corner of the whole image.

item {stfaf blcf} - the formatted absolute world coordinate of the bottom left corner of the bounding box of the region of interest.

item {stfaf trc} - the absolute pixel coordinate of the top right corner of the bounding box of the region of interest.

item {stfaf trcf} - the formatted absolute world coordinate of the top right corner of the bounding box of the region of interest.

item {stfaf minpos} - absolute pixel coordinate of minimum pixel value

item {stfaf maxpos} - absolute pixel coordinate of maximum pixel value

item {stfaf minposf} - formatted string of the world coordinate of the minimum pixel value

item {stfaf maxposf} - formatted string of the world coordinate of the maximum pixel value

end{itemize}

The last four fields only appear if you evaluate the statistics over all axes in the image. As an example, if the returned record is captured in {stfaf ‘mystats’}, then you could access the ‘mean’ field via {cf print mystats[‘mean’]}.

If there are no good points (e.g. all pixels are masked bad in the region), then the length of these fields will be 0 (e.g. {cf len(mystats[‘npts’])==0}).

You have no control over which statistics are listed to the logger, you always get the same selection. You can choose to list the statistics or not (argument {stfaf list}).

As well as the simple (and faster to calculate) statistics like means and sums, you can also compute some robust (quantile-like) statistics. Currently these are the median, median absolute deviations from the median, the first and third quartiles, and the inner-quartile range. Because these are computationally expensive, they are only computed if robust=True.

Note that if the axes are set to all of the axes in the image (which is the default) there is just one value per statistic.

You have control over which pixels are included in the statistics computations via the {stfaf includepix} and {stfaf excludepix} arguments. These vectors specify a range of pixel values for which pixels are either included or excluded. They are mutually exclusive; you can specify one or the other, but not both. If you only give one value for either of these, say {stfaf includepix=b}, then this is interpreted as {stfaf includepix=[-abs(b),abs(b)]}.

This function generates a ‘storage’ lattice, into which the statistics are written as a means to improve performance on successive identical runs. After the initial execution of this method, it is only regenerated as necessary. For example, if you run the method twice with identical arguments, the statistics will be directly retrieved from the storage lattice the second time, and the statistics will not be recomputed. However, you can force regeneration of the statistics if you set force=True. VERY IMPORTANT NOTE. If you have an open image tool on which you’ve run statistics(), and change the pixel values of the opened image via another tool or a task, the opened image tool has no knowledge that pixel values have been changed, and so if you run ia.statistics() again with that tool, you will very likely get incorrect results since they come from the statistics stored from a previous run. First of all, it is highly discouraged that you do anything like this. Tools maintain state that only they know about internally; if you change an image that is already opened with another image tool or task, the general outcome will be undefined when you run methods on the already opened tool. However, if for some reason you must do this, first consider changing your algorithm because there should be no reason to have to do this. But, if you decide to continue down that path that will likely lead you over the edge of a cliff, then you can set force=True to ensure statistics are always recomputed.

The storage medium is either in memory or on disk, depending upon its size. You can force it to disk if you set {stfaf disk=T}, otherwise it decides for itself.

CURSOR AXES The axes parameter allows one to set the cursor axes over which statistics are computed. For example, consider a 3-dimensional image for which axes=[0,2]. The statistics would be computed for each XZ (axes 0 and 2) plane in the image. One could then examine those statistics as a function of the Y (axis 1) axis.

Each statistic is stored in an array in its own field in the returned dictionary. The dimensionality of these arrays is equal to the number of axes over which the statistics were not evaluated (called the display axes). For example, if the input image has four axes, and axes=[0], the output statistic arrays will have three dimensions. If axes=[0, 1], the output statistic arrays will have two dimensions.

The shape of the output arrays when axes has a positive number of elements is based on the region selection. If there is no region selection, the shape of the statistic arrays is just the shape of the image along the display (non-cursor) axes. For example, if the input image has dimensions of 300x400x4x80 and axes=[0, 1], in the absence of a region selection, the shape of the output statistic arrays will be 4x80. If there is a region selection, the shape of the output statistic arrays will be determined by the number of planes along the display axes chosen in the region selection. For example, continuing with our example, if axes=[0,1] and region=rg.box([0, 0, 1, 20], [299,399, 2, 60]), the output statistic arrays will have shapes of 2x41. Only the selected planes will be displayed in the logger output if verbose=True.

In the case where the image has a pixel mask, and/or the mask parameter is specified, and because of this specification a plane is entirely masked, this element is included in the statistic arrays (usually with a value of 0). It is not included in the logger output if verbose=True. One can exclude such elements from computations on the output arrays by using the numpy.extract() method. For example, to compute the minimum rms value, not including any fully masked planes, one could use

stats = ia.statistics(…) rmsmin = numpy.min(numpy.extract(stats[‘npts’]>0, stats[‘rms’]))

Thus in the computation of rmsmin, only the rms elements are considered which have associated values of npts that are not zero.

ALGORITHMS

Several types of statistical algorithms are supported:

  • classic: This is the familiar algorithm, in which all unmasked pixels, subject to any specified pixel ranges, are used. One may choose one of two methods, which vary only by performance, for computing classic statistics, via the clmethod parameter. The “tiled” method is the old method and is fastest in cases where there are a large number of individual sets of statistics to be computed and a small number of data points per set. This can occur when one sets the axes parameter, which causes several individual sets of statistics to be computed. The “framework” method uses the new statistics framework to compute statistics. This method is fastest in the regime where one has a small number of individual sets of statistics to calculate, and each set has a large number of points. For example, this method is fastest when computing statistics over an entire image in one go (no axes specified). A third option, “auto”, chooses which method to use by predicting which be faster based on the number of pixels in the image and the choice of the axes parameter.

  • fit-half: This algorithm calculates statistics on a dataset created from real and virtual pixel values. The real values are determined by the input parameters center and lside. The parameter center tells the algorithm where the center value of the combined real+virtual dataset should be. Options are the mean or the median of the input image’s pixel values, or at zero. The lside parameter tells the algorithm on which side of this center the real pixel values are located. True indicates that the real pixel values to be used are <= center. False indicates the real pixel values to be used are >= center. The virtual part of the dataset is then created by reflecting all the real values through the center value, to create a perfectly symmetric dataset composed of a real and a virtual component. Statistics are then calculated on this resultant dataset. These two parameters are ignored if algorithm is not “fit-half”. Because the maximum value is virtual if lside is True and the minimum value is virtual if lside is False, the value of the maximum position (if lside=True) or minimum position (if lside=False) is not reported in the returned record.

  • hinges-fences: This algorithm calculates statistics by including data in a range between Q1 - f*D and Q3 + f*D, inclusive, where Q1 is the first quartile of the distribution of unmasked data, subject to any specified pixel ranges, Q3 is the third quartile, D = Q3 - Q1 (the inner quartile range), and f is the user-specified fence factor. Negative values of f indicate that the full distribution is to be used (ie, the classic algorithm is used). Sufficiently large values of f will also be equivalent to using the classic algorithm. For f = 0, only data in the inner quartile range is used for computing statistics. The value of fence is silently ignored if algortihm is not “hinges-fences”.

  • chauvenet: The idea behind this algorithm is to eliminate outliers based on a maximum z-score value. A z-score is the number of standard deviations a point is from the mean of a distribution. This method thus is meant to be used for (nearly) normal distributions. In general, this is an iterative process, with successive iterations discarding additional outliers as the remaining points become closer to forming a normal distribution. Iterating stops when no additional points lie beyond the specified zscore value, or, if zscore is negative, when Chauvenet’s criterion is met (see below). The parameter maxiter can be set to a non-negative value to prematurely abort this iterative process. When verbose=True, the “N iter” column in the table that is logged represents the number of iterations that were executed.

    Chauvenet’s criterion allows the target z-score to decrease as the number of points in the distribution decreases on subsequent iterations. Essentially, the criterion is that the probability of having one point in a normal distribution at a maximum z-score of z_max must be at least 0.5. z_max is therefore a function of (only) the number of points in the distrbution and is given by

    npts = 0.5/erfc(z_max/sqrt(2))

    where erfc() is the complementary error function. As iterating proceeds, the number of remaining points decreases as outliers are discarded, and so z_max likewise decreases. Convergence occurs when all remaining points fall within a z-score of z_max. Below is an illustrative table of z_max values and their corresponding npts values. For example, it is likely that there will be a 5-sigma “noise bump” in a perfectly noisy image with one million independent elements.

    z_max npts 1.0 1 1.5 3 2.0 10 2.5 40 3.0 185 3.5 1,074 4.0 7,893 4.5 73,579 5.0 872,138 5.5 13,165,126 6.0 253,398,672 6.5 6,225,098,696 7.0 195,341,107,722

  • biweight: The biweight algorithm is a robust iterative algorithm that computes two quantities called the “location” and the “scale”, which are analogous to the mean and the standard deviation. In this case, the only keys present in the returned dictionary are ‘mean’ (location), ‘sigma’ (scale), ‘npts’, ‘min’, and ‘max’. The last three represent the values using the entire distribution. Note that the biweight algorithm does not support computation of quantile-like values (median, madm, q1, q3, and iqr), so setting robust=True will cause a warning message to be logged regarding that, and the computation will proceed.

    Important equations for the biweight algorithm are

    1. How to compute u_i values, which are related to the weights w_i = (1 - u_i*u_i), using the equation

      u_i = (x_i - c_bi)/(6.0*s_bi)                  (1)

      where x_i are the data values, c_bi is the biweight location and s_bi is the biweight scale. For the initial computation of the u_i values, c_bi is set equal to the median of the distribution and s_bi is set equal to the normalized median of the absolute deviation about the median (that is the median of the absolute deviation about the median multiplied by the value of the probit function at 0.75).

    2. The location, c_bi, is computed from

      c_bi = sum(x_i * w_i^2)/sum(w_i^2)             (2)

      where only values of u_i which satisfy abs(u_i) < 1 (w_i > 0) are used in the sums.

    3. The scale value is computed using:

      .             n * sum((x_i - c_bi)^2 * w_i^4)
      . s_bi^2 =    _______________________________    (3)
      .                    p * max(1, p - 1)               
      

    where n is the number of points for the entire distribution (which includes all the data, even for which abs(u_i) >= 1) and p is given by:

    p = abs(sum((w_i) \* (5\*w_i - 4)))
    

    Again, the sums include only data for which abs(u_i) < 1.

    The algorithm proceeds as follows.
    1. Compute initial u_i values (and hence w_i values) from equation (1), setting c_bi equal to the median of the distribution and s_bi equal to the normalized median of the absolute deviation about the median.

    2. Compute the initial value of the scale using the w_i values computed in step 1. using equation 3.

    3. Recompute u_i/w_i values using the most recent previous scale and location values.

    4. Compute the location using the u_i.w_i values from step 3 and equation (2).

    5. Recompute u_i/w_i values using the most recent previous scale and location values.

    6. Compute the new scale value using the the u_i/w_i values computed in step 5 and the value of the location computed in step 4.

    7. Steps 3. - 6. are repeated until convergence occurs or the maximum number of iterations (specified in the niter parameter) is reached. The convergence criterion is given by

      abs(s_bi - s_bi,prev) < 0.03 * sqrt(0.5/(n - 1))

      where s_bi,prev is the value of the scale computed in the previous iteration.

    In the special case where niter is specified to be negative, the faster, non-iterative algorithm proceeds as follows:

    1. Compute u_i/w_i values using the median for the location and the normalized median of the absolute deviation about the median as the scale

    2. Compute the location and scale (which can be carried out simultaneously) using the u_i/w_i values computed in step 1. The value of the location is just the median that is used in equation (3) to compute the scale

NOTES ON FLUX DENSITIES AND FLUXES

Fluxes and flux densities are not computed if any of the following conditions is met:

  1. The image does not have a direction coordinate

  2. The image does not have a intensity-like brightness unit. Examples of such units are Jy/beam (in which case the image must also have a beam) and K.

  3. There are no direction axes in the cursor axes that are used.

  4. If the (specified region of the) image has a non-degenerate spectral axis, and the image has a tablular spectral axis (axis with varying increments)

  5. Any axis that is not a direction nor a spectral axis that is included in the cursor axes is not degenerate within in the specified region

Note that condition 4 may be removed in the future.

In cases where none of the above conditions is met, the flux density(ies) (intensities integrated over direction planes) will be computed if any of the following conditions are met:

  1. The image has no spectral coordinate

  2. The cursor axes do not include the spectral axis

  3. The spectral axis in the chosen region is degenerate

In the case where there is a nondegenerate spectral axis that is included in the cursor axes, the flux (flux density integrated over spectral planes) will be computed. In this case, the spectral portion of the flux unit will be the velocity unit of the spectral coordinate if it has one (eg, if the brightness unit is Jy/beam and the velocity unit is km/s, the flux will have units of Jy.km/s). If not, the spectral portion of the flux unit will be the frequency unit of the spectral axis (eg, if the brightness unit is K and the frequency unit is Hz, the resulting flux unit will be K.arcsec2.Hz).

In both cases of flux density or flux being computed, the resulting numerical value is assigned to the “flux” key in the output dictionary.

If the image has units of Jy/beam, the flux density is just the mean intensity multiplied by the number of beam areas included in the region. The beam area is defined as the volume of the elliptical Gaussian defined by the synthesized beam, divided by the maximum of that function, which is equivalent to

pi/(4*ln(2)) * major * minor

where ln() is the natural logarithm and major and minor are the major and minor FWHM axes of the beam, respectively.

Parameters

  • axes (intVec=[-1]) - List of axes to evaluate statistics over. Default is all axes.

  • region ({string, string, record}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • includepix (doubleVec=[-1]) - Range of pixel values to include. Vector of 1 or 2 doubles. Default is to include all pixels.

  • excludepix (doubleVec=[-1]) - Range of pixel values to exclude. Vector of 1 or 2 doubles. Default is exclude no pixels.

  • list (bool=False) - If True print bounding box and statistics to logger.

  • force (bool=False) - If True then force the stored statistical accumulations to be regenerated

  • disk (bool=False) - If T then force the storage image to disk

  • robust (bool=False) - If T then compute robust statistics as well

  • verbose (bool=False) - If T then log statistics

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

  • logfile (string='') - Name of file to which to write statistics.

  • append (bool=True) - Append results to logfile? Logfile must be specified. Default is to append. False means overwrite existing file if it exists.

  • algorithm (string='classic') - Algorithm to use. Supported values are “biweight”, “chauvenet”, “classic”, “fit-half”, and “hinges-fences”. Minimum match is supported.

  • fence (double=-1) - Fence value for hinges-fences. A negative value means use the entire data set (ie default to the “classic” algorithm). Ignored if algorithm is not “hinges-fences”.

  • center (string='mean') - Center to use for fit-half. Valid choices are “mean”, “median”, and “zero”. Ignored if algorithm is not “fit-half”.

  • lside (bool=True) - For fit-half, real data are <=; center? If False, real data are >= center. Ignored if algorithm is not “fit-half”.

  • zscore (double=-1) - For chauvenet, this is the target maximum number of standard deviations data may have to be included. If negative, use Chauvenet’s criterion. Ignored if algorithm is not “chauvenet”.

  • maxiter (int=-1) - For chauvenet, this is the maximum number of iterations to attempt. Iterating will stop when either this limit is reached, or the zscore criterion is met. If negative, iterate until the zscore criterion is met. Ignored if algorithm is not “chauvenet”.

  • clmethod (string='auto') - Method to use for calculating classical statistics. Supported methods are “auto”, “tiled”, and “framework”. Ignored if algorithm is not “classic”.

  • niter (int=3) - For biweight, this is the maximum number of iterations to attempt. Iterating will stop when either this limit is reached, or the zscore criterion is met. If negative, do a fast, simple computation (see description). Ignored if the algorithm is not “chauvenet”.

Returns

record

Examples

#
print "\t----\t statistics Ex 1 \t----"
ia.maketestimage()
ia.statistics()
ia.close()
#

# evaluate statistics for each spectral plane in an ra x dec x frequency image
ia.fromshape("", [20,30,40])
# give pixels non-zero values
ia.addnoise()
# These are the display axes, the calculation of statistics occurs
# for each (hyper)plane along axes not listed in the axes parameter,
# in this case axis 2 (the frequency axis)
# display the rms for each frequency plane (your mileage will vary with
# the values).
stats = ia.statistics(axes=[0,1])
 stats["rms"]
  Out[10]: 
array([ 0.99576014,  1.03813124,  0.97749186,  0.97587883,  1.04189885,
        1.03784776,  1.03371549,  1.03153074,  1.00841606,  0.947155  ,
        0.97335404,  0.94389403,  1.0010221 ,  0.97151822,  1.03942156,
        1.01158476,  0.96957082,  1.04212773,  1.00589049,  0.98696715,
        1.00451481,  1.02307892,  1.03102005,  0.97334671,  0.95209879,
        1.02088714,  0.96999902,  0.98661619,  1.01039267,  0.96842754,
        0.99464947,  1.01536798,  1.02466023,  0.96956468,  0.98090756,
        0.9835844 ,  0.95698935,  1.05487967,  0.99846411,  0.99634868])




In this example, we ask to see statistics evaluated over the
entire image. 
subimage(outfile='', region='', mask='', dropdeg=False, overwrite=False, list=True, stretch=False, wantreturn=True, keepaxes='')[source]

This function copies all or part of the image to another on-the-fly Image tool. Both float and complex valued images are supported.

If outfile is given, the subimage is written to the specified disk file. If outfile is unset, the returned Image  actually references the input image file (i.e. that associated with the Image  to which you are applying this function). So if you deleted the input image disk file, it would render this  useless. When you destroy this  (with the done function) the reference connection is broken.

Sometimes it is useful to drop axes of length one (degenerate axes). Use the dropdeg argument if you want to do this. Further control is provided via the keepaxes parameter. If dropdeg=True, you may specify a list of degenerate axes to keep in the keep axes parameter. This allows you to drop only a subset of degenerate axes. This parameter is ignored if dropdeg=False. If dropdeg=True, all degenerate axes are dropped if keepaxes is set to an empty list (this is the default behavior). Nondegenerate axes are implicitly kept, even if they are included in the keepaxes list.

The output mask is the combination (logical OR) of the default input  (if any) and the OTF mask. Any other input  will not be copied. Use function maskhandler if you need to copy other masks too.

If the mask has fewer dimensions than the image and if the shape of the dimensions the mask and image have in common are the same, the mask will automatically have the missing dimensions added so it conforms to the image.

If stretch is True and if the number of mask dimensions is less than or equal to the number of image dimensions and some axes in the mask are degenerate while the corresponding axes in the image are not, the mask will be stetched in the degenerate dimensions. For example, if the input image has shape [100, 200, 10] and the input mask has shape [100, 200, 1] and stretch is True, the mask will be stretched along the third dimension to shape [100, 200, 10]. However if the mask is shape [100, 200, 2], stretching is not possible and an error will result.

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • region (variant='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • dropdeg (bool=False) - Drop degenerate axes

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • list (bool=True) - List informative messages to the logger

  • stretch (bool=False) - Stretch the mask if necessary and possible?

  • wantreturn (bool=True) - Return an image analysis tool attached to the created subimage

  • keepaxes (intVec='') - If dropdeg=True, these are the degenerate axes to keep. Nondegenerate axes are implicitly always kept.

Returns

image

Examples

#
print "\t----\t subimage Ex 1 \t----"
ia.maketestimage('myfile',overwrite=True)
im2 = ia.subimage()                # a complete copy
r1 = rg.box([10,10],[30,40],[5,5]) # A strided pixel box region
im3 = ia.subimage(outfile='/tmp/foo', region=r1, overwrite=True)
                                   # Explicitly named subimage
im2.done()
im3.done()
ia.close()
#


# As an example of the usage of the keepaxes parameter, consider an image
# that has axes RA, Dec, Stokes, and Freq. The Stokes and Freq axes are both
# degenerate while the RA and Dec axes are not, and it is desired to make a
# subimage in which the Stokes axis is discarded. The following command will
# accomplish that.
ia.open("my.im")
subim = ia.subimage(outfile="discarded_stokes.im", dropdeg=True, keepaxes=[3])
ia.done()
subim.done()
summary(doppler='RADIO', list=True, pixelorder=True, verbose=False)[source]

This function summarizes various metadata such as shape, Coordinate System, restoring beams, and masks.

If called without any arguments, this function displays a summary of the image metadata to the logger; where appropriate, values will be formatted nicely (e.g. HH:MM:SS.SS for the reference value of RA axes).

For spectral axes, the information is listed in both the velocity and frequency domains. The doppler parameter allows one to specify what velocity doppler convention it is listed in. Supported values are: “radio”, “optical”, and “True”. Alternative values are “z” for “optical”, and “beta” or “relativistic” for True. The default is “radio”. The definitions are

  • radio: \(1 - F\)

  • optical: \(-1 + 1/F\)

  • True: \((1 - F^2)/(1 + F^2)\)

where \(F = \nu/\nu_0\) and \(\nu_0\) is the rest frequency. If the rest frequency has not been set in your image, you can set it via a coordinate system (cs) tool using the setrestfrequency() method().

The keys of the returned dictionary are

ndim Dimension of the image. shape Length of each axis in the image. tileshape Shape of the chunk which is most efficient for I/O. axisnames Name of each axis. refpix Reference pixel for each axis (0-relative) refval Reference value for each axis. incr Increment for each axis. axisunits Unit name for each axis. unit Brightness units for the pixels. hasmask True if the image has a mask. defaultmask The name of the mask which is applied by default. masks The names of all the masks stored in this image. restoringbeam The restoring beam(s) if present. imagetype The image type.

For an image with multiple beams, the restoringbeam field is a dictionary of dictionaries with keys of names “” followed by the channel number, if the image has a spectral coordinate, or the polarization number if it does not. That is, the keys have names “”, “”, “”, etc. If the image has both a spectral and a polarization coordinate, each of these dictionaries is a dictionary with keys of the same form which range from 0 to the number of polarizations minus 1; “”, “”, … The dictionaries pointed to by the channel and/or polarization number contain information for the beam at that position.

If the list parameter is set to False, then the summary will not be written to the logger. The return value of the method, in the “header” field is a vector string containing the formatted output that would have been logged in the list=True case.

If verbose is True and the image contains multiple beams, the formatted output, whether it is written to the logger or placed in the output record, will have information on every beam in the dataset. If verbose=False and the image has multiple beams, only a summary of beams for each polarization is listed. In this case, the beams with the maximum area, the minimum area, and the median area for each polarization are listed. However, all the beams can still be found in the restoringbeam field of the returned dictionary. If the image does not have multiple beams, verbose is not used.

Parameters

  • doppler (string='RADIO') - If there is a spectral axis, list velocity too, with this doppler definition

  • list (bool=True) - List the summary to the logger

  • pixelorder (bool=True) - List axis descriptors in pixel or world axis order

  • verbose (bool=False) - Give a full listing of beams or just a short summary? Only used when the image has multiple beams.

Returns

record

Examples

#
print "\t----\t summary Ex 1 \t----"
ia.maketestimage('myim1', overwrite=True)
ia.summary()                 # summarize to logging only
s = ia.summary(list=False)       # store header in record
if s['header']['ndim'] == 2: # program using header values
  print s['header']['axisnames']
ia.close()
#
tofits(outfile='', velocity=False, optical=True, bitpix=-32, minpix=1, maxpix=-1, region='', mask='', overwrite=False, dropdeg=False, deglast=False, dropstokes=False, stokeslast=True, wavelength=False, airwavelength=False, stretch=False, history=True)[source]

This function converts the image into a  file.

If the image has a rest frequency associated with it, it will always write velocity information into the  file. By default the frequency information will be primary as it is the internal native format. If you select velocity=T then by default the velocity is written in the optical convention, but if optical=F it will use the radio convention instead. Alternatively, if you use velocity=F and wavelength=T, the spectral axis will be written in wavelength.

The  definition demands equal increment pixels. Therefore, if you write wavelength or optical velocity information as primary, the increment is computed at the spectral reference pixel. If the bandwidth is large, this may incur non-negligible coordinate calculation errors far from the reference pixel if the spectral bins are not originally equidistant in wavelength. Images generated by the CASA clean task have spectral axes which are always equidistant in frequency.

By default the image is written as a floating point  file (bitpix= -32). Under rare circumstances you might want to save space and write it as scaled 16 bit integers (bitpix = 16). You can have tofits calculate the scaling factors by using the default minpix and maxpix. If you set minpix and maxpix, values outside of that range will be truncated. This can be useful if all of the  images dynamic range is being used by a few high or low values and you are not interested in preserving those values exactly. Besides the factor of two space savings you get by using 16 instead of 32 bits, integer images usually also compress well (for example, with the standard GNU software facility gzip).

If the specified  extends beyond the image, it is truncated.

The output mask is the combination (logical OR) of the default input  (if any) and the OTF mask.

Sometimes it is useful to drop axes of length one (degenerate axes) because not all FITS readers can handle them. Use the dropdeg argument if you want to do this. If you want to specifically only drop a degenerate Stokes axis, use the dropstokes argument.

If you want to place degenerate axes last in the FITS header, use the deglast argument. If you want to make sure that the Stokes axis is placed last in the FITS header, use the stokeslast argument.

Parameters

  • outfile (string='') - FITS file name. Default is input name + ‘.fits’

  • velocity (bool=False) - prefer velocity (rather than frequency) as primary spectral axis?

  • optical (bool=True) - use the optical (rather than radio) velocity convention?

  • bitpix (int=-32) - Bits per pixel, -32 (floating point) or 16 (integer)

  • minpix (double=1) - Minimum pixel value for BITPIX=16, Default is to autoscale if minpix > maxpix.

  • maxpix (double=-1) - Maximum pixel value for BITPIX=16, Default is to autoscale if maxpix < minpix.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • dropdeg (bool=False) - Drop degenerate axes?

  • deglast (bool=False) - Put degenerate axes last in header?

  • dropstokes (bool=False) - Drop Stokes axis?

  • stokeslast (bool=True) - Put Stokes axis last in header?

  • wavelength (bool=False) - Write spectral axis in units of wavelength (instead of velocity or frequency)?

  • airwavelength (bool=False) - When wirting the spectral axis in units of wavelength, use air wavelength instead of vacuum wavelength?

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

  • history (bool=True) - Write the image history to the FITS file? Default True

Returns

bool

Examples

#
print "\t----\t tofits Ex 1 \t----"
ia.maketestimage()
ok = ia.tofits('MYFILE.FITS',overwrite=True)
                       # write FITS image file
ok = ia.tofits('MYFILE2.FITS', bitpix=16, overwrite=True)
                       # Write as scaled 16 bit integers
ia.close()
#
topixel(value='')[source]

This method converts from world to pixel coordinates. The world coordinate can be provided in many formats (numeric, string, quantum etc.) via the value parameter. These match the output of the toworld() method.

This function is just a wrapper for the coordsys tool method of the same name, so see that documentation for a description and more examples.

Parameters

  • value (variant='') - Absolute world coordinate, Numeric vector, vector of strings representing quantities, or record of format analogous to that produced by ia.toworld(). Default is reference value.

Returns

record

Examples

#
print "\t----\t topixel Ex 1 \t----"
ia.maketestimage();
w = ia.toworld([10,10], 'n')        # Numeric vector
ia.topixel(w)
#{'ar_type': 'absolute',
# 'numeric': array([10., 10.]), 'pw_type': 'pixel'}
w = ia.toworld([10,10], 'm')        # Record of measures
ia.topixel(w)
#{'ar_type': 'absolute',
# 'numeric': array([10., 10.]), 'pw_type': 'pixel'}
ia.close()
#


Convert a pixel coordinate to world as floats and then
back to pixel.  Do the same with the world coordinate
formatted as measures instead.
torecord()[source]

You can convert an associated image to a record for manipulation or passing it to inputs of other methods of other tools. This method and fromrecord() are used for serialization and deserialization.

toworld(value='', format='n', dovelocity=True)[source]

This method converts between pixel and world coordinates. A variety of return formats is supported. If format=’n’, numerical values are returned. If format=’q’, values formatted as quantities are returned. If format=’s’, values formatted as strings are returned. If format=’m’, values formatted as measures are returned. If format=’m’, one can choose to have the corresponding velcocites of an extant spectral coordinate computed as well by specifyting dovelocity=True (dovelocity is ignored if format is not equal to ’m’ or if the image does not have a spectral coordinate).

Parameters

  • value (variant='') - Absolute pixel coordinate. Numeric vector is allowed. Default is reference pixel.

  • format (string='n') - What type of formatting? String from combination of ‘n’ (numeric), ‘q’ (quantity), ‘m’ (measure), ‘s’ (string).

  • dovelocity (bool=True) - Compute corresponding spectral velocities if format=’m’?

Returns

record

Examples

#
print "\t----\t toworld Ex 1 \t----"
ia.maketestimage('hcn',overwrite=True)
w = ia.toworld([10,10], 'n')
print w
#{'numeric': array([ 0.00174533, -0.0015708 ])}
w = ia.toworld([10,10], 'nmq')
print w
#{'measure': {'direction': {'m0': {'unit': 'rad',
#                                  'value': 0.0017453323593185704},
#                           'm1': {'unit': 'rad',
#                                  'value': -0.0015707969259645381},
#                           'refer': 'J2000',
#                           'type': 'direction'}},
# 'numeric': array([ 0.00174533, -0.0015708 ]),
# 'quantity': {'\*1': {'unit': 'rad', 'value': 0.0017453323593185704},
#              '\*2': {'unit': 'rad', 'value': -0.0015707969259645381}}}
ia.close()
#


Convert to a vector of floats and then to a record
holding a vector of floats, a vector of quantities
and a subrecord of measures.
transpose(outfile='', order='')[source]

This method transposes the axes in the input image to the specified order. The associated pixel and mask values and coordinate system are transposed.

If the outfile parameter is empty, only a temporary image is created; no output image is written to disk.

The order parameter describes the mapping of the input axes to the output axes. It can be one of three types: a non-negative integer, a string, or a list of strings. If a string or non-negative integer, it should contain zero-based digits describing the new order of the input axes. It must contain the same number of (unique) digits as the number of input axes. For example, specifying reorder=”1032” or reorder=1032 for a four axes image maps input axes 1, 0, 3, 2 to output axes 0, 1, 2, 3. In the case of order being a nonnegative integer and the zeroth axis in the input being mapped to zeroth axis in the output, the zeroth digit is implicitly understood to be 0 so that to transpose an image where one would use a string order=”0321”, one could equivalently specify an int order=321. IMPORTANT: When specifying a non-negative integer and mapping the zeroth axis of the input to the zeroth axis of the output, do *not* explicitly specify the leading 0; eg, specify order=321 rather than order=0321. Python interprets an integer with a leading 0 as an octal number.

Because of ambiguity for axes numbers greater than nine, using string or integer order specifications cannot handle images containing more than 10 axes. The order parameter can also be specified as a list of strings which uniquely minimally match, ignoring case, the image axis names (ia.coordsys().names()). So to reorder an image with right ascension, declination, and frequency axes, one could specify order=[“d”, “f”, “r”] or equivalently [“decl”, “frequ”, “right a”]. Note that specifying “ra” for the right ascension axis will result in an error because “ra” does not match the first two characters of right ascension. Axes can be simultaneously inverted in cases where order is a string or an array of strings by specifying negative signs in front of the axis/axes to be inverted. So, in a 4-D image, order=”-10-3-2” maps input axes 1, 0, 3, 2 to output axes 0, 1, 2, 3 and reverses the direction and values of input axes 1, 3, and 2.

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • order (variant='') - Zero-based order of axes in output image (eg “120” $=>$ input$->$ output 0-$>$2, 1-$>$0, 2-$>$1))

Returns

image

Examples

# swap stokes (axis 2) and spectral (axis 3) axes in a 4 dimensional image
ia.open("myimage.fits")
reordim = ia.transpose(outfile="my_reordered_image.im", order="0132")
ia.done()
twopointcorrelation(outfile='', region='', mask='', axes=[-1], method='structurefunction', overwrite=False, stretch=False)[source]

This function computes two-point auto-correlation functions from an image.

By default, the auto-correlation function is computed for the Sky axes. If there is no sky in the image, then the first two axes are used. Otherwise you can specify which axes the auto-correlation function lags are computed over with the axes argument (must be of length 2).

Presently, only the Structure Function is implemented. This is defined as :

\[S(lx,ly) = \< (data(i,j) - data(i+lx,j+ly))^2 \>\]

where \(lx, ly\) are integer lags in the x (0-axis) and y (1-axis) directions. The ensemble average is over all the values at the same lag pair. This process is extremely compute intensive and so you may have to be patient.

In an auto-correlation function image there are some symmetries. The first and third quadrants are symmetric, and the second and fourth are symmetric. So in principle, all the information is in the top or bottom half of the image. We just write it all out to look nice. The long lags don’t have a lot of contributing values of course.

Parameters

  • outfile (string='') - Output image file name. Default is unset.

  • region ({record, string}='') - Region selection. Default is to use the full image.

  • mask (variant='') - Mask to use. Default is none.

  • axes (intVec=[-1]) - The pixel axes to compute structure function over. The default is sky or first two axes.

  • method (string='structurefunction') - The method of computation. String from ‘structurefunction’.

  • overwrite (bool=False) - Overwrite (unprompted) pre-existing output file?

  • stretch (bool=False) - Stretch the mask if necessary and possible? Default False

Returns

bool

Examples

ia.maketestimage("test.image")
ia.twopointcorrelation("2pt.image")
ia.done()
type()[source]

This function returns the string ’image’. It can be used in a script to make sure this variable is an Image .

unlock()[source]

This function releases any lock set on the  (and also flushes any outstanding I/O to disk). It is not of general user interest. It can be useful in scripts when a file is being shared between more than one process. See also functions lock and haslock.