From 193c6ca2d3b09445b5196e54554957a9af2ec31c Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Wed, 2 Dec 2020 19:36:58 -0300 Subject: [PATCH 01/20] Set NDData.unit when reading FITS files --- astrodata/fits.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/astrodata/fits.py b/astrodata/fits.py index 5116fd61fd..55471d47c0 100644 --- a/astrodata/fits.py +++ b/astrodata/fits.py @@ -30,6 +30,12 @@ NO_DEFAULT = object() LOGGER = logging.getLogger(__name__) +known_invalid_fits_unit_strings = { + 'ELECTRONS/S': u.electron/u.s, + 'ELECTRONS': u.electron, + 'electrons': u.electron, +} + class FitsHeaderCollection: """Group access to a list of FITS Header-like objects. @@ -391,7 +397,7 @@ def _prepare_hdulist(hdulist, default_extension='SCI', extname_parser=None): return HDUList(sorted(new_list, key=fits_ext_comp_key)) -def read_fits(cls, source, extname_parser=None): +def read_fits(cls, source, extname_parser=None, default_unit='adu'): """ Takes either a string (with the path to a file) or an HDUList as input, and tries to return a populated AstroData (or descendant) instance. @@ -477,12 +483,24 @@ def associated_extensions(ver): not isinstance(parts['uncertainty'], FitsLazyLoadable)): parts['uncertainty'] = ADVarianceUncertainty(parts['uncertainty']) + bunit = header.get('BUNIT') + if bunit: + if bunit in known_invalid_fits_unit_strings: + bunit = known_invalid_fits_unit_strings[bunit] + try: + unit = u.Unit(bunit, format='fits') + except ValueError: + unit = u.Unit(default_unit) + else: + unit = u.Unit(default_unit) + # Create the NDData object nd = NDDataObject( data=parts['data'], uncertainty=parts['uncertainty'], mask=parts['mask'], meta={'header': header}, + unit=unit, ) if parts['wcs'] is not None: From 7bdfee601fecc0bc96ef0d595de3f624586f92ab Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Mon, 18 Jan 2021 19:24:35 -0300 Subject: [PATCH 02/20] Add setter for unit --- astrodata/nddata.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/astrodata/nddata.py b/astrodata/nddata.py index c3ac1d25aa..e76e892b93 100644 --- a/astrodata/nddata.py +++ b/astrodata/nddata.py @@ -8,8 +8,8 @@ from copy import deepcopy from functools import reduce +import astropy.units as u import numpy as np - from astropy.io.fits import ImageHDU from astropy.modeling import Model, models from astropy.nddata import (NDArithmeticMixin, NDData, NDSlicingMixin, @@ -369,6 +369,18 @@ def variance(self, value): self.uncertainty = (ADVarianceUncertainty(value) if value is not None else None) + # Override unit so that we can add a setter. + @property + def unit(self): + return self._unit + + @unit.setter + def unit(self, value): + if value is None: + self._unit = None + else: + self._unit = u.Unit(value) + def set_section(self, section, input): """ Sets only a section of the data. This method is meant to prevent From 4c5ea655c749bf1456a9d52860090caad96a667a Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Mon, 18 Jan 2021 20:08:53 -0300 Subject: [PATCH 03/20] Write unit to BUNIT --- astrodata/fits.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/astrodata/fits.py b/astrodata/fits.py index 55471d47c0..b31d4d1e5d 100644 --- a/astrodata/fits.py +++ b/astrodata/fits.py @@ -588,9 +588,16 @@ def ad_to_hdulist(ad): if 'APPROXIMATE' not in wcs_dict.get('FITS-WCS', ''): wcs = None # There's no need to create a WCS extension - hdul.append(new_imagehdu(ext.data, header, 'SCI')) + data_hdu = new_imagehdu(ext.data, header, 'SCI') + + if ext.unit is not None and ext.unit is not u.dimensionless_unscaled: + data_hdu.header['BUNIT'] = ext.unit.to_string() + + hdul.append(data_hdu) + if ext.uncertainty is not None: hdul.append(new_imagehdu(ext.uncertainty.array, header, 'VAR')) + if ext.mask is not None: hdul.append(new_imagehdu(ext.mask, header, 'DQ')) From b430235cac88890f77ae5cd18aa7c75a331a41ed Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Mon, 15 Mar 2021 18:05:37 -0300 Subject: [PATCH 04/20] Fix getting unit from AstroData for Spek1D --- gempy/library/spectral.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/gempy/library/spectral.py b/gempy/library/spectral.py index 362a494f89..34efc99334 100644 --- a/gempy/library/spectral.py +++ b/gempy/library/spectral.py @@ -42,12 +42,16 @@ def __init__(self, spectrum=None, spectral_axis=None, wcs=None, **kwargs): # Unit handling try: # for NDData-like - flux_unit = spectrum.unit + if isinstance(spectrum, AstroData): + flux_unit = spectrum.nddata.unit + else: + flux_unit = spectrum.unit except AttributeError: try: # for AstroData flux_unit = u.Unit(spectrum.hdr.get('BUNIT')) except (TypeError, ValueError): # unknown/missing flux_unit = None + if flux_unit is None: flux_unit = u.dimensionless_unscaled try: From dcf3ae2efa4e8ee33a69ed49acdc3de5d4d16a31 Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Mon, 15 Mar 2021 18:10:35 -0300 Subject: [PATCH 05/20] Update unit in ADUToElectrons --- geminidr/core/primitives_preprocess.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/geminidr/core/primitives_preprocess.py b/geminidr/core/primitives_preprocess.py index f6a25b849d..0bfd85822e 100644 --- a/geminidr/core/primitives_preprocess.py +++ b/geminidr/core/primitives_preprocess.py @@ -9,11 +9,11 @@ from copy import deepcopy from functools import partial +import astropy.units as u import astrodata import gemini_instruments # noqa import matplotlib.pyplot as plt import numpy as np -from astrodata import NDAstroData from astrodata.provenance import add_provenance from astropy.table import Table from geminidr import PrimitivesBASE @@ -107,7 +107,7 @@ def ADUToElectrons(self, adinputs=None, suffix=None): "the gain".format(ad.filename)) for ext, gain in zip(ad, gain_list): log.stdinfo(f" gain for extension {ext.id} = {gain}") - ext.multiply(gain) + ext.multiply(gain * u.electron / u.adu) # Update the headers of the AstroData Object. The pixel data now # has units of electrons so update the physical units keyword. From 90b96889b9ecf8f79379abc7258e351c48a963df Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Mon, 15 Mar 2021 18:50:48 -0300 Subject: [PATCH 06/20] Set default unit for .append and fix unit tests --- astrodata/core.py | 8 +++++++- astrodata/tests/test_core.py | 32 +++++++++++++++----------------- 2 files changed, 22 insertions(+), 18 deletions(-) diff --git a/astrodata/core.py b/astrodata/core.py index d55b538118..a87423f406 100644 --- a/astrodata/core.py +++ b/astrodata/core.py @@ -891,7 +891,7 @@ def _append_raw_nddata(self, raw_nddata, name, header, add_to): custom_header=header) return self._append_nddata(processed_nddata, name=name, add_to=add_to) - def _append_nddata(self, new_nddata, name, add_to): + def _append_nddata(self, new_nddata, name, add_to, unit='adu'): # NOTE: This method is only used by others that have constructed NDData # according to our internal format. We don't accept new headers at this # point, and that's why it's missing from the signature. 'name' is @@ -900,6 +900,12 @@ def _append_nddata(self, new_nddata, name, add_to): raise TypeError("You can only append NDData derived instances " "at the top level") + if new_nddata.unit is None: + # FIXME: setting unit to adu by default but we should allow + # passing a unit in .append() + new_nddata.unit = 'adu' + warnings.warn('setting unit to ADU') + hd = new_nddata.meta['header'] hname = hd.get('EXTNAME', DEFAULT_EXTENSION) if hname == DEFAULT_EXTENSION: diff --git a/astrodata/tests/test_core.py b/astrodata/tests/test_core.py index ab80e42fcd..537049f65a 100644 --- a/astrodata/tests/test_core.py +++ b/astrodata/tests/test_core.py @@ -6,6 +6,7 @@ from numpy.testing import assert_array_equal import astrodata +import astropy.units as u from astropy.io import fits from astropy.nddata import NDData, VarianceUncertainty from astropy.table import Table @@ -48,7 +49,7 @@ def test_attributes(ad1): (operator.truediv, 0.5, 2) ]) def test_arithmetic(op, res, res2, ad1, ad2): - for data in (ad2, ad2.data): + for data in (ad2, ad2.data * u.adu): result = op(ad1, data) assert_array_equal(result.data, res) assert isinstance(result, astrodata.AstroData) @@ -56,10 +57,10 @@ def test_arithmetic(op, res, res2, ad1, ad2): assert isinstance(result[0].data, np.ndarray) assert isinstance(result[0].hdr, fits.Header) - result = op(data, ad1) - assert_array_equal(result.data, res2) + result = op(ad2, ad1) + assert_array_equal(result.data, res2) - for data in (ad2[0], ad2[0].data): + for data in (ad2[0], ad2[0].data * u.adu): result = op(ad1[0], data) assert_array_equal(result.data, res) assert isinstance(result, astrodata.AstroData) @@ -67,7 +68,6 @@ def test_arithmetic(op, res, res2, ad1, ad2): assert isinstance(result[0].data, np.ndarray) assert isinstance(result[0].hdr, fits.Header) - # FIXME: should work also with ad2[0].data, but crashes result = op(ad2[0], ad1[0]) assert_array_equal(result.data, res2) @@ -79,7 +79,7 @@ def test_arithmetic(op, res, res2, ad1, ad2): (operator.itruediv, 0.5, 2) ]) def test_arithmetic_inplace(op, res, res2, ad1, ad2): - for data in (ad2, ad2.data): + for data in (ad2, ad2.data * u.adu): ad = deepcopy(ad1) op(ad, data) assert_array_equal(ad.data, res) @@ -88,11 +88,7 @@ def test_arithmetic_inplace(op, res, res2, ad1, ad2): assert isinstance(ad[0].data, np.ndarray) assert isinstance(ad[0].hdr, fits.Header) - # data2 = deepcopy(ad2[0]) - # op(data2, ad1) - # assert_array_equal(data2, res2) - - for data in (ad2[0], ad2[0].data): + for data in (ad2[0], ad2[0].data * u.adu): ad = deepcopy(ad1) op(ad[0], data) assert_array_equal(ad.data, res) @@ -111,13 +107,13 @@ def test_arithmetic_inplace(op, res, res2, ad1, ad2): def test_arithmetic_multiple_ext(op, res, ad1): ad1.append(np.ones(SHAPE, dtype=np.uint16) + 4) - result = op(ad1, 2) + result = op(ad1, 2 * u.adu) assert len(result) == 2 assert_array_equal(result[0].data, res[0]) assert_array_equal(result[1].data, res[1]) for i, ext in enumerate(ad1): - result = op(ext, 2) + result = op(ext, 2 * u.adu) assert len(result) == 1 assert_array_equal(result[0].data, res[i]) @@ -132,7 +128,7 @@ def test_arithmetic_inplace_multiple_ext(op, res, ad1): ad1.append(np.ones(SHAPE, dtype=np.uint16) + 4) ad = deepcopy(ad1) - result = op(ad, 2) + result = op(ad, 2 * u.adu) assert len(result) == 2 assert_array_equal(result[0].data, res[0]) assert_array_equal(result[1].data, res[1]) @@ -141,13 +137,13 @@ def test_arithmetic_inplace_multiple_ext(op, res, ad1): # as it is now independant from its parent for i, ext in enumerate(ad1): ext = deepcopy(ext) - result = op(ext, 2) + result = op(ext, 2 * u.adu) assert len(result) == 1 assert_array_equal(result[0].data, res[i]) # Now work directly on the input object, will creates single ad objects for i, ext in enumerate(ad1): - result = op(ext, 2) + result = op(ext, 2 * u.adu) assert len(result) == 1 assert_array_equal(result.data, res[i]) @@ -157,7 +153,7 @@ def test_arithmetic_inplace_multiple_ext(op, res, ad1): ('multiply', 3, 3), ('divide', 2, 0.5)]) def test_operations(ad1, op, arg, res): - result = getattr(ad1, op)(arg) + result = getattr(ad1, op)(arg * u.adu) assert_array_equal(result.data, res) assert isinstance(result, astrodata.AstroData) assert isinstance(result[0].data, np.ndarray) @@ -170,6 +166,7 @@ def test_operate(): nd = NDData(data=[[1, 2], [3, 4]], uncertainty=VarianceUncertainty(np.ones((2, 2))), mask=np.identity(2), + unit='adu', meta={'header': fits.Header()}) ad.append(nd) @@ -184,6 +181,7 @@ def test_write_and_read(tmpdir, capsys): nd = NDData(data=[[1, 2], [3, 4]], uncertainty=VarianceUncertainty(np.ones((2, 2))), mask=np.identity(2), + unit='adu', meta={'header': fits.Header()}) ad.append(nd) From 16cbd3472505bdbf0bbbb46afe3ed2cad9c8a4d7 Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Mon, 15 Mar 2021 19:21:35 -0300 Subject: [PATCH 07/20] Fix test_stack --- geminidr/core/tests/test_stack.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/geminidr/core/tests/test_stack.py b/geminidr/core/tests/test_stack.py index 8561e1ea9b..7d8828cc49 100644 --- a/geminidr/core/tests/test_stack.py +++ b/geminidr/core/tests/test_stack.py @@ -4,6 +4,7 @@ import numpy as np import pytest +import astropy.units as u from astropy.io import fits from astropy.table import Table from numpy.testing import assert_array_equal @@ -90,7 +91,7 @@ def test_stackframes_refcat_propagation(niri_adinputs): def test_rejmap(niri_adinputs): for i in (2, 3, 4): - niri_adinputs.append(niri_adinputs[0] + i) + niri_adinputs.append(niri_adinputs[0] + i * u.adu) p = NIRIImage(niri_adinputs) p.prepare() From 7ec18cd814ff0b72d7b1911feeffbc78830975ac Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Mon, 15 Mar 2021 19:25:24 -0300 Subject: [PATCH 08/20] Fix test_spectral --- gempy/library/tests/test_spectral.py | 1 + 1 file changed, 1 insertion(+) diff --git a/gempy/library/tests/test_spectral.py b/gempy/library/tests/test_spectral.py index a26128739d..e30948d366 100644 --- a/gempy/library/tests/test_spectral.py +++ b/gempy/library/tests/test_spectral.py @@ -33,6 +33,7 @@ def _create_fake_data(): _ad.add_extension(hdu, pixel_scale=1.0) _ad[0].wcs = None # or else imaging WCS will be added + _ad[0].nddata.unit = units.electron _ad[0].data = _ad[0].data.ravel() + 1. _ad[0].mask = np.zeros(_ad[0].data.size, dtype=np.uint16) # ToDo Requires mask _ad[0].variance = np.ones_like(_ad[0].data) # ToDo Requires Variance From e2f48c69551b72525093e20808ceecc122941e3b Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Mon, 15 Mar 2021 19:29:37 -0300 Subject: [PATCH 09/20] Use AstroFaker branch --- tox.ini | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tox.ini b/tox.ini index 3b8a5a7636..753413a3ce 100644 --- a/tox.ini +++ b/tox.ini @@ -47,7 +47,7 @@ commands = which pytest pip freeze -l pip install git+https://github.com/astropy/astroscrappy.git#egg=astroscrappy - pip install git+https://github.com/GeminiDRSoftware/AstroFaker#egg=AstroFaker + pip install git+https://github.com/GeminiDRSoftware/AstroFaker@handle-units#egg=AstroFaker unit: python -m coverage run -m pytest -v --dragons-remote-data --durations=20 -m "not integration_test and not gmosls and not regression and not slow" {posargs:astrodata geminidr gemini_instruments gempy recipe_system} integ: python -m coverage run -m pytest -v --dragons-remote-data --durations=20 -m "integration_test and not slow" {posargs:astrodata geminidr gemini_instruments gempy recipe_system} gmosls: python -m coverage run -m pytest -v --dragons-remote-data --durations=20 -m "gmosls and not slow" {posargs:astrodata geminidr gemini_instruments gempy recipe_system} From 513d0ad9ff1f85c8ae9320ce24b2d8334e40a18e Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Tue, 16 Mar 2021 18:36:20 -0300 Subject: [PATCH 10/20] Add unit argument for .append --- astrodata/core.py | 45 ++++++++++++--------- astrodata/fits.py | 1 + astrodata/tests/test_object_construction.py | 26 ++++++++++-- 3 files changed, 48 insertions(+), 24 deletions(-) diff --git a/astrodata/core.py b/astrodata/core.py index a87423f406..9516ec1767 100644 --- a/astrodata/core.py +++ b/astrodata/core.py @@ -849,7 +849,8 @@ def _process_pixel_plane(self, pixim, name=None, top_level=False, return nd - def _append_array(self, data, name=None, header=None, add_to=None): + def _append_array(self, data, name=None, header=None, add_to=None, + unit=None): if name in {'DQ', 'VAR'}: raise ValueError(f"'{name}' need to be associated to a " f"'{DEFAULT_EXTENSION}' one") @@ -866,21 +867,22 @@ def _append_array(self, data, name=None, header=None, add_to=None): hdu = fits.ImageHDU(data, header=header) hdu.header['EXTNAME'] = hname ret = self._append_imagehdu(hdu, name=hname, header=None, - add_to=None) + add_to=None, unit=unit) else: ret = add_to.meta['other'][name] = data return ret - def _append_imagehdu(self, hdu, name, header, add_to): + def _append_imagehdu(self, hdu, name, header, add_to, unit): if name in {'DQ', 'VAR'} or add_to is not None: - return self._append_array(hdu.data, name=name, add_to=add_to) + return self._append_array(hdu.data, name=name, add_to=add_to, + unit=unit) else: nd = self._process_pixel_plane(hdu, name=name, top_level=True, custom_header=header) - return self._append_nddata(nd, name, add_to=None) + return self._append_nddata(nd, name, add_to=None, unit=unit) - def _append_raw_nddata(self, raw_nddata, name, header, add_to): + def _append_raw_nddata(self, raw_nddata, name, header, add_to, unit): # We want to make sure that the instance we add is whatever we specify # as NDDataObject, instead of the random one that the user may pass top_level = add_to is None @@ -889,9 +891,10 @@ def _append_raw_nddata(self, raw_nddata, name, header, add_to): processed_nddata = self._process_pixel_plane(raw_nddata, top_level=top_level, custom_header=header) - return self._append_nddata(processed_nddata, name=name, add_to=add_to) + return self._append_nddata(processed_nddata, name=name, add_to=add_to, + unit=unit) - def _append_nddata(self, new_nddata, name, add_to, unit='adu'): + def _append_nddata(self, new_nddata, name, add_to, unit): # NOTE: This method is only used by others that have constructed NDData # according to our internal format. We don't accept new headers at this # point, and that's why it's missing from the signature. 'name' is @@ -901,10 +904,7 @@ def _append_nddata(self, new_nddata, name, add_to, unit='adu'): "at the top level") if new_nddata.unit is None: - # FIXME: setting unit to adu by default but we should allow - # passing a unit in .append() - new_nddata.unit = 'adu' - warnings.warn('setting unit to ADU') + new_nddata.unit = unit hd = new_nddata.meta['header'] hname = hd.get('EXTNAME', DEFAULT_EXTENSION) @@ -916,7 +916,7 @@ def _append_nddata(self, new_nddata, name, add_to, unit='adu'): return new_nddata - def _append_table(self, new_table, name, header, add_to): + def _append_table(self, new_table, name, header, add_to, unit): tb = _process_table(new_table, name, header) hname = tb.meta['header'].get('EXTNAME') @@ -948,7 +948,7 @@ def find_next_num(tables): add_to.meta['other'][hname] = tb return tb - def _append_astrodata(self, ad, name, header, add_to): + def _append_astrodata(self, ad, name, header, add_to, unit): if not ad.is_single: raise ValueError("Cannot append AstroData instances that are " "not single slices") @@ -960,9 +960,10 @@ def _append_astrodata(self, ad, name, header, add_to): if header is not None: new_nddata.meta['header'] = deepcopy(header) - return self._append_nddata(new_nddata, name=None, add_to=None) + return self._append_nddata(new_nddata, name=None, add_to=None, + unit=unit) - def _append(self, ext, name=None, header=None, add_to=None): + def _append(self, ext, name=None, header=None, add_to=None, unit='adu'): """ Internal method to dispatch to the type specific methods. This is called either by ``.append`` to append on top-level objects only or @@ -978,12 +979,14 @@ def _append(self, ext, name=None, header=None, add_to=None): for bases, method in dispatcher: if isinstance(ext, bases): - return method(ext, name=name, header=header, add_to=add_to) + return method(ext, name=name, header=header, add_to=add_to, + unit=unit) # Assume that this is an array for a pixel plane - return self._append_array(ext, name=name, header=header, add_to=add_to) + return self._append_array(ext, name=name, header=header, add_to=add_to, + unit=unit) - def append(self, ext, name=None, header=None): + def append(self, ext, name=None, header=None, unit='adu'): """ Adds a new top-level extension. @@ -1004,6 +1007,8 @@ def append(self, ext, name=None, header=None): It can consist in a combination of numbers and letters, with the restriction that the letters have to be all capital, and the first character cannot be a number ("[A-Z][A-Z0-9]*"). + header : `astropy.io.fits.Header` + FITS Header to be associated with the NDData or Table object. Returns -------- @@ -1038,7 +1043,7 @@ def append(self, ext, name=None, header=None): UserWarning) name = name.upper() - return self._append(ext, name=name, header=header) + return self._append(ext, name=name, header=header, unit=unit) @classmethod def read(cls, source, extname_parser=None): diff --git a/astrodata/fits.py b/astrodata/fits.py index b31d4d1e5d..c4fadb1a00 100644 --- a/astrodata/fits.py +++ b/astrodata/fits.py @@ -34,6 +34,7 @@ 'ELECTRONS/S': u.electron/u.s, 'ELECTRONS': u.electron, 'electrons': u.electron, + 'electron': u.electron, } diff --git a/astrodata/tests/test_object_construction.py b/astrodata/tests/test_object_construction.py index 3b464173a1..0c3a757601 100644 --- a/astrodata/tests/test_object_construction.py +++ b/astrodata/tests/test_object_construction.py @@ -94,16 +94,19 @@ def test_create_invalid(): def test_append_image_hdu(): + shape = (4, 5) ad = astrodata.create(fits.PrimaryHDU()) - ad.append(fits.ImageHDU(data=np.zeros((4, 5)))) - ad.append(fits.ImageHDU(data=np.zeros((4, 5))), name='SCI') + ad.append(fits.ImageHDU(data=np.zeros(shape))) + ad.append(fits.ImageHDU(data=np.zeros(shape)), name='SCI', unit='electron') with pytest.raises(ValueError, match="Arbitrary image extensions can only be added " "in association to a 'SCI'"): - ad.append(fits.ImageHDU(data=np.zeros((4, 5))), name='SCI2') + ad.append(fits.ImageHDU(data=np.zeros(shape)), name='SCI2') assert len(ad) == 2 + assert ad[0].nddata.unit == 'adu' + assert ad[1].nddata.unit == 'electron' def test_append_lowercase_name(): @@ -113,12 +116,25 @@ def test_append_lowercase_name(): ad.append(NDData(np.zeros((4, 5))), name='sci') +def test_append_nddata(): + shape = (4, 5) + ad = astrodata.create({}) + ad.append(NDData(np.zeros(shape))) + ad.append(NDData(np.zeros(shape), unit='electron')) + ad.append(NDData(np.zeros(shape)), unit='electron') + + assert ad[0].nddata.unit == 'adu' + assert ad[1].nddata.unit == 'electron' + assert ad[2].nddata.unit == 'electron' + + def test_append_arrays(tmp_path): testfile = tmp_path / 'test.fits' ad = astrodata.create({}) ad.append(np.zeros(10)) ad[0].ARR = np.arange(5) + ad.append(np.zeros(10), unit='electron') with pytest.raises(AttributeError): ad[0].SCI = np.arange(5) @@ -138,7 +154,9 @@ def test_append_arrays(tmp_path): ad.write(testfile) ad = astrodata.open(testfile) - assert len(ad) == 1 + assert len(ad) == 2 + assert ad[0].nddata.unit == 'adu' + assert ad[1].nddata.unit == 'electron' assert ad[0].nddata.meta['header']['EXTNAME'] == 'SCI' assert_array_equal(ad[0].ARR, np.arange(5)) From a146920a2f22efc29ff2eefcc94dac5daee34f4f Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Tue, 16 Mar 2021 18:41:26 -0300 Subject: [PATCH 11/20] Revert to use master branch of AstroFaker --- tox.ini | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tox.ini b/tox.ini index 753413a3ce..3b8a5a7636 100644 --- a/tox.ini +++ b/tox.ini @@ -47,7 +47,7 @@ commands = which pytest pip freeze -l pip install git+https://github.com/astropy/astroscrappy.git#egg=astroscrappy - pip install git+https://github.com/GeminiDRSoftware/AstroFaker@handle-units#egg=AstroFaker + pip install git+https://github.com/GeminiDRSoftware/AstroFaker#egg=AstroFaker unit: python -m coverage run -m pytest -v --dragons-remote-data --durations=20 -m "not integration_test and not gmosls and not regression and not slow" {posargs:astrodata geminidr gemini_instruments gempy recipe_system} integ: python -m coverage run -m pytest -v --dragons-remote-data --durations=20 -m "integration_test and not slow" {posargs:astrodata geminidr gemini_instruments gempy recipe_system} gmosls: python -m coverage run -m pytest -v --dragons-remote-data --durations=20 -m "gmosls and not slow" {posargs:astrodata geminidr gemini_instruments gempy recipe_system} From d9797d42da1056af5daa51c04fe3663fd729b128 Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Tue, 16 Mar 2021 18:58:52 -0300 Subject: [PATCH 12/20] Add test for ADUToElectrons --- geminidr/core/tests/test_preprocess.py | 31 +++++++++++++++++--------- 1 file changed, 21 insertions(+), 10 deletions(-) diff --git a/geminidr/core/tests/test_preprocess.py b/geminidr/core/tests/test_preprocess.py index d51d448b27..0f9ee6b7f6 100644 --- a/geminidr/core/tests/test_preprocess.py +++ b/geminidr/core/tests/test_preprocess.py @@ -309,6 +309,27 @@ def test_fixpixels_multiple_ext(niriprim2): assert_almost_equal(ad[1].data[sy, sx].max(), 60.333, decimal=2) +def test_adu_to_electrons(astrofaker, caplog): + caplog.set_level('INFO') + ad = astrofaker.create("NIRI", "IMAGE") + ad.init_default_extensions() + ad[0].data[:] = 1 + p = NIRIImage([ad]) + + ad = p.ADUToElectrons()[0] + assert ad[0].nddata.unit == 'electron' + assert_array_almost_equal(ad[0].data, ad[0].gain()) + assert caplog.messages[2] == ('Converting N20010101S0001.fits from ADU to ' + 'electrons by multiplying by the gain') + caplog.clear() + + ad = p.ADUToElectrons()[0] + assert ad[0].nddata.unit == 'electron' + assert_array_almost_equal(ad[0].data, ad[0].gain()) + assert ('No changes will be made to N20010101S0001_ADUToElectrons.fits, ' + 'since it has already been processed') in caplog.messages[2] + + # TODO @bquint: clean up these tests # @pytest.fixture @@ -357,16 +378,6 @@ def test_fixpixels_multiple_ext(niriprim2): # assert all(ext.mask[ext.OBJMASK == 1] == ext_orig.mask[ext.OBJMASK == 1] | 1) -# @pytest.mark.xfail(reason="Test needs revision", run=False) -# def test_adu_to_electrons(astrofaker): -# ad = astrofaker.create("NIRI", "IMAGE") -# # astrodata.open(os.path.join(TESTDATAPATH, 'NIRI', 'N20070819S0104_dqAdded.fits')) -# p = NIRIImage([ad]) -# ad = p.ADUToElectrons()[0] -# assert ad_compare(ad, os.path.join(TESTDATAPATH, 'NIRI', -# 'N20070819S0104_ADUToElectrons.fits')) - - # @pytest.mark.xfail(reason="Test needs revision", run=False) # def test_associateSky(): # filenames = ['N20070819S{:04d}_flatCorrected.fits'.format(i) From ffd90ec3b0b8f557aa4da021b7376f1de2578fb1 Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Tue, 16 Mar 2021 19:29:17 -0300 Subject: [PATCH 13/20] Add .unit attribute for AstroData --- astrodata/core.py | 14 ++++++++++++++ astrodata/tests/test_object_construction.py | 11 +++++++---- 2 files changed, 21 insertions(+), 4 deletions(-) diff --git a/astrodata/core.py b/astrodata/core.py index 9516ec1767..f3be858ce7 100644 --- a/astrodata/core.py +++ b/astrodata/core.py @@ -444,6 +444,20 @@ def wcs(self): def wcs(self, value): self.nddata.wcs = value + @property + @returns_list + def unit(self): + """ + A list of `astropy.units.Unit` objects (or a single object if this + is a single slice) attached to the science data, for each extension. + """ + return [nd.unit for nd in self._nddata] + + @unit.setter + @assign_only_single_slice + def unit(self, value): + self.nddata.unit = value + def __iter__(self): if self.is_single: yield self diff --git a/astrodata/tests/test_object_construction.py b/astrodata/tests/test_object_construction.py index 0c3a757601..2903cae86b 100644 --- a/astrodata/tests/test_object_construction.py +++ b/astrodata/tests/test_object_construction.py @@ -116,16 +116,19 @@ def test_append_lowercase_name(): ad.append(NDData(np.zeros((4, 5))), name='sci') -def test_append_nddata(): +def test_append_nddata_and_units(): shape = (4, 5) ad = astrodata.create({}) ad.append(NDData(np.zeros(shape))) ad.append(NDData(np.zeros(shape), unit='electron')) ad.append(NDData(np.zeros(shape)), unit='electron') - assert ad[0].nddata.unit == 'adu' - assert ad[1].nddata.unit == 'electron' - assert ad[2].nddata.unit == 'electron' + assert ad[0].unit == 'adu' + assert ad[1].unit == 'electron' + assert ad[2].unit == 'electron' + assert ad.unit == ['adu', 'electron', 'electron'] + ad[1].unit = 'adu' + assert ad.unit == ['adu', 'adu', 'electron'] def test_append_arrays(tmp_path): From 6a1c010e99eae932d1ecc76b9f19d2a4129e846a Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Tue, 16 Mar 2021 19:44:07 -0300 Subject: [PATCH 14/20] Use .unit where possible --- astrodata/fits.py | 1 + geminidr/core/primitives_spect.py | 32 ++++++++----------- geminidr/core/primitives_visualize.py | 4 +-- geminidr/core/tests/test_preprocess.py | 4 +-- geminidr/core/tests/test_spect.py | 2 +- .../tests/spect/test_calculate_sensitivity.py | 4 +-- .../gmos/tests/spect/test_flux_calibration.py | 5 ++- gempy/library/spectral.py | 16 ++-------- gempy/library/tests/test_spectral.py | 3 +- 9 files changed, 27 insertions(+), 44 deletions(-) diff --git a/astrodata/fits.py b/astrodata/fits.py index c4fadb1a00..4d8783974a 100644 --- a/astrodata/fits.py +++ b/astrodata/fits.py @@ -33,6 +33,7 @@ known_invalid_fits_unit_strings = { 'ELECTRONS/S': u.electron/u.s, 'ELECTRONS': u.electron, + 'ELECTRON': u.electron, 'electrons': u.electron, 'electron': u.electron, } diff --git a/geminidr/core/primitives_spect.py b/geminidr/core/primitives_spect.py index 201ae73ba6..e010501be5 100644 --- a/geminidr/core/primitives_spect.py +++ b/geminidr/core/primitives_spect.py @@ -2092,7 +2092,6 @@ def fluxCalibrate(self, adinputs=None, **params): log.stdinfo(f"{ad.filename}: Correcting for airmass of " f"{delta_airmass:5.3f}") - for index, ext in enumerate(ad): ext_std = std[min(index, len_std-1)] extname = f"{ad.filename} extension {ext.id}" @@ -2107,20 +2106,17 @@ def fluxCalibrate(self, adinputs=None, **params): std_physical_unit = (std_flux_unit.physical_unit if isinstance(std_flux_unit, u.LogUnit) else std_flux_unit) - try: - sci_flux_unit = u.Unit(ext.hdr.get('BUNIT')) - except: - sci_flux_unit = None - if not (std_physical_unit is None or sci_flux_unit is None): - unit = sci_flux_unit * std_physical_unit / flux_units + + if not (std_physical_unit is None or ext.unit is None): + unit = ext.unit * std_physical_unit / flux_units if unit.is_equivalent(u.s): log.fullinfo("Dividing {} by exposure time of {} s". format(extname, exptime)) ext /= exptime - sci_flux_unit /= u.s + ext.unit /= u.s elif not unit.is_equivalent(u.dimensionless_unscaled): log.warning(f"{extname} has incompatible units ('" - f"{sci_flux_unit}' and '{std_physical_unit}'" + f"{ext.unit}' and '{std_physical_unit}'" "). Cannot flux calibrate") continue else: @@ -2163,14 +2159,15 @@ def fluxCalibrate(self, adinputs=None, **params): else: sens_factor *= 10**(0.4 * delta_airmass * extinction_correction) - final_sens_factor = (sci_flux_unit * sens_factor / pixel_sizes).to( + final_sens_factor = (ext.unit * sens_factor / pixel_sizes).to( final_units, equivalencies=u.spectral_density(waves)).value if ndim == 2 and dispaxis == 0: ext *= final_sens_factor[:, np.newaxis] else: ext *= final_sens_factor - ext.hdr['BUNIT'] = final_units + + ext.unit = final_units # Timestamp and update the filename gt.mark_history(ad, primname=self.myself(), keyword=timestamp_key) @@ -3335,8 +3332,8 @@ def conserve_or_interpolate(ext, user_conserve=None, flux_calibrated=False, bool : whether or not to conserve the flux """ ext_str = f"{ext.filename} extension {ext.id}" - ext_unit = ext.hdr["BUNIT"] - if ext_unit in (None, ""): + + if ext.unit is u.dimensionless_unscaled: if user_conserve is None: this_conserve = not flux_calibrated log.stdinfo(f"{ext_str} has no units but " @@ -3349,13 +3346,12 @@ def conserve_or_interpolate(ext, user_conserve=None, flux_calibrated=False, f"been flux calibrated but conserve={user_conserve}") return this_conserve - ext_unit = u.Unit(ext_unit) # Test for units like flux density units_imply_conserve = True for unit1 in ("W", "photon", "electron", "adu"): for unit2 in ("m2", ""): try: - ext_unit.to(u.Unit(f"{unit1} / ({unit2} nm)"), + ext.unit.to(u.Unit(f"{unit1} / ({unit2} nm)"), equivalencies=u.spectral_density(1. * u.m)) except u.UnitConversionError: pass @@ -3365,15 +3361,15 @@ def conserve_or_interpolate(ext, user_conserve=None, flux_calibrated=False, if flux_calibrated and units_imply_conserve: log.warning(f"Possible unit mismatch for {ext_str}. File has been " - f"flux calibrated but units are {ext_unit}") + f"flux calibrated but units are {ext.unit}") if user_conserve is None: this_conserve = units_imply_conserve log.stdinfo(f"Setting conserve={this_conserve} for {ext_str} since " - f"units are {ext_unit}") + f"units are {ext.unit}") else: if user_conserve != units_imply_conserve: log.warning(f"conserve is set to {user_conserve} but the " - f"units of {ext_str} are {ext_unit}") + f"units of {ext_str} are {ext.unit}") this_conserve = user_conserve # but do what we're told return this_conserve diff --git a/geminidr/core/primitives_visualize.py b/geminidr/core/primitives_visualize.py index 122179f164..723a330e61 100644 --- a/geminidr/core/primitives_visualize.py +++ b/geminidr/core/primitives_visualize.py @@ -605,8 +605,6 @@ def plotSpectraForQA(self, adinputs=None, **params): [float(w), float(s)] for w, s in zip(wavelength, stddev)] - _units = ext.hdr["BUNIT"] - center = np.round(ext.hdr["XTRACTED"]) lower = np.round(ext.hdr["XTRACTLO"]) upper = np.round(ext.hdr["XTRACTHI"]) @@ -619,7 +617,7 @@ def plotSpectraForQA(self, adinputs=None, **params): "wavelength_units": w_units, "id": np.round(center + offset), "intensity": _intensity, - "intensity_units": _units, + "intensity_units": ext.unit.to_string(), "slices": _slices, "stddev": _stddev, } diff --git a/geminidr/core/tests/test_preprocess.py b/geminidr/core/tests/test_preprocess.py index 0f9ee6b7f6..0aa8d02f45 100644 --- a/geminidr/core/tests/test_preprocess.py +++ b/geminidr/core/tests/test_preprocess.py @@ -317,14 +317,14 @@ def test_adu_to_electrons(astrofaker, caplog): p = NIRIImage([ad]) ad = p.ADUToElectrons()[0] - assert ad[0].nddata.unit == 'electron' + assert ad[0].unit == 'electron' assert_array_almost_equal(ad[0].data, ad[0].gain()) assert caplog.messages[2] == ('Converting N20010101S0001.fits from ADU to ' 'electrons by multiplying by the gain') caplog.clear() ad = p.ADUToElectrons()[0] - assert ad[0].nddata.unit == 'electron' + assert ad[0].unit == 'electron' assert_array_almost_equal(ad[0].data, ad[0].gain()) assert ('No changes will be made to N20010101S0001_ADUToElectrons.fits, ' 'since it has already been processed') in caplog.messages[2] diff --git a/geminidr/core/tests/test_spect.py b/geminidr/core/tests/test_spect.py index 39b7d39523..25e964b0bf 100644 --- a/geminidr/core/tests/test_spect.py +++ b/geminidr/core/tests/test_spect.py @@ -309,7 +309,7 @@ def test_flux_conservation_consistency(astrofaker, caplog, unit, ad = astrofaker.create("NIRI") ad.init_default_extensions() p = NIRIImage([ad]) # doesn't matter but we need a log object - ad.hdr["BUNIT"] = unit + ad[0].unit = unit conserve = primitives_spect.conserve_or_interpolate(ad[0], user_conserve=user_conserve, flux_calibrated=flux_calibrated, log=p.log) diff --git a/geminidr/gmos/tests/spect/test_calculate_sensitivity.py b/geminidr/gmos/tests/spect/test_calculate_sensitivity.py index 0f4e49f3c0..a662ad20d4 100644 --- a/geminidr/gmos/tests/spect/test_calculate_sensitivity.py +++ b/geminidr/gmos/tests/spect/test_calculate_sensitivity.py @@ -49,13 +49,13 @@ def _create_fake_data(): _ad = astrofaker.create('GMOS-S') _ad.add_extension(hdu, pixel_scale=1.0) + _ad[0].unit = "electron" _ad[0].data = _ad[0].data.ravel() + 1. _ad[0].mask = np.zeros(_ad[0].data.size, dtype=np.uint16) # ToDo Requires mask _ad[0].variance = np.ones_like(_ad[0].data) # ToDo Requires Variance _ad[0].phu.set('OBJECT', "DUMMY") _ad[0].phu.set('EXPTIME', 1.) - _ad[0].hdr.set('BUNIT', "electron") _ad[0].hdr.set('CTYPE1', "Wavelength") _ad[0].hdr.set('CUNIT1', "nm") _ad[0].hdr.set('CRPIX1', 1) @@ -280,4 +280,4 @@ def create_inputs_recipe(): if "--create-inputs" in sys.argv[1:]: create_inputs_recipe() else: - pytest.main() \ No newline at end of file + pytest.main() diff --git a/geminidr/gmos/tests/spect/test_flux_calibration.py b/geminidr/gmos/tests/spect/test_flux_calibration.py index f49612d92c..23b4e90468 100644 --- a/geminidr/gmos/tests/spect/test_flux_calibration.py +++ b/geminidr/gmos/tests/spect/test_flux_calibration.py @@ -67,7 +67,6 @@ def _get_spectrophotometric_data(object_name): return wavelength, flux def _create_fake_data(object_name): - from astropy.table import Table astrofaker = pytest.importorskip('astrofaker') wavelength, flux = _get_spectrophotometric_data(object_name) @@ -89,6 +88,7 @@ def _create_fake_data(object_name): _ad = astrofaker.create('GMOS-S') _ad.add_extension(hdu, pixel_scale=1.0) + _ad[0].unit = "electron" _ad[0].data = _ad[0].data.ravel() _ad[0].mask = np.zeros(_ad[0].data.size, dtype=np.uint16) # ToDo Requires mask _ad[0].variance = np.ones_like(_ad[0].data) # ToDo Requires Variance @@ -102,7 +102,6 @@ def _create_fake_data(object_name): _ad[0].hdr.set('NAXIS', 1) _ad[0].phu.set('OBJECT', object_name) _ad[0].phu.set('EXPTIME', 1.) - _ad[0].hdr.set('BUNIT', "electron") assert _ad.object() == object_name assert _ad.exposure_time() == 1 @@ -286,4 +285,4 @@ def create_inputs_recipe(): if "--create-inputs" in sys.argv[1:]: create_inputs_recipe() else: - pytest.main() \ No newline at end of file + pytest.main() diff --git a/gempy/library/spectral.py b/gempy/library/spectral.py index 34efc99334..10b5ddc102 100644 --- a/gempy/library/spectral.py +++ b/gempy/library/spectral.py @@ -41,19 +41,9 @@ def __init__(self, spectrum=None, spectral_axis=None, wcs=None, **kwargs): raise TypeError("Input spectrum must be a single AstroData slice") # Unit handling - try: # for NDData-like - if isinstance(spectrum, AstroData): - flux_unit = spectrum.nddata.unit - else: - flux_unit = spectrum.unit - except AttributeError: - try: # for AstroData - flux_unit = u.Unit(spectrum.hdr.get('BUNIT')) - except (TypeError, ValueError): # unknown/missing - flux_unit = None + if spectrum.unit is None: + spectrum.unit = u.dimensionless_unscaled - if flux_unit is None: - flux_unit = u.dimensionless_unscaled try: kwargs['mask'] = spectrum.mask except AttributeError: @@ -65,7 +55,7 @@ def __init__(self, spectrum=None, spectral_axis=None, wcs=None, **kwargs): # If spectrum was a Quantity, it already has units so we'd better # not multiply them in again! if not isinstance(flux, u.Quantity): - flux *= flux_unit + flux *= spectrum.unit # If no wavelength information is included, get it from the input if spectral_axis is None and wcs is None: diff --git a/gempy/library/tests/test_spectral.py b/gempy/library/tests/test_spectral.py index e30948d366..a07e8b7ab5 100644 --- a/gempy/library/tests/test_spectral.py +++ b/gempy/library/tests/test_spectral.py @@ -33,14 +33,13 @@ def _create_fake_data(): _ad.add_extension(hdu, pixel_scale=1.0) _ad[0].wcs = None # or else imaging WCS will be added - _ad[0].nddata.unit = units.electron + _ad[0].unit = units.electron _ad[0].data = _ad[0].data.ravel() + 1. _ad[0].mask = np.zeros(_ad[0].data.size, dtype=np.uint16) # ToDo Requires mask _ad[0].variance = np.ones_like(_ad[0].data) # ToDo Requires Variance _ad[0].phu.set('OBJECT', "DUMMY") _ad[0].phu.set('EXPTIME', 1.) - _ad[0].hdr.set('BUNIT', "electron") _ad[0].hdr.set('CTYPE1', "Wavelength") _ad[0].hdr.set('CUNIT1', "nm") _ad[0].hdr.set('CRPIX1', 1) From 20d6b43717b3206cd7d1eff806defcd2d314a517 Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Tue, 16 Mar 2021 19:45:26 -0300 Subject: [PATCH 15/20] Add comment when setting BUNIT --- astrodata/fits.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/astrodata/fits.py b/astrodata/fits.py index 4d8783974a..bea2ae4312 100644 --- a/astrodata/fits.py +++ b/astrodata/fits.py @@ -593,7 +593,8 @@ def ad_to_hdulist(ad): data_hdu = new_imagehdu(ext.data, header, 'SCI') if ext.unit is not None and ext.unit is not u.dimensionless_unscaled: - data_hdu.header['BUNIT'] = ext.unit.to_string() + data_hdu.header['BUNIT'] = (ext.unit.to_string(), + "Physical units of the array values") hdul.append(data_hdu) From b99352c65e5ce6e8d409c5d6ee7575f5847b09a4 Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Tue, 16 Mar 2021 19:51:38 -0300 Subject: [PATCH 16/20] Test that BUNIT is correctly written after ADUToElectrons --- geminidr/core/primitives_preprocess.py | 1 - geminidr/core/tests/test_preprocess.py | 13 ++++++++----- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/geminidr/core/primitives_preprocess.py b/geminidr/core/primitives_preprocess.py index 0bfd85822e..104eb8dd66 100644 --- a/geminidr/core/primitives_preprocess.py +++ b/geminidr/core/primitives_preprocess.py @@ -111,7 +111,6 @@ def ADUToElectrons(self, adinputs=None, suffix=None): # Update the headers of the AstroData Object. The pixel data now # has units of electrons so update the physical units keyword. - ad.hdr.set('BUNIT', 'electron', self.keyword_comments['BUNIT']) gt.mark_history(ad, primname=self.myself(), keyword=timestamp_key) ad.update_filename(suffix=suffix, strip=True) return adinputs diff --git a/geminidr/core/tests/test_preprocess.py b/geminidr/core/tests/test_preprocess.py index 0aa8d02f45..747ed98c20 100644 --- a/geminidr/core/tests/test_preprocess.py +++ b/geminidr/core/tests/test_preprocess.py @@ -1,16 +1,13 @@ -# import os -# from copy import deepcopy - import os import astrodata import gemini_instruments import numpy as np import pytest +from astropy.io import fits from astrodata.testing import download_from_archive from geminidr.core.primitives_preprocess import Preprocess from geminidr.gemini.lookups import DQ_definitions as DQ -# from geminidr.gmos.primitives_gmos_image import GMOSImage from geminidr.niri.primitives_niri_image import NIRIImage from gempy.library.astrotools import cartesian_regions_to_slices from numpy.testing import (assert_almost_equal, assert_array_almost_equal, @@ -309,7 +306,7 @@ def test_fixpixels_multiple_ext(niriprim2): assert_almost_equal(ad[1].data[sy, sx].max(), 60.333, decimal=2) -def test_adu_to_electrons(astrofaker, caplog): +def test_adu_to_electrons(astrofaker, caplog, tmp_path): caplog.set_level('INFO') ad = astrofaker.create("NIRI", "IMAGE") ad.init_default_extensions() @@ -329,6 +326,12 @@ def test_adu_to_electrons(astrofaker, caplog): assert ('No changes will be made to N20010101S0001_ADUToElectrons.fits, ' 'since it has already been processed') in caplog.messages[2] + testfile = tmp_path / 'test.fits' + ad.write(testfile) + assert fits.getval(testfile, 'BUNIT', extname='SCI') == 'electron' + ad = astrodata.read(testfile) + assert ad[0].unit == 'electron' + # TODO @bquint: clean up these tests From aadc1b3c070c48e104aace0893f34d6a5e49e475 Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Tue, 16 Mar 2021 20:15:45 -0300 Subject: [PATCH 17/20] Fix reading file --- geminidr/core/tests/test_preprocess.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geminidr/core/tests/test_preprocess.py b/geminidr/core/tests/test_preprocess.py index 747ed98c20..efc0ec4657 100644 --- a/geminidr/core/tests/test_preprocess.py +++ b/geminidr/core/tests/test_preprocess.py @@ -329,7 +329,7 @@ def test_adu_to_electrons(astrofaker, caplog, tmp_path): testfile = tmp_path / 'test.fits' ad.write(testfile) assert fits.getval(testfile, 'BUNIT', extname='SCI') == 'electron' - ad = astrodata.read(testfile) + ad = astrodata.open(testfile) assert ad[0].unit == 'electron' From 393e9386170be233c223ed3718196c62227744dc Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Wed, 17 Mar 2021 14:54:06 -0300 Subject: [PATCH 18/20] Fix fluxCalibrate to have correct unit for uncertainty --- astrodata/nddata.py | 49 +++++++++++++++++++++++++++++++ geminidr/core/primitives_spect.py | 23 ++++++++------- 2 files changed, 62 insertions(+), 10 deletions(-) diff --git a/astrodata/nddata.py b/astrodata/nddata.py index e76e892b93..69879a1eed 100644 --- a/astrodata/nddata.py +++ b/astrodata/nddata.py @@ -381,6 +381,55 @@ def unit(self, value): else: self._unit = u.Unit(value) + def convert_unit_to(self, unit, equivalencies=[]): + """ + Returns a new `NDData` object whose values have been converted + to a new unit. + + Copied from Astropy's NDDataArray. + + Parameters + ---------- + unit : `astropy.units.UnitBase` instance or str + The unit to convert to. + + equivalencies : list of equivalence pairs, optional + A list of equivalence pairs to try if the units are not + directly convertible. See :ref:`unit_equivalencies`. + + Returns + ------- + result : `~astropy.nddata.NDData` + The resulting dataset + + Raises + ------ + UnitsError + If units are inconsistent. + + """ + if self.unit is None: + raise ValueError("No unit specified on source data") + data = self.unit.to(unit, self.data, equivalencies=equivalencies) + if self.uncertainty is not None: + uncertainty_values = self.unit.to(unit, self.uncertainty.array, + equivalencies=equivalencies) + # should work for any uncertainty class + uncertainty = self.uncertainty.__class__(uncertainty_values) + else: + uncertainty = None + if self.mask is not None: + new_mask = self.mask.copy() + else: + new_mask = None + # Call __class__ in case we are dealing with an inherited type + result = self.__class__(data, uncertainty=uncertainty, + mask=new_mask, + wcs=self.wcs, + meta=self.meta, unit=unit) + + return result + def set_section(self, section, input): """ Sets only a section of the data. This method is meant to prevent diff --git a/geminidr/core/primitives_spect.py b/geminidr/core/primitives_spect.py index e010501be5..f50eb5ee06 100644 --- a/geminidr/core/primitives_spect.py +++ b/geminidr/core/primitives_spect.py @@ -2112,8 +2112,7 @@ def fluxCalibrate(self, adinputs=None, **params): if unit.is_equivalent(u.s): log.fullinfo("Dividing {} by exposure time of {} s". format(extname, exptime)) - ext /= exptime - ext.unit /= u.s + ext /= exptime * u.s elif not unit.is_equivalent(u.dimensionless_unscaled): log.warning(f"{extname} has incompatible units ('" f"{ext.unit}' and '{std_physical_unit}'" @@ -2159,15 +2158,19 @@ def fluxCalibrate(self, adinputs=None, **params): else: sens_factor *= 10**(0.4 * delta_airmass * extinction_correction) - final_sens_factor = (ext.unit * sens_factor / pixel_sizes).to( - final_units, equivalencies=u.spectral_density(waves)).value - + sens_factor /= pixel_sizes if ndim == 2 and dispaxis == 0: - ext *= final_sens_factor[:, np.newaxis] - else: - ext *= final_sens_factor - - ext.unit = final_units + sens_factor = sens_factor[:, np.newaxis] + + new_nddata = ext.nddata.multiply(sens_factor).convert_unit_to( + final_units, equivalencies=u.spectral_density(waves)) + # Replace manually data and uncertainty... Not really optimal + # but 1) we don't have a way to replace a NDData object in an + # AstroData one (as there may be other references), and 2) + ext.data = new_nddata.data + ext.unit = new_nddata.unit + ext.uncertainty.array = new_nddata.uncertainty.array + ext.uncertainty.unit = new_nddata.uncertainty.unit # Timestamp and update the filename gt.mark_history(ad, primname=self.myself(), keyword=timestamp_key) From c72a7c29303e9d591b6cebd8ae1295bd4f696b7b Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Wed, 17 Mar 2021 14:57:53 -0300 Subject: [PATCH 19/20] Avoid running test twice --- geminidr/gmos/tests/spect/test_flux_calibration.py | 1 - 1 file changed, 1 deletion(-) diff --git a/geminidr/gmos/tests/spect/test_flux_calibration.py b/geminidr/gmos/tests/spect/test_flux_calibration.py index 23b4e90468..e5853585f4 100644 --- a/geminidr/gmos/tests/spect/test_flux_calibration.py +++ b/geminidr/gmos/tests/spect/test_flux_calibration.py @@ -119,7 +119,6 @@ def _create_fake_data(object_name): @pytest.mark.gmosls @pytest.mark.preprocessed_data -@pytest.mark.regression @pytest.mark.parametrize("ad", test_datasets, indirect=True) def test_regression_on_flux_calibration(ad, ref_ad_factory, change_working_dir): """ From ae59e963acd6396b0a8c487cc8839f035f5ef921 Mon Sep 17 00:00:00 2001 From: Simon Conseil Date: Thu, 18 Mar 2021 16:08:39 -0300 Subject: [PATCH 20/20] Allow to set a new NDData --- astrodata/core.py | 9 +++++++++ geminidr/core/primitives_spect.py | 9 +-------- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/astrodata/core.py b/astrodata/core.py index f3be858ce7..8ea3208c8a 100644 --- a/astrodata/core.py +++ b/astrodata/core.py @@ -321,6 +321,15 @@ def nddata(self): """ return self._nddata[0] if self.is_single else self._nddata + @nddata.setter + @assign_only_single_slice + def nddata(self, new_nddata): + self.data = new_nddata.data + self.unit = new_nddata.unit + self.uncertainty = new_nddata.uncertainty + self.mask = new_nddata.mask + self.wcs = new_nddata.wcs + def table(self): # FIXME: do we need this in addition to .tables ? return self._tables.copy() diff --git a/geminidr/core/primitives_spect.py b/geminidr/core/primitives_spect.py index f50eb5ee06..f600913462 100644 --- a/geminidr/core/primitives_spect.py +++ b/geminidr/core/primitives_spect.py @@ -2162,15 +2162,8 @@ def fluxCalibrate(self, adinputs=None, **params): if ndim == 2 and dispaxis == 0: sens_factor = sens_factor[:, np.newaxis] - new_nddata = ext.nddata.multiply(sens_factor).convert_unit_to( + ext.nddata = ext.nddata.multiply(sens_factor).convert_unit_to( final_units, equivalencies=u.spectral_density(waves)) - # Replace manually data and uncertainty... Not really optimal - # but 1) we don't have a way to replace a NDData object in an - # AstroData one (as there may be other references), and 2) - ext.data = new_nddata.data - ext.unit = new_nddata.unit - ext.uncertainty.array = new_nddata.uncertainty.array - ext.uncertainty.unit = new_nddata.uncertainty.unit # Timestamp and update the filename gt.mark_history(ad, primname=self.myself(), keyword=timestamp_key)