pipeline.extern.findContinuum

This file implements the algorithm to determine continuum channel ranges to use from an ALMA image cube (dirty or clean) in CASA format. All dependencies on what were originally analysisUtils functions have been pasted into this file for convenience. This function is meant to be run inside CASA. Simple usage:

import findContinuum as fc fc.findContinuum(‘my_dirty_cube.image’)

This file can be found in a typical pipeline distribution directory, e.g.: /lustre/naasc/sciops/comm/rindebet/pipeline/branches/trunk/pipeline/extern As of March 7, 2019 (version 3.36), it is compatible with both python 2 and 3.

Code changes for Pipeline2020: (as of August 9, 2020) 0) fix for PIPE-554 (bug found in cube imaging of calibrators by ARI-L project) 1) fix for PIPE-525 (divide by zero in a 130-target MOUS) 2) new feature PIPE-702 (expand mask if mom8fc image has emission outside it) 3) Add warning messages for too little bandwidth or too little spread. 4) if outdir is specified and does not exist, stop and print error message. 5) Modified the is_binary() function to work in both python 2 and 3. 6) Add LSRK frequency ranges in a new final line of the .dat file 7) Increase to 2 significant digits on sigmaFC in the png name 8) Allow manually specified narrow value to be used everywhere instead of ‘auto’ 9) Print Final selection message to casa log 10) fix bug in convertSelectionIntoChannelList (only relevant if avoidance used) 11) allow single value in convertSelection 12) fix bug in avoidance range when gaps were trimmed (not relevant to PL) 13) remove warning about no beam in cube header 14) Add totalChannels & medianWidthOfRanges parameters+logic to pickAutoTrimChannels() 15) Insure that channel ranges are in increasing order if a second range is added to a solo range 16) add imstatListit=False to all imstat calls 17) add minor ticks to lower x-axis (channel number) of the plot 18) added 21 new control parameters to findContinuum (default values listed) for a total of 98:

  • returnWarnings=False (PL 2020 might set this to True if Dirk has time)

  • sigmaFindContinuumMode=’auto’ (split from sigmaFindContinuum)

  • returnSigmaFindContinuum=False

  • enableRejectNarrowInnerWindows=True (to be able to disable it manually)

  • avoidExtremaInNoiseCalcForJointMask=False

  • buildMom8fc=True

  • momentdir=’’

  • amendMaskIterations=’auto’

  • skipchan=1

  • checkIfMaskWouldHaveBeenAmended=False

  • fontsize=10

  • thresholdForSame=0.01 (to control the new 4-letter codes)

  • keepIntermediatePngs=True

  • enableOnlyExtraMask=True

  • useMomentDiff=True

  • smallBandwidthFraction=0.05

  • smallSpreadFraction=0.33

  • skipAmendMask=False (for manual usage)

  • useAnnulus=True

  • cubeSigmaThreshold=7.5

  • npixThreshold=7 (for onlyExtraMask)

  1. Added 27 new functions:

  • round_half_up (to support the same rounding in python 3 as 2)

  • byteDecode (to convert bytes to string for python 3)

  • imageSNR

  • imageSNRAnnulus

  • plotChannelSelections

  • compute4LetterCodeAndUpdateLegend

  • compute4LetterCode

  • cubeNoiseLevel

  • updateChannelRangesOnPlot

  • computeNpixCubeMedian

  • computeNpixMom8Median

  • computeNpixMom8MedianBadAtm

  • replaceLineFullRangesWithNoise

  • gatherWarnings

  • tooLittleBandwidth

  • tooLittleSpread

  • computeSpread

  • amendMaskYesOrNo

  • extraMaskYesOrNo

  • onlyExtraMaskYesOrNo

  • invertChannelRanges

  • robustMADofContinuumRanges

  • findWidestContiguousListInChannelRange

  • imagePercentileNoMask (not used by PL because avoidExtremaInNoiseCalcForJointMask=False)

  • combineContDat (not used by pipeline)

  • getSpwFromPipelineImageName (not used by pipeline)

  • getFieldnameFromPipelineImageName (not used by pipeline)

-Todd Hunter

Functions

CalcAtmTransmissionForImage(img, imageInfo)

This function is called by atmosphereVariation.

ExpandYLimitsForLegend()

Called by runFindContinuum.

MAD(a[, c, axis])

This function is called by removeInitialQuadraticIfNeeded, findContinuumChannels, meanSpectrum, computeStatisticalSpectrumFromMask, plotStatisticalSpectrumFromMask and runFindContinuum. Median Absolute Deviation along given axis of an array: median(abs(a - median(a))) / c c = 0.6745 is the constant to convert from MAD to std.

aboveBelow(avgSpectrumNansReplaced, threshold)

This function is called by runFindContinuum. Given an array of values (i.e. a spectrum) and a threshold value, this function computes and returns 6 items: * the number of channels above that threshold (group 1) * the number of channels below that threshold (group 2) * the ratio of these numbers (i.e., #ChanAboveThreshold / #ChanBelowThreshold)) * sum of the intensities in group (1) * sum of the intensities in group (2) * the ratio of these sums (i.e., FluxAboveThreshold / FluxBelowThreshold) (This value is not currently used by the calling function.) Example with 24 channels where 3 have positive spikes, and rest are zero: 3 2 3 channelsAbove = 3 +threshold=+1 …………………… sumAbove = (3-1) + (2-1) + (3-1) = 5 zero ————————— channelsBelow = 21 -threshold=-1 …………………… sumBelow = 21*(1-0) = 21.

allContinuumSelected(selection, nchan[, …])

Returns True if only one range is selected and that range covers sufficient bandwidth (default fraction = 92.5%) nchan: number of channels in the spectrum fraction: default = 0.925 which is 92.5%

amendMaskYesOrNo(badAtmosphere, median, …)

Used in Pipeline2020 going forward.

atmosphereVariation(img, imageInfo, chanInfo)

This function is called by findContinuum. Computes the absolute and percentage variation in atmospheric transmission and sky temperature across an image cube. Returns 5 values: max(Trans)-min(Trans), and as percentage of mean, Max(Tsky)-min(Tsky), and as percentage of mean, standard devation of Tsky.

avgOverCube(pixels[, useAbsoluteValue, …])

This function was used by Cycle 4 and 5 pipeline from meanSpectrum(), but will not be used in the Cycle 6 default heuristic of ‘mom0mom8jointMask’. Computes the average spectrum across a multi-dimensional array read from an image cube, ignoring any NaN values. Inputs: pixels: a 3D array (with degenerate Stokes axis dropped) threshold: value in Jy median: if True, compute the median instead of the mean mask: use pixels inside this spatial mask (i.e. with value==1 or True) to compute the average (the spectral mask is taken as the union of all channels) This is the mask located inside the image cube specified in findContinuum(img=’’) usermask: this array results from the mask specified in findContinuum(mask=’’) parameter If threshold is specified, then it only includes pixels with an intensity above that value. Returns: * average spectrum * percentage of pixels not masked Note: This function assumes that the spectral axis is the final axis. If it is not, there is no single setting of the axis parameter that can make it work.

buildMeanSpectrumFilename(img, …[, …])

This function is called by findContinuum and runFindContinuum.

byteDecode(buffer)

casaVersionCompare(comparitor, versionString)

This function is called by meanSpectrum.

casalogPost(mystring[, debug])

Generates an INFO message prepended with the version number of findContinuum.

centralArcsecArgumentMismatch(centralArcsec, …)

This function is called by checkForMismatch.

channelSelectionRangesToIndexArray(selection)

This function is called by runFindContinuum.

channelsInLargestGroup(selection)

This function is called by runFindContinuum.

checkForMismatch(meanSpectrumFile, img, …)

This function is called by runFindContinuum.

checkForTriangularWavePattern(img[, …])

+++++++ This function is not used for meanSpectrumMethod == ‘mom0mom8jointMask’.

combineContDat(contdatlist[, outputfile, …])

Takes a list of *.dat files created by findContinuum for a single field/spw and pulls the line that contains LSRK frequency ranges from each one, and writes it to a new cont.dat file suitable for the pipeline that contains an entry for each spw of that field. contdatlist: a list or comma-delimited string of the .dat files, or of the cube names (to which ‘_findContinuum.dat’ will be appended) fieldname: use this value if specified, otherwise parse from filename spw: integer or string integer, use this value if specified, otherwise parse from the filename Note: It assumes that none of the files originate from an ALLcont spw, so it will never write the “ALL” line that the pipeline triggers on. But this feature could be added by calling the following function if there is only one range of channels read: allContinuumSelected(selection, nchan, fraction=ALL_CONTINUUM_CRITERION) We would only need to know nchan.

compute4LetterCode(mom8fcPeak, mom8fcPeak0, …)

Used in Pipeline2020 going forward.

compute4LetterCodeAndUpdateLegend(…[, …])

Used in Pipeline2020 going forward.

computeBandwidth(selection, channelWidth[, loc])

This function is called by runFindContinuum and findContinuum.

computeMadSpectrum(img[, box])

+++++++ This function is not used for meanSpectrumMethod == ‘mom0mom8jointMask’.

computeMedianCorrectionFactor(baselineMode, …)

This function is called by findContinuumChannels and medianCorrected.

computeNpixCubeMedian(mom8fcMedian, …)

Used in Pipeline2020 going forward.

computeNpixMom8Median(mom8fcMedian, …[, …])

Used in Pipeline2020 going forward.

computeNpixMom8MedianBadAtm(mom8fcMedian, …)

Used in Pipeline2020 going forward.

computeSpread(selection, channelWidth)

Used in Pipeline2020 going forward.

computeStatisticalSpectrumFromMask(cube, …)

New heuristic for Cycle 6 pipeline. It is called by meanSpectrumFromMom0Mom8JointMask Uses ia.getprofile to compute the mean spectrum of a cube within a masked area. jointmask: a 2D or 3D mask image indicating which pixels shall be used or an expression with < or > along with a 2D or 3D mask image jointMaskForNormalize: a 2D or 3D mask image to be used when normalizeByMAD==True pbimg: the path to the primary beam of this cube, or a 2D primary beam with an expression statistic: passed to ia.getprofile via the ‘function’ parameter normalizeByMAD: if True, then create the inverse of the jointmask and normalize the spectrum by the spectrum of ‘xmadm’ (scaled MAD) computed on the inverse mask Returns: three arrays: channels, freqs(Hz), intensities, and a Boolean which says if normalization was applied.

convertChannelListIntoSelection(channels[, …])

This function is called by findContinuumChannels.

convertColonDelimitersToHMSDMS(mystring[, …])

This function is called by rad2radec. Converts HH:MM:SS.SSS, +DD:MM:SS.SSS to HHhMMmSS.SSSs, +DDdMMmSS.SSSs or HH:MM:SS.SSS +DD:MM:SS.SSS to HHhMMmSS.SSSs +DDdMMmSS.SSSs or HH:MM:SS.SSS, +DD:MM:SS.SSS to HHhMMmSS.SSSs, +DD.MM.SS.SSS or HH:MM:SS.SSS +DD:MM:SS.SSS to HHhMMmSS.SSSs +DD.MM.SS.SSS s: whether or not to include the trailing ‘s’ in both axes.

convertSelectionIntoChannelList(selection)

Converts an arbitrary CASA selection string into an array of channel numbers ‘0~1;4~6’ – > [0,4,5] Ignores any initial spw number that preceeds a colon.

countChannels(channels)

This function is called by runFindContinuum. Counts the number of channels in a CASA channel selection string. If multiple spws are found, then it returns a dictionary of counts: e.g. “1:5~20;30~40” yields 27; or ‘6~30’ yields 25 “1:5~20;30~40,2:6~30” yields {1:27, 2:25}.

countChannelsInRanges(channels[, separator])

This function is called by runFindContinuum.

countPixelsAboveZero(img[, useImstat, value])

This function is called by meanSpectrumFromMom0Mom8JointMask.

countUnmaskedPixels(img[, useImstat])

This function is called by meanSpectrumFromMom0Mom8JointMask.

createNoisyQuadratic([n, nsigma])

Function meant simply to test removeStatContQuadratic.

create_casa_quantity(myqatool, value, unit)

This function is called by CalcAtmTransmissionForImage, frames, and lsrkToRest.

cubeLSRKToTopo(img, imageInfo[, freqrange, …])

This function is called by CalcAtmTransmissionForImage and writeContDat. Reads the date of observation, central RA and Dec, and observatory from an image cube and then calls lsrkToTopo to return the specified frequency range in TOPO. freqrange: desired range of frequencies (empty string or list = whole cube) floating point list of two frequencies, or a delimited string (delimiter = ‘,’, ‘~’ or space) prec: in fractions of Hz (only used to display the value when verbose=True) vis: read date of observation from the specified measurement set (instead of from img).

cubeNoiseLevel(cube[, pbcube, mask, …])

Used in Pipeline2020 going forward. Computes median of the per-channel scaled MAD within the standard noise annulus and outside the specified mask chans: channel selection string (passed to imstat) mask: either mask image to use (e.g. the string passed to imstat will be ‘“mask” == 0’) or a proper mask string containing >, < or == percentile: compute this percentile value of the per-channel scaled MAD (examples: 50 = median, 90 = close to the highest value) if it is a list, then a list of values is returned Returns: median of MAD, median or [list of MADs], median.

drawYlabel(img, typeOfMeanSpectrum, …[, …])

This function is called by runFindContinuum.

extraMaskYesOrNo(badAtmosphere, median, …)

Used in Pipeline2020 going forward.

findContinuum([img, pbcube, psfcube, …])

This function calls functions to: 1) compute a representative ‘mean’ spectrum of a dirty cube 2) find the continuum channels 3) plot the results and writes a text file with channel ranges and frequency ranges It calls runFindContinuum up to 2 times.

