From 4f7477db801bdc265a05bbfa889195a18473760d Mon Sep 17 00:00:00 2001 From: Tony Wong Date: Sun, 12 Nov 2023 19:56:17 -0600 Subject: [PATCH] Allow SNR images to be generated when an external mask is supplied. --- README.md | 6 ++++-- maskmoment/maskmoment.py | 42 ++++++++++++++++++++-------------------- 2 files changed, 25 insertions(+), 23 deletions(-) diff --git a/README.md b/README.md index 37b57cc..96c7eb6 100644 --- a/README.md +++ b/README.md @@ -114,12 +114,14 @@ How to use: Output the moment maps in K units if the cube is in Jy/beam units. Default: True vmin : float or :class:`~astropy.units.Quantity`, optional - Minimum channel number or velocity to use for moment calculation. + Minimum channel number or velocity to use for all calculations. Note + that channels are discarded before any rms estimation using edgech. If given as astropy quantity, should be provided in velocity units. If NOT given as astropy quantity, interpreted as channel number. Default: Do not impose a velocity cut. vmax : float or :class:`~astropy.units.Quantity`, optional - Maximum channel number or velocity to use for moment calculation. + Maximum channel number or velocity to use for all calculations. Note + that channels are discarded before any rms estimation using edgech. If given as astropy quantity, should be provided in velocity units. If NOT given as astropy quantity, interpreted as channel number. Default: Do not impose a velocity cut. diff --git a/maskmoment/maskmoment.py b/maskmoment/maskmoment.py index c8fccc3..da90eda 100644 --- a/maskmoment/maskmoment.py +++ b/maskmoment/maskmoment.py @@ -33,7 +33,7 @@ def maskmoment(img_fits, gain_fits=None, rms_fits=None, mask_fits=None, outdir=' This should have the same dimensions and units as the image cube. If it is a 2D array it will be replicated along the velocity axis. NOTE: If rms_fits is not given, a noise cube is generated from the - image cube, after removing any gain variation using the gain cube. + image cube, taking into account any gain variation using the gain cube. rms : float, optional Global estimate of the rms noise, to be used instead of trying to estimate it from the data. It should have the units of the input cube. @@ -264,29 +264,29 @@ def maskmoment(img_fits, gain_fits=None, rms_fits=None, mask_fits=None, outdir=' # # --- GENERATE AND OUTPUT SNR CUBE, PEAK SNR IMAGE # + snr_cube = image_cube / rms_cube + print('SNR cube:\n',snr_cube) + if output_snr_cube: + hd3d['datamin'] = snr_cube.min().value + hd3d['datamax'] = snr_cube.max().value + hd3d['bunit'] = ' ' + fits.writeto(pth+basename+'.snrcube.fits.gz', snr_cube._data.astype(np.float32), + hd3d, overwrite=True) + print('Wrote', pth+basename+'.snrcube.fits.gz') + if output_snr_peak: + snr_peak = snr_cube.max(axis=0) + hd2d['datamin'] = np.nanmin(snr_peak.value) + hd2d['datamax'] = np.nanmax(snr_peak.value) + hd2d['bunit'] = ' ' + fits.writeto(pth+basename+'.snrpk.fits.gz', snr_peak.astype(np.float32), + hd2d, overwrite=True) + print('Wrote', pth+basename+'.snrpk.fits.gz') + # + # --- GENERATE AND OUTPUT DILATED MASK + # if mask_fits is not None: dilatedmask = fits.getdata(mask_fits) else: - snr_cube = image_cube / rms_cube - print('SNR cube:\n',snr_cube) - if output_snr_cube: - hd3d['datamin'] = snr_cube.min().value - hd3d['datamax'] = snr_cube.max().value - hd3d['bunit'] = ' ' - fits.writeto(pth+basename+'.snrcube.fits.gz', snr_cube._data.astype(np.float32), - hd3d, overwrite=True) - print('Wrote', pth+basename+'.snrcube.fits.gz') - if output_snr_peak: - snr_peak = snr_cube.max(axis=0) - hd2d['datamin'] = np.nanmin(snr_peak.value) - hd2d['datamax'] = np.nanmax(snr_peak.value) - hd2d['bunit'] = ' ' - fits.writeto(pth+basename+'.snrpk.fits.gz', snr_peak.astype(np.float32), - hd2d, overwrite=True) - print('Wrote', pth+basename+'.snrpk.fits.gz') - # - # --- GENERATE AND OUTPUT DILATED MASK - # if fwhm is not None or vsm is not None: sm_snrcube = smcube(snr_cube, fwhm=fwhm, vsm=vsm, vsm_type=vsm_type, edgech=edgech, huge=huge_operations)