Loading etl/mdl/mdl_extract_ch3oh_spectra.py +8 −7 Changes for etl/mdl/mdl_extract_ch3oh_spectra.py: 8 added lines, 7 removed lines. Original line number Diff line number Diff line Loading @@ -457,9 +457,9 @@ def run_and_plot_specfit(apspec, tied, line_guesses, limits, limited, # plot fit components apspec.plotter(title=f"{source} {_region}") apspec.specfit.plot_fit(show_components=True, annotate=False) plt.savefig(os.path.join( plotdir, f"{_name_prefix}_{source}_{species}{_region}.svg" ), bbox_inches='tight') # plt.savefig(os.path.join( # plotdir, f"{_name_prefix}_{source}_{species}{_region}.svg" # ), bbox_inches='tight') # return the RMS of the residuals return float(np.nanstd(apspec.specfit.residuals)) Loading Loading @@ -504,7 +504,7 @@ def overlay_region_on_map(apspec, ax_map.set_xlabel("RA [J2000]") ax_map.set_ylabel("Dec [J2000]") fig_map.tight_layout() fig_map.savefig(plotfile, bbox_inches='tight') # fig_map.savefig(plotfile, bbox_inches='tight') plt.close(fig_map) def setup_spectral_analysis(source, species): Loading Loading @@ -709,7 +709,8 @@ def process_source(source, species, guesses, prob_threshold, samples): densities=densities, moment_map=moment_map, wcs=wcs) save_results(densities, source) print(source, densities) # save_results(densities, source) def plot_density_vs_fwhm(densities, fwhms, source, species, regions): Loading @@ -720,7 +721,7 @@ def plot_density_vs_fwhm(densities, fwhms, source, species, regions): color = 'blue' if region.startswith('core') else 'cyan' plt.scatter(densities[region]['best_fit'], fwhms[region], label=region, edgecolors=color, figure=fig) plt.legend() plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_density_FWHM_{source}_{species}.svg'), bbox_inches='tight') # plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_density_FWHM_{source}_{species}.svg'), bbox_inches='tight') plt.close(fig) Loading Loading @@ -758,7 +759,7 @@ def overlay_regions_density(source, regions, densities, moment_map, wcs): plt.xlabel('RA') plt.colorbar(pc, label=r'log($n($H$_2)$ [cm$^{-3}$])', shrink=0.68) plt.tight_layout() plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_{source}_{species}.svg'), bbox_inches='tight') # plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_{source}_{species}.svg'), bbox_inches='tight') plt.close() if __name__ == '__main__': Loading etl/prs/prs_2d_kinematics.py +3 −1 Changes for etl/prs/prs_2d_kinematics.py: 3 added lines, 1 removed line. Original line number Diff line number Diff line Loading @@ -6,6 +6,8 @@ from sklearn.mixture import BayesianGaussianMixture from scipy.optimize import curve_fit from mpl_toolkits.axes_grid1 import make_axes_locatable from assets.constants import config_vis3d as config, source_properties, datadir import matplotlib matplotlib.use('svg') fontsize_small = 14 fontsize_medium = 16 Loading Loading @@ -140,7 +142,7 @@ def plot_2d_clustered_kinematics(source, species): if source == 'G332.96': ax.legend(loc='upper right', fontsize=fontsize_small, framealpha=0.9) ax.grid(alpha=0.3) ax.grid(False) plt.tight_layout() output_filename = f'projected_2d_kinematics_{source}_{species}.png' Loading etl/stg/stg_check_feather.py +22 −1 Changes for etl/stg/stg_check_feather.py: 22 added lines, 1 removed line. Original line number Diff line number Diff line import os import glob import numpy as np import matplotlib # Force headless backend to prevent the Wayland/Qt crash matplotlib.use('Agg') import matplotlib.pyplot as plt from spectral_cube import SpectralCube from regions import CircleSkyRegion Loading Loading @@ -39,7 +43,6 @@ def compare_tp_feather_spectralcube(source, species, ap_radius_arcsec=40.0, data subcube_f = cube_f.subcube_from_regions([aperture]) subcube_tp = cube_tp.subcube_from_regions([aperture]) # Enable huge ops just in case the subcube is still large subcube_f.allow_huge_operations = True subcube_tp.allow_huge_operations = True Loading @@ -59,6 +62,24 @@ def compare_tp_feather_spectralcube(source, species, ap_radius_arcsec=40.0, data spec_tp_k_interp = np.interp(freq_f, freq_tp[sort_tp], spec_tp_k[sort_tp]) diff_k = spec_f_k - spec_tp_k_interp # --- Compute Maximum Relative Residual Correctly --- peak_tp = np.nanmax(spec_tp_k_interp) # We only care about residuals in the velocity range where the line actually exists # So we mask out the pure noise baseline (e.g., regions > 5% of peak) signal_mask = spec_tp_k_interp > (0.05 * peak_tp) if np.any(signal_mask): # Calculate max absolute residual within the line profile max_abs_res = np.nanmax(np.abs(diff_k[signal_mask])) # Scale it against the peak flux to get the true percentage max_rel_res = (max_abs_res / peak_tp) * 100 print(f" [Info] Peak TP signal: {peak_tp:.3f} K") print(f" [Info] Max absolute residual: {max_abs_res:.3f} K") print(f" [Info] Max residual relative to peak: {max_rel_res:.1f}%") else: print(" [Info] Could not compute relative residual reliably.") # --- Plotting --- if source == 'G332.96': source = 'G332' Loading Loading
etl/mdl/mdl_extract_ch3oh_spectra.py +8 −7 Changes for etl/mdl/mdl_extract_ch3oh_spectra.py: 8 added lines, 7 removed lines. Original line number Diff line number Diff line Loading @@ -457,9 +457,9 @@ def run_and_plot_specfit(apspec, tied, line_guesses, limits, limited, # plot fit components apspec.plotter(title=f"{source} {_region}") apspec.specfit.plot_fit(show_components=True, annotate=False) plt.savefig(os.path.join( plotdir, f"{_name_prefix}_{source}_{species}{_region}.svg" ), bbox_inches='tight') # plt.savefig(os.path.join( # plotdir, f"{_name_prefix}_{source}_{species}{_region}.svg" # ), bbox_inches='tight') # return the RMS of the residuals return float(np.nanstd(apspec.specfit.residuals)) Loading Loading @@ -504,7 +504,7 @@ def overlay_region_on_map(apspec, ax_map.set_xlabel("RA [J2000]") ax_map.set_ylabel("Dec [J2000]") fig_map.tight_layout() fig_map.savefig(plotfile, bbox_inches='tight') # fig_map.savefig(plotfile, bbox_inches='tight') plt.close(fig_map) def setup_spectral_analysis(source, species): Loading Loading @@ -709,7 +709,8 @@ def process_source(source, species, guesses, prob_threshold, samples): densities=densities, moment_map=moment_map, wcs=wcs) save_results(densities, source) print(source, densities) # save_results(densities, source) def plot_density_vs_fwhm(densities, fwhms, source, species, regions): Loading @@ -720,7 +721,7 @@ def plot_density_vs_fwhm(densities, fwhms, source, species, regions): color = 'blue' if region.startswith('core') else 'cyan' plt.scatter(densities[region]['best_fit'], fwhms[region], label=region, edgecolors=color, figure=fig) plt.legend() plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_density_FWHM_{source}_{species}.svg'), bbox_inches='tight') # plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_density_FWHM_{source}_{species}.svg'), bbox_inches='tight') plt.close(fig) Loading Loading @@ -758,7 +759,7 @@ def overlay_regions_density(source, regions, densities, moment_map, wcs): plt.xlabel('RA') plt.colorbar(pc, label=r'log($n($H$_2)$ [cm$^{-3}$])', shrink=0.68) plt.tight_layout() plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_{source}_{species}.svg'), bbox_inches='tight') # plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_{source}_{species}.svg'), bbox_inches='tight') plt.close() if __name__ == '__main__': Loading
etl/prs/prs_2d_kinematics.py +3 −1 Changes for etl/prs/prs_2d_kinematics.py: 3 added lines, 1 removed line. Original line number Diff line number Diff line Loading @@ -6,6 +6,8 @@ from sklearn.mixture import BayesianGaussianMixture from scipy.optimize import curve_fit from mpl_toolkits.axes_grid1 import make_axes_locatable from assets.constants import config_vis3d as config, source_properties, datadir import matplotlib matplotlib.use('svg') fontsize_small = 14 fontsize_medium = 16 Loading Loading @@ -140,7 +142,7 @@ def plot_2d_clustered_kinematics(source, species): if source == 'G332.96': ax.legend(loc='upper right', fontsize=fontsize_small, framealpha=0.9) ax.grid(alpha=0.3) ax.grid(False) plt.tight_layout() output_filename = f'projected_2d_kinematics_{source}_{species}.png' Loading
etl/stg/stg_check_feather.py +22 −1 Changes for etl/stg/stg_check_feather.py: 22 added lines, 1 removed line. Original line number Diff line number Diff line import os import glob import numpy as np import matplotlib # Force headless backend to prevent the Wayland/Qt crash matplotlib.use('Agg') import matplotlib.pyplot as plt from spectral_cube import SpectralCube from regions import CircleSkyRegion Loading Loading @@ -39,7 +43,6 @@ def compare_tp_feather_spectralcube(source, species, ap_radius_arcsec=40.0, data subcube_f = cube_f.subcube_from_regions([aperture]) subcube_tp = cube_tp.subcube_from_regions([aperture]) # Enable huge ops just in case the subcube is still large subcube_f.allow_huge_operations = True subcube_tp.allow_huge_operations = True Loading @@ -59,6 +62,24 @@ def compare_tp_feather_spectralcube(source, species, ap_radius_arcsec=40.0, data spec_tp_k_interp = np.interp(freq_f, freq_tp[sort_tp], spec_tp_k[sort_tp]) diff_k = spec_f_k - spec_tp_k_interp # --- Compute Maximum Relative Residual Correctly --- peak_tp = np.nanmax(spec_tp_k_interp) # We only care about residuals in the velocity range where the line actually exists # So we mask out the pure noise baseline (e.g., regions > 5% of peak) signal_mask = spec_tp_k_interp > (0.05 * peak_tp) if np.any(signal_mask): # Calculate max absolute residual within the line profile max_abs_res = np.nanmax(np.abs(diff_k[signal_mask])) # Scale it against the peak flux to get the true percentage max_rel_res = (max_abs_res / peak_tp) * 100 print(f" [Info] Peak TP signal: {peak_tp:.3f} K") print(f" [Info] Max absolute residual: {max_abs_res:.3f} K") print(f" [Info] Max residual relative to peak: {max_rel_res:.1f}%") else: print(" [Info] Could not compute relative residual reliably.") # --- Plotting --- if source == 'G332.96': source = 'G332' Loading