findContinuumChannels(spectrum[, …])

This function is called by runFindContinuum. Trys to find continuum channels in a spectrum, based on a threshold or some number of edge channels and their median and standard deviation. Inputs: spectrum: a one-dimensional array of intensity values nBaselineChannels: number of channels over which to compute standard deviation/MAD sigmaFindContinuum: value to multiply the standard deviation by then add to median to get threshold. Default = 3. narrow: the minimum number of channels in a contiguous block to accept if 0<narrow<1, then it is interpreted as the fraction of all channels within identified blocks nanmin: the value that NaNs were replaced by in previous steps baselineMode: ‘min’ or ‘edge’, ‘edge’ will use nBaselineChannels/2 from each edge trimChannels: if integer, use that number of channels. If float between 0 and 1, then use that fraction of channels in each contiguous list (rounding up). If ‘auto’, then use 0.1 but not more than maxTrim channels. The ‘auto’ mode can also cause trimChannels to be set to 13 or 6 in certain circumstances (when moderately strong lines are present). maxTrim: if trimChannels=’auto’, this is the max channels to trim per group for TDM spws; it is automatically scaled upward for FDM spws. maxTrimFraction: in trimChannels=’auto’, the max fraction of channels to trim per group separator: the character to use to separate groups of channels in the string returned negativeThresholdFactor: scale the nominal negative threshold by this factor (to adjust sensitivity to absorption features: smaller values=more sensitive) dropBaselineChannels: percentage of extreme values to drop in baseline mode ‘min’ madRatioUpperLimit, madRatioLowerLimit: if ratio of MADs is between these values, then apply dropBaselineChannels when defining the MAD of the baseline range fitResult: coefficients from linear or quadratic fit, only relevant when meanSpectrumMethod is not ‘mom0mom8jointMask’ signalRatioTier1: threshold for signalRatio, above which we desensitize the level to consider line emission in order to avoid small differences in noise levels from affecting the result (e.g. that occur between serial and parallel tclean cubes) signalRatio=1 means: no lines seen, while closer to zero: more lines seen signalRatioTier2: second threshold for signalRatio, used for FDM spws (nchan>192) and for cases of peakOverMad < 5. Should be < signalRatioTier1. sigmaFindContinuumMode: if ‘auto’, then do not desensitize with signalRatioTier1.

findOuterAnnulusForPBCube(pbcube[, …])

Given a PB cube, finds the minimum sensitivity value, then computes the corresponding higher value to form an annulus.

findSpectralAxis(img)

This function is called by computeStatisticalSpectrumFromMask, getImageInfo, findContinuum and numberOfChannelsInCube.

findWidestContiguousListInChannelRange(…)

Given a list of channels, find the longest contiguous list in the first half or second half example: CASA <43>: fc.findWidestContiguousListInOtherHalf([2,3,7,8,9,12,13,14,15], True, 20) Out[43]: [12, 13, 14, 15] channels: list of blue points channel numbers channelRange: range of channels to look for widest set of blue points continuumChannels: current selection of continuuum channels

findWidestContiguousListInOtherHalf(…)

Given a list of channels, find the longest contiguous list in the first half or second half example: CASA <43>: fc.findWidestContiguousListInOtherHalf([2,3,7,8,9,12,13,14,15], True, 20) Out[43]: [12, 13, 14, 15] channels: list of blue points channel numbers lowerHalf: a Boolean (True for lower half, False for upper half) nchan: total channels in spectrum

flattenMask(mask[, outfile, overwrite])

