April 2024 updates (ACES-CMZ#418)

* add a split step for continuum cleaning (even though we've now imaged
95% of the data without taking this step...)

* remove split approach

* verboser text

* allow for target_line.ms to be in filenames but not exist

* pre-split continuum data

* fix up some critical bookkeeping

* missed import, more flexible retrieval script

* reduce requirements for installation so we can use ACES within CASA

* fix typos/errors

* dammit, I deleted a bunch of critical .ms files.

* more safety checks

* bump up cyclefactor for field aa

* pblimit for aa back to 0.1, plus extra logging commands

* add debug verbosity in imstats, which keeps "failing" when it succeeds, and switch pipeline to casa 6.5.4-9

* minor typo

* some mosaic stuff

* more sanity checks in make_mosaics and allow for another exception for field am in cleaning

* remove bad check

* mosaics need a relativethreshold

* fix manual cases

* pbcor things

* giant mosaics of feathered

* try to make mosaics work more...

* add a velocity of peak map.  Try forcing field aj to use original MSes

* whitespace

* revert aj ms selection: just changing to datacolumn=data worked, which doesn't make a whole lot of sense but 🤷

* fix undefined variable error

* add workaround for CI - raise an error at runtime if not in CASA environment

* revise mosaicing file finder

* try to make product-based images work again

* increase verbosity and add manual cubes to max/mom calculations

* p: use datacolumn=data

* ah, v, and l all need to switch to datacolumn='data' now too.

* force datacolumn='data' if that's specified in the tclean pars

* for field i, low-frequency, turn up cyclefactor and limit niter

* same thing for ag

* b was using contsubd data

* whitespace

* mosaicing & giant cube analysis updates

* mosaicing file path fixes, RMS folder choice, and attempt to fix analysis issue

aces/analysis/giantcube_cuts.py

@@ -77,6 +77,16 @@
     makepng(data=mx.value, wcs=mx.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_max.png",
             stretch='asinh', min_percent=0.1, max_percent=99.9)
 
+    print(f"argmax.  dt={time.time() - t0}")
+    argmx = cube.argmax(axis=0, **howargs)
+    vmax = cube.spectral_axis[argmx]
+    hdu = mx.hdu
+    hdu.data = vmax.value
+    hdu.writeto(f"{mompath}/{molname}_CubeMosaic_vpeak.fits", overwrite=True)
+    # use mx.wcs
+    makepng(data=vmax.value, wcs=mx.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_vpeak.png",
+            stretch='asinh', min_percent=0.1, max_percent=99.9)
+
     if dopv:
         print(f"PV max 2.  dt={time.time() - t0}")
         pv_max = cube.max(axis=2, **howargs)
@@ -96,15 +106,21 @@
         make_downsampled_cube(f'{cubepath}/{molname}_CubeMosaic.fits', f'{cubepath}/{molname}_CubeMosaic_downsampled9.fits')
 
     print(f"masked mom0.  dt={time.time() - t0}")
-    std = cube.mad_std()
+    try:
+        std = cube.mad_std()
+    except ValueError:
+        # mad_std requires whole cube in memory; we can't afford that
+        # instead, do a cheap version of sigma clipping
+        std = cube.std()
+        std = cube.with_mask(cube < std * 5).std()
     mom0 = cube.with_mask(cube > std).moment0(axis=0, **howargs)
     mom0.write(f"{mompath}/{molname}_CubeMosaic_masked_mom0.fits", overwrite=True)
     makepng(data=mom0.value, wcs=mom0.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_masked_mom0.png",
             stretch='asinh', min_cut=-0.1, max_percent=99.5)
 
     from dask_image import ndmorph
     signal_mask = cube > std
-    signal_mask = ndmorph.binary_dilation(signal_mask, structure=np.ones([3, 3, 3]), iterations=1)
+    signal_mask = ndmorph.binary_dilation(signal_mask.include(), structure=np.ones([3, 3, 3]), iterations=1)
     mom0 = cube.with_mask(signal_mask).moment0(axis=0, **howargs)
     mom0.write(f"{mompath}/{molname}_CubeMosaic_masked_dilated_mom0.fits", overwrite=True)
     makepng(data=mom0.value, wcs=mom0.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_masked_dilated_mom0.png",

aces/analysis/imstats.py

@@ -497,18 +497,25 @@ def default(self, obj):
 
 
 def savestats(basepath=basepath,
-              suffix='image.tt0*', filetype=".fits"):
+              suffix='image.tt0*', filetype=".fits",
+              verbose=True
+             ):
     """
     filtype : str
         ".fits" or "" for CASA image
     """
 
+    if verbose:
+        print("Assembling statistics")
     stats = assemble_stats(
         f"{basepath}/data/2021.1.00172.L/science_goal.uid___A001_X1590_X30a8/group.uid___A001_X1590_X30a9/*/calibrated/working/*.cont.I.iter1.{suffix}{filetype}",
         ditch_suffix=f".{suffix[:-1]}")
     with open(f'{basepath}/tables/metadata_{suffix}.json', 'w') as fh:
         json.dump(stats, fh, cls=MyEncoder)
 
+    if verbose:
+        print("Done assembling statistics & dumping them to file")
+
     requested = get_requested_sens()
 
     meta_keys = ['region', 'band', 'array', 'robust', 'suffix',
@@ -521,6 +528,9 @@ def savestats(basepath=basepath,
     req_keys = ['B3_res', 'B3_sens', ]
     req_keys_head = ['Req_Res', 'Req_Sens']
 
+    if verbose:
+        print("Assembling table rows")
+
     rows = []
     for entry in stats:
         band = entry['meta']['band']
@@ -540,15 +550,20 @@ def savestats(basepath=basepath,
     tbl.add_column(Column(name='BeamVsReq', data=(tbl['bmaj'] * tbl['bmin'])**0.5 / tbl['Req_Res']))
     tbl.add_column(Column(name='BmajVsReq', data=tbl['bmaj'] / tbl['Req_Res']))
 
+    if verbose:
+        print("Writing tables")
     tbl.write(f'{basepath}/tables/metadata_{suffix.strip("*")}.ecsv', overwrite=True)
     tbl.write(f'{basepath}/tables/metadata_{suffix.strip("*")}.html',
               format='ascii.html', overwrite=True)
     tbl.write(f'{basepath}/tables/metadata_{suffix.strip("*")}.tex', overwrite=True)
     tbl.write(f'{basepath}/tables/metadata_{suffix.strip("*")}.js.html',
               format='jsviewer', overwrite=True)
+    if verbose:
+        print("Done writing tables")
 
     return tbl
 
 
 def main():
     return savestats()
+    print("Finished with imstats.main()")

aces/analysis/parse_contdotdat.py

@@ -191,7 +191,8 @@ def cont_channel_selection_to_contdotdat(cont_channel_selection, msname,
 
 except ModuleNotFoundError:
     # if using this not in a casa6 environment
-    pass
+    def contchannels_to_linechannels(*args, **kwargs):
+        raise NotImplementedError
 
 # flagchannels='0:0~60;180~300;2790~2880;3280~3360;3460~3490;3830~3839,1:60~130;200~250;320~420;580~650;1000~1040;1200~1360;1420~1460;1720~1790;1860~1919,2:40~300;630~700;800~1000;1440~1640;1780~1919,3:100~150;470~540;640~820;920~980;1220~1260;1370~1420;1710~1780,4:0~60;180~300;2790~2880;3280~3360;3460~3490;3830~3839,5:60~130;200~250;320~420;580~650;1000~1040;1200~1360;1420~1460;1720~1790;1860~1919,6:40~300;630~700;800~1000;1440~1640;1780~1919,7:100~150;470~540;640~820;920~980;1220~1260;1370~1420;1710~1780,8:0~60;180~300;2790~2880;3280~3360;3460~3490;3830~3839,9:60~130;200~250;320~420;580~650;1000~1040;1200~1360;1420~1460;1720~1790;1860~1919,10:40~300;630~700;800~1000;1440~1640;1780~1919,11:100~150;470~540;640~820;920~980;1220~1260;1370~1420;1710~1780,12:0~60;180~300;900~1050;1860~1950;2100~2140;2230~2280;2790~2880;3050~3100;3280~3360;3460~3490;3590~3650;3830~3839,13:60~130;200~250;265~285;320~420;435~460;580~650;670~700;760~810;1000~1040;1200~1360;1420~1460;1720~1790;1800~1840;1860~1919,14:40~300;630~700;800~1000;1440~1640;1780~1919,15:100~150;470~540;640~820;920~980;1170~1190;1220~1260;1370~1420;1710~1780'
 # freqsel = cont_channel_selection_to_contdotdat(flagchannels, ('./science_goal.uid___A001_X1290_X44/group.uid___A001_X1290_X45/member.uid___A001_X1290_X46/calibrated/calibrated_final.ms/'), spw_mapping={0:25,1:27,2:29,3:31})

aces/analysis/statcont_cubes.py

@@ -261,9 +261,11 @@ def main():
                 try:
                     SpectralCube.read(fn).beam
                 except NoBeamError:
+                    print(f"Neither {fn} nor {contsubfn} have a beam")
                     raise ValueError(f"Neither {fn} nor {contsubfn} have a beam")
                 except AttributeError:
                     SpectralCube.read(fn).beams
+                    print(f"{fn} is a multi-beam cube")
                     raise AttributeError(f"{fn} is a multi-beam cube")
                 redo = True
             except AttributeError:

aces/hipergator_scripts/job_runner.py

@@ -299,6 +299,7 @@ def main():
                     if continue_started_only:
                         msfiles = glob.glob(f'{workdir}/{tempdir_name}/*.ms')
                         if len(msfiles) == 0:
+                            print("No MS files found; skipping")
                             continue
                         else:
                             msfstr = '\n'.join(msfiles)
@@ -333,11 +334,16 @@ def main():
                                          else os.path.join(datapath, os.path.basename(vis)).replace("targets", "target")
                                          for vis in tcpars['vis']]
                         for ii, vis in enumerate(tcpars['vis']):
-                            try:
-                                assert os.path.exists(vis), vis
-                            except AssertionError:
-                                assert os.path.exists(vis.replace("_line", ""))
-                                tcpars['vis'][ii] = vis.replace("_line", "")
+                            if not os.path.exists(vis):
+                                if os.path.exists(vis.replace("_line", "")):
+                                    tcpars['vis'][ii] = vis.replace("_line", "")
+                                elif os.path.exists(vis.replace("_target_line", "")):
+                                    tcpars['vis'][ii] = vis.replace("_target_line", "")
+                                elif os.path.exists(vis.replace("_target", "")):
+                                    tcpars['vis'][ii] = vis.replace("_target", "")
+                                else:
+                                    matches = glob.glob(vis.replace(".ms", "*"))
+                                    raise IOError(f"MS {vis} does not exist, nor do any of its variants.  glob={matches}")
                         if 'nchan' not in tcpars:
                             print(tcpars)
                             print(sous, gous, mous)

aces/hipergator_scripts/run_mosaic_12m.sh

@@ -8,8 +8,8 @@
 #SBATCH --time=96:00:00               # Time limit hrs:min:sec
 #SBATCH --qos=astronomy-dept-b
 #SBATCH --account=astronomy-dept
-#SBATCH --output=/blue/adamginsburg/adamginsburg/ACES/logs/ACES_mosaic12m_cont_%j.log
-#SBATCH --job-name=ACES_mosaic12m_cont
+#SBATCH --output=/blue/adamginsburg/adamginsburg/ACES/logs/ACES_mosaic12m_%j.log
+#SBATCH --job-name=ACES_mosaic12m
 #SBATCH --export=ALL
 
 
@@ -42,7 +42,7 @@ echo "test import"
 #echo "Make TP mosaic"
 #/orange/adamginsburg/miniconda3/envs/python39/bin/aces_mosaic_TP || exit 1
 echo "Make 12m mosaic"
-/orange/adamginsburg/miniconda3/envs/python39/bin/aces_mosaic_12m --contonly || exit 1
+/orange/adamginsburg/miniconda3/envs/python39/bin/aces_mosaic_12m || exit 1
 
 # technically shouldn't need to be re-run, but as I add new mosaics, it will
 # (but this causes a 'failed' error)

aces/hipergator_scripts/run_pipeline_slurm_burst_mpi.sh

@@ -29,6 +29,8 @@ export ACES_ROOTDIR="/orange/adamginsburg/ACES/"
 CASAVERSION=casa-6.2.1-7-pipeline-2021.2.0.128
 CASAVERSION=casa-6.4.3-2-pipeline-2021.3.0.17
 CASAVERSION=casa-6.4.1-12-pipeline-2022.2.0.68
+# required for QA3 version of field af
+CASAVERSION=casa-6.5.4-9-pipeline-2023.1.0.125
 
 export CASAPATH=/orange/adamginsburg/casa/${CASAVERSION}
 export MPICASA=${CASAPATH}/bin/mpicasa

aces/imaging/make_mosaic.py

@@ -79,16 +79,21 @@ def read_as_2d(fn, minval=None):
 
 
 def get_peak(fn, slab_kwargs=None, rest_value=None, suffix="", save_file=True,
-             folder=None, threshold=None):
+             folder=None, threshold=None, rel_threshold=None, fail_on_zeros=True):
     print(".", end='', flush=True)
     outfn = fn.replace(".fits", "") + f"{suffix}_max.fits"
     if folder is not None:
         outfn = os.path.join(folder, os.path.basename(outfn))
     if os.path.exists(outfn):
         hdu = fits.open(outfn)
         proj = Projection.from_hdu(hdu)
+        if rel_threshold is not None:
+            threshold = rel_threshold * proj.max()
+            print(f"Set threshold to {threshold} based on rel_threshold={rel_threshold}")
         if threshold is not None:
             proj[proj < threshold] = 0
+        if fail_on_zeros and np.nansum(proj) == 0:
+            raise ValueError(f"File {fn} reduced to all zeros")
         return proj
     else:
         ft = 'fits' if fn.endswith(".fits") else "casa_image"
@@ -116,8 +121,14 @@ def get_peak(fn, slab_kwargs=None, rest_value=None, suffix="", save_file=True,
                         mxjy = mxjy.with_beam(beam, raise_error_jybm=False)
 
                     mx = mxjy.to(u.K, equivalencies=equiv)
+        if fail_on_zeros and np.nansum(mx.value) == 0:
+            raise ValueError(f"File {fn} reduced to all zeros")
         if save_file:
             mx.hdu.writeto(outfn)
+
+        if rel_threshold is not None:
+            threshold = rel_threshold * mx.max()
+            print(f"Set threshold to {threshold} based on rel_threshold={rel_threshold}")
         if threshold is not None:
             mx[mx.value < threshold] = 0
         return mx
@@ -153,6 +164,20 @@ def get_m0(fn, slab_kwargs=None, rest_value=None, suffix="", save_file=True, fol
         return moment0
 
 
+def check_hdus(hdus):
+    bad = 0
+    for hdu in hdus:
+        if isinstance(hdu, fits.PrimaryHDU):
+            ttl = np.nansum(hdu.data)
+        else:
+            ttl = np.nansum(hdu)
+        if ttl == 0:
+            bad = bad + 1
+
+    if bad > 0:
+        raise ValueError(f"Found {bad} bad HDUs (all zero)")
+
+
 def makepng(data, wcs, imfn, footprint=None, **norm_kwargs):
     import pylab as pl
     from astropy import visualization
@@ -192,6 +217,9 @@ def make_mosaic(twod_hdus, name, norm_kwargs={}, slab_kwargs=None,
                 array='7m',
                 folder=None,  # must be specified though
                 basepath='./'):
+    """
+    Given a long list of 2D HDUs and an output name, make a giant mosaic.
+    """
 
     if target_header is None:
         log.info(f"Finding WCS for {len(twod_hdus)} HDUs")
@@ -257,6 +285,10 @@ def repr_function(*args, **kwargs):
     if commonbeam is not None:
         header.update(cb.to_header_keywords())
 
+    assert not np.all(np.isnan(prjarr))
+
+    log.info(f"DEBUG: footprint[prjarr==0] = {footprint[prjarr == 0]}")
+
     outfile = f'{basepath}/mosaics/{folder}/{array}_{name}_mosaic.fits'
     log.info(f"Writing reprojected data to {outfile}")
     fits.PrimaryHDU(data=prjarr, header=header).writeto(outfile, overwrite=True)
@@ -367,7 +399,7 @@ def repr_function(*args, **kwargs):
         return cb
 
 
-def all_lines(header, parallel=False, array='12m', glob_suffix='cube.I.iter1.image.pbcor',
+def all_lines(header, parallel=False, array='12m', glob_suffixes=('cube.I.iter1.image.pbcor.statcont.contsub.fits', 'cube.I.manual.image.pbcor.statcont.contsub.fits'),
               lines='all', folder='',
               globdir='calibrated/working', use_weights=True):
 
@@ -388,10 +420,14 @@ def all_lines(header, parallel=False, array='12m', glob_suffix='cube.I.iter1.ima
 
         log.info(f"{array} {line} {restf}")
 
-        filelist = glob.glob(f'{basepath}/rawdata/2021.1.00172.L/s*/g*/m*/{globdir}/*spw{spwn}.{glob_suffix}')
+        filelist = []
+        for glob_suffix in glob_suffixes:
+            for spwpre in ('spw', 'sci'):
+                filelist += glob.glob(f'{basepath}/rawdata/2021.1.00172.L/s*/g*/m*/{globdir}/*{spwpre}{spwn}.{glob_suffix}')
+        log.debug(f"Filelist={filelist}")
         assert len(filelist) > 0
         if use_weights:
-            weightfiles = [fn.replace("image.pbcor", "weight") for fn in filelist]
+            weightfiles = [fn.replace("image.pbcor.statcont.contsub.fits", "weight") for fn in filelist]
             for ii, (ifn, wfn) in enumerate(zip(list(filelist), list(weightfiles))):
                 if not os.path.exists(wfn):
                     log.error(f"Missing file {wfn}")
@@ -408,36 +444,42 @@ def all_lines(header, parallel=False, array='12m', glob_suffix='cube.I.iter1.ima
                             filelist,
                             )
             hdus = [x.hdu for x in hdus]
+            check_hdus(hdus)
             if use_weights:
                 wthdus = pool.map(partial(get_peak,
                                           **{'slab_kwargs': {'lo': -2 * u.km / u.s, 'hi': 2 * u.km / u.s},
                                              'rest_value': restf},
                                           suffix=f'_{line}',
-                                          threshold=0.5,  # pb limit
+                                          rel_threshold=0.25,  # pb limit
                                           ),
                                   weightfiles)
                 wthdus = [x.hdu for x in wthdus]
+                check_hdus(wthdus)
         else:
             hdus = [get_peak(fn, slab_kwargs={'lo': -200 * u.km / u.s, 'hi': 200 * u.km / u.s},
                              rest_value=restf, suffix=f'_{line}').hdu for fn in filelist]
+            check_hdus(hdus)
             print(flush=True)
             if use_weights:
                 wthdus = [get_peak(fn, slab_kwargs={'lo': -2 * u.km / u.s, 'hi': 2 * u.km / u.s},
                                    rest_value=restf, suffix=f'_{line}',
-                                   threshold=0.5,  # pb limit
+                                   rel_threshold=0.25,  # pb limit
                                    ).hdu for fn in weightfiles]
+                check_hdus(wthdus)
                 print(flush=True)
 
         if parallel:
-            print(f"Starting function make_mosaic for {line} peak intensity")
+            print(f"Starting function make_mosaic for {line} peak intensity (parallel mode)")
             proc = Process(target=make_mosaic, args=(hdus, ),
                            kwargs=dict(name=f'{line}_max', cbar_unit='K',
                            norm_kwargs=dict(max_cut=5, min_cut=-0.01, stretch='asinh'),
                            array=array, basepath=basepath, weights=wthdus if use_weights else None,
+                           folder=folder,
                            target_header=header))
             proc.start()
             processes.append(proc)
         else:
+            print(f"Starting function make_mosaic for {line} peak intensity")
             make_mosaic(hdus, name=f'{line}_max', cbar_unit='K',
                         norm_kwargs=dict(max_cut=5, min_cut=-0.01, stretch='asinh'),
                         array=array, basepath=basepath, weights=wthdus if use_weights else None,
@@ -458,13 +500,16 @@ def all_lines(header, parallel=False, array='12m', glob_suffix='cube.I.iter1.ima
             proc = Process(target=make_mosaic, args=(m0hdus, ),
                            kwargs=dict(name=f'{line}_m0', cbar_unit='K km/s',
                            norm_kwargs=dict(max_cut=20, min_cut=-1, stretch='asinh'),
+                           folder=folder,
                            array=array, basepath=basepath, weights=wthdus if use_weights else None, target_header=header))
             proc.start()
             processes.append(proc)
         else:
             make_mosaic(m0hdus, name=f'{line}_m0', cbar_unit='K km/s',
                         norm_kwargs=dict(max_cut=20, min_cut=-1, stretch='asinh'),
-                        array=array, basepath=basepath, weights=wthdus if use_weights else None, target_header=header)
+                        array=array, basepath=basepath,
+                        folder=folder,
+                        weights=wthdus if use_weights else None, target_header=header)
             del m0hdus
             if use_weights:
                 del wthdus