2014年8月6日星期三

【SDSS】SQL search through SDSS website

This is just a note from my work.

SDSS search website:
      http://www.sdss3.org/dr10/
      http://classic.sdss.org/dr7/
  DR10包括DR8和DR9, 但是DR7以及之前的结果都需要在DR7里面搜索。
(1) search particular objects
  http://data.sdss3.org/bulkSpectra/
(2) search objects with certain criterion
  http://skyservice.pha.jhu.edu/casjobs/default.aspx


********************************************************************************
          DR10
********************************************************************************
commands
search at website http://skyserver.sdss3.org/CasJobs/MyDB.aspx or from individual sites.
********************************************************************************
http://skyserver.sdss3.org/dr10/en/help/docs/realquery.aspx
-- a) Finding galaxies by their emission lines:
-- This query selects galaxy spectra with high internal reddening,
-- as measured by the standard Balmer decrement technique. It
-- makes use of the galSpec tables for the measurements of
-- galaxy lines. In this case we use galSpecLine, which has
-- emission line measurements.

SELECT
s.specObjID, s.plate, s.fiberid, s.mjd, s.z, s.zwarning,
g.h_beta_flux, g.h_beta_flux_err,
g.h_alpha_flux, g.h_alpha_flux_err,
g.oiii_5007_flux, g.oiii_5007_flux_err,
g.oiii_5007_eqw
FROM GalSpecLine AS g
  JOIN SpecObj AS s ON s.specobjid = g.specobjid
WHERE
  --h_alpha_flux > h_alpha_flux_err*5
  --AND h_beta_flux > h_beta_flux_err*5
  --AND h_beta_flux_err > 0
  --AND h_alpha_flux > 10.*h_beta_flux
  oiii_5007_flux > oiii_5007_flux_err*5
  AND oiii_5007_eqw > 600.
  AND s.class = 'GALAXY'
  AND s.zwarning = 0

********************************************************************************
          DR7
********************************************************************************
http://cas.sdss.org/dr7/en/tools/search/sql.asp

SELECT
G.objId, S.plate, S.fiberID, S.mjd, S.z,  L1.ew,L1.ewerr, L2.ew,L2.ewerr,L1.ew*L1.continuum,L1.ewerr*L1.continuum

FROM Galaxy as G
   JOIN SpecObj as S ON G.objId=S.bestObjId
   JOIN SpecLine as hL1 ON S.specObjId=L1.specObjId
   JOIN SpecLine as L2 ON S.specObjId=L2.specObjId  
   JOIN SpecLine as Lb ON S.specObjId=Lb.specObjId
WHERE
     L1.LineId = 5008
     and L1.ew > 1000.
     and L1.ew/L1.ewerr > 3.
     and L2.LineId =6565
     and Lb.LineId =4863
     --and log10( (L1.ew*L1.continuum)/(Lb.ew*Lb.continuum) ) < 0.75


********************************************************************************
********************************************************************************
Check image:
   http://cas.sdss.org/dr7/en/tools/chart/list.asp
   http://skyserver.sdss3.org/public/en/tools/chart/listinfo.aspx

********************************************************************************
# Example script for looking at BOSS spectra and redshift fits via Python.
#
# Written by Adam S. Bolton, University of Utah, Oct. 2009
#

# Imports:
import numpy as n
import pyfits as pf
import matplotlib as mpl
mpl.use('TkAgg')
mpl.interactive(True)
from matplotlib import pyplot as p

# Set topdir:
topdir = '/data/BOSS/spectro/redux/v5_5_12/'

# Pick your plate/mjd and read the data:
plate = '3621'
mjd = '55104'
spfile = topdir + plate + '/spPlate-' + plate + '-' + mjd + '.fits'
zbfile = topdir + plate + '/spZbest-' + plate + '-' + mjd + '.fits'
hdulist = pf.open(spfile)
c0 = hdulist[0].header['coeff0']
c1 = hdulist[0].header['coeff1']
npix = hdulist[0].header['naxis1']
wave = 10.**(c0 + c1 * n.arange(npix))
# Following commented-out bit was needed for some of the early redux:
#bzero = hdulist[0].header['bzero']
bunit = hdulist[0].header['bunit']
flux = hdulist[0].data
ivar = hdulist[1].data
hdulist.close()
hdulist = 0
hdulist = pf.open(zbfile)
synflux = hdulist[2].data
zstruc = hdulist[1].data
hdulist.close()
hdulist = 0

i = 499
# Set starting fiber point (above), then copy and paste
# the following repeatedly to loop over spectra:
i+=1
# Following commented-out bit was needed for some of the early redux:
#p.plot(wave, (flux[i,:]-bzero) * (ivar[i,:] > 0), 'k', hold=False)
p.plot(wave, flux[i,:] * (ivar[i,:] > 0), 'k', hold=False)
p.plot(wave, synflux[i,:], 'g', hold=True)
p.xlabel('Angstroms')
p.ylabel(bunit)
p.title(zstruc[i].field('class') + ', z = ' + str(zstruc[i].field('z')))




********************************************************************************
  SDSS get the images and the spectra
********************************************************************************
One example from Peter Erwin:   http://www.mpe.mpg.de/~erwin/code/telarchive/
# This is the URL of the SDSS images ( from astroML)
http://www.astroml.org/sklearn_tutorial/auto_examples/plot_sdss_images.html
url = ("http://casjobs.sdss.org/ImgCutoutDR7/"
           "getjpeg.aspx?ra=%.8f&dec=%.8f&scale=%.2f&width=%i&height=%i"
           % (RA, DEC, scale, width, height))
However, this link does not work now. Please use :
SDSS_URL = ('http://skyservice.pha.jhu.edu/DR7/ImgCutout/getjpeg.aspx?'
            'ra=%(RA).8f&dec=%(DEC).8f&scale=%(scale).2f&width=%(width)i&height=%(height)i&opt=&query=')

# This is the URL of the sdss fits spectra (from astroML)
FITS_FILENAME = 'spSpec-%(mjd)05i-%(plate)04i-%(fiber)03i.fit'
SDSS_URL = ('http://das.sdss.org/spectro/1d_26/%(plate)04i/1d/spSpec-%(mjd)05i-%(plate)04i-%(fiber)03i.fit')



********************************************************************************
  DR7   about the spectra files 
********************************************************************************
http://classic.sdss.org/dr7/products/spectra/read_spSpec.html
The spSpec file's header contains several parameters such as object identification and coordinates, observing information, spectral classification, redshift, etc..
The fits image contains a 4x(roughly 4000) image whose
  first row is the spectrum,
  second the continuum subtracted spectrum,
  third the error,
  forth a bitmask.
  fifth ?   sky light ?
Notice that the wavelength vector is not contained in the image, but must be generated from parameters in the header.
SDSS spectra are binned in constant Log(Λ ) and the wavelength can be obtained from the header parameters COEFF0 and COEFF1 (or alternatively CRVAL1 and CD1_1) as follows:

lambda = 10**(COEFF0 + COEFF1*i), where i denotes the (zero indexed) pixel number.

Six binary tables are included, most importantly,
 HDU 2 which contains a table of gaussian fits to line positions,
 HDU 3 which contains a table on line index parameters. Also included, but of less general interest, are:
 HDU 1 which contains line parameters used in the emission-line redshift determination;
 HDU 3 which contains information on emission-line redshifts;
 HDU 4 which contains a table of cross-correlation redshifts with all the templates, and
 HDU 6 which has additional mask information as well as the spectral resolution as a function of wavelength.

2014年8月4日星期一

【PYTHON】Python中filter、map、reduce、lambda 的用法

