From d9b04b2e246b173ee4333292dc63d6863921abc4 Mon Sep 17 00:00:00 2001 From: Takeshi Nakazato Date: Mon, 4 Aug 2025 07:44:03 +0000 Subject: [PATCH 1/5] Fix wrong frequency label of output image - fixed wrong reference channel index - convert TOPO frequency into LSRK --- python/priism/alma/imager.py | 56 +++++++++++++++++++++++-------- python/priism/alma/imagewriter.py | 29 +++++++++++++--- 2 files changed, 66 insertions(+), 19 deletions(-) diff --git a/python/priism/alma/imager.py b/python/priism/alma/imager.py index 0c0e012..8a74b3c 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 @@ -214,7 +216,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,23 +232,46 @@ 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, + + if (isinstance(self.imparam.start, str) and self.imparam.start.isdigit()) \ + or isinstance(self.imparam.start, int): + start_index = self.imparam.start + assert self.imparam.nchan == 1 + chan_width = self.imparam.width + spw = int(mssel_index['spw'][0]) + spw_chan = mssel_index['channel'][0] + chan_list = np.arange(spw_chan[1], spw_chan[2] + 1, spw_chan[3]) + start = chan_list[start_index] + nchan = chan_width + end = chan_list[start_index + nchan - 1] + channel = (start + end) / 2.0 + print(f'Use Freuquency for channel {channel} spw {spw}') + cf_vis, cw = imagewriter.ImageWriter.frequency_setup_for_spw( + vis=vis, + spw_id=spw, + channel=channel + ) + reference_time = imagemeta.observing_date + phase_direction = imparam_for_writer.phasecenter + observatory_position = imagemeta.telescope_position + cf_lsrk = imagewriter.ImageWriter.to_lsrk( + cf_vis, + reference_time, + phase_direction, + observatory_position + ) + cf_new = cf_lsrk['m0'] + imparam_for_writer.start = f'{cf_new["value"]:16.12f}{cf_new["unit"]}' + width = nchan * cw + imparam_for_writer.width = f'{width:16.12f}Hz' + + writer = imagewriter.ImageWriter(imparam_for_writer, self.imagearray.data, imagemeta) writer.write(overwrite=overwrite) diff --git a/python/priism/alma/imagewriter.py b/python/priism/alma/imagewriter.py index 85c185b..5b9fe0c 100644 --- a/python/priism/alma/imagewriter.py +++ b/python/priism/alma/imagewriter.py @@ -34,11 +34,30 @@ def phase_direction_for_field(vis, 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 frequency_setup_for_spw(vis: str, spw_id: int, channel: float) -> tuple[dict, float]: + me = casa.CreateCasaMeasure() + qa = casa.CreateCasaQuantity() + with casa.OpenMSMetaData(vis) as msmd: + chan_freq = msmd.chanfreqs(spw_id) + chan_width = msmd.chanwidths(spw_id) + ref_freq = msmd.reffreq(spw_id) + freq_frame = me.getref(ref_freq) + xp = np.arange(len(chan_freq)) + cf = np.interp(channel, xp, chan_freq) + cw = np.interp(channel, xp, chan_width) + freq_quantity = qa.quantity(cf, 'Hz') + cf_vis = me.frequency(rf=freq_frame, v0=freq_quantity) + return cf_vis, cw + + @staticmethod + def to_lsrk(cf: dict, reference_time: dict, phase_direction: dict, observatory_position: dict) -> dict: + me = casa.CreateCasaMeasure() + me.done() + me.doframe(reference_time) + me.doframe(phase_direction) + me.doframe(observatory_position) + cf_lsrk = me.measure(cf, 'LSRK') + return cf_lsrk def __init__(self, imageparam, imagearray, imagemeta=None): self.imageparam = imageparam From 1fa1ee7b83429a3c3f115f56d965a0c7d78d5ca0 Mon Sep 17 00:00:00 2001 From: Takeshi Nakazato Date: Fri, 29 Aug 2025 05:34:11 +0000 Subject: [PATCH 2/5] Fail if nchan != 1 --- python/priism/alma/imager.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/python/priism/alma/imager.py b/python/priism/alma/imager.py index 8a74b3c..5ad49d7 100644 --- a/python/priism/alma/imager.py +++ b/python/priism/alma/imager.py @@ -118,7 +118,7 @@ def defineimage(self, imsize=100, cell='1arcsec', phasecenter='', projection='SI | |-------|-------|-------| nchan=3 |<----->| - width= + width= Parameters: @@ -135,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): From 1f1185210ddaf09f7ce5de14a73be95aea98b2fb Mon Sep 17 00:00:00 2001 From: Takeshi Nakazato Date: Fri, 29 Aug 2025 07:35:29 +0000 Subject: [PATCH 3/5] support multiple MS inputs and multiple data selections to single MS --- python/priism/alma/imager.py | 48 +++++++++++----------------- python/priism/alma/imagewriter.py | 26 --------------- python/priism/alma/paramcontainer.py | 41 ++++++++++++++++++++++++ 3 files changed, 60 insertions(+), 55 deletions(-) diff --git a/python/priism/alma/imager.py b/python/priism/alma/imager.py index 5ad49d7..2520ee3 100644 --- a/python/priism/alma/imager.py +++ b/python/priism/alma/imager.py @@ -247,38 +247,28 @@ def exportimage(self, imagename, overwrite=False): if (isinstance(self.imparam.start, str) and self.imparam.start.isdigit()) \ or isinstance(self.imparam.start, int): - start_index = self.imparam.start assert self.imparam.nchan == 1 - chan_width = self.imparam.width - spw = int(mssel_index['spw'][0]) - spw_chan = mssel_index['channel'][0] - chan_list = np.arange(spw_chan[1], spw_chan[2] + 1, spw_chan[3]) - start = chan_list[start_index] - nchan = chan_width - end = chan_list[start_index + nchan - 1] - channel = (start + end) / 2.0 - print(f'Use Freuquency for channel {channel} spw {spw}') - cf_vis, cw = imagewriter.ImageWriter.frequency_setup_for_spw( - vis=vis, - spw_id=spw, - channel=channel - ) - reference_time = imagemeta.observing_date - phase_direction = imparam_for_writer.phasecenter - observatory_position = imagemeta.telescope_position - cf_lsrk = imagewriter.ImageWriter.to_lsrk( - cf_vis, - reference_time, - phase_direction, - observatory_position - ) - cf_new = cf_lsrk['m0'] - imparam_for_writer.start = f'{cf_new["value"]:16.12f}{cf_new["unit"]}' - width = nchan * cw + 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) + imparam_for_writer.start = f'{lsrk_freq_min:16.12f}Hz' + width = lsrk_freq_max - lsrk_freq_min imparam_for_writer.width = f'{width:16.12f}Hz' - writer = imagewriter.ImageWriter(imparam_for_writer, self.imagearray.data, - imagemeta) + 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 5b9fe0c..9ec7650 100644 --- a/python/priism/alma/imagewriter.py +++ b/python/priism/alma/imagewriter.py @@ -33,32 +33,6 @@ def phase_direction_for_field(vis, field_id): pdir = msmd.phasecenter(field_id) return pdir - @staticmethod - def frequency_setup_for_spw(vis: str, spw_id: int, channel: float) -> tuple[dict, float]: - me = casa.CreateCasaMeasure() - qa = casa.CreateCasaQuantity() - with casa.OpenMSMetaData(vis) as msmd: - chan_freq = msmd.chanfreqs(spw_id) - chan_width = msmd.chanwidths(spw_id) - ref_freq = msmd.reffreq(spw_id) - freq_frame = me.getref(ref_freq) - xp = np.arange(len(chan_freq)) - cf = np.interp(channel, xp, chan_freq) - cw = np.interp(channel, xp, chan_width) - freq_quantity = qa.quantity(cf, 'Hz') - cf_vis = me.frequency(rf=freq_frame, v0=freq_quantity) - return cf_vis, cw - - @staticmethod - def to_lsrk(cf: dict, reference_time: dict, phase_direction: dict, observatory_position: dict) -> dict: - me = casa.CreateCasaMeasure() - me.done() - me.doframe(reference_time) - me.doframe(phase_direction) - me.doframe(observatory_position) - cf_lsrk = me.measure(cf, 'LSRK') - return cf_lsrk - def __init__(self, imageparam, imagearray, imagemeta=None): self.imageparam = imageparam self.imagearray = imagearray diff --git a/python/priism/alma/paramcontainer.py b/python/priism/alma/paramcontainer.py index b950dfc..68c3ded 100644 --- a/python/priism/alma/paramcontainer.py +++ b/python/priism/alma/paramcontainer.py @@ -84,6 +84,47 @@ 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) + position = msmd.observatoryposition() + freqs = (chan_freq_for_spw[i] for i in (chan_start, chan_end)) + 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): """ From 58528b0044daaa0ef0b7e57b5b36e67e6ed8b90d Mon Sep 17 00:00:00 2001 From: Takeshi Nakazato Date: Thu, 6 Nov 2025 00:59:02 +0000 Subject: [PATCH 4/5] Additional fix to frequency shift issue --- python/priism/alma/imager.py | 2 +- python/priism/alma/imagewriter.py | 24 +++++++++++------------- python/priism/alma/paramcontainer.py | 6 +++++- 3 files changed, 17 insertions(+), 15 deletions(-) diff --git a/python/priism/alma/imager.py b/python/priism/alma/imager.py index 2520ee3..7b74adf 100644 --- a/python/priism/alma/imager.py +++ b/python/priism/alma/imager.py @@ -260,8 +260,8 @@ def exportimage(self, imagename, overwrite=False): ) lsrk_freq_min = min(lsrk_freq_min, _freq_min) lsrk_freq_max = max(lsrk_freq_max, _freq_max) - imparam_for_writer.start = f'{lsrk_freq_min:16.12f}Hz' 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( diff --git a/python/priism/alma/imagewriter.py b/python/priism/alma/imagewriter.py index 9ec7650..d3fe80e 100644 --- a/python/priism/alma/imagewriter.py +++ b/python/priism/alma/imagewriter.py @@ -98,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 68c3ded..8892075 100644 --- a/python/priism/alma/paramcontainer.py +++ b/python/priism/alma/paramcontainer.py @@ -108,8 +108,12 @@ def to_lsrk_range(self, spw_chan: list[int], phasecenter: dict) -> tuple[float, 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[i] for i in (chan_start, chan_end)) + freqs = ( + chan_freq_for_spw[chan_start] - chan_width_for_spw[spw_id] / 2, + chan_freq_for_spw[chan_end] + chan_width_for_spw[spw_id] / 2 + ) me.done() me.doframe(position) me.doframe(start_time) From 9e3fa1407f19c306f465bd268e397c182031b9c1 Mon Sep 17 00:00:00 2001 From: Takeshi Nakazato Date: Thu, 6 Nov 2025 02:50:42 +0000 Subject: [PATCH 5/5] fix indexing bug --- python/priism/alma/paramcontainer.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/priism/alma/paramcontainer.py b/python/priism/alma/paramcontainer.py index 8892075..fc02a9f 100644 --- a/python/priism/alma/paramcontainer.py +++ b/python/priism/alma/paramcontainer.py @@ -111,8 +111,8 @@ def to_lsrk_range(self, spw_chan: list[int], phasecenter: dict) -> tuple[float, chan_width_for_spw = msmd.chanwidths(spw_id) position = msmd.observatoryposition() freqs = ( - chan_freq_for_spw[chan_start] - chan_width_for_spw[spw_id] / 2, - chan_freq_for_spw[chan_end] + chan_width_for_spw[spw_id] / 2 + 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)