Skip to content

Commit

Permalink
Array Quantities (#111)
Browse files Browse the repository at this point in the history 
* Rough test of array quanitties.

* Changed Quantity object instantiation in order to preserve datatype.

* Add support for a roman unit class

* Create asdf extension for `Unit`

* Fix remaining unit tests to support `Unit`

* Temporarily use rad branch with the new schema

* Update changes

* Test unit serialization as part of quantity

* Units updated.

* Add support for composite units

* Add force roman unit function

* Save for rebase.

* Rebased and fixed level 1 maker utility.

* Code clean up.

* Formatting fixes.

* Fix roman_datamodels tests

* Update to new asdf-astropy version

* Fix memmap test

* Rebase and updated Changelog.

* Fixed formatting.

* Pointed RAD to development version.

* Rough test of array quanitties.

* Changed Quantity object instantiation in order to preserve datatype.

* Units updated.

* Save for rebase.

* Rebased and fixed level 1 maker utility.

* Code clean up.

* Formatting fixes.

* Fix roman_datamodels tests

* Update to new asdf-astropy version

* Fix memmap test

* Rebase and updated Changelog.

* Fixed formatting.

* Pointed RAD to development version.

* Updated units.

* Temporary change to rad dependency to pass tests.

* Pointed RAD to main repo.

* Removed note on magnitude units.

Co-authored-by: William Jamieson <[email protected]>
  • Loading branch information
@PaulHuwe @WilliamJamieson
PaulHuwe and WilliamJamieson authored Dec 15, 2022
1 parent e6deb8e commit fa016af
Show file tree
Hide file tree
Showing 6 changed files with 200 additions and 89 deletions.
2 changes: 2 additions & 0 deletions CHANGES.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@
- Remove the unused ``stnode_test`` module. [#110]
- Add support for non-VOUnits to be used by Roman. [#109]

- Changed science arrays to quantities. [#111]


0.13.0 (2022-08-23)
===================
Expand Down
6 changes: 4 additions & 2 deletions setup.cfg
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,12 +22,14 @@ setup_requires =
setuptools_scm
install_requires =
asdf>=2.13.0
asdf-astropy>=0.2.0
asdf-astropy>=0.3.0
gwcs>=0.18.1
psutil>=5.7.2
numpy
astropy>=5.0.4
rad>=0.14.0
#rad>=0.14.0
# Temporary path for github PR tests
rad @git+https://github.com/spacetelescope/rad.git@main
asdf-standard>=1.0.3
package_dir =
=src
Expand Down
2 changes: 2 additions & 0 deletions src/roman_datamodels/stnode.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@
from .stuserdict import STUserDict as UserDict
from .units import Unit, CompositeUnit

from roman_datamodels.units import Unit, CompositeUnit, force_roman_unit

if sys.version_info < (3, 9):
import importlib_resources
else:
Expand Down
97 changes: 59 additions & 38 deletions src/roman_datamodels/testing/factories.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@
from astropy.modeling import models
import numpy as np

from roman_datamodels import units as ru

from .. import stnode
# from .. import table_definitions

Expand Down Expand Up @@ -136,35 +138,48 @@ def _random_bool():
return _random_choice(*[True, False])


def _random_array_float32(size=(4096, 4096), min=None, max=None):
def _random_array_float32(size=(4096, 4096), min=None, max=None, units=None):
if min is None:
min = np.finfo("float32").min
if max is None:
max = np.finfo("float32").max

array = np.random.default_rng().random(size=size, dtype=np.float32)
return min + max * array - min * array
array = min + max * array - min * array
if units:
array = u.Quantity(array, units, dtype=np.float32)
return array

def _random_array_uint8(size=(4096, 4096), min=None, max=None):
def _random_array_uint8(size=(4096, 4096), min=None, max=None, units=None):
if min is None:
min = np.iinfo("uint8").min
if max is None:
max = np.iinfo("uint8").max
return np.random.randint(min, high=max, size=size, dtype=np.uint8)
array = np.random.randint(min, high=max, size=size, dtype=np.uint8)
if units:
array = u.Quantity(array, units, dtype=np.uint8)
return array

def _random_array_uint16(size=(4096, 4096), min=None, max=None):
def _random_array_uint16(size=(4096, 4096), min=None, max=None, units=None):
if min is None:
min = np.iinfo("uint16").min
if max is None:
max = np.iinfo("uint16").max
return np.random.randint(min, high=max, size=size, dtype=np.uint16)
array = np.random.randint(min, high=max, size=size, dtype=np.uint16)
if units:
array = u.Quantity(array, units, dtype=np.uint16)
return array


def _random_array_uint32(size=(4096, 4096), min=None, max=None):
def _random_array_uint32(size=(4096, 4096), min=None, max=None, units=None):
if min is None:
min = np.iinfo("uint32").min
if max is None:
max = np.iinfo("uint32").max
return np.random.randint(min, high=max, size=size, dtype=np.uint32)
array = np.random.randint(min, high=max, size=size, dtype=np.uint32)
if units:
array = u.Quantity(array, units, dtype=np.uint32)
return array


def _random_exposure_type():
Expand Down Expand Up @@ -1030,15 +1045,17 @@ def create_ramp(**kwargs):

raw = {
"meta": create_meta(),
"data": _random_array_float32((2, 4096, 4096)),

"data": _random_array_float32((2, 4096, 4096), units=ru.electron),
"pixeldq": _random_array_uint32((4096, 4096)),
"groupdq": _random_array_uint8((2, 4096, 4096)),
"err": _random_array_float32(size=(2, 4096, 4096),min=0.0),
"amp33": _random_array_uint16((2, 4096, 128)),
"border_ref_pix_right": _random_array_float32((2, 4096, 4)),
"border_ref_pix_left": _random_array_float32((2, 4096, 4)),
"border_ref_pix_top": _random_array_float32((2, 4, 4096)),
"border_ref_pix_bottom": _random_array_float32((2, 4, 4096)),
"err": _random_array_float32(size=(2, 4096, 4096),min=0.0, units=ru.electron),

"amp33": _random_array_uint16((2, 4096, 128), units=ru.DN),
"border_ref_pix_right": _random_array_float32((2, 4096, 4), units=ru.DN),
"border_ref_pix_left": _random_array_float32((2, 4096, 4), units=ru.DN),
"border_ref_pix_top": _random_array_float32((2, 4, 4096), units=ru.DN),
"border_ref_pix_bottom": _random_array_float32((2, 4, 4096), units=ru.DN),
"dq_border_ref_pix_right": _random_array_uint32((4096, 4)),
"dq_border_ref_pix_left": _random_array_uint32((4096, 4)),
"dq_border_ref_pix_top": _random_array_uint32((4, 4096)),
Expand Down Expand Up @@ -1067,15 +1084,17 @@ def create_ramp_fit_output(**kwargs):

raw = {
"meta": create_meta(),
"slope": _random_array_float32(seg_shape),
"sigslope": _random_array_float32(seg_shape),
"yint": _random_array_float32(seg_shape),
"sigyint": _random_array_float32(seg_shape),
"pedestal": _random_array_float32(seg_shape[1:]),

"slope": _random_array_float32(seg_shape, units=ru.electron / u.s),
"sigslope": _random_array_float32(seg_shape, units=ru.electron / u.s),
"yint": _random_array_float32(seg_shape, units=ru.electron),
"sigyint": _random_array_float32(seg_shape, units=ru.electron),
"pedestal": _random_array_float32(seg_shape[1:], units=ru.electron),
"weights": _random_array_float32(seg_shape),
"crmag": _random_array_float32(seg_shape),
"var_poisson": _random_array_float32(seg_shape),
"var_rnoise": _random_array_float32(seg_shape),
"crmag": _random_array_float32(seg_shape, units=ru.electron),
"var_poisson": _random_array_float32(seg_shape, units=ru.electron**2 / u.s**2),
"var_rnoise": _random_array_float32(seg_shape, units=ru.electron**2 / u.s**2),

}
raw.update(kwargs)

Expand All @@ -1101,9 +1120,9 @@ def create_guidewindow(**kwargs):

raw = {
"meta": create_meta(),
"pedestal_frames": _random_array_uint16(seg_shape),
"signal_frames": _random_array_uint16(seg_shape),
"amp33": _random_array_uint16(seg_shape),
"pedestal_frames": _random_array_uint16(seg_shape, units=ru.DN),
"signal_frames": _random_array_uint16(seg_shape, units=ru.DN),
"amp33": _random_array_uint16(seg_shape, units=ru.DN),
}
raw.update(kwargs)

Expand Down Expand Up @@ -1284,18 +1303,20 @@ def create_wfi_image(**kwargs):
"""

raw = {
"data": _random_array_float32((4088, 4088)),
"data": _random_array_float32((4088, 4088), units=ru.electron / u.s),
"dq": _random_array_uint32((4088, 4088)),
"err": _random_array_float32((4088, 4088), min=0.0),
"err": _random_array_float32((4088, 4088), min=0.0, units=ru.electron / u.s),
"meta": create_meta(),
"var_flat": _random_array_float32((4088, 4088)),
"var_poisson": _random_array_float32((4088, 4088)),
"var_rnoise": _random_array_float32((4088, 4088)),
"amp33": _random_array_uint16((2, 4096, 128)),
"border_ref_pix_right": _random_array_float32((2, 4096, 4)),
"border_ref_pix_left": _random_array_float32((2, 4096, 4)),
"border_ref_pix_top": _random_array_float32((2, 4, 4096)),
"border_ref_pix_bottom": _random_array_float32((2, 4, 4096)),

"var_flat": _random_array_float32((4088, 4088), units=ru.electron**2 / u.s**2),
"var_poisson": _random_array_float32((4088, 4088), units=ru.electron**2 / u.s**2),
"var_rnoise": _random_array_float32((4088, 4088), units=ru.electron**2 / u.s**2),

"amp33": _random_array_uint16((2, 4096, 128), units=ru.DN),
"border_ref_pix_right": _random_array_float32((2, 4096, 4), units=ru.DN),
"border_ref_pix_left": _random_array_float32((2, 4096, 4), units=ru.DN),
"border_ref_pix_top": _random_array_float32((2, 4, 4096), units=ru.DN),
"border_ref_pix_bottom": _random_array_float32((2, 4, 4096), units=ru.DN),
"dq_border_ref_pix_right": _random_array_uint32((4096, 4)),
"dq_border_ref_pix_left": _random_array_uint32((4096, 4)),
"dq_border_ref_pix_top": _random_array_uint32((4, 4096)),
Expand Down Expand Up @@ -1348,8 +1369,8 @@ def create_wfi_science_raw(**kwargs):
"""
raw = {
# TODO: What should this shape be?
"data": _random_array_uint16((2, 4096, 4096)),
"amp33": _random_array_uint16((2, 4096, 128)),
"data": _random_array_uint16((2, 4096, 4096), units=ru.DN),
"amp33": _random_array_uint16((2, 4096, 128), units=ru.DN),
"meta": create_meta(),
}
raw.update(kwargs)
Expand Down
99 changes: 60 additions & 39 deletions src/roman_datamodels/testing/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
import numpy as np
from astropy import units as u
from astropy.modeling import models
from roman_datamodels import units as ru

from .factories import _random_positive_int, _random_string
from .. import stnode
Expand Down Expand Up @@ -450,10 +451,12 @@ def mk_level1_science_raw(shape=None, filepath=None):
else:
n_ints = shape[0]

wfi_science_raw['data'] = np.zeros(shape, dtype=np.uint16)
wfi_science_raw['data'] = u.Quantity(np.zeros(shape, dtype=np.uint16),
ru.DN, dtype=np.uint16)

# add amp 33 ref pix
wfi_science_raw['amp33'] = np.zeros((n_ints, 4096, 128), dtype=np.uint16)
wfi_science_raw['amp33'] = u.Quantity(np.zeros((n_ints, 4096, 128), dtype=np.uint16),
ru.DN, dtype=np.uint16)

if filepath:
af = asdf.AsdfFile()
Expand Down Expand Up @@ -501,14 +504,15 @@ def mk_level2_image(shape=None, n_ints=None, filepath=None):
n_ints = 8

# add border reference pixel arrays
wfi_image['border_ref_pix_left'] = np.zeros((n_ints, shape[0] + 8, 4),
dtype=np.float32)
wfi_image['border_ref_pix_right'] = np.zeros((n_ints, shape[0] + 8, 4),
dtype=np.float32)
wfi_image['border_ref_pix_top'] = np.zeros((n_ints, 4, shape[1] + 8),
dtype=np.float32)
wfi_image['border_ref_pix_bottom'] = np.zeros((n_ints, 4, shape[1] + 8),
dtype=np.float32)
wfi_image['border_ref_pix_left'] = u.Quantity(np.zeros((n_ints, shape[0] + 8, 4), dtype=np.float32),
ru.DN, dtype=np.float32)
wfi_image['border_ref_pix_right'] = u.Quantity(np.zeros((n_ints, shape[0] + 8, 4), dtype=np.float32),
ru.DN, dtype=np.float32)
wfi_image['border_ref_pix_top'] = u.Quantity(np.zeros((n_ints, shape[0] + 8, 4), dtype=np.float32),
ru.DN, dtype=np.float32)
wfi_image['border_ref_pix_bottom'] = u.Quantity(np.zeros((n_ints, shape[0] + 8, 4), dtype=np.float32),
ru.DN, dtype=np.float32)


# and their dq arrays
wfi_image['dq_border_ref_pix_left'] = np.zeros((shape[0] + 8, 4), dtype=np.uint32)
Expand All @@ -518,14 +522,21 @@ def mk_level2_image(shape=None, n_ints=None, filepath=None):

# add amp 33 ref pixel array
amp33_size = (n_ints, 4096, 128)
wfi_image['amp33'] = np.zeros(amp33_size, dtype=np.uint16)
wfi_image['amp33'] = u.Quantity(np.zeros(amp33_size, dtype=np.uint16),
ru.DN, dtype=np.uint16)

wfi_image['data'] = np.zeros(shape, dtype=np.float32)
wfi_image['data'] = u.Quantity(np.zeros(shape, dtype=np.float32),
ru.electron / u.s, dtype=np.float32)
wfi_image['dq'] = np.zeros(shape, dtype=np.uint32)
wfi_image['err'] = np.zeros(shape, dtype=np.float32)
wfi_image['var_poisson'] = np.zeros(shape, dtype=np.float32)
wfi_image['var_rnoise'] = np.zeros(shape, dtype=np.float32)
wfi_image['var_flat'] = np.zeros(shape, dtype=np.float32)
wfi_image['err'] = u.Quantity(np.zeros(shape, dtype=np.float32),
ru.electron / u.s, dtype=np.float32)

wfi_image['var_poisson'] = u.Quantity(np.zeros(shape, dtype=np.float32),
ru.electron**2 / u.s**2, dtype=np.float32)
wfi_image['var_rnoise'] = u.Quantity(np.zeros(shape, dtype=np.float32),
ru.electron**2 / u.s**2, dtype=np.float32)
wfi_image['var_flat'] = u.Quantity(np.zeros(shape, dtype=np.float32),
ru.electron**2 / u.s**2, dtype=np.float32)
wfi_image['cal_logs'] = mk_cal_logs()

if filepath:
Expand Down Expand Up @@ -946,14 +957,14 @@ def mk_ramp(shape=None, n_ints=None, filepath=None):
shape = (8, 4096, 4096)

# add border reference pixel arrays
ramp['border_ref_pix_left'] = np.zeros((shape[0], shape[1], 4),
dtype=np.float32)
ramp['border_ref_pix_right'] = np.zeros((shape[0], shape[1], 4),
dtype=np.float32)
ramp['border_ref_pix_top'] = np.zeros((shape[0], 4, shape[2]),
dtype=np.float32)
ramp['border_ref_pix_bottom'] = np.zeros((shape[0], 4, shape[2]),
dtype=np.float32)
ramp['border_ref_pix_left'] = u.Quantity(np.zeros((shape[0], shape[1], 4), dtype=np.float32),
ru.DN, dtype=np.float32)
ramp['border_ref_pix_right'] = u.Quantity(np.zeros((shape[0], shape[1], 4), dtype=np.float32),
ru.DN, dtype=np.float32)
ramp['border_ref_pix_top'] = u.Quantity(np.zeros((shape[0], 4, shape[2]), dtype=np.float32),
ru.DN, dtype=np.float32)
ramp['border_ref_pix_bottom'] = u.Quantity(np.zeros((shape[0], 4, shape[2]), dtype=np.float32),
ru.DN, dtype=np.float32)

# and their dq arrays
ramp['dq_border_ref_pix_left'] = np.zeros((shape[1], 4), dtype=np.uint32)
Expand All @@ -962,13 +973,12 @@ def mk_ramp(shape=None, n_ints=None, filepath=None):
ramp['dq_border_ref_pix_bottom'] = np.zeros((4, shape[2]), dtype=np.uint32)

# add amp 33 ref pixel array
amp33_size = (shape[0], 4096, 128)
ramp['amp33'] = np.zeros(amp33_size, dtype=np.uint16)
ramp['amp33'] = u.Quantity(np.full(shape, 1.0, dtype=np.uint16), ru.DN, dtype=np.uint16)

ramp['data'] = np.full(shape, 1.0, dtype=np.float32)
ramp['data'] = u.Quantity(np.full(shape, 1.0, dtype=np.float32), ru.DN, dtype=np.float32)
ramp['pixeldq'] = np.zeros(shape[1:], dtype=np.uint32)
ramp['groupdq'] = np.zeros(shape, dtype=np.uint8)
ramp['err'] = np.zeros(shape, dtype=np.float32)
ramp['err'] = u.Quantity(np.zeros(shape, dtype=np.float32), ru.DN, dtype=np.float32)

if filepath:
af = asdf.AsdfFile()
Expand Down Expand Up @@ -1001,15 +1011,23 @@ def mk_rampfitoutput(shape=None, filepath=None):
if not shape:
shape = (8, 4096, 4096)

rampfitoutput['slope'] = np.zeros(shape, dtype=np.float32)
rampfitoutput['sigslope'] = np.zeros(shape, dtype=np.float32)
rampfitoutput['yint'] = np.zeros(shape, dtype=np.float32)
rampfitoutput['sigyint'] = np.zeros(shape, dtype=np.float32)
rampfitoutput['pedestal'] = np.zeros(shape[1:], dtype=np.float32)
rampfitoutput['slope'] = u.Quantity(np.zeros(shape, dtype=np.float32),
ru.electron / u.s, dtype=np.float32)
rampfitoutput['sigslope'] = u.Quantity(np.zeros(shape, dtype=np.float32),
ru.electron / u.s, dtype=np.float32)
rampfitoutput['yint'] = u.Quantity(np.zeros(shape, dtype=np.float32),
ru.electron, dtype=np.float32)
rampfitoutput['sigyint'] = u.Quantity(np.zeros(shape, dtype=np.float32),
ru.electron, dtype=np.float32)
rampfitoutput['pedestal'] = u.Quantity(np.zeros(shape[1:], dtype=np.float32),
ru.electron, dtype=np.float32)
rampfitoutput['weights'] = np.zeros(shape, dtype=np.float32)
rampfitoutput['crmag'] = np.zeros(shape, dtype=np.float32)
rampfitoutput['var_poisson'] = np.zeros(shape, dtype=np.float32)
rampfitoutput['var_rnoise'] = np.zeros(shape, dtype=np.float32)
rampfitoutput['crmag'] = u.Quantity(np.zeros(shape, dtype=np.float32),
ru.electron, dtype=np.float32)
rampfitoutput['var_poisson'] = u.Quantity(np.zeros(shape, dtype=np.float32),
ru.electron**2 / u.s**2, dtype=np.float32)
rampfitoutput['var_rnoise'] = u.Quantity(np.zeros(shape, dtype=np.float32),
ru.electron**2 / u.s**2, dtype=np.float32)

if filepath:
af = asdf.AsdfFile()
Expand Down Expand Up @@ -1072,9 +1090,12 @@ def mk_guidewindow(shape=None, filepath=None):
if not shape:
shape = (2, 8, 16, 32, 32)

guidewindow['pedestal_frames'] = np.zeros(shape, dtype=np.uint16)
guidewindow['signal_frames'] = np.zeros(shape, dtype=np.uint16)
guidewindow['amp33'] = np.zeros(shape, dtype=np.uint16)
guidewindow['pedestal_frames'] = u.Quantity(np.zeros(shape, dtype=np.uint16),
ru.DN, dtype=np.uint16)
guidewindow['signal_frames'] = u.Quantity(np.zeros(shape, dtype=np.uint16),
ru.DN, dtype=np.uint16)
guidewindow['amp33'] = u.Quantity(np.zeros(shape, dtype=np.uint16),
ru.DN, dtype=np.uint16)

if filepath:
af = asdf.AsdfFile()
Expand Down
Loading

0 comments on commit fa016af

Please sign in to comment.