Takes a multi-channel CASA image mask and propagates all pixels to a new single plane image mask outfile: default name = <mask>.flattened - Todd Hunter

frames([velocity, datestring, ra, dec, …])

This function is called by restToTopo.

gatherWarnings(selection, chanInfo, …)

Used in Pipeline2020 going forward.

gaussianBeamOffset([response, fwhm])

++++ This function is used by pipeline in computeStatisticalSpectrumFromMask().

gaussianBeamResponse(radius, fwhm)

++++ This function is used by pipeline in computeStatisticalSpectrumFromMask().

getDateObs(img[, myia])

This function is called by cubeLSRKToTopo.

getEquinox(img[, myia])

This function is called by cubeLSRKToTopo.

getFieldnameFromPipelineImageName(img[, verbose])

Extracts the field name from the pipeline image file name.

getFreqType(img)

This function is called by runFindContinuum and cubeLSRKToTopo.

getImageInfo(img[, returnBeamAreaPerChannel])

This function is called by findContinuum and meanSpectrum. Extract the beam and pixel information from a CASA image. This was copied from getFitsBeam in analysisUtils.py. Returns: A list of 12 things: bmaj, bmin, bpa, cdelt1, cdelt2, naxis1, naxis2, frequency, shape, crval1, crval2, maxBaseline Beam angles are in arcseconds (bpa in degrees), crvals are in radians Frequency is in GHz and is the central frequency MaxBaseline is in meters, inferred from freq and beamsize (will be zero for .residual images).

getMaxpixpos(img)

This function is called by findContinuum, but only in Cycle 4+5 heuristic.

getMemorySize()

This function is called by findContinuum and runFindContinuum.

getObservationStart(vis[, obsid, verbose])

Reads the start time of the observation from the OBSERVATION table (using tb tool) and reports it in MJD seconds.

getObservationStartDate(vis[, obsid, …])

Uses the tb tool to read the start time of the observation and reports the date.

getObservatoryName(vis)

+++++++ This function is not used by pipeline, because it is only used by getScienceSpws.

getScienceSpws(vis[, intent, returnString, …])

+++++++ This function is not used by pipeline, because it is only used by isSingleContinuum. Return a list of the each spw with the specified intent. For ALMA data, it ignores channel-averaged and SQLD spws. returnString: if True, it returns: ‘1,2,3’ if False, it returns: [1,2,3].

getSpwFromPipelineImageName(img[, verbose])

Extracts the spw ID from the pipeline image file name.

getTelescope(img[, myia])

This function is called by CalcAtmTransmissionForImage and cubeLSRKToTopo.

getcube(i[, filename])

++++ This function is not used by pipeline.

grep(filename, arg)

This function is called only by maskArgumentMismatch and centralArcsecArgumentMismatch.

help([match, debug])

Print an alphabetized list of all the defined functions at the top level in this python file.

imagePercentileNoMask(img, percentiles)

percentiles: a single value or a list Returns: a single value or a list

imageSNR(img[, axes, mask, maskWithAnnulus, …])

Uses imstat (which automatically applies the internal mask, i.e avoids the black edges of the image) to compute (max-median)/scaledMAD axes: passed to imstat (use [0,1,2] to get a vector SNR per channel) img: a 2D image or a cube mask: additional mask image (or expression) to apply when computing everything but max maskWithAnnulus: additional mask image (or expression) to use for median and MAD if useAnnulus==True (but always take the smaller of the 2 possible MADs) applyMaskToAll: also apply the provided mask when computing the max (default=False=anywhere) returnAllStats: if True, then also return the peak, median and scaledMAD where the peak has not had the median subtracted (in contrast to the SNR) returnFractionalDifference: if True, then also return the fraction of negative pixels usepb: (optional) name of a single-plane pb image to include in the mask (at >0.23 level) -ToddHunter

imageSNRAnnulus(mymomDiff, jointMaskTest, …)

invertChannelRanges(invertstring[, nchan, …])

Takes a CASA channel selection string and inverts it.

isSingleContinuum(vis[, spw, source, …])

+++++++ This function is not used by pipeline, rather it is an expected input parameter. Checks whether an spw was defined as single continuum in the OT by looking at the transition name in the SOURCE table of a measurement set. vis: a single measurement set, a comma-delimited list, or a python list (only the first one will be used) Optional parameters: spw: can be integer ID or string integer ID; if not specified, then it will use the first science spw source: passed to transition(); if not specified, it will use the first one intent: if source is blank then use first one with matching intent and spw mymsmd: an existing instance of the msmd tool.

is_binary(filename)

This function is called by runFindContinuum.

linfit(x, y, yerror[, pinit])

This function is called by atmosphereVariation in Cycle 4+5+6 and by runFindContinuum in Cycle 4+5.

locatePBCube(cube)

lsrkToRest(lsrkFrequency, velocityLSRK, …)

This function is called by lsrkToTopo.

lsrkToTopo(lsrkFrequency, datestring, ra, dec)

This function is called by cubeLSRKToTopo.

makeUvcontsub([files, fitorder, useFrequency])

+++++++ This function is not used anywhere else in this file. Takes a list of output files from findContinuum and builds an spw selection for uvcontsub, then prints the resulting commands to the screen. files: a list of files, either a comma-delimited string, a python list, or a wildcard string fitorder: passed through to the uvcontsub useFrequency: if True, then produce selection string in frequency; otherwise use channels Note: frequencies in the findContinuum .dat file are topo, which is what uvcontsub wants.

maskArgumentMismatch(mask, meanSpectrumFile, …)

This function is called by checkForMismatch.

maxLengthOfLists(lists)

This function is called by findContinuumChannels.

meanSpectrum(img[, nBaselineChannels, …])

