A Community Python Library for Astronomy

Astropy is a bundle of several kinds of functionality for Astronomy. It consists of the core packages, providing base functionality for everyone, and affiliated packages which ontain more specialized functionality. The core packages have a minimum of dependencies, while affiliated packages are allowed to depend on external libraries.

A selection of core packages:

  • astropy.constants
  • astropy.units
  • astropy.nddata - Numpy arrays on steroids
  • astropy.io - file input/output, i.e. FITS and ASCII
  • astropy.convolution - smoothing, filtering etc.
  • astropy.cosmology
  • astropy.modeling - modelling and fitting

Examples of affiliated packages:

  • Ginga - FITS image viewer based on Astropy
  • Montage-wrapper - Advanced image handling
  • APLPy ("Apple-pie") - Astronomical image plotting
  • astroquery - Get data from online databases

...etc. etc.

Obviously, there is way too much functionality in Astropy to do more than just scratch the paint. We have already looked at astropy.io.fits. We will now take a little look at:

  • astropy.units
  • astropy.cosmology

astropy.units

In [33]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from astropy import units as u
from astropy import constants as con 
# Only when running in n IPython notebook:
%matplotlib inline  
In [34]:
# Just for good looks, no need to worry about this now.
# You are welcome to do this but it is not necessary.
c1 = cm.Reds(0.6)
c2 = cm.Blues(0.6)
c3 = cm.PuOr(0.3)
plt.rc('xtick', labelsize=13)
plt.rc('ytick', labelsize=13)
plt.rc('font', size = 12)
plt.rc('lines', linewidth=3)
plt.rc('figure', figsize=(10.0, 7.0))
plt.rcParams['axes.formatter.limits'] = (-3,3) #
In [38]:
dist = 3.4159 * u.km
dist
Out[38]:
$3.4159 \; \mathrm{km}$
In [44]:
time = 2.7183 * u.h
time
Out[44]:
$2.7183 \; \mathrm{h}$
In [45]:
dist / time
Out[45]:
$1.25663 \; \mathrm{\frac{km}{h}}$
In [47]:
(_).to(u.km / u.s)
Out[47]:
$0.000349064 \; \mathrm{\frac{km}{s}}$

We can of course define custom units for our convenience:

In [49]:
kmph = u.Unit(u.km / u.hour)
kms = u.Unit(u.km / u.s)
mph = u.Unit(u.imperial.mile / u.hour)
(200 * kmph).to(mph)
Out[49]:
$124.274 \; \mathrm{\frac{mi}{h}}$

It handles NumPy arrays as well as single scalars:

In [36]:
speeds = np.arange(100, 200) * kmph
print(speeds.to(mph))
[  62.13711922   62.75849042   63.37986161   64.0012328    64.62260399
   65.24397518   65.86534638   66.48671757   67.10808876   67.72945995
   68.35083115   68.97220234   69.59357353   70.21494472   70.83631592
   71.45768711   72.0790583    72.70042949   73.32180068   73.94317188
   74.56454307   75.18591426   75.80728545   76.42865665   77.05002784
   77.67139903   78.29277022   78.91414141   79.53551261   80.1568838
   80.77825499   81.39962618   82.02099738   82.64236857   83.26373976
   83.88511095   84.50648214   85.12785334   85.74922453   86.37059572
   86.99196691   87.61333811   88.2347093    88.85608049   89.47745168
   90.09882287   90.72019407   91.34156526   91.96293645   92.58430764
   93.20567884   93.82705003   94.44842122   95.06979241   95.6911636
   96.3125348    96.93390599   97.55527718   98.17664837   98.79801957
   99.41939076  100.04076195  100.66213314  101.28350433  101.90487553
  102.52624672  103.14761791  103.7689891   104.3903603   105.01173149
  105.63310268  106.25447387  106.87584506  107.49721626  108.11858745
  108.73995864  109.36132983  109.98270103  110.60407222  111.22544341
  111.8468146   112.46818579  113.08955699  113.71092818  114.33229937
  114.95367056  115.57504176  116.19641295  116.81778414  117.43915533
  118.06052653  118.68189772  119.30326891  119.9246401   120.54601129
  121.16738249  121.78875368  122.41012487  123.03149606  123.65286726] mi / h
In [5]:
print speeds.to(u.lightyear/u.millisecond)
[  2.93611343e-18   2.96547456e-18   2.99483570e-18   3.02419683e-18
   3.05355796e-18   3.08291910e-18   3.11228023e-18   3.14164137e-18
   3.17100250e-18   3.20036364e-18   3.22972477e-18   3.25908590e-18
   3.28844704e-18   3.31780817e-18   3.34716931e-18   3.37653044e-18
   3.40589158e-18   3.43525271e-18   3.46461384e-18   3.49397498e-18
   3.52333611e-18   3.55269725e-18   3.58205838e-18   3.61141952e-18
   3.64078065e-18   3.67014178e-18   3.69950292e-18   3.72886405e-18
   3.75822519e-18   3.78758632e-18   3.81694746e-18   3.84630859e-18
   3.87566972e-18   3.90503086e-18   3.93439199e-18   3.96375313e-18
   3.99311426e-18   4.02247540e-18   4.05183653e-18   4.08119766e-18
   4.11055880e-18   4.13991993e-18   4.16928107e-18   4.19864220e-18
   4.22800334e-18   4.25736447e-18   4.28672560e-18   4.31608674e-18
   4.34544787e-18   4.37480901e-18   4.40417014e-18   4.43353128e-18
   4.46289241e-18   4.49225354e-18   4.52161468e-18   4.55097581e-18
   4.58033695e-18   4.60969808e-18   4.63905922e-18   4.66842035e-18
   4.69778148e-18   4.72714262e-18   4.75650375e-18   4.78586489e-18
   4.81522602e-18   4.84458716e-18   4.87394829e-18   4.90330942e-18
   4.93267056e-18   4.96203169e-18   4.99139283e-18   5.02075396e-18
   5.05011510e-18   5.07947623e-18   5.10883736e-18   5.13819850e-18
   5.16755963e-18   5.19692077e-18   5.22628190e-18   5.25564304e-18
   5.28500417e-18   5.31436530e-18   5.34372644e-18   5.37308757e-18
   5.40244871e-18   5.43180984e-18   5.46117098e-18   5.49053211e-18
   5.51989324e-18   5.54925438e-18   5.57861551e-18   5.60797665e-18
   5.63733778e-18   5.66669892e-18   5.69606005e-18   5.72542118e-18
   5.75478232e-18   5.78414345e-18   5.81350459e-18   5.84286572e-18] lyr / ms
In [6]:
type(speeds)
Out[6]:
astropy.units.quantity.Quantity

Astropy's Quantity objects plays nicely with the rest of the SciPy ecosystem by disguising themselves as numpy.ndarray's. No need to convert before sending them to e.g. Matplotlib. They also play nicely with most simple numpy functions:

In [7]:
plt.plot(speeds, np.sqrt(speeds.to(mph)))
Out[7]:
[<matplotlib.lines.Line2D at 0x7f8e6a278f90>]

...But not all, which makes sense because log(km/h) is meaningless:

In [8]:
np.sqrt(speeds)
Out[8]:
<Quantity [ 10.        , 10.04987562, 10.09950494, 10.14889157,
            10.19803903, 10.24695077, 10.29563014, 10.34408043,
            10.39230485, 10.44030651, 10.48808848, 10.53565375,
            10.58300524, 10.63014581, 10.67707825, 10.72380529,
            10.77032961, 10.81665383, 10.86278049, 10.90871211,
            10.95445115, 11.        , 11.04536102, 11.09053651,
            11.13552873, 11.18033989, 11.22497216, 11.26942767,
            11.3137085 , 11.35781669, 11.40175425, 11.44552314,
            11.48912529, 11.53256259, 11.5758369 , 11.61895004,
            11.66190379, 11.70469991, 11.74734012, 11.78982612,
            11.83215957, 11.87434209, 11.91637529, 11.95826074,
            12.        , 12.04159458, 12.08304597, 12.12435565,
            12.16552506, 12.20655562, 12.24744871, 12.28820573,
            12.32882801, 12.36931688, 12.40967365, 12.4498996 ,
            12.489996  , 12.52996409, 12.56980509, 12.60952021,
            12.64911064, 12.68857754, 12.72792206, 12.76714533,
            12.80624847, 12.84523258, 12.88409873, 12.92284798,
            12.9614814 , 13.        , 13.03840481, 13.07669683,
            13.11487705, 13.15294644, 13.19090596, 13.22875656,
            13.26649916, 13.3041347 , 13.34166406, 13.37908816,
            13.41640786, 13.45362405, 13.49073756, 13.52774926,
            13.56465997, 13.60147051, 13.6381817 , 13.67479433,
            13.7113092 , 13.74772708, 13.78404875, 13.82027496,
            13.85640646, 13.89244399, 13.92838828, 13.96424004,
            14.        , 14.03566885, 14.07124728, 14.10673598] km(1/2) / h(1/2)>
In [9]:
np.log10(speeds)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-9-abd8bbabc22f> in <module>()
----> 1 np.log10(speeds)

/home/trive/anaconda/lib/python2.7/site-packages/astropy/units/quantity.pyc in __array_prepare__(self, obj, context)
    281         # the unit the output from the ufunc will have.
    282         if function in UFUNC_HELPERS:
--> 283             scales, result_unit = UFUNC_HELPERS[function](function, *units)
    284         else:
    285             raise TypeError("Unknown ufunc {0}.  Please raise issue on "

/home/trive/anaconda/lib/python2.7/site-packages/astropy/units/quantity_helper.pyc in helper_dimensionless_to_dimensionless(f, unit)
     87         raise TypeError("Can only apply '{0}' function to "
     88                         "dimensionless quantities"
---> 89                         .format(f.__name__))
     90     return [scale], dimensionless_unscaled
     91 

TypeError: Can only apply 'log10' function to dimensionless quantities

The solution to this is a little annoying but there's no way around it:

In [10]:
logspeeds = np.log10(speeds.to(mph).value)  # No units anymore!
print(logspeeds)
plt.plot(speeds, logspeeds)
plt.show()
[ 1.79335111  1.79767249  1.80195129  1.80618834  1.81038445  1.81454041
  1.81865698  1.82273489  1.82677487  1.83077761  1.8347438   1.83867409
  1.84256914  1.84642956  1.85025597  1.85404896  1.8578091   1.86153698
  1.86523312  1.86889808  1.87253236  1.87613649  1.87971095  1.88325623
  1.8867728   1.89026113  1.89372166  1.89715484  1.90056108  1.90394083
  1.90729447  1.91062241  1.91392505  1.91720276  1.92045591  1.92368488
  1.92689002  1.93007168  1.9332302   1.93636592  1.93947915  1.94257023
  1.94563946  1.94868715  1.95171361  1.95471912  1.95770397  1.96066845
  1.96361283  1.96653738  1.96944237  1.97232806  1.9751947   1.97804255
  1.98087184  1.98368281  1.98647571  1.98925077  1.9920082   1.99474824
  1.9974711   2.00017699  2.00286613  2.00553872  2.00819496  2.01083506
  2.0134592   2.01606759  2.0186604   2.02123782  2.02380004  2.02634723
  2.02887956  2.03139722  2.03390036  2.03638916  2.03886378  2.04132438
  2.04377112  2.04620415  2.04862362  2.05102969  2.0534225   2.0558022
  2.05816894  2.06052284  2.06286406  2.06519272  2.06750896  2.06981292
  2.07210472  2.07438448  2.07665234  2.07890842  2.08115284  2.08338573
  2.08560719  2.08781734  2.09001631  2.09220419]

astropy.cosmology

In [11]:
from astropy.cosmology import Planck13, LambdaCDM
In [12]:
pl13 = Planck13
redshifts = np.linspace(0.001, 100, 100)
In [13]:
pl13.lookback_time(7.8)
Out[13]:
$13.1355 \; \mathrm{Gyr}$
In [14]:
pl13.angular_diameter_distance(redshifts)
Out[14]:
<Quantity [    4.41823263, 1701.89483502, 1768.15666645, 1622.24519705,
            1461.51845089, 1319.11084429, 1198.20127448, 1096.10255984,
            1009.43135739,  935.23098014,  871.1271297 ,  815.25727217,
             766.16458236,  722.70378618,  683.96706805,  649.22804468,
             617.89990425,  589.50413756,  563.64706083,  540.00205354,
             518.29600201,  498.2988584 ,  479.81552582,  462.67949537,
             446.74781434,  431.89707476,  418.02019063,  405.02378999,
             392.82608989,  381.35515352,  370.54745206,  360.34667099,
             350.70271389,  341.5708667 ,  332.91109316,  324.68743805,
             316.86751948,  309.42209513,  302.32469024,  295.55127716,
             289.07999861,  282.8909275 ,  276.96585816,  271.28812395,
             265.84243768,  260.61475154,  255.59213367,  250.76265928,
             246.11531429,  241.6399098 ,  237.32700606,  233.16784473,
             229.15428827,  225.27876584,  221.53422466,  217.91408632,
             214.4122075 ,  211.02284444,  207.74062096,  204.56049941,
             201.47775434,  198.48794871,  195.58691206,  192.77072083,
             190.03568023,  187.37830781,  184.79531837,  182.28361017,
             179.84025223,  177.46247274,  175.14764838,  172.89329443,
             170.6970558 ,  168.55669856,  166.47010232,  164.43525302,
             162.45023638,  160.51323176,  158.62250648,  156.77641056,
             154.97337186,  153.2118915 ,  151.49053965,  149.80795155,
             148.16282391,  146.55391143,  144.9800236 ,  143.44002173,
             141.93281615,  140.4573636 ,  139.01266476,  137.59776197,
             136.21173707,  134.85370939,  133.52283384,  132.21829911,
             130.93932603,  129.68516599,  128.45509941,  127.24843438] Mpc>
In [15]:
plt.plot(redshifts, pl13.angular_diameter_distance(redshifts))
Out[15]:
[<matplotlib.lines.Line2D at 0x7f8e668c3150>]
In [16]:
pl13.comoving_distance(redshifts).to(u.lightyear)
Out[16]:
<Quantity [  1.44247579e+07,  1.11632409e+10,  1.74230218e+10,
             2.13296904e+10,  2.40313519e+10,  2.60355611e+10,
             2.75965894e+10,  2.88561678e+10,  2.98999955e+10,
             3.07832293e+10,  3.15431515e+10,  3.22059726e+10,
             3.27907228e+10,  3.33115945e+10,  3.37794155e+10,
             3.42026086e+10,  3.45878370e+10,  3.49404502e+10,
             3.52647993e+10,  3.55644639e+10,  3.58424199e+10,
             3.61011641e+10,  3.63428092e+10,  3.65691562e+10,
             3.67817516e+10,  3.69819313e+10,  3.71708567e+10,
             3.73495423e+10,  3.75188792e+10,  3.76796537e+10,
             3.78325625e+10,  3.79782253e+10,  3.81171958e+10,
             3.82499701e+10,  3.83769943e+10,  3.84986709e+10,
             3.86153640e+10,  3.87274038e+10,  3.88350910e+10,
             3.89386994e+10,  3.90384795e+10,  3.91346608e+10,
             3.92274539e+10,  3.93170526e+10,  3.94036356e+10,
             3.94873677e+10,  3.95684017e+10,  3.96468788e+10,
             3.97229304e+10,  3.97966785e+10,  3.98682366e+10,
             3.99377106e+10,  4.00051995e+10,  4.00707957e+10,
             4.01345857e+10,  4.01966505e+10,  4.02570663e+10,
             4.03159045e+10,  4.03732323e+10,  4.04291130e+10,
             4.04836063e+10,  4.05367683e+10,  4.05886522e+10,
             4.06393081e+10,  4.06887835e+10,  4.07371235e+10,
             4.07843704e+10,  4.08305648e+10,  4.08757449e+10,
             4.09199472e+10,  4.09632062e+10,  4.10055549e+10,
             4.10470246e+10,  4.10876451e+10,  4.11274448e+10,
             4.11664510e+10,  4.12046894e+10,  4.12421849e+10,
             4.12789610e+10,  4.13150403e+10,  4.13504443e+10,
             4.13851939e+10,  4.14193088e+10,  4.14528080e+10,
             4.14857096e+10,  4.15180312e+10,  4.15497894e+10,
             4.15810003e+10,  4.16116794e+10,  4.16418416e+10,
             4.16715010e+10,  4.17006714e+10,  4.17293660e+10,
             4.17575976e+10,  4.17853783e+10,  4.18127199e+10,
             4.18396339e+10,  4.18661311e+10,  4.18922222e+10,
             4.19179173e+10] lyr>
In [17]:
plt.plot(redshifts, pl13.comoving_distance(redshifts).to(u.lightyear))
Out[17]:
[<matplotlib.lines.Line2D at 0x7f8e667dded0>]
In [18]:
pl13.angular_diameter_distance(8.4)
Out[18]:
$984.829 \; \mathrm{Mpc}$
In [19]:
plt.plot(redshifts, pl13.luminosity_distance(redshifts))
Out[19]:
[<matplotlib.lines.Line2D at 0x7f8e6678fb50>]
In [20]:
redshift = 0.3
pl13.comoving_distance(redshift)
Out[20]:
$1231.38 \; \mathrm{Mpc}$

Question: What is the current recession velocity at this distance?

In [21]:
def vel_now(z):
    return pl13.H0 * pl13.comoving_distance(z)

vel_now(redshift)
Out[21]:
$83450.9 \; \mathrm{\frac{km}{s}}$

But not to worry, the more interesting question is: What was the recession velocity at the time of emission:

In [22]:
def vel_em(z): 
    return pl13.H(z) * pl13.scale_factor(z) * pl13.comoving_distance(z)

vel_em(redshift)
Out[22]:
$75121.1 \; \mathrm{\frac{km}{s}}$
In [23]:
pl13.H(redshift)
Out[23]:
$79.307 \; \mathrm{\frac{km}{Mpc\,s}}$
In [24]:
pl13.scale_factor(redshift)
Out[24]:
0.7692307692307692
In [25]:
pl13.comoving_distance(redshift)
Out[25]:
$1231.38 \; \mathrm{Mpc}$
In [31]:
plt.plot((redshifts), vel_em(redshifts), label='Velovity now')
plt.plot((redshifts), vel_now(redshifts), label='Velocity at emission')
plt.axhline(con.c.to(kms).value, color='k', ls='--')
plt.legend(loc='best')
#plt.axis((-.1, 3, 0, 4e5))
#plt.plot(redshifts, vel_em(redshifts))
Out[31]:
<matplotlib.legend.Legend at 0x7f8e66455810>
In [177]:
MatDom = LambdaCDM(0.7, 0.3, 0)
MatDom
Out[177]:
LambdaCDM(H0=0.7 km / (Mpc s), Om0=0.3, Ode0=0, Tcmb0=2.725 K, Neff=3.04, m_nu=[ 0.  0.  0.] eV)
In [30]:
 
In [ ]: