A look at GAMA single-sersic colour profiles
Code collected from BD_analysis.ipynb and plot_illustrations.py
In [1]:
%matplotlib inline
import matplotlib
from matplotlib import pyplot as plt
# better-looking plots
plt.rcParams['font.family'] = 'serif'
plt.rcParams['figure.figsize'] = (10.0*1.3, 8*1.3)
plt.rcParams['font.size'] = 18*1.3
In [2]:
import numpy as np
In [3]:
import pandas as pd
#from galapagos_to_pandas import galapagos_to_pandas
## convert the GALAPAGOS data
#galapagos_to_pandas()
## convert the SIGMA data
#galapagos_to_pandas('/home/ppzsb1/projects/gama/qc/raw/StructureCat_SersicExp.fits',
# '/home/ppzsb1/quickdata/StructureCat_SersicExp.h5')
In [4]:
from astropy.io import fits as pyfits
In [5]:
data_all = pyfits.getdata('/Users/spb/Work/projects/MegaMorph/benedetta/gama_verylast.fits', 1)
In [6]:
len(data_all)
Out[6]:
In [7]:
data_all = pd.DataFrame(data_all)
In [8]:
data_all = data_all+0
data_all
Out[8]:
In [9]:
## band information
allbands = list('ugrizYJHK')
#band_wl = pd.Series([3543,4770,6231,7625,9134,10395,12483,16313,22010], index=allbands)
normband = 'r'
bands = list('ugrizYJH')
band_labels = ['${}$'.format(i) for i in bands]
band_wl = pd.Series([3543,4770,6231,7625,9134,10395,12483,16313], index=bands)
#normband = 'Z'
#bands = list('ugriYJHK')
#band_wl = numpy.array([3543,4770,6231,7625,10395,12483,16313,22010])
In [10]:
## restrict to Vulcani selection objects
#data = data_all[(data_all.r_abs < -21.2) & (data_all.r_abs > -50) &
# (data_all.Z < 0.3) & (data_all.N_GALFIT_RFBAND_1 < 7.75)]
#len(data)
In [11]:
## restrict to low-z selection objects
data = data_all[(data_all.r_abs < -19.48) & (data_all.r_abs > -50) &
(data_all.Z < 0.15) & (data_all.N_GALFIT_RFBAND_1 < 7.75)]
len(data)
Out[11]:
In [12]:
## extract magnitudes and use consistent column labels
mag_s = pd.DataFrame(data[['{}_abs'.format(b) for b in list('rugizYJHK')]])
mag_s.columns = list('rugizYJHK')
In [13]:
re_s = pd.DataFrame(data[['RE_GALFIT_RFBAND_{}_corr'.format(i) for i in range(1, 10)]])
re_s.columns = list('rugizYJHK')
n_s = pd.DataFrame(data[['N_GALFIT_RFBAND_{}'.format(i) for i in range(1, 10)]])
n_s.columns = list('rugizYJHK')
In [14]:
from scipy.special import gamma, gammaincinv
from numpy import log, log10, exp, power, pi
class Sersic:
# currently doesn't handle uncertainties
def __init__(self, mag, re, n, ar=1.0, pa=0.0,
mag_err=None, re_err=None, n_err=None, ar_err=None, pa_err=None, xc_err=None, yc_err=None):
self.mag = mag
self.re = re
self.n = n
self.ar = ar
self.pa = pa
self.mag_err = mag_err
self.re_err = re_err
self.n_err = n_err
self.ar_err = ar_err
self.pa_err = pa_err
def __call__(self, r):
return self.mu_r(r)
def mu_r(self, r):
# Returns the surface brightess at specified major axis radius,
# within annular ellipses corresponding to the shape of each component individualy
# Taking, e.g. colours, this currently assumes major axes of components align
# to be more generally correct need to account for AR, PA, XC, YC,
# and either select specific vector, or properly compute azimuthal average
mag = self.mag
re = self.re
n = self.n
bn = self.bn()
mue = mag + 5.0*log10(re) + 2.5*log10(2.0*pi*n*gamma(2.0*n)*exp(bn)/power(bn, 2.0*n))
mu = mue + 2.5 * bn / log(10) * (power(r/re, 1.0/n) - 1.0)
return mu
def bn(self):
return gammaincinv(2.0*self.n, 0.5)
# These need testing
def I_el(self, r_m, ar_m, pa_m=0):
return quad(self.I_el_theta, 0, 2*pi, args=(r_m, ar_m, pa_m))[0] / (2*pi)
def mu_el(self, r_m, ar_m, pa_m=0):
return -2.5*numpy.log10(self.I_el(r_m, ar_m, pa_m))
def mu_el_theta(self, theta, r_m, ar_m, pa_m=0):
x = r_m * numpy.cos(theta - pa_m)
y = ar_m * r_m * numpy.sin(theta - pa_m)
r_c = numpy.sqrt(x**2 + self.ar**2 * y**2)
return self.mu_r(r_c)
def I_el_theta(self, theta, r_m, ar_m, pa_m=0):
return 10**(-0.4*self.mu_el_theta(theta, r_m, ar_m, pa_m))
In [15]:
def make_funcs(ids):
func = {}
for i in ids:
func[i] = {}
remax = 0
for band in bands:
func[i][band] = Sersic(mag_s[band].iloc[i], re_s[band].iloc[i]/re_s[normband].iloc[i], n_s[band].iloc[i])
return func
In [16]:
f = make_funcs(np.arange(len(data)))
In [17]:
def plot_colour_profile(prof, band1='g', band2='r', *args, **kwargs):
lograd = np.arange(-2, log10(5)+0.1, 0.1)
rad = 10**lograd
b1 = prof[band1](rad)
b2 = prof[band2](rad)
col = b1 - b2
plt.plot(rad, col, *args, **kwargs)
plt.semilogx()
plt.gca().set_xticks([0.01, 0.03, 0.1, 0.3, 1.0, 3.0])
plt.gca().xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
plt.xlim((0.01, 5))
plt.ylim((-2, 5))
In [18]:
from scipy.odr import *
def gradfit(rad, col):
# use orthogonal regression a la La Barbera
def f(B, x):
return B[0]*x + B[1]
linear = Model(f)
mydata = RealData(rad, col, sx=0.00001, sy=0.01)
myodr = ODR(mydata, linear, beta0=[0., col.mean()])
myoutput = myodr.run()
return myoutput.beta
In [19]:
def colour_gradient(prof, band1='g', band2='r', rad1=0.1, rad2=1.0, fit=True):
if not fit:
rad = [rad1, rad2]
else:
rad = np.arange(rad1, rad2+0.00001, 0.01)
b1 = prof[band1](rad)
b2 = prof[band2](rad)
col = b1 - b2
if not fit:
grad = col[1] - col[0]
else:
grad = gradfit(rad, col)
return grad
In [20]:
def plot_colour_gradient_distribs(profs, hsel=None, lsel=None, fit=False):
fig, ax = plt.subplots(len(bands[2:]), 1, sharex=True, sharey=True)
for i, b in enumerate(bands[2:]):
g = np.array([colour_gradient(profs[id], 'g', b, fit=fit) for id in profs.keys()])
g = g[~np.isnan(g)].clip(-4, 4)
w = 1.0/len(g)
label = '${}\quad all:\quad \mu={:.2f}\quad \sigma={:.2f}$'.format(b, g.mean(), g.std())
ax[i].hist(g, range=(-4, 4), bins=100, histtype='step', label=label,
weights=np.ones_like(g)*w)
if hsel is not None:
g = np.array([colour_gradient(profs[id], 'g', b, fit=fit) for id in profs.keys() if hsel[id]])
g = g[~np.isnan(g)].clip(-4, 4)
label = '${}\ high-n:\quad \mu={:.2f}\quad \sigma={:.2f}$'.format(b, g.mean(), g.std())
ax[i].hist(g, range=(-4, 4), bins=100, histtype='step', label=label,
weights=np.ones_like(g)*w)
if lsel is not None:
g = np.array([colour_gradient(profs[id], 'g', b, fit=fit) for id in profs.keys() if lsel[id]])
g = g[~np.isnan(g)].clip(-4, 4)
label = '${}\ \ \ low-n:\quad \mu={:.2f}\quad \sigma={:.2f}$'.format(b, g.mean(), g.std())
ax[i].hist(g, range=(-4, 4), bins=100, histtype='step', label=label,
weights=np.ones_like(g)*w)
ax[i].legend(fontsize='x-small',loc='upper right')
plt.xlim((-2.5, 3.5))
plt.ylim((0, 0.16))
ax[0].set_yticks((0.0, 0.1))
ax[-1].set_xlabel('$\\nabla(g-x)$')
fig.text(0.05, 0.5, 'norm. freq.', va='center', rotation='vertical')
plt.subplots_adjust(hspace=0)
In [21]:
plot_colour_gradient_distribs(f, hsel=(n_s.r >= 2.5).values, lsel=(n_s.r < 2.5).values)
In [22]:
ids = np.random.choice(len(data), 2000, replace=False)
In [23]:
coltot = (mag_s.u - mag_s.r)
plt.subplot(211)
for id in ids:
if n_s.r.iloc[id] < 2.5:
if coltot.iloc[id] > 2.1:
c = 'r'
elif coltot.iloc[id] > 1.6:
c = 'g'
else:
c = 'b'
plot_colour_profile(f[id], 'g', 'r', c=c, alpha=0.05)
plot_colour_profile(f[id], 'g', 'J', c=c, alpha=0.05)
plt.hlines([1.0], 0.01, 5.0)
plt.ylim((-0.5, 5))
plt.subplot(212)
for id in ids:
if n_s.r.iloc[id] > 2.5:
if coltot.iloc[id] > 2.1:
c = 'r'
elif coltot.iloc[id] > 1.6:
c = 'g'
else:
c = 'b'
plot_colour_profile(f[id], 'g', 'r', c=c, alpha=0.05)
plot_colour_profile(f[id], 'g', 'J', c=c, alpha=0.05)
plt.hlines([1.0], 0.01, 5.0)
_=plt.ylim((-0.5, 5))
In [24]:
def colour_at_rad(profs, rad=1.0, band1='g', band2='r'):
return np.array([f[id][band1](rad) - f[id][band2](rad) for id in f.keys()])
In [25]:
def colour_at_rad_hist(profs, band1='g', band2='r'):
col1re = colour_at_rad(profs, 1.0, band1, band2)
col0p1re = colour_at_rad(profs, 0.1, band1, band2)
col2re = colour_at_rad(profs, 2.0, band1, band2)
col3re = colour_at_rad(profs, 3.0, band1, band2)
coltot = (mag_s[band1] - mag_s[band2]).values
plt.hist((coltot, col0p1re, col1re, col2re, col3re), range=(-1, 4), bins=50, histtype='step',
label=('total', '0.1Re', '1Re', '2Re', '3Re'), normed=True)
plt.legend()
plt.xlabel('${}-{}$'.format(band1, band2))
plt.ylabel('norm. freq.')
In [26]:
colour_at_rad_hist(f, 'g', 'r')
In [27]:
from scipy.stats import binned_statistic
In [28]:
def plot_colour_gradient_vs_mag(profs, ax=None, sel=None, fit=False,
labellines=False, labeloffset={}, label=None):
if ax is None:
fig, ax = plt.subplots()
for i, b in enumerate(bands[2:]):
if sel is None:
g = np.array([colour_gradient(profs[id], 'g', b, fit=fit) for id in profs.keys()])
m = mag_s['r']
else:
g = np.array([colour_gradient(profs[id], 'g', b, fit=fit) for id in profs.keys() if sel[id]])
m = mag_s['r'][sel]
ok = ~np.isnan(g) & (g > -4) & (g < 4)
g = g[ok]
m = m[ok]
# HALF-MAGNITUDE BINS
counts, edges, bins = binned_statistic(m, g, bins=8, statistic='count', range=(-23.5,-19.5))
meds, edges, bins = binned_statistic(m, g, bins=8, statistic='median', range=(-23.5,-19.5))
centres = (edges[:-1]+edges[1:])/2.0
ok = counts > 10
centres = centres[ok]
meds = meds[ok]
#offset = -0.1*(i+1) - np.median(meds)
offset = 0 # NO OFFSETS
meds += offset
g += offset
col = plt.cm.jet((0.8*i/5.0 + 0.1))
ax.plot(m, g, '.', alpha=0.1, c=col, zorder=1)
ax.plot(centres, meds, c='w', lw=6, zorder=2)
ax.plot(centres, meds, c=col, lw=2, zorder=3)
if labellines:
y = meds[0]
if labeloffset.get(b) is not None:
y += labeloffset[b]
ax.text(centres[0]-0.1, y, 'g-{}'.format(b), horizontalalignment='right', color=col,
fontsize='x-small')
if label is not None:
ax.text(0.05, 0.95, label, verticalalignment='top',
transform=ax.transAxes, fontsize='small')
In [29]:
# all types
plot_colour_gradient_vs_mag(f, labellines=True, label='all galaxies')
plt.xlim((-23.99, -19.51))
plt.ylim((-1.75, 0.49))
plt.xlabel('$M_r$')
plt.ylabel('$\\nabla_{g-x}$')
plt.locator_params(nbins=6)
In [30]:
highn = (n_s.r >= 2.5).values
lown = ~highn
red = (coltot >= 2.1).values
green = (coltot < 2.1).values & (coltot >= 1.6).values
blue = (coltot < 1.6).values
fig, ax = plt.subplots(2, 3, sharex=True, sharey=True)
plot_colour_gradient_vs_mag(f, ax[0,0], sel=red & highn, label='red, high-$n$')
plot_colour_gradient_vs_mag(f, ax[0,1], sel=green & highn, label='green, high-$n$')
plot_colour_gradient_vs_mag(f, ax[0,2], sel=blue & highn, label='blue, high-$n$')
plot_colour_gradient_vs_mag(f, ax[1,0], sel=red & lown, label='red, low-$n$')
plot_colour_gradient_vs_mag(f, ax[1,1], sel=green & lown, label='green, low-$n$',
labellines=True, labeloffset={'z':0.01, 'H':-0.05})
plot_colour_gradient_vs_mag(f, ax[1,2], sel=blue & lown, label='blue, low-$n$')
plt.xlim((-23.49, -19.51))
plt.ylim((-1.75, 0.49))
ax[1,1].set_xlabel('$M_r$')
fig.text(0.02, 0.5, '$\\nabla_{g-x}$', va='center', rotation='vertical')
plt.locator_params(nbins=6)
plt.subplots_adjust(hspace=0, wspace=0)