This function is not used by Cycle 6 pipeline, but remains as manual user option. Computes a threshold and then uses it to compute the average spectrum across a CASA image cube, via the specified method. Inputs: nBaselineChannels: number of channels to use as the baseline in the ‘meanAboveThreshold’ methods. baselineMode: how to select the channels to use as the baseline: ‘edge’: use an equal number of channels on each edge of the spw ‘min’: use the percentile channels with the lowest absolute intensity sigmaCube: multiply this value by the rms to get the threshold above which a pixel is included in the mean spectrum nanBufferChannels: when removing or replacing NaNs, do this many extra channels beyond their actual extent useAbsoluteValue: passed to avgOverCube percentile: used with baselineMode=’min’ continuumThreshold: if specified, only use pixels above this intensity level meanSpectrumFile: name of ASCII file to produce to speed up future runs centralArcsec: default=-1 means whole field, or a floating point value in arcsec mask: a mask image (e.g. from clean); restrict calculations to pixels with mask=1 chanInfo: the list returned by numberOfChannelsInCube meanSpectrumMethod: ‘peakOverMad’, ‘peakOverRms’, ‘meanAboveThreshold’, ‘meanAboveThresholdOverRms’, or ‘meanAboveThresholdOverMad’ * ‘meanAboveThreshold’: uses a selection of baseline channels to compute the rms to be used as a threshold value (similar to constructing a moment map). * ‘meanAboveThresholdOverRms/Mad’: scale spectrum by RMS or MAD * ‘peakOverRms/Mad’ computes the ratio of the peak in a spatially-smoothed version of cube to the RMS or MAD. Smoothing is set by peakFilterFWHM. peakFilterFWHM: value used by ‘peakOverRms’ and ‘peakOverMad’ to presmooth each plane with a Gaussian kernel of this width (in pixels) Set to 1 to not do any smoothing. useIAGetProfile: if True, then for peakOverMad, or meanAboveThreshold with baselineMode=’min’, then use ia.getprofile instead of ia.getregion and the subsequent arithmetic (because it should be faster) useThresholdWithMask: if False, then just take mean rather than meanAboveThreshold when a mask has been specified Returns 8 items: * avgspectrum (vector) * avgspectrumAboveThresholdNansRemoved (vector) * avgspectrumAboveThresholdNansReplaced (vector) * threshold (scalar) * edgesUsed: 0=lower, 1=upper, 2=both * nchan (scalar) * nanmin (the value used to replace NaNs) * percentagePixelsNotMasked.

meanSpectrumFromMom0Mom8JointMask(cube, …)

This function is called by runFindContinuum when meanSpectrumMethod=’mom0mom8jointMask’. This is the new heuristic for Cycle 6 which creates the moment 0 and moment 8 images for a cube, takes their union and determines the mean spectrum by calling computeStatisticalSpectrumFromMask(), which uses ia.getprofile. pbcube: if not specified, then assume ‘.residual’ should be replaced in the name by ‘.pb’ overwriteMoments: rebuild the mom0 and mom8 images even if they already exist overwriteMasks: rebuild the mom0 and mom8 mask images even if they already exist phase2: if True, then run a second phase if SNR is high minPixelsInJointMask: if fewer than these pixels are found, then use all pixels above pb=0.3 userJointMask: if specified, use this joint mask instead of computing one; if it is a cube mask, as from tclean, then this function will form a flattened version first, and then use that snrThreshold: if SNR is higher than this, and phase2==True, then run a phase 2 mask calculation mom0minsnr: sets the threshold for the mom0 image (i.e. median + mom0minsnr*scaledMAD of the whole image) mom8minsnr: sets the threshold for the mom8 image (i.e. median + mom8minsnr*scaledMAD of thw whole image) rmStatContQuadratic: if True, then do not remove the old-style quadratic in this function in favor of the new-style removal elsewhere bidirectionalMaskPhase2: True=extend mask to negative values beyond threshold; False=Cycle6 behavior outdir: directory to write the mom0, mom8, masks, .dat and mean spectrum text files avoidExtremaInNoiseCalcForJointMask: experimental Boolean to avoid pixels <5%ile and >95%ile in the chauvenet imstat of mom0 and mom8 images momentdir: alternate directory to look for existing mom0 and mom8 images Returns: 13 things: * 1) the meanSpectrum * 2) a Boolean which states whether normalization was applied * 3) the number of pixels in the mask * 4) a Boolean for whether the mask reverted to pb-based * 5) a Boolean for whether a quadratic was removed * 6) the initialQuadraticImprovementRatio * 7) the list of 3 mom0snrs * 8) the list of 3 mom8snrs * 9) the number of regions pruned * 10) number of pixels in moment 8 mask * 11) mom0peak (Jy*km/s) * 12) mom8peak (Jy) * 13) name of jointMask file.

medianCorrected(baselineMode, percentile, …)

This function is called by findContinuumChannels. Computes the true median of a datastream from the observed median and MAD of the lowest Nth percentile points. signalRatio: when there is a lot of signal, we need to reduce the correction factor because it is less like Gaussian noise It is 1.0 if no lines are present, 0.25 if half the channels have signal, etc.

mjdSecondsToMJDandUT(mjdsec[, debug, prec, …])

This function is called by mjdToUT.

mjdToUT(mjd[, use_metool, prec])

This function is called by getDateObs. Converts an MJD value to a UT date and time string such as ‘2012-03-14 00:00:00 UT’ use_metool: whether or not to use the CASA measures tool if running from CASA. This parameter is simply for testing the non-casa calculation. -Todd Hunter.

nanmean(a[, axis])

This function is called by findContinuumChannels, runFindContinuum, and avgOverCube.

nanmedian(x[, axis, preop])

This function is called by MAD, avgOverCube, and meanSpectrum.

numberOfChannelsInCube(img[, returnFreqs, …])

This function is called by findContinuum, cubeLSRKToTopo, computeStatisticalSpectrumFromMask, and meanSpectrum. Finds the number of channels in a CASA image cube. returnFreqs: if True, then also return the frequency of the first and last channel (in Hz) returnChannelWidth: if True, then also return the channel width (in Hz) verbose: if True, then print the frequencies of first and last channel -Todd Hunter.

oneEvent(npts[, events, positive, verbose])

This function is called by meanSpectrumFromMom0Mom8JointMask.

onlyExtraMaskYesOrNo(badAtmosphere, median, …)

Used in Pipeline2020 going forward. badAtmosphere: either a boolean or a string (‘goodAtm’ or ‘badAtm’) momDiffSNR: SNR where peak is measured over whole image (instead of outside mask), (mom scaledMAD is in denominator) momDiffSNRCube: SNR where peak is measured over whole image (instead of outside mask), (cube scaledMAD is in denominator) NpixCubeAnywhere: unused now.

parseFrequencyArgument(bandwidth)

This function is called by parseFrequencyArgumentToGHz, topoFreqToChannel and cubeLSRKToTopo.

parseFrequencyArgumentToGHz(bandwidth)

This function is called by frames and lsrkToRest.

pickAutoTrimChannels(totalChannels, length, …)

This function is called by splitListIntoContiguousListsAndTrim and pickTrimString.

pickNarrow(length)

This function is called by pickNarrowString and findContinuumChannels. Automatically picks a setting for the narrow parameter of findContinuumChannels() based on the number of channels in the spectrum. Returns: an integer Examples: This formula results in the following values: length: 64,128,240,480,960,1920,3840,7680: return: 2, 3, 3, 3, 3, 4, 4, 4 (ceil(log10)) **** current function **** 2, 2, 2, 3, 3, 3, 4, 4 (round_half_up(log10)) 1, 2, 2, 2, 2, 3, 3, 3 (floor(log10)) 5, 5, 6, 7, 7, 8, 9, 9 (ceil(log)) 4, 5, 5, 6, 7, 8, 8, 9 (round_half_up(log)) 4, 4, 5, 6, 6, 7, 8, 8 (floor(log)).

pickNarrowString(narrow, length[, …])

This function is called by runFindContinuum.

pickRandomError([seed])

Picks a random value from a Gaussian distribution with mean 0 and standard deviation = 1.

pickRandomErrors([nvalues])

Picks a series of random values from a Gaussian distribution with mean 0 and standard deviation = 1 and returns it as an array.

pickTrimString(totalChannels, trimChannels, …)

This function is called by runFindContinuum.

pick_sFC_TDM(meanSpectrumMethod, singleContinuum)

This function is called by runFindContinuum.

plotChannelSelections(ax1, selection, …[, …])

avgspectrumAboveThreshold: computed in runFindContinuum medianTrue, threshold: computed by findContinuumChannels Returns labelDescs: list of plot artists that may be removed and replaced later

plotMeanSpectrum(meanSpectrumFile[, …])

Reads a *.findcont.residual.meanSpectrum.<method> file and plots it as intensity vs. channel. selection: channel selection string to draw in cyan horizontal lines.

plotStatisticalSpectrumFromMask(cube[, …])

Simple wrapper to call computeStatisticalSpectrumFromMask and plot it. jointMask: either a mask image, or an expression with < or > on an mask image.

polyfit(x, y, yerror[, pinit])

This function is called by removeInitialQuadraticIfNeeded in Cycle 6, and by runFindContinuum in Cycle 4+5. Basic second-order polynomial function fitter with error bars in y-axis data points. Uses scipy.optimize.leastsq(). Accepts either lists or arrays. Input: x, y: x and y axis values yerror: uncertainty in the y-axis values (vector or scalar) pinit contains the initial guess of [slope, intercept] y = order2coeff*(x-xoffset)**2 + slope*(x-xoffset) + y-intercept Output: The fit result as: [xoffset, order2coeff, slope, y-intercept].

propagateMaskToAllChannels(mask[, axis])

This function is called by meanSpectrum. Takes a spectral mask array, and propagates every masked spatial pixel to all spectral channels of the array. Returns: a 2D mask (of the same spatial dimension as the input, but with only 1 channel).

pruneMask(mymask[, psfcube, minbeamfrac, …])

available in CASA >= 5.4.0-X (see CAS-11335) Removes any regions smaller than prunesize contiguous pixels. psfcube: if specified, then use minbeamfrac * pixelsPerBeam, otherwise use prunesize prunesize: value used if psfcube is None; use 6 because pipeline size mitigation could produce images that have only 3 pixels across the beam, meaning only ~7-8 pts/beam nchan: number of channels to process, will always be 1 in findContinuum’s use case overwrite: if True, replaces mymask with pruned mask, putting original in mymask.prepruned if False, is simply creates: mymask.pruned, but only if it pruned anything If all regions are pruned, it checks if it is ACA7m (maxBaseline < 60m) and if so, tries again with minbeamfrac*0.5.

rad2radec([ra, dec, prec, verbose, …])

This function is called by cubeLSRKToTopo.

readContDat(filename)

Reads the channel selection for a _findContinuum.dat file

readContDatAggregateContinuum(filename)

Reads the aggregate continuum from a _findContinuum.dat file.

readMeanSpectrumFITSFile(meanSpectrumFile[, …])

++++++ This function is not used by the pipeline.

readPreviousMeanSpectrum(meanSpectrumFile[, …])

++++++ This function is not used by the pipeline. Read a previous calculated mean spectrum and its parameters, or an ASCII file created by the CASA viewer (or tt.ispec). Note: only the intensity column(s) are returned, not the channel/frequency columns. This function will not typically be used by the pipeline, only manual users. Returns: 11 things: avgspectrum, avgSpectrumNansReplaced, threshold, edgesUsed, nchan, nanmin, firstFreq, lastFreq, centralArcsec, mask, percentagePixelsNotMasked.

readViewerOutputFile(lines[, debug])

++++++ This function is not used by the pipeline.

rejectNarrowInnerWindowsChannels(channels[, …])

This function is called by findContinuumChannels.

removeIfNecessary(imgname)

removeInitialQuadraticIfNeeded(avgSpectrum)

This function is called by runFindContinuum when meanSpectrumMethod=’mom0mom8jointMask’.

removeNaNs(a[, replaceWithMin, verbose, …])

+++++++ This function is not used for meanSpectrumMethod == ‘mom0mom8jointMask’ by pipeline,

removeStatContQuadratic(y[, nsigma, …])

Iteratively removes outlier points down to nsigma, then fits a quadratic to the remaining points, removes this quadratic from the full initial spectrum, and returns it, along with the list of channels used for the fit and the fit evaluated at all channels.

removeZeros(a)

+++++ This function is not used for meanSpectrumMethod=’mom0mom8jointMask’.

replaceLineFullRangesWithNoise(intensity, …)

Used in Pipeline2020 going forward.

restToTopo(restFrequency, velocityLSRK, …)

This function is called by lsrkToTopo.

robustMADofContinuumRanges(myChannelList, …)

Used in Pipeline2020 going forward.

roundFigures(value, digits)

This function is called by runFindContinuum and drawYlabel.

round_half_up(x)

runFindContinuum([img, pbcube, psfcube, …])

This function is called by findContinuum. It calls functions that: 1) compute the mean spectrum of a dirty cube 2) find the continuum channels 3) plot the results Inputs: channelWidth: in Hz Returns: 23 items *1 A channel selection string suitable for the spw parameter of clean. *2 The name of the png produced. *3 The slope of the linear fit (or None if no fit attempted). *4 The channel width in Hz (only necessary for the fitsTable option). *5 nchan in the cube (only necessary to return this for the fitsTable option, will be computed otherwise). *6 Boolean: True if it used low values as the baseline (as opposed to high values) *7 SNRs in the moment0 image (raw, outside mask, outside phase2 mask) *8 SNRs in the moment8 image (raw, outside mask, outside phase2 mask) *9 Boolean: True if the middle-valued channels were used as the baseline (as opposed to low or high) *10 list: selectionPreBluePruning *11 float: finalSigmaFindContinuum *12 string: name of jointMask *13 float array: avgspectrumAboveThreshold *14 float: the true median of the spectrum *15 list: of plot artist descriptors (things to potentially remove and replot) *16 axis instance: ax1 (lower x-axis) *17 axis instance: ax2 (upper x-axis) *18 float: threshold (the positive threshold for line detection) *19 string: final line of upper text legend *20 int: number of narrow central groups that were dropped (or would have been dropped if enableRejectNarrowWindows had been True, but wasn’t), summed over all runs of findContinuumChannels() *21 float: effectiveSigma (product of finalSigmaFindContinuum*correctionFactor) *22 float: baselineMAD *23 string: upper x-axis label Inputs: img: the image cube to operate upon spw: the spw name or number to put in the x-axis label transition: the name of the spectral transition (for the plot title) baselineModeA: ‘min’ or ‘edge’, method to define the baseline in meanSpectrum() sigmaCube: multiply this value by the MAD to get the threshold above which a pixel is included in the mean spectrum nBaselineChannels: if integer, then the number of channels to use in: a) computing the mean spectrum with the ‘meanAboveThreshold’ methods. b) computing the MAD of the lowest/highest channels in findContinuumChannels if float, then the fraction of channels to use (i.e. the percentile) default = 0.19, which is 24 channels (i.e. 12 on each side) of a TDM window sigmaFindContinuum: passed to findContinuumChannels, ‘auto’ starts with 3.5 sigmaFindContinuumMode: ‘auto’, ‘autolower’, or ‘fixed’ verbose: if True, then print additional information during processing png: the name of the png to produce (‘’ yields default name) pngBasename: if True, then remove the directory from img name before generating png name nanBufferChannels: when removing or replacing NaNs, do this many extra channels beyond their extent source: the name of the source, to be shown in the title of the spectrum. if None, then use the filename, up to the first underscore. findContinuum: if True, then find the continuum channels before plotting overwrite: if True, or ASCII file does not exist, then recalculate the mean spectrum writing it to <img>.meanSpectrum if False, then read the mean spectrum from the existing ASCII file trimChannels: after doing best job of finding continuum, remove this many channels from each edge of each block of channels found (for margin of safety) If it is a float between 0..1, then trim this fraction of channels in each group (rounding up). If it is ‘auto’, use 0.1 but not more than maxTrim channels and not more than maxTrimFraction percentile: control parameter used with baselineMode=’min’ continuumThreshold: if specified, only use pixels above this intensity level narrow: the minimum number of channels that a group of channels must have to survive if 0<narrow<1, then it is interpreted as the fraction of all channels within identified blocks if ‘auto’, then use int(ceil(log10(nchan))) titleText: default is img name and transition and the control parameter values meanSpectrumFile: an alternative ASCII text file to use for the mean spectrum. Must also set img=’’. meanSpectrumMethod: ‘peakOverMad’, ‘peakOverRms’, ‘meanAboveThreshold’, ‘meanAboveThresholdOverRms’, ‘meanAboveThresholdOverMad’, where ‘peak’ refers to the peak in a spatially-smoothed version of cube peakFilterFWHM: value used by ‘peakOverRms’ and ‘peakOverMad’ to presmooth each plane with a Gaussian kernel of this width (in pixels) Set to 1 to not do any smoothing. centralArcsec: radius of central box within which to compute the mean spectrum default=’auto’ means start with whole field, then reduce to 1/10 if only one window is found (unless mask is specified); or ‘fwhm’ for the central square with the same radius as the PB FWHM; or a floating point value in arcseconds mask: a mask image equal in shape to the parent image that is used to determine the region to compute the noise (outside the mask) and the mean spectrum (inside the mask) option ‘auto’ will look for the <filename>.mask that matches the <filename>.image and if it finds it, it will use it; otherwise it will use no mask plotAtmosphere: ‘’, ‘tsky’, or ‘transmission’ airmass: for plot of atmospheric transmission pwv: in mm (for plot of atmospheric transmission) channelFractionForSlopeRemoval: if at least this many channels are initially selected, or if there are only 1-2 ranges found and the widest range has > nchan*0.3 channels, then fit and remove a linear slope and re-identify continuum channels. Set to 1 to turn off. quadraticFit: if True, fit quadratic polynomial to the noise regions when appropriate; otherwise fit only a linear slope megapixels: simply used to label the plot triangularPatternSeen: if True, then slightly boost sigmaFindContinuum if it is in ‘auto’ mode maxMemory: only used to name the png if it is set applyMaskToMask: if True, apply the mask inside the user mask image to set its masked pixels to 0 plotBaselinePoints: if True, then plot the baseline-defining points as black dots dropBaselineChannels: percentage of extreme values to drop in baseline mode ‘min’ useIAGetProfile: if True, then for meanAboveThreshold and baselineMode=’min’, then use ia.getprofile instead of ia.getregion and subsequent arithmetic (faster) initialQuadraticImprovementThreshold: if removal of a quadratic from the raw meanSpectrum reduces the MAD by this factor or more, then proceed with removing this quadratic (new Cycle 6 heuristic) maxMadRatioForSFCAdjustment: madRatio must be below this value (among other things) in order for automatic lowering of sigmaFindContinuum value to occur avoidance: a CASA channel selection string to avoid selecting (for manual use).

