Week 14 - Astropy

Today's Agenda

  • Useful functions of Astropy
  • Units
  • Time
  • Coordinates
  • FITS files
  • Analytic functions
  • AstroPy Tables and different formats

Astropy is a package that is meant to provide a lot of basic functionality for astronomy work in Python

This can be roughly broken up into two areas. One is astronomical calculations:

  • unit and physical quantity conversions
  • physical constants specific to astronomy
  • celestial coordinate and time transformations

The other is file type and structures:

  • FITS files, implementing the former standalone PyFITS interface
  • Virtual Observatory (VO) tables
  • common ASCII table formats, e.g. for online catalogues or data supplements of scientific publications
  • Hierarchical Data Format (HDF5) files

AstroPy normallly comes with the Anaconda installation. But in case you happen to not have it installed it on your computer, you can simply do a

pip install --no-deps astropy

You can always update it via

conda update astropy

This is just a glimpse of all the features that AstroPy has:

For purposes of today, we'll focus just on what astropy can do for units, time, coordinates, image manipulation, and more.


In [85]:
# Importing Modules
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
sns.set_context("notebook")

In [86]:
import astropy
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy import constants as const

Units

Astropy.units introduces units and allows for unit conversions. It doesn't, however, correctly handle spherical coordinates, but the astropy.coordinates package will address this later.

These units can be used to create objects that are made up of both a value and a unit, and basic math can be easily carried out with these. We can add the .unit and .value properties to get the units and numerical values, respectively.


In [87]:
d=42*u.meter
t=6*u.second
v=d/t
print v
print v.unit


7.0 m / s
m / s

Astropy includes a large number of units, and this can include imperial units as well if desired by importing and enabling imperial units. The .find_equivalent_units() function will also return all the other units that are already defined in astropy. Below we do a quick list of the units that are defined for time and length units


In [88]:
from astropy.units import imperial
imperial.enable()
print( u.s.find_equivalent_units() )
print( u.m.find_equivalent_units() )


  Primary name | Unit definition | Aliases 
[
  a            | 3.15576e+07 s   | annum    ,
  d            | 86400 s         | day      ,
  fortnight    | 1.2096e+06 s    |          ,
  h            | 3600 s          | hour, hr ,
  min          | 60 s            | minute   ,
  s            | irreducible     | second   ,
  sday         | 86164.1 s       |          ,
  wk           | 604800 s        | week     ,
  yr           | 3.15576e+07 s   | year     ,
]
  Primary name | Unit definition | Aliases                         
[
  AU           | 1.49598e+11 m   | au, astronomical_unit            ,
  Angstrom     | 1e-10 m         | AA, angstrom                     ,
  cm           | 0.01 m          | centimeter                       ,
  earthRad     | 6.37814e+06 m   | R_earth, Rearth                  ,
  ft           | 0.3048 m        | foot                             ,
  fur          | 201.168 m       | furlong                          ,
  inch         | 0.0254 m        |                                  ,
  jupiterRad   | 7.1492e+07 m    | R_jup, Rjup, R_jupiter, Rjupiter ,
  lyr          | 9.46073e+15 m   | lightyear                        ,
  m            | irreducible     | meter                            ,
  mi           | 1609.34 m       | mile                             ,
  micron       | 1e-06 m         |                                  ,
  mil          | 2.54e-05 m      | thou                             ,
  nmi          | 1852 m          | nauticalmile, NM                 ,
  pc           | 3.08568e+16 m   | parsec                           ,
  solRad       | 6.95508e+08 m   | R_sun, Rsun                      ,
  yd           | 0.9144 m        | yard                             ,
]

The package also provides constants, with the units included. The full list of units can be found here. We can take a quick look at c and G below, and see that these are objects which have value, uncertainty, and units.


In [89]:
print const.c
print const.G


  Name   = Speed of light in vacuum
  Value  = 299792458.0
  Uncertainty  = 0.0
  Unit  = m / s
  Reference = CODATA 2010
  Name   = Gravitational constant
  Value  = 6.67384e-11
  Uncertainty  = 8e-15
  Unit  = m3 / (kg s2)
  Reference = CODATA 2010

Astropy has an aditional function that will allow for unit conversions. So we can, for example, create an object that is the distance to Mars, and then convert that to kilometers or miles. A brief note is that if you try to convert a pure unit (like the 4th line below) into another unit, you'll get a unitless value representing the conversion between the two.

This can also be used to convert constants into other units, so we can convert the speed of light to the somewhat useful pc/yr or the entirely unuseful furlong/fortnight


In [90]:
Mars=1.5*u.AU
print Mars.to('kilometer')
print Mars.to('mile')
print u.AU.to('kilometer')
print const.c.to('pc/yr')
print const.c.to('fur/fortnight')


224396806.05 km
139433710.91 mi
149597870.7
0.306601393788 pc / yr
1.80261749979e+12 fur / fortnight

To use this more practically, we can calculate the time it will take for light to reach the earth just by dividing 1 AU by the speed of light, as done below. Since AU is a unit, and c is in m/s, we end up with an answer that is (AU*m/s). By using .decompose() we can simplify that expression, which in this case will end up with an answer that is just in seconds. Finally, we can then convert that answer to minutes to get the answer of about 8 1/3 minutes that is commonly used. None of this required our doing the conversions where we might've slipped up.


In [91]:
time=1*u.AU/const.c
print(time)

time_s=time.decompose()
print(time_s)

time_min=time_s.to(u.minute)
print(time_min)


3.33564095198e-09 AU s / m
499.004783836 s
8.31674639727 min

Time

Astropy handles time in a similar way to units, with creating Time objects. These objects have two main properties.
The format is simply how the time is displayed. This is the difference between, for example, Julian Date, Modified Julian Date, and ISO time (YYYY-MM-DD HH:MM:SS). The second is the scale, and is the difference between terrestrial time vs time at the barycenter of the solar system.

We can start off by changing a time from one format to many others. We can also subtract times and we will get a timedelta unit.


In [92]:
from astropy.time import Time
t=Time(57867.346424, format='mjd', scale='utc')
t1=Time(58867.346424, format='mjd', scale='utc')
print t.mjd
print t.iso
print t.jyear
t1-t


57867.346424
2017-04-24 08:18:51.034
2017.31101006
Out[92]:
<TimeDelta object: scale='tai' format='jd' value=1000.0>

Coordinates

Coordinates again work by using an object time defined for this purpose. We can establish a point in the ICRS frame (this is approximately the equatorial coordinate) by defining the ra and dec. Note that here we are using u.degree in specifying the coordinates.

We can then print out the RA and dec, as well as change the units displayed. In the last line, we can also convert from ICRS equatorial coordinates to galactic coordinates.


In [93]:
c = SkyCoord(ra=10.68458*u.degree, dec=41.26917*u.degree, frame='icrs')
print c
print c.ra
print c.dec
print c.ra.hour
print c.ra.hms
print c.galactic


<SkyCoord (ICRS): (ra, dec) in deg
    ( 10.68458,  41.26917)>
10d41m04.488s
41d16m09.012s
0.712305333333
hms_tuple(h=0.0, m=42.0, s=44.299200000000525)
<SkyCoord (Galactic): (l, b) in deg
    ( 121.17424181, -21.57288557)>

Slightly practical application of this

Using some of these astropy functions, we can do some fancier applications. Starting off, we import a listing of stars with RA and dec from the attached table, and store them in the coordinate formats that are used by astropy. We then use matplotlib to plot this, and are able to easily convert them into radians thanks to astropy. This plot is accurate, but it lacks reference for where these points are.


In [94]:
hosts={}
data=np.loadtxt('./data/planets.tab', dtype='str', delimiter='\t')
print data[0]
hosts['ra_hours']=data[1:,9].astype(float)
hosts['ra']=data[1:,6].astype(float)
hosts['dec']=data[1:,8].astype(float)
#print hosts['ra_hours']
#print hosts['dec']

import astropy.units as u
import astropy.coordinates as coord
from astropy.coordinates import SkyCoord
ra = coord.Angle(hosts['ra']*u.degree)
ra = ra.wrap_at(180*u.degree)
dec = coord.Angle(hosts['dec']*u.degree)

fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111, projection="mollweide")
plt.title('Map of Exoplanets')
ax.scatter(ra.radian, dec.radian)
ax.set_xticklabels(['14h','16h','18h','20h','22h','0h','2h','4h','6h','8h','10h'])
ax.grid(True)
plt.show()


['rowid' 'pl_hostname' 'pl_letter' 'pl_discmethod' 'pl_pnum' 'ra_str' 'ra'
 'dec_str' 'dec' 'st_rah']

To fix this, we will add some references to this by adding a few more sets of data points. The first is relatively simple, we put in a line at the celestial equator. This just has to be a set of points that are all at declination of 0, and from -180 to +180 degrees in RA. These are a and b in the below.

We also want to add the planes of the ecliptic and the galaxy on this. For both, we use coordinate objects and provide numpy arrays where one coordinate is at zero, and the other goes from 0 to 360. With astropy we can then easily convert from each coordinate system to ICRS. There's some for loops to modify the plotting, but the important thing is that this will give us a plot that has not just the locations of all the planets that we've plotted, but will also include the celestial equator, galactic plane, and ecliptic plane on it.


In [95]:
a=coord.Angle((np.arange(361)-180)*u.degree)
b=coord.Angle(np.zeros(len(a))*u.degree)
numpoints=360
galaxy=SkyCoord(l=coord.Angle((np.arange(numpoints))*u.degree), b=coord.Angle(np.zeros(numpoints)*u.degree), frame='galactic')
ecliptic=SkyCoord(lon=coord.Angle((np.arange(numpoints))*u.degree), lat=coord.Angle(np.zeros(numpoints)*u.degree), frame='geocentrictrueecliptic')
ecl_eq=ecliptic.icrs
gal_eq=galaxy.icrs
#print gal_eq
fixed_ra=[]
for item in gal_eq.ra.radian:
   if item < np.pi:
      fixed_ra.append(item)
   else:
      fixed_ra.append(item-2*np.pi)
i=np.argmin(fixed_ra)
fixed_dec=[x for x in gal_eq.dec.radian]

fixed_ra_eq=[]
for item in ecl_eq.ra.radian:
   if item < np.pi:
      fixed_ra_eq.append(item)
   else:
      fixed_ra_eq.append(item-2*np.pi)
j=np.argmin(fixed_ra_eq)
fixed_dec_eq=[x for x in ecl_eq.dec.radian]

fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111, projection="mollweide")
plt.title('Map of Exoplanets')
ax.scatter(ra.radian, dec.radian)
ax.plot(a.radian, b.radian, color='r', lw=2)
#ax.scatter(gal_eq.ra.radian, gal_eq.dec.radian, color='g')
ax.plot(fixed_ra[i:]+fixed_ra[:i], fixed_dec[i:]+fixed_dec[:i], color='g', lw=2)
ax.plot(fixed_ra_eq[j:]+fixed_ra_eq[:j], fixed_dec_eq[j:]+fixed_dec_eq[:j], color='m', lw=2)
ax.set_xticklabels(['14h','16h','18h','20h','22h','0h','2h','4h','6h','8h','10h'])
ax.grid(True)
plt.show()


Reading in FITS files

One of the useful things with Astropy is that you can use it for reading in FITS files, and extracting info such as bands, exposure times, intrument information, etc.

In this example, we will read in a FITS image file, and extract its information


In [96]:
# We will use  `wget` to download the necessary file to the `data` folder.
!wget 'http://star.herts.ac.uk/~gb/python/656nmos.fits' -O ./data/hst_image.fits


--2017-04-24 09:54:34--  http://star.herts.ac.uk/~gb/python/656nmos.fits
Resolving star.herts.ac.uk... 147.197.221.254
Connecting to star.herts.ac.uk|147.197.221.254|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 10298880 (9.8M) [image/fits]
Saving to: `./data/hst_image.fits'

100%[======================================>] 10,298,880  4.86M/s   in 2.0s    

2017-04-24 09:54:36 (4.86 MB/s) - `./data/hst_image.fits' saved [10298880/10298880]

Now we can extract some of the information stored in the FITS file.


In [97]:
from astropy.io import fits
filename = './data/hst_image.fits'
hdulist = fits.open(filename)

The returned object, hdulist, (an instance of the HDUList class) behaves like a Python list, and each element maps to a Header-Data Unit (HDU) in the FITS file. You can view more information about the FITS file with:


In [98]:
hdulist.info()


Filename: ./data/hst_image.fits
No.    Name         Type      Cards   Dimensions   Format
  0  PRIMARY     PrimaryHDU     290   (1600, 1600)   float32   
  1  656nmos_cvt.tab  TableHDU       353   1R x 49C   [D25.17, D25.17, E15.7, E15.7, E15.7, E15.7, E15.7, E15.7, E15.7, E15.7, A1, E15.7, I12, I12, D25.17, D25.17, A8, A8, I12, E15.7, E15.7, E15.7, E15.7, E15.7, E15.7, I12, I12, I12, I12, I12, I12, I12, I12, A48, E15.7, E15.7, E15.7, E15.7, E15.7, E15.7, E15.7, E15.7, E15.7, E15.7, E15.7, E15.7, E15.7, E15.7, E15.7]   

