#
# LOFAR Transients Key Project
"""
General purpose astronomical coordinate handling routines.
"""
import sys
import math
import pywcs
import logging
import datetime
import pytz
from casacore.measures import measures
from casacore.quanta import quantity
logger = logging.getLogger(__name__)
# Note that we take a +ve longitude as WEST.
CORE_LAT = 52.9088
CORE_LON = -6.8689
# ITRF position of CS002
# Should be a good approximation for anything refering to the LOFAR core.
ITRF_X = 3826577.066110000
ITRF_Y = 461022.947639000
ITRF_Z = 5064892.786
# Useful constants
SECONDS_IN_HOUR = 60**2
SECONDS_IN_DAY = 24 * SECONDS_IN_HOUR
[docs]def julian_date(time=None, modified=False):
"""
Calculate the Julian date at a given timestamp.
Args:
time (datetime.datetime): Timestamp to calculate JD for.
modified (bool): If True, return the Modified Julian Date:
the number of days (including fractions) which have elapsed between
the start of 17 November 1858 AD and the specified time.
Returns:
float: Julian date value.
"""
if not time:
time = datetime.datetime.now(pytz.utc)
mjdstart = datetime.datetime(1858, 11, 17, tzinfo=pytz.utc)
mjd = time - mjdstart
mjd_daynumber = (mjd.days + mjd.seconds / (24. * 60**2) +
mjd.microseconds / (24. * 60**2 * 1000**2))
if modified:
return mjd_daynumber
return 2400000.5 + mjd_daynumber
[docs]def mjd2datetime(mjd):
"""
Convert a Modified Julian Date to datetime via 'unix time' representation.
NB 'unix time' is defined by the casacore/casacore package.
"""
q = quantity("%sd" % mjd)
return datetime.datetime.fromtimestamp(q.to_unix_time())
[docs]def mjd2lst(mjd, position=None):
"""
Converts a Modified Julian Date into Local Apparent Sidereal Time in
seconds at a given position. If position is None, we default to the
reference position of CS002.
mjd -- Modified Julian Date (float, in days)
position -- Position (casacore measure)
"""
dm = measures()
position = position or dm.position(
"ITRF", "%fm" % ITRF_X, "%fm" % ITRF_Y, "%fm" % ITRF_Z
)
dm.do_frame(position)
last = dm.measure(dm.epoch("UTC", "%fd" % mjd), "LAST")
fractional_day = last['m0']['value'] % 1
return fractional_day * 24 * SECONDS_IN_HOUR
[docs]def mjds2lst(mjds, position=None):
"""
As mjd2lst(), but takes an argument in seconds rather than days.
Args:
mjds (float):Modified Julian Date (in seconds)
position (casacore measure): Position for LST calcs
"""
return mjd2lst(mjds/SECONDS_IN_DAY, position)
[docs]def jd2lst(jd, position=None):
"""
Converts a Julian Date into Local Apparent Sidereal Time in seconds at a
given position. If position is None, we default to the reference position
of CS002.
Args:
jd (float): Julian Date
position (casacore measure): Position for LST calcs.
"""
return mjd2lst(jd - 2400000.5, position)
# NB: datetime is not sensitive to leap seconds.
# However, leap seconds were first introduced in 1972.
# So there are no leap seconds between the start of the
# Modified Julian epoch and the start of the Unix epoch,
# so this calculation is safe.
#julian_epoch = datetime.datetime(1858, 11, 17)
#unix_epoch = datetime.datetime(1970, 1, 1, 0, 0)
#delta = unix_epoch - julian_epoch
#deltaseconds = delta.total_seconds()
#unix_epoch = 3506716800
# The above is equivalent to this:
unix_epoch = quantity("1970-01-01T00:00:00").get_value('s')
[docs]def julian2unix(timestamp):
"""
Convert a modifed julian timestamp (number of seconds since 17 November
1858) to Unix timestamp (number of seconds since 1 January 1970).
Args:
timestamp (numbers.Number): Number of seconds since the Unix epoch.
Returns:
numbers.Number: Number of seconds since the modified Julian epoch.
"""
return timestamp - unix_epoch
[docs]def unix2julian(timestamp):
"""
Convert a Unix timestamp (number of seconds since 1 January 1970) to a
modified Julian timestamp (number of seconds since 17 November 1858).
Args:
timestamp (numbers.Number): Number of seconds since the modified
Julian epoch.
Returns:
numbers.Number: Number of seconds since the Unix epoch.
"""
return timestamp + unix_epoch
[docs]def sec2deg(seconds):
"""Seconds of time to degrees of arc"""
return 15.0 * seconds / 3600.0
[docs]def sec2days(seconds):
"""Seconds to number of days"""
return seconds / (24.0 * 3600)
[docs]def sec2hms(seconds):
"""Seconds to hours, minutes, seconds"""
hours, seconds = divmod(seconds, 60**2)
minutes, seconds = divmod(seconds, 60)
return (int(hours), int(minutes), seconds)
[docs]def altaz(mjds, ra, dec, lat=CORE_LAT):
"""Calculates the azimuth and elevation of source from time and position
on sky. Takes MJD in seconds and ra, dec in degrees. Returns (alt, az) in
degrees."""
# compute hour angle in degrees
ha = mjds2lst(mjds) - ra
if (ha < 0):
ha = ha + 360
# convert degrees to radians
ha, dec, lat = [math.radians(value) for value in (ha, dec, lat)]
# compute altitude in radians
sin_alt = (math.sin(dec) * math.sin(lat) +
math.cos(dec) * math.cos(lat) * math.cos(ha))
alt = math.asin(sin_alt)
# compute azimuth in radians
# divide by zero error at poles or if alt = 90 deg
cos_az = ((math.sin(dec) - math.sin(alt) * math.sin(lat)) /
(math.cos(alt) * math.cos(lat)))
az = math.acos(cos_az)
# convert radians to degrees
hrz_altitude, hrz_azimuth = [math.degrees(value) for value in (alt, az)]
# choose hemisphere
if (math.sin(ha) > 0):
hrz_azimuth = 360 - hrz_azimuth
return hrz_altitude, hrz_azimuth
[docs]def ratohms(radegs):
"""Convert RA in decimal degrees format to hours, minutes,
seconds format.
Keyword arguments:
radegs -- RA in degrees format
Return value:
ra -- tuple of 3 values, [hours,minutes,seconds]
"""
radegs %= 360
raseconds = radegs * 3600 / 15.0
return sec2hms(raseconds)
[docs]def dectodms(decdegs):
"""Convert Declination in decimal degrees format to hours, minutes,
seconds format.
Keyword arguments:
decdegs -- Dec. in degrees format
Return value:
dec -- list of 3 values, [degrees,minutes,seconds]
"""
sign = -1 if decdegs < 0 else 1
decdegs = abs(decdegs)
if decdegs > 90:
raise ValueError("coordinate out of range")
decd = int(decdegs)
decm = int((decdegs - decd) * 60)
decs = (((decdegs - decd) * 60) - decm) * 60
# Necessary because of potential roundoff errors
if decs - 60 > -1e-7:
decm += 1
decs = 0
if decm == 60:
decd += 1
decm = 0
if decd > 90:
raise ValueError("coordinate out of range")
if sign == -1:
if decd == 0:
if decm == 0:
decs = -decs
else:
decm = -decm
else:
decd = -decd
return (decd, decm, decs)
[docs]def propagate_sign(val1, val2, val3):
"""
casacore (reasonably enough) demands that a minus sign (if required)
comes at the start of the quantity. Thus "-0D30M" rather than "0D-30M".
Python regards "-0" as equal to "0"; we need to split off a separate sign
field.
If more than one of our inputs is negative, it's not clear what the user
meant: we raise.
Args:
val1(float): (,val2,val3) input values (hour/min/sec or deg/min/sec)
Returns:
tuple: "+" or "-" string denoting sign,
val1, val2, val3 (numeric) denoting absolute values of inputs.
"""
signs = [x<0 for x in (val1, val2, val3)]
if signs.count(True) == 0:
sign = "+"
elif signs.count(True) == 1:
sign, val1, val2, val3 = "-", abs(val1), abs(val2), abs(val3)
else:
raise ValueError("Too many negative coordinates")
return sign, val1, val2, val3
[docs]def hmstora(rah, ram, ras):
"""Convert RA in hours, minutes, seconds format to decimal
degrees format.
Keyword arguments:
rah,ram,ras -- RA values (h,m,s)
Return value:
radegs -- RA in decimal degrees
"""
sign, rah, ram, ras = propagate_sign(rah, ram, ras)
ra = quantity("%s%dH%dM%f" % (sign, rah, ram, ras)).get_value()
if abs(ra) >= 360:
raise ValueError("coordinates out of range")
return ra
[docs]def dmstodec(decd, decm, decs):
"""Convert Dec in degrees, minutes, seconds format to decimal
degrees format.
Keyword arguments:
decd, decm, decs -- list of Dec values (d,m,s)
Return value:
decdegs -- Dec in decimal degrees
"""
sign, decd, decm, decs = propagate_sign(decd, decm, decs)
dec = quantity("%s%dD%dM%f" % (sign, decd, decm, decs)).get_value()
if abs(dec) > 90:
raise ValueError("coordinates out of range")
return dec
[docs]def angsep(ra1, dec1, ra2, dec2):
"""Find the angular separation of two sources, in arcseconds,
using the proper spherical trig formula
Keyword arguments:
ra1,dec1 - RA and Dec of the first source, in decimal degrees
ra2,dec2 - RA and Dec of the second source, in decimal degrees
Return value:
angsep - Angular separation, in arcseconds
"""
b = (math.pi / 2) - math.radians(dec1)
c = (math.pi / 2) - math.radians(dec2)
temp = (math.cos(b) * math.cos(c)) + (math.sin(b) * math.sin(c) * math.cos(math.radians(ra1 - ra2)))
# Truncate the value of temp at +- 1: it makes no sense to do math.acos()
# of a value outside this range, but occasionally we might get one due to
# rounding errors.
if abs(temp) > 1.0:
temp = 1.0 * cmp(temp, 0)
return 3600 * math.degrees(math.acos(temp))
[docs]def alphasep(ra1, ra2, dec1, dec2):
"""Find the angular separation of two sources in RA, in arcseconds
Keyword arguments:
ra1,dec1 - RA and Dec of the first source, in decimal degrees
ra2,dec2 - RA and Dec of the second source, in decimal degrees
Return value:
angsep - Angular separation, in arcseconds
"""
return 3600 * (ra1 - ra2) * math.cos(math.radians((dec1 + dec2) / 2.0))
[docs]def deltasep(dec1, dec2):
"""Find the angular separation of two sources in Dec, in arcseconds
Keyword arguments:
dec1 - Dec of the first source, in decimal degrees
dec2 - Dec of the second source, in decimal degrees
Return value:
angsep - Angular separation, in arcseconds
"""
return 3600 * (dec1 - dec2)
# Find angular separation in Dec of 2 positions, in arcseconds
[docs]def alpha(l, m, alpha0, delta0):
"""Convert a coordinate in l,m into an coordinate in RA
Keyword arguments:
l,m -- direction cosines, given by (offset in cells) x cellsi (radians)
alpha_0, delta_0 -- centre of the field
Return value:
alpha -- RA in decimal degrees
"""
return (alpha0 + (math.degrees(math.atan2(l, (
(math.sqrt(1 - (l*l) - (m*m)) * math.cos(math.radians(delta0))) -
(m * math.sin(math.radians(delta0))))))))
[docs]def alpha_inflate(theta, decl):
"""Compute the ra expansion for a given theta at a given declination
Keyword arguments:
theta, decl are both in decimal degrees.
Return value:
alpha -- RA inflation in decimal degrees
For a derivation, see MSR TR 2006 52, Section 2.1
http://research.microsoft.com/apps/pubs/default.aspx?id=64524
"""
if abs(decl) + theta > 89.9:
return 180.0
else:
return math.degrees(abs(math.atan(math.sin(math.radians(theta)) / math.sqrt(abs(math.cos(math.radians(decl - theta)) * math.cos(math.radians(decl + theta)))))))
# Find the RA of a point in a radio image, given l,m and field centre
[docs]def delta(l, m, delta0):
"""Convert a coordinate in l, m into an coordinate in Dec
Keyword arguments:
l, m -- direction cosines, given by (offset in cells) x cellsi (radians)
alpha_0, delta_0 -- centre of the field
Return value:
delta -- Dec in decimal degrees
"""
return math.degrees(math.asin(m * math.cos(math.radians(delta0)) +
(math.sqrt(1 - (l*l) - (m*m)) *
math.sin(math.radians(delta0)))))
[docs]def l(ra, dec, cra, incr):
"""Convert a coordinate in RA,Dec into a direction cosine l
Keyword arguments:
ra,dec -- Source location
cra -- RA centre of the field
incr -- number of degrees per pixel (negative in the case of RA)
Return value:
l -- Direction cosine
"""
return ((math.cos(math.radians(dec)) * math.sin(math.radians(ra - cra))) /
(math.radians(incr)))
[docs]def m(ra, dec, cra, cdec, incr):
"""Convert a coordinate in RA,Dec into a direction cosine m
Keyword arguments:
ra,dec -- Source location
cra,cdec -- centre of the field
incr -- number of degrees per pixel
Return value:
m -- direction cosine
"""
return ((math.sin(math.radians(dec)) * math.cos(math.radians(cdec))) -
(math.cos(math.radians(dec)) * math.sin(math.radians(cdec)) *
math.cos(math.radians(ra-cra)))) / math.radians(incr)
[docs]def lm_to_radec(ra0, dec0, l, m):
"""
Find the l direction cosine in a radio image, given an RA and Dec and the
field centre
"""
# This function should be the inverse of radec_to_lmn, but it is
# not. There is likely an error here.
sind0 = math.sin(dec0)
cosd0 = math.cos(dec0)
dl = l
dm = m
d0 = dm * dm * sind0 * sind0 + dl * dl - 2 * dm * cosd0 * sind0
sind = math.sqrt(abs(sind0 * sind0 - d0))
cosd = math.sqrt(abs(cosd0 * cosd0 + d0))
if (sind0 > 0):
sind = abs(sind)
else:
sind = -abs(sind)
dec = math.atan2(sind, cosd)
if l != 0:
ra = math.atan2(-dl, (cosd0 - dm * sind0)) + ra0
else:
ra = math.atan2((1e-10), (cosd0 - dm * sind0)) + ra0
# Calculate RA,Dec from l,m and phase center. Note: As done in
# Meqtrees, which seems to differ from l, m functions above. Meqtrees
# equation may have problems, judging from my difficulty fitting a
# fringe to L4086 data. Pandey's equation is now used in radec_to_lmn
return (ra, dec)
[docs]def radec_to_lmn(ra0, dec0, ra, dec):
l = math.cos(dec) * math.sin(ra - ra0)
sind0 = math.sin(dec0)
if sind0 != 0:
# from pandey; gives same results for casa and cyga
m = (math.sin(dec) * math.cos(dec0) -
math.cos(dec) * math.sin(dec0) * math.cos(ra - ra0))
else:
m = 0
n = math.sqrt(1 - l**2 - m**2)
return (l, m, n)
[docs]def eq_to_gal(ra, dec):
"""Find the Galactic co-ordinates of a source given the equatorial
co-ordinates
Keyword arguments:
(alpha,delta) -- RA, Dec in decimal degrees
Return value:
(l,b) -- Galactic longitude and latitude, in decimal degrees
"""
dm = measures()
result = dm.measure(
dm.direction("J200", "%fdeg" % ra, "%fdeg" % dec),
"GALACTIC"
)
lon_l = math.degrees(result['m0']['value']) % 360 # 0 < ra < 360
lat_b = math.degrees(result['m1']['value'])
return lon_l, lat_b
[docs]def gal_to_eq(lon_l, lat_b):
"""Find the Galactic co-ordinates of a source given the equatorial
co-ordinates
Keyword arguments:
(l, b) -- Galactic longitude and latitude, in decimal degrees
Return value:
(alpha, delta) -- RA, Dec in decimal degrees
"""
dm = measures()
result = dm.measure(
dm.direction("GALACTIC", "%fdeg" % lon_l, "%fdeg" % lat_b),
"J2000"
)
ra = math.degrees(result['m0']['value']) % 360 # 0 < ra < 360
dec = math.degrees(result['m1']['value'])
return ra, dec
[docs]def eq_to_cart(ra, dec):
"""Find the cartesian co-ordinates on the unit sphere given the eq. co-ords.
ra, dec should be in degrees.
"""
return (math.cos(math.radians(dec)) * math.cos(math.radians(ra)), # Cartesian x
math.cos(math.radians(dec)) * math.sin(math.radians(ra)), # Cartesian y
math.sin(math.radians(dec))) # Cartesian z
[docs]class CoordSystem(object):
"""A container for constant strings representing different coordinate
systems."""
FK4 = "B1950 (FK4)"
FK5 = "J2000 (FK5)"
[docs]def coordsystem(name):
"""Given a string, return a constant from class CoordSystem."""
mappings = {
'j2000': CoordSystem.FK5,
'fk5': CoordSystem.FK5,
CoordSystem.FK5.lower(): CoordSystem.FK5,
'b1950': CoordSystem.FK4,
'fk4': CoordSystem.FK4,
CoordSystem.FK4.lower(): CoordSystem.FK4
}
return mappings[name.lower()]
[docs]def convert_coordsystem(ra, dec, insys, outsys):
"""
Convert RA & dec (given in decimal degrees) between equinoxes.
"""
dm = measures()
if insys == CoordSystem.FK4:
insys = "B1950"
elif insys == CoordSystem.FK5:
insys = "J2000"
else:
raise Exception("Unknown Coordinate System")
if outsys == CoordSystem.FK4:
outsys = "B1950"
elif outsys == CoordSystem.FK5:
outsys = "J2000"
else:
raise Exception("Unknown Coordinate System")
result = dm.measure(
dm.direction(insys, "%fdeg" % ra, "%fdeg" % dec),
outsys
)
ra = math.degrees(result['m0']['value']) % 360 # 0 < ra < 360
dec = math.degrees(result['m1']['value'])
return ra, dec
[docs]class WCS(object):
"""
Wrapper around pywcs.WCS.
This is primarily to preserve API compatibility with the earlier,
home-brewed python-wcslib wrapper. It includes:
* A fix for the reference pixel lying at the zenith;
* Raises ValueError if coordinates are invalid.
"""
# ORIGIN is the upper-left corner of the image. pywcs supports both 0
# (NumPy, C-style) or 1 (FITS, Fortran-style). The TraP uses 1.
ORIGIN = 1
# We can set these attributes on the pywcs.WCS().wcs object to configure
# the coordinate system.
WCS_ATTRS = ("crpix", "cdelt", "crval", "ctype", "cunit", "crota")
def __init__(self):
# Currently, we only support two dimensional images.
self.wcs = pywcs.WCS(naxis=2)
def __setattr__(self, attrname, value):
if attrname in self.WCS_ATTRS:
# Account for arbitrary coordinate rotations in images pointing at
# the North Celestial Pole. We set the reference direction to
# infintesimally less than 90 degrees to avoid any ambiguity. See
# discussion at #4599.
if attrname == "crval" and (value[1] == 90 or value[1] == math.pi/2):
value = (value[0], value[1] * (1 - sys.float_info.epsilon))
self.wcs.wcs.__setattr__(attrname, value)
else:
super(WCS, self).__setattr__(attrname, value)
def __getattr__(self, attrname):
if attrname in self.WCS_ATTRS:
return getattr(self.wcs.wcs, attrname)
else:
super(WCS, self).__getattr__(attrname)
[docs] def p2s(self, pixpos):
"""
Pixel to Spatial coordinate conversion.
Args:
pixpos (list): [x, y] pixel position
Returns:
ra (float): Right ascension corresponding to position [x, y]
dec (float): Declination corresponding to position [x, y]
"""
[ra], [dec] = self.wcs.wcs_pix2sky(pixpos[0], pixpos[1], self.ORIGIN)
if math.isnan(ra) or math.isnan(dec):
raise RuntimeError("Spatial position is not a number")
return ra, dec
[docs] def s2p(self, spatialpos):
"""
Spatial to Pixel coordinate conversion.
Args:
pixpos (list): [ra, dec] spatial position
Returns:
x (float): X pixel value corresponding to position [ra, dec]
y (float): Y pixel value corresponding to position [ra, dec]
"""
[x], [y] = self.wcs.wcs_sky2pix(spatialpos[0], spatialpos[1], self.ORIGIN)
if math.isnan(x) or math.isnan(y):
raise RuntimeError("Pixel position is not a number")
return x, y