Fix wrong frequency label of output image

python/priism/alma/imager.py

@@ -16,6 +16,8 @@
 # along with PRIISM.  If not, see <https://www.gnu.org/licenses/>.
 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=<constant channel width of image channel>
+          width=<number of visibility channels to be accumulated into one image channel>
 
 
         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):

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']

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):
     """