setYLimitsAvoidingEdgeChannels(…[, chan])

Called by runFindContinuum to set the y-axis limits.. 1) Avoid spikes in edge channels from skewing the plot limits, by setting plot limits on the basis of channels chan..-chan 2) Enforce a maximum dynamic range (peak-median)/mad so that weak features can be seen, allowing strongest ones to overflow the top.

sigmaCorrectionFactor(baselineMode, npts, …)

This function is called by findContinuumChannels. Computes the correction factor for the fact that the measured rms (or MAD) on the N%ile lowest points of a datastream will be less than the true rms (or MAD) of all points. Returns: a value between 0 and 1, which will be used to reduce the estimate of the population sigma based on the sample sigma.

splitListIntoContiguousLists(mylist)

This function is called by findContinuumChannels.

splitListIntoContiguousListsAndRejectNarrow(…)

This function is called by findContinuumChannels. Split a list of channels into contiguous lists, and reject those that have a small number of channels. narrow: if >=1, then interpret it as the minimum number of channels that a group must have in order to survive if <1, then interpret it as the minimum fraction of channels that a group must have relative to the total number of channels Returns: a new single list (as an array) Example: [1,2,3,4,6,7,9,10,11] –> [ 1, 2, 3, 4, 9, 10, 11].

splitListIntoContiguousListsAndRejectZeroStd(…)

This function is called by findContinuumChannels.

splitListIntoContiguousListsAndTrim(…[, …])

This function is called by findContinuumChannels. Split a list of channels into contiguous lists, and trim some number of channels from each edge. totalChannels: number of channels in the spectrum channels: a list of channels selected for potential continuum trimChannels: if integer, use that number of channels. If float between 0 and 1, then use that fraction of channels in each contiguous list (rounding up). If ‘auto’, then use 0.1 but not more than maxTrim channels and not more than maxTrimFraction of channels. maxTrim: used in ‘auto’ mode Returns: a new single list.

splitListIntoHomogeneousLists(mylist)

This function is called only by removeNaNs, and hence is not used in Cycle 6-onward pipeline.

submask(mask, region)

This function is called by meanSpectrum, but only in Cycle 4+5 heuristic.

tdmSpectrum(channelWidth, nchan)

This function is called by runFindContinuum and findContinuum.

tooLittleBandwidth(selection, chanInfo[, …])

Used in Pipeline2020 going forward.

tooLittleSpread(selection, chanInfo[, fraction])

Used in Pipeline2020 going forward.

topoFreqRangeListToChannel([contdatline, …])

++++ This function is used by pipeline in writeContDat() Converts a semicolon-delimited string list of topocentric frequency ranges to channel ranges in the specified ms.

topoFreqToChannel(freqlist, vis, spw[, mymsmd])

++++ This function is used by pipeline in writeContDat() Converts a python floating point list of topocentric frequency ranges to channel ranges in the specified ms.

transition(vis, spw[, source, intent, …])

+++++++ This function is not used by pipeline, because it is only used by isSingleContinuum.

updateChannelRangesOnPlot(labelDescs, …[, …])

Used in Pipeline2020 going forward.

version([showfile])

Returns the CVS revision number.

versionStringToArray(versionString)

This function is called by casaVersionCompare.

widthOfMaskArcsec(mask, maskInfo)

++++ This function is not used by pipeline when meanSpectrumMethod=’mom0mom8jointMask’ or if mask=’’ Finds width of the smallest central square that circumscribes all masked pixels.

writeContDat(meanSpectrumFile, selection, …)

This function is called by findContinuum. Writes out an ASCII file called <meanSpectrumFile>_findContinuum.dat that contains the selected channels, the png name and the aggregate bandwidth in GHz. Returns: None meanSpectrumFile: only used to generate the name of the .dat file, unless firstFreq <= 0, in which case the readPreviousMeanSpectrum is called on it (It splits off name before ‘.meanSpectrum’ and appends _findContinuum.dat) vis: if specified, then also write a line with the topocentric velocity ranges spw: integer or string ID number; if specified (along with vis), then also write a final line with the topocentric channel ranges for this spw.

writeMeanSpectrum(meanSpectrumFile, …[, …])

This function is called by meanSpectrum.