from casatasks.private.casa_transition import is_CASA6
from casatasks import casalog
from casatools import quanta
from casatools import image
from taskinit import casalog
from taskinit import qatool as quanta
from taskinit import iatool as image
def plotprofilemap(imagename=None, figfile=None, overwrite=None, transparent=None,
pol=None, spectralaxis=None, restfreq=None, plotrange=None, title=None,
linecolor=None, linestyle=None, linewidth=None,
separatepanel=None, plotmasked=None, maskedcolor=None,
showaxislabel=None, showtick=None, showticklabel=None,
figsize=None, numpanels=None):
casalog.origin('plotprofilemap')
if len(figfile) > 0 and os.path.exists(figfile) and overwrite == False:
raise RuntimeError('overwrite is False and output file exists: \'%s\''%(figfile))
image = SpectralImage(imagename)
if pol < 0 or pol > npol - 1:
raise RuntimeError('pol {pol} is out of range (Stokes axis {stokes})'.format(pol=pol,stokes=image.stokes))
parsed_size = parse_figsize(figsize)
nx, ny = parse_numpanels(numpanels)
if (not isinstance(restfreq, str)) or len(restfreq) == 0:
plot_profile_map(image, figfile, pol, spectralaxis, restfreq, title, linecolor, linestyle, linewidth,
separatepanel, plotmasked, maskedcolor,
showaxislabel, showtick, showticklabel, parsed_size,
nx, ny, transparent, plotrange)
casalog.post('Error: %s' % (str(e)), priority='SEVERE')
casalog.post(traceback.format_exc(), priority='DEBUG')
NoDataThreshold = NoData + 10000.0
LightSpeedQuantity = qa.constants('c')
LightSpeed = qa.convert(LightSpeedQuantity, 'km/s')['value'] # speed of light in km/s
def Deg2HMS(x, arrowance):
# Transform degree to HHMMSS.sss format
s = ((xx % 15) * 4 - m) * 60.0
(h, m, s) = Deg2HMS(x, 1 / 240000.0)
#return '%02dh%02dm%05.2fs' % (h, m, s)
return '%02d%s%02d%s%05.2f' % (h, hsyb, m, msyb, s)
def Deg2DMS(x, arrowance):
# Transform degree to +ddmmss.ss format
xxx = (x + 90) % 180 - 90
xx = abs(xxx) + arrowance
s = ((xx % 1) * 60 - m) * 60.0
(d, m, s) = Deg2DMS(x, 1 / 360000.0)
#return '%+02dd%02dm%04.1fs' % (d, m, s)
sstr = ('%3.1f'% ( s - int(s))).lstrip('0')
return '%+02d%s%02d\'%02d\"%s' % (d, dsyb, m, sint, sstr)
def parse_figsize(figsize):
casalog.post('parse_figsize input: {input}'.format(input=figsize), priority='DEBUG')
if figsize is not None and isinstance(figsize, str) and len(figsize) > 0:
size_list = figsize.split(',')
size_inch_list = [x / 25.4 for s in size_list for x in qa.getvalue(qa.convert(qa.quantity(s),outunit='mm'))]
if len(size_inch_list) == 1:
parsed = tuple(size_inch_list[:2])
casalog.post('parse_figsize output: {output}'.format(output=parsed), priority='DEBUG')
def parse_numpanels(numpanels):
if isinstance(numpanels, str) and len(numpanels) > 0:
n = map(int, numpanels.split(','))
casalog.post('parse_numpanels output: {output}'.format(output=parsed), priority='DEBUG')
class ProfileMapAxesManager(object):
label_map = {'Right Ascension': 'RA',
def __init__(self, nh, nv, brightnessunit, direction_label, direction_reference,
spectral_label, spectral_unit, ticksize, title='',
showaxislabel=False, showtick=False, showticklabel=False,
self.brightnessunit = brightnessunit
self.numeric_formatter = pl.FormatStrFormatter('%.2f')
self.direction_label = direction_label
self.direction_reference = direction_reference
self.separatepanel = separatepanel
self.spectral_label = spectral_label
self.spectral_unit = spectral_unit
self.showaxislabel = showaxislabel
self.showticklabel = showticklabel
casalog.post('figsize={figsize}'.format(figsize=self.figsize), priority='DEBUG')
self.normalization_factor = 1
# to resize matplotlib window to specified size
pl.figure(self.MATPLOTLIB_FIGURE_ID)
pl.figure(self.MATPLOTLIB_FIGURE_ID)
pl.figure(self.MATPLOTLIB_FIGURE_ID, figsize=self.figsize)
def MATPLOTLIB_FIGURE_ID(self):
if self._axes_spmap is None:
self._axes_spmap = list(self.__axes_spmap())
return 0.01 + 0.2 / self.ncolumn
return 0.01 + 0.5 / self.nrow
def horizontal_space(self):
return self.horizontal_subplot_size * 0.1
def vertical_space(self):
return self.vertical_subplot_size * 0.1
if self.showaxislabel or (self.showtick and self.showticklabel):
if self.showaxislabel or (self.showtick and self.showticklabel):
if isinstance(self.title, str) and len(self.title) > 0:
return 0.04 * (self.title.count('\n') + 1)
def horizontal_subplot_size(self):
return (1.0 - self.left_margin - self.right_margin - self.ylabel_area) / self.ncolumn
def vertical_subplot_size(self):
return (1.0 - self.bottom_margin - self.top_margin - self.xlabel_area - self.title_area) / self.nrow
def set_normalization_factor(self, factor):
self.normalization_factor = factor
w = self.horizontal_subplot_size
h = self.vertical_subplot_size
l = 1.0 - self.right_margin - w * (x + 1) + 0.5 * self.horizontal_space
b = self.bottom_margin + self.ylabel_area + h * y + 0.5 * self.vertical_space
axes = pl.axes([l, b, w - self.horizontal_space, h - self.vertical_space])
if self.showaxislabel and y == 0 and x == self.nh - 1:
casalog.post('label "{label}" unit "{unit}"'.format(label=self.spectral_label, unit=self.spectral_unit), priority='DEBUG')
if self.spectral_unit is not None and len(self.spectral_unit) > 0:
spectral_label = '%s [%s]' % (self.spectral_label, self.spectral_unit)
spectral_label = self.spectral_label
axes.xaxis.set_label_text(spectral_label,
if self.normalization_factor < 100 and self.normalization_factor > 0.01:
label_text = 'Intensity [%s]' % self.brightnessunit
label_text = 'Intensity [1e%d x %s]' % (int(numpy.log10(self.normalization_factor)),
axes.yaxis.set_label_text(label_text,
size=self.ticksize, rotation='vertical')
xlocator = pl.AutoLocator()
xlocator.set_params(nbins=6, prune='upper')
axes.xaxis.set_major_locator(xlocator)
ylocator = pl.AutoLocator()
ylocator.set_params(nbins=4)
axes.yaxis.set_major_locator(ylocator)
xformatter = pl.ScalarFormatter(useOffset=False)
axes.xaxis.set_major_formatter(xformatter)
axes.xaxis.set_tick_params(labelsize=max(self.ticksize - 1, 1))
axes.yaxis.set_tick_params(labelsize=max(self.ticksize - 1, 1))
if y != 0 or x != self.nh - 1:
axes.xaxis.set_major_formatter(pl.NullFormatter())
axes.yaxis.set_major_formatter(pl.NullFormatter())
axes.xaxis.set_major_formatter(pl.NullFormatter())
axes.yaxis.set_major_formatter(pl.NullFormatter())
axes.yaxis.set_major_locator(pl.NullLocator())
axes.xaxis.set_major_locator(pl.NullLocator())
def setup_labels(self, label_ra, label_dec):
w = self.horizontal_subplot_size
l = 1.0 - self.right_margin - w * (x + 1)
h = self.bottom_margin * 0.5
b = self.bottom_margin - h
a1 = pl.axes([l, b, w, h])
pl.text(0.5, 0.2, HHMMSSss((label_ra[x][0] + label_ra[x][1]) / 2.0, 0),
horizontalalignment='center', verticalalignment='center', size=self.ticksize)
a1.texts[0].set_text(HHMMSSss((label_ra[x][0] + label_ra[x][1]) / 2.0, 0))
w = self.horizontal_subplot_size
h = self.vertical_subplot_size
b = self.bottom_margin + y * h
a1 = pl.axes([l, b, w, h])
pl.text(0.5, 0.5, DDMMSSs((label_dec[y][0] + label_dec[y][1]) / 2.0, 0),
horizontalalignment='center', verticalalignment='center', size=self.ticksize)
a1.texts[0].set_text(DDMMSSs((label_dec[y][0] + label_dec[y][1]) / 2.0, 0))
l = self.left_margin + self.xlabel_area
h = self.bottom_margin * 0.5
w = 1.0 - l - self.right_margin
a1 = pl.axes([l, b, w, h])
xpos = (1.0 + 0.5 * self.nh) / self.ncolumn
casalog.post('xpos=%s' % (xpos), priority='DEBUG')
pl.text(xpos, 0.5, '%s (%s)' % (self.direction_label[0], self.direction_reference),
horizontalalignment='center', verticalalignment='center',
size=(self.ticksize + 2))
h = self.vertical_subplot_size
b = self.bottom_margin + 0.5 * (h * self.nrow - self.vertical_subplot_size)
a1 = pl.axes([l, b, w, h])
pl.text(1.0, 0.5, '%s (%s)' % (self.direction_label[1], self.direction_reference),
horizontalalignment='right', verticalalignment='center',
rotation='vertical', size=(self.ticksize + 2))
left = self.left_margin + self.xlabel_area
bottom = 1.0 - self.title_area - self.top_margin
width = 1.0 - left - self.right_margin
a1 = pl.axes([left, bottom, width, height])
xpos = (1.0 + 0.5 * self.nh) / self.ncolumn
pl.text(xpos, 0.1, self.title,
horizontalalignment='center', verticalalignment='bottom',
def plot_profile_map(image, figfile, pol, spectralaxis='', restfreq=None, title=None,
linecolor='b', linestyle='-', linewidth=0.2,
separatepanel=True, plotmasked=None, maskedcolor=None,
showaxislabel=False, showtick=False, showticklabel=False,
figsize=None, nx=-1, ny=-1, transparent=False, user_xlim=None):
plotmasked -- 'none': show neither lines nor panels
'empty': show empty panel
'text': show 'NO DATA' symbol
'plot': plot masked data with different color
if separatepanel is None:
maskedcolor = 'gray' if linecolor != 'gray' else 'black'
if showaxislabel is None:
if showticklabel is None:
# user-specified range for horizontal (spectral) axis
if isinstance(user_xlim, str) and user_xlim.find('~') != -1:
user_range = user_xlim.split('~')
if len(user_range[0]) > 0:
user_xmin = float(user_range[0])
if len(user_range[1]) > 0:
user_xmax = float(user_range[1])
num_panel = min(max(x_max - x_min + 1, y_max - y_min + 1), MaxPanel)
STEP = int((max(x_max - x_min + 1, y_max - y_min + 1) - 1) / num_panel) + 1
NH = (x_max - x_min) / STEP + 1
NV = (y_max - y_min) / STEP + 1
xSTEP = image.nx / NH + (1 if image.nx % NH > 0 else 0)
ySTEP = image.ny / NV + (1 if image.ny % NV > 0 else 0)
casalog.post('num_panel=%s, xSTEP=%s, ySTEP=%s, NH=%s, NV=%s' % (num_panel, xSTEP, ySTEP, NH, NV))
direction_label = image.direction_label
direction_reference = image.direction_reference
default_spectral_label = image.spectral_label
default_spectral_unit = image.spectral_unit
if default_spectral_label == 'Frequency':
default_spectral_unit = 'GHz'
default_spectral_data = image.spectral_data[chan0:chan1]
if spectralaxis is None or spectralaxis == '' or spectralaxis == default_spectral_label.lower():
spectral_label = default_spectral_label
spectral_unit = default_spectral_unit
spectral_data = default_spectral_data
elif spectralaxis == 'channel':
spectral_label = 'Channel'
spectral_data = numpy.arange(chan0, chan1, dtype=numpy.int32)
spectral_label = spectralaxis.capitalize()
if spectralaxis == 'frequency':
spectral_data = numpy.fromiter((image.to_frequency(v, freq_unit='GHz')
for v in default_spectral_data), dtype=numpy.float64)
elif spectralaxis == 'velocity':
casalog.post('User-specified rest frequency %s' % (restfreq))
casalog.post('Default rest frequency from input image %s' % (image.coordsys.restfrequency()))
spectral_data = numpy.fromiter((image.to_velocity(f, freq_unit='GHz', restfreq=restfreq)
for f in default_spectral_data), dtype=numpy.float64)
plotter = SDProfileMapPlotter(NH, NV, xSTEP, ySTEP, image.brightnessunit,
direction_label, direction_reference,
spectral_label, spectral_unit,
separatepanel=separatepanel,
showaxislabel=showaxislabel,
showticklabel=showticklabel,
masked_data = image.data * image.mask
masked_data[numpy.logical_not(numpy.isfinite(masked_data))] = 0.0
refpix[0], refval[0], increment[0] = image.direction_axis(0, unit='deg')
refpix[1], refval[1], increment[1] = image.direction_axis(1, unit='deg')
plotter.setup_labels(refpix, refval, increment)
stokes = image.stokes[pol]
casalog.post('Generate profile map for pol {stokes}'.format(stokes=stokes))
casalog.post('masked_data.shape=%s id_stokes=%s' % (list(masked_data.shape), image.id_stokes), priority='DEBUG')
masked_data_p = masked_data.take([pol], axis=image.id_stokes).squeeze(axis=image.id_stokes)
Plot = numpy.zeros((NH, NV, (chan1 - chan0)), numpy.float32) + NoData
mask_p = image.mask.take([pol], axis=image.id_stokes).squeeze(axis=image.id_stokes)
isvalid = numpy.any(mask_p, axis=2)
Nsp = sum(isvalid.flatten())
casalog.post('Nsp=%s' % (Nsp))
valid_index = isvalid[x0:x1, y0:y1].nonzero()
chunk = masked_data_p[x0:x1, y0:y1]
valid_sp = chunk[valid_index[0], valid_index[1], :]
Plot[x][y] = valid_sp.mean(axis=0)
plot_mask = numpy.logical_and(numpy.isfinite(Plot), Plot != NoData)
max_data = numpy.abs(Plot[plot_mask]).max()
casalog.post('max_data = %s' % (max_data), priority='DEBUG')
normalization_factor = numpy.power(10.0, int(numpy.log10(max_data)))
if normalization_factor < 1.0:
normalization_factor /= 10.
casalog.post('normalization_factor = %s' % (normalization_factor), priority='DEBUG')
plotter.set_normalization_factor(normalization_factor)
Plot[plot_mask] /= normalization_factor
status = plotter.plot(figfile, Plot,
class SDProfileMapPlotter(object):
def __init__(self, nh, nv, xstep, ystep, brightnessunit, direction_label, direction_reference,
spectral_label, spectral_unit, title=None, separatepanel=True,
showaxislabel=False, showtick=False, showticklabel=False,
if self.xstep > 1 or self.ystep > 1:
step = max(self.xstep, self.ystep)
ticksize = 10 - int(max(nh, nv) * step / (step - 1)) / 2
ticksize = 10 - int(max(nh, nv)) / 2
self.axes = ProfileMapAxesManager(nh, nv, brightnessunit,
direction_label, direction_reference,
spectral_label, spectral_unit,
separatepanel=separatepanel,
showaxislabel=showaxislabel,
showticklabel=showticklabel,
self.lines_averaged = None
self.reference_level = None
self.global_scaling = True
self.deviation_mask = None
return self.axes.ticksize
def setup_labels(self, refpix_list, refval_list, increment_list):
LabelRA = numpy.zeros((self.nh, 2), numpy.float32) + NoData
LabelDEC = numpy.zeros((self.nv, 2), numpy.float32) + NoData
increment = increment_list[0]
#casalog.post('axis 0: refpix,refval,increment=%s,%s,%s'%(refpix,refval,increment))
x0 = (self.nh - x - 1) * self.xstep
x1 = (self.nh - x - 2) * self.xstep + 1
LabelRA[x][0] = refval + (x0 - refpix) * increment
LabelRA[x][1] = refval + (x1 - refpix) * increment
increment = increment_list[1]
#casalog.post('axis 1: refpix,refval,increment=%s,%s,%s'%(refpix,refval,increment))
y1 = (y + 1) * self.ystep - 1
LabelDEC[y][0] = refval + (y0 - refpix) * increment
LabelDEC[y][1] = refval + (y1 - refpix) * increment
self.axes.setup_labels(LabelRA, LabelDEC)
def setup_reference_level(self, level=0.0):
self.reference_level = level
def set_global_scaling(self):
self.global_scaling = True
def unset_global_scaling(self):
self.global_scaling = False
def set_normalization_factor(self, factor):
self.axes.set_normalization_factor(factor)
def plot(self, figfile, map_data, frequency,
linecolor='b', linestyle='-', linewidth=0.2,
plotmasked='none', maskedcolor='gray', transparent=False,
user_xmin=None, user_xmax=None):
casalog.post('Number of spectral channels is 1. Resulting profile map will be useless.',
global_xmin = frequency[0] - 1.
global_xmax = frequency[0] + 1.
global_xmin = min(frequency[0], frequency[-1])
global_xmax = max(frequency[0], frequency[-1])
casalog.post('full spectral range: global_xmin=%s, global_xmax=%s' % (global_xmin, global_xmax))
if user_xmin is not None:
casalog.post('user-specified global_xmin %s' % (user_xmin))
if user_xmax is not None:
casalog.post('user-specified global_xmax %s' % (user_xmax))
if global_xmin == global_xmax:
raise RuntimeError('Wrong spectral axis range: no range to plot.')
elif global_xmin > global_xmax:
raise RuntimeError('Wrong spectral axis range: minimum > maximum')
# to eliminate max/min value due to bad pixel or bad fitting,
# 1/10-th value from max and min are used instead
valid_index = numpy.where(map_data.min(axis=2) > NoDataThreshold)
valid_data = map_data[valid_index[0], valid_index[1], :]
ListMax = valid_data.max(axis=1)
ListMin = valid_data.min(axis=1)
casalog.post('ListMax=%s' % (list(ListMax)))
casalog.post('ListMin=%s' % (list(ListMin)))
global_ymax = numpy.sort(ListMax)[len(ListMax) - len(ListMax) / 10 - 1]
global_ymin = numpy.sort(ListMin)[len(ListMin) / 10]
global_ymax = global_ymax + (global_ymax - global_ymin) * 0.2
global_ymin = global_ymin - (global_ymax - global_ymin) * 0.1
casalog.post('global_ymin=%s, global_ymax=%s' % (global_ymin, global_ymax))
no_data = numpy.zeros(len(frequency), dtype=numpy.float32)
if self.global_scaling is True:
if map_data[x][y].min() > NoDataThreshold:
median = numpy.median(map_data[x][y])
mad = numpy.median(map_data[x][y] - median)
sigma = map_data[x][y].std()
ymin = median - 2.0 * sigma
ymax = median + 5.0 * sigma
casalog.post('Per panel scaling turned on: ymin=%s, ymax=%s (global ymin=%s, ymax=%s)' % (ymin, ymax, global_ymin, global_ymax))
pl.gcf().sca(self.axes.axes_spmap[y + (self.nh - x - 1) * self.nv])
if map_data[x][y].min() > NoDataThreshold:
pl.plot(frequency, map_data[x][y], color=linecolor, linestyle=linestyle,
if plotmasked == 'empty':
elif plotmasked == 'text':
pl.text((xmin + xmax) / 2.0, (ymin + ymax) / 2.0, 'NO DATA', ha='center', va='center',
size=(self.TickSize + 1))
elif plotmasked == 'zero':
pl.plot(frequency, no_data,
color=maskedcolor, linestyle=linestyle, linewidth=linewidth)
elif plotmasked == 'none':
if self.axes.xlabel_area > 0. and y == 0 and x == 0:
elif plotmasked == 'plot':
m = map_data[x][y] == NoDataThreshold
ma = pl.ma.masked_array(map_data[x][y], m)
color=maskedcolor, linestyle=linestyle, linewidth=linewidth)
pl.axis((xmin, xmax, ymin, ymax))
casalog.post('figfile=\'%s\'' % (figfile), priority='DEBUG')
if figfile is not None and len(figfile) > 0:
casalog.post('Output profile map to %s' % (figfile))
pl.savefig(figfile, format='png', dpi=DPIDetail, transparent=transparent)
class SpectralImage(object):
def __init__(self, imagename):
self.imagename = imagename
self.image_shape = ia.shape()
self.coordsys = ia.coordsys()
coord_types = self.coordsys.axiscoordinatetypes()
self.units = self.coordsys.units()
self.names = self.coordsys.names()
self.direction_reference = self.coordsys.referencecode('dir')[0]
self.spectral_reference = self.coordsys.referencecode('spectral')[0]
self.id_direction = coord_types.index('Direction')
self.id_direction = [self.id_direction, self.id_direction + 1]
self.id_spectral = coord_types.index('Spectral')
self.id_stokes = coord_types.index('Stokes')
stokes_axis_exists = True
# if no Stokes axis exists, set dummy axis at the end
self.id_stokes = len(coord_types)
stokes_axis_exists = False
casalog.post('id_direction=%s' % (self.id_direction), priority='DEBUG')
casalog.post('id_spectral=%s' % (self.id_spectral), priority='DEBUG')
casalog.post('id_stokes=%s (%s stokes axis)' % (self.id_stokes, ('real' if stokes_axis_exists else 'dummy')),
self.data = ia.getchunk()
self.mask = ia.getchunk(getmask=True)
if not stokes_axis_exists:
# put degenerate Stokes axis at the end
self.data = numpy.expand_dims(self.data, axis=-1)
self.mask = numpy.expand_dims(self.mask, axis=-1)
casalog.post('add dummy Stokes axis to data and mask since input image doesn\'t have it.')
casalog.post('data.shape={}'.format(self.data.shape), priority='DEBUG')
casalog.post('mask.shape={}'.format(self.mask.shape), priority='DEBUG')
bottom = ia.toworld(numpy.zeros(len(self.image_shape), dtype=int), 'q')['quantity']
top = ia.toworld(self.image_shape - 1, 'q')['quantity']
key = lambda x: '*%s' % (x + 1)
ra_min = bottom[key(self.id_direction[0])]
ra_max = top[key(self.id_direction[0])]
ra_min, ra_max = ra_max, ra_min
self.dec_min = bottom[key(self.id_direction[1])]
self.dec_max = top[key(self.id_direction[1])]
self._brightnessunit = ia.brightnessunit()
if self.spectral_label == 'Frequency':
refpix, refval, increment = self.spectral_axis(unit='GHz')
self.spectral_data = numpy.array([refval + increment * (i - refpix) for i in range(self.nchan)])
elif self.spectral_label == 'Velocity':
refpix, refval, increment = self.spectral_axis(unit='km/s')
self.spectral_data = numpy.array([refval + increment * (i - refpix) for i in range(self.nchan)])
self.stokes = self.coordsys.stokes()
# if no Stokes axis exists, set ['I']
return self.image_shape[self.id_direction[0]]
return self.image_shape[self.id_direction[1]]
return self.image_shape[self.id_spectral]
return self.image_shape[self.id_stokes]
def brightnessunit(self):
return self._brightnessunit
def direction_label(self):
return [self.names[i] for i in self.id_direction]
def spectral_label(self):
return self.names[self.id_spectral]
return self.units[self.id_spectral]
def to_velocity(self, frequency, freq_unit='GHz', restfreq=None):
rest_frequency = self.coordsys.restfrequency()
# user-defined rest frequency takes priority
vrf = qa.convert(qa.quantity(restfreq), freq_unit)['value']
elif rest_frequency['unit'] != freq_unit:
vrf = qa.convert(rest_frequency, freq_unit)['value']
vrf = rest_frequency['value']
return (1.0 - (frequency / vrf)) * LightSpeed
def to_frequency(self, velocity, freq_unit='km/s'):
rest_frequency = self.coordsys.restfrequency()
if rest_frequency['unit'] != freq_unit:
vrf = qa.convert(rest_frequency, freq_unit)['value']
vrf = rest_frequency['value']
return (1.0 - velocity / LightSpeed) * vrf
def spectral_axis(self, unit='GHz'):
return self.__axis(self.id_spectral, unit=unit)
def direction_axis(self, idx, unit='deg'):
return self.__axis(self.id_direction[idx], unit=unit)
def __axis(self, idx, unit):
refpix = self.coordsys.referencepixel()['numeric'][idx]
refval = self.coordsys.referencevalue()['numeric'][idx]
increment = self.coordsys.increment()['numeric'][idx]
refval = qa.convert(qa.quantity(refval, _unit), unit)['value']
increment = qa.convert(qa.quantity(increment, _unit), unit)['value']
#return numpy.array([refval+increment*(i-refpix) for i in xrange(self.nchan)])
return (refpix, refval, increment)