import astroML.datasets as datasets

# datasets.

from IPython.display import Image
Image("http://www.sdss3.org/images/dr8_spectro_coverage.png", width=800)

sdss_data = datasets.fetch_imaging_sample()

figure(figsize=(4,4))
plot(sdss_data['ra'], sdss_data['dec'],
     marker='.', linestyle='none', alpha=0.1)
xlim(360, 0)
ylim(-90, 90)

figure(figsize=(12,6))
plot(sdss_data['ra'], sdss_data['dec'],
     marker='.', linestyle='none', alpha=0.05)
title("Source density")
xlim(10.5, -0.5)
ylim(-2, 2)

stars = sdss_data[sdss_data['type'] == 6]
galaxies = sdss_data[sdss_data['type'] == 3]
stars.size, galaxies.size

fig,axes = subplots(1,2,figsize=(12,4),
                    sharex=True,sharey=True)

axes[0].plot(stars['gRaw']-stars['rRaw'], 
             stars['rRaw']-stars['iRaw'],
             alpha=0.02, marker='.', linestyle='none')
axes[0].set_xlabel(r"$g-r$")
axes[0].set_ylabel(r"$r-i$")
axes[0].set_title("Stars")

axes[1].plot(galaxies['gRaw']-galaxies['rRaw'], 
             galaxies['rRaw']-galaxies['iRaw'],
             alpha=0.02, marker='.', linestyle='none')
axes[1].set_xlabel(r"$g-r$")
axes[1].set_title("Galaxies")

axes[0].set_xlim(-0.5, 3)
axes[0].set_ylim(-1, 2)

from astropy.coordinates import ICRS
import astropy.units as u

eq = ICRS(stars['ra']*u.deg, stars['dec']*u.deg)
l,b = eq.galactic.l, eq.galactic.b

plot(l.degree, b.degree,
     marker='.',linestyle='none',alpha=0.05)
xlim(95,117)
ylim(-64,-59)
xlabel('$l$ [deg]')
ylabel('$b$ [deg]')

spectrum = datasets.fetch_sdss_spectrum(2872, 54533, 71)

from scipy.ndimage import gaussian_filter

figure(figsize=(12,6))
plot(spectrum.wavelength(), gaussian_filter(spectrum.spectrum,2), 
     drawstyle='steps', alpha=0.75, color='k')
xlim(min(spectrum.wavelength()), max(spectrum.wavelength()))
xlabel(r"Wavelength [$\AA$]")
ylabel(r"Flux")

qso = datasets.fetch_dr7_quasar()

figure(figsize=(8,5))
n,bins,patches = hist(qso['redshift'], bins=100, 
                      alpha=0.75)
xlabel('Redshift')

segue = datasets.fetch_sdss_sspp()

eq = ICRS(segue['ra']*u.deg, segue['dec']*u.deg)
l,b = eq.galactic.l, eq.galactic.b

l = l.wrap_at(180*u.deg)

figure(figsize=(12,8))
subplot(projection='hammer')
plot(l.radian, b.radian,
     marker='.',linestyle='none',alpha=0.05)
xlabel('$l$ [deg]')
ylabel('$b$ [deg]')

repr(segue.dtype.names)

good = ((segue['radVelErr']/segue['radVel']) < 0.1) & \
       ((segue['FeHErr']/segue['FeH']) < 0.1) & \
       ((segue['alphFeErr']/segue['alphFe']) < 0.1) & \
        (segue['rpsf'] > 15) & (segue['rpsf'] < 19)

good_segue = segue[good]
good_segue.size

from streams.coordinates import vhel_to_vgsr

eq = ICRS(good_segue['ra']*u.deg, good_segue['dec']*u.deg)
l,b = eq.galactic.l, eq.galactic.b
vgsr = vhel_to_vgsr(l, b, good_segue['radVel']*u.km/u.s)

plot(good_segue['FeH'], vgsr.value,
     marker='.',alpha=0.02,linestyle='none')
xlim(-3, 0.1)
ylim(-400, 400)
xlabel(r'[Fe/H]')
ylabel(r'$V_{\rm gsr}$')

plot(good_segue['FeH'], good_segue['alphFe'],
     marker='.',alpha=0.03,linestyle='none')
xlabel(r'[Fe/H]')
ylabel(r'[$\alpha$/Fe]')
xlim(-3,0)
ylim(0.05,0.6)

linear = datasets.fetch_LINEAR_sample()

trg = linear.targets

sesar2013 = ((trg['raLIN'] > 309) | (trg['raLIN'] < 60)) & \
            (np.abs(trg['decLIN']) < 1.23) & \
            (trg['ug'] > 0.75) & (trg['ug'] < 1.45) & \
            (trg['gr'] > -0.25) & (trg['gr'] < 0.4) & \
            (trg['ri'] > -0.2) & (trg['ri'] < 0.2) & \
            (trg['iz'] > -0.3) & (trg['iz'] < 0.3)

sesar2013 = trg[sesar2013]
sesar2013.size

rrlyr = sesar2013[0]
P = 10**rrlyr['LP1']
t,mag,err = linear[rrlyr['objectID']].T
phase = (t / P + 0.1) % 1

figure(figsize=(12,6))
errorbar(phase, mag, err, linestyle='none',
         marker='o', ecolor='#888888')
xlabel("Phase")
ylabel("mag")