Fixed case in figure titles

import numpy as np

from astropy.io import fits

from astropy.wcs import WCS

from astropy.visualization.wcsaxes import SphericalCircle

from astropy import units as u

import matplotlib.pyplot as plt

from astropy.coordinates import SkyCoord

from astropy.visualization.wcsaxes import add_scalebar

plt.rcParams['font.size'] = 16

def scale_image(data, min_val=None, max_val=None):

    """Linearly scale image data to [0, 1]."""

    if min_val is None:

        min_val = np.nanpercentile(data, 1)

    if max_val is None:

        max_val = np.nanpercentile(data, 99)

    scaled = (data - min_val) / (max_val - min_val)

    return np.clip(scaled, 0, 1)

# The distances to the sources in parsecs

distances = {

    'AGAL014.492-00.139': 3100 * u.pc,

    'AGAL028.273-00.167': 4500 * u.pc,

    'AGAL332.969-00.029': 4000 * u.pc

}

for source in ['AGAL014.492-00.139', 'AGAL028.273-00.167', 'AGAL332.969-00.029']:

    r_data = fits.getdata(f'{source}_GLM_4_reprojected.fits')

    g_data = fits.getdata(f'{source}_GLM_2_reprojected.fits')

    b_data = fits.getdata(f'{source}_GLM_1_reprojected.fits')

    wcs = WCS(fits.getheader(f'{source}_GLM_4_reprojected.fits'))

    r_scaled = scale_image(r_data)

    g_scaled = scale_image(g_data)

    b_scaled = scale_image(b_data)

    # Combine into RGB

    rgb_image = np.dstack([r_scaled, g_scaled, b_scaled])

    # Load contour image

    contour_hdu = fits.open(f'{source}_ATLASGAL_reprojected.fits')[0]

    contour_data = contour_hdu.data

    contour_wcs = WCS(contour_hdu.header)

    # Plot with WCSAxes

    fig = plt.figure(figsize=(6, 6))

    ax = plt.subplot(projection=wcs)

    # Show RGB image

    ax.imshow(rgb_image, origin='lower')

    ax.set_xlim(40, 205)

    ax.set_ylim(40, 205)

    # Add contours

    levels = np.linspace(0.1, np.nanmax(contour_data), 7)[1:]

    ax.contour(contour_data, levels=levels, colors='white',

               transform=ax.get_transform(contour_wcs), linewidths=1)

    # Add 19" circle at the center

    ny, nx = r_data.shape

    center_sky = wcs.pixel_to_world(nx / 2, ny / 2)

    circle = SphericalCircle(center_sky, 50 * u.arcsec, edgecolor='cyan',

                             facecolor='none', transform=ax.get_transform('world'), linewidth=1.5)

    ax.add_patch(circle)

    # Add 0.1 pc ruler

    dist = distances[source]

    ang_size_01pc = (0.5 * u.pc / dist * u.rad).to(u.arcsec, equivalencies=u.parallax())

    add_scalebar(ax, length=ang_size_01pc, label='0.5 pc',

                 corner='bottom left', frame=False, pad=0.5, color='white')

    # Labels

    ax.set_xlabel("RA")

    ax.set_ylabel("Dec")

    ax.set_title(source)

    plt.tight_layout()

    plt.savefig(f"{source}_rgb_atlasgal.png", dpi=300)

    plt.show()

 No newline at end of file

import glob

import os

from astropy.io import fits

from astropy.wcs import WCS

import numpy as np

from reproject import reproject_interp

for filename in glob.glob('*_reprojected.fits'):

    os.remove(filename)

lookup_sources = {

    'AGAL014.492-00.139': 'G14',

    'AGAL028.273-00.167': 'G28',

    'AGAL332.969-00.029': 'G332.96'

}

for source in lookup_sources.keys():

    for input_fits in glob.glob(f'{source}*.fits'):

        # Input and output paths

        alma_hdu = fits.open(f'../../../etl/assets/data/{lookup_sources[source]}_CH3OH_3mm_feather.fits')[0]

        alma_wcs = WCS(alma_hdu.header)

        alma_shape = alma_hdu.data.shape  # (ny, nx)

        # Determine center coordinate of ALMA image

        alma_nstokes, alma_nchans, alma_ny, alma_nx = alma_shape

        alma_center_x, alma_center_y = alma_nx / 2, alma_ny / 2

        alma_center_icrs = alma_wcs.pixel_to_world(alma_center_x, alma_center_y, 1, 1)[0]

        output_fits = input_fits.replace('.fits', '_reprojected.fits')

        # Load the original FITS and WCS

        hdu = fits.open(input_fits)[0]

        data = hdu.data

        wcs_in = WCS(hdu.header)

        # Define target WCS: Equatorial (RA/Dec), same shape as input

        shape_out = data.shape

        # Get the pixel scale from the original WCS

        cdelt = wcs_in.wcs.cd

        # Create a new WCS object in ICRS (RA/Dec)

        target_wcs = WCS(naxis=2)

        target_wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN']

        ny, nx = data.shape  # note: shape is (ny, nx)

        x_center, y_center = nx / 2, ny / 2

        # Get sky coord of the center in Galactic

        center_coord_gal = alma_wcs.pixel_to_world(1,1,x_center, y_center)[0]

        # Convert to ICRS (RA/Dec)

        center_coord_icrs = center_coord_gal.transform_to('icrs')

        # Set the reference value in RA/Dec

        target_wcs.wcs.crval = [alma_center_icrs.ra.deg, alma_center_icrs.dec.deg]

        target_wcs.wcs.crpix = [x_center, y_center]

        target_wcs.wcs.cd = cdelt

        target_wcs.wcs.cunit = ['deg', 'deg']

        # Reproject the image

        reprojected_data, footprint = reproject_interp((data, wcs_in), target_wcs, shape_out=shape_out)

        # Create new HDU and write to file

        new_hdu = fits.PrimaryHDU(data=reprojected_data, header=target_wcs.to_header())

        new_hdu.writeto(output_fits, overwrite=True)

        print(f"Reprojected FITS saved to: {output_fits}")

\documentclass[class=minimal,border=0pt]{standalone}

\usepackage{graphicx}

\usepackage{tikz}

\usetikzlibrary{calc}

\usetikzlibrary{arrows.meta}

\usepackage{svg}

% \usepackage[margin=1in]{geometry}

\definecolor{cy}{HTML}{00ffff}

\begin{document}

  \begin{tikzpicture}

    \draw[draw=black,

    line width=1.5mm,              % thickness of the shaft

    -{Latex[length=8mm,width=5mm]} % arrow tip style and size

  ] (-1,20.5) -- node[midway, below, scale=1.5, rotate=90] {$L/M$}  (-1,2);

    \begin{scope}[shift={(0,14)}]

        % Main image

        \node[anchor=south west, inner sep=0] (main) at (0,0)

          {\includegraphics[width=0.45\textwidth]{AGAL332.969-00.029_rgb_atlasgal.png}}; % Replace with your large-scale figure

        \node[anchor=south west] (inset) at (17, 0.25) % adjust as needed

          {\includesvg[width=0.615\textwidth]{G332.96_cont.svg}}; % Replace with your inset image

        \begin{scope}

            % Define the clipping area.

            % We start from the bottom-left of the intended node (7, 0.15)

            % The width of the SVG is 0.75\textwidth, which is roughly 11.4 units in your coordinate system.

            % We clip a rectangle but "cut off" the bottom right.

%             \clip (7, 0.15) -- (18.4, 0.15) -- (18.4, 1.5) -- (16.5, 0.15) -- cycle;

            % Note: The path above is an example. For a simple rectangular clip with a corner missing:

%             \clip (7, 0.15) rectangle (18.4, 8); % Standard box

            \newsavebox{\myimage}

            \sbox{\myimage}{\includesvg[width=0.75\textwidth]{aperture_spectra_extraction_G332.96_CH3OH.svg}}

            \begin{scope}[shift={(7.11, 0.27)}]

                % Define the clip path relative to (0,0) which is now (7, 0.15)

                \clip (0,0) -- (0, \ht\myimage) -- (\wd\myimage, \ht\myimage)

                      -- (\wd\myimage, 1.5)              % Vertical offset up from bottom-right

                      -- ({\wd\myimage-1.5cm}, 0)        % Horizontal offset left from bottom-right

                      -- cycle;

                % Place the node with ZERO padding to prevent shifting

                \node[anchor=south west, inner sep=0pt, outer sep=0pt] (dens) at (0,0) {\usebox{\myimage}};

            \end{scope}

%             \clip (dens.south west) -- (dens.north west) -- (dens.north east) -- (18.4, 8) -- (7, 8) -- cycle; % Clipped corner

        \end{scope}

%         \node[anchor=south west] (dens) at (7, 0.15) % adjust as needed

%           {\includesvg[width=0.75\textwidth]{aperture_spectra_extraction_G332.96_CH3OH.svg}}; % Replace with your inset image

        % Coordinate system for overlay

        \begin{scope}[x={(main.south east)}, y={(main.north west)}]

          % Draw a line to the inset

          \draw[cy, thick] (0.61,0.396) -- (1.305,0.19);

          \draw[cy, thick] (0.61,0.743) -- (1.305,0.94);

        \end{scope}

        % Inset image (positioned absolutely)

    %     \node[anchor=south west, draw=blue!80, thick, inner sep=0pt] at (7, 5) % adjust as needed

    %       {\includesvg[width=0.5\textwidth]{G14_CH3OH_cont.svg}}; % Replace with your inset image

    %     \node[anchor=south west, draw=blue!80, thick, inner sep=0pt] at (7, -2) % adjust as needed

    %       {\includesvg[width=0.5\textwidth]{G14_CH3OH_mom_zero.svg}}; % Replace with your inset image

    %     \node[anchor=south west] at (9,-8) % adjust as needed

    %       {\includesvg[width=0.5\textwidth]{G14_N2H+_mom_zero.svg}}; % Replace with your inset image

    %     \node[anchor=south west] at (4.4, -5.3) % adjust as needed

    %       {\includesvg[width=0.5\textwidth]{G14_CH3OH_mom_zero.svg}}; % Replace with your inset image

    %     \node[anchor=south west] at (-3.8, -5.3) % adjust as needed

    %       {\includesvg[width=0.5\textwidth]{G14_CH3OH_cont.svg}}; % Replace with your inset image

    \end{scope}

    \begin{scope}[shift={(0,7)}]

        % Main image

        \node[anchor=south west, inner sep=0] (main2) at (0.1,0)

          {\includegraphics[width=0.45\textwidth]{AGAL028.273-00.167_rgb_atlasgal.png}}; % Replace with your large-scale figure

        \node[anchor=south west] (inset) at (17.3, 0.25) % adjust as needed

          {\includesvg[width=0.593\textwidth]{G28_cont.svg}}; % Replace with your inset image

        \node[anchor=south west] (dens) at (7.23, 0.15) % adjust as needed

          {\includesvg[width=0.735\textwidth]{aperture_spectra_extraction_G28_CH3OH.svg}}; % Replace with your inset image

            % Coordinate system for overlay

        \begin{scope}[

                        x={($(main2.south east)-(main2.south west)$)},

                        y={($(main2.north west)-(main2.south west)$)},

                        shift={(main2.south west)}

                     ]

          % Draw a line to the inset

          \draw[cy, thick] (0.61,0.396) -- (1.292,0.189);

          \draw[cy, thick] (0.61,0.743) -- (1.292,0.941);

        \end{scope}

    \end{scope}

    \begin{scope}[shift={(0,0)}]

        % Main image

        \node[anchor=south west, inner sep=0] (main3) at (0,0)

          {\includegraphics[width=0.45\textwidth]{AGAL014.492-00.139_rgb_atlasgal.png}}; % Replace with your large-scale figure

        \node[anchor=south west] (inset) at (17, 0.25) % adjust as needed

          {\includesvg[width=0.63\textwidth]{G14_cont.svg}}; % Replace with your inset image

        \node[anchor=south west] (dens) at (7, 0.15) % adjust as needed

          {\includesvg[width=0.75\textwidth]{aperture_spectra_extraction_G14_CH3OH.svg}}; % Replace with your inset image

        % Coordinate system for overlay

        \begin{scope}[x={(main3.south east)}, y={(main3.north west)}]

          % Draw a line to the inset

          \draw[cy, thick] (0.61,0.396) -- (1.305,0.19);

          \draw[cy, thick] (0.61,0.743) -- (1.305,0.94);

        \end{scope}

    \end{scope}

  \end{tikzpicture}

\end{document}

import matplotlib.pyplot as plt

import matplotlib

import numpy as np

import yaml

import os

matplotlib.use('svg')

matplotlib.rc('font', size=16)

def compute_density_contrast(data):

    core_values, core_weights = [], []

    envelope_values, envelope_weights = [], []

    for key, item in data.items():

        best_fit_list = item.get("best_fit", [])

        reliability_list = item.get("reliability_score", [])

        uncertainty_list = item.get("uncertainty", [])

        if not reliability_list:

            continue

        highest_reliability_index = np.argmax(reliability_list)

        best_fit = best_fit_list[highest_reliability_index]

        uncertainty_range = uncertainty_list[highest_reliability_index]

        sigma_log = (uncertainty_range[1] - uncertainty_range[0]) / 2

        if sigma_log == 0:

            continue

        weight = 1 / (sigma_log ** 2)

        val_linear = 10 ** best_fit

        if key.startswith("core"):

            core_values.append(val_linear)

            core_weights.append(weight)

        elif key.startswith("envelope"):

            envelope_values.append(val_linear)

            envelope_weights.append(weight)

    np_core_values = np.array(core_values)

    np_core_weights = np.array(core_weights)

    np_envelope_values = np.array(envelope_values)

    np_envelope_weights = np.array(envelope_weights)

    if np_core_weights.sum() == 0 or np_envelope_weights.sum() == 0:

        return None, None

    mean_core = np.average(np_core_values, weights=np_core_weights)

    mean_envelope = np.average(np_envelope_values, weights=np_envelope_weights)

    # Use the weighted sample variance to capture intrinsic physical scatter

    if len(np_core_values) > 1:

        var_core = np.average((np_core_values - mean_core) ** 2, weights=np_core_weights)

        err_core = np.sqrt(var_core / len(np_core_values))

    else:

        err_core = mean_core * 0.1

    if len(np_envelope_values) > 1:

        var_env = np.average((np_envelope_values - mean_envelope) ** 2, weights=np_envelope_weights)

        err_envelope = np.sqrt(var_env / len(np_envelope_values))

    else:

        err_envelope = mean_envelope * 0.1

    density_contrast = mean_core / mean_envelope

    contrast_err = density_contrast * np.sqrt((err_core / mean_core) ** 2 + (err_envelope / mean_envelope) ** 2)

    return density_contrast, contrast_err

def normalize_and_get_errors(val_array, err_array, err_type='linear', ref_idx=0):

    """

    Normalizes variables to a reference index and computes asymmetric errors.

    Uses log-space propagation for highly uncertain variables to preserve strictly positive bounds.

    """

    ref_val = val_array[ref_idx]

    rel_val = val_array / ref_val

    if err_type == 'linear':

        # err_array contains absolute linear errors

        rel_err = err_array / ref_val

        rel_err[ref_idx] = 0.0  # Reference point has 0 error in normalization

        return rel_val, np.vstack((rel_err, rel_err))

    elif err_type == 'log':

        # err_array contains log10 errors (dex)

        factors = 10 ** err_array

        lower_err = rel_val - (rel_val / factors)

        upper_err = (rel_val * factors) - rel_val

        lower_err[ref_idx] = 0.0

        upper_err[ref_idx] = 0.0

        return rel_val, np.vstack((lower_err, upper_err))

# --- 1. Base Data ---

nc_a_raw = np.array([5, 13, 23])

nc_a = np.array([4, 8, 28])

cfe = np.array([5.9, 9.8, 23.4])

m_cl_m_j = np.array([444, 1065, 1136])

ic_vel_disp = np.array([0.8, 1.4, 1.9])

ch3oh_abund = np.array([2.3e-9, 3.8e-9, 4.3e-9])

fwhm = np.array([3.66, 4.33, 4.13])

T_dust = np.array([12.2, 11.8, 16.6])

sound_speed = np.array([0.22, 0.22, 0.258])

average_clump_density = np.array([48000, 53000, 108000])

mach_number = fwhm / 2.355 / sound_speed

nc_a_fwhm = nc_a / mach_number

# --- 2. Error Definitions (Linear & Log-Space) ---

# Linear Errors

T_err = np.array([2.0, 2.0, 2.0])

fwhm_err_frac = np.array([0.05, 0.05, 0.05])

ch3oh_abund_err_frac = np.array([0.30, 0.30, 0.30])

nc_err_frac = 1.0 / np.sqrt(nc_a_raw)

ic_vel_disp_err_frac = 1.0 / np.sqrt(2 * (nc_a_raw - 1))

cfe_err_frac = nc_err_frac

sound_speed_err_frac = 0.5 * (T_err / T_dust)

mach_err_frac = np.sqrt(fwhm_err_frac ** 2 + sound_speed_err_frac ** 2)

fragmentation_err_frac = np.sqrt(nc_err_frac ** 2 + mach_err_frac ** 2)

# Log Errors (dex)

clump_dens_log_err = np.array([

    (4.95 - 4.26) / 2.0,

    (4.88 - 4.57) / 2.0,

    (5.18 - 4.90) / 2.0

])

mass_log_err = 0.5

T_log_err = (1.0 / np.log(10)) * (T_err / T_dust)

# Log-space propagation for M_cl / M_J

m_cl_m_j_log_err = np.sqrt(mass_log_err ** 2 + (0.5 * clump_dens_log_err) ** 2 + (1.5 * T_log_err) ** 2)

# --- 3. Compute Density Contrast ---

density_contrast, density_contrast_err = [], []

for source in ['G332.96', 'G28', 'G14']:

    file_path = os.path.join('mdl', 'inferred_densities', f'densities_{source}.yml')

    if os.path.exists(file_path):

        with open(file_path) as infile:

            data = yaml.load(infile, Loader=yaml.FullLoader)

        val, err = compute_density_contrast(data=data)

        density_contrast.append(val)

        density_contrast_err.append(err)

    else:

        density_contrast.append(1.0)

        density_contrast_err.append(0.1)

density_contrast = np.array(density_contrast)

density_contrast_err = np.array(density_contrast_err)

# --- 4. Normalize and Retrieve Asymmetric Error Arrays ---

density_contrast_rel, d_contrast_rel_err = normalize_and_get_errors(density_contrast, density_contrast_err,

                                                                    err_type='linear')

nc_a_rel, nc_a_rel_err = normalize_and_get_errors(nc_a, nc_a * nc_err_frac, err_type='linear')

cfe_rel, cfe_rel_err = normalize_and_get_errors(cfe, cfe * cfe_err_frac, err_type='linear')

ch3oh_abund_rel, ch3oh_abund_rel_err = normalize_and_get_errors(ch3oh_abund, ch3oh_abund * ch3oh_abund_err_frac,

                                                                err_type='linear')

fragmentation_rel, fragmentation_rel_err = normalize_and_get_errors(nc_a_fwhm, nc_a_fwhm * fragmentation_err_frac,

                                                                    err_type='linear')

ic_vel_disp_rel, ic_vel_disp_rel_err = normalize_and_get_errors(ic_vel_disp, ic_vel_disp * ic_vel_disp_err_frac,

                                                                err_type='linear')

# Variables natively requiring Log-Space propagation

average_clump_density_rel, avg_dens_rel_err = normalize_and_get_errors(average_clump_density, clump_dens_log_err,

                                                                       err_type='log')

m_cl_m_j_rel, m_cl_m_j_rel_err = normalize_and_get_errors(m_cl_m_j, m_cl_m_j_log_err, err_type='log')

# --- 5. Single-Panel Linear Plotting with Jitter ---

plt.style.use('seaborn-v0_8-whitegrid')

fig, ax = plt.subplots(figsize=(15, 7))

x_vals = np.array([1, 2, 3])

x_labels = ['G332', 'G28', 'G14']

arrays_to_plot = {

    r'Density contrast ($\frac{\overline{n_{core}\mathrm{(H_2)}}}{\overline{n_{env}\mathrm{(H_2)}}}$)': (

        density_contrast_rel, d_contrast_rel_err),

    'Cores per unit area ($N_c/A$)': (nc_a_rel, nc_a_rel_err),

    'CFE': (cfe_rel, cfe_rel_err),

    'Jeans mass count ($M_{cl}/M_J$)': (m_cl_m_j_rel, m_cl_m_j_rel_err),

    'Methanol abundance ($X_{CH_3OH}$)': (ch3oh_abund_rel, ch3oh_abund_rel_err),

    '$N_c/A$ per Mach number': (fragmentation_rel, fragmentation_rel_err),

    'Clump $\overline{n_{cl}\mathrm{(H_2)}}$': (average_clump_density_rel, avg_dens_rel_err),

    'Core-to-core $\sigma_{V_{LSR}}$': (ic_vel_disp_rel, ic_vel_disp_rel_err)

}

colors = plt.cm.viridis(np.linspace(0.05, 0.95, len(arrays_to_plot)))

markers = ['o', 's', '^', 'D', '*', 'h', 'P', 'v']

num_vars = len(arrays_to_plot)

jitter_offsets = np.linspace(-0.2, 0.2, num_vars)

# Plotting Loop

for idx, (qty_label, (y_vals, y_errs)) in enumerate(arrays_to_plot.items()):

    x_jittered = x_vals + jitter_offsets[idx]

    # Faint lines connecting points

    ax.plot(x_jittered, y_vals, color=colors[idx], linestyle='--', alpha=0.3, linewidth=1.5)

    # Bold markers and asymmetric error bars natively handled by matplotlib

    ax.errorbar(x_jittered, y_vals, yerr=y_errs, fmt=markers[idx], color=colors[idx],

                label=qty_label, markersize=10, fillstyle='none',

                capsize=4, elinewidth=1.8, markeredgewidth=1.8)

# Formatting Linear Scale

ax.axhline(1.0, color='grey', linestyle='-', linewidth=2)

ax.text(3.3, 0.4, 'Turbulence\nBaseline', color='grey', ha='right', va='center', fontsize=14)

# Draw the green arrow

ax.annotate('', xy=(2.9, 8.5), xytext=(1.1, 3.0),

            arrowprops=dict(facecolor='green', edgecolor='none', shrink=0.05, width=2.5, headwidth=12, headlength=12))

# Rotated arrow text

ax.text(1.5, 4.0, 'Gravity-dominated', color='green', fontweight='bold',

        fontsize=15, rotation=17, ha='center', va='bottom')

ax.set_title('Evolutionary discriminator: gravity vs. turbulence', fontsize=20, fontweight='bold', pad=20)

ax.set_xlabel('Source (proposed evolutionary sequence)', fontsize=18, labelpad=15)

ax.set_ylabel('Relative change (normalized to G332)', fontsize=18, labelpad=15)

ax.set_xticks(x_vals)

ax.set_xticklabels(x_labels, fontsize=16)

ax.set_yscale('linear')

ax.set_ylim(-0.5, 10.0)

# ax.grid(True, which="major", ls="--", lw=0.5, c='0.8')

ax.grid()

# Legend placed in the empty top-left quadrant

ax.legend(loc='upper left', frameon=True, fontsize=14, ncol=2)

plt.tight_layout()

output_dir = os.path.join('prs', 'relative_variation_figure_input')

os.makedirs(output_dir, exist_ok=True)

plt.savefig(os.path.join(output_dir, 'gravity_vs_turbulence_discriminator_plot.png'), dpi=300, bbox_inches='tight')

plt.show()

 No newline at end of file