As we can see, this file contains two HDUs. The first contains the image, the second a data table. To access the primary HDU, which contains the main data, you can then do:


In [99]:
hdu = hdulist[0]

To read the header of the FITS file, you can read hdulist. The following shows the different keys for the header


In [100]:
np.sort(hdulist[0].header.keys())


Out[100]:
array(['', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '',
       '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '',
       '', '', '', '', '', '', '', 'ALLG-MAX', 'ALLG-MIN', 'ASTR_I',
       'ATODCORR', 'ATODFILE', 'ATODGAIN', 'ATODSAT', 'BACKGRND',
       'BADPIXEL', 'BIASCORR', 'BIASDFIL', 'BIASEVEN', 'BIASFILE',
       'BIASODD', 'BITPIX', 'BLEVCORR', 'BLEVDFIL', 'BLEVFILE', 'BSCALE',
       'BZERO', 'CALIBDEF', 'CD1_1', 'CD1_2', 'CD2_1', 'CD2_2', 'CDBSFILE',
       'COMPTAB', 'CRPIX1', 'CRPIX2', 'CRVAL1', 'CRVAL2', 'CTYPE1',
       'CTYPE2', 'DADSCLAS', 'DADSDATE', 'DADSFILE', 'DARKCORR',
       'DARKDFIL', 'DARKFILE', 'DARKTIME', 'DATALOST', 'DATAMAX',
       'DATAMEAN', 'DATAMIN', 'DATE', 'DATE-OBS', 'DEC_SUN', 'DEC_TARG',
       'DETECTOR', 'DEZERO', 'DOHISTOS', 'DOPHOTOM', 'DOSATMAP',
       'EPLONGPM', 'EQNX_SUN', 'EQRADTRG', 'EQUINOX', 'ERRCNT', 'EXPEND',
       'EXPFLAG', 'EXPSTART', 'EXPTIME', 'EXTEND', 'FGSLOCK', 'FILENAME',
       'FILETYPE', 'FILLCNT', 'FILTER1', 'FILTER2', 'FILTNAM1', 'FILTNAM2',
       'FILTROT', 'FITSDATE', 'FLATCORR', 'FLATDFIL', 'FLATFILE',
       'FLATNTRG', 'FPKTTIME', 'FR_GAL1', 'FR_GAL2', 'FR_GAL3', 'FR_GAL4',
       'FR_GALS', 'FR_LOWF', 'FR_LOWF1', 'FR_LOWF2', 'FR_LOWF3',
       'FR_LOWF4', 'FR_STAR1', 'FR_STAR2', 'FR_STAR3', 'FR_STAR4',
       'FR_STARS', 'FTOTAL', 'FTOTAL1', 'FTOTAL2', 'FTOTAL3', 'FTOTAL4',
       'GOODMAX', 'GOODMIN', 'GPIXELS', 'GRAPHTAB', 'HISTORY', 'HISTORY',
       'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY',
       'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY',
       'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY',
       'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY',
       'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY',
       'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY',
       'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY',
       'HISTORY', 'HISTORY', 'HISTORY', 'HISTORY', 'HISTWIDE', 'IMAGETYP',
       'IMAG_I', 'INSTRUME', 'IRAF_I', 'KSPOTS', 'LINENUM', 'LONGPMER',
       'LPKTTIME', 'LRFWAVE', 'MASKCORR', 'MASKFILE', 'MEANC10',
       'MEANC100', 'MEANC200', 'MEANC25', 'MEANC300', 'MEANC50', 'MEDIAN',
       'MEDSHADO', 'MIR_REVR', 'MODE', 'MOONANGL', 'MTFLAG', 'NAXIS',
       'NAXIS1', 'NAXIS2', 'NCOMBINE', 'NPDECTRG', 'NPRATRG', 'NSHUTA17',
       'OCRVL1', 'OCRVL2', 'ODATTYPE', 'OPSIZE', 'ORIENTAT', 'ORIGIN',
       'OUTDTYPE', 'OVERLAP', 'PA_V3', 'PEP_EXPO', 'PHOTBW', 'PHOTFLAM',
       'PHOTMODE', 'PHOTPLAM', 'PHOTTAB', 'PHOTZPT', 'PKTFMT', 'PODPSFF',
       'PROPOSID', 'PSTPTIME', 'PSTRTIME', 'PYS_HEXP', 'PYS_ITER',
       'PYS_VERS', 'RA_SUN', 'RA_TARG', 'READTIME', 'ROOTNAME', 'ROTRTTRG',
       'RSDPFILL', 'SATURATE', 'SDASMGNU', 'SEQLINE', 'SEQNAME', 'SERIALS',
       'SEXC_I', 'SHADCORR', 'SHADFILE', 'SHUTTER', 'SIMPLE', 'SKEWNESS',
       'SKYSUB1', 'SKYSUB2', 'SKYSUB3', 'SKYSUB4', 'SOFTERRS', 'STAC_I',
       'STATICD', 'STDCFFF', 'STDCFFP', 'SUNANGLE', 'SUN_ALT', 'SURFALTD',
       'SURFLATD', 'SURFLONG', 'TARGNAME', 'TIME-OBS', 'UBAY3TMP',
       'UCH1CJTM', 'UCH2CJTM', 'UCH3CJTM', 'UCH4CJTM', 'UEXPODUR',
       'UEXPOTIM', 'USCALE', 'UZERO', 'XSHIFT', 'YSHIFT'], 
      dtype='|S8')

As we can see, this file contains two HDUs. The first contains the image, the second a data table.

Let's look at the image of the FITS file. The hdu object then has two important attributes: data, which behaves like a Numpy array, can be used to access the data, and header, which behaves like a dictionary, can be used to access the header information. First, we can take a look at the data:


In [101]:
hdu.data.shape


Out[101]:
(1600, 1600)

This tells us that it is a 1600-by-1600 pixel image. We can now take a peak at the header. To access the primary HDU, which contains the main data, you can then do:


In [102]:
hdu.header


Out[102]:
SIMPLE  =                    T / FITS STANDARD                                  
BITPIX  =                  -32 / FITS BITS/PIXEL                                
NAXIS   =                    2 / NUMBER OF AXES                                 
NAXIS1  =                 1600 /                                                
NAXIS2  =                 1600 /                                                
EXTEND  =                    T / There maybe standard extensions                
BSCALE  =                1.0E0 / REAL = TAPE*BSCALE + BZERO                     
BZERO   =                0.0E0 /                                                
OPSIZE  =                 2112 / PSIZE of original image                        
ORIGIN  = 'STScI-STSDAS'       / Fitsio version 21-Feb-1996                     
FITSDATE= '2005-07-01'         / Date FITS file was created                     
FILENAME= '656nmos_cvt.hhh'    / Original filename                              
ALLG-MAX=           0.000000E0 / Data max in all groups                         
ALLG-MIN=           0.000000E0 / Data min in all groups                         
ODATTYPE= 'FLOATING'           / Original datatype: Single precision real       
SDASMGNU=                    1 / Number of groups in original image             
CRVAL1  =      274.71149247724                                                  
CRVAL2  =     -13.816384007184                                                  
CRPIX1  =                386.5                                                  
CRPIX2  =                 396.                                                  
CD1_1   =          1.878013E-5                                                  
CD1_2   =         -2.031193E-5                                                  
CD2_1   =         -2.029358E-5                                                  
CD2_2   =         -1.879711E-5                                                  
DATAMIN =           0.000000E0 / DATA MIN                                       
DATAMAX =           0.000000E0 / DATA MAX                                       
MIR_REVR=                    T                                                  
ORIENTAT=            -131.9115                                                  
FILLCNT =                    0                                                  
ERRCNT  =                    0                                                  
FPKTTIME=     49808.7324543649                                                  
LPKTTIME=      49808.732622189                                                  
CTYPE1  = 'RA---TAN'                                                            
CTYPE2  = 'DEC--TAN'                                                            
DETECTOR=                    4                                                  
DEZERO  =             311.3041                                                  
BIASEVEN=             311.3347                                                  
BIASODD =             311.2736                                                  
GOODMIN =             20.02589                                                  
GOODMAX =             5584.467                                                  
DATAMEAN=              89.2569                                                  
GPIXELS =               559905                                                  
SOFTERRS=                    0                                                  
CALIBDEF=                80081                                                  
STATICD =                    0                                                  
ATODSAT =                   15                                                  
DATALOST=                    0                                                  
BADPIXEL=                    0                                                  
OVERLAP =                    0                                                  
PHOTMODE= 'WFPC2,4,A2D7,F656N,,CAL                         '                    
PHOTFLAM=         1.528990E-16                                                  
PHOTZPT =                -21.1                                                  
PHOTPLAM=             6563.585                                                  
PHOTBW  =             55.16187                                                  
MEDIAN  =             93.54508                                                  
MEDSHADO=             66.10242                                                  
HISTWIDE=             62.41553                                                  
SKEWNESS=           -0.1202522                                                  
MEANC10 =             126.9503                                                  
MEANC25 =             128.0693                                                  
MEANC50 =             134.1882                                                  
MEANC100=               132.23                                                  
MEANC200=             114.3746                                                  
MEANC300=             100.2731                                                  
BACKGRND=             42.37756                                                  
DADSFILE= 'U2LX0501T.D0F'      /                                                
DADSCLAS= 'CAL     '           /                                                
DADSDATE= '02-APR-1995 16:09:04' /                                              
                                                                                
                 / GROUP PARAMETERS: OSS                                        
                                                                                
                 / GROUP PARAMETERS: PODPS                                      
                                                                                
                 / GROUP PARAMETERS: DATA QUALITY FILE SUMMARY                  
                                                                                
                 / GROUP PARAMETERS: PHOTOMETRY                                 
                                                                                
                 / GROUP PARAMETERS: IMAGE STATISTICS                           
                                                                                
                 / WFPC-II DATA DESCRIPTOR KEYWORDS                             
INSTRUME= 'WFPC2             ' / instrument in use                              
ROOTNAME= 'U2LX0501T         ' / rootname of the observation set                
FILETYPE= 'SCI               ' / shp, ext, edq, sdq, sci                        
                                                                                
                 / SCIENCE INSTRUMENT CONFIGURATION                             
MODE    = 'FULL              ' / instr. mode: FULL (full res.), AREA (area int.)
SERIALS = 'ON                ' / serial clocks: ON, OFF                         
                                                                                
                 / IMAGE TYPE CHARACTERISTICS                                   
IMAGETYP= 'EXT               ' / DARK/BIAS/IFLAT/UFLAT/VFLAT/KSPOT/EXT/ECAL     
CDBSFILE= 'NO                ' / GENERIC/BIAS/DARK/FLAT/MASK/NO                 
PKTFMT  =                  104 / packet format code                             
DATE    = '02/04/95          ' / date file written (dd/mm/yy)                   
                                                                                
                 / FILTER CONFIGURATION                                         
FILTNAM1= 'F656N             ' / first filter name                              
FILTNAM2= '                  ' / second filter name                             
FILTER1 =                   31 / first filter number (0-48)                     
FILTER2 =                    0 / second filter number (0-48)                    
FILTROT =                  0.0 / partial filter rotation angle (degrees)        
LRFWAVE =                  0.0 / linear ramp filter wavelength                  
                                                                                
                 / INSTRUMENT STATUS USED IN DATA PROCESSING                    
UCH1CJTM=      -88.34862518311 / TEC cold junction #1 temperature (Celsius)     
UCH2CJTM=      -88.80734252930 / TEC cold junction #2 temperature (Celsius)     
UCH3CJTM=      -88.39449310303 / TEC cold junction #3 temperature (Celsius)     
UCH4CJTM=      -88.90411376953 / TEC cold junction #4 temperature (Celsius)     
UBAY3TMP=       13.61565399170 / bay 3 A1 temperature (deg C)                   
KSPOTS  = 'OFF               ' / Status of Kelsall spot lamps: ON, OFF          
SHUTTER = 'B                 ' / Shutter in place at the beginning of the exposu
ATODGAIN=       7.000000000000 / Analog to Digital Gain (Electrons/DN)          
                                                                                
                 / RSDP CONTROL KEYWORDS                                        
MASKCORR= 'COMPLETE'           / Do mask correction: PERFORM, OMIT, COMPLETE    
ATODCORR= 'COMPLETE'           / Do A-to-D correction: PERFORM, OMIT, COMPLETE  
BLEVCORR= 'COMPLETE'           / Do bias level correction: PERFORM, OMIT, COMPLE
BIASCORR= 'COMPLETE'           / Do bias correction: PERFORM, OMIT, COMPLETE    
DARKCORR= 'COMPLETE'           / Do dark correction: PERFORM, OMIT, COMPLETE    
FLATCORR= 'COMPLETE'           / Do flat field correction: PERFORM, OMIT, COMPLE
SHADCORR= 'COMPLETE'           / Do shaded shutter correction: PERFORM, OMIT, CO
DOSATMAP= 'OMIT              ' / Output saturated pixel map: PERFORM, OMIT, COMP
DOPHOTOM= 'COMPLETE'           / Fill photometry keywords: PERFORM, OMIT, COMPLE
DOHISTOS= 'OMIT    '           / Make histograms: PERFORM, OMIT, COMPLETE       
OUTDTYPE= 'REAL    '           / Output image datatype: REAL, LONG, SHORT       
                                                                                
                 / CALIBRATION REFERENCE FILES                                  
MASKFILE= 'uref$f8213081u.r0h' / name of the input DQF of known bad pixels      
ATODFILE= 'uref$dbu1405iu.r1h' / name of the A-to-D conversion file             
BLEVFILE= 'ucal$u2lx0501t.x0h' / Engineering file with extended register data   
BLEVDFIL= 'ucal$u2lx0501t.q1h' / Engineering file DQF                           
BIASFILE= 'uref$fcb1503du.r2h' / name of the bias frame reference file          
BIASDFIL= 'uref$fcb1503du.b2h' / name of the bias frame reference DQF           
DARKFILE= 'uref$f461210tu.r3h' / name of the dark reference file                
DARKDFIL= 'uref$f461210tu.b3h' / name of the dark reference DQF                 
FLATFILE= 'uref$g6h0945cu.r4h' / name of the flat field reference file          
FLATDFIL= 'uref$g6h0945cu.b4h' / name of the flat field reference DQF           
SHADFILE= 'uref$e371355iu.r5h' / name of the reference file for shutter shading 
PHOTTAB = '                  ' / name of the photometry calibration table       
GRAPHTAB= 'mtab$m1p1255om.tmg' / the HST graph table                            
COMPTAB = 'mtab$m4i1217km.tmc' / the HST components table                       
                                                                                
                 / DEFAULT KEYWORDS SET BY STSCI                                
SATURATE=                 4095 / Data value at which saturation occurs          
USCALE  =                  1.0 / Scale factor for output image                  
UZERO   =                  0.0 / Zero point for output image                    
                                                                                
                 / READOUT DURATION INFORMATION                                 
READTIME=                  464 / Length of time for CCD readout in clock ticks  
                                                                                
                 / PLANETARY SCIENCE KEYWORDS                                   
PA_V3   =        0.9299789E+02 / position angle of v3 axis of HST               
RA_SUN  =  0.1065970464620E+02 / right ascension of the sun (deg)               
DEC_SUN =  0.4585364087424E+01 / declination of the sun (deg)                   
EQNX_SUN=               2000.0 / equinox of the sun                             
MTFLAG  =                    F / moving target flag                             
EQRADTRG=                  0.0 / equatorial radius of target                    
FLATNTRG=                  0.0 / flattening of target                           
NPDECTRG=                  0.0 / north pole declination of target               
NPRATRG =                  0.0 / north pole right ascension of target           
ROTRTTRG=                  0.0 / rotation rate of target                        
LONGPMER=                  0.0 / longitude of prime meridian                    
EPLONGPM=                  0.0 / epoch of longitude of prime meridian           
SURFLATD=                  0.0 / surface feature latitude                       
SURFLONG=                  0.0 / surface feature longitude                      
SURFALTD=                  0.0 / surface feature altitude                       
                                                                                
                 / PODPS FILL VALUES                                            
PODPSFF =                    0 / 0=(no podps fill), 1=(podps fill present)      
STDCFFF =                    0 / 0=(no st dcf fill), 1=(st dcf fill present)    
STDCFFP = '0000              ' / st dcf fill pattern (hex)                      
RSDPFILL=                 -100 / bad data fill value for calibrated images      
                                                                                
                 / EXPOSURE TIME AND RELATED INFORMATION                        
UEXPODUR=                 1100 / Commanded duration of exposure (seconds)       
NSHUTA17=                    1 / Number of AP17 shutter B closes                
DARKTIME=       1100.000000000 / Dark time (seconds)                            
UEXPOTIM=                23408 / Major frame pulse time preceding exposure start
PSTRTIME= '1995.091:17:32:39 ' / Predicted obs. start time  (yyyy:ddd:hh:mm:ss) 
PSTPTIME= '1995.091:17:52:39 ' / Predicted obs. stop  time  (yyyy:ddd:hh:mm:ss) 
                                                                                
                 / EXPOSURE INFORMATION                                         
SUNANGLE=       96.85704040527 / angle between sun and V1 axis (deg)            
MOONANGL=       115.2464447021 / angle between moon and V1 axis (deg)           
SUN_ALT =       36.93915557861 / altitude of the sun above Earth's limb (deg)   
FGSLOCK = 'FINE              ' / commanded FGS lock (FINE,COARSE,GYROS,UNKNOWN) 
                                                                                
DATE-OBS= ' 1/04/95          ' / UT date of start of observation (dd/mm/yy)     
TIME-OBS= '17:15:17          ' / UT time of start of observation (hh:mm:ss)     
EXPSTART=       49808.71894742 / exposure start time (Modified Julian Date)     
EXPEND  =       49808.73167890 / exposure end time (Modified Julian Date)       
EXPTIME =                1100. / exposure duration (seconds)--calculated        
EXPFLAG = 'NORMAL            ' / Exposure interruption indicator                
                                                                                
                 / TARGET & PROPOSAL ID                                         
TARGNAME= 'M16-A             ' / proposer's target name                         
RA_TARG =  0.2747039166667E+03 / right ascension of the target (deg) (J2000)    
DEC_TARG= -0.1383091666667E+02 / declination of the target (deg) (J2000)        
                                                                                
PROPOSID=                05773 / PEP proposal identifier                        
PEP_EXPO= '4.2000017#001     ' / PEP exposure identifier including sequence     
LINENUM = '4.010             ' / PEP proposal line number                       
SEQLINE = '4.200             ' / PEP line number of defined sequence            
SEQNAME = 'M16SEQ            ' / PEP define/use sequence name                   
                                                                                
