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:
The other is file type and structures:
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.
# Importing Modules
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
sns.set_context("notebook")
import astropy
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy import constants as const
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.
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
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.
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
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.
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
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.
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
<TimeDelta object: scale='tai' format='jd' value=1000.0>
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.
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)>
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.
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.
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()
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
# 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.
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:
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:
hdu = hdulist[0]
To read the header of the FITS file, you can read hdulist
. The following shows the different keys for the header
np.sort(hdulist[0].header.keys())
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:
hdu.data.shape
(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:
hdu.header
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:
hdu.header['INSTRUME']
'WFPC2'
hdu.header['EXPTIME']
1100.0
We can plot the image using matplotlib:
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
<matplotlib.image.AxesImage at 0x118a9ae10>
You can also add new fields to the FITS file
hdu.header['MODIFIED'] = '2014-12-01' # adds a new keyword
and we can also change the data, for example subtracting a background value:
hdu.data = hdu.data - 0.5
This only changes the FITS file in memory. You can write to a file with:
hdu.writeto('./data/hubble-image-background-subtracted.fits', overwrite=True)
!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
Astropy
comes with some built-in analytic functions, e.g. the blackbody radiation
function.
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}$$from astropy.analytic_functions import blackbody_lambda, blackbody_nu
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$
lam_arr = np.arange(1e2, 2e4)
nu_arr = (const.c/(lam_arr * u.AA)).to(1./u.s).value
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})
!head ./data/sources.dat
obsid redshift X Y object 3102 0.32 4167 4085 Q1250+568-A 877 0.22 4378 3892 "Source 82"
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
You can also write directoy to a file using the data in the AstroPy
table.
Let's create a new AstroPy
Table:
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
ascii.write(data, './data/astropy_data.tb', overwrite=True)
Let's see what's in the astropy_data.tb
file
!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.
ascii.write(data, './data/astropy_data_2.tb', delimiter=',', overwrite=True)
!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
The AstroPy
tables can also be converted to multiple formats
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
df = data.to_pandas()
df.head()
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
data
x | y |
---|---|
float64 | float64 |
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 |
... | ... |
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 |
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
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:
ascii.write(data, './data/astropy_data_latex.tex', format='latex')
# 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 \\
ascii.write(data, './data/astropy_data_csv.csv', format='csv', fast_writer=False)
!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
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.
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.
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')
t
x | y |
---|---|
m | |
float32 | bool |
1.0 | False |
-- | 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
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
Table.read(table_string, format='ascii')
x | y |
---|---|
m | |
float32 | bool |
1.0 | False |
-- | True |
Or you can dump it into a file
t.write('./data/astropy_data_ecsv.ecsv', format='ascii.ecsv', overwrite=True)
And you can now read it in
data_ecsv = ascii.read('./data/astropy_data_ecsv.ecsv', format='ecsv')
data_ecsv
x | y |
---|---|
m | |
float32 | bool |
1.0 | False |
-- | True |
data_ecsv['x']
1.0 |
-- |
For further reading and exercises, you can check out: