Remove units from roman data model outputs (#165)

* Remove units from romanisim.

docs/romanisim/l2.rst

@@ -1,7 +1,7 @@
 L2 images
 =========
 
-L2 images are constructed from :doc:`L1 </romanisim/l1>` images, which are in turn constructed from idealized :doc:`count </romanisim/image>` images.  This means that even when constructing L2 images, one must go through the process of simulating how electrons counts get apportioned among the various reads.  It is challenging to realistically model the statistics in the noise of L2 images without going through this process.
+L2 images are constructed from :doc:`L1 </romanisim/l1>` images, which are in turn constructed from idealized :doc:`count </romanisim/image>` images.  This means that even when constructing L2 images, one must go through the process of simulating how electrons get apportioned among the various reads.  It is challenging to realistically model the statistics in the noise of L2 images without going through this process.
 
 L2 images are constructed by doing "ramp fitting" on the level 1 images.  Each pixel of a level 1 image is a series of "resultants", giving the measured value of that pixel averaged over a series on non-destructive reads as the exposure is being observed.  A simple model for a pixel is that its flux as a function of time is simply two numbers: a pedestal and a linear ramp representing the rate at which photons are detected by the pixel.  Ramp fitting turns these one-dimensional series of resultants into a "slope" image that is of interest astronomically.  Due to details of the H4RG detectors, the pedestals of the ramp vary widely from exposure to exposure, and so current fitting completely throws away as non-astronomical any information in the ramp that is sensitive to the pedestal.
 
@@ -20,10 +20,11 @@ This is a fairly faithful representation of how level two image construction wor
 L2 source injection
 ===================
 We support L2 source injection through :meth:`romanisim.image.inject_sources_into_l2`.  This routine is intended for cases when existing L2 images are available, and a small number of sources need to be added on top of those images.  The L2 source injection is intended to be a fairly faithful simulation of the process.  It includes the following steps:
+
 * The model WCS is used to generate the pixel locations for sources.
 * The number of additional electrons in each pixel is simulated.
 * A virtual ramp is created that evenly apportions the original L2 flux along the ramp.
-* The new counts are apportioned randomly along the ramp following a Poisson distribution in the same way as for the usual L1 simulations.
+* The new electrons are apportioned randomly along the ramp following a Poisson distribution in the same way as for the usual L1 simulations.
 * The new ramp is then fit using the ramp fitting code, and the new fit replaces the old fit in the L2 file.
 
 Note that this procedure does not include any new cosmic rays, or any special treatment of pixels that had cosmic rays in them and now have additional source flux.  It is hard to do this right since we do not know which pixel the cosmic ray landed in on the basis of the L2 file.  If you need more accuracy than this you probably need to directly inject sources into the L1 file.

pyproject.toml

@@ -17,16 +17,16 @@ classifiers = [
     "Programming Language :: Python :: 3.10",
 ]
 dependencies = [
-    "asdf >=2.15.1",
+    "asdf >=3.3.0",
     "astropy >=5.3",
     "crds >=11.16.16",
     "defusedxml >=0.5.0",
     "galsim >=2.5.1",
-    "rad >=0.19.2",
-    "roman_datamodels >=0.19.1",
-    # "rad @ git+https://github.com/spacetelescope/rad.git",
-    # "roman_datamodels @ git+https://github.com/spacetelescope/roman_datamodels.git",
-    "gwcs >=0.18.1",
+    # "rad >=0.19.2",
+    # "roman_datamodels >=0.19.1",
+    "rad @ git+https://github.com/spacetelescope/rad.git",
+    "roman_datamodels @ git+https://github.com/spacetelescope/roman_datamodels.git",
+    "gwcs >=0.19.0",
     "jsonschema >=4.8",
     "numpy >=1.22",
     "webbpsf >=1.2.1",

romanisim/image.py

@@ -860,8 +860,6 @@ def make_asdf(slope, slopevar_rn, slopevar_poisson, metadata=None,
               filepath=None, persistence=None, dq=None, imwcs=None,
               gain=None):
     """Wrap a galsim simulated image with ASDF/roman_datamodel metadata.
-
-    Eventually this needs to get enough info to reconstruct a refit WCS.
     """
 
     out = maker_utils.mk_level2_image(
@@ -898,13 +896,13 @@ def make_asdf(slope, slopevar_rn, slopevar_poisson, metadata=None,
 
     util.update_photom_keywords(out, gain=gain)
 
-    out['data'] = slope
+    out['data'] = slope.value
     out['dq'] = np.zeros(slope.shape, dtype='u4')
     if dq is not None:
         out['dq'][:, :] = dq
-    out['var_poisson'] = slopevar_poisson
-    out['var_rnoise'] = slopevar_rn
-    out['var_flat'] = slopevar_rn * 0
+    out['var_poisson'] = slopevar_poisson.value
+    out['var_rnoise'] = slopevar_rn.value
+    out['var_flat'] = slopevar_rn.value * 0
     out['err'] = np.sqrt(out['var_poisson'] + out['var_rnoise'] + out['var_flat'])
     extras = dict()
     if persistence is not None:
@@ -1030,10 +1028,10 @@ def inject_sources_into_l2(model, cat, x=None, y=None, psf=None, rng=None,
     # create injected source ramp resultants
     resultants, dq = romanisim.l1.apportion_counts_to_resultants(
         sourcecounts.array[m], tij, rng=rng)
-    resultants *= u.electron
+    resultants = resultants * u.electron
 
     # Inject source to original image
-    newramp = model.data[None, :] * tbar[:, None, None] * u.s
+    newramp = model.data[None, :] * tbar[:, None, None] * u.DN
     newramp[:, m] += resultants / gain
     # newramp has units of DN
 

romanisim/l1.py

@@ -416,7 +416,7 @@ def make_asdf(resultants, dq=None, filepath=None, metadata=None, persistence=Non
     if metadata is not None:
         out['meta'].update(metadata)
     extras = dict()
-    out['data'][:, nborder:-nborder, nborder:-nborder] = resultants
+    out['data'][:, nborder:-nborder, nborder:-nborder] = resultants.value
     if dq is not None:
         extras['dq'] = np.zeros(out['data'].shape, dtype='i4')
         extras['dq'][:, nborder:-nborder, nborder:-nborder] = dq

romanisim/l3.py

@@ -72,7 +72,7 @@ def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos, ypos, psf,
     # Create Image canvas to add objects to
     if isinstance(l3_mos, (rdm.MosaicModel, WfiMosaic)):
         sourcecountsall = galsim.ImageF(
-            l3_mos.data.value, wcs=romanisim.wcs.GWCS(l3_mos.meta.wcs),
+            l3_mos.data, wcs=romanisim.wcs.GWCS(l3_mos.meta.wcs),
             xmin=0, ymin=0)
     else:
         sourcecountsall = l3_mos
@@ -87,7 +87,7 @@ def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos, ypos, psf,
 
     # Save array with added sources
     if isinstance(l3_mos, (rdm.MosaicModel, WfiMosaic)):
-        l3_mos.data = sourcecountsall.array * l3_mos.data.unit
+        l3_mos.data = sourcecountsall.array
 
     return outinfo
 
@@ -162,26 +162,26 @@ def inject_sources_into_l3(model, cat, x=None, y=None, psf=None, rng=None,
         # Set scaling factor for injected sources
         # Flux / sigma_p^2
         xidx, yidx = int(np.round(x0)), int(np.round(y0))
-        if model.var_poisson[yidx, xidx].value != 0:
+        if model.var_poisson[yidx, xidx] != 0:
             Ct.append(math.fabs(
-                model.data[yidx, xidx].value /
-                model.var_poisson[yidx, xidx].value))
+                model.data[yidx, xidx] /
+                model.var_poisson[yidx, xidx]))
         else:
             Ct.append(1.0)
     Ct = np.array(Ct)
     # etomjysr = 1/C; C converts fluxes to electrons
     exptimes = Ct * etomjysr
 
-    Ct_all = (model.data.value /
-              (model.var_poisson.value + (model.var_poisson.value == 0)))
+    Ct_all = (model.data /
+              (model.var_poisson + (model.var_poisson == 0)))
 
     # compute the total number of counts we got from the source
     res = add_objects_to_l3(
         model, cat, exptimes, x, y, psf, etomjysr=etomjysr,
         maggytoes=maggytoes, filter_name=filter_name, bandpass=None,
         rng=rng)
 
-    model.var_poisson = (model.data.value / Ct_all) * model.var_poisson.unit
+    model.var_poisson = (model.data / Ct_all)
 
     return res
 
@@ -655,7 +655,7 @@ def make_l3(image, metadata, efftimes, var_poisson=None,
     """
 
     # Create mosaic data object
-    mosaic = image.array.copy() * u.MJy / u.sr
+    mosaic = image.array.copy()
 
     # Ensure that effective times are an array
     if isinstance(efftimes, np.ndarray):
@@ -677,11 +677,11 @@ def make_l3(image, metadata, efftimes, var_poisson=None,
     context = (np.ones((1,) + mosaic.shape, dtype=np.uint32)
                if context is None else context)
 
-    mosaic_node.var_poisson[...] = var_poisson * mosaic.unit ** 2
-    mosaic_node.var_rnoise[...] = var_rnoise * mosaic.unit ** 2
-    mosaic_node.var_flat[...] = var_flat * mosaic.unit ** 2
+    mosaic_node.var_poisson[...] = var_poisson
+    mosaic_node.var_rnoise[...] = var_rnoise
+    mosaic_node.var_flat[...] = var_flat
     mosaic_node.err[...] = np.sqrt(
-        var_poisson + var_rnoise + var_flat) * mosaic.unit
+        var_poisson + var_rnoise + var_flat)
     mosaic_node.context = context
 
 

romanisim/tests/test_image.py

@@ -658,15 +658,15 @@ def test_inject_source_into_image():
     iminj = image.inject_sources_into_l2(im, catinj, x=[xpos], y=[ypos])
 
     # Test that all pixels near the PSF are different from the original values
-    assert np.all((im.data.value[ypos - 1:ypos + 2, xpos - 1:xpos + 2] !=
-                   iminj.data.value[ypos - 1:ypos + 2, xpos - 1: xpos + 2]))
+    assert np.all((im.data[ypos - 1:ypos + 2, xpos - 1:xpos + 2] !=
+                   iminj.data[ypos - 1:ypos + 2, xpos - 1: xpos + 2]))
 
     # Test that pixels far from the injected source are close to the original image
-    assert np.all(im.data.value[-10:, -10:] == iminj.data.value[-10:, -10:])
+    assert np.all(im.data[-10:, -10:] == iminj.data[-10:, -10:])
 
     # Test that the amount of added flux makes sense
-    fluxeps = flux * romanisim.bandpass.get_abflux('F158') * u.electron / u.s
-    assert np.abs(np.sum(iminj.data - im.data) * parameters.reference_data['gain'] /
+    fluxeps = flux * romanisim.bandpass.get_abflux('F158')  # u.electron / u.s
+    assert np.abs(np.sum(iminj.data - im.data) * parameters.reference_data['gain'].value /
                   fluxeps - 1) < 0.1
 
     # Create log entry and artifacts
@@ -733,7 +733,7 @@ def test_image_input(tmpdir):
     res = image.simulate(meta, tab, usecrds=False, webbpsf=False)
 
     # did we get all the flux?
-    totflux = np.sum(res[0].data.value - np.median(res[0].data.value))
+    totflux = np.sum(res[0].data - np.median(res[0].data))
     expectedflux = (romanisim.bandpass.get_abflux('F087') * np.sum(tab['F087'])
                     / parameters.reference_data['gain'].value)
     assert np.abs(totflux / expectedflux - 1) < 0.1
@@ -743,7 +743,7 @@ def test_image_input(tmpdir):
         x, y = imwcs.toImage(r, d, units=galsim.degrees)
         x = int(x)
         y = int(y)
-        assert res[0].data.value[y, x] > np.median(res[0].data.value) * 5
+        assert res[0].data[y, x] > np.median(res[0].data) * 5
     log.info('DMS228: Successfully added rendered sources from image input; '
              'sources are present and flux matches.')
 

romanisim/tests/test_l3.py

@@ -67,16 +67,16 @@ def test_inject_sources_into_mosaic():
     unit_factor = romanisim.bandpass.etomjysr(filter_name)
 
     # Populate the mosaic data array with gaussian noise from generators
-    g1.generate(l3_mos.data.value[0:100, 0:100])
-    g2.generate(l3_mos.data.value[0:100, 100:200])
-    g3.generate(l3_mos.data.value[100:200, 0:100])
-    g4.generate(l3_mos.data.value[100:200, 100:200])
+    g1.generate(l3_mos.data[0:100, 0:100])
+    g2.generate(l3_mos.data[0:100, 100:200])
+    g3.generate(l3_mos.data[100:200, 0:100])
+    g4.generate(l3_mos.data[100:200, 100:200])
 
     # Define Poisson Noise of mosaic
-    l3_mos.var_poisson.value[0:100, 0:100] = 0.01**2 * meanflux**2
-    l3_mos.var_poisson.value[0:100, 100:200] = 0.02**2 * meanflux**2
-    l3_mos.var_poisson.value[100:200, 0:100] = 0.05**2 * meanflux**2
-    l3_mos.var_poisson.value[100:200, 100:200] = 0.1**2 * meanflux**2
+    l3_mos.var_poisson[0:100, 0:100] = 0.01**2 * meanflux**2
+    l3_mos.var_poisson[0:100, 100:200] = 0.02**2 * meanflux**2
+    l3_mos.var_poisson[100:200, 0:100] = 0.05**2 * meanflux**2
+    l3_mos.var_poisson[100:200, 100:200] = 0.1**2 * meanflux**2
 
     # Create normalized psf source catalog (same source in each quadrant)
     mag_flux = 1e-9
@@ -103,18 +103,18 @@ def test_inject_sources_into_mosaic():
 
     # Ensure that every data pixel value has increased or
     # remained the same with the new sources injected
-    assert np.all(l3_mos.data.value >= l3_mos_orig.data.value)
+    assert np.all(l3_mos.data >= l3_mos_orig.data)
 
     # Ensure that every pixel's poisson variance has increased or
     # remained the same with the new sources injected
     # Numpy isclose is needed to determine equality, due to float precision issues
-    close_mask = np.isclose(l3_mos.var_poisson.value, l3_mos_orig.var_poisson.value, rtol=1e-06)
+    close_mask = np.isclose(l3_mos.var_poisson, l3_mos_orig.var_poisson, rtol=1e-06)
     assert False in close_mask
-    assert np.all(l3_mos.var_poisson.value[~close_mask] > l3_mos_orig.var_poisson.value[~close_mask])
+    assert np.all(l3_mos.var_poisson[~close_mask] > l3_mos_orig.var_poisson[~close_mask])
 
     # Ensure total added flux matches expected added flux
     total_rec_flux = np.sum(l3_mos.data - l3_mos_orig.data)  # MJy / sr
-    total_theo_flux = 4 * mag_flux * cps_conv * unit_factor * u.MJy / u.sr
+    total_theo_flux = 4 * mag_flux * cps_conv * unit_factor  # u.MJy / u.sr
     assert np.isclose(total_rec_flux, total_theo_flux, rtol=4e-02)
 
     # Create log entry and artifacts
@@ -190,11 +190,11 @@ def test_sim_mosaic():
     for x, y in zip(x_all, y_all):
         x = int(x)
         y = int(y)
-        assert mosaic.data.value[y, x] > (np.median(mosaic.data.value) * 5)
+        assert mosaic.data[y, x] > (np.median(mosaic.data) * 5)
 
     # Did we get all the flux?
     etomjysr = romanisim.bandpass.etomjysr(filter_name)
-    totflux = np.sum(mosaic.data.value - np.median(mosaic.data.value)) / etomjysr
+    totflux = np.sum(mosaic.data - np.median(mosaic.data)) / etomjysr
 
     # Flux to counts
     cps_conv = romanisim.bandpass.get_abflux(filter_name)
@@ -205,7 +205,7 @@ def test_sim_mosaic():
 
     # Is the noise about right?
     assert np.abs(
-        mad_std(mosaic.data.value) / np.median(mosaic.err.value) - 1) < 0.5
+        mad_std(mosaic.data) / np.median(mosaic.err) - 1) < 0.5
     # note large ~50% error bar there; value in initial test run is 0.25
     # a substantial number of source pixels have flux, so the simple medians
     # and mads aren't terribly right.
@@ -351,8 +351,8 @@ def test_simulate_vs_cps():
                                )
 
     # Ensure that the simulate and simulate_cps output matches for each type
-    assert np.allclose(im1.array, im3['data'].value)
-    assert np.allclose(im2.array, im4['data'].value)
+    assert np.allclose(im1.array, im3['data'])
+    assert np.allclose(im2.array, im4['data'])
 
 
 def test_simulate_cps():
@@ -500,12 +500,12 @@ def test_exptime_array():
                                )
 
     # Ensure that the poisson variance scales with exposure time difference
-    assert np.isclose((np.median(im1['var_poisson'][0:50, :].value) / np.median(im1['var_poisson'][50:, :].value)), expfactor, rtol=0.02)
-    assert np.isclose((np.median(im2['var_poisson'][0:50, :].value) / np.median(im2['var_poisson'][50:, :].value)), expfactor, rtol=0.02)
+    assert np.isclose((np.median(im1['var_poisson'][0:50, :]) / np.median(im1['var_poisson'][50:, :])), expfactor, rtol=0.02)
+    assert np.isclose((np.median(im2['var_poisson'][0:50, :]) / np.median(im2['var_poisson'][50:, :])), expfactor, rtol=0.02)
 
     # Ensure that the data remains consistent across the exposure times
-    assert np.isclose(np.median(im1['data'][0:50, :].value), np.median(im1['data'][50:, :].value), rtol=0.02)
-    assert np.isclose(np.median(im2['data'][0:50, :].value), np.median(im2['data'][50:, :].value), rtol=0.02)
+    assert np.isclose(np.median(im1['data'][0:50, :]), np.median(im1['data'][50:, :]), rtol=0.02)
+    assert np.isclose(np.median(im2['data'][0:50, :]), np.median(im2['data'][50:, :]), rtol=0.02)
 
 
 def test_scaling():
@@ -533,13 +533,13 @@ def test_scaling():
         imdict['tabcatalog'], seed=rng_seed, effreadnoise=0)
 
     # check that sky level doesn't depend on pixel scale (in calibrated units!)
-    skyfracdiff = np.median(im1.data.value) / np.median(im2.data.value) - 1
+    skyfracdiff = np.median(im1.data) / np.median(im2.data) - 1
     log.info(f'skyfracdiff: {skyfracdiff:.3f}')
     assert np.abs(skyfracdiff) < 0.1
 
     # check that uncertainties match observed standard deviations
-    err1fracdiff = mad_std(im1.data.value) / np.median(im1.err.value) - 1
-    err2fracdiff = mad_std(im2.data.value) / np.median(im2.err.value) - 1
+    err1fracdiff = mad_std(im1.data) / np.median(im1.err) - 1
+    err2fracdiff = mad_std(im2.data) / np.median(im2.err) - 1
     log.info(f'err1fracdiff: {err1fracdiff:.3f}, '
              f'err2fracdiff: {err2fracdiff:.3f}')
     assert np.abs(err1fracdiff) < 0.1
@@ -551,18 +551,18 @@ def test_scaling():
         imdict['tabcatalog'], seed=rng_seed, effreadnoise=0)
 
     # check that sky level doesn't depend on exposure time (in calibrated units!)
-    sky3fracdiff = np.median(im1.data.value) / np.median(im3.data.value) - 1
+    sky3fracdiff = np.median(im1.data) / np.median(im3.data) - 1
     log.info(f'sky3fracdiff: {sky3fracdiff:.4f}')
     assert np.abs(sky3fracdiff) < 0.1
 
     # check that variances still work out
-    err3fracdiff = mad_std(im3.data.value) / np.median(im3.err.value) - 1
+    err3fracdiff = mad_std(im3.data) / np.median(im3.err) - 1
     log.info(f'err3fracdiff: {err3fracdiff:.3f}')
     assert np.abs(err3fracdiff) < 0.1
 
     # check that new variances are smaller than old ones by an appropriate factor
     errfracdiff = (
-        np.median(im1.err.value) / np.median(im3.err.value) - np.sqrt(10))
+        np.median(im1.err) / np.median(im3.err) - np.sqrt(10))
     log.info(f'err3 ratio diff from 1/sqrt(10) err1: {errfracdiff:0.3f}')
     assert np.abs(errfracdiff) < 0.1
 
@@ -574,7 +574,7 @@ def test_scaling():
             imdict['tabcatalog']['ra'][0], imdict['tabcatalog']['dec'][0])
         pind = [int(x) for x in pix]
         margin = 30 * fac
-        flux = np.sum(im.data.value[pind[1] - margin: pind[1] + margin,
+        flux = np.sum(im.data[pind[1] - margin: pind[1] + margin,
                                     pind[0] - margin: pind[0] + margin])
         fluxes.append(flux / fac ** 2)
         # division by fac ** 2 accounts for different pixel scale