diff --git a/python/priism/alma/imager.py b/python/priism/alma/imager.py
index 0c0e012..7b74adf 100644
--- a/python/priism/alma/imager.py
+++ b/python/priism/alma/imager.py
@@ -16,6 +16,8 @@
# along with PRIISM. If not, see .
from __future__ import absolute_import
+import copy
+
import numpy as np
from . import paramcontainer
@@ -116,7 +118,7 @@ def defineimage(self, imsize=100, cell='1arcsec', phasecenter='', projection='SI
|
|-------|-------|-------| nchan=3
|<----->|
- width=
+ width=
Parameters:
@@ -133,6 +135,12 @@ def defineimage(self, imsize=100, cell='1arcsec', phasecenter='', projection='SI
stokes stokes parameter (fixed to 'I')
"""
self.imparam = paramcontainer.ImageParamContainer.CreateContainer(**locals())
+ if self.imparam.nchan != 1:
+ raise RuntimeError(
+ 'Only nchan=1 is supported in this version. '
+ 'If you mean to accumulate multiple visibility channels into '
+ 'one image channel, please use width instead of nchan.'
+ )
def configuregrid(self, convsupport, convsampling, gridfunction):
if isinstance(gridfunction, str):
@@ -214,7 +222,10 @@ def exportimage(self, imagename, overwrite=False):
"""
if self.imparam is None:
raise RuntimeError('You have to define image configuration before export!')
- self.imparam.imagename = imagename
+
+ imparam_for_writer = copy.deepcopy(self.imparam)
+
+ imparam_for_writer.imagename = imagename
if self.imagearray is None:
raise RuntimeError('You don\'t have an image array!')
@@ -227,24 +238,37 @@ def exportimage(self, imagename, overwrite=False):
field_id = int(self.imparam.phasecenter)
phase_direction = imagewriter.ImageWriter.phase_direction_for_field(vis=vis,
field_id=field_id)
- self.imparam.phasecenter = phase_direction
- if (isinstance(self.imparam.start, str) and self.imparam.start.isdigit()) \
- or isinstance(self.imparam.start, int):
- # TODO: we need LSRK frequency
- start = self.imparam.start
- spw = int(self.visparams[0].as_msindex()['spw'][0])
- print('Use Freuquency for channel {0} spw {1}'.format(start, spw))
- cf, cw = imagewriter.ImageWriter.frequency_setup_for_spw(vis=vis,
- spw_id=spw,
- chan=start)
- self.imparam.start = cf
- self.imparam.width = cw
+ imparam_for_writer.phasecenter = phase_direction
+
mssel_index = self.visparams[0].as_msindex()
field = None if len(mssel_index['field']) == 0 else mssel_index['field'][0]
spw_selected = None if len(mssel_index['spw']) == 0 else mssel_index['spw']
imagemeta = paramcontainer.ImageMetaInfoContainer.fromvis(vis, field, spw_selected)
- writer = imagewriter.ImageWriter(self.imparam, self.imagearray.data,
- imagemeta)
+
+ if (isinstance(self.imparam.start, str) and self.imparam.start.isdigit()) \
+ or isinstance(self.imparam.start, int):
+ assert self.imparam.nchan == 1
+ lsrk_freq_min, lsrk_freq_max = 1e100, 0
+ for visparam in self.visparams:
+ _mssel_index = visparam.as_msindex()
+ spw_chan_list = _mssel_index['channel']
+ print(f'spw_chan_list: {spw_chan_list}')
+ for spw_chan in spw_chan_list:
+ _freq_min, _freq_max = visparam.to_lsrk_range(
+ spw_chan,
+ imparam_for_writer.phasecenter
+ )
+ lsrk_freq_min = min(lsrk_freq_min, _freq_min)
+ lsrk_freq_max = max(lsrk_freq_max, _freq_max)
+ width = lsrk_freq_max - lsrk_freq_min
+ imparam_for_writer.start = f'{lsrk_freq_min + width / 2:16.12f}Hz'
+ imparam_for_writer.width = f'{width:16.12f}Hz'
+
+ writer = imagewriter.ImageWriter(
+ imparam_for_writer,
+ self.imagearray.data,
+ imagemeta
+ )
writer.write(overwrite=overwrite)
def getimage(self, imagename):
diff --git a/python/priism/alma/imagewriter.py b/python/priism/alma/imagewriter.py
index 85c185b..d3fe80e 100644
--- a/python/priism/alma/imagewriter.py
+++ b/python/priism/alma/imagewriter.py
@@ -33,13 +33,6 @@ def phase_direction_for_field(vis, field_id):
pdir = msmd.phasecenter(field_id)
return pdir
- @staticmethod
- def frequency_setup_for_spw(vis, spw_id, chan):
- with casa.OpenTableForRead(os.path.join(vis, 'SPECTRAL_WINDOW')) as tb:
- chan_freq = tb.getcell('CHAN_FREQ', spw_id)
- chan_width = tb.getcell('CHAN_WIDTH', spw_id)
- return '{0:16.12f}Hz'.format(chan_freq[chan]), '{0:16.12f}Hz'.format(chan_width[chan])
-
def __init__(self, imageparam, imagearray, imagemeta=None):
self.imageparam = imageparam
self.imagearray = imagearray
@@ -105,29 +98,27 @@ def _setup_coordsys(self):
# spectral coordinate
refframe = 'LSRK'
- print('start {0} width {1}'.format(self.imageparam.start,
- self.imageparam.width))
+ print(f'start {self.imageparam.start} width {self.imageparam.width}')
start = qa.convert(self.imageparam.start, 'Hz')
width = qa.convert(self.imageparam.width, 'Hz')
nchan = self.imageparam.nchan
- f = np.fromiter((start['value'] + i * width['value'] for i in range(nchan)), dtype=np.float64)
- print('f = {0}'.format(f))
- frequencies = qa.quantity(f, 'Hz')
veldef = 'radio'
- if len(f) > 1:
+ if nchan > 1:
+ f = np.fromiter((start['value'] + i * width['value'] for i in range(nchan)), dtype=np.float64)
+ print('f = {0}'.format(f))
+ frequencies = qa.quantity(f, 'Hz')
c.setspectral(refcode=refframe,
frequencies=frequencies,
doppler=veldef)
else:
print('set increment for spectral axis: {0}'.format(width))
- #c.setreferencepixel(value=0, type='spectral')
- #c.setreferencevalue(value=start, type='spectral')
- #c.setincrement(value=width, type='spectral')
r = c.torecord()
- if 'spectral2' in r:
- key = 'spectral2'
- elif 'spectral1' in r:
- key = 'spectral1'
+ for k in r:
+ if k.startswith('spectral'):
+ key = k
+ break
+ else:
+ raise RuntimeError('spectral axis not found in coordinate system')
r[key]['wcs']['crpix'] = 0.0
r[key]['wcs']['crval'] = start['value']
r[key]['wcs']['cdelt'] = width['value']
diff --git a/python/priism/alma/paramcontainer.py b/python/priism/alma/paramcontainer.py
index b950dfc..fc02a9f 100644
--- a/python/priism/alma/paramcontainer.py
+++ b/python/priism/alma/paramcontainer.py
@@ -84,6 +84,51 @@ def as_msindex(self):
idx = ms.msselectedindices()
return idx
+ def to_lsrk_range(self, spw_chan: list[int], phasecenter: dict) -> tuple[float, float]:
+ """
+ Convert channel selection specified by spw_chan
+ to LSRK frequency range.
+
+ Parameters:
+ spw_chan: spw channel selection in the form of
+ [spw_id, start_chan, end_chan, chan_increment].
+ end_chan should be inclusive.
+ phasecenter: phase center direction measure.
+
+ Returns:
+ Minimum and maximum LSRK frequencies in Hz.
+ """
+ me = casa.CreateCasaMeasure()
+ qa = casa.CreateCasaQuantity()
+ spw_id, chan_start, chan_end, _ = spw_chan
+ with casa.OpenMSMetaData(self.vis) as msmd:
+ time_list_for_spw = msmd.timesforspws(spw_id)
+ start_time = me.epoch(
+ 'UTC',
+ qa.quantity(time_list_for_spw.min(), 's')
+ )
+ chan_freq_for_spw = msmd.chanfreqs(spw_id)
+ chan_width_for_spw = msmd.chanwidths(spw_id)
+ position = msmd.observatoryposition()
+ freqs = (
+ chan_freq_for_spw[chan_start] - chan_width_for_spw[chan_start] / 2,
+ chan_freq_for_spw[chan_end] + chan_width_for_spw[chan_end] / 2
+ )
+ me.done()
+ me.doframe(position)
+ me.doframe(start_time)
+ me.doframe(phasecenter)
+ topo_freqs = map(lambda f: qa.quantity(f, 'Hz'), freqs)
+ converted = map(
+ lambda q: me.measure(me.frequency('TOPO', q), 'LSRK')['m0'],
+ topo_freqs
+ )
+ lsrk_freqs = map(
+ lambda q: qa.convert(q, 'Hz')['value'],
+ converted
+ )
+ return tuple(sorted(lsrk_freqs))
+
class ImageParamContainer(base_container.ParamContainer):
"""