Feb 2025 edits from me

aces/analysis/downsample_all_spectrally.py

@@ -0,0 +1,69 @@
+import os
+import shutil
+from aces.imaging.make_mosaic import downsample_spectrally
+from aces import conf
+
+basepath = conf.basepath
+
+dsfactor_dict = {
+    "HNCO_7m12mTP": 7,
+    "HCOP_noTP": 7,
+    #"CH3CHO": 1,
+    #"NSplus": 1,
+    #"H40a": 1,
+    "HC15N": 2,
+    "SO21": 2,
+    "H13CN": 2,
+    "HN13C": 2,
+    "H13COp": 2,
+    #"CS21": 1,
+    #"HC3N": 1,
+    "HCOP": 7,
+    "HCOP_mopra": 3,
+    "SiO21": 2,
+    #"SO32": 1,
+}
+
+
+def main():
+    if os.getenv('MOLNAME'):
+        molnames = [os.getenv('MOLNAME')]
+    else:
+        molnames = dsfactor_dict.keys()
+    if os.getenv('USE_DASK'):
+        use_dask = os.getenv('USE_DASK').lower() == 'true'
+    else:
+        use_dask = False
+
+    if os.getenv('NUM_CORES'):
+        numcores = int(os.getenv('NUM_CORES'))
+    elif os.getenv('SLURM_NTASKS') and os.getenv('SLURM_JOB_CPUS_PER_NODE'):
+        numcores = int(os.getenv('SLURM_NTASKS')) * int(os.getenv('SLURM_JOB_CPUS_PER_NODE'))
+    else:
+        numcores = 1
+
+    print(f"Using Dask: {use_dask}")
+    print(f"Number of cores: {numcores}")
+    print(f"molnames: {molnames}")
+
+    for molname in molnames:
+        print(molname)
+        target_path = '/red/adamginsburg/workdir/mosaics/'
+        cubename = f"{basepath}/mosaics/cubes/{molname}_CubeMosaic.fits"
+        outname = f"{target_path}/mosaics/cubes/{molname}_CubeMosaic_spectrally.fits"
+        downsample_spectrally(cubename,
+                              outname,
+                              factor=dsfactor_dict[molname],
+                              use_dask=use_dask,
+                              num_cores=numcores,
+                              )
+        if os.path.exists(outname):
+            destination = f"{basepath}/mosaics/cubes/{os.path.basename(outname)}"
+            if os.path.exists(destination):
+                print(f"Overwriting destionation {destination}")
+            print(f"Moving {outname} to {destination}")
+            shutil.move(outname, destination)
+
+
+if __name__ == "__main__":
+    main()

aces/analysis/giantcube_cuts.py

@@ -116,6 +116,28 @@ def get_prunemask_space_dask(mask, npix=10):
     return da.stack(masks)
 
 
+def get_noedge_mask(cube, iterations=40):
+    """
+    Erode the edges on a channel-by-channel basis
+
+    Empirically, there are ~35ish pixels along the edge of fields with visibily inflated noise
+    """
+    good = cube.mask.include()
+
+    if hasattr(cube, 'rechunk'):
+        ndi = ndmorph
+    else:
+        ndi = ndimage
+
+    struct = ndi.generate_binary_structure(3, 1)
+    # make cubic binary structure, but ignore the spectral dimension
+    struct[0, :, :] = False
+    struct[-1, :, :] = False
+    good = ndi.binary_erosion(good, structure=struct, iterations=iterations)
+
+    return good
+
+
 def get_pruned_mask(cube, noise, threshold1=1.0, threshold2=5.0):
     beam_area_pix = get_beam_area_pix(cube)
 
@@ -153,7 +175,7 @@ def get_beam_area_pix(cube):
     return beam_area_pix
 
 
-def do_pvs(cube, molname, mompath=f'{basepath}/mosaics/cubes/moments/',
+def do_pvs(cube, molname, mask=None, mompath=f'{basepath}/mosaics/cubes/moments/',
            howargs={}):
 
     t0 = time.time()
@@ -193,20 +215,24 @@ def do_pvs(cube, molname, mompath=f'{basepath}/mosaics/cubes/moments/',
     makepng(data=pv_mean.value, wcs=pv_mean.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_PV_b_mean.png",
             stretch='asinh', min_percent=1, max_percent=99.5)
 
-    # PHANGS-like cuts
-    noise = fits.getdata(f"{mompath}/{molname}_CubeMosaic_noisemap.fits")
-    mcube = cube.with_mask(cube > noise * cube.unit)
+    if mask is None:
+        print("Loading noise map & calculating mask")
+        noise = fits.getdata(f"{mompath}/{molname}_CubeMosaic_noisemap.fits")
+        mcube = cube.with_mask(cube > noise * cube.unit)
+    else:
+        print("Using provided mask")
+        mcube = cube.with_mask(mask)
 
-    print(f"PV mean.  dt={time.time() - t0}")
+    print(f"PV mean with mask.  dt={time.time() - t0}")
     pv_mean_masked = mcube.mean(axis=1, **howargs)
     pv_mean_masked.write(f"{mompath}/{molname}_CubeMosaic_PV_mean_masked.fits", overwrite=True)
-    makepng(data=pv_mean.value, wcs=pv_mean.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_PV_mean_masked.png",
+    makepng(data=pv_mean_masked.value, wcs=pv_mean_masked.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_PV_mean_masked.png",
             stretch='asinh', min_percent=1, max_percent=99.5)
 
-    print(f"PV mean 2.  dt={time.time() - t0}")
+    print(f"PV mean with mask 2.  dt={time.time() - t0}")
     pv_mean_masked = mcube.mean(axis=2, **howargs)
     pv_mean_masked.write(f"{mompath}/{molname}_CubeMosaic_PV_b_mean_masked.fits", overwrite=True)
-    makepng(data=pv_mean.value, wcs=pv_mean.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_PV_b_mean_masked.png",
+    makepng(data=pv_mean_masked.value, wcs=pv_mean_masked.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_PV_b_mean_masked.png",
             stretch='asinh', min_percent=1, max_percent=99.5)
 
     # 2.5 sigma cut
@@ -216,13 +242,13 @@ def do_pvs(cube, molname, mompath=f'{basepath}/mosaics/cubes/moments/',
     print(f"PV mean.  dt={time.time() - t0}")
     pv_mean_masked = mdcube_2p5.mean(axis=1, **howargs)
     pv_mean_masked.write(f"{mompath}/{molname}_CubeMosaic_PV_mean_masked_2p5.fits", overwrite=True)
-    makepng(data=pv_mean.value, wcs=pv_mean.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_PV_mean_masked_2p5.png",
+    makepng(data=pv_mean_masked.value, wcs=pv_mean_masked.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_pv_mean_masked_2p5.png",
             stretch='asinh', min_percent=1, max_percent=99.5)
 
     print(f"PV mean 2.  dt={time.time() - t0}")
     pv_mean_masked = mdcube_2p5.mean(axis=2, **howargs)
     pv_mean_masked.write(f"{mompath}/{molname}_CubeMosaic_PV_b_mean_masked_2p5.fits", overwrite=True)
-    makepng(data=pv_mean.value, wcs=pv_mean.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_PV_b_mean_masked_2p5.png",
+    makepng(data=pv_mean_masked.value, wcs=pv_mean_masked.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_PV_b_mean_masked_2p5.png",
             stretch='asinh', min_percent=1, max_percent=99.5)
 
 
@@ -383,6 +409,32 @@ def do_all_stats(cube, molname, mompath=f'{basepath}/mosaics/cubes/moments/',
             stretch='asinh', vmin=-0.1, max_percent=99.5)
     del mom0
 
+    print(f"Computing no-edge mask. dt={time.time() - t0}")
+    noedge = get_noedge_mask(cube)
+
+    print(f"no-edge max.  dt={time.time() - t0}", flush=True)
+    nemx = cube.with_mask(noedge).max(axis=0, **howargs)
+    print(f"Done computing no-edge max, now writing. dt={time.time() - t0}")
+    nemx.write(f"{mompath}/{molname}_CubeMosaic_edgelessmax.fits", overwrite=True)
+    print(f"Done writing no-edge max, now pnging. dt={time.time() - t0}")
+    makepng(data=nemx.value, wcs=nemx.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_edgelessmax.png",
+            stretch='asinh', min_percent=0.1, max_percent=99.9)
+
+    print(f"No-edge Noisemap. dt={time.time() - t0}")
+
+    if hasattr(cube, 'rechunk'):
+        print(f"no-edge madstd (dask).  dt={time.time() - t0}")
+        noise_madstd = cube.with_mask(noedge).mad_std(axis=0)
+    else:
+        print(f"no-edge madstd (numpy).  dt={time.time() - t0}")
+        noise_madstd = cube.with_mask(noedge).mad_std(axis=0, how='ray', progressbar=True)
+    print(f"Done with no-edge madstd. Writing noisemap to {mompath}/{molname}_CubeMosaic_madstd.fits.  dt={time.time() - t0}")
+    noise_madstd.write(f"{mompath}/{molname}_CubeMosaic_edgelessmadstd.fits", overwrite=True)
+    makepng(data=noise_madstd.value, wcs=noise_madstd.wcs, imfn=f"{mompath}/{molname}_CubeMosaic_edgelessmadstd.png",
+            stretch='asinh', min_percent=0.5, max_percent=99.5)
+
+    return BooleanArrayMask(signal_mask_both, cube.wcs)
+
 
 def main():
     dodask = os.getenv('USE_DASK')
@@ -425,6 +477,8 @@ def main():
         t0 = time.time()
 
         cube = SpectralCube.read(cubefilename)
+        if cube.shape[0] > 500 and 'H40a' not in cubefilename:
+            cube = SpectralCube.read(cubefilename.replace(".fits", "_spectrally.fits"))
 
         howargs = {'how': 'slice', 'progressbar': True}
         print(f"Non-Dask: how are we computing? {howargs}")
@@ -448,6 +502,8 @@ def main():
         print(f"Memory limit: {memory_limit}")
 
         cube = SpectralCube.read(cubefilename, use_dask=True)
+        if cube.shape[0] > 500 and 'H40a' not in cubefilename:
+            cube = SpectralCube.read(cubefilename.replace(".fits", "_spectrally.fits"), use_dask=True)
 
         t0 = time.time()
 
@@ -499,10 +555,10 @@ def main():
             print(f"Dask client number of workers: {len(client.scheduler_info()['workers'])}")
             print(f"Using scheduler {scheduler}", flush=True)
 
-    do_all_stats(cube, molname=molname, howargs=howargs)
+    signal_mask_both = do_all_stats(cube, molname=molname, howargs=howargs)
 
     if do_pv:
-        do_pvs(cube, molname=molname, howargs=howargs)
+        do_pvs(cube, molname=molname, howargs=howargs, mask=signal_mask_both)
 
     if dods:
         print(f"Downsampling.  dt={time.time() - t0}")

aces/analysis/latex_table.py

@@ -228,6 +228,7 @@ def make_observation_table(access_token=None):
     metadata['Observation Start Time'].format = 'fits'
     metadata['Observation End Time'] = Time(metadata['t_max'], format='mjd')
     metadata['Observation End Time'].format = 'fits'
+    metadata = metadata[metadata['target_name'] == 'Sgr_A_star']
 
     sub_meta = metadata['schedblock_name', 'Observation Start Time',
                         'Observation End Time', 'pwv', 't_exptime',
@@ -279,26 +280,34 @@ def make_observation_table(access_token=None):
                 path = f'{basepath}/data/2021.1.00172.L/science_goal.uid___A001_X1590_X30a8/group.uid___A001_X1590_X30a9/member.{mousstr}/calibrated/working/{xidstr}.ms'
                 try:
                     tbl = Table.read(path + "/ASDM_CALWVR")
-                    these_pwvs.append(np.median(tbl['water']) * 1000)
+                    med_pwv = np.median(tbl['water'])
+                    if np.isnan(med_pwv):
+                        raise ValueError("WTF?")
+                    these_pwvs.append(med_pwv * 1000)
                 except ValueError:
                     lltbl = casa_formats_io.table_reader.CASATable.read(path + "/ASDM_CALWVR")
                     tbl = lltbl.as_astropy_table(include_columns=['water'])
-                    these_pwvs.append(np.median(tbl['water']) * 1000)
+                    med_pwv = np.median(tbl['water'])
+                    if np.isnan(med_pwv):
+                        raise ValueError("WTF?")
+                    these_pwvs.append(med_pwv * 1000)
                 except IOError as ex:
                     if 'TP' not in sbname:
                         path = f'{basepath}/data/2021.1.00172.L/science_goal.uid___A001_X1590_X30a8/group.uid___A001_X1590_X30a9/member.{mousstr}/calibrated/{xidstr}.ms'
                         if os.path.exists(path):
                             lltbl = casa_formats_io.table_reader.CASATable.read(path + "/ASDM_CALWVR")
                             tbl = lltbl.as_astropy_table(include_columns=['water'])
-                            these_pwvs.append(np.median(tbl['water']) * 1000)
+                            med_pwv = np.median(tbl['water'])
+                            if np.isnan(med_pwv):
+                                raise ValueError("WTF?")
+                            these_pwvs.append(med_pwv * 1000)
                         else:
                             these_pwvs.append(np.nan)
-                            if 'TM1' in sbname:
-                                print(sbname, mous, xid, ex)
+                            print("No CALWVR table found for ", sbname, mous, xid, ex)
+                    print(f"Skipped TP {path}")
                 except Exception as ex:
                     these_pwvs.append(np.nan)
-                    if 'TM1' in sbname:
-                        print(sbname, mous, xid, ex)
+                    print("Unknown exception for ", sbname, mous, xid, ex)
             pwv_measurements[sbname] = these_pwvs
         pwvs.append(fr"\makecell{{{',\\\\ '.join(map(pwvfmt, these_pwvs))}}}".replace('nan', '-'))
     with open(pwv_cache_fn, 'w') as fh:
@@ -376,9 +385,11 @@ def try_read_table(index):
     usub_meta['Time'].unit = u.min
     usub_meta['PWV'].unit = u.mm
 
+    colnames_noconfig = colnames[:2] + colnames[3:]
+
     usub_meta[colnames][tm].write(f'{basepath}/tables/observation_metadata_12m.ecsv', format='ascii.ecsv', overwrite=True)
-    usub_meta[colnames][sm].write(f'{basepath}/tables/observation_metadata_7m.ecsv', format='ascii.ecsv', overwrite=True)
-    usub_meta[colnames][tp].write(f'{basepath}/tables/observation_metadata_TP.ecsv', format='ascii.ecsv', overwrite=True)
+    usub_meta[colnames_noconfig][sm].write(f'{basepath}/tables/observation_metadata_7m.ecsv', format='ascii.ecsv', overwrite=True)
+    usub_meta[colnames_noconfig][tp].write(f'{basepath}/tables/observation_metadata_TP.ecsv', format='ascii.ecsv', overwrite=True)
 
     usub_meta['Execution Dates'] = [fr'\makecell{{{",\\\\ ".join([x for x in ed.split(", ")])}.}}' for ed in usub_meta['Execution Dates']]
 
@@ -477,12 +488,13 @@ def make_spw_table():
     latexdict['tabletype'] = 'table*'
     latexdict['tablefoot'] = (
         "}\\par\n"
-        "ACES Spectral Configuration, including a non-exhaustive lists of prominent, "
+        "ACES Spectral Configuration, including a non-exhaustive list of prominent, "
         "potentially continuum-affecting, lines.  The included lines are those that are, "
-        "in at least some portion of the survey, masked out (see Section \\ref{sec:continuum_selection})."
+        "in at least some portion of the survey, masked out (see Section \\ref{sec:continuum_selection}).  "
+        "The rest frequencies of the targeted lines are given in GHz in the row below their names."
     )
 
-    ftbl.write(f"{basepath}/papers/continuum_data/spectral_setup.tex", formats=formats,
+    ftbl.write(f"{basepath}/papers/continuum_data/tables/spectral_setup.tex", formats=formats,
                overwrite=True, latexdict=latexdict)
 
     return ftbl

aces/data/regions/zoom_regions.reg

@@ -2,7 +2,7 @@
 global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1
 galactic
 box(0.096993759, 0.024355304, 504.8203", 358.8801", 0) # color=#2EE6D6 width=2 text={North Arches} font="helvetica 10 normal roman"
-box(0.163727690, -0.070544532, 326.2236", 284.0760", 0) # color=#2EE6D6 width=2 text={Pistol Quintuplet} font="helvetica 10 normal roman"
+box(0.163727690, -0.070544532, 326.2236", 284.0760", 0) # color=#2EE6D6 width=2 text={Pistol Sickle Quintuplet} font="helvetica 10 normal roman"
 box(0.514917217, -0.052354790, 498.0334", 412.1990", 0) # color=#2EE6D6 width=2 text={Sgr B1} font="helvetica 10 normal roman"
 box(359.442350761, -0.096457228, 208.4775", 166.0009", 0) # color=#2EE6D6 width=2 text={Sgr C} font="helvetica 10 normal roman"
 box(359.974945550, -0.078751871, 98.2085", 82.4607", 0) # color=#2EE6D6 width=2 text={Sgr A East} font="helvetica 10 normal roman"
@@ -13,7 +13,7 @@ box(359.611683297, -0.239006506, 131.4033", 107.9750", 0) # color=#2EE6D6 width=
 box(0.254526363, 0.019530785, 393.3003", 297.6328", 0) # color=#2EE6D6 width=2 text={Brick} font="helvetica 10 normal roman"
 box(0.375781645, 0.017482982, 451.7147", 390.1172", 0) # color=#2EE6D6 width=2 text={Dust Ridge C} font="helvetica 10 normal roman"
 box(0.675977373, -0.045253042, 660.7727", 696.2378", 0) # color=#2EE6D6 width=2 text={Sgr B2} font="helvetica 10 normal roman"
-box(0.815018274, -0.190889216, 382.5156", 396.7885", 0) # color=#2EE6D6 width=2 text={Nonhebel Bubble} font="helvetica 10 normal roman"
+box(0.815018274, -0.190889216, 382.5156", 396.7885", 0) # color=#2EE6D6 width=2 text={M0.8 Ring} font="helvetica 10 normal roman"
 box(359.944870512, -0.046021049, 178.8130", 162.1015", 0) # color=#2EE6D6 width=2 text={Sgr A*} font="helvetica 10 normal roman"
-box(0.063474920, -0.079908407, 370.9951", 342.5856", 0) # color=#2EE6D6 width=2 text={TLP} font="helvetica 10 normal roman"
+box(0.063474920, -0.079908407, 370.9951", 342.5856", 0) # color=#2EE6D6 width=2 text={South Dust Ridge} font="helvetica 10 normal roman"
 box(359.632471104, -0.054839112, 279.0820", 269.0550", 0) # color=#2EE6D6 width=2 text={G359.63} font="helvetica 10 normal roman"

aces/hipergator_scripts/run_allgiantcube_analysis.sh

@@ -35,7 +35,15 @@ else
 fi
 
 #for MOLNAME in SO21 H13CN HN13C H40a CH3CHO NSplus; do # H13COp CS21 HC3N HCOP SiO21 HNCO_7m12mTP; do
-for MOLNAME in CH3CHO NSplus H40a HC15N SO21 H13CN HN13C H13COp CS21 HC3N HCOP SiO21 SO32 HNCO_7m12mTP; do
+for MOLNAME in HNCO_7m12mTP HCOP_noTP CH3CHO NSplus H40a HC15N SO21 H13CN HN13C H13COp CS21 HC3N HCOP SiO21 SO32; do
+
+    if [[ $MOLNAME == *"HNCO"* ]]; then
+        export mem=512
+    elif [[ $MOLNAME == *"HCOP"* ]]; then
+        export mem=512
+    else
+        export mem=256
+    fi
 
     if [ -e /orange/adamginsburg/ACES/mosaics/cubes/${MOLNAME}_CubeMosaic.fits ]; then
 
@@ -45,18 +53,21 @@ for MOLNAME in CH3CHO NSplus H40a HC15N SO21 H13CN HN13C H13COp CS21 HC3N HCOP S
         export MOLNAME
         # optional
         export DOWNSAMPLE=False
-        export DO_PV=False
+        export DO_PV=True
 
         echo "Giant ${MOLNAME} cube"
         #/red/adamginsburg/miniconda3/envs/python312/bin/python /orange/adamginsburg/ACES/reduction_ACES/aces/analysis/giantcube_cuts.py || exit 1
         sbatch --job-name=aces_giantcube_analysis_${MOLNAME}${DASK}${DASK_CLIENT} \
                --output=/red/adamginsburg/ACES/logs/aces_giantcube_analysis_${MOLNAME}${DASK}${DASK_CLIENT}_%j.log \
                --account=astronomy-dept --qos=astronomy-dept-b --ntasks=32 --nodes=1 \
-               --mem=256gb --time=96:00:00 \
+               --mem=${mem}gb --time=96:00:00 \
                --export=MOLNAME=${MOLNAME},USE_DASK=${USE_DASK},USE_LOCAL=${USE_LOCAL},DASK_CLIENT=${DASK_CLIENT},DOWNSAMPLE=${DOWNSAMPLE},DO_PV=${DO_PV} \
                --wrap /red/adamginsburg/miniconda3/envs/python312/bin/aces_giantcube_analysis
     else
         echo "${MOLNAME}_CubeMosaic.fits does not exist"
         ls -lh ${MOLNAME}*fits
     fi
 done
+
+# one=line version
+# export MOLNAME=SiO21 DOWNSAMPLE=True DOPV=True mem=256 USE_LOCAL=True USE_DASK=False;            sbatch --job-name=aces_giantcube_analysis_${MOLNAME}${DASK}${DASK_CLIENT}                --output=/red/adamginsburg/ACES/logs/aces_giantcube_analysis_${MOLNAME}${DASK}${DASK_CLIENT}_%j.log                --account=astronomy-dept --qos=astronomy-dept-b --ntasks=32 --nodes=1                --mem=${mem}gb --time=96:00:00                --export=MOLNAME=${MOLNAME},USE_DASK=${USE_DASK},USE_LOCAL=${USE_LOCAL},DASK_CLIENT=${DASK_CLIENT},DOWNSAMPLE=${DOWNSAMPLE},DO_PV=${DO_PV}                --wrap /red/adamginsburg/miniconda3/envs/python312/bin/aces_giantcube_analysis