filter(function, sequence):对sequence中的item依次执行function(item),将执行结果为True的item组成一个List/String/Tuple(取决于sequence的类型)返回:
>>> def f(x): return x % 2 != 0 and x % 3 != 0 
>>> filter(f, range(2, 25)) 
[5, 7, 11, 13, 17, 19, 23]
>>> def f(x): return x != 'a' 
>>> filter(f, "abcdef") 
'bcdef'

map(function, sequence) :对sequence中的item依次执行function(item),见执行结果组成一个List返回:
>>> def cube(x): return x*x*x 
>>> map(cube, range(1, 11)) 
[1, 8, 27, 64, 125, 216, 343, 512, 729, 1000]
>>> def cube(x) : return x + x 
... 
>>> map(cube , "abcde") 
['aa', 'bb', 'cc', 'dd', 'ee']
另外map也支持多个sequence,这就要求function也支持相应数量的参数输入:
>>> def add(x, y): return x+y 
>>> map(add, range(8), range(8)) 
[0, 2, 4, 6, 8, 10, 12, 14]

reduce(function, sequence, starting_value):对sequence中的item顺序迭代调用function,如果有starting_value,还可以作为初始值调用,例如可以用来对List求和:
>>> def add(x,y): return x + y 
>>> reduce(add, range(1, 11)) 
55 (注:1+2+3+4+5+6+7+8+9+10)
>>> reduce(add, range(1, 11), 20) 
75 (注:1+2+3+4+5+6+7+8+9+10+20)

lambda:这是Python支持一种有趣的语法,它允许你快速定义单行的最小函数,类似与C语言中的宏,这些叫做lambda的函数,是从LISP借用来的,可以用在任何需要函数的地方: 
>>> g = lambda x: x * 2 
>>> g(3) 

>>> (lambda x: x * 2)(3) 
6


我们也可以把filter map reduce 和lambda结合起来用,函数就可以简单的写成一行。
例如
kmpathes = filter(lambda kmpath: kmpath,                  
map(lambda kmpath: string.strip(kmpath),
string.split(l, ':')))              
看起来麻烦,其实就像用语言来描述问题一样,非常优雅。
对 l 中的所有元素以':'做分割,得出一个列表。对这个列表的每一个元素做字符串strip,形成一个列表。对这个列表的每一个元素做直接返回操作(这个地方可以加上过滤条件限制),最终获得一个字符串被':'分割的列表,列表中的每一个字符串都做了strip,并可以对特殊字符串过滤。



---------------------------------------------------------------

lambda表达式返回一个函数对象
例子:
func = lambda x,y:x+y
func相当于下面这个函数
def func(x,y):
    return x+y
 
注意def是语句而lambda是表达式
下面这种情况下就只能用lambda而不能用def
[(lambda x:x*x)(x) for x in range(1,11)]
 
map,reduce,filter中的function都可以用lambda表达式来生成!
 
map(function,sequence)
把sequence中的值当参数逐个传给function,返回一个包含函数执行结果的list。
如果function有两个参数,即map(function,sequence1,sequence2)。
 
例子:
求1*1,2*2,3*3,4*4
map(lambda x:x*x,range(1,5))
返回值是[1,4,9,16]
 
reduce(function,sequence)
function接收的参数个数只能为2
先把sequence中第一个值和第二个值当参数传给function,再把function的返回值和第三个值当参数传给
function,然后只返回一个结果。
 
例子:
求1到10的累加
reduce(lambda x,y:x+y,range(1,11))
返回值是55。
 
filter(function,sequence)
function的返回值只能是True或False
把sequence中的值逐个当参数传给function,如果function(x)的返回值是True,就把x加到filter的返回值里面。一般来说filter的返回值是list,特殊情况如sequence是string或tuple,则返回值按照sequence的类型。
 
例子:
找出1到10之间的奇数
filter(lambda x:x%2!=0,range(1,11))
返回值
[1,3,5,7,9]
 
如果sequence是一个string
filter(lambda x:len(x)!=0,'hello')返回'hello'
filter(lambda x:len(x)==0,'hello')返回''