Python and Astronomy

In [2]:

    

import this


    

In [3]:

    

import antigravity


    

In [4]:

    

from __future__ import division, print_function, absolute_import, unicode_literals


    

In [5]:

    

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import rcParams


    

In [6]:

    

%matplotlib inline


    

In [7]:

    

rcParams['figure.figsize'] = 13, 8.1


    

In [8]:

    

import astropy.units as u
import astropy.constants as c


    

In [10]:

    

658.28 * u.nm + 10 * u.Angstrom


    

    
Out[10]:





$659.28 \; \mathrm{nm}$

In [11]:

    

658.29 * u.nm + 10


    

    




---------------------------------------------------------------------------
UnitsError                                Traceback (most recent call last)
<ipython-input-11-bbe3a116b0ed> in <module>()
----> 1 658.29 * u.nm + 10

/mnt/data/ajh/Documents/virtualenvs/python3.4.3/lib/python3.4/site-packages/astropy/units/quantity.py in __array_prepare__(self, obj, context)
    344                                      "argument is not a quantity (unless the "
    345                                      "latter is all zero/infinity/nan)"
--> 346                                      .format(function.__name__))
    347             except TypeError:
    348                 # _can_have_arbitrary_unit failed: arg could not be compared

UnitsError: Can only apply 'add' function to dimensionless quantities when other argument is not a quantity (unless the latter is all zero/infinity/nan)

In [12]:

    

sky_brightness_gr = (22.5, 21.5) * u.ABmag
pivot_wavelength_gr = (4702, 6175) * u.Angstrom
band_width_gr = (138.7, 137.3) * u.nm
pixel_scale = 1.0 * u.arcsecond
sky_solid_angle = pixel_scale**2 / u.pixel
aperture_area = np.pi * (3.9**2 - 1.5**2) / 4 * u.m**2
throughput = 0.7
QE = 0.9 * u.electron / u.photon


    

In [13]:

    

sky_sfd_gr = sky_brightness_gr.to(u.Watt / (u.m**2 * u.micron), \
                                equivalencies=u.equivalencies.spectral_density(pivot_wavelength_gr))
sky_sfd_gr


    

    
Out[13]:





$[4.9232904 \times 10^{-17},~7.1704561 \times 10^{-17}] \; \mathrm{\frac{W}{\mu m\,m^{2}}}$

In [14]:

    

focal_plane_flux = sky_sfd_gr * band_width_gr * aperture_area * throughput * sky_solid_angle / u.arcsecond**2
focal_plane_flux.to(u.Watt / u.pixel)


    

    
Out[14]:





$[4.8654704 \times 10^{-17},~7.0147184 \times 10^{-17}] \; \mathrm{\frac{W}{pix}}$

In [15]:

    

photon_energy_gr = c.h * c.c / (pivot_wavelength_gr * u.photon)
photon_energy_gr.to(u.J / u.photon)


    

    
Out[15]:





$[4.2246824 \times 10^{-19},~3.2169161 \times 10^{-19}] \; \mathrm{\frac{J}{ph}}$

In [16]:

    

electron_rate = focal_plane_flux / photon_energy_gr * QE
electron_rate.to(u.electron / (u.second * u.pixel))


    

    
Out[16]:





$[103.65095,~196.25151] \; \mathrm{\frac{e^{-}}{pix\,s}}$

In [17]:

    

from astropy.table import Table


    

In [18]:

    

with open('resources/nonsense.csv') as f:
    print(f.read())


    

    




ID,RA,dec,g_mag
1,04 23 15.3,-22 36 19,26.5
2,22 17 34.7,+02 22 59,23.4
3,10 44 29.1,-48 29 51,24.0
4,09 14 23.1,-10 10 41,27.3
5,15 51 11.6,+13 43 17,25.1

In [19]:

    

nonsense = Table.read('resources/nonsense.csv')
nonsense


    

    
Out[19]:




<Table length=5>

ID
RA
dec
g_mag
int64
str10
str9
float64
1
04 23 15.3
-22 36 19
26.5
2
22 17 34.7
+02 22 59
23.4
3
10 44 29.1
-48 29 51
24.0
4
09 14 23.1
-10 10 41
27.3
5
15 51 11.6
+13 43 17
25.1

In [20]:

    

nonsense.write('resources/nonsense.tex', format='ascii.latex')


    

In [21]:

    

with open('resources/nonsense.tex') as f:
    print(f.read())


    

    




\begin{table}
\begin{tabular}{cccc}
ID & RA & dec & g_mag \\
1 & 04 23 15.3 & -22 36 19 & 26.5 \\
2 & 22 17 34.7 & +02 22 59 & 23.4 \\
3 & 10 44 29.1 & -48 29 51 & 24.0 \\
4 & 09 14 23.1 & -10 10 41 & 27.3 \\
5 & 15 51 11.6 & +13 43 17 & 25.1 \\
\end{tabular}
\end{table}

In [22]:

    

import os, glob, subprocess


    

In [23]:

    

import ccdproc
from astroML.plotting import hist
from matplotlib.colors import LogNorm
from astropy.wcs import WCS


    

In [24]:

    

raw_path = os.path.abspath('resources/raw/')
calib_path = os.path.abspath('resources/calibrations/')
temp_path = os.path.abspath('resources/temp/')
red_path = os.path.abspath('resources/reduced/')


    

In [25]:

    

os.makedirs(temp_path, mode=0o775, exist_ok=True)
os.makedirs(red_path, mode=0o775, exist_ok=True)


    

In [26]:

    

os.chdir(temp_path)


    

In [27]:

    

os.listdir(calib_path)


    

    
Out[27]:





['master_unbiased_dark.fits', 'master_flat.fits', 'master_bias.fits']

In [28]:

    

bias = ccdproc.CCDData.read(os.path.join(calib_path, 'master_bias.fits'))
dark = ccdproc.CCDData.read(os.path.join(calib_path, 'master_unbiased_dark.fits'))
flat = ccdproc.CCDData.read(os.path.join(calib_path, 'master_flat.fits'))


    

In [29]:

    

def isbadpixel(data, min_value, max_value):
    isgoodpixel = np.logical_and(np.isfinite(data), \
                                np.logical_and(data >= min_value, data <= max_value))
    return np.logical_not(isgoodpixel)


    

In [30]:

    

bias.mask = np.where(isbadpixel(bias.data, 1,60000), True, False)
bias.data = np.where(isbadpixel(bias.data, 1,60000), 0.0, bias.data)
dark.mask = np.where(isbadpixel(dark.data, -1000, 60000), True, False)
dark.data = np.where(isbadpixel(dark.data, -1000, 60000), 0.0, dark.data)
flat.mask = np.where(isbadpixel(flat.data, 0.7, 1.4), True, False)
flat.data = np.where(isbadpixel(flat.data, 0.7, 1.4), 0.0, flat.data)


    

    




/mnt/data/ajh/Documents/virtualenvs/python3.4.3/lib/python3.4/site-packages/ipykernel/__main__.py:2: RuntimeWarning: invalid value encountered in greater_equal
  from ipykernel import kernelapp as app
/mnt/data/ajh/Documents/virtualenvs/python3.4.3/lib/python3.4/site-packages/ipykernel/__main__.py:2: RuntimeWarning: invalid value encountered in less_equal
  from ipykernel import kernelapp as app

In [31]:

    

plt.imshow(bias, interpolation='none', aspect='equal', cmap='cubehelix', vmin=1025, vmax=1050)
plt.colorbar()


    

    
Out[31]:





<matplotlib.colorbar.Colorbar at 0x7f64a3b8db70>






    

In [32]:

    

plt.imshow(dark, interpolation='none', aspect='equal', cmap='cubehelix', vmax=40, vmin=0)
plt.colorbar()


    

    
Out[32]:





<matplotlib.colorbar.Colorbar at 0x7f64a02ab0f0>






    

In [33]:

    

plt.imshow(flat, interpolation='none', aspect='equal', cmap='cubehelix', vmax=1.4, vmin=0.7)
plt.colorbar()


    

    
Out[33]:





<matplotlib.colorbar.Colorbar at 0x7f64a01e87f0>






    

In [34]:

    

os.listdir(raw_path)


    

    
Out[34]:





['83F011167_127_light.fits',
 '83F011167_126_light.fits',
 '83F011167_129_light.fits']

In [35]:

    

raw_filenames = glob.glob(os.path.join(raw_path, '*_light.fits'))


    

In [36]:

    

images = [ccdproc.CCDData.read(filename, unit='adu') for filename in raw_filenames]


    

    




WARNING: FITSFixedWarning: EPOCH = 'JNOW ' / epoch of coordinates 
a floating-point value was expected. [astropy.wcs.wcs]

In [37]:

    

images[0].header


    

    
Out[37]:





SIMPLE  =                    T / file does conform to FITS standard             
BITPIX  =                   16 / number of bits per data pixel                  
NAXIS   =                    2 / number of data axes                            
NAXIS1  =                 3352 / length of data axis 1                          
NAXIS2  =                 2532 / length of data axis 2                          
EXTEND  =                    T / FITS dataset may contain extensions            
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
BZERO   =                32768 / offset data range to that of unsigned short    
BSCALE  =                    1 / default scaling factor                         
EXPTIME =                 300. / exposure time (seconds)                        
TEMPERAT=                 -20. / temperature (C)                                
IMAGETYP= 'light   '           / image type                                     
FILTNUM =                    0                                                  
DATE    = '2014-09-23T18:55:31' / file creation date (YYYY-MM-DDThh:mm:ss UT)   
SERIALNO= '83F011167'          / serial number                                  
TARGET  = 'M42     '           / target name                                    
RA      = 'Succeeded.'         / right ascension                                
DEC     = 'RA:     '           / declination                                    
EPOCH   = 'JNOW    '           / epoch of coordinates                           
OBJCTRA = 'Succeeded.'         / right ascension                                
OBJCTDEC= 'RA:     '           / declination                                    
ALTITUDE=                   5. / Altitude (deg)                                 
AZIMUTH =                  36. / Azimuth (deg)                                  

In [38]:

    

[print(image.header['EXPTIME']) for image in images]


    

    




300.0
10.0
60.0






    
Out[38]:





[None, None, None]

In [39]:

    

plt.imshow(images[0], norm=LogNorm(vmin=1000,vmax=64000), interpolation='none', aspect='equal', cmap='cubehelix')
plt.colorbar()


    

    
Out[39]:





<matplotlib.colorbar.Colorbar at 0x7f64a0142588>






    

In [40]:

    

plt.imshow(images[1], norm=LogNorm(vmin=250,vmax=64000), interpolation='none', aspect='equal', cmap='cubehelix')
plt.colorbar()


    

    
Out[40]:





<matplotlib.colorbar.Colorbar at 0x7f649ff590b8>






    

In [41]:

    

plt.imshow(images[2], norm=LogNorm(vmin=10,vmax=64000), interpolation='none', aspect='equal', cmap='cubehelix')
plt.colorbar()


    

    
Out[41]:





<matplotlib.colorbar.Colorbar at 0x7f649fee5080>






    

In [42]:

    

images = [ccdproc.trim_image(image[30:,:]) for image in images]


    

    




WARNING: VerifyWarning: Keyword name 'trim_image' is greater than 8 characters or contains characters not allowed by the FITS standard; a HIERARCH card will be created. [astropy.io.fits.card]

In [43]:

    

for image in images:
    image.mask = np.where(isbadpixel(image.data, 0, 64000), True, False)


    

In [44]:

    

images = [ccdproc.subtract_bias(image, bias) for image in images]


    

    




WARNING: VerifyWarning: Keyword name 'subtract_bias' is greater than 8 characters or contains characters not allowed by the FITS standard; a HIERARCH card will be created. [astropy.io.fits.card]

In [45]:

    

images = [ccdproc.subtract_dark(image, dark, exposure_time='EXPTIME', \
                                exposure_unit=u.second, scale=True) for image in images]


    

    




WARNING: VerifyWarning: Keyword name 'subtract_dark' is greater than 8 characters or contains characters not allowed by the FITS standard; a HIERARCH card will be created. [astropy.io.fits.card]

In [46]:

    

images = [ccdproc.flat_correct(image, flat) for image in images]


    

    




/mnt/data/ajh/Documents/virtualenvs/python3.4.3/lib/python3.4/site-packages/astropy/nddata/mixins/ndarithmetic.py:99: RuntimeWarning: divide by zero encountered in true_divide
  data = operation(self.data * self_unit, operand.data * operand_unit)
WARNING: VerifyWarning: Keyword name 'flat_correct' is greater than 8 characters or contains characters not allowed by the FITS standard; a HIERARCH card will be created. [astropy.io.fits.card]

In [47]:

    

images[0].write('m42_300.fits', clobber=True)
images[1].write('m42_060.fits', clobber=True)
images[2].write('m42_010.fits', clobber=True)


    

In [46]:

    

subprocess.call(['solve-field', '--use-sextractor', '--no-plot', '--overwrite', \
                '--ra', '83.81', '--dec', '-5.389', '--radius', '1', \
                'm42_010.fits', 'm42_060.fits', 'm42_300.fits'])


    

    
Out[46]:





0

In [48]:

    

images[0].wcs = WCS('m42_300.wcs')
images[1].wcs = WCS('m42_060.wcs')
images[2].wcs = WCS('m42_010.wcs')


    

    




WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]

In [49]:

    

images[1] = ccdproc.wcs_project(images[1], images[0].wcs)
images[2] = ccdproc.wcs_project(images[2], images[0].wcs)


    

    




/mnt/data/ajh/Documents/virtualenvs/python3.4.3/lib/python3.4/site-packages/ccdproc/core.py:829: RuntimeWarning: invalid value encountered in greater
  reprojected_mask = reprojected_mask > 1e-8
WARNING: VerifyWarning: Keyword name 'wcs_project' is greater than 8 characters or contains characters not allowed by the FITS standard; a HIERARCH card will be created. [astropy.io.fits.card]

In [50]:

    

images = [image.divide(image.header['EXPTIME']) for image in images]


    

In [51]:

    

combiner = ccdproc.Combiner(images)
weights = np.zeros((3, bias.data.shape[0], bias.data.shape[1]))
weights[0] = 1.0
weights[1] = 60/300.0
weights[2] = 10/300.0
combiner.weights = weights
m42 = combiner.average_combine()
m42.wcs = images[0].wcs


    

In [52]:

    

plt.subplot(111, projection=m42.wcs)
plt.imshow(m42, norm=LogNorm(vmin=3,vmax=3000), cmap='cubehelix')
plt.colorbar()


    

    
Out[52]:





<matplotlib.colorbar.Colorbar at 0x7f64a017fc18>






    

In [53]:

    

plt.subplot(111, projection=m42.wcs)
plt.imshow(m42, norm=LogNorm(vmin=3,vmax=3000), cmap='cubehelix')
plt.xlim(2400,1100)
plt.ylim(610,1640)
plt.colorbar()


    

    
Out[53]:





<matplotlib.colorbar.Colorbar at 0x7f64998140b8>