diff --git a/tests/test_database/test_subroutines.py b/tests/test_database/test_subroutines.py index cd97183d..585a9a47 100644 --- a/tests/test_database/test_subroutines.py +++ b/tests/test_database/test_subroutines.py @@ -2,10 +2,10 @@ import unittest +from sourcefinder.utility.coordinates import eq_to_cart as py_cartesian import tkp.db from tkp.testutil.decorators import requires_database from tkp.testutil.db_queries import convert_to_cartesian as db_cartesian -from tkp.utility.coordinates import eq_to_cart as py_cartesian """Test miscellaneous minor database functions""" diff --git a/tests/test_utility/test_coordinates.py b/tests/test_utility/test_coordinates.py deleted file mode 100644 index 1096911a..00000000 --- a/tests/test_utility/test_coordinates.py +++ /dev/null @@ -1,313 +0,0 @@ -""" -Test functions used for manipulating coordinates in the TKP pipeline. -""" - -import pytz - -import unittest -from casacore.measures import measures - -import datetime -from tkp.utility import coordinates - - -class mjd2lstTest(unittest.TestCase): - # Note that the precise conversion can vary slightly depending on the - # casacore measures database, so we allow for a 1s uncertainty in the - # results - MJD = 56272.503491875 - MJDs = coordinates.SECONDS_IN_DAY * MJD - - def testDefaultLocation(self): - last = 64496.73666805029 - self.assertTrue(abs(coordinates.mjd2lst(self.MJD) - last) < 1.0) - self.assertTrue(abs(coordinates.mjds2lst(self.MJDs) - last) < 1.0) - - def testSpecificLocation(self): - last = 73647.97547170892 - dm = measures() - position = dm.position("itrf", "1m", "1m", "1m") - self.assertTrue(abs(coordinates.mjd2lst(self.MJD, position) - last) < 1.0) - self.assertTrue(abs(coordinates.mjds2lst(self.MJDs, position) - last) < 1.0) - - -class jd2lsttest(unittest.TestCase): - # Note that the precise conversion can vary slightly depending on the - # casacore measures database, so we allow for a 1s uncertainty in the - # results - JD = 56272.503491875 + 2400000.5 - - def testdefaultlocation(self): - last = 64496.73666805029 - self.assertTrue(abs(coordinates.jd2lst(self.JD) - last) < 1.0) - - def testspecificlocation(self): - last = 73647.97547170892 - dm = measures() - position = dm.position("itrf", "1m", "1m", "1m") - self.assertTrue(abs(coordinates.jd2lst(self.JD, position) - last) < 1.0) - - -class gal_to_eq_and_eq_to_gal_Test(unittest.TestCase): - knownValues = ( - [(10, 10), (262.89288972929216, -15.207965232900214)], - [(350, -10), (271.0071629114459, -42.544903877104645)], - [(1, 1), (266.0298496879791, -27.561144848129747)], - [(118.27441982715264, -52.7682822322771), (10, 10)], - [(67.05481764677926, -62.47515197465108), (350, -10)], - [(98.9410284460564, -59.6437963999095), (1, 1)] - ) - - def testKnownValues(self): - for coord_pair in self.knownValues: - gal_l, gal_b = coord_pair[0] - ra, dec = coord_pair[1] - - self.assertAlmostEqual(coordinates.gal_to_eq(gal_l, gal_b)[0], ra, 3) - self.assertAlmostEqual(coordinates.gal_to_eq(gal_l, gal_b)[1], dec, 3) - - self.assertAlmostEqual(coordinates.eq_to_gal(ra, dec)[0], gal_l, 3) - self.assertAlmostEqual(coordinates.eq_to_gal(ra, dec)[1], gal_b, 3) - - -class ratohmsTest(unittest.TestCase): - knownValues = ((0, (0, 0, 0)), - (90, (6, 0, 0)), - (180, (12, 0, 0)), - (-180, (12, 0, 0)), - (-540, (12, 0, 0)), - (270, (18, 0, 0)), - (360, (0, 0, 0)), - (450, (6, 0, 0)), - (0.1, (0, 0, 24)), - (1.0, (0, 4, 0)), - ) - - def testknownValues(self): - for input, output in self.knownValues: - result = coordinates.ratohms(input) - self.assertEqual(result, output) - - def testTypes(self): - self.assertRaises(TypeError, coordinates.ratohms, 'a') - - -class dectodmsTest(unittest.TestCase): - knownValues = ((0, (0, 0, 0)), - (45, (45, 0, 0)), - (-45, (-45, 0, 0)), - (90, (90, 0, 0)), - (-90, (-90, 0, 0)), - (0.5, (0, 30, 0)), - (1.0/60**2, (0, 0, 1)), - ) - - def testknownValues(self): - for input, output in self.knownValues: - result = coordinates.dectodms(input) - self.assertEqual(result, output) - - def testTypes(self): - self.assertRaises(TypeError, coordinates.dectodms, 'a') - - def testRange(self): - self.assertRaises(ValueError, coordinates.dectodms, 91) - -class propagate_signTest(unittest.TestCase): - knownValues = ( - ((0, 0, 0), ("+", 0, 0, 0)), - ((1, 0, 0), ("+", 1, 0, 0)), - ((0, 1, 0), ("+", 0, 1, 0)), - ((0, 0, 1), ("+", 0, 0, 1)), - ((0, 0, -1), ("-", 0, 0, 1)), - ((0, -1, 0), ("-", 0, 1, 0)), - ((-1, 0, 0), ("-", 1, 0, 0)) - ) - - def testknownValues(self): - for input, output in self.knownValues: - result = coordinates.propagate_sign(*input) - self.assertEqual(result, output) - - def testCrazyInput(self): - self.assertRaises(ValueError, coordinates.propagate_sign, -1, -1, -1) - -class hmstoraTest(unittest.TestCase): - knownValues = (((0, 0, 0), 0), - ((6, 0, 0), 90), - ((12, 0, 0), 180), - ((18, 0, 0), 270), - ((0, 0, 24), 0.1), - ((0, 4, 0), 1.0), - ((0, 0, 59.1), 0.24625), - ) - - def testknownValues(self): - for input, output in self.knownValues: - result = coordinates.hmstora(*input) - self.assertAlmostEqual(result, output) - - def testTypes(self): - self.assertRaises(TypeError, coordinates.hmstora, 'a') - - def testSanity(self): - for h in range(0, 24): - for m in range(0, 60): - for s in range(0, 60): - ra = coordinates.dmstodec(h, m, s) - ch, cm, cs = coordinates.dectodms(ra) - self.assertEqual(h, ch) - self.assertEqual(m, cm) - self.assertAlmostEqual(s, cs) - - def testRange(self): - self.assertRaises(ValueError, coordinates.hmstora, 24, 0, 1) - self.assertRaises(ValueError, coordinates.hmstora, 25, 0, 0) - -class dmstodecTest(unittest.TestCase): - knownValues = (((0, 0, 0), 0), - ((45, 0, 0), 45), - ((-45, 0, 0), -45), - ((90, 0, 0), 90), - ((-90, 0, 0), -90), - ((0, 30, 0), 0.5), - ((0, 0, 1), 1.0/60**2), - ((0, 0, 59.004), 59.004/60**2) - ) - - def testknownValues(self): - for input, output in self.knownValues: - result = coordinates.dmstodec(*input) - self.assertAlmostEqual(result, output) - - def testTypes(self): - self.assertRaises(TypeError, coordinates.dmstodec, 'a') - - def testSanity(self): - for d in range(-89, 90): - for m in range(0, 60): - for s in range(0, 60): - dec = coordinates.dmstodec(d, m, s) - cd, cm, cs = coordinates.dectodms(dec) - self.assertEqual(d, cd) - self.assertEqual(m, cm) - self.assertAlmostEqual(s, cs) - dec = coordinates.dmstodec(90, 0, 0) - cd, cm, cs = coordinates.dectodms(dec) - self.assertEqual(90, cd) - self.assertEqual(0, cm) - self.assertAlmostEqual(0, cs) - dec = coordinates.dmstodec(0, -30, 0) - cd, cm, cs = coordinates.dectodms(dec) - self.assertEqual(0, cd) - self.assertEqual(-30, cm) - self.assertAlmostEqual(0, cs) - - def testRange(self): - self.assertRaises(ValueError, coordinates.dmstodec, 91, 0, 0) - self.assertRaises(ValueError, coordinates.dmstodec, 90, 0, 1) - self.assertRaises(ValueError, coordinates.dmstodec, -90, 0, 1) - - -class juliandate(unittest.TestCase): - knownValues = ( - (datetime.datetime(2008, 4, 17, 12, 30, tzinfo=pytz.utc), 2454574.0208333), - (datetime.datetime(2007, 1, 14, 13, 18, 59, 900000, tzinfo=pytz.utc), 2454115.05486), - (datetime.datetime(1858, 11, 17, tzinfo=pytz.utc), 2400000.5), - ) - - def testknownValues(self): - for date, jd in self.knownValues: - calc_jd = coordinates.julian_date(date) - self.assertAlmostEqual(calc_jd, jd) - calc_mjd = coordinates.julian_date(date, modified=True) - self.assertAlmostEqual(calc_mjd, jd - 2400000.5) - - def testNow(self): - now = coordinates.julian_date() - self.assertTrue(now > 2454574) - -class coordsystemTest(unittest.TestCase): - knownValues = ( - {"fk4": (10.0, 10.0), "fk5": (10.64962347, 10.273829)}, - {"fk4": (9.351192, 9.725562), "fk5": (10.0, 10.0)}, - {"fk4": (179.357791, 45.278329), "fk5": (180, 45)}, - {"fk4": (180, 45), "fk5": (180.63913, 44.721730)}, - {"fk4": (349.35037, -10.27392), "fk5": (-10, -10)}, - {"fk4": (-10, -10), "fk5": (350.648870, -9.725590)} - ) - - def testKnownValues(self): - # Note that RA is always positive in the range 0 < RA < 360 - for coord_pair in self.knownValues: - ra, dec = coordinates.convert_coordsystem( - coord_pair["fk4"][0], coord_pair["fk4"][1], - coordinates.CoordSystem.FK4, - coordinates.CoordSystem.FK4 - ) - self.assertAlmostEqual(ra, coord_pair["fk4"][0] % 360, 3) - self.assertAlmostEqual(dec, coord_pair["fk4"][1], 3) - - ra, dec = coordinates.convert_coordsystem( - coord_pair["fk5"][0], coord_pair["fk5"][1], - coordinates.CoordSystem.FK4, - coordinates.CoordSystem.FK4 - ) - self.assertAlmostEqual(ra, coord_pair["fk5"][0] % 360, 3) - self.assertAlmostEqual(dec, coord_pair["fk5"][1], 3) - - ra, dec = coordinates.convert_coordsystem( - coord_pair["fk4"][0], coord_pair["fk4"][1], - coordinates.CoordSystem.FK4, - coordinates.CoordSystem.FK5 - ) - self.assertAlmostEqual(ra, coord_pair["fk5"][0] % 360, 3) - self.assertAlmostEqual(dec, coord_pair["fk5"][1], 3) - - ra, dec = coordinates.convert_coordsystem( - coord_pair["fk5"][0], coord_pair["fk5"][1], - coordinates.CoordSystem.FK5, - coordinates.CoordSystem.FK4 - ) - self.assertAlmostEqual(ra, coord_pair["fk4"][0] % 360, 3) - self.assertAlmostEqual(dec, coord_pair["fk4"][1], 3) - - -class angsepTest(unittest.TestCase): - def testZeroSeparation(self): - self.assertEqual(coordinates.angsep(0, 0, 0, 0), 0) - - def testSmallSeparation(self): - # Due to float rounding issues, feeding these values to - # coordinates.angsep() can cause us to attempt math.acos(x) for x - # greater than 1, which throws a ValueError. - ra1 = 33.655858149999993145229382207617 - dec1 = 87.061899808796610500394308473915 - ra2 = 33.655860050872931310550484340638 - dec2 = 87.061899872535235545001341961324 - coordinates.angsep(ra1, dec1, ra2, dec2) - - -class altazTest(unittest.TestCase): - def testExecute(self): - # Simply testing the function is usable, because earlier changes had - # broken it completely. This does not check it gives the right - # results! - coordinates.altaz(55000 * coordinates.SECONDS_IN_DAY, 1, 1) - - -class UnixToMJDTestCase(unittest.TestCase): - def test_unix_epoch(self): - # By definition; this is the number of seconds between - # 1858-11-17 00:00 and 1970-01-01 00:00:00. - self.assertEqual(coordinates.unix_epoch, 3506716800) - - def test_julian2unix(self): - self.assertEqual(coordinates.julian2unix(coordinates.unix_epoch), 0) - - def test_unix2julian(self): - self.assertEqual(coordinates.unix2julian(0), coordinates.unix_epoch) - - -if __name__ == '__main__': - unittest.main() diff --git a/tests/test_utility/test_uncertain.py b/tests/test_utility/test_uncertain.py deleted file mode 100644 index 8e9f49aa..00000000 --- a/tests/test_utility/test_uncertain.py +++ /dev/null @@ -1,12 +0,0 @@ -import unittest -from tkp.utility.uncertain import Uncertain - -class UncertainTestCase(unittest.TestCase): - def test_castable_to_floating_point(self): - """ - float(Uncertain(int(x))) should not raise. - """ - x = Uncertain(int(0), int(0)) - self.assertFalse(isinstance(x.value, float)) - self.assertFalse(isinstance(x.error, float)) - float(x) diff --git a/tests/test_utility/test_wcs.py b/tests/test_utility/test_wcs.py deleted file mode 100644 index be151bdf..00000000 --- a/tests/test_utility/test_wcs.py +++ /dev/null @@ -1,120 +0,0 @@ -import unittest - -from tkp.utility import coordinates -from sourcefinder import extract -from tkp.utility.uncertain import Uncertain - - -# Specify the number of digits you want to include when checking if positions are equal. -nod=12 - -class TestNCP(unittest.TestCase): - """ - Check that we retrieve the correct position for objects in images centred - on the North Celestial Pole. - - At the NCP (dec=90), our coordinate system becomes ambiguous. We avoid the - problem by subtracting an infintesimal quantity from the reference dec, - such that it is never quite 90 degrees. - - See discussion at issue #4599. - """ - def test_3c61(self): - # Coordinate system and position of 3C 61.1 based on an image of the - # NCP provided by Adam Stewart. - wcs = coordinates.WCS() - wcs.ctype = ('RA---SIN', 'DEC--SIN') - wcs.crval = (15.0, 90.0) - wcs.cdelt = (-0.01111111111111, 0.01111111111111) - wcs.crpix = (1025.0, 1025.0) - wcs.unit = ("deg", "deg") - calculated_position = wcs.p2s([908, 715]) - self.assertAlmostEqual(calculated_position[0], 35.7, 1) - self.assertAlmostEqual(calculated_position[1], 86.3, 1) - - -class Sanity(unittest.TestCase): - """Some sanity checks because of issue #2787. - - Previously, 1 was added or subtracted in the - pixel <-> sky conversion, resulting in positions - being one pixel off in both x and y - """ - - def setUp(self): - self.wcs = coordinates.WCS() - self.wcs.crota = [0, 0] - self.wcs.cdelt = [1, 1] - self.wcs.crpix = [10, 10] - self.wcs.crval = [10, 10] - - def tests2p(self): - x, y = self.wcs.s2p([10, 10]) - self.assertAlmostEqual(x, 10) - self.assertAlmostEqual(y, 10) - x, y = self.wcs.s2p([0, 0]) - self.assertAlmostEqual(x, 0) - self.assertAlmostEqual(y, 0) - x, y = self.wcs.s2p([1, 1]) - self.assertAlmostEqual(x, 1) - self.assertAlmostEqual(y, 1) - - def testp2s(self): - r, d = self.wcs.p2s([10, 10]) - self.assertAlmostEqual(r, 10) - self.assertAlmostEqual(d, 10) - r, d = self.wcs.p2s([0, 0]) - self.assertAlmostEqual(r, 0) - self.assertAlmostEqual(d, 0) - r, d = self.wcs.p2s([1, 1]) - self.assertAlmostEqual(r, 1) - self.assertAlmostEqual(d, 1) - - -class wcsSin(unittest.TestCase): - known_values = ( - ([1442.0, 1442.0], [350.785563529949, 58.848317299883],), - ([1441.0, 1501.0], [350.85000000000002, 60.815406379421418],), - ([1441.0, 1381.0], [350.85000000000002, 56.814593620578577],), - ([1381.0, 1441.0], [354.70898191482644, 58.757358624100824],), - ([1501.0, 1441.0], [346.99101808517361, 58.757358624100824],), - ) - - # "header" parameters for this test - # (Taken from arbitrary FITS file) - crval = (3.508500000000E+02, 5.881500000000E+01) - crpix = (1.441000000000E+03, 1.441000000000E+03) - cdelt = (-3.333333333333E-02, 3.333333333333E-02) - ctype = ('RA---SIN', 'DEC--SIN') - crota = (0, 0) - cunit = ('deg', 'deg') - - def setUp(self): - # Initialise a wcs object with the above parameters - self.wcs = coordinates.WCS() - self.wcs.crval = self.crval - self.wcs.crpix = self.crpix - self.wcs.cdelt = self.cdelt - self.wcs.ctype = self.ctype - self.wcs.crota = self.crota - self.wcs.cunit = self.cunit - - def testPixelToSpatial(self): - for pixel, spatial in self.known_values: - result = self.wcs.p2s(pixel) - self.assertEqual([round(result[0],nod),round(result[1],nod)], [round(spatial[0],nod),round(spatial[1],nod)]) - - def testSpatialToPixel(self): - for pixel, spatial in self.known_values: - result = list(map(round, self.wcs.s2p(spatial))) - self.assertEqual(result, pixel) - - def testSanity(self): - import random - pixel = [random.randrange(500, 1500), random.randrange(500, 1500)] - result = list(map(round, self.wcs.s2p(self.wcs.p2s(pixel)))) - self.assertEqual(result, pixel) - - -if __name__ == '__main__': - unittest.main() diff --git a/tkp/db/alchemy/image.py b/tkp/db/alchemy/image.py index 1b0f09c6..2a1aa223 100644 --- a/tkp/db/alchemy/image.py +++ b/tkp/db/alchemy/image.py @@ -1,7 +1,7 @@ import math from datetime import datetime +from sourcefinder.utility.coordinates import eq_to_cart from tkp.db.model import Frequencyband, Skyregion, Image, Dataset -from tkp.utility.coordinates import eq_to_cart from sqlalchemy import func, cast from sqlalchemy.dialects.postgresql import DOUBLE_PRECISION as Double diff --git a/tkp/db/general.py b/tkp/db/general.py index 049c0751..129fe479 100644 --- a/tkp/db/general.py +++ b/tkp/db/general.py @@ -11,11 +11,11 @@ import tkp.db from datetime import datetime +from sourcefinder.utility.coordinates import alpha_inflate +from sourcefinder.utility.coordinates import eq_to_cart from tkp.db.alchemy.image import insert_dataset as alchemy_insert_dataset from tkp.db.generic import columns_from_table from tkp.utility import substitute_inf -from tkp.utility.coordinates import alpha_inflate -from tkp.utility.coordinates import eq_to_cart logger = logging.getLogger(__name__) diff --git a/tkp/quality/brightsource.py b/tkp/quality/brightsource.py index 50e54bde..f9d80f3a 100644 --- a/tkp/quality/brightsource.py +++ b/tkp/quality/brightsource.py @@ -2,7 +2,7 @@ import logging from io import BytesIO from casacore.measures import measures -from tkp.utility.coordinates import unix2julian +from sourcefinder.utility.coordinates import unix2julian from tkp.utility.redirect_stream import redirect_stream logger = logging.getLogger(__name__) diff --git a/tkp/testutil/db_subs.py b/tkp/testutil/db_subs.py index 4c173549..3fd187af 100644 --- a/tkp/testutil/db_subs.py +++ b/tkp/testutil/db_subs.py @@ -2,6 +2,8 @@ from collections import namedtuple import datetime, math + +import sourcefinder.utility.coordinates as coords import tkp.db from tkp.db.generic import get_db_rows_as_dicts from tkp.db.database import Database @@ -11,8 +13,6 @@ from tkp.db import nulldetections import tkp.testutil.data as testdata -import tkp.utility.coordinates as coords - ExtractedSourceTuple = namedtuple("ExtractedSourceTuple", ['ra', 'dec' , 'ra_fit_err' , 'dec_fit_err' , diff --git a/tkp/testutil/mock.py b/tkp/testutil/mock.py index 8277e64d..a8fc9474 100644 --- a/tkp/testutil/mock.py +++ b/tkp/testutil/mock.py @@ -3,7 +3,7 @@ """ import numpy as np from sourcefinder.accessors.dataaccessor import DataAccessor -from tkp.utility.coordinates import WCS +from sourcefinder.utility.coordinates import WCS import datetime class Mock(object): @@ -102,4 +102,4 @@ def __init__(self, def calculate_phase_centre(self): x, y = self.data.shape centre_ra, centre_decl = self.wcs.p2s((x / 2, y / 2)) - return centre_ra, centre_decl \ No newline at end of file + return centre_ra, centre_decl diff --git a/tkp/utility/containers.py b/tkp/utility/containers.py deleted file mode 100644 index c1d43e23..00000000 --- a/tkp/utility/containers.py +++ /dev/null @@ -1,88 +0,0 @@ -""" -Container classes for the TKP pipeline. - -These provide convenient means of marshalling the various types of data -- -lightcurves, detections, sources, etc -- that the pipeline must handle. -""" - -import logging - -logger = logging.getLogger(__name__) - - -class ObjectContainer(list): - """A container class for objects. - - What sort of objects? Well, anything that has a position and we - want to keep lists of, really. So detections (ie, an individual - source measurement on an image), sources (ie all the detections of - a given object in a given image stack) and lightcurves (ie, all - the sources associated with a given object through time). - - You probably don't want to use this on it's own: see ExtractionResults, - TargetList or source for more useful derived classes. - """ - - def closest_to(self, pix_x, pix_y): - distance, target = False, False - logger.debug("Beginning a search for objects near %.1f, %.1f: ", - pix_x, pix_y) - logger.debug("%s contains %d objects", str(self), len(self)) - for obj in self: - tmpdist = (pix_x - obj.x) ** 2 + (pix_y - obj.y) ** 2 - logger.debug("Object at %f, %f", obj.x, obj.y) - logger.debug("Has distance %f", tmpdist) - if not distance: - distance = tmpdist - target = obj - else: - if tmpdist < distance: - target = obj - distance = tmpdist - logger.debug("Best distance is now %f", distance) - logger.debug("From object %s", str(target)) - if not distance: - return (target, distance) - else: - return (target, distance ** 0.5) - - def __setslice__(self, section, items): - """ - Not implemented. - """ - raise NotImplementedError - - def __iadd__(self, y): - """ - Not implemented. - """ - raise NotImplementedError - - def __imul__(self, y): - """ - Not implemented. - """ - raise NotImplementedError - - def __mul__(self, y): - """ - Not implemented. - """ - raise NotImplementedError - - def __rmul__(self, y): - """ - Not implemented. - """ - raise NotImplementedError - - def __str__(self): - return 'Container: ' + str(len(self)) + ' object(s).' - - -class ExtractionResults(ObjectContainer): - """Container for the results of running source extraction on an - ImageData object""" - - def __str__(self): - return 'ExtractionResults: ' + str(len(self)) + ' detection(s).' diff --git a/tkp/utility/coordinates.py b/tkp/utility/coordinates.py deleted file mode 100644 index 595b8406..00000000 --- a/tkp/utility/coordinates.py +++ /dev/null @@ -1,698 +0,0 @@ -# -# LOFAR Transients Key Project -""" -General purpose astronomical coordinate handling routines. -""" - -import datetime -import logging -import math -import sys - -import pytz -from astropy import wcs as pywcs -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 - - -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 - - -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()) - - -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 - - -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) - - -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') - - -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 - - -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 - - -def sec2deg(seconds): - """Seconds of time to degrees of arc""" - return 15.0 * seconds / 3600.0 - - -def sec2days(seconds): - """Seconds to number of days""" - return seconds / (24.0 * 3600) - - -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) - - -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 - - -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) - - -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) - - -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 - - -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 - - -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 - -def cmp(a, b): - return (a > b) - (a < b) - -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)) - - -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)) - - -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 -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)))))))) - - -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 -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))))) - - -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))) - - -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) - - -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) - - -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) - - -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 - - -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 - - -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 - - -class CoordSystem(object): - """A container for constant strings representing different coordinate - systems.""" - FK4 = "B1950 (FK4)" - FK5 = "J2000 (FK5)" - - -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()] - - -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 - - -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) - - def p2s(self, pixpos): - """ - Pixel to Spatial coordinate conversion. - - Args: - pixpos (tuple): [x, y] pixel position - - Returns: - tuple: ra (float) Right ascension corresponding to position [x, y] - dec (float) Declination corresponding to position [x, y] - """ - ra, dec = self.wcs.wcs_pix2world(pixpos[0], pixpos[1], self.ORIGIN) - if math.isnan(ra) or math.isnan(dec): - raise RuntimeError("Spatial position is not a number") - return float(ra), float(dec) - - def s2p(self, spatialpos): - """ - Spatial to Pixel coordinate conversion. - - Args: - pixpos (tuple): [ra, dec] spatial position - - Returns: - tuple: X pixel value corresponding to position [ra, dec], - Y pixel value corresponding to position [ra, dec] - """ - x, y = self.wcs.wcs_world2pix(spatialpos[0], spatialpos[1], self.ORIGIN) - if math.isnan(x) or math.isnan(y): - raise RuntimeError("Pixel position is not a number") - return float(x), float(y) diff --git a/tkp/utility/memoize.py b/tkp/utility/memoize.py deleted file mode 100644 index 876636e5..00000000 --- a/tkp/utility/memoize.py +++ /dev/null @@ -1,37 +0,0 @@ -# -# LOFAR Transients Key Project -# -# Memoization. -# -from functools import update_wrapper -from weakref import WeakKeyDictionary - - -class Memoize(object): - """Decorator to cache the results of methods. - - Examples in e.g. image.py:: - - @Memoize - def _grids(self): - return self.__grids() - grids = property(fget=_grids, fdel=_grids.delete) - - """ - - def __init__(self, funct): - self.funct = funct - self.memo = WeakKeyDictionary() - update_wrapper(self, self.funct) - - def __call__(self, instance): - if instance not in self.memo: - self.memo[instance] = self.funct(instance) - return self.memo[instance] - - def delete(self, instance): - """Forget a memoized value""" - try: - del (self.memo[instance]) - except KeyError: - pass diff --git a/tkp/utility/uncertain.py b/tkp/utility/uncertain.py deleted file mode 100644 index 08fbebbc..00000000 --- a/tkp/utility/uncertain.py +++ /dev/null @@ -1,140 +0,0 @@ -# -# LOFAR Transients Key Project -# -# Uncertain quantities. -# Based on: -# (c) Robert Jordens -# Made available freely under the Python license -# http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/535164 -import math - - -class Uncertain(object): - """Represents a numeric value with a known small uncertainty (error, - standard deviation...). - - Numeric operators are overloaded to work with other Uncertain or - numeric objects. The uncertainty (error) must be small. Otherwise - the linearization employed here becomes wrong. - """ - - def __init__(self, value=0., error=0., *a, **t): - self.value = value - self.error = abs(error) - super(Uncertain, self).__init__(*a, **t) - - def __str__(self): - return "%g+-%g" % (self.value, self.error) - - def __repr__(self): - return "Uncertain(%s, %s)" % (self.value, self.error) - - def __float__(self): - return float(self.value) - - def assign(self, other): - if isinstance(other, Uncertain): - self.value = other.value - self.error = other.error - else: - self.value = other - self.error = 0. - - def __abs__(self): - return Uncertain(abs(self.value), self.error) - - def __add__(self, other): - if isinstance(other, Uncertain): - v = self.value + other.value - e = (self.error ** 2 + other.error ** 2) ** .5 - return Uncertain(v, e) - else: - return Uncertain(self.value + other, self.error) - - def __radd__(self, other): - return self + other # __add__ - - def __sub__(self, other): - return self + (-other) # other.__neg__ and __add__ - - def __rsub__(self, other): - return -self + other # __neg__ and __add__ - - def __mul__(self, other): - if isinstance(other, Uncertain): - v = self.value * other.value - e = ((self.error * other.value) ** 2 + - (other.error * self.value) ** 2) ** .5 - return Uncertain(v, e) - else: - return Uncertain(self.value * other, self.error * other) - - def __rmul__(self, other): - return self * other # __mul__ - - def __neg__(self): - return self * -1 # __mul__ - - def __pos__(self): - return self - - def __round__(self, pos): - return round(self.value, pos) - - def __div__(self, other): - return self * (1. / other) # other.__div__ and __mul__ - - def __truediv__(self, other): - return self * (1. / other) - - def __rdiv__(self, other): - return (self / other) ** -1. # __pow__ and __div__ - - def __pow__(self, other): - if isinstance(other, Uncertain): - v = self.value ** other.value - e = ((self.error * other.value * - self.value ** (other.value - 1.0)) ** 2 + - (other.error * math.log(self.value) * - self.value ** other.value) ** 2 - ) ** .5 - return Uncertain(v, e) - else: - return Uncertain(self.value ** other, - self.error * other * self.value ** (other - 1)) - - def __rpow__(self, other): - assert not isinstance(other, Uncertain) - # otherwise other.__pow__ would have been called - return Uncertain(other ** self.value, - self.error * math.log(other) * other ** self.value) - - def __cmp__(self, compare): - try: - return cmp(self.value, compare.value) - except AttributeError: - return cmp(self.value, compare) - - def __lt__(self, compare): - try: - return self.value < compare.value - except AttributeError: - return self.value < compare - - def __gt__(self, compare): - try: - return self.value > compare.value - except AttributeError: - return self.value > compare - - def exp(self): - return math.e ** self - - def log(self): - return Uncertain(math.log(self.value), self.error / self.value) - - def max(self): - return self.value + self.error - - def min(self): - return self.value - self.error