read_vanderwel_2012.py.old

# This is obsolete.  Real version is in  /Volumes/Apps_and_Docs/jrrigby1/Python/NB/MRD to Pandas to Seaborn plot.ipynb
from astropy.io import ascii
from astropy.table import Table
from matplotlib import pyplot as plt
import seaborn as sns
import numpy as np

download = False  # Download the catlaog from the ApJ, or use local copy?
if download : 
    file = "http://iopscience.iop.org/0067-0049/203/2/24/suppdata/apjs449744t2_mrt.txt"
else :
    file = "apjs449744t2_mrt.txt"  # MRT version of Table 2 of VanderWel et al 2012
temp_table = ascii.read(file)

print "Have read in van der Wel table as an Astropy table"
table = temp_table.to_pandas()  # Convert from astropy Table to Pandas Data Frame.  Needs astropy 1.2
print "Have converted to a Pandas Data frame"
print table.keys()
print table.info()
print table.head()

good = table[table['Flag'].eq(0)]  # filter out bad Galfit fits
H = good[good['Filter'].str.contains("H")]  # filter just H-band
J = good[good['Filter'].str.contains("J")]  # filter just J-band

small = H[H['ID'].lt(3000)]  # make a small subset of H, for quick plotting
plotted = small  #Set to small to run fast, H to run slow but more accurate
print "Plotting some stuff"
sns.set(font_scale=2)
sns.set_style("white")
ax = sns.kdeplot(plotted.r, plotted.n, n_levels=22)
ax.set_xlim(0,1.5)
ax.set_ylim(0,5)
ax.set_xlabel("effective radius (\")")
ax.set_ylabel("Sersic index")

kpc_per_arcsec = 8.085
# add in s1110            _e(")          n     d(r_e("))        d(n)
s1110_candelized= (2.3/kpc_per_arcsec, 1.0, 0.3/kpc_per_arcsec, 0.4)   # average over all bands
s1110_F160cand = (0.32, 0.68, 0.029, 0.26)   #F160W only, from Table 1 of paper draft
plt.scatter( s1110_candelized[0], s1110_candelized[1], color='k', s=100)
plt.errorbar(s1110_candelized[0],s1110_candelized[1], xerr=s1110_candelized[2],yerr=s1110_candelized[3], ecolor='k')
plt.scatter( s1110_F160cand[0], s1110_F160cand[1], color='b', s=100)
plt.errorbar(s1110_F160cand[0],s1110_F160cand[1], xerr=s1110_F160cand[2],yerr=s1110_F160cand[3], ecolor='b')

# this looked powerful, but i can't get the histograms to bin right, and i need to move on
#sns.jointplot(H.r, H.n, kind="scatter", space=0, color='g', xlim=(0,2), ylim=(0,5), stat_func=None, marginal_kws=dict(bins=15, joint_kws=dict(n_levels=30)))
plt.show()
# only plotting subset, for speed to get graph pretty.  then replace w H for all datapoints.

#For an extra cookie, let's try to plot w thie histograms in the margins
# code borrowed from https://github.com/mwaskom/seaborn/issues/469
g = sns.JointGrid(plotted.r, plotted.n)
g.set_axis_labels("effective radius (\")", "Sersic index")
g.ax_marg_x.hist(plotted.r, bins=np.arange(0, 1.4, 0.05))
g.ax_marg_y.hist(plotted.n, bins=np.arange(0, 5, 0.2), orientation="horizontal")
g.plot_joint(plt.hexbin, gridsize=30, extent=[0, 1.4, 0, 5], cmap="Blues")