HISTORY   MASKFILE=uref$f8213081u.r0h  MASKCORR=COMPLETED                       
HISTORY   PEDIGREE=INFLIGHT 01/01/1994 - 15/05/1995                             
HISTORY   DESCRIP=STATIC MASK - INCLUDES CHARGE TRANSFER TRAPS                  
HISTORY   BIASFILE=uref$fcb1503du.r2h  BIASCORR=COMPLETED                       
HISTORY   PEDIGREE=INFLIGHT 13/01/1995 - 11/04/1995                             
HISTORY   DESCRIP=not signif. different from f1j16* but generated from new data 
HISTORY   DARKFILE=uref$f461210tu.r3h  DARKCORR=COMPLETED                       
HISTORY   PEDIGREE=INFLIGHT 27/03/1995 - 04/04/1995                             
HISTORY   DESCRIP=dark,full mode,serials on,gain=7,-88C                         
HISTORY   FLATFILE=uref$g6h0945cu.r4h  FLATCORR=COMPLETED                       
HISTORY   PEDIGREE=INFLIGHT 01/05/1994 - 01/10/1994                             
HISTORY   DESCRIP=Improved Cyc4 flat, fixed errors at CCD edges - now < 0.5% RMS
HISTORY   SHADFILE=uref$e371355iu.r5h  SHADCORR=COMPLETED                       
HISTORY   PEDIGREE=GROUND                                                       
HISTORY   DESCRIP=                                                              
HISTORY   PC1: bias jump level ~0.155 DN.                                       
HISTORY   The following throughput tables were used:                            
HISTORY   crotacomp$hst_ota_007_syn.fits, crwfpc2comp$wfpc2_optics_006_syn.fits,
HISTORY   crwfpc2comp$wfpc2_f656n_005_syn.fits,                                 
HISTORY   crwfpc2comp$wfpc2_dqepc1_005_syn.fits,                                
HISTORY   crwfpc2comp$wfpc2_a2d7pc1_004_syn.fits,                               
HISTORY   crwfpc2comp$wfpc2_flatpc1_003_syn.fits                                
HISTORY   WF2: bias jump level ~0.180 DN.                                       
HISTORY   The following throughput tables were used:                            
HISTORY   crotacomp$hst_ota_007_syn.fits, crwfpc2comp$wfpc2_optics_006_syn.fits,
HISTORY   crwfpc2comp$wfpc2_f656n_005_syn.fits,                                 
HISTORY   crwfpc2comp$wfpc2_dqewfc2_005_syn.fits,                               
HISTORY   crwfpc2comp$wfpc2_a2d7wf2_004_syn.fits,                               
HISTORY   crwfpc2comp$wfpc2_flatwf2_003_syn.fits                                
HISTORY   WF3: bias jump level ~0.274 DN.                                       
HISTORY   The following throughput tables were used:                            
HISTORY   crotacomp$hst_ota_007_syn.fits, crwfpc2comp$wfpc2_optics_006_syn.fits,
HISTORY   crwfpc2comp$wfpc2_f656n_005_syn.fits,                                 
HISTORY   crwfpc2comp$wfpc2_dqewfc3_005_syn.fits,                               
HISTORY   crwfpc2comp$wfpc2_a2d7wf3_004_syn.fits,                               
HISTORY   crwfpc2comp$wfpc2_flatwf3_003_syn.fits                                
HISTORY   WF4: bias jump level ~0.202 DN.                                       
HISTORY   The following throughput tables were used:                            
HISTORY   crotacomp$hst_ota_007_syn.fits, crwfpc2comp$wfpc2_optics_006_syn.fits,
HISTORY   crwfpc2comp$wfpc2_f656n_005_syn.fits,                                 
HISTORY   crwfpc2comp$wfpc2_dqewfc4_005_syn.fits,                               
HISTORY   crwfpc2comp$wfpc2_a2d7wf4_004_syn.fits,                               
HISTORY   crwfpc2comp$wfpc2_flatwf4_003_syn.fits                                
HISTORY   Ran task WARMPIX, version 1.1 (Sep 28, 1999) at Mon 00:16:49          
HISTORY   29-Jul-2002, rej_thresh=0.10000, fix_thresh=0.00300,                  
HISTORY   var_thresh=0.00300,                                                   
HISTORY   fix_dqval=1024, rej_val=INDEF                                         
SKYSUB1 =                 13.7 / Sky value for chip 1                           
PYS_VERS=                  1.7 / PyStack.py: Version                            
IRAF_I  = 'IRAF V2.11 May 1997 release:2.11.3b'                                 
STAC_I  = '$RCSfile: imstack_imshift_imcomb.cl,v $ $Revision: 1.1 $'            
PYS_ITER=                    5 / PyStack.py: Number of iterations               
PYS_HEXP=                    2 / PyStack.py: Value of exponent                  
SKYSUB2 =                86.37 / Sky value for chip 2                           
SKYSUB3 =                88.54 / Sky value for chip 3                           
SKYSUB4 =                44.66 / Sky value for chip 4                           
NCOMBINE=                    2 / Number of exposures for this association       
ASTR_I  = '$RCSfile: image_type.cl,v $ $Revision: 1.1 $'                        
XSHIFT  =         -2.136230E-4 / Y shift in degrees from USNO stars             
YSHIFT  =         -2.765656E-5 / X shift in degrees from USNO stars             
OCRVL1  =             274.7039 / Old CRVAL1 value before USNO correction        
OCRVL2  =            -13.83092 / Old CRVAL2 value before USNO correction        
FR_LOWF1=            0.9928686 / Fraction of flux in low frequency for chip 1   
FR_STAR1=          5.451115E-4 / Fraction of flux in stars for chip 1           
FR_GAL1 =          0.006586314 / Fraction of flux in stars for chip 1           
FTOTAL1 =             3886.887 / Total flux in chip 1                           
FR_LOWF2=            0.9486716 / Fraction of flux in low frequency for chip 2   
FR_STAR2=          0.004499074 / Fraction of flux in stars for chip 2           
FR_GAL2 =           0.04682927 / Fraction of flux in stars for chip 2           
FTOTAL2 =             12345.12 / Total flux in chip 2                           
FR_LOWF3=            0.9776544 / Fraction of flux in low frequency for chip 3   
FR_STAR3=          6.581650E-4 / Fraction of flux in stars for chip 3           
FR_GAL3 =           0.02168746 / Fraction of flux in stars for chip 3           
FTOTAL3 =             16654.92 / Total flux in chip 3                           
FR_LOWF4=            0.9892846 / Fraction of flux in low frequency for chip 4   
FR_STAR4=          0.001161799 / Fraction of flux in stars for chip 4           
FR_GAL4 =           0.00955356 / Fraction of flux in stars for chip 4           
FTOTAL4 =             24094.24 / Total flux in chip 4                           
FR_LOWF =            0.9773308 / Fraction of flux in low frequency for all chips
FR_STARS=          0.001695556 / Fraction of flux in stars for all chips        
FR_GALS =           0.02097363 / Fraction of flux in stars for all chips        
FTOTAL  =             56981.16 / Total flux in all chips                        
IMAG_I  = '$RCSfile$ $Revision$'                                                
EQUINOX =                2000.                                                  
SEXC_I  = '$RCSfile$ $Revision$'                                                
HISTORY   Ran task WMOSAIC, version 2.1 (Jun 1995), at Fri 16:17:34 01-Jul-2005 

We can access individual header keywords using standard item notation:


In [103]:
hdu.header['INSTRUME']


Out[103]:
'WFPC2'

In [104]:
hdu.header['EXPTIME']


Out[104]:
1100.0

We can plot the image using matplotlib:


In [105]:
plt.figure(figsize=(10,10))
plt.imshow(np.log10(hdu.data), origin='lower', cmap='gray', vmin=1.5, vmax=3)


/Users/victor2/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:2: RuntimeWarning: divide by zero encountered in log10
  from ipykernel import kernelapp as app
Out[105]:
<matplotlib.image.AxesImage at 0x118a9ae10>

You can also add new fields to the FITS file


In [106]:
hdu.header['MODIFIED'] = '2014-12-01'  # adds a new keyword

and we can also change the data, for example subtracting a background value:


In [107]:
hdu.data = hdu.data - 0.5

This only changes the FITS file in memory. You can write to a file with:


In [108]:
hdu.writeto('./data/hubble-image-background-subtracted.fits', overwrite=True)

In [109]:
!ls ./data


README.md                               astropy_data_latex.tex
astropy_data.tb                         hst_image.fits
astropy_data_2.tb                       hubble-image-background-subtracted.fits
astropy_data_csv.csv                    planets.tab
astropy_data_ecsv.ecsv                  sources.dat

Analytic Functions

Astropy comes with some built-in analytic functions, e.g. the blackbody radiation function.

Blackbody Radiation

Blackbody flux is calculated with Planck law (Rybicki & Lightman 1979)

$$B_{\lambda}(T) = \frac{2 h c^{2} / \lambda^{5}}{exp(h c / \lambda k T) - 1}$$$$B_{\nu}(T) = \frac{2 h \nu^{3} / c^{2}}{exp(h \nu / k T) - 1}$$

In [110]:
from astropy.analytic_functions import blackbody_lambda, blackbody_nu

In [111]:
def Planck_func(temp, lam_arr, opt='lam'):
    """
    Computes the Blackbody radiation curve of a blackbody of a given temperature `temp`.
    
    Parameters
    ----------
    temp: float or array-like
        temperature(s) of the blackbody
    
    lam_arr: float or array_like
        aray of wavelenths to evaluate the Planck function.
    
    opt: str, optional (default = 'lam')
        Option for returning either the flux of `lambda` (wavelength) or `nu` (frequency).
        Options:
            - `lam`: Return flux for `lambda' or wavelength
            - `nu` : Returns flux for `nu` (frequency)
    """
    wavelengths = lam_arr * u.AA
    temperature = temp * u.K
    with np.errstate(all='ignore'):
        flux_lam = blackbody_lambda(wavelengths, temperature)
        flux_nu  = blackbody_nu(wavelengths, temperature)
    
    if opt=='lam':
        return flux_lam
    if opt=='nu':
        return flux_nu

Let's plot the Planck function for two bodies with temperatures $T_1 = 8000\ K$ and $T_2 = 6000\ K$


In [112]:
lam_arr = np.arange(1e2, 2e4)
nu_arr  = (const.c/(lam_arr * u.AA)).to(1./u.s).value

In [113]:
fig = plt.figure(figsize=(15,8))
ax1  = fig.add_subplot(121, axisbg='white')
ax2  = fig.add_subplot(122, axisbg='white')
ax1.set_xlabel(r'$\lambda$ (Ansgstrom)', fontsize=25)
ax1.set_ylabel(r'$B_{\lambda}(T)$', fontsize=25)
ax2.set_xlabel(r'$\nu$ (\textrm{s}^{-1})', fontsize=25)
ax2.set_ylabel(r'$B_{\nu}(T)$', fontsize=25)
ax2.set_xscale('log')

temp_arr = [6e3, 8e3, 1e4, 1.2e4]
for temp in temp_arr:
    ax1.plot(lam_arr, Planck_func(temp, lam_arr=lam_arr, opt='lam'), label='T = {0} K'.format(int(temp)))
    ax2.plot(nu_arr , Planck_func(temp, lam_arr=lam_arr, opt='nu' ), label='T = {0} K'.format(int(temp)))
    ax1.legend(loc=1, prop={'size':20})


AstroPy Tables

Read files

You can use Astropy to read tables from data files. We'll use it to read the sources.dat file, which contains columns and rows of data


In [114]:
!head ./data/sources.dat


obsid redshift  X      Y     object
3102  0.32      4167  4085   Q1250+568-A
877   0.22      4378  3892   "Source 82"


In [115]:
from astropy.io import ascii
sources_tb = ascii.read('./data/sources.dat')

print( sources_tb )


obsid redshift  X    Y      object  
----- -------- ---- ---- -----------
 3102     0.32 4167 4085 Q1250+568-A
  877     0.22 4378 3892   Source 82

Write to files

You can also write directoy to a file using the data in the AstroPy table. Let's create a new AstroPy Table:


In [116]:
from astropy.table import Table, Column, MaskedColumn
x = np.random.uniform(low=10, high=20, size=(1000,))
y = np.random.uniform(low=100, high=50, size=(x.size,))
z = np.random.uniform(low=30, high=50, size=(x.size,))
data = Table([x, y], names=['x', 'y'])
print(data)


      x             y      
------------- -------------
15.2414976634 87.3902978029
14.2858102529 63.4778235534
12.6837892045 95.1038332346
15.7433361157 84.0554583013
10.6276318955 96.0939150693
13.6078233437 90.8772263701
16.7426639018 53.5680565112
18.5854346321 63.4464077849
16.8978153196 53.8207361937
11.3037331537 62.9365050838
          ...           ...
 18.942591823 76.8984462197
14.3360864888 93.3963823352
19.5899858716 54.7713213972
14.8603557394 65.0842136962
10.2954772658 86.2165081125
11.9427749909  98.733046694
10.8330331892 55.5301147739
16.8250678957 69.3254817244
 16.792578986 59.2414887155
12.4648161361 70.4455159566
10.8867007288  99.185886158
Length = 1000 rows

In [117]:
ascii.write(data, './data/astropy_data.tb', overwrite=True)

Let's see what's in the astropy_data.tb file


In [118]:
!head ./data/astropy_data.tb


x y
15.241497663441134 87.3902978028505
14.285810252905431 63.47782355343405
12.68378920454439 95.10383323455497
15.743336115743206 84.05545830126047
10.627631895531712 96.09391506928101
13.607823343675943 90.87722637008545
16.742663901777668 53.568056511155795
18.585434632105862 63.44640778487338
16.897815319630578 53.820736193723626

You can also specify the delimiter of the file. For example, we can separate it with a comma.


In [119]:
ascii.write(data, './data/astropy_data_2.tb', delimiter=',', overwrite=True)

In [120]:
!head ./data/astropy_data_2.tb


x,y
15.241497663441134,87.3902978028505
14.285810252905431,63.47782355343405
12.68378920454439,95.10383323455497
15.743336115743206,84.05545830126047
10.627631895531712,96.09391506928101
13.607823343675943,90.87722637008545
16.742663901777668,53.568056511155795
18.585434632105862,63.44640778487338
16.897815319630578,53.820736193723626

AstroPy Tables to other Formats

The AstroPy tables can also be converted to multiple formats

to Pandas DataFrames

A nice feature of AstroPy Tables is that you can export your data into different formats. For example, you can export it as a Pandas Dataframe.

See here for more info on how to use pandas with Astropy: http://docs.astropy.org/en/stable/table/pandas.html


In [121]:
df = data.to_pandas()
df.head()


Out[121]:
x y
0 15.241498 87.390298
1 14.285810 63.477824
2 12.683789 95.103833
3 15.743336 84.055458
4 10.627632 96.093915

And to compare, let's see the AstroPy Tables format


In [122]:
data


Out[122]:
<Table length=1000>
xy
float64float64
15.241497663487.3902978029
14.285810252963.4778235534
12.683789204595.1038332346
15.743336115784.0554583013
10.627631895596.0939150693
13.607823343790.8772263701
16.742663901853.5680565112
18.585434632163.4464077849
16.897815319653.8207361937
11.303733153762.9365050838
......
14.336086488893.3963823352
19.589985871654.7713213972
14.860355739465.0842136962
10.295477265886.2165081125
11.942774990998.733046694
10.833033189255.5301147739
16.825067895769.3254817244
16.79257898659.2414887155
12.464816136170.4455159566
10.886700728899.185886158

to LaTeX tables

A nice thing about AstroPy is that you can convert your data into LaTeX tables. This is easily done with writing it to a file. You can then copy it and use it on your next publication


In [123]:
import sys
ascii.write(data[0:10], sys.stdout, format='latex')


\begin{table}
\begin{tabular}{cc}
x & y \\
15.2414976634 & 87.3902978029 \\
14.2858102529 & 63.4778235534 \\
12.6837892045 & 95.1038332346 \\
15.7433361157 & 84.0554583013 \\
10.6276318955 & 96.0939150693 \\
13.6078233437 & 90.8772263701 \\
16.7426639018 & 53.5680565112 \\
18.5854346321 & 63.4464077849 \\
16.8978153196 & 53.8207361937 \\
11.3037331537 & 62.9365050838 \\
\end{tabular}
\end{table}

To save it as a file, you can do this:


In [124]:
ascii.write(data, './data/astropy_data_latex.tex', format='latex')

In [125]:
# I'm only showing the first 10 lines
!head ./data/astropy_data_latex.tex


\begin{table}
\begin{tabular}{cc}
x & y \\
15.2414976634 & 87.3902978029 \\
14.2858102529 & 63.4778235534 \\
12.6837892045 & 95.1038332346 \\
15.7433361157 & 84.0554583013 \\
10.6276318955 & 96.0939150693 \\
13.6078233437 & 90.8772263701 \\
16.7426639018 & 53.5680565112 \\

to CSV files


In [126]:
ascii.write(data, './data/astropy_data_csv.csv', format='csv', fast_writer=False)

In [127]:
!head ./data/astropy_data_csv.csv


x,y
15.2414976634,87.3902978029
14.2858102529,63.4778235534
12.6837892045,95.1038332346
15.7433361157,84.0554583013
10.6276318955,96.0939150693
13.6078233437,90.8772263701
16.7426639018,53.5680565112
18.5854346321,63.4464077849
16.8978153196,53.8207361937

Other formats

AstroPy tables come with a great support for many different types of files. This is a list of the supported files that you can import/export AstroPy tables.

Data tables and Column types

You can also use AstroPy tables to preserve the metadata of a column. For example, you can keep the units of each column, so that you use the data later on, and still be able to use unit conversions, etc. for this.


In [128]:
t = Table(masked=True)
t['x'] = MaskedColumn([1.0, 2.0], unit='m', dtype='float32')
t['x'][1] = np.ma.masked
t['y'] = MaskedColumn([False, True], dtype='bool')

In [129]:
t


Out[129]:
<Table masked=True length=2>
xy
m
float32bool
1.0False
--True

Now we can save it into a ecsv file. This type of file will preserve the type of units, and more, for each of the columns


In [130]:
from astropy.extern.six.moves import StringIO
fh = StringIO()
t.write(fh, format='ascii.ecsv')  
table_string = fh.getvalue()      
print(table_string)


# %ECSV 0.9
# ---
# datatype:
# - {name: x, unit: m, datatype: float32}
# - {name: y, datatype: bool}
x y
1.0 False
"" True


In [131]:
Table.read(table_string, format='ascii')


Out[131]:
<Table masked=True length=2>
xy
m
float32bool
1.0False
--True

Or you can dump it into a file


In [132]:
t.write('./data/astropy_data_ecsv.ecsv', format='ascii.ecsv', overwrite=True)

And you can now read it in


In [133]:
data_ecsv = ascii.read('./data/astropy_data_ecsv.ecsv', format='ecsv')
data_ecsv


Out[133]:
<Table masked=True length=2>
xy
m
float32bool
1.0False
--True

In [134]:
data_ecsv['x']


Out[134]:
<MaskedColumn name='x' dtype='float32' unit='m' length=2>
1.0
--

Resources

For further reading and exercises, you can check out: