Python tutorial 2

1. random number generators: numpy.random

https://docs.scipy.org/doc/numpy/reference/routines.random.html


In [1]:
import numpy as np
print(dir(np.random))


['Lock', 'RandomState', '__RandomState_ctor', '__all__', '__builtins__', '__cached__', '__doc__', '__file__', '__loader__', '__name__', '__package__', '__path__', '__spec__', '_numpy_tester', 'absolute_import', 'bench', 'beta', 'binomial', 'bytes', 'chisquare', 'choice', 'dirichlet', 'division', 'exponential', 'f', 'gamma', 'geometric', 'get_state', 'gumbel', 'hypergeometric', 'info', 'laplace', 'logistic', 'lognormal', 'logseries', 'mtrand', 'multinomial', 'multivariate_normal', 'negative_binomial', 'noncentral_chisquare', 'noncentral_f', 'normal', 'np', 'operator', 'pareto', 'permutation', 'poisson', 'power', 'print_function', 'rand', 'randint', 'randn', 'random', 'random_integers', 'random_sample', 'ranf', 'rayleigh', 'sample', 'seed', 'set_state', 'shuffle', 'standard_cauchy', 'standard_exponential', 'standard_gamma', 'standard_normal', 'standard_t', 'test', 'triangular', 'uniform', 'vonmises', 'wald', 'warnings', 'weibull', 'zipf']

how to draw samples from a gaussian distribution


In [2]:
%pylab inline
import matplotlib.pyplot as plt

from matplotlib import rcParams
rcParams.update({'font.size': 20})

rdata = np.random.randn(1000)
fig = plt.figure(figsize=(6, 4))
plt.hist(rdata)
print(np.mean(rdata), np.median(rdata), np.std(rdata))


Populating the interactive namespace from numpy and matplotlib
0.0131778852121 -0.0127231840331 1.00130532414

In [3]:
np.std(np.random.randn(1000) + np.random.randn(1000))


Out[3]:
1.4343155445498037

other distributions ...


In [4]:
randexp = np.random.exponential(2., size=(1000))
hist(randexp, np.linspace(0,10,50));



In [5]:
randps = np.random.poisson(10, size=(10000,))
hist(randps, np.arange(20));



In [6]:
M = 4.
m = 15.
merr = 0.1

rand_m = np.random.randn(1000)*0.1+m
hist(rand_m);


$\log_{10}(d) = 1 + \mu /5 $


In [7]:
rand_d = 10.**(1+0.2*(rand_m-M))
hist(rand_d, np.linspace(1300, 1900, 30));



In [ ]:

2. plotting


In [8]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams.update({'font.size':20})

In [9]:
# fig = plt.figure(figsize=(10,10))
x = np.linspace(0, 6*np.pi, 100)
plt.plot(x, np.cos(x), 'rv--');
plt.plot(x, np.sin(x), 'bs-.', alpha=.1);



In [10]:
plt.scatter(x, np.cos(x)+0.2, s=np.random.rand(*x.shape)*80, c=np.sin(x)+1)


Out[10]:
<matplotlib.collections.PathCollection at 0x7f58a6191940>

3. IO (text files & fits files)


In [11]:
# use numpy.savetxt & numpy.loadtxt
a = np.random.randn(4, 5)
print(a)


[[-0.97016615 -1.16455194  1.12048518  0.50861574 -0.10792383]
 [ 2.86641519 -0.65292932 -0.14310716  0.01853606 -0.73349613]
 [-0.14018446  1.31599772  0.36529062  0.02056543  0.49264984]
 [ 0.61179898  0.01365724  0.38080106 -0.03712658 -1.17633178]]

In [12]:
np.savetxt('./data/text/rdata.dat', a)

In [13]:
!gedit ./data/text/rdata.dat

In [14]:
b = np.loadtxt('./data/text/rdata.dat')
print(b)


[[-0.97016615 -1.16455194  1.12048518  0.50861574 -0.10792383]
 [ 2.86641519 -0.65292932 -0.14310716  0.01853606 -0.73349613]
 [-0.14018446  1.31599772  0.36529062  0.02056543  0.49264984]
 [ 0.61179898  0.01365724  0.38080106 -0.03712658 -1.17633178]]

In [15]:
a==b.reshape(4, 5)


Out[15]:
array([[ True,  True,  True,  True,  True],
       [ True,  True,  True,  True,  True],
       [ True,  True,  True,  True,  True],
       [ True,  True,  True,  True,  True]], dtype=bool)

In [16]:
impath = "./data/image_data/G178_final.850.fits"

%pylab inline
from matplotlib import rcParams
rcParams.update({'font.size': 20})

from aplpy import FITSFigure
fig = FITSFigure(impath)
fig.show_colorscale()


Populating the interactive namespace from numpy and matplotlib
INFO: Setting slices=[0] [aplpy.core]
INFO: Auto-setting vmin to -9.026e-02 [aplpy.core]
INFO: Auto-setting vmax to  1.985e-01 [aplpy.core]

In [ ]:


In [17]:
impath = "./data/wise_image/w1_cut.fits"

%pylab inline
%matplotlib inline
from matplotlib import rcParams
rcParams.update({'font.size': 20})
from aplpy import FITSFigure
fig = FITSFigure(impath)
fig.show_colorscale()
fig.show_grayscale()


Populating the interactive namespace from numpy and matplotlib
INFO: Auto-setting vmin to -3.489e+01 [aplpy.core]
INFO: Auto-setting vmax to  4.563e+02 [aplpy.core]
INFO: Auto-setting vmin to -3.489e+01 [aplpy.core]
INFO: Auto-setting vmax to  4.563e+02 [aplpy.core]

In [ ]:

LAMOST spectra


In [18]:
ls ./data/lamost_dr2_spectra/


LAMOST DR2 - LAMOST Data Release Two.pdf*  spec-55911-B91105_sp09-181.fits*
spec-55892-F9205_sp09-174.fits*            spec-55915-F5591504_sp05-140.fits*
spec-55899-B9903_sp07-046.fits*

In [19]:
specpath = "./data/lamost_dr2_spectra/spec-55892-F9205_sp09-174.fits"

In [20]:
from astropy.io import fits

In [21]:
hl = fits.open(specpath)

In [22]:
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   

In [23]:
hl


Out[23]:
[<astropy.io.fits.hdu.image.PrimaryHDU object at 0x7f58a31abf98>]

In [24]:
hl[0]


Out[24]:
<astropy.io.fits.hdu.image.PrimaryHDU at 0x7f58a31abf98>

In [25]:
hl[0].header


Out[25]:
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                                

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 ' /


In [26]:
wave = 10.**(hl[0].header['CRVAL1']+np.arange(hl[0].header['NAXIS1'])*hl[0].header['CD1_1'])

In [27]:
wave


Out[27]:
array([ 3699.98531173,  3700.83736292,  3701.68961033, ...,  9082.3869326 ,
        9084.47847027,  9086.57048958])

In [28]:
np.log10(wave)


Out[28]:
array([ 3.5682,  3.5683,  3.5684, ...,  3.9582,  3.9583,  3.9584])

In [29]:
flux = hl[0].data # [flux, ivar, wave, and_mask, or_mask]

In [30]:
flux


Out[30]:
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 [31]:
%pylab
%matplotlib inline
fig = figure(figsize=(10, 5))
plt.plot(wave, flux[0, :])


Using matplotlib backend: TkAgg
Populating the interactive namespace from numpy and matplotlib
Out[31]:
[<matplotlib.lines.Line2D at 0x7f58cab8e2e8>]

In [32]:
# fig.savefig("here goes the file path")

In [ ]:


In [ ]: