
from spectral_cube import SpectralCube
import astropy.units as u
from astropy.coordinates import SkyCoord
import numpy as np
import matplotlib.pyplot as p

hi_cube = SpectralCube.read(
    '/media/eric/MyRAID/M33/14B-088/HI/imaging/south_arm_800_1200.image.fits')

hi_cube = hi_cube.to(u.K, hi_cube.beam.jtok_equiv(1420.40575177*u.MHz))

co21_cube = SpectralCube.read('/media/eric/Data_3/M33/IRAM/m33.co21_iram.fits')
co21_cube = co21_cube.spectral_slab(*hi_cube.spectral_extrema)

# Only keep the overlap
lat_extrema = co21_cube.latitude_extrema
long_extrema = co21_cube.longitude_extrema
hi_cube = hi_cube.subcube(xlo=long_extrema[1],
                          xhi=long_extrema[0],
                          ylo=lat_extrema[0],
                          yhi=lat_extrema[1])

# Now find the position in each cube corresponding to a given
# coordinate

posn = SkyCoord("1h33m25.684", "+30d32m17.009", frame='icrs')


def get_closest_posn(posn, spatial_footprint):
    '''
    '''

    spatial_coords = SkyCoord(ra=spatial_footprint[1],
                              dec=spatial_footprint[0],
                              frame='icrs')

    min_posn = spatial_coords.separation(posn).argmin()

    twod_posn = np.unravel_index(min_posn, spatial_footprint[0].shape)

    return twod_posn

hi_posn = get_closest_posn(posn, hi_cube.spatial_coordinate_map)

co21_posn = get_closest_posn(posn, co21_cube.spatial_coordinate_map)

hi_spectrum = hi_cube[:, hi_posn[0], hi_posn[1]]

co21_spectrum = co21_cube[:, co21_posn[0], co21_posn[1]]


p.plot(hi_spectrum.spectral_axis, hi_spectrum.value)
p.plot(co21_spectrum.spectral_axis, co21_spectrum.value)
@@ -873,14 +873,14 @@
 h33m
-25.684
+14.905
 %22, %22
@@ -888,15 +888,15 @@
 30d3
-2
+3
 m17.
-00
+47
 9%22, 
@@ -1526,18 +1526,38 @@
 sn%5B1%5D%5D%0A%0A
-%0Ap
+ax = p.subplot(131)%0Aax
 .plot(hi
@@ -1579,16 +1579,22 @@
 ral_axis
+/1000.
 , hi_spe
@@ -1606,17 +1606,174 @@
 .value)%0A
-p
+ax.set_xlabel('Velocity (km/s)')%0Aax.set_ylabel('HI Surface Brightness (K)')%0Afor tick in ax.get_xticklabels():%0A    tick.set_rotation(45)%0Aax_2 = ax.twinx()%0Aax_2
 .plot(co
@@ -1797,16 +1797,22 @@
 ral_axis
+/1000.
 , co21_s
@@ -1812,21 +1812,659 @@
  co21_spectrum.value
-)
+, 'g')%0Aax_2.set_ylabel('CO(2-1) Surface Brightness (K)')%0A%0Ap.subplot(233)%0Ahi_mom0 = np.nansum(hi_cube.filled_data%5B:%5D.value, axis=0)%0Ahi_slice = %5Bslice(hi_posn%5B0%5D-30, hi_posn%5B0%5D+30),%0A            slice(hi_posn%5B1%5D-30, hi_posn%5B1%5D+30)%5D%0Ap.imshow(hi_mom0%5Bhi_slice%5D, origin='lower')%0Ap.scatter(%5B30%5D, %5B30%5D, c='r', s=40)%0A%0Ap.subplot(236)%0Aco21_mom0 = co21_cube.moment0().value%0Aco21_slice = %5Bslice(co21_posn%5B0%5D-30, co21_posn%5B0%5D+30),%0A              slice(co21_posn%5B1%5D-30, co21_posn%5B1%5D+30)%5D%0Ap.imshow(co21_mom0%5Bco21_slice%5D, origin='lower')%0Ap.scatter(%5B30%5D, %5B30%5D, c='r', s=40)%0A%0Ap.savefig(%22hi_co_overlay_%22+new_coord.to_string('hmsdms').replace(' ', '_')+%22.pdf%22)%0A

Add plotting