<jupyter_start><jupyter_text>Figure 4: H-Alpha Map Requirerd python packages:* numpy* matplotlib* scipy* [seaborn](https://seaborn.pydata.org/index.html)* [whampy](https://whampy.readthedocs.io/en/latest/index.html)* [astropy](https://docs.astropy.org/en/stable/index.html)* [bettermoments](https://bettermoments.readthedocs.io/en/latest/)* [dustmaps](https://dustmaps.readthedocs.io/en/latest/)* [extinction](https://extinction.readthedocs.io/en/latest/)<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
%matplotlib notebook

from whampy import SkySurvey

from scipy.interpolate import interp1d


import seaborn as sns
pal = sns.color_palette("colorblind")


from astropy.table import Table
import astropy.units as u
import glob

from bettermoments.collapse_cube import collapse_width
from bettermoments.collapse_cube import collapse_quadratic
from bettermoments.collapse_cube import collapse_zeroth


from astropy.coordinates import SkyCoord

from mpl_toolkits.axes_grid1.inset_locator import inset_axes

from dustmaps.bayestar import BayestarQuery
from extinction import fm07 as extinction_law<jupyter_output><empty_output><jupyter_text>Read in Data<jupyter_code>vel_range = [-270, -200] * u.km/u.s

ha_pds = SkySurvey(filename = "HA_PDS_BLOCKS_R_WITHIP.sav", idl_var = "ha_pds_ae_r")
ha_pds["VELOCITY"][ha_pds["VELOCITY"] < -310] = np.nan

# Remove Problematic pointings
vel_mask = np.array([(ha_pds[ell]["VELOCITY"] < vel_range[1].value) & 
            (ha_pds[ell]["VELOCITY"] > vel_range[0].value) 
            for ell in range(len(ha_pds))])

pos_mask = np.array([ha_pds[ell]["DATA"] > 0
            for ell in range(len(ha_pds))])

remove_row_inds = np.array([ell for ell,(vm,pm) in enumerate(zip(vel_mask, pos_mask)) if np.sum(vm&pm) < 5])
if len(remove_row_inds) > 0:
    ha_pds.remove_rows(remove_row_inds)
def smooth(data, width = 3):
    from scipy.signal import savgol_filter
    data = savgol_filter(data, width, polyorder=2,
                         axis=0, mode = "wrap")
    return data<jupyter_output><empty_output><jupyter_text>Intensity Measurements<jupyter_code>vel_range = [-270, -200] * u.km/u.s



vel_mask = np.array([(ha_pds[ell]["VELOCITY"] < vel_range[1].value) & 
            (ha_pds[ell]["VELOCITY"] > vel_range[0].value) 
            for ell in range(len(ha_pds))])

pos_mask = np.array([ha_pds[ell]["DATA"] > 0
            for ell in range(len(ha_pds))])

remove_row_inds = np.array([ell for ell,(vm,pm) in enumerate(zip(vel_mask, pos_mask)) if np.sum(vm&pm) < 5])
if len(remove_row_inds) > 0:
    ha_pds.remove_rows(remove_row_inds)

vel_mask = np.array([(ha_pds[ell]["VELOCITY"] < vel_range[1].value) & 
            (ha_pds[ell]["VELOCITY"] > vel_range[0].value) 
            for ell in range(len(ha_pds))])

pos_mask = np.array([ha_pds[ell]["DATA"] > 0
            for ell in range(len(ha_pds))])




ha_zeroth = np.vstack([collapse_zeroth(ha_pds[ell]["VELOCITY"][vel_mask[ell] & pos_mask[ell]], 
                             ha_pds[ell]["DATA"][vel_mask[ell] & pos_mask[ell]], 
                             rms = 0.015, 
                                       threshold = 0.0) 
             for ell in range(len(ha_pds))])

ha_quadratic = np.vstack([collapse_quadratic(ha_pds[ell]["VELOCITY"][vel_mask[ell] & pos_mask[ell]], 
                                   smooth(ha_pds[ell]["DATA"][vel_mask[ell] & pos_mask[ell]], width = 3), 
                                   rms = 0.015)
                for ell in range(len(ha_pds))])

mom0_value = ha_zeroth[:,0]
mom0_std = ha_zeroth[:,1]

sig1 = mom0_value > 1 * mom0_std
sig2 = mom0_value > 2 * mom0_std
sig3 = mom0_value > 3 * mom0_std

ha_vel = ha_quadratic[:,0]
ha_vel_err = ha_quadratic[:,1]<jupyter_output><empty_output><jupyter_text>Prepare Extintion Map<jupyter_code>bayestar = BayestarQuery()
# Try Deredden
distance = 6.5 * u.kpc

pds_map_coords = ha_pds.get_SkyCoord(distance = np.full(len(ha_pds), fill_value = distance.value) * u.kpc)


Av_bayestar = 2.742 * bayestar(pds_map_coords)
wave_ha = np.array([6562.8])
wave_nii = np.array([6584.])
A_V_to_A_ha = extinction_law(wave_ha, 1.)
A_V_to_A_nii = extinction_law(wave_nii, 1.)<jupyter_output><empty_output><jupyter_text>Plot Figure 4: Done in Two Parts<jupyter_code>fig = plt.figure(figsize = (7.3,4))

ax = fig.add_subplot(121)
ax2 = fig.add_subplot(122, sharey = ax)

plt.subplots_adjust(hspace = 0.0, wspace = 0)

cmap = "Reds"

s, fig = ha_pds[sig1].intensity_map(s_factor = 6, 
                                                        c = mom0_value[sig1], 
                                    cmap = cmap, 
                                    vmin = .1, vmax = 1., fig = fig, ax = ax, 
                                        alpha = 1, zorder = 2, return_sc = True)
pl_cl = sns.mpl_palette(cmap, 256)
ax.set_facecolor(pl_cl[0])
# ax.set_alpha(0.1)
ax.scatter(ha_pds["GAL-LON"][~sig1], ha_pds["GAL-LAT"][~sig1], 
           facecolors='none', edgecolors = "k", 
           s = 140, alpha = 0.5, zorder = 1)

ax.scatter(ha_pds["GAL-LON"][sig3], ha_pds["GAL-LAT"][sig3], 
           facecolors='none', edgecolors = "b",linestyles = "-", 
           s = 140, alpha = 0.9, zorder = 2, label = r"$>3 \sigma$", linewidths = 1.5)

ax.scatter(ha_pds["GAL-LON"][sig2 & ~sig3], ha_pds["GAL-LAT"][sig2 & ~sig3], 
           facecolors='none', edgecolors = "b", linestyles = ":",
           s = 140, alpha = 0.9, zorder = 2, label = r"$>2 \sigma$", linewidths = 2)


ax.scatter(10.4, 11.2, marker = "*", color = pal[0], s= 400, label = "PDS 456", alpha = 0.7, zorder = 4)
ax.invert_xaxis()
xlim = ax.get_xlim()
ylim = ax.get_ylim()

xx = np.linspace(xlim[0]-5, xlim[1]+5, 500)
yy = np.linspace(ylim[0]-5, ylim[1]+5, 500)

cgrid_lon, cgrid_lat = np.meshgrid(xx, yy)

grid_coords = SkyCoord(l = cgrid_lon * u.deg, 
                       b = cgrid_lat * u.deg, 
                       distance = np.full_like(cgrid_lat, 
                                          fill_value = distance.value), 
                       frame = "galactic")

Av_grid = 2.742 * bayestar(grid_coords)

Av_med_grid = np.median(Av_grid, axis = 2)

# im = ax.imshow(Av_med_grid, 
#                 cmap = "Greys", 
#                 extent = [xlim[1]+5, xlim[0]-5, ylim[0]-5, ylim[1]+5], 
#                 zorder = 0, 
#                 alpha = 1, 
#                 origin = "lower", 
#                vmin = 0.5, 
#                vmax = 7)

ax.text(16.7, 1.5, 
        "Observed", 
        fontsize = 14, ha = "center", va = "bottom", color = 'k', fontweight = "bold")
ax.text(16.7, 0.5, 
        r"H$\mathbf{\alpha}$", 
        fontsize = 14, ha = "center", va = "bottom", color = 'k', fontweight = "bold")

# ax.text(18, 14, 
#         r"($A_V$ greyscale)", 
#         fontsize = 14, ha = "center", va = "bottom")

# ax.set_xlabel(None)
# ax.set_xticklabels([])

# Bubble Boundaries

north_bubble_l = [-9.9, 0, 12, 16, 18.4]
north_bubble_b = [ 5, 0,5, 10, 15]

boundary_int = interp1d(north_bubble_l, north_bubble_b, bounds_error=False, kind = "quadratic", fill_value = 20)
bcol = pal[2]
xx = np.linspace(xlim[0], xlim[1], 100)
yy = boundary_int(xx)
ax.plot(xx, yy, lw = 3, ls = "--", color = bcol, alpha = 0.9)
ax2.plot(xx, yy, lw = 3, ls = "--", color = bcol, label = "Su+ (2010) Boundary", alpha = 0.9)

xlim = ax.set_xlim(xlim)
ylim = ax.set_ylim(ylim)
ax.invert_xaxis()



lg = ax.legend(fontsize = 12, loc = 2, facecolor = "w", framealpha = 0.5)
lg2 = ax2.legend(fontsize = 12, loc = 2, facecolor = "w", framealpha = 0.8)
ax.set_aspect("equal")




ha_intrinsic_grid = np.full_like(Av_med_grid, fill_value = 1.0)
ha_reddened_grid = ha_intrinsic_grid * 10**(-0.4 * A_V_to_A_ha  * Av_med_grid)

im2 = ax2.imshow(Av_med_grid, 
                cmap = "Greys", 
                extent = [xlim[1]+5, xlim[0]-5, ylim[0]-5, ylim[1]+5], 
                zorder = 0, 
                alpha = 1, 
                origin = "lower", 
               norm = None, vmin = 0, vmax = 5)

ax2.scatter(ha_pds["GAL-LON"][~sig2], ha_pds["GAL-LAT"][~sig2], 
           facecolors='none', edgecolors = "k", 
           s = 140, alpha = 0.5, zorder = 1)

ax2.scatter(ha_pds["GAL-LON"][sig3], ha_pds["GAL-LAT"][sig3], 
           facecolors='none', edgecolors = "b", linestyles = "-", 
           s = 140, alpha = 0.9, zorder = 2, linewidths = 1.5)

ax2.scatter(ha_pds["GAL-LON"][sig2 & ~sig3], ha_pds["GAL-LAT"][sig2 & ~sig3], 
           facecolors='none', edgecolors = "b", linestyles = ":",
           s = 140, alpha = 0.9, zorder = 2, linewidths = 2)

ax2.scatter(10.4, 11.2, marker = "*", color = pal[0], s= 400, label = "PDS 456", alpha = 1)

xlim = ax2.set_xlim(xlim)
ylim = ax2.set_ylim(ylim)
ax2.invert_xaxis()

ax2.text(16.7, 1.5, 
        r"$\mathbf{{A_V}}$ at", 
        fontsize = 14, ha = "center", va = "bottom", color = "w", fontweight = "bold")
ax2.text(16.7, .3, 
        "D = 6.5 kpc", 
        fontsize = 14, ha = "center", va = "bottom", color = "w", fontweight = "bold")


# plt.colorbar(im)
ax2.yaxis.tick_right()

ax2.set_aspect("equal")

ax2.set_xlabel("Galactic Longitude (deg)", fontsize = 12)
ax.set_xlabel("Galactic Longitude (deg)", fontsize = 12)
ax2.yaxis.set_label_position("right")
ax2.set_ylabel("Galactic Latitude (deg)", fontsize = 12)
ax.set_ylabel(None)





plt.tight_layout()

plt.subplots_adjust(left = 0.12)

cax = fig.add_axes([0.065, 0.165, .013, 0.79])


cb = plt.colorbar(s, cax = cax, orientation = "vertical", extend = "both")
cax.tick_params(axis='y', rotation=90)
cax_label = cax.set_ylabel(r"$\log_{10}(H\alpha/R)$", fontsize = 12)
cax.yaxis.tick_left()
cax.yaxis.set_label_position("left")

# plt.savefig("/Users/dk/Desktop/H-Alpha_MapArea.png", dpi = 300, transparent = 300)
# plt.savefig("/Users/dk/Desktop/H-Alpha_MapArea.svg", transparent = 300)
peak_mask = ha_vel_err < 12



fig = plt.figure(figsize = (7.3,4))
ax = fig.add_subplot(212)

lat_jitter = 0.05 * np.random.randn(len(ha_pds))
lon_jitter = 0.05 * np.random.randn(len(ha_pds))

ax.scatter(ha_pds["GAL-LAT"][peak_mask] + lat_jitter[peak_mask], 
           ha_vel[peak_mask], 
           c = mom0_value[peak_mask], 
           vmin = .1, vmax = 1., alpha = 0.9, cmap = cmap, norm = LogNorm(), s = 150)
ax.scatter(ha_pds["GAL-LAT"][(sig3) & (peak_mask)] + lat_jitter[(sig3) & (peak_mask)], 
           ha_vel[(sig3) & (peak_mask)],
           vmin = .1, vmax = 2.5, alpha = 0.9, linewidths = 1.5, 
           facecolors = "none", edgecolors = "b", 
           norm = LogNorm(), s = 150)
ax.scatter(ha_pds["GAL-LAT"][(sig2) & (peak_mask) & (~sig3)] + lat_jitter[(sig2) & (peak_mask) & (~sig3)], 
           ha_vel[(sig2) & (peak_mask) & (~sig3)],
           vmin = .1, vmax = 2.5, alpha = 0.9, linewidths = 1.5,
           facecolors = "none", edgecolors = "b", linestyles = ":", 
           norm = LogNorm(), s = 150)

    
plotter = sns.regression._RegressionPlotter(ha_pds["GAL-LAT"][peak_mask], 
                                            ha_vel[peak_mask], robust = True)
sns.regplot(ha_pds["GAL-LAT"][peak_mask], 
                                            ha_vel[peak_mask], robust = True, 
            line_kws = {"color":pal[0], "alpha":0.5, "zorder":0, "lw":2},scatter = False, ax = ax)

grid, yhat, err_bands = plotter.fit_regression(ax)
g1,y1,e1 = grid, yhat, err_bands 
slope = (yhat[-1] - yhat[0]) / (grid[-1] - grid[0])
intercept = yhat[0] - slope * grid[0]
err_slope_0 = (err_bands[0][-1] - err_bands[1][0]) / (grid[-1]-grid[0])
err_slope_1 = (err_bands[1][-1] - err_bands[0][0]) / (grid[-1]-grid[0])

# ax.set_ylim(-250, -235)

ax.set_xlabel("Galactic Latitude (deg)", fontsize = 12)

# ax.set_ylabel("Centroid Velocity (km/s)", fontsize = 12)


ax.text(11, -250, r"Slope = ${0:.2f}^{{+{2:.2f}}}_{{-{1:.2f}}}$ $\frac{{km/s}}{{deg}}$ ".format(slope, 
                                                                       err_slope_1 - slope, slope - err_slope_0), 
        fontsize = 14, va = "center", ha = "left")




ax = fig.add_subplot(211)


sc = ax.scatter(ha_pds["GAL-LON"][peak_mask] + lon_jitter[peak_mask], 
                ha_vel[peak_mask], 
                c = mom0_value[peak_mask], 
           vmin = .1, vmax = 1., alpha = 0.9, cmap = cmap, norm = LogNorm(), s = 150)
ax.scatter(ha_pds["GAL-LON"][(sig3) & (peak_mask)] + lon_jitter[(sig3) & (peak_mask)], 
           ha_vel[(sig3) & (peak_mask)], 
           alpha = 0.9, 
           facecolors = "none", linewidths = 1.5,
           edgecolors = "b", 
           label = r"$> 3 \sigma$", s = 150)
ax.scatter(ha_pds["GAL-LON"][(sig2) & (peak_mask) & (~sig3)] + lon_jitter[(sig2) & (peak_mask) & (~sig3)], 
           ha_vel[(sig2) & (peak_mask) & (~sig3)], 
           alpha = 0.9, 
           facecolors = "none", linewidths = 1.5, 
           edgecolors = "b", linestyles = ":",
           label = r"$> 2 \sigma$", s = 150)

    
plotter = sns.regression._RegressionPlotter(ha_pds["GAL-LON"][peak_mask], 
                                            ha_vel[peak_mask], robust = True)
sns.regplot(ha_pds["GAL-LON"][peak_mask], 
                                            ha_vel[peak_mask], robust = True, 
            line_kws = {"color":pal[0], "alpha":0.5, "zorder":0, "lw":2},
            scatter = False, ax =ax)

grid, yhat, err_bands = plotter.fit_regression(ax)
slope = (yhat[-1] - yhat[0]) / (grid[-1] - grid[0])
intercept = yhat[0] - slope * grid[0]
err_slope_0 = (err_bands[0][-1] - err_bands[1][0]) / (grid[-1]-grid[0])
err_slope_1 = (err_bands[1][-1] - err_bands[0][0]) / (grid[-1]-grid[0])


# ax.yaxis.tick_right()

ax.set_xlabel("Galactic Longitude (deg)", fontsize = 12)
# ax.set_ylabel("Centroid Velocity (km/s)", fontsize = 12)

ax.text(20, -250, r"Slope = ${0:.2f}^{{+{1:.2f}}}_{{-{2:.2f}}}$ $\frac{{km/s}}{{deg}}$ ".format(slope, 
                                                                       slope-err_slope_0, err_slope_1-slope), 
        fontsize = 14, ha = "left", va = "center")

ax.invert_xaxis()


fig.text(0, 0.55, "Centroid Velocity (km/s)", fontsize = 12, ha = 'left', va = "center", rotation = 90)


plt.tight_layout()
plt.subplots_adjust(left = 0.095, right = 0.98)
# plt.savefig("../Figures/Velocity_Gradient_LatLon.png", dpi = 300, transparent = True)
# plt.savefig("../Figures/Velocity_Gradient_LatLon.svg", transparent = True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17425464749210803, 'fraction_numerical': 0.041388986320589266, 'mean_word_length': 2.8046437149719776, 'pattern_counts': {'":': 12, '<': 28, '<?xml version=': 0, '>': 36, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '158224', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6165, 'n_tokens_neox': 5628, 'n_words': 1018, 'path': 'Figure_4.ipynb', 'repo_name': 'Deech08/WHAM_PDS456', 'size': '14255'}