diff --git a/aces/analysis/downsample_all_spectrally.py b/aces/analysis/downsample_all_spectrally.py new file mode 100644 index 00000000..b685b2f1 --- /dev/null +++ b/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() diff --git a/aces/analysis/giantcube_cuts.py b/aces/analysis/giantcube_cuts.py index debf12a2..00891163 100644 --- a/aces/analysis/giantcube_cuts.py +++ b/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}") diff --git a/aces/analysis/latex_table.py b/aces/analysis/latex_table.py index 9b929c1c..1de71386 100644 --- a/aces/analysis/latex_table.py +++ b/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 diff --git a/aces/data/regions/zoom_regions.reg b/aces/data/regions/zoom_regions.reg index 605feb11..c7c45709 100644 --- a/aces/data/regions/zoom_regions.reg +++ b/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" diff --git a/aces/hipergator_scripts/run_allgiantcube_analysis.sh b/aces/hipergator_scripts/run_allgiantcube_analysis.sh index bbe2d755..f4b62798 100644 --- a/aces/hipergator_scripts/run_allgiantcube_analysis.sh +++ b/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,14 +53,14 @@ 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 @@ -60,3 +68,6 @@ for MOLNAME in CH3CHO NSplus H40a HC15N SO21 H13CN HN13C H13COp CS21 HC3N HCOP S 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 \ No newline at end of file diff --git a/aces/hipergator_scripts/run_downsample_spectrally.sh b/aces/hipergator_scripts/run_downsample_spectrally.sh new file mode 100644 index 00000000..213a5dd5 --- /dev/null +++ b/aces/hipergator_scripts/run_downsample_spectrally.sh @@ -0,0 +1,39 @@ +date + +. ~/.gh_token +echo $GITHUB_TOKEN + +cd /red/adamginsburg/ACES/workdir/ +pwd + + +#for MOLNAME in SO21 H13CN HN13C H40a CH3CHO NSplus; do # H13COp CS21 HC3N HCOP SiO21 HNCO_7m12mTP; do +for MOLNAME in HNCO_7m12mTP HCOP_noTP HC15N SO21 H13CN HN13C H13COp HCOP SiO21; do + + if [ -e /orange/adamginsburg/ACES/mosaics/cubes/${MOLNAME}_CubeMosaic.fits ]; then + + #echo "test import" + #/orange/adamginsburg/miniconda3/envs/python312/bin/python -c "import zipfile" || exit 1 + + export MOLNAME + # optional + export DOWNSAMPLE=False + export DO_PV=True + + echo "Giant ${MOLNAME} cube - spectrally downsampling" + #/red/adamginsburg/miniconda3/envs/python312/bin/python /orange/adamginsburg/ACES/reduction_ACES/aces/analysis/giantcube_cuts.py || exit 1 + sbatch --job-name=aces_downsample_spectrally_${MOLNAME} \ + --output=/red/adamginsburg/ACES/logs/aces_downsample_spectrally_${MOLNAME}_%j.log \ + --account=astronomy-dept --qos=astronomy-dept-b --ntasks=32 --nodes=1 \ + --cpus-per-task=1 \ + --mem=256gb --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/python -c "from aces.analysis import downsample_all_spectrally; downsample_all_spectrally.main()"' + 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 \ No newline at end of file diff --git a/aces/hipergator_scripts/slurm_hcop_mopra_mosaic_jobs.sh b/aces/hipergator_scripts/slurm_hcop_mopra_mosaic_jobs.sh new file mode 100644 index 00000000..13977b0b --- /dev/null +++ b/aces/hipergator_scripts/slurm_hcop_mopra_mosaic_jobs.sh @@ -0,0 +1,16 @@ + +jobid=$(sbatch --job-name=aces_hcop_mopra_mos_arr \ + --output=/red/adamginsburg/ACES/logs/aces_hcop_mopra_mosaic_%j_%A_%a.log \ + --array=0-54 \ + --account=astronomy-dept --qos=astronomy-dept-b \ + --ntasks=8 --nodes=1 --mem=64gb --time=96:00:00 --parsable \ + --wrap "/red/adamginsburg/miniconda3/envs/python310/bin/python -c \"from aces.imaging.mosaic_12m import make_giant_mosaic_cube_hcop_mopra; make_giant_mosaic_cube_hcop_mopra(channels='slurm', skip_final_combination=True, verbose=True,)\"") + +echo "Job IDs are ${jobid}" + +sbatch --job-name=aces_hcop_mopra_mosaic_merge \ + --output=/red/adamginsburg/ACES/logs/aces_hcop_mopra_mosaic_merge_%j.log \ + --account=astronomy-dept --qos=astronomy-dept-b \ + --dependency=afterok:$jobid \ + --ntasks=16 --nodes=1 --mem=64gb --time=96:00:00 \ + --wrap "/red/adamginsburg/miniconda3/envs/python310/bin/python -c \"from aces.imaging.mosaic_12m import make_giant_mosaic_cube_hcop_mopra; make_giant_mosaic_cube_hcop_mopra(channels='all', skip_channel_mosaicing=True, verbose=True,)\"" diff --git a/aces/imaging/make_mosaic.py b/aces/imaging/make_mosaic.py index 84705d63..5256dc12 100644 --- a/aces/imaging/make_mosaic.py +++ b/aces/imaging/make_mosaic.py @@ -933,6 +933,7 @@ def combine_channels_into_mosaic_cube(header, cubename, nchan, channels, pbar = tqdm(channels, desc='Channels') else: pbar = channels + status = {'filled': [], 'skipped': []} for chan in pbar: chanfn = f'{channelmosaic_directory}/{cubename}_CubeMosaic_channel{chan}.fits' if os.path.exists(chanfn): @@ -942,10 +943,13 @@ def combine_channels_into_mosaic_cube(header, cubename, nchan, channels, if verbose: pbar.set_description('Channels (flushing)') hdu.flush() + status['filled'].append(chan) else: pbar.set_description("Channels (skip)") + status['skipped'].append(chan) if verbose: + print(f"Status of filling: {status}") print(f"Moving {output_working_file} to {output_file}") shutil.move(output_working_file, output_file) assert os.path.exists(output_file), f"Failed to move {output_working_file} to {output_file}" @@ -1030,6 +1034,38 @@ def make_downsampled_cube(cubename, outcubename, factor=9, overwrite=True, dscube_s.write(outcubename.replace(".fits", "_spectrally.fits"), overwrite=True) +def downsample_spectrally(cubename, outcubename, factor=9, overwrite=True, + smooth=True, num_cores=1, + verbose=True, + use_dask=True, spectrally_too=True): + assert outcubename.endswith('.fits') + from astropy.convolution import Gaussian1DKernel + + cube = SpectralCube.read(cubename, use_dask=use_dask) + + kernel = Gaussian1DKernel(factor / 2) + print(f"Downsampling cube {cubename} -> {outcubename}") + print(cube) + + if use_dask: + import dask.array as da + from dask.diagnostics import ProgressBar as DaskProgressBar + + with DaskProgressBar(): + # smooth with sigma width half the downsampling + dscube = cube.spectral_smooth(kernel) + + dscube_s = dscube.downsample_axis(factor=factor, axis=0) + dscube_s.write(outcubename, overwrite=True) + else: + cube.allow_huge_operations = True + + dscube = cube.spectral_smooth(kernel, use_memmap=True, num_cores=num_cores, verbose=verbose) + print("Done smoothing, now downsampling") + dscube_s = dscube.downsample_axis(factor=factor, axis=0) + dscube_s.write(outcubename, overwrite=True) + + def rms_map(img, kernel=Gaussian2DKernel(10)): """ Gaussian2DKernel should be larger than the beam, probably at least 2x larger diff --git a/aces/imaging/mosaic_12m.py b/aces/imaging/mosaic_12m.py index d7aa9db4..4cd83630 100644 --- a/aces/imaging/mosaic_12m.py +++ b/aces/imaging/mosaic_12m.py @@ -16,7 +16,8 @@ get_m0, make_downsampled_cube, make_giant_mosaic_cube, - rms + rms, + downsample_spectrally ) from aces.imaging.make_mosaic import make_mosaic as make_mosaic_, all_lines as all_lines_ # import os @@ -707,6 +708,10 @@ def make_giant_mosaic_cube_sio21(**kwargs): make_downsampled_cube(f'{basepath}/mosaics/cubes/SiO21_CubeMosaic.fits', f'{basepath}/mosaics/cubes/SiO21_CubeMosaic_downsampled9.fits', ) + downsample_spectrally(f'{basepath}/mosaics/cubes/SiO21_CubeMosaic.fits', + f'{basepath}/mosaics/cubes/SiO21_CubeMosaic_spectrally.fits', + factor=2, + ) def make_giant_mosaic_cube_hnco(**kwargs): @@ -745,6 +750,10 @@ def make_giant_mosaic_cube_hnco(**kwargs): make_downsampled_cube(f'{basepath}/mosaics/cubes/HNCO_CubeMosaic.fits', f'{basepath}/mosaics/cubes/HNCO_CubeMosaic_downsampled9.fits', ) + downsample_spectrally(f'{basepath}/mosaics/cubes/HNCO_CubeMosaic.fits', + f'{basepath}/mosaics/cubes/HNCO_CubeMosaic_spectrally.fits', + factor=7, + ) def make_giant_mosaic_cube_nsplus(**kwargs): @@ -892,6 +901,10 @@ def make_giant_mosaic_cube_hnco_TP7m12m(**kwargs): make_downsampled_cube(f'{basepath}/mosaics/cubes/HNCO_7m12mTP_CubeMosaic.fits', f'{basepath}/mosaics/cubes/HNCO_7m12mTP_CubeMosaic_downsampled9.fits', ) + downsample_spectrally(f'{basepath}/mosaics/cubes/HNCO_7m12mTP_CubeMosaic.fits', + f'{basepath}/mosaics/cubes/HNCO_7m12mTP_CubeMosaic_spectrally.fits', + factor=7, + ) def make_giant_mosaic_cube_hcop(**kwargs): @@ -929,6 +942,10 @@ def make_giant_mosaic_cube_hcop(**kwargs): make_downsampled_cube(f'{basepath}/mosaics/cubes/HCOP_CubeMosaic.fits', f'{basepath}/mosaics/cubes/HCOP_CubeMosaic_downsampled9.fits', ) + downsample_spectrally(f'{basepath}/mosaics/cubes/HCOP_CubeMosaic.fits', + f'{basepath}/mosaics/cubes/HCOP_CubeMosaic_spectrally.fits', + factor=7, + ) def make_giant_mosaic_cube_hnco_TP7m(**kwargs): @@ -968,6 +985,10 @@ def make_giant_mosaic_cube_hnco_TP7m(**kwargs): f'{basepath}/mosaics/cubes/HNCO_7mTP_CubeMosaic_downsampled9.fits', overwrite=True ) + downsample_spectrally(f'{basepath}/mosaics/cubes/HNCO_7mTP_CubeMosaic.fits', + f'{basepath}/mosaics/cubes/HNCO_7mTP_CubeMosaic_spectrally.fits', + factor=7, + ) def make_giant_mosaic_cube_ch3cho(**kwargs): @@ -1070,6 +1091,10 @@ def make_giant_mosaic_cube_h13cn(**kwargs): make_downsampled_cube(f'{basepath}/mosaics/cubes/H13CN_CubeMosaic.fits', f'{basepath}/mosaics/cubes/H13CN_CubeMosaic_downsampled9.fits', ) + downsample_spectrally(f'{basepath}/mosaics/cubes/H13CN_CubeMosaic.fits', + f'{basepath}/mosaics/cubes/H13CN_CubeMosaic_spectrally.fits', + factor=2, + ) def make_giant_mosaic_cube_h13cop(**kwargs): @@ -1102,6 +1127,10 @@ def make_giant_mosaic_cube_h13cop(**kwargs): make_downsampled_cube(f'{basepath}/mosaics/cubes/H13COp_CubeMosaic.fits', f'{basepath}/mosaics/cubes/H13COp_CubeMosaic_downsampled9.fits', ) + downsample_spectrally(f'{basepath}/mosaics/cubes/H13COp_CubeMosaic.fits', + f'{basepath}/mosaics/cubes/H13COp_CubeMosaic_spectrally.fits', + factor=2, + ) def make_giant_mosaic_cube_hn13c(**kwargs): @@ -1135,6 +1164,10 @@ def make_giant_mosaic_cube_hn13c(**kwargs): make_downsampled_cube(f'{basepath}/mosaics/cubes/HN13C_CubeMosaic.fits', f'{basepath}/mosaics/cubes/HN13C_CubeMosaic_downsampled9.fits', ) + downsample_spectrally(f'{basepath}/mosaics/cubes/HN13C_CubeMosaic.fits', + f'{basepath}/mosaics/cubes/HN13C_CubeMosaic_spectrally.fits', + factor=2, + ) def make_giant_mosaic_cube_so21(**kwargs): @@ -1169,6 +1202,10 @@ def make_giant_mosaic_cube_so21(**kwargs): make_downsampled_cube(f'{basepath}/mosaics/cubes/SO21_CubeMosaic.fits', f'{basepath}/mosaics/cubes/SO21_CubeMosaic_downsampled9.fits', ) + downsample_spectrally(f'{basepath}/mosaics/cubes/SO21_CubeMosaic.fits', + f'{basepath}/mosaics/cubes/SO21_CubeMosaic_spectrally.fits', + factor=2, + ) def make_giant_mosaic_cube_hc15n(**kwargs): @@ -1201,6 +1238,10 @@ def make_giant_mosaic_cube_hc15n(**kwargs): make_downsampled_cube(f'{basepath}/mosaics/cubes/HC15N_CubeMosaic.fits', f'{basepath}/mosaics/cubes/HC15N_CubeMosaic_downsampled9.fits', ) + downsample_spectrally(f'{basepath}/mosaics/cubes/HC15N_CubeMosaic.fits', + f'{basepath}/mosaics/cubes/HC15N_CubeMosaic_spectrally.fits', + factor=2, + ) def make_giant_mosaic_cube_h40a(**kwargs): @@ -1241,7 +1282,7 @@ def make_giant_mosaic_cube_hcop_noTP(**kwargs): #filelist += glob.glob(f'{basepath}/rawdata/2021.1.00172.L/s*/g*/m*/calibrated/working/*spw29.cube.I.manual*image.pbcor.statcont.contsub.fits') #filelist += glob.glob(f'{basepath}/rawdata/2021.1.00172.L/s*/g*/m*/calibrated/working/*sci29.cube.I.manual*image.pbcor.statcont.contsub.fits') #filelist += glob.glob(f'{basepath}/rawdata/2021.1.00172.L/s*/g*/m*/manual/*29.cube.I.manual.pbcor.fits') - filelist = glob.glob('/orange/adamginsburg/ACES/upload//upload/HCOp_feather_images/TM_7M_only_feather/Sgr_A_st_*.7M_12M_feather.SPW_29.image.statcont.contsub.fits') + filelist = glob.glob('/orange/adamginsburg/ACES/upload/HCOp_feather_images/MOPRA_12M_7M_featherTM_7M_only_feather/Sgr_A_st_*.7M_12M_feather.SPW_29.image.statcont.contsub.fits') print(f"Found {len(filelist)} HCOP-containing spw29 files") @@ -1269,4 +1310,49 @@ def make_giant_mosaic_cube_hcop_noTP(**kwargs): if not kwargs.get('skip_final_combination') and not kwargs.get('test'): make_downsampled_cube(f'{basepath}/mosaics/cubes/HCOP_noTP_CubeMosaic.fits', f'{basepath}/mosaics/cubes/HCOP_noTP_CubeMosaic_downsampled9.fits', + ) + downsample_spectrally(f'{basepath}/mosaics/cubes/HCOP_noTP_CubeMosaic.fits', + f'{basepath}/mosaics/cubes/HCOP_noTP_CubeMosaic_spectrally.fits', + factor=7, + ) + + +def make_giant_mosaic_cube_hcop_mopra(**kwargs): + + #filelist = glob.glob(f'{basepath}/rawdata/2021.1.00172.L/s*/g*/m*/calibrated/working/*spw29.cube.I.iter1.image.pbcor.statcont.contsub.fits') + #filelist += glob.glob(f'{basepath}/rawdata/2021.1.00172.L/s*/g*/m*/calibrated/working/*spw29.cube.I.manual*image.pbcor.statcont.contsub.fits') + #filelist += glob.glob(f'{basepath}/rawdata/2021.1.00172.L/s*/g*/m*/calibrated/working/*sci29.cube.I.manual*image.pbcor.statcont.contsub.fits') + #filelist += glob.glob(f'{basepath}/rawdata/2021.1.00172.L/s*/g*/m*/manual/*29.cube.I.manual.pbcor.fits') + filelist = glob.glob('/orange/adamginsburg/ACES/upload/HCOp_feather_images/MOPRA_12M_7M_feather/Sgr_A_st_*.SD_7M_12M_feather.hco*.image.statcont.contsub.fits') + + print(f"Found {len(filelist)} HCOP-containing spw29 files (MOPRA)") + + check_files(filelist) + + weightfilelist = [get_weightfile(fn, spw=29) for fn in filelist] + for fn in weightfilelist: + assert os.path.exists(fn) + + restfrq = 89.188526e9 + cdelt_kms = 0.5401059211273 + make_giant_mosaic_cube(filelist, + weightfilelist=weightfilelist, + reference_frequency=restfrq, + cdelt_kms=cdelt_kms, + cubename='HCOP_mopra', + nchan=550, + beam_threshold=3.3 * u.arcsec, + channelmosaic_directory=f'{basepath}/mosaics/HCOP_Channels/', + fail_if_cube_dropped=False, + #use_reproject_cube=True, + parallel=os.getenv('SLURM_NTASKS'), + **kwargs,) + + if not kwargs.get('skip_final_combination') and not kwargs.get('test'): + make_downsampled_cube(f'{basepath}/mosaics/cubes/HCOP_mopra_CubeMosaic.fits', + f'{basepath}/mosaics/cubes/HCOP_mopra_CubeMosaic_downsampled9.fits', + ) + downsample_spectrally(f'{basepath}/mosaics/cubes/HCOP_mopra_CubeMosaic.fits', + f'{basepath}/mosaics/cubes/HCOP_mopra_CubeMosaic_spectrally.fits', + factor=3, ) \ No newline at end of file