import numpy as np
import scipy.interpolate as interp
import matplotlib.pyplot as plt
from astroML.plotting import setup_text_plots
#Lets text in plots use latex
setup_text_plots(usetex=True,fontsize=18)

from astroquery.sdss import SDSS
from astropy import coordinates as coords
import astropy.units as u
#Find a spectrum using astroquery
pos = coords.SkyCoord('120.01976d +45.916684d', frame='icrs')
xid = SDSS.query_region(pos, spectro=True,radius=2*u.arcsec)

/wandajune_home/jakub/PYTHON/loopy/lib/python2.7/site-packages/astroquery/sdss/__init__.py:33: UserWarning: Experimental: SDSS has not yet been refactored to have its API match the rest of astroquery (but it's nearly there).
  warnings.warn("Experimental: SDSS has not yet been refactored to have its API "

#Find matches based on object id assignment from SDSS
#It turns out this object was observed with two different intruments in Sloan
#BOSS and SDSS spectragraphs have different responses so it might be interesting to see how well we
#can recover the flux for the different spectrographs
print(xid)
sp = SDSS.get_spectra(matches=xid)

WARNING: AstropyDeprecationWarning: Since 0.4, config parameter 'astropy.utils.data.REMOTE_TIMEOUT' is deprecated. Use 'astropy.utils.data.conf.remote_timeout' instead. [astropy.config.configuration]
WARNING:astropy:AstropyDeprecationWarning: Since 0.4, config parameter 'astropy.utils.data.REMOTE_TIMEOUT' is deprecated. Use 'astropy.utils.data.conf.remote_timeout' instead.

      ra           dec             objid        ... run2d  instrument
------------- ------------- ------------------- ... ------ ----------
120.019735505 45.9167053402 1237651272957165652 ...     26       SDSS
120.019735505 45.9167053402 1237651272957165652 ... v5_7_0       BOSS

#Create a plot using the two spectra
#loop over spectra and plot
fig, ax = plt.subplots(figsize=(14.,8.5))
for i in np.arange(xid['ra'].size):
    ax.plot(10.**sp[i][1].data['loglam'],sp[i][1].data['flux'],label=xid['instrument'][i])
ax.set_ylabel('Flux [10$^{-17}$ ergs/cm$^2$/s/\AA]')
ax.set_xlabel('Wavelength [\AA]')
ax.legend(loc='upper right')

<matplotlib.legend.Legend at 0x7fcbe7255350>

#manually create a figure
def create_figure():

   # definitions for the axes
   left, width = 0.1, 0.65
   bottom, height = 0.1, 0.65
   bottom_h = left_h = left+width+0.001

   rect_line = [left, bottom, width, height]
   rect_scatterx = [left, bottom_h, width, 0.2]

   # start with a Figure
   plt.figure(figsize=(18,12))
   axline = plt.axes(rect_line)
   axscatter = plt.axes(rect_scatterx)
#return the axis objects
   return axline,axscatter

#add information to subplot
def plotter(x,y1,y2,axline,axscatter):
   from matplotlib.ticker import NullFormatter
#remove default formatting from subplot
   nullfmt   = NullFormatter()         # no labels
   setup_text_plots(fontsize=24,usetex=True)
   # no labels on xaxis of subplot
   axscatter.xaxis.set_major_formatter(nullfmt)
   print((y1-y2)/y2)
#Plot a line of zero
   axscatter.plot(axline.get_xlim(),[0,0],'--',color='gray')
   axscatter.scatter(x,(y1-y2)/y2, color='black')

   axscatter.set_xlim( axline.get_xlim() )
   axscatter.set_ylim([-1,1])
   axscatter.set_ylabel(u'$\delta$ Flux/Flux')

#We see the fluxes are not exactly the same for the two objects especially in the blue end
#Now it might be useful to visualize how the flux at a given wavelength changes 
#for that a more complicated graph is needed.
sdss_spec = sp[0][1].data['flux']
wav0 = sp[0][1].data['loglam']
#interpolate the boss spectrum in to the same wavelength scale as the SDSS
boss_spec = np.interp(wav0,sp[1][1].data['loglam'],sp[1][1].data['flux'])
 
#Create the figures    
axline, axscatter = create_figure()
#add info to the main plot
axline.plot(10.**wav0,sdss_spec,label='SDSS')
axline.plot(10.**wav0,boss_spec,label='BOSS')
axline.set_ylabel('Flux [10$^{-17}$ ergs/cm$^2$/s/\AA]')
axline.set_xlabel('Wavelength [\AA]')
axline.legend(loc='upper right')
#add subplot info
plotter(10.**wav0,boss_spec,sdss_spec,axline,axscatter)

[-0.22289521 -0.0908732   0.04027061 ...,  3.17577691 -0.31567347
 -0.62209734]