Welcome to my lessons

Bo Zhang (NAOC, mailto:bozhang@nao.cas.cn) will have a few lessons on python.

These are very useful knowledge, skills and code styles when you use python to process astronomical data.
All materials can be found on my github page.
jupyter notebook (formerly named ipython notebook) is recommeded to use

These lectures are organized as below:

install python
basic syntax
numerical computing
scientific computing
plotting
astronomical data processing
high performance computing
version control

Astronomical data format standard - .fits files

HDUs
- Primary HDU
- Image HDU
- BinTable HDU (Binary Table)
HDU header
- keywords
- column information for tables
HDU data

reference:

Usually, you will meet these:

1d image stored using WCS (world coordinate system)
- SDSS DR7 spectra
- IRAF spectra fits
2d image stored using WCS (world coordinate system)
- all photometry images (2MASS, WISE, ...)
table / catalog
- LAMOST catalog
- SDSS DR10/12 spectra



In [1]:

    
%pylab inline
from astropy.io import fits

# read hdu list
hl = fits.open('./data/lamost_dr2_spectra/spec-55892-F9205_sp09-174.fits')









    



Populating the interactive namespace from numpy and matplotlib



In [2]:

    
# show fits info
print hl.info()









    



Filename: ./data/lamost_dr2_spectra/spec-55892-F9205_sp09-174.fits
No.    Name         Type      Cards   Dimensions   Format
0    Flux        PrimaryHDU     132   (3903, 5)    float32   
None



In [3]:

    
# read HDU header
hl[0].header









    Out[3]:





SIMPLE  =                    T /Primary Header created by MWRFITS v1.11b        
BITPIX  =                  -32 /                                                
NAXIS   =                    2 / Number of array dimensions                     
NAXIS1  =                 3903 /                                                
NAXIS2  =                    5 /                                                
EXTEND  =                    T /                                                
                                                                                
COMMENT --------FILE INFORMATION                                                
FILENAME= 'spec-55892-F9205_sp09-174.fits' /                                    
OBSID   =                20566 / Unique number ID of this spectrum              
AUTHOR  = 'LAMOST Pipeline'    / Who compiled the information                   
DATA_V  = 'LAMOST DR3'         / Data release version                           
EXTEN0  = 'Flux, Inverse, Wavelength, Andmask, Ormask' /                        
N_EXTEN =                    1 / The extension number                           
EXTNAME = 'Flux    '           / The extension name                             
ORIGIN  = 'NAOC-LAMOST'        / Organization responsible for creating this file
DATE    = '2015-12-16T13:29:15' / Time when this HDU is created (UTC)           
COMMENT --------TELESCOPE PARAMETERS                                            
TELESCOP= 'LAMOST  '           / GuoShouJing Telescope                          
LONGITUD=               117.58 / [deg] Longitude of site                        
LATITUDE=                40.39 / [deg] Latitude of site                         
FOCUS   =                19964 / [mm] Telescope focus                           
CAMPRO  = 'NEWCAM  '           / Camera program name                            
CAMVER  = 'v2.0    '           / Camera program version                         
COMMENT --------OBSERVATION PARAMETERS                                          
DATE-OBS= '2011-11-26T19:40:26.84' / The observation median UTC                 
DATE-BEG= '2011-11-27T03:07:13.0' / The observation start local time            
DATE-END= '2011-11-27T04:55:45.0' / The observation end local time              
LMJD    =                55892 / Local Modified Julian Day                      
MJD     =                55891 / Modified Julian Day                            
LMJMLIST= '80484661-80484699-80484739' / Local Modified Julian Minute list      
PLANID  = 'F9205   '           / Plan ID in use                                 
RA      =           123.396690 / [deg] Right ascension of object                
DEC     =             6.963130 / [deg] Declination of object                    
RA_OBS  =           123.396690 / [deg] Right ascension during observing         
DEC_OBS =             6.963130 / [deg] Declination during observing             
OFFSET  =                    F / Whether there's a offset during observing      
OFFSET_V=                 0.00 / Offset value in arcsecond                      
DESIG   = 'LAMOST J081335.20+065747.2' / Designation of LAMOST target           
FIBERID =                  174 / Fiber ID of Object                             
CELL_ID = 'G2826   '           / Fiber Unit ID on the focal plane               
X_VALUE =             -171.982 / [mm] X coordinate of object on the focal plane 
Y_VALUE =             -273.869 / [mm] Y coordinate of object on the focal plane 
OBJNAME = '257233284699743'    / Name of object                                 
OBJTYPE = 'gal     '           / Object type from input catalog                 
TSOURCE = 'PILOT   '           / Name of input catalog                          
TCOMMENT= 'No      '           / Target information                             
TFROM   = '-       '           / Target catalog                                 
FIBERTYP= 'Obj     '           / Fiber type of object                           
MAGTYPE = 'ugriz   '           / Magnitude type of object                       
MAG1    =                24.62 / [mag] Mag1 of object                           
MAG2    =                16.40 / [mag] Mag2 of object                           
MAG3    =                17.78 / [mag] Mag3 of object                           
MAG4    =                15.26 / [mag] Mag4 of object                           
MAG5    =                22.84 / [mag] Mag5 of object                           
MAG6    =                10.00 / [mag] Mag6 of object                           
MAG7    =                10.00 / [mag] Mag7 of object                           
OBS_TYPE= 'OBJ     '           / The type of target (OBJ, FLAT, ARC or BIAS)    
OBSCOMM = 'Science '           / Science or Test                                
RADECSYS= 'FK5     '           / Equatorial coordinate system                   
EQUINOX =              2000.00 / Equinox in years                               
LAMPLIST= 'lamphgcdne.dat'     / Arc lamp emission line list                    
SKYLIST = 'skylines.dat'       / Sky emission line list                         
EXPID01 = '09b-20111127030708-8-80484661' / ID string for blue exposure 1       
EXPID02 = '09b-20111127034533-9-80484699' / ID string for blue exposure 2       
EXPID03 = '09b-20111127042507-10-80484739' / ID string for blue exposure 3      
EXPID04 = '09r-20111127030742-8-80484661' / ID string for red exposure 1        
EXPID05 = '09r-20111127034606-9-80484699' / ID string for red exposure 2        
EXPID06 = '09r-20111127042540-10-80484739' / ID string for red exposure 3       
NEXP    =                    3 / Number of valid exposures                      
NEXP_B  =                    3 / Number of valid blue exposures                 
NEXP_R  =                    3 / Number of valid red exposures                  
EXPT_B  =              5400.00 / [s] Blue exposure duration time                
EXPT_R  =              5400.00 / [s] Red exposure duration time                 
EXPTIME =              5400.00 / [s] Minimum of exposure time for all cameras   
BESTEXP =             80484661 / MJM of the best exposure                       
SCAMEAN =                 3.10 / [ADU] Mean level of scatter light              
COMMENT --------SPECTROGRAPH PARAMETERS                                         
SPID    =                    9 / Spectrograph ID                                
SPRA    =          123.7406162 / [deg] Average RA of this spectrograph          
SPDEC   =            6.8539852 / [deg] Average DEC of this spectrograph         
SLIT_MOD= 'x2/3    '           / Slit mode, x1, x2/3 or x1/2                    
COMMENT --------WEATHER CONDITION                                               
TEMPCCDB=              -113.10 / [deg] The temperature of blue CCD              
TEMPCCDR=              -122.20 / [deg] The temperature of red CCD               
SEEING  =                 3.00 / [arcsec] Seeing during exposure                
MOONPHA =                 0.01 / [day] Moon phase for a 29.53 days period       
TEMP_AIR=                 1.60 / [deg] Temperature outside dome                 
TEMP_FP =                 6.30 / [degree celsius] Temprature of the focalplane  
DEWPOINT=                -8.60 / [deg]                                          
DUST    = '        '           / Reservation                                    
HUMIDITY=                 0.50 /                                                
WINDD   =                18.70 / [deg] Wind direction                           
WINDS   =                 1.40 / [m/s] Wind speed                               
SKYLEVEL= '        '           / Reservation                                    
COMMENT --------DATA REDUCTION PARAMETERS                                       
EXTRACT = 'aperture'           / Extraction method                              
SFLATTEN=                    T / Super flat has been applied                    
PCASKYSB=                    T / PCA sky-subtraction has been applied           
NSKIES  =                   20 / Sky fiber number                               
SKYCHI2 =                  2.4 / Mean chi^2 of sky-subtraction                  
SCHI2MIN=                  1.6 / Minimum chi^2 of sky-subtraction               
SCHI2MAX=                  3.2 / Maximum chi^2 of sky-subtraction               
NSTD    =                    5 / Number of (good) standard stars                
FSTAR   = '60-76-89-135-152'   / FiberID of flux standard stars                 
FCBY    = 'catalog '           / Standard stars origin (auto, manual or catalog)
HELIO   =                    T / Heliocentric correction                        
HELIO_RV=            -25.38969 / [km/s] Heliocentric correction                 
VACUUM  =                    T / Wavelengths are in vacuum                      
NWORDER =                    2 / Number of linear-log10 coefficients            
WFITTYPE= 'LOG-LINEAR'         / Linear-log10 dispersion                        
COEFF0  =               3.5682 / Central wavelength (log10) of first pixel      
COEFF1  =               0.0001 / Log10 dispersion per pixel                     
WAT0_001= 'system=linear'      /                                                
WAT1_001= 'wtype=linear label=Wavelength units=Angstroms' /                     
CRVAL1  =               3.5682 / Central wavelength (log10) of first pixel      
CD1_1   =               0.0001 / Log10 dispersion per pixel                     
CRPIX1  =                    1 / Starting pixel (1-indexed)                     
CTYPE1  = 'LINEAR  '           /                                                
DC-FLAG =                    1 / Log-linear flag                                
COMMENT --------SPECTRA ANALYSIS RESULTS                                        
VERSPIPE= 'v2.7.5  '           / Version of Pipeline                            
CLASS   = 'STAR    '           / Class of object                                
SUBCLASS= 'F0      '           / Subclass of object                             
Z       =          -0.00001500 / Redshift of object                             
Z_ERR   =           0.00007800 / Redshift error of object                       
ZFLAG   = 'LASP    '           / Which method computes the redshift             
SNRU    =               164.85 / SNR of u filter                                
SNRG    =               446.69 / SNR of g filter                                
SNRR    =               590.62 / SNR of r filter                                
SNRI    =               570.01 / SNR of i filter                                
SNRZ    =               376.18 / SNR of z filter



In [4]:

    
# read HDU data
hl[0].data









    Out[4]:





array([[  1.73385332e+04,   1.78001680e+04,   1.74559727e+04, ...,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
       [  5.20499998e-05,   5.11758699e-05,   5.23977396e-05, ...,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
       [  3.69998633e+03,   3.70083838e+03,   3.70169067e+03, ...,
          9.08238965e+03,   9.08448047e+03,   9.08657324e+03],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00]], dtype=float32)



In [5]:

    
# transform data to WVS coordinates
naxis1 = hl[0].header['NAXIS1']
crval1 = hl[0].header['CRVAL1']
crpix1 = hl[0].header['CRPIX1']
cdelt1 = hl[0].header['CD1_1']

wave = 10**( crval1 + (np.arange(naxis1)+1-crpix1)*cdelt1)
print wave

flux = hl[0].data[0, :]
print flux









    



[ 3699.98531173  3700.83736292  3701.68961033 ...,  9082.3869326
  9084.47847027  9086.57048958]
[ 17338.53320312  17800.16796875  17455.97265625 ...,      0.              0.
      0.        ]



In [6]:

    
fig = plt.figure(figsize=(15,8 ))
plt.plot(wave, flux)









    Out[6]:





[<matplotlib.lines.Line2D at 0x7eff25e74bd0>]



In [7]:

    
from astropy.table import Table, Column
spec_table = Table([Column(wave, 'wave'),
                    Column(flux, 'flux')])
# Table.write('./data/lamost_dr2_spectra/test_spec_table.fits')
spec_table.pprint









    Out[7]:





<bound method Table.pprint of <Table length=3903>
     wave       flux 
   float64    float32
------------- -------
3699.98531173 17338.5
3700.83736292 17800.2
3701.68961033 17456.0
  3702.542054 17002.9
3703.39469398 15677.3
 3704.2475303 14652.4
3705.10056302 14114.7
3705.95379218 13921.9
3706.80721783 14344.3
   3707.66084 14902.9
          ...     ...
9067.75964584     0.0
9069.84781506     0.0
9071.93646515     0.0
9074.02559623     0.0
 9076.1152084     0.0
9078.20530178     0.0
9080.29587648     0.0
 9082.3869326     0.0
9084.47847027     0.0
9086.57048958     0.0>



In [ ]:

how to convert (ra, dec) to (l, b)



In [8]:

    
%pylab
from astropy.coordinates import SkyCoord, Longitude, Latitude
import astropy.units as u









    



Using matplotlib backend: TkAgg
Populating the interactive namespace from numpy and matplotlib



In [9]:

    
ra = Longitude([1, 2, 3], unit=u.deg)  # Could also use Angle
dec = np.array([4.5, 5.2, 6.3]) * u.deg  # Astropy Quantity
c = SkyCoord(ra, dec, frame='icrs')
c









    Out[9]:





<SkyCoord (ICRS): (ra, dec) in deg
    [(1.0, 4.5), (2.0, 5.2), (3.0, 6.3)]>



In [10]:

    
c.info()









    



dtype = object
unit = deg,deg
class = SkyCoord
n_bad = 0
length = 3



In [11]:

    
c.galactic









    Out[11]:





<SkyCoord (Galactic): (l, b) in deg
    [(101.22678436, -56.35949807), (103.32241344, -56.00907284),
     (105.56483231, -55.23635591)]>



In [12]:

    
c.galactic.l, c.galactic.b









    Out[12]:





(<Longitude [ 101.22678436, 103.32241344, 105.56483231] deg>,
 <Latitude [-56.35949807,-56.00907284,-55.23635591] deg>)



In [ ]:

2-D data (images)

recommended packages: aplpy, healpy

temporally these will not be introduced in our course

HOMEWORK

download a LAMOST spectra, and write a function to read data in the fits file.

requirements:
1. usage
```
spec = read_spectrum_lamost('filepath')
```
2. spec should be a 5-column table, including these data:
  
  wavelength, flux, inverse variance, and-mask, or-mask

hint: the data structure of LAMOST DR2 spectra fits files can be found here



In [ ]: