Astrometry in Python with PyEphem
[Last updated: 2010/12/04]
- Date and time
- Positions in the sky and coordinate systems
- Constants and utility functions
PyEphem, by Brandon Craig Rhodes, is a python package for performing astrometric calculations. PyEphem can be used to perform calculations such as, rising and setting times for astronomical objects and Earth satellites, Julian dates, sidereal time, astronomical coordinate of objects at various epochs, conversion between different coordinate systems, and many others. It uses routines from the XEphem package.
PyEphem can be downloaded from the PyEphem home page. It is available for both python 2.x and 3.x releases.
PyEphem comes with very good documentation that includes a large collection of examples. This document aims to complement the official PyEphem documentation.
In this document, we start by introducing the fundamental concepts used in PyEphem: angle, date and time, representation of positions in the sky and coordinate systems, body, and observer. We then look at the various utility methods provided in the package. Each of the above sections also provide examples that illustrate the concepts introduced. We end with a few examples that illustrate the usage of PyEphem in performing astronomical calculations.
PyEphem uses object oriented programming techniques and as such, provides most of its facilities as methods of relevant objects. For example, a particular equatorial coordinate can be represented using an ephem.Equatorial object. To find the right ascension of this coordinate we can use the ra property of the object, as shown below:
>>> import ephem >>> # ra = 12 hours, dec = 20 degrees and 30 minutes >>> eq = ephem.Equatorial("12:0:0", "20:30:00") >>> print eq.ra, eq.dec 12:00:00.00 20:30:00.0
Angles are internally represented in radians. But when converted into strings, using print or str(), the angle is converted into degrees or hours. The latter unit is used for right ascension, where 24 hours equals 360 degrees.
We can create an Angle object using the methods ephem.hours() and ephem.degrees(). We can specify an input angle, for both of the above methods, either as a floating point number or as a string. In the former case, the input is considered to be in radians. In the latter case, data in the string is interpreted as a colon separated string of the format dd:mm:ss by the ephem.degrees() method, and as a colon separated string of the format hh:mm:ss by the ephem.hours() method. Only the dd part of the degree string and hh part of the hour string are mandatory, and they can be floating point numbers.
Some examples are given below:
>>> import ephem >>> # input in radians; π radians == 12 hours >>> a = ephem.hours(ephem.pi) >>> # input in radians; π radians == 180 degrees >>> b = ephem.degrees(ephem.pi) >>> # input in degrees >>> c = ephem.degrees("180.0") >>> c1 = ephem.degrees("180:00") >>> c2 = ephem.degrees("180:00:00") >>> # input in hours >>> d = ephem.hours("12.0") >>> d1 = ephem.hours("12:00") >>> d2 = ephem.hours("12:00:00") >>> # print internal representation; shows that internally all >>> # the above represent the same value. >>> print repr(a), repr(b) 3.1415926535897931 3.1415926535897931 >>> print repr(c), repr(c1), repr(c2) 3.1415926535897931 3.1415926535897931 3.1415926535897931 >>> print repr(d), repr(d1), repr(d2) 3.1415926535897931 3.1415926535897931 3.1415926535897931 >>> # print interpreted values in a nice format; shows that the values >>> # are interpreted in different ways. >>> print a, b 12:00:00.00 180:00:00.0 >>> print c, c1, c2 180:00:00.0 180:00:00.0 180:00:00.0 >>> print d, d1, d2 12:00:00.00 12:00:00.00 12:00:00.00
The result of mathematical operations on Angle objects is always a floating point number in radians. Use one of the above two methods to convert the results into hours or degrees, as appropriate.
>>> # input in radians; pi radians == 12 hours >>> a = ephem.hours(ephem.pi) >>> # ephem.pi/4.0 == 3 hours == 45 degrees >>> print a + ephem.degrees(ephem.pi/4.0) 3.92699081699 >>> print ephem.hours(a+ephem.degrees(ephem.pi/4.0)) 15:00:00.00 >>> print ephem.degrees(a+ephem.degrees(ephem.pi/4.0)) 225:00:00.0
To get angles restricted to the range [0,2π) and [-π, π), we can use the norm and znorm properties, respectively, of a date object.
>>> print ephem.degrees(3*ephem.pi).norm 180:00:00.0 >>> print ephem.degrees(3*ephem.pi).znorm -180:00:00.0
The constant ephem.degree gives the number of radians per degree, so that ephem.degree * 180.0 gives π (3.14…) radians.
Date and time in PyEphem are stored using floating point numbers, that denote the number of days passed since 1899/12/31 12:00:00 UTC. So 1899/12/31 12:00:00 UTC is 0.0, 1900/1/1 00:00:00 UTC is 0.5, 2010/01/24 08:30:16.4 UTC is 40200.854356481483 and so on. Negative numbers are used to represent time before 1899/12/31 12:00:00 UTC.
PyEphem represents date and time using a Date object. PyEphem always deals with UTC, i.e., Coordinated Universal Time, and never the time in any particular time zone. A Date object, though, can be converted into the local time zone using the method ephem.localtime(), which takes a Date object as input. An example is show later in this document.
A Date object can be created using the method ephem.Date(). This method takes an argument that can be of any of the following types:
- A string of the format yyyy/mm/dd hh:mm:ss.ss, where all but the
year part is optional; BC dates: 1 BC is -1, 344 BC is -344.
- A floating point number corresponding to the number of days,
including fractional part, that has elapsed since 1899/12/31
12:00:00 UTC; negative numbers for dates older than this.
- A python datetime object.
- Any format that is accepted by the python datetime.datetime()
A few examples follow.
>>> a_date = ephem.Date("2010") # Set to 2010/1/1 00:00:00 UTC >>> a_date = ephem.Date("2010/03/27.5)" # 2010/3/27 12:00:00 UTC >>> a_date = ephem.Date("2010/3/27 12.56") # 2010/3/27 12:33:35 UTC >>> a_date = ephem.Date(-365) # 1898/12/31 12:00:00 UTC >>> a_date = ephem.Date((2010,1,10.5)) # 2010/1/10 12:00:00 UTC >>> # 2010/1/10 12:20:30.56 UTC >>> a_date = ephem.Date((2010,1,10,12,20,30.56)) >>> print ephem.Date(0.5) 1900/1/1 00:00:00
Since dates are stored as floating point numbers, arithmetic operations can be performed on Date objects, the results of which are floating point numbers. These can then be converted into Date objects.
>>> # Get the current UTC date, and time. >>> a_date = ephem.now() >>> a_date.real # the number stored internally 40511.953101851854 >>> # Move forward by two days >>> print a_date + 2 40513.9531019 >>> # Move back 2 days and 3 hours >>> print a_date - (2 + 3.0/24.0) 40509.8281019 >>> # Create a date object corresponding to 2 days and 3 hours before >>> # the current time. >>> b = ephem.Date(a_date - (2 + 3.0/24.0)) >>> print a_date 2010/12/1 10:52:28 >>> print b 2010/11/29 07:52:28
PyEphem provides three constants that represent an hour, a minute and a second, as a fraction of a day. These are ephem.hour, ephem.minute and ephem.second, respectively. So, for example, to add one hour to a particular date we can use the expression date + ephem.hour.
>>> a_date = ephem.now() # current UTC date and time >>> print a_date # note that fractional part is missing from seconds 2010/12/1 14:03:59 >>> # Add 2 days, 3 hours and 5 minutes >>> print ephem.Date(a_date + 2 + 3*ephem.hour + 5*ephem.minute) 2010/12/3 17:08:59
PyEphem also provides the constants ephem.J2000, ephem.B1950 and ephem.B1900, that store the day number of the corresponding epochs.
>>> print ephem.Date(ephem.J2000) 2000/1/1 12:00:00 >>> print ephem.Date(ephem.B1950) # This is NOT 1950 but Besselian 1950 1949/12/31 22:09:50
Information in Date objects can be extracted in several ways, in addition to obtaining it as a string using str() function, or as a floating point number in arithmetic calculations. The method tuple() gives a 6 element tuple containing the year, month, day, hour, minute and second, including the fractional part, of a Date object, in that order. The method triple() gives a three element tuple containing the year, month and day, including fractional part, of a Date object. The former can be passed to datetime.datetime() to create a python datetime.datetime object.
>>> a_date = ephem.now() >>> (year, month, day) = a_date.triple() >>> print year, month, day 2010 12 1.58611111111 >>> (y, mn, d, h, min, s) = a_date.tuple() >>> print y, mn, d, h, min, s 2010 12 1 14 3 59.9999997346
The local time corresponding to a Date object can be obtained using the ephem.localtime() method. This method returns a python datetime.datetime object, with date and time provided as argument converted into the local time zone. Note that time zone properties such as, tzname, tzinfo etc., of this object are not set.
>>> a_date = ephem.now() >>> print ephem.localtime(a_date) datetime.datetime(2010, 12, 1, 19, 33, 59, 2)
PyEphem can also calculate the Julian Day number / Julian date for a given date and time. PyEphem uses the Proleptic Julian Day.
>>> # Current Julian date. >>> jd = ephem.julian_date(ephem.now()) >>> print jd 2455532.09928
When specifying BC dates prefix the year with a negative sign. So, 1 BC is -1, 344 BC is -344 and so on. This is different from the usual practice in astronomy, where BC dates are represented by subtracting 1 from the year and prefixing a negative sign.
Positions in the sky are represented using the classes, ephem.Equatorial, ephem.Galactic and ephem.Ecliptic. Each of these represent a sky position in the appropriate coordinate system. We create instances of these by passing string representation of right ascension or longitude and declination or latitude, along with an optional epoch, to the above class constructors. The default epoch is 2000.0.
A position is defined by creating an object, of the type of the coordinate system required, by passing to the appropriate class the coordinates and an optional epoch. Positions can also be created from Body objects as described in the Body section.
Equatorial positions have the properties, ra and dec while the galactic and ecliptic positions have the properties, lat and long. Both galactic and ecliptic longitudes must be specified in degrees and not hours.
>>> # RA and DEC for equatorial coordinates; epoch = 2000.0 >>> eq = ephem.Equatorial("00:42:44.31","+41:16:09.4") >>> print eq.ra, eq.dec, eq.epoch # All three are needed to identify position (0:42:44.31, 41:16:09.4, 2000/1/1 00:00:00) >>> # Epoch 2010.5 i.e., the same ra, dec but at epoch 2010.5 >>> # NOT same as the previous coordinate >>> eq = ephem.Equatorial("00:42:44.31","+41:16:09.4",epoch='2010.5') print(eq.ra, eq.dec, eq.epoch) (0:42:44.31, 41:16:09.4, 2010/7/2 12:00:00) >>> # Latitude and longitude for galactic coordinates; >>> ga = ephem.Galactic("0","90",epoch=ephem.J2000) >>> print ga.long, ga.lat 0:00:00.0 90:00:00.0 >>> # Latitude and longitude for ecliptic coordinates >>> ec = ephem.Ecliptic("0","90",epoch=ephem.J2000) >>> print ec.long, ec.lat 0:00:00.0 90:00:00.0
Coordinates can be converted from one system to another by passing an object of the first type, as input to the class of the second type. This returns an object of the latter type, with coordinates converted to the latter system.
>>> # Galactic pole; epoch 2000.0 >>> ga = ephem.Galactic("0","90",epoch=ephem.J2000) >>> # Find ra and dec of Galactic pole >>> eq = ephem.Equatorial(ga) >>> print eq.ra, eq.dec 12:51:26.28 27:07:41.7 >>> # or use the ``get()`` method >>> print eq.get() (12:51:26.28, 27:07:41.7) >>> # Convert equatorial system from one epoch to another >>> eq = ephem.Equatorial("00:42:44.31","+41:16:09.4") #J2000.0 >>> eq1 = ephem.Equatorial(eq,epoch=ephem.B1950) >>> eq2 = ephem.Equatorial(eq,epoch="2010.0") >>> print eq1.get() (12:49:00.03, 27:24:00.0) >>> print eq2.get() (12:51:55.50, 27:04:26.3)
Objects in the sky are represented by instances of the Body class. Each Body object stores various properties, in addition to its coordinates. The PyEphem documentation lists properties of various types of bodies.
For a given body, PyEphem calculates three types of equatorial positions. These are :
- Astrometric geocentric position represented by the properties
a_ra and a_dec.
- Apparent geocentric position represented by the properties
g_ra and g_dec
- Apparent topocentric position represented by the properties
ra and dec, and alt and az.
When a Body is initialized, it needs a date and time, an optional epoch, and positions and velocities at a reference date and time. For built-in bodies PyEphem stores the positions and velocities i.e., proper motion, at the reference time J2000 and in a particular reference coordinate system.
If a date and time is provided when initializing a Body, PyEphem will calculate where the body will be at the specified date and time, using the stored reference position and velocities. If an epoch is provided, then PyEphem will precess the coordinate system and calculates the position of the body, at the given date and time, in the precessed coordinate system. This is then stored as the Astrometric geocentric position. The epoch defaults to J2000 and hence no precession gets applied if epoch is not given. These values are available as a_ra and a_dec. This is the position that will appear in a star atlas for the given epoch.
To calculate where the object will "appear" for an observer at the center of the Earth, at the given date and time, PyEphem takes the position found, before applying precession, and then applies corrections for light-travel time (for solar system objects), gravitational deflection of light, nutation and aberration. Aberration effects are not considered for the Moon. This position refers to the orientation of the Earth’s axes at the given date and time and hence are referred to as "epoch-of-date" coordinates; this will not be same as star atlas position. This position is stored as Apparent geocentric position and are available as g_ra and g_dec.
If an Observer is passed to the constructor, instead of a date and time, then PyEphem will use the date and time, and epoch if none is provided, stored in the Observer to perform the above calculations. In addition, PyEphem will apply parallax corrections, to g_ra and g_dec, to get the apparent position of the body as seen by the Observer on the surface of the Earth. This position is called the Apparent topocentric position and are stored in the members ra and dec. If no Observer is given, i.e., only date and time, and optional epoch, then ra and dec are set to the same values as g_ra and g_dec, respectively. In addition to ra and dec, PyEphem also calculates the alt and az properties of the Body: ra and dec and alt and az both define where in the sky the object will appear for the given Observer at the Observer‘s date and time.
After creating a Body we can calculate/re-calculate these positions by calling the method compute() of a body. In fact PyEphem calls this method if appropriate information is available while creating the object. As mentioned above, if the information passed to compute is date and time, with optional epoch, then the properties a_ra, a_dec, g_ra and g_dec are set. If the information passed is an Observer, then, the above are set for the date and time of the observer, and ra, dec, alt and az are also set. If no explicit epoch is given then the epoch property of the Observer is used.
In the following sections we describe several Body types used in PyEphem. For all the bodies discussed below, see the Observer section for setting topocentric properties using an Observer.
>>> sun = ephem.Sun() >>> # Position are not set. >>> print sun.a_ra, sun.a_dec # raises RuntimeError Traceback (most recent call last): File "<stdin>", line 1, in <module> RuntimeError: field astrora undefined until first compute() >>> # a_ra, a_dec, g_ra, g_dec are set to values for time 2010/2/1 >>> # epoch 1900 >>> sun.compute("2010/2/1",epoch=1900) >>> print sun.a_ra, sun.a_dec 20:51:58.70 -17:35:39.8 >>> print sun.g_ra, sun.g_dec 20:57:52.72 -17:11:29.5 >>> # Apparent topocentric positions are not set >>> print sun.ra , sun.dec # same as sun.g_ra and sun.g_dec since no 20:57:52.72 -17:11:29.5 # Observer is specified. >>> print sun.alt, sun.az # raises RuntimeError Traceback (most recent call last): File "<stdin>", line 1, in <module> RuntimeError: field alt undefined because the most recent compute() was supplied a date rather than an Observer >>> moon = ephem.Moon(ephem.now()) # current time, epoch 2000.0 >>> # Other properties of Moon. Earth distance is in AU. >>> print moon.name, moon.earth_distance, moon.a_epoch Moon 0.00247154687531 2000/1/1 12:00:00
PyEphem provides the following planets and moons:
ephem.Io, ephem.Europa, ephem.Ganymede, ephem.Callisto
ephem.Mimas, ephem.Enceladus, ephem.Tethys, ephem.Dione, ephem.Rhea, ephem.Titan, ephem.Iapetus, ephem.Hyperion
ephem.Miranda, ephem.Ariel, ephem.Umbriel, ephem.Titania, ephem.Oberon
In addition to the properties described above, planets and moons have many other properties. Solar system bodies have hlong, the Heliocentric longitude, hlat, the Heliocentric latitude, phase, percentage of surface illuminated etc., . Saturn bodies have information on their rings. See PyEphem documentation on Bodies for more details.
PyEphem comes with a collection of named stars. These can be created by passing the name of the star to the method ephem.star().
>>> altair = ephem.star("Altair") >>> print altair.name Altair
The following code lists all named stars available in PyEphem:
>>> import ephem.stars >>> for star in ephem.stars.db.split("\n"): ... print star.split(",") Sirrah Caph Algenib Schedar Mirach Achernar Almach Hamal Polaris . . . Fomalhaut Scheat Markab
Bodies can also be defined by specifying their properties in the XEphem database format or the TLE format and reading it in using the ephem.readdb() and ephem.readtle() methods, respectively. This is needed for solar system bodies such as asteroids, comets, and Earth satellites; their orbits need to be frequently updated using latest orbital elements. A leading underscore is used to distinguish properties read from database from those set by the compute() method; for example, _ra and _dec.
For more information on the XEphem database format, see section 126.96.36.199 of the XEphem manual. The NORAD Two-Line Element (TLE) format is explained in section 7.14 of the XEphem manual. Also see the NASA spaceflight FAQ page.
See, sections titled "catalog formats" and "bodies with orbital elements" in the PyEphem quick-start page, for more on catalog formats. This page also has documentation listing properties of various types of objects.
An observer in PyEphem is identified by, date, the date of observation, epoch, the epoch of observation, lat, the geographic latitude, long, the geographic longitude, elevation, the elevation in meters, temp, the temperature in degree Celsius and pressure, the atmospheric pressure in milli-bar.
The property epoch is the epoch that will be used to calculate astrometric coordinates, when the position of an object with respect to an observer is to be calculated. This defaults to 2000.0. Date defaults to ephem.now(), pressure to 1010 milli-bars and temperature to 10 degree Celsius. If pressure is non-zero then correction for atmospheric refraction is included while calculating ra, dec, alt and az. Set pressure = 0, to prevent the application of this correction.
We can create an observer and set various properties as shown below.
>>> kpno = ephem.Observer() >>> # The following data are from the Skycalc program >>> # by J. R. Thorstensen >>> # longitude is 7.44111 hours or 111.61665 degrees west >>> # same as kpno.long = "-111.61665" (degrees) >>> # positive values are used for eastern longitudes. >>> kpno.long = ephem.hours("-7.44111") >>> # positive values are for northern latitudes. >>> kpno.lat = ephem.degrees("31.9533") >>> kpno.elevation = 1925.0 + 700.0 # meters >>> kpno.date = "2010/1/1"
PyEphem has a number of built-in cities, that can be used to create Observer instances for these places.
>>> london = ephem.city("London") >>> london.date = "2010/1/1" >>> print london.sidereal_time() 6:41:42.04
The following code will list all the cities defined in PyEphem.
>>> from ephem import cities >>> for i in cities._city_data.keys(): ... print i Kiev Lille Paris Oslo Rio de Janeiro Johannesburg San Francisco Sao Paulo Santiago Abu Dhabi . . . Colombo Lisbon Amsterdam Copenhagen Warsaw
To calculate the position of an object as seen from the location of the observer, we pass the Observer instance to the Body instance. The epoch property of the Observer is used in calculating astrometric positions; default is 2000.0.
>>> kpno = ephem.Observer() >>> kpno.long = ephem.hours("-7.44111") >>> kpno.lat = ephem.degrees("31.9533") >>> kpno.elevation = 1925.0 + 700.0 >>> kpno.date = "2010/1/1" >>> jupiter = ephem.Jupiter(kpno) >>> print jupiter.alt, jupiter.az 40:35:18.8 205:57:25.3 >>> print jupiter.ra, jupiter.dec 21:55:41.10 -13:36:33.5
An Observer instance has several methods to find the time of occurrence of various events at its location.
>>> # Set date at kpno to the current date and time >>> kpno.date = ephem.now() >>> # All dates and time are in UTC; KPNO is UTC - 7 >>> print ephem.Date(kpno.date - 7*ephem.hour) # local time 2010/12/1 23:29:45 >>> # We want to find the rising and setting time for Sun at KPNO >>> sun = ephem.Sun() >>> print kpno.previous_rising(sun) # UTC time 2010/12/1 14:08:44 >>> print ephem.Date(kpno.previous_rising(sun) - 7*ephem.hour) # local time 2010/12/1 07:08:44 >>> print kpno.next_setting(sun) # UTC time 2010/12/3 00:22:09 >>> print ephem.Date(kpno.next_setting(sun) - 7*ephem.hour) # local time 2010/12/2 17:22:09
PyEphem uses the horizon property of an Observer while calculating the setting and rising time of an Object from the Observer‘s location. The horizon is the altitude of the upper limb of a body at the moment it is considered to be rising and setting. So, for example, we can set it to a value greater than 0 if we want to define rising as the event when a body is high above obstructions.
To get values that are consistent with the United States Naval Observatory definition of rising and setting, we should set pressure to 0 and then set horizon to -34 arcminutes.
The following is an example taken directly from the PyEphem documentation.
>>> sun = ephem.Sun() >>> greenwich = ephem.Observer() >>> greenwich.lat = '51:28:38' >>> print greenwich.horizon 0:00:00.0 >>> greenwich.date = '2007/10/1' >>> r1 = greenwich.next_rising(sun) >>> # Using USNO definition >>> greenwich.pressure = 0 >>> greenwich.horizon = '-0:34' >>> greenwich.date = '2007/10/1' >>> r2 = greenwich.next_rising(sun) >>> print 'Visual sunrise:', r1 Visual sunrise: 2007/10/1 05:59:29 >>> print 'Naval Observatory sunrise:', Naval Observatory sunrise: 2007/10/1 05:59:50
- ephem.c = 299792458
Exact speed of light in meters/second.
- ephem.earth_radius = 6378160
Earth’s equatorial radius in meters
- ephem.moon_radius = 1740000
Moon’s equatorial radius in meters
- ephem.meters_per_au = 1.4959787e11
1 astronomical unit in meters.
Gives the difference between Terrestrial Time and Coordinated Universal Time (UTC), in seconds.
Given an argument that can be converted into an ephem.Date object, this function returns the Julian date corresponding to the given date and time. If no input is given, then the Julian date of the current time is returned.
This method converts the given PyEphem date into local time and returns a python datetime object representing this value.
These methods take as input a date and returns the date of the relevant event closest to the input date.
The following give dates of relevant equinox:
The following give dates of relevant solstice:
These methods take as input a date and returns the date of the relevant event closest to the input date.
This method takes as input an object or coordinate, and returns the constellation in which the given object or coordinate lie.
Returns the angular separation between two objects or positions in degrees.
ephem.uranometria and ephem.uranometria2000
Given right ascension and declination in radians, these methods return the page of the Uranometria and Uranometria 2000 catalogs, repectively, displaying that location.
Given a one-dimensional function, a lower limit, an upper limit and an optional precision, returns the value of the independent variable at which the function reaches zero, to within the given precision. The default precision is a half-second of clock time.
Which constellation is Mars in?
>>> import ephem as ep >>> mars = ep.Mars(ep.now()) >>> print ep.constellation(mars) ('Oph', 'Ophiuchus')
Is Mars in the sky above London right now? What about Saturn?
>>> import ephem as ep >>> # create an observer object for London >>> london = ep.city("London") # London is built-in. >>> london.date = ep.now() >>> london.pressure = 0 # Let's neglect atmospheric refraction >>> # Create a Mars object and a Saturn object, for the London observer. >>> mars = ep.Mars(london) >>> saturn = ep.Saturn(london) >>> # print altitude >>> print "Mars:", mars.alt Mars: -19:00:30.9 >>> print "Saturn:", saturn.alt Saturn: 32:20:45.1
The angular separation between Mars and Saturn?
>>> # continuing from the above example... >>> print ep.separation(mars,saturn) # result in degrees. 71:04:39.1