<commit_before>
import astropy.units as u
from spectral_cube import SpectralCube, OneDSpectrum
import numpy as np
from astropy.io import fits
from os.path import join as osjoin
import os
from astropy.coordinates import Angle

from cube_analysis.spectral_stacking import radial_stacking

from paths import (fourteenB_wGBT_HI_file_dict,
                   iram_co21_14B088_data_path,
                   fourteenB_HI_data_wGBT_path)
from constants import hi_freq
from galaxy_params import gal_feath as gal

'''
Create profiles of HI and CO after subtracting velocity surfaces
'''

co_stackpath = lambda x: osjoin(iram_co21_14B088_data_path("", no_check=True), "stacked_spectra", x)
if not os.path.exists(co_stackpath("")):
    os.mkdir(co_stackpath(""))

dr = 500 * u.pc
max_radius = (7.0 * u.kpc).to(u.pc)
wstring = "{0}{1}".format(int(dr.value), dr.unit)
maxrad_string = "{0}{1}".format(int(max_radius.value), max_radius.unit)

pa_bounds_n = Angle([0.5 * np.pi * u.rad, -0.5 * np.pi * u.rad])
pa_bounds_n = pa_bounds_n.wrap_at(0.5 * np.pi * u.rad)

pa_bounds_s = Angle([-0.5 * np.pi * u.rad, 0.5 * np.pi * u.rad])

# CO stacking first
co_cube = SpectralCube.read(iram_co21_14B088_data_path("m33.co21_iram.14B-088_HI_feather.rotation_corrected.fits"))

bin_centers, total_spectrum_co_radial, num_pixels = \
    radial_stacking(gal, co_cube, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=None,
                    verbose=True,
                    how='cube')

bin_centers, total_spectrum_co_radial_n, num_pixels_n = \
    radial_stacking(gal, co_cube, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=pa_bounds_n,
                    verbose=True,
                    how='cube')

bin_centers, total_spectrum_co_radial_s, num_pixels_s = \
    radial_stacking(gal, co_cube, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=pa_bounds_s,
                    verbose=True,
                    how='cube')

spec_shape = co_cube.shape[0]

rot_stack = SpectralCube(data=total_spectrum_co_radial.T.reshape((spec_shape, bin_centers.size, 1)),
                         wcs=co_cube.wcs)
rot_stack.write(co_stackpath("rotation_stacked_radial_{}.fits".format(wstring)),
                overwrite=True)
rot_stack_n = SpectralCube(data=total_spectrum_co_radial_n.T.reshape((spec_shape, bin_centers.size, 1)),
                           wcs=co_cube.wcs)
rot_stack_n.write(co_stackpath("rotation_stacked_radial_north_{}.fits".format(wstring)),
                  overwrite=True)
rot_stack_s = SpectralCube(data=total_spectrum_co_radial_s.T.reshape((spec_shape, bin_centers.size, 1)),
                           wcs=co_cube.wcs)
rot_stack_s.write(co_stackpath("rotation_stacked_radial_south_{}.fits".format(wstring)),
                  overwrite=True)

# Separately save the number of pixels in each bin
np.save(co_stackpath("radial_stacking_pixelsinbin_{}.npy").format(wstring), num_pixels)
np.save(co_stackpath("radial_stacking_pixelsinbin_north_{}.npy").format(wstring), num_pixels_n)
np.save(co_stackpath("radial_stacking_pixelsinbin_south_{}.npy").format(wstring), num_pixels_s)


# Save the total profiles over the inner 7 kpc

total_spectrum_co = total_spectrum_co_radial.sum(0)

oned_wcs = co_cube[:, 0, 0].wcs
OneDSpectrum(total_spectrum_co.value, unit=total_spectrum_co.unit,
             wcs=oned_wcs).write(co_stackpath("rotation_stacked_{}.fits".format(maxrad_string)),
                                 overwrite=True)

del co_cube

co_cube_cent = \
    SpectralCube.read(iram_co21_14B088_data_path("m33.co21_iram.14B-088_HI_feather.centroid_corrected.fits"))

total_spectrum_co_radial_cent = \
    radial_stacking(gal, co_cube_cent, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=None,
                    verbose=True,
                    how='cube')[1]

total_spectrum_co_radial_n_cent = \
    radial_stacking(gal, co_cube_cent, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=pa_bounds_n,
                    verbose=True,
                    how='cube')

total_spectrum_co_radial_s_cent = \
    radial_stacking(gal, co_cube_cent, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=pa_bounds_s,
                    verbose=True,
                    how='cube')

spec_shape = co_cube_cent.shape[0]

cent_stack = SpectralCube(data=total_spectrum_co_radial_cent.T.reshape((spec_shape, bin_centers.size, 1)),
                          wcs=co_cube_cent.wcs)
cent_stack.write(co_stackpath("centroid_stacked_radial_{}.fits".format(wstring)),
                 overwrite=True)

cent_stack_n = SpectralCube(data=total_spectrum_co_radial_cent_n.T.reshape((spec_shape, bin_centers.size, 1)),
                            wcs=co_cube_cent.wcs)
cent_stack_n.write(co_stackpath("centroid_stacked_radial_north_{}.fits".format(wstring)),
                   overwrite=True)
cent_stack_s = SpectralCube(data=total_spectrum_co_radial_cent_s.T.reshape((spec_shape, bin_centers.size, 1)),
                            wcs=co_cube_cent.wcs)
cent_stack_s.write(co_stackpath("centroid_stacked_radial_south_{}.fits".format(wstring)),
                   overwrite=True)

total_spectrum_co_cent = total_spectrum_co_radial_cent.sum(0)
oned_wcs = co_cube_cent[:, 0, 0].wcs
OneDSpectrum(total_spectrum_co_cent.value, unit=total_spectrum_co_cent.unit,
             wcs=oned_wcs).write(co_stackpath("centroid_stacked_{}.fits".format(maxrad_string)),
                                 overwrite=True)

del co_cube_cent

co_cube_peakvel = \
    SpectralCube.read(iram_co21_14B088_data_path("m33.co21_iram.14B-088_HI_feather.peakvels_corrected.fits"))

total_spectrum_co_radial_peakvel = \
    radial_stacking(gal, co_cube_peakvel, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=None,
                    verbose=True,
                    how='cube')[1]

total_spectrum_co_radial_n_peakvel = \
    radial_stacking(gal, co_cube_peakvel, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=pa_bounds_n,
                    verbose=True,
                    how='cube')

total_spectrum_co_radial_s_peakvel = \
    radial_stacking(gal, co_cube_peakvel, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=pa_bounds_s,
                    verbose=True,
                    how='cube')
spec_shape = co_cube_peakvel.shape[0]<commit_msg>Fix the N/S masking for CO stacking<commit_after>
import astropy.units as u
from spectral_cube import SpectralCube, OneDSpectrum
import numpy as np
from astropy.io import fits
from os.path import join as osjoin
import os
from astropy.coordinates import Angle

from cube_analysis.spectral_stacking import radial_stacking

from paths import (fourteenB_wGBT_HI_file_dict,
                   iram_co21_14B088_data_path,
                   fourteenB_HI_data_wGBT_path)
from constants import hi_freq
from galaxy_params import gal_feath as gal

'''
Create profiles of HI and CO after subtracting velocity surfaces
'''

co_stackpath = lambda x: osjoin(iram_co21_14B088_data_path("", no_check=True), "stacked_spectra", x)
if not os.path.exists(co_stackpath("")):
    os.mkdir(co_stackpath(""))

dr = 500 * u.pc
max_radius = (7.0 * u.kpc).to(u.pc)
wstring = "{0}{1}".format(int(dr.value), dr.unit)
maxrad_string = "{0}{1}".format(int(max_radius.value), max_radius.unit)

pa_bounds_n = Angle([0.5 * np.pi * u.rad, -0.5 * np.pi * u.rad])

pa_bounds_s = Angle([-0.5 * np.pi * u.rad, 0.5 * np.pi * u.rad])

# CO stacking first
co_cube = SpectralCube.read(iram_co21_14B088_data_path("m33.co21_iram.14B-088_HI_feather.rotation_corrected.fits"))

bin_centers, total_spectrum_co_radial, num_pixels = \
    radial_stacking(gal, co_cube, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=None,
                    verbose=True,
                    how='cube')

bin_centers, total_spectrum_co_radial_n, num_pixels_n = \
    radial_stacking(gal, co_cube, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=pa_bounds_n,
                    verbose=True,
                    how='cube')

bin_centers, total_spectrum_co_radial_s, num_pixels_s = \
    radial_stacking(gal, co_cube, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=pa_bounds_s,
                    verbose=True,
                    how='cube')

spec_shape = co_cube.shape[0]

rot_stack = SpectralCube(data=total_spectrum_co_radial.T.reshape((spec_shape, bin_centers.size, 1)),
                         wcs=co_cube.wcs)
rot_stack.write(co_stackpath("rotation_stacked_radial_{}.fits".format(wstring)),
                overwrite=True)
rot_stack_n = SpectralCube(data=total_spectrum_co_radial_n.T.reshape((spec_shape, bin_centers.size, 1)),
                           wcs=co_cube.wcs)
rot_stack_n.write(co_stackpath("rotation_stacked_radial_north_{}.fits".format(wstring)),
                  overwrite=True)
rot_stack_s = SpectralCube(data=total_spectrum_co_radial_s.T.reshape((spec_shape, bin_centers.size, 1)),
                           wcs=co_cube.wcs)
rot_stack_s.write(co_stackpath("rotation_stacked_radial_south_{}.fits".format(wstring)),
                  overwrite=True)

# Separately save the number of pixels in each bin
np.save(co_stackpath("radial_stacking_pixelsinbin_{}.npy").format(wstring), num_pixels)
np.save(co_stackpath("radial_stacking_pixelsinbin_north_{}.npy").format(wstring), num_pixels_n)
np.save(co_stackpath("radial_stacking_pixelsinbin_south_{}.npy").format(wstring), num_pixels_s)


# Save the total profiles over the inner 7 kpc

total_spectrum_co = total_spectrum_co_radial.sum(0)

oned_wcs = co_cube[:, 0, 0].wcs
OneDSpectrum(total_spectrum_co.value, unit=total_spectrum_co.unit,
             wcs=oned_wcs).write(co_stackpath("rotation_stacked_{}.fits".format(maxrad_string)),
                                 overwrite=True)

del co_cube

co_cube_cent = \
    SpectralCube.read(iram_co21_14B088_data_path("m33.co21_iram.14B-088_HI_feather.centroid_corrected.fits"))

total_spectrum_co_radial_cent = \
    radial_stacking(gal, co_cube_cent, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=None,
                    verbose=True,
                    how='cube')[1]

total_spectrum_co_radial_cent_n = \
    radial_stacking(gal, co_cube_cent, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=pa_bounds_n,
                    verbose=True,
                    how='cube')[1]

total_spectrum_co_radial_cent_s = \
    radial_stacking(gal, co_cube_cent, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=pa_bounds_s,
                    verbose=True,
                    how='cube')[1]

spec_shape = co_cube_cent.shape[0]


cent_stack = SpectralCube(data=total_spectrum_co_radial_cent.T.reshape((spec_shape, bin_centers.size, 1)),
                          wcs=co_cube_cent.wcs)
cent_stack.write(co_stackpath("centroid_stacked_radial_{}.fits".format(wstring)),
                 overwrite=True)

cent_stack_n = SpectralCube(data=total_spectrum_co_radial_cent_n.T.reshape((spec_shape, bin_centers.size, 1)),
                            wcs=co_cube_cent.wcs)
cent_stack_n.write(co_stackpath("centroid_stacked_radial_north_{}.fits".format(wstring)),
                   overwrite=True)
cent_stack_s = SpectralCube(data=total_spectrum_co_radial_cent_s.T.reshape((spec_shape, bin_centers.size, 1)),
                            wcs=co_cube_cent.wcs)
cent_stack_s.write(co_stackpath("centroid_stacked_radial_south_{}.fits".format(wstring)),
                   overwrite=True)

total_spectrum_co_cent = total_spectrum_co_radial_cent.sum(0)
oned_wcs = co_cube_cent[:, 0, 0].wcs
OneDSpectrum(total_spectrum_co_cent.value, unit=total_spectrum_co_cent.unit,
             wcs=oned_wcs).write(co_stackpath("centroid_stacked_{}.fits".format(maxrad_string)),
                                 overwrite=True)

del co_cube_cent

co_cube_peakvel = \
    SpectralCube.read(iram_co21_14B088_data_path("m33.co21_iram.14B-088_HI_feather.peakvels_corrected.fits"))

total_spectrum_co_radial_peakvel = \
    radial_stacking(gal, co_cube_peakvel, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=None,
                    verbose=True,
                    how='cube')[1]

total_spectrum_co_radial_peakvel_n = \
    radial_stacking(gal, co_cube_peakvel, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=pa_bounds_n,
                    verbose=True,
                    how='cube')[1]

total_spectrum_co_radial_peakvel_s = \
    radial_stacking(gal, co_cube_peakvel, dr=dr,
                    max_radius=max_radius,
                    pa_bounds=pa_bounds_s,
                    verbose=True,
                    how='cube')[1]

spec_shape = co_cube_peakvel.shape[0]
{'fraction_non_alphanumeric': 0.10068933467585779, 'fraction_numerical': 0.015490666873208893, 'mean_word_length': 2.853178155774396, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7587487', 'n_tokens_mistral': 4829, 'n_tokens_neox': 4517, 'n_words': 569}