From b4a9166e8b5325d1607c29950cbc1c029914602d Mon Sep 17 00:00:00 2001 From: Marc Pound <22331890+mpound@users.noreply.github.com> Date: Fri, 1 Dec 2023 12:52:45 -0500 Subject: [PATCH 1/4] new core file, veldef function and test --- src/dysh/fits/core.py | 78 ++++++++++++++++++++++++++++++++ src/dysh/fits/tests/test_core.py | 48 ++++++++++++++++++++ 2 files changed, 126 insertions(+) create mode 100644 src/dysh/fits/core.py create mode 100644 src/dysh/fits/tests/test_core.py diff --git a/src/dysh/fits/core.py b/src/dysh/fits/core.py new file mode 100644 index 00000000..9202f917 --- /dev/null +++ b/src/dysh/fits/core.py @@ -0,0 +1,78 @@ +""" +Core functions for FITS/SDFITS +""" + +# Velocity frame conventions, partially stolen from pyspeckit. +# See also Section6.2.5.2 of the GBT observer's guide https://www.gb.nrao.edu/scienceDocs/GBTog.pdf + +frame_dict = { + "VLSR": "LSRK", + "VRAD": "LSRK", + "VELO": "LSRK", + "VOPT": "LSRK", + "LSRD": "LSRD", + "LSRK": "LSRK", + "-LSR": "LSRK", + "-HEL": "heliocentric", + "-BAR": "barycentric", + "BAR": "barycentric", + "BARY": "barycentric", + "-OBS": "obs", + "VHEL": "heliocentric", + "VGEO": "topocentric", + "TOPO": "topocentric", + "VREST": "rest", + "Z": "rest", + "FREQ": "rest", + "WAV": "rest", + "WAVE": "rest", + "CMB": "cmb", + "GALAC": "galactic", + "GALA": "galactic", + "ALAC": "galactic", +} + +# Dictionary to convert from FITS velocity convention to specutils string. +# At GBT, VELO was written by sdfits filler for some unknown amount of +# time instead of RELA, so allow for it here +vconv_dict = { + "OPTI": "doppler_optical", + "RADI": "doppler_radio", + "RELA": "doppler_relativistic", + "VELO": "doppler_relativistic", +} + + +def decode_veldef(veldef): + """ + Parse the SDFITS VELDEF value into its two components, the velocity + definition and velocity reference frame. This value must contain + no more than 8 characters where the first 4 characters describe the velocity + definition and the last 4 characters describe the reference frame. + + Parameters + ---------- + veldef : str + The definition string, consisting of a velocity convention and a velocity frame, e.g., 'OPTI-LSR' + + Returns + ------- + A str tuple of velocity convention and velocity frame type, e.g., ('doppler_radio', 'LSRK') + """ + if len(veldef) > 8: + # in future, possibly relax this requirement + # if string not coming from FITS + raise ValueError(f"VELDEF string {veldef} must be no more than 8 characters.") + vconv = veldef[:4] + try: + velocity_convention = vconv_dict[vconv] + except KeyError: + raise KeyError(f"Velocity convention {vconv} not recognized.") + + frame = veldef[4:] + try: + frame_type = frame_dict[frame] + except KeyError: + raise KeyError(f"Velocity frame {frame} not recognized.") + + return velocity_convention, frame_type diff --git a/src/dysh/fits/tests/test_core.py b/src/dysh/fits/tests/test_core.py new file mode 100644 index 00000000..bfb9060f --- /dev/null +++ b/src/dysh/fits/tests/test_core.py @@ -0,0 +1,48 @@ +from dysh.fits import decode_veldef + + +class TestCore: + """Test dysh.fits core functions""" + + def test_veldef(self): + # first make sure we get correct answers for normal inputs + inputs = [ + "RADILSRK", + "RADI-LSR", + "RADILSRD", + "OPTICMB", + "VELO-BAR", + "OPTIBARY", + "RELATOPO", + "RADIGALA", + "OPTI-HEL", + ] + outputs = [ + ("doppler_radio", "LSRK"), + ("doppler_radio", "LSRK"), + ("doppler_radio", "LSRD"), + ("doppler_optical", "cmb"), + ("doppler_relativistic", "barycentric"), + ("doppler_optical", "barycentric"), + ("doppler_relativistic", "topocentric"), + ("doppler_radio", "galactic"), + ("doppler_optical", "heliocentric"), + ] + for i, j in zip(inputs, outputs): + assert decode_veldef(i) == j + + # Now test that bad input raises an exception + try: + decode_veldef("This is more than 8 chars") + except ValueError: + assert True + try: + # frame fails + decode_veldef("OPTI-LRS") + except KeyError: + assert True + try: + # convention fails + decode_veldef("MAXILSRK") + except KeyError: + assert True From 1f3b836b3cf74b24f74a79d097a42112f553082e Mon Sep 17 00:00:00 2001 From: Marc Pound <22331890+mpound@users.noreply.github.com> Date: Fri, 1 Dec 2023 13:23:29 -0500 Subject: [PATCH 2/4] rename because pytest doesn't like files with the same name --- src/dysh/fits/tests/test_fits_core.py | 48 ++++++++++++++++ src/dysh/spectra/tests/test_spectra_core.py | 62 +++++++++++++++++++++ 2 files changed, 110 insertions(+) create mode 100644 src/dysh/fits/tests/test_fits_core.py create mode 100644 src/dysh/spectra/tests/test_spectra_core.py diff --git a/src/dysh/fits/tests/test_fits_core.py b/src/dysh/fits/tests/test_fits_core.py new file mode 100644 index 00000000..bfb9060f --- /dev/null +++ b/src/dysh/fits/tests/test_fits_core.py @@ -0,0 +1,48 @@ +from dysh.fits import decode_veldef + + +class TestCore: + """Test dysh.fits core functions""" + + def test_veldef(self): + # first make sure we get correct answers for normal inputs + inputs = [ + "RADILSRK", + "RADI-LSR", + "RADILSRD", + "OPTICMB", + "VELO-BAR", + "OPTIBARY", + "RELATOPO", + "RADIGALA", + "OPTI-HEL", + ] + outputs = [ + ("doppler_radio", "LSRK"), + ("doppler_radio", "LSRK"), + ("doppler_radio", "LSRD"), + ("doppler_optical", "cmb"), + ("doppler_relativistic", "barycentric"), + ("doppler_optical", "barycentric"), + ("doppler_relativistic", "topocentric"), + ("doppler_radio", "galactic"), + ("doppler_optical", "heliocentric"), + ] + for i, j in zip(inputs, outputs): + assert decode_veldef(i) == j + + # Now test that bad input raises an exception + try: + decode_veldef("This is more than 8 chars") + except ValueError: + assert True + try: + # frame fails + decode_veldef("OPTI-LRS") + except KeyError: + assert True + try: + # convention fails + decode_veldef("MAXILSRK") + except KeyError: + assert True diff --git a/src/dysh/spectra/tests/test_spectra_core.py b/src/dysh/spectra/tests/test_spectra_core.py new file mode 100644 index 00000000..5882cd64 --- /dev/null +++ b/src/dysh/spectra/tests/test_spectra_core.py @@ -0,0 +1,62 @@ +import os + +import numpy as np +import pytest +from astropy.io import fits + +from dysh import util +from dysh.spectra import core + +LOCALDIR = os.path.dirname(os.path.realpath(__file__)) + + +class TestMeanTsys: + """ + Tests for `dysh.spectra.core.dcmeantsys` function. + """ + + def setup_method(self): + self.root_dir = util.get_project_root() + self.data_dir = f"{self.root_dir}/testdata" + + def test_tsys(self): + expected = np.array([17.24000345, 17.17140405, 17.15663698]) + + path_to_file = f"{self.data_dir}/TGBT21A_501_11" + filename = "TGBT21A_501_11_ifnum_0_int_0-2.fits" + sdf_file = f"{path_to_file}/{filename}" + + # Open and select data. + hdu_sdf = fits.open(sdf_file) + table = hdu_sdf[1].data + table_pl0 = table[table["PLNUM"] == 0] + table_pl0_off = table_pl0[table_pl0["SCAN"] == 153] + tcal = table_pl0_off["TCAL"][0] + tsys_dysh = np.empty(table_pl0_off["DATA"].shape[0] // 2, dtype=float) + for i in range(len(tsys_dysh)): + tsys_dysh[i] = core.mean_tsys( + calon=table_pl0_off["DATA"][1::2][i], caloff=table_pl0_off["DATA"][0::2][i], tcal=tcal + ) + # Compare. + assert tsys_dysh == pytest.approx(expected) + + def test_tsys2(self): + path_to_file = f"{self.data_dir}/TGBT21A_501_11" + filein = f"{path_to_file}/TGBT21A_501_11.raw.vegas.fits" + gbtidl_file = f"{path_to_file}/TGBT21A_501_11_getps_scan_152_intnum_0_ifnum_0_plnum_0.fits" + + hdu = fits.open(filein) + table = hdu[1].data + mask = (table["SCAN"] == 153) & (table["IFNUM"] == 0) & (table["PLNUM"] == 0) + mask_on = table[mask]["CAL"] == "T" + mask_off = table[mask]["CAL"] == "F" + table_on = table[mask][mask_on] + table_off = table[mask][mask_off] + nchan = table["DATA"].shape[1] + tsys_dysh = core.mean_tsys(table_on["DATA"][0], table_off["DATA"][0], table_on["TCAL"][0]) + + hdu = fits.open(gbtidl_file) + gbtidl_table = hdu[1].data + gbtidl_tsys = gbtidl_table["TSYS"] + + assert tsys_dysh == gbtidl_tsys From e90d6fed21374443a579b7dac51a90240e342df5 Mon Sep 17 00:00:00 2001 From: Marc Pound <22331890+mpound@users.noreply.github.com> Date: Fri, 1 Dec 2023 13:23:53 -0500 Subject: [PATCH 3/4] add core --- src/dysh/fits/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/dysh/fits/__init__.py b/src/dysh/fits/__init__.py index b3597907..4b1b1ab4 100644 --- a/src/dysh/fits/__init__.py +++ b/src/dysh/fits/__init__.py @@ -1,4 +1,5 @@ """Classes and functions for importing SDFITS files""" +from dysh.fits.core import * # noqa from dysh.fits.gb20mfitsload import GB20MFITSLoad # noqa from dysh.fits.gbtfitsload import GBTFITSLoad # noqa from dysh.fits.sdfitsload import SDFITSLoad # noqa From f931368b2dab709c82c0d02100c23e5919501905 Mon Sep 17 00:00:00 2001 From: Marc Pound <22331890+mpound@users.noreply.github.com> Date: Mon, 4 Dec 2023 09:46:12 -0500 Subject: [PATCH 4/4] delete old names --- src/dysh/fits/tests/test_core.py | 48 ---------------------- src/dysh/spectra/tests/test_core.py | 62 ----------------------------- 2 files changed, 110 deletions(-) delete mode 100644 src/dysh/fits/tests/test_core.py delete mode 100644 src/dysh/spectra/tests/test_core.py diff --git a/src/dysh/fits/tests/test_core.py b/src/dysh/fits/tests/test_core.py deleted file mode 100644 index bfb9060f..00000000 --- a/src/dysh/fits/tests/test_core.py +++ /dev/null @@ -1,48 +0,0 @@ -from dysh.fits import decode_veldef - - -class TestCore: - """Test dysh.fits core functions""" - - def test_veldef(self): - # first make sure we get correct answers for normal inputs - inputs = [ - "RADILSRK", - "RADI-LSR", - "RADILSRD", - "OPTICMB", - "VELO-BAR", - "OPTIBARY", - "RELATOPO", - "RADIGALA", - "OPTI-HEL", - ] - outputs = [ - ("doppler_radio", "LSRK"), - ("doppler_radio", "LSRK"), - ("doppler_radio", "LSRD"), - ("doppler_optical", "cmb"), - ("doppler_relativistic", "barycentric"), - ("doppler_optical", "barycentric"), - ("doppler_relativistic", "topocentric"), - ("doppler_radio", "galactic"), - ("doppler_optical", "heliocentric"), - ] - for i, j in zip(inputs, outputs): - assert decode_veldef(i) == j - - # Now test that bad input raises an exception - try: - decode_veldef("This is more than 8 chars") - except ValueError: - assert True - try: - # frame fails - decode_veldef("OPTI-LRS") - except KeyError: - assert True - try: - # convention fails - decode_veldef("MAXILSRK") - except KeyError: - assert True diff --git a/src/dysh/spectra/tests/test_core.py b/src/dysh/spectra/tests/test_core.py deleted file mode 100644 index 5882cd64..00000000 --- a/src/dysh/spectra/tests/test_core.py +++ /dev/null @@ -1,62 +0,0 @@ -import os - -import numpy as np -import pytest -from astropy.io import fits - -from dysh import util -from dysh.spectra import core - -LOCALDIR = os.path.dirname(os.path.realpath(__file__)) - - -class TestMeanTsys: - """ - Tests for `dysh.spectra.core.dcmeantsys` function. - """ - - def setup_method(self): - self.root_dir = util.get_project_root() - self.data_dir = f"{self.root_dir}/testdata" - - def test_tsys(self): - expected = np.array([17.24000345, 17.17140405, 17.15663698]) - - path_to_file = f"{self.data_dir}/TGBT21A_501_11" - filename = "TGBT21A_501_11_ifnum_0_int_0-2.fits" - sdf_file = f"{path_to_file}/{filename}" - - # Open and select data. - hdu_sdf = fits.open(sdf_file) - table = hdu_sdf[1].data - table_pl0 = table[table["PLNUM"] == 0] - table_pl0_off = table_pl0[table_pl0["SCAN"] == 153] - tcal = table_pl0_off["TCAL"][0] - tsys_dysh = np.empty(table_pl0_off["DATA"].shape[0] // 2, dtype=float) - for i in range(len(tsys_dysh)): - tsys_dysh[i] = core.mean_tsys( - calon=table_pl0_off["DATA"][1::2][i], caloff=table_pl0_off["DATA"][0::2][i], tcal=tcal - ) - # Compare. - assert tsys_dysh == pytest.approx(expected) - - def test_tsys2(self): - path_to_file = f"{self.data_dir}/TGBT21A_501_11" - filein = f"{path_to_file}/TGBT21A_501_11.raw.vegas.fits" - gbtidl_file = f"{path_to_file}/TGBT21A_501_11_getps_scan_152_intnum_0_ifnum_0_plnum_0.fits" - - hdu = fits.open(filein) - table = hdu[1].data - mask = (table["SCAN"] == 153) & (table["IFNUM"] == 0) & (table["PLNUM"] == 0) - mask_on = table[mask]["CAL"] == "T" - mask_off = table[mask]["CAL"] == "F" - table_on = table[mask][mask_on] - table_off = table[mask][mask_off] - nchan = table["DATA"].shape[1] - tsys_dysh = core.mean_tsys(table_on["DATA"][0], table_off["DATA"][0], table_on["TCAL"][0]) - - hdu = fits.open(gbtidl_file) - gbtidl_table = hdu[1].data - gbtidl_tsys = gbtidl_table["TSYS"] - - assert tsys_dysh == gbtidl_tsys