Basic astronomical calculations (`pwkit.astutil`)¶

This module provides functions and constants for doing a variety of basic calculations and conversions that come up in astronomy.

This topics covered in this module are:

Useful Constants

Sexagesimal Notation

Working with Angles

Simple Operations on 2D Gaussians

Basic Astrometry

Miscellaneous Astronomical Computations

Angles are always measured in radians, whereas some other astronomical codebases prefer degrees.

Useful Constants¶

The following useful constants are provided:

pi

Mathematical π.

twopi

Mathematical 2π.

halfpi

Mathematical π/2.

R2A

A constant for converting radians to arcseconds by multiplication:

arcsec = radians * astutil.R2A

Equal to 3600 * 180 / pi or about 206265.

A2R

A constant for converting arcseconds to radians by multiplication:

radians = arcsec * astutil.A2R

R2D

Analogous to R2A: a constant for converting radians to degrees

D2R

Analogous to A2R: a constant for converting degrees to radians

R2H

Analogous to R2A: a constant for converting radians to hours

H2R

Analogous to A2R: a constant for converting hours to radians

F2S

A constant for converting a Gaussian FWHM (full width at half maximum) to a standard deviation (σ) value by multiplication:

sigma = fwhm * astutil.F2S

Equal to (8 * ln(2))**-0.5 or about 0.425.

S2F

A constant for converting a Gaussian standard deviation (σ) value to a FWHM (full width at half maximum) by multiplication.

J2000

The astronomical J2000.0 epoch as a MJD (modified Julian Date). Precisely equal to 51544.5.

Sexagesimal Notation¶

`fmthours`(radians[, norm, precision, seps])	Format an angle as sexagesimal hours in a string.
`fmtdeglon`(radians[, norm, precision, seps])	Format a longitudinal angle as sexagesimal degrees in a string.
`fmtdeglat`(radians[, norm, precision, seps])	Format a latitudinal angle as sexagesimal degrees in a string.
`fmtradec`(rarad, decrad[, precision, raseps, ...])	Format equatorial coordinates in a single sexagesimal string.
`parsehours`(hrstr)	Parse a string formatted as sexagesimal hours into an angle.
`parsedeglat`(latstr)	Parse a latitude formatted as sexagesimal degrees into an angle.
`parsedeglon`(lonstr)	Parse a longitude formatted as sexagesimal degrees into an angle.

pwkit.astutil.fmthours(radians, norm='wrap', precision=3, seps='::')[source]¶

Format an angle as sexagesimal hours in a string.

Arguments are:

radians: The angle, in radians.
norm (default “wrap”): The normalization mode, used for angles outside of the standard range of 0 to 2π. If “none”, the value is formatted ignoring any potential problems. If “wrap”, it is wrapped to lie within the standard range. If “raise”, a ValueError is raised.
precision (default 3): The number of decimal places in the “seconds” place to use in the formatted string.
seps (default “::”): A two- or three-item iterable, used to separate the hours, minutes, and seconds components. If a third element is present, it appears after the seconds component. Specifying “hms” yields something like “12h34m56s”; specifying ['', ''] yields something like “123456”.

Returns a string.

pwkit.astutil.fmtdeglon(radians, norm='wrap', precision=2, seps='::')[source]¶

Format a longitudinal angle as sexagesimal degrees in a string.

Arguments are:

radians: The angle, in radians.
norm (default “wrap”): The normalization mode, used for angles outside of the standard range of 0 to 2π. If “none”, the value is formatted ignoring any potential problems. If “wrap”, it is wrapped to lie within the standard range. If “raise”, a ValueError is raised.
precision (default 2): The number of decimal places in the “arcseconds” place to use in the formatted string.
seps (default “::”): A two- or three-item iterable, used to separate the degrees, arcminutes, and arcseconds components. If a third element is present, it appears after the arcseconds component. Specifying “dms” yields something like “12d34m56s”; specifying ['', ''] yields something like “123456”.

Returns a string.

pwkit.astutil.fmtdeglat(radians, norm='raise', precision=2, seps='::')[source]¶

Format a latitudinal angle as sexagesimal degrees in a string.

Arguments are:

radians: The angle, in radians.
norm (default “raise”): The normalization mode, used for angles outside of the standard range of -π/2 to π/2. If “none”, the value is formatted ignoring any potential problems. If “wrap”, it is wrapped to lie within the standard range. If “raise”, a ValueError is raised.
precision (default 2): The number of decimal places in the “arcseconds” place to use in the formatted string.
seps (default “::”): A two- or three-item iterable, used to separate the degrees, arcminutes, and arcseconds components. If a third element is present, it appears after the arcseconds component. Specifying “dms” yields something like “+12d34m56s”; specifying ['', ''] yields something like “123456”.

Returns a string. The return value always includes a plus or minus sign. Note that the default of norm is different than in fmthours() and fmtdeglon() since it’s not so clear what a “latitude” of 110 degrees (e.g.) means.

pwkit.astutil.fmtradec(rarad, decrad, precision=2, raseps='::', decseps='::', intersep=' ')[source]¶

Format equatorial coordinates in a single sexagesimal string.

Returns a string of the RA/lon coordinate, formatted as sexagesimal hours, then intersep, then the Dec/lat coordinate, formatted as degrees. This yields something like “12:34:56.78 -01:23:45.6”. Arguments are:

rarad: The right ascension coordinate, in radians. More generically, this is the longitudinal coordinate; note that the ordering in this function differs than the other spherical functions, which generally prefer coordinates in “lat, lon” order.
decrad: The declination coordinate, in radians. More generically, this is the latitudinal coordinate.
precision (default 2): The number of decimal places in the “arcseconds” place of the latitudinal (declination) coordinate. The longitudinal (right ascension) coordinate gets one additional place, since hours are bigger than degrees.
raseps (default “::”): A two- or three-item iterable, used to separate the hours, minutes, and seconds components of the RA/lon coordinate. If a third element is present, it appears after the seconds component. Specifying “hms” yields something like “12h34m56s”; specifying ['', ''] yields something like “123456”.
decseps (default “::”): A two- or three-item iterable, used to separate the degrees, arcminutes, and arcseconds components of the Dec/lat coordinate.
intersep (default ” “): The string separating the RA/lon and Dec/lat coordinates

pwkit.astutil.parsehours(hrstr)[source]¶

Parse a string formatted as sexagesimal hours into an angle.

This function converts a textual representation of an angle, measured in hours, into a floating point value measured in radians. The format of hrstr is very limited: it may not have leading or trailing whitespace, and the components of the sexagesimal representation must be separated by colons. The input must therefore resemble something like "12:34:56.78". A ValueError will be raised if the input does not resemble this template. Hours greater than 24 are not allowed, but negative values are.

pwkit.astutil.parsedeglat(latstr)[source]¶

Parse a latitude formatted as sexagesimal degrees into an angle.

This function converts a textual representation of a latitude, measured in degrees, into a floating point value measured in radians. The format of latstr is very limited: it may not have leading or trailing whitespace, and the components of the sexagesimal representation must be separated by colons. The input must therefore resemble something like "-00:12:34.5". A ValueError will be raised if the input does not resemble this template. Latitudes greater than 90 or less than -90 degrees are not allowed.

pwkit.astutil.parsedeglon(lonstr)[source]¶

Parse a longitude formatted as sexagesimal degrees into an angle.

This function converts a textual representation of a longitude, measured in degrees, into a floating point value measured in radians. The format of lonstr is very limited: it may not have leading or trailing whitespace, and the components of the sexagesimal representation must be separated by colons. The input must therefore resemble something like "270:12:34.5". A ValueError will be raised if the input does not resemble this template. Values of any sign and magnitude are allowed, and they are not normalized (e.g. to lie within the range [0, 2π]).

Working with Angles¶

`angcen`(a)
`orientcen`(a)
`sphdist`(lat1, lon1, lat2, lon2)	Calculate the distance between two locations on a sphere.
`sphbear`(lat1, lon1, lat2, lon2[, tol])	Calculate the bearing between two locations on a sphere.
`sphofs`(lat1, lon1, r, pa[, tol, rmax])	Offset from one location on the sphere to another.
`parang`(hourangle, declination, latitude)	Calculate the parallactic angle of a sky position.

pwkit.astutil.angcen(a)¶

“Center” an angle a to be between -π and +π.

This is done by adding or subtracting multiples of 2π as necessary. Both a and the return value are in radians. The argument may be a vector.

pwkit.astutil.orientcen(a)¶

“Center” an orientation a to be between -π/2 and +π/2.

This is done by adding or subtract multiples of π as necessary. Both a and the return value are in radians. The argument may be a vector.

An “orientation” is different than an angle because values that differ by just π, not 2π, are considered equivalent. Orientations can come up in the discussion of linear polarization, for example.

pwkit.astutil.sphdist(lat1, lon1, lat2, lon2)[source]¶

Calculate the distance between two locations on a sphere.

lat1: The latitude of the first location.
lon1: The longitude of the first location.
lat2: The latitude of the second location.
lon2: The longitude of the second location.

Returns the separation in radians. All arguments are in radians as well. The arguments may be vectors.

Note that the ordering of the arguments maps to the nonstandard ordering (Dec, RA) in equatorial coordinates. In a spherical projection it maps to (Y, X) which may also be unexpected.

The distance is computed with the “specialized Vincenty formula”. Faster but more error-prone formulae are possible; see Wikipedia on Great-circle Distance.

pwkit.astutil.sphbear(lat1, lon1, lat2, lon2, tol=1e-15)[source]¶

Calculate the bearing between two locations on a sphere.

lat1: The latitude of the first location.
lon1: The longitude of the first location.
lat2: The latitude of the second location.
lon2: The longitude of the second location.
tol: Tolerance for checking proximity to poles and rounding to zero.

The bearing (AKA the position angle, PA) is the orientation of point 2 with regards to point 1 relative to the longitudinal axis. Returns the bearing in radians. All arguments are in radians as well. The arguments may be vectors.

Note that the ordering of the arguments maps to the nonstandard ordering (Dec, RA) in equatorial coordinates. In a spherical projection it maps to (Y, X) which may also be unexpected.

The sign convention is astronomical: bearings range from -π to π, with negative values if point 2 is in the western hemisphere with regards to point 1, positive if it is in the eastern. (That is, “east from north”.) If point 1 is very near the pole, the bearing is undefined and the result is NaN.

The tol argument is used for checking proximity to the poles and for rounding the bearing to precisely zero if it’s extremely small.

Derived from bear() in angles.py from Prasanth Nair. His version is BSD licensed. This one is sufficiently different that I think it counts as a separate implementation.

pwkit.astutil.sphofs(lat1, lon1, r, pa, tol=0.01, rmax=None)[source]¶

Offset from one location on the sphere to another.

This function is given a start location, expressed as a latitude and longitude, a distance to offset, and a direction to offset (expressed as a bearing, AKA position angle). It uses these to compute a final location. This function mirrors sphdist() and sphbear() such that:

# If:
r = sphdist (lat1, lon1, lat2a, lon2a)
pa = sphbear (lat1, lon1, lat2a, lon2a)
lat2b, lon2b = sphofs (lat1, lon1, r, pa)
# Then lat2b = lat2a and lon2b = lon2a

Arguments are:

lat1: The latitude of the start location.
lon1: The longitude of the start location.
r: The distance to offset by.
pa: The position angle (“PA” or bearing) to offset towards.
tol: The tolerance for the accuracy of the calculation.
rmax: The maximum allowed offset distance.

Returns a pair (lat2, lon2). All arguments and the return values are measured in radians. The arguments may be vectors. The PA sign convention is astronomical, measuring orientation east from north.

Note that the ordering of the arguments and return values maps to the nonstandard ordering (Dec, RA) in equatorial coordinates. In a spherical projection it maps to (Y, X) which may also be unexpected.

The offset is computed naively as:

lat2 = lat1 + r * cos (pa)
lon2 = lon1 + r * sin (pa) / cos (lat2)

This will fail for large offsets. Error checking can be done in two ways. If tol is not None, sphdist() is used to calculate the actual distance between the two locations, and if the magnitude of the fractional difference between that and r is larger than tol, ValueError is raised. This will add an overhead to the computation that may be significant if you’re going to be calling this function a lot.

Additionally, if rmax is not None, magnitudes of r greater than rmax are rejected. For reference, an r of 0.2 (~11 deg) gives a maximum fractional distance error of ~3%.

pwkit.astutil.parang(hourangle, declination, latitude)[source]¶

Calculate the parallactic angle of a sky position.

This computes the parallactic angle of a sky position expressed in terms of an hour angle and declination. Arguments:

hourangle: The hour angle of the location on the sky.
declination: The declination of the location on the sky.
latitude: The latitude of the observatory.

Inputs and outputs are all in radians. Implementation adapted from GBTIDL parangle.pro.

Simple Operations on 2D Gaussians¶

pwkit.astutil.gaussian_convolve(maj1, min1, pa1, maj2, min2, pa2)[source]¶

Convolve two Gaussians analytically.

Given the shapes of two 2-dimensional Gaussians, this function returns the shape of their convolution.

Arguments:

maj1: Major axis of input Gaussian 1.
min1: Minor axis of input Gaussian 1.
pa1: Orientation angle of input Gaussian 1, in radians.
maj2: Major axis of input Gaussian 2.
min2: Minor axis of input Gaussian 2.
pa2: Orientation angle of input Gaussian 2, in radians.

The return value is (maj3, min3, pa3), with the same format as the input arguments. The axes can be measured in any units, so long as they’re consistent.

Implementation copied from MIRIAD’s gaufac.

pwkit.astutil.gaussian_deconvolve(smaj, smin, spa, bmaj, bmin, bpa)[source]¶

Deconvolve two Gaussians analytically.

Given the shapes of 2-dimensional “source” and “beam” Gaussians, this returns a deconvolved “result” Gaussian such that the convolution of “beam” and “result” is “source”.

Arguments:

smaj: Major axis of source Gaussian.
smin: Minor axis of source Gaussian.
spa: Orientation angle of source Gaussian, in radians.
bmaj: Major axis of beam Gaussian.
bmin: Minor axis of beam Gaussian.
bpa: Orientation angle of beam Gaussian, in radians.

The return value is (rmaj, rmin, rpa, status). The first three values have the same format as the input arguments. The status result is one of “ok”, “pointlike”, or “fail”. A “pointlike” status indicates that the source and beam shapes are difficult to distinguish; a “fail” status indicates that the two shapes seem to be mutually incompatible (e.g., source and beam are very narrow and orthogonal).

The axes can be measured in any units, so long as they’re consistent.

Ideally if:

rmaj, rmin, rpa, status = gaussian_deconvolve (smaj, smin, spa, bmaj, bmin, bpa)

then:

smaj, smin, spa = gaussian_convolve (rmaj, rmin, rpa, bmaj, bmin, bpa)

Implementation derived from MIRIAD’s gaudfac. This function currently doesn’t do a great job of dealing with pointlike sources, i.e. ones where “source” and “beam” are nearly indistinguishable.

Basic Astrometry¶

The AstrometryInfo class can be used to perform basic astrometric calculations that are nonetheless fairly accurate.

class pwkit.astutil.AstrometryInfo(simbadident=None, **kwargs)[source]¶

Holds astrometric data and their uncertainties, and can predict positions with uncertainties.

The attributes encoding the astrometric data are as follows. Values of None will be treated as unknown. Most of this information can be automatically filled in from the fill_from_simbad() function, if you trust Simbad.

`ra`	The J2000 right ascension of the object, measured in radians.
`dec`	The J2000 declination of the object, measured in radians.
`pos_u_maj`	Major axis of the error ellipse for the object position, in radians.
`pos_u_min`	Minor axis of the error ellipse for the object position, in radians.
`pos_u_pa`	Position angle (really orientation) of the error ellipse for the object position, east from north, in radians.
`pos_epoch`	The epoch of position, that is, the date when the position was measured, in MJD[TT].
`promo_ra`	The proper motion in right ascension, in milliarcsec per year.
`promo_dec`	The object's proper motion in declination, in milliarcsec per year.
`promo_u_maj`	Major axis of the error ellipse for the object's proper motion, in milliarcsec per year.
`promo_u_min`	Minor axis of the error ellipse for the object's proper motion, in milliarcsec per year.
`promo_u_pa`	Position angle (really orientation) of the error ellipse for the object proper motion, east from north, in radians.
`parallax`	The object's parallax, in milliarcsec.
`u_parallax`	Uncertainty in the object's parallax, in milliarcsec.
`vradial`	The object's radial velocity, in km/s.
`u_vradial`	The uncertainty in the object's radial velocity, in km/s.

Methods are:

`verify`([complain])	Validate that the attributes are self-consistent.
`predict`(mjd[, complain, n])	Predict the object position at a given MJD.
`print_prediction`(ptup[, precision])	Print a summary of a predicted position.
`predict_without_uncertainties`(mjd[, complain])	Predict the object position at a given MJD.
`fill_from_simbad`(ident[, debug])	Fill in astrometric information using the Simbad web service.
`fill_from_allwise`(ident[, catalog_ident])	Fill in astrometric information from the AllWISE catalog using Astroquery.

The stringification of an AstrometryInfo class formats its fields in a human-readable, multiline format that uses Unicode characters.

AstrometryInfo.ra = None¶: The J2000 right ascension of the object, measured in radians.

AstrometryInfo.dec = None¶: The J2000 declination of the object, measured in radians.

AstrometryInfo.pos_u_maj = None¶: Major axis of the error ellipse for the object position, in radians.

AstrometryInfo.pos_u_min = None¶: Minor axis of the error ellipse for the object position, in radians.

AstrometryInfo.pos_u_pa = None¶: Position angle (really orientation) of the error ellipse for the object position, east from north, in radians.

AstrometryInfo.pos_epoch = None¶: The epoch of position, that is, the date when the position was measured, in MJD[TT].

AstrometryInfo.promo_ra = None¶: The proper motion in right ascension, in milliarcsec per year. XXX: cos(dec) confusion!

AstrometryInfo.promo_dec = None¶: The object’s proper motion in declination, in milliarcsec per year.

AstrometryInfo.promo_u_maj = None¶: Major axis of the error ellipse for the object’s proper motion, in milliarcsec per year.

AstrometryInfo.promo_u_min = None¶: Minor axis of the error ellipse for the object’s proper motion, in milliarcsec per year.

AstrometryInfo.promo_u_pa = None¶: Position angle (really orientation) of the error ellipse for the object proper motion, east from north, in radians.

AstrometryInfo.parallax = None¶: The object’s parallax, in milliarcsec.

AstrometryInfo.u_parallax = None¶: Uncertainty in the object’s parallax, in milliarcsec.

AstrometryInfo.vradial = None¶: The object’s radial velocity, in km/s. NOTE: not relevant in our usage.

AstrometryInfo.u_vradial = None¶: The uncertainty in the object’s radial velocity, in km/s. NOTE: not relevant in our usage.

AstrometryInfo.verify(complain=True)[source]¶

Validate that the attributes are self-consistent.

This function does some basic checks of the object attributes to ensure that astrometric calculations can legally be performed. If the complain keyword is true, messages may be printed to sys.stderr if non-fatal issues are encountered.

Returns self.

AstrometryInfo.predict(mjd, complain=True, n=20000)[source]¶

Predict the object position at a given MJD.

The return value is a tuple (ra, dec, major, minor, pa), all in radians. These are the predicted position of the object and its uncertainty at mjd. If complain is True, print out warnings for incomplete information. n is the number of Monte Carlo samples to draw for computing the positional uncertainty.

The uncertainty ellipse parameters are sigmas, not FWHM. These may be converted with the S2F constant.

This function relies on the external skyfield package.

AstrometryInfo.print_prediction(ptup, precision=2)[source]¶

Print a summary of a predicted position.

The argument ptup is a tuple returned by predict(). It is printed to sys.stdout in a reasonable format that uses Unicode characters.

AstrometryInfo.predict_without_uncertainties(mjd, complain=True)[source]¶

Predict the object position at a given MJD.

The return value is a tuple (ra, dec), in radians, giving the predicted position of the object at mjd. Unlike predict(), the astrometric uncertainties are ignored. This function is therefore deterministic but potentially misleading.

If complain is True, print out warnings for incomplete information.

This function relies on the external skyfield package.

AstrometryInfo.fill_from_simbad(ident, debug=False)[source]¶

Fill in astrometric information using the Simbad web service.

This uses the CDS Simbad web service to look up astrometric information for the source name ident and fills in attributes appropriately. Values from Simbad are not always reliable.

Returns self.

AstrometryInfo.fill_from_allwise(ident, catalog_ident='II/328/allwise')[source]¶

Fill in astrometric information from the AllWISE catalog using Astroquery.

This uses the astroquery module to query the AllWISE (2013wise.rept….1C) source catalog through the Vizier (2000A&AS..143…23O) web service. It then fills in the instance with the relevant information. Arguments are:

ident: The AllWISE catalog identifier of the form "J112254.70+255021.9".
catalog_ident: The Vizier designation of the catalog to query. The default is “II/328/allwise”, the current version of the AllWISE catalog.

Raises PKError if something unexpected happens that doesn’t itself result in an exception within astroquery.

You should probably prefer fill_from_simbad() for objects that are known to the CDS Simbad service, but not all objects in the AllWISE catalog are so known.

If you use this function, you should acknowledge AllWISE and Vizier.

Returns self.

A few helper functions may also be of interest:

`load_skyfield_data`()	Load data files used in Skyfield.
`get_2mass_epoch`(tmra, tmdec[, debug])	Given a 2MASS position, look up the epoch when it was observed.
`get_simbad_astrometry_info`(ident[, items, debug])	Fetch astrometric information from the Simbad web service.

pwkit.astutil.load_skyfield_data()[source]¶

Load data files used in Skyfield. This will download files from the internet if they haven’t been downloaded before.

Skyfield downloads files to the current directory by default, which is not ideal. Here we abuse astropy and use its cache directory to cache the data files per-user. If we start downloading files in other places in pwkit we should maybe make this system more generic. And the dep on astropy is not at all necessary.

Skyfield will print out a progress bar as it downloads things.

Returns (planets, ts), the standard Skyfield ephemeris and timescale data files.

pwkit.astutil.get_2mass_epoch(tmra, tmdec, debug=False)[source]¶

Given a 2MASS position, look up the epoch when it was observed.

This function uses the CDS Vizier web service to look up information in the 2MASS point source database. Arguments are:

tmra: The source’s J2000 right ascension, in radians.
tmdec: The source’s J2000 declination, in radians.
debug: If True, the web server’s response will be printed to sys.stdout.

The return value is an MJD. If the lookup fails, a message will be printed to sys.stderr (unconditionally!) and the J2000 epoch will be returned.

pwkit.astutil.get_simbad_astrometry_info(ident, items=['COO(d;A)', 'COO(d;D)', 'COO(E)', 'COO(B)', 'PM(A)', 'PM(D)', 'PM(E)', 'PLX(V)', 'PLX(E)', 'RV(V)', 'RV(E)'], debug=False)[source]¶

Fetch astrometric information from the Simbad web service.

Given the name of a source as known to the CDS Simbad service, this function looks up its positional information and returns it in a dictionary. In most cases you should use an AstrometryInfo object and its fill_from_simbad() method instead of this function.

Arguments:

ident: The Simbad name of the source to look up.
items: An iterable of data items to look up. The default fetches position, proper motion, parallax, and radial velocity information. Each item name resembles the string COO(d;A) or PLX(E). The allowed formats are defined on this CDS page.
debug: If true, the response from the webserver will be printed.

The return value is a dictionary with a key corresponding to the textual result returned for each requested item.

Miscellaneous Astronomical Computations¶

These functions don’t fit under the other rubrics very well.

`abs2app`(abs_mag, dist_pc)	Convert an absolute magnitude to an apparent magnitude, given a source's (luminosity) distance in parsecs.
`app2abs`(app_mag, dist_pc)	Convert an apparent magnitude to an absolute magnitude, given a source's (luminosity) distance in parsecs.

pwkit.astutil.abs2app(abs_mag, dist_pc)[source]¶

Convert an absolute magnitude to an apparent magnitude, given a source’s (luminosity) distance in parsecs.

Arguments:

abs_mag: Absolute magnitude.
dist_pc: Distance, in parsecs.

Returns the apparent magnitude. The arguments may be vectors.

pwkit.astutil.app2abs(app_mag, dist_pc)[source]¶

Convert an apparent magnitude to an absolute magnitude, given a source’s (luminosity) distance in parsecs.

Arguments:

app_mag: Apparent magnitude.
dist_pc: Distance, in parsecs.

Returns the absolute magnitude. The arguments may be vectors.

Basic astronomical calculations (pwkit.astutil)¶

Useful Constants¶

Sexagesimal Notation¶

Working with Angles¶

Simple Operations on 2D Gaussians¶

Basic Astrometry¶

Miscellaneous Astronomical Computations¶

Basic astronomical calculations (`pwkit.astutil`)¶