Samples of Python and SQL code written for my graduate school research.

ApogeeQueries.sql

-- Get parameters for all stars with [Fe/H] < -2 with no BAD FLAGS set:
-- You can also select a subset of the stars based on their
-- properties. This example finds a set of metal-poor stars, without any
-- flags set indicating that the observations or analysis is bad.

SELECT TOP 300
   s.apogee_id, s.ra, s.dec, s.glon, s.glat,
   s.vhelio_avg, s.vscatter,
   a.teff, a.logg, a.param_m_h, a.param_alpha_m,
   dbo.fApogeeAspcapFlagN(a.aspxcapflag),
   dbo.fApogeeStarFlagN(s.starflag)
FROM apogeeStar s
   JOIN aspcapStar a on a.apstar_id = s.apstar_id
WHERE (a.aspxcapflag & dbo.fApogeeAspcapFlag('STAR_BAD')) = 0
   AND a.teff > 0 AND a.param_m_h < -2;
   

-- Get individual RVs from individual visits, the ASPCAP parameters for the combined
spectra for stars which have more than 6 visits:
-- Each star is visited several times, and in some case many times, in
-- order to build up signal-to-noise and to detect radial velocity
-- variations. The information about each visit to each star is in the
-- apogeeVisit table. One could join this table with apogeeStar on
-- apogee_id in order to literally find all visits to each star. However,
-- in this example we are interested in just finding those visits that
-- actually contributed to each combined spectrum. In this case, bad visits
-- are excluded and commissioning data and survey data are kept separate
-- (not combined). To find these stars, one may use the apogeeStarVisit
-- table in CAS, or the array visit_pk which exists for each star in the
-- allStar file. Alternatively, if you wanted to find all visits to a
-- particular star, one could replace in the code below apogeeStarVisit
-- with apogeeStarAllVisit and visit_pk with all_visit_pk.

SELECT TOP 300
visit.*, aspcap.teff, aspcap.logg, aspcap.param_m_h
FROM apogeeVisit visit
   JOIN apogeeStarVisit starvisit ON visit.visit_id = starvisit.visit_id
   JOIN aspcapStar aspcap ON aspcap.apstar_id = starvisit.apstar_id
   JOIN apogeeStar star ON star.apstar_id = starvisit.apstar_id
WHERE (aspcap.aspxcapflag & dbo.fApogeeAspcapFlag('STAR_BAD')) = 0 and aspcap.teff > 0
    AND (star.apogee_target1 & dbo.fApogeeTarget1('APOGEE_LONG')) > 0
    AND star.nvisits > 6
ORDER BY visit.apogee_id

RadialBinAbundances.py

from astropy.io import fits
from astropy.coordinates import SkyCoord as sc
from astropy.coordinates import Longitude as lon
from astropy.coordinates import Latitude as lat
from astropy.coordinates import Distance as dis
from math import cos, sin, sqrt, radians
import matplotlib.colors as colors
import matplotlib.cm as cm
from sdss.apogee import apload as ap


hdudata = fits.open('/home/regulus/tpicard/data/allStar+-t30e.2.fits')
table = hdudata[1].data

hdulist = ap.allStar()
e_table = hdulist[3].data
e_list = e_table[0][1].upper()

# get chemical abundances and other stellar parameters 
e_abun = []
for e in range(len(e_list)):
  element = e_list[e]
  if (element != 'FE') and 
     (element != 'CE') and 
     (element != 'GE') and 
     (element != 'RB') and 
     (element != 'Y') and 
     (element != 'ND'):
    e_abun.append(table[element+'_FE'])
FEH = table['FE_H']  # metal abundances
teff = table['TEFF'] # effective temperatures
ASPCAPFLAG = table['ASPCAPFLAG'] # data insights
logg = table['LOGG'] # surface gravities
SNR = table['SNR']   # signal-to-noise ratios
diso = table['diso'][:,1] / 1000 # kpc
l = table['GLON']    # galactic longitudes
b = table['GLAT']    # galactic latitudes

# filter out "bad" stars
good = np.where((ASPCAPFLAG & 2**23 == 0) & 
                (teff <= 5500) & 
                (teff >= 3500) & 
                (FEH >= -5) &
                (logg > 1) & 
                (logg < 3.8) & 
                (SNR > 80) & 
                (diso > 1) & 
                (diso != np.nan))[0] 
for e in range(len(e_abun)):
  e_abun[e] = e_abun[e][good]

# map "good" stars using spherical galactic coordinates
coordinates = sc(lon(l[good], unit='deg'), 
                 lat(b[good], unit='deg'), 
                 dis(diso[good], unit='kpc'), 
                 frame='galactic')
z = np.array(abs(coordinates.galactocentric.cartesian.z))
R = np.array(coordinates.galactocentric.spherical.distance)

# put stellar parameters into radial bins
Rbins, FEHbins, e_abunbins, teffbins = ([],[],[],[])
edge = 3
for bin in range(6):
  Ri = np.where((ASPCAPFLAG[good] & 2**23 == 0) & 
                (teff[good] <= 5500) & 
                (teff[good] >= 3500) & 
                (FEH[good] >= -5) &
                (R > edge) & 
                (R < edge + 2) & 
                (logg[good] > 1) & 
                (logg[good] < 3.8) & 
                (SNR[good] > 80) & 
                (diso[good] > 1) & 
                (diso[good] != np.nan))[0]
  Rbins.append(R[Ri])
  FEHbins.append(FEH[good][Ri])
  teffbins.append(teff[good][Ri])
  edge += 2
  for e in range(len(e_abun)):
    e_abunbins.append(e_abun[e][Ri])
edge = 3
Rbins = np.array(Rbins)
FEHbins = np.array(FEHbins)
e_abunbins = np.array(e_abunbins)

# make a plot for each radial bin with a subplot for each element
e_labels = ["C", "CI", "N", "O", "Na", "Mg", "Al", "Si", "P", "S",
	        "K", "Ca", "Ti", "TiII", "V", "Cr", "Mn", "Co", "Ni", "Cu"]
row, col, i = (0,)*3
for bin in range(6):
  f, axarr = plt.subplots(4, 5)
  for e in range(len(e_abunbins[bin])):
    if (row == 4):
      row = 0
      break
    if (e % 5 == 0) and (e != 0):
      row += 1
      col = 0
    if (row < 4) and (col < 5):
      axarr[row, col].scatter(FEHbins[bin], e_abunbins[i], 
                              c=teffbins[bin], vmin=3500, vmax=5500)
      axarr[row, col].set_xlim([-2.25, 0.75])
      axarr[row, col].set_ylim([-4, 2])
      axarr[row, col].set_xticks(axarr[row, col].get_xticks()[1:6])
      axarr[row, col].set_yticks(axarr[row, col].get_yticks()[1:7])
      axarr[row, col].text(-2, -3, e_labels[e], weight='bold', size='large')
      col += 1
      i += 1

  # adjust plot settings
  f.subplots_adjust(hspace=0, wspace=0)    #smush plots together
  f.text(0.5, 0.04, '[Fe/H]', ha='center', size='large')
  f.text(0.04, 0.5, '[X/Fe]', va='center', size='large')
  sm = cm.ScalarMappable(norm=colors.Normalize(), cmap='jet')
  sm.set_array(teff)
  cb_ax = f.add_axes([0.9, 0.1, 0.025, 0.8])
  cb = f.colorbar(sm, orientation='vertical', cax=cb_ax)
  cb.set_label('$T_{eff}$ (K)', size='large')
  plt.suptitle(np.str(edge) + ' kpc < R < ' + np.str(edge+2) + ' kpc')
  edge += 2
  plt.show()