B3_spectra.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b7cefb02-66f0-4f11-8dad-636abcb90e03",
   "metadata": {},
   "source": [
    "### For Juan and Walker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1a3ff657-2dff-4582-9d66-fcc4d7317aa1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Based on plot_all_spectra_NEW.ipynb\n",
    "# Also based on modeling_core.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8755420c-78af-411b-84d4-caa882967ddf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For parallelization, which helps to make convolution faster\n",
    "import dask\n",
    "dask.config.set(scheduler = 'threads', num_workers = 8)\n",
    "from dask.diagnostics import ProgressBar\n",
    "ProgressBar().register()\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7fb6cbcc-4161-4718-993f-955230b412bd",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "make: *** No targets specified and no makefile found.  Stop.\n"
     ]
    }
   ],
   "source": [
    "from spectral_cube import OneDSpectrum\n",
    "import pyspeckit\n",
    "from astropy.io import fits\n",
    "from astropy import units as u\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b07cbd33-a8b9-4b0c-976d-4f33acaab7bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = '/blue/adamginsburg/abulatek/brick/symlinks/imaging_results/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "24f93510-26ed-412c-8c39-f8b783e96364",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get all spws for B3\n",
    "freq_spws = [\n",
    "'87_spw25',\n",
    "'89_spw27',\n",
    "'91_spw25', \n",
    "'93_spw27',\n",
    "'95_spw25', \n",
    "'97_spw27',\n",
    "'98_spw29',\n",
    "'99_spw31',\n",
    "'102_spw23',\n",
    "'102_spw29',\n",
    "'103_spw31',\n",
    "'106_spw29',\n",
    "'107_spw31',\n",
    "'110_spw29',\n",
    "'111_spw31',\n",
    "'112_spw27',\n",
    "'114_spw29']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9f397855-3336-4e47-b5a3-721dcdfcbc51",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[########################################] | 100% Completed | 93.77 s\n",
      "[########################################] | 100% Completed | 93.65 s\n",
      "[########################################] | 100% Completed | 93.40 s\n",
      "[########################################] | 100% Completed | 92.94 s\n",
      "[########################################] | 100% Completed | 93.32 s\n",
      "[########################################] | 100% Completed | 92.86 s\n",
      "[########################################] | 100% Completed | 93.83 s\n",
      "[########################################] | 100% Completed | 93.71 s\n",
      "[########################################] | 100% Completed | 94.13 s\n",
      "[########################################] | 100% Completed | 94.28 s\n",
      "[########################################] | 100% Completed | 93.74 s\n",
      "[########################################] | 100% Completed | 93.90 s\n",
      "[########################################] | 100% Completed | 94.13 s\n",
      "[########################################] | 100% Completed | 93.78 s\n",
      "[########################################] | 100% Completed | 93.80 s\n",
      "[########################################] | 100% Completed | 93.87 s\n",
      "[########################################] | 100% Completed | 93.98 s\n",
      "[########################################] | 100% Completed | 94.07 s\n",
      "[########################################] | 100% Completed | 95.46 s\n",
      "[########################################] | 100% Completed | 95.34 s\n",
      "[########################################] | 100% Completed | 95.26 s\n",
      "[########################################] | 100% Completed | 94.55 s\n",
      "[########################################] | 100% Completed | 110.87 s\n",
      "[########################################] | 100% Completed | 94.91 s\n",
      "[########################################] | 100% Completed | 98.12 s\n",
      "[########################################] | 100% Completed | 96.37 s\n",
      "[########################################] | 100% Completed | 96.46 s\n",
      "[########################################] | 100% Completed | 101.99 s\n",
      "[########################################] | 100% Completed | 103.71 s\n",
      "[########################################] | 100% Completed | 101.30 s\n",
      "[########################################] | 100% Completed | 99.99 s\n",
      "[########################################] | 100% Completed | 99.67 s\n",
      "[########################################] | 100% Completed | 96.54 s\n",
      "[########################################] | 100% Completed | 95.01 s\n"
     ]
    }
   ],
   "source": [
    "from spectral_cube import SpectralCube\n",
    "from astropy import coordinates\n",
    "from astropy import wcs\n",
    "\n",
    "for spw_name in freq_spws:\n",
    "    # Get cube based on spw name and smooth it to a common beam\n",
    "    ind = freq_spws.index(spw_name) # Silly way to do this, but easiest at the moment\n",
    "    cube = SpectralCube.read(results+'source_ab_'+freq_spws[ind]+'_clean_2sigma_n50000_masked_3sigma_pbmask0p18.image', use_dask = True)\n",
    "    cube_common_beam = cube.beams.common_beam(max_iter = 20, max_epsilon = 0.01)\n",
    "    smoothed_cube = cube.convolve_to(cube_common_beam) # Convert from VaryingResolution to regular\n",
    "    # Extract spectrum for a coordinate in central core region\n",
    "    crd = coordinates.SkyCoord(\"17:46:10.63 -28:42:17.8\", frame='fk5', unit=(u.h, u.deg)) # This is the central core, about\n",
    "    x, y = map(int, smoothed_cube.wcs.celestial.world_to_pixel(crd))\n",
    "    spectrum = smoothed_cube[:, y, x].to(u.K)\n",
    "    spectrum.write(f'B3_spectra_for_juan/{spw_name}.fits', format='fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "e0f4992e-fae2-47e1-a6eb-093803961d74",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "from astropy import wcs\n",
    "from astropy import units as u\n",
    "from radio_beam import Beam\n",
    "from spectral_cube import OneDSpectrum\n",
    "\n",
    "# Verify that you can read in and plot the spectrum in this form\n",
    "# Paraphrased from spectral_cube source code\n",
    "\n",
    "spw_name_example = '106_spw29'\n",
    "\n",
    "with fits.open(f'{spw_name_example}.fits') as hdul:\n",
    "    # spectrum_in = OneDSpectrum.from_hdu(hdul) # Doesn't work, gives \"beam must be a radio_beam.Beams object;\" Adam made bugfix on 2022/12/01\n",
    "    mywcs = wcs.WCS(hdul[0].header)\n",
    "    unit = u.Unit(hdul[0].header[\"BUNIT\"])\n",
    "    beam = Beam.from_fits_header(hdul[0].header)\n",
    "    spectrum_in = OneDSpectrum(hdul[0].data, unit = unit, wcs = mywcs, header = hdul[0].header, beam = beam)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}