Update romanisim to use the new L1 / L2 metadata (#166)

* Remove units.

* Update pyproject.toml.

* Remove ggsaa from file names.

* change to lowercase

* Point rad/rdm to main.

* Point rad/rdm to main better?

* Increase minimum asdf version to match rdm.

* Increase minimum gwcs version to match rdm.

* Updates to support L1/L2 metadata refactoring

* Make unit changes less dumb.

* Metadata and filename changes.

* Point rad/rdm to L1/L2 refactor branches.

* Fix rebase failure.

* fix pyproject.toml syntax.

* Make program an integer.

* Update wfi lowercase and docs.

* Misc L3 updates.

* Repoint rad/rdm to main.

Author: schlafly

docs/romanisim/running.rst

@@ -15,37 +15,45 @@ The ``romanisim-make-image`` command line interface (CLI) has a number of argume
 this functionality::
 
     romanisim-make-image -h
-    usage: romanisim-make-image [-h] [--catalog CATALOG] [--radec RADEC RADEC] [--bandpass BANDPASS]
-                                [--sca SCA] [--usecrds] [--webbpsf] [--date DATE [DATE ...]]
-                                [--level LEVEL] [--ma_table_number MA_TABLE_NUMBER] [--seed SEED]
-                                [--nobj NOBJ] [--boresight] [--previous PREVIOUS]
-                                filename
-    
+    usage: romanisim-make-image [-h] [--bandpass BANDPASS] [--boresight] [--catalog CATALOG] [--config CONFIG]
+				[--date DATE] [--level LEVEL] [--ma_table_number MA_TABLE_NUMBER] [--nobj NOBJ]
+				[--previous PREVIOUS] [--radec RADEC RADEC] [--rng_seed RNG_SEED] [--roll ROLL]
+				[--sca SCA] [--usecrds] [--webbpsf] [--truncate TRUNCATE]
+				[--pretend-spectral PRETEND_SPECTRAL] [--drop-extra-dq]
+				filename
+
     Make a demo image.
-    
+
     positional arguments:
-      filename              output image (fits)
-    
-    optional arguments:
+      filename              output image (asdf)
+
+    options:
       -h, --help            show this help message and exit
-      --catalog CATALOG     input catalog (csv) (default: None)
-      --radec RADEC RADEC   ra and dec (deg) (default: None)
       --bandpass BANDPASS   bandpass to simulate (default: F087)
-      --sca SCA             SCA to simulate (default: 7). Use -1 to generate images for all SCAs.
-      --usecrds             Use CRDS for distortion map (default: False)
-      --webbpsf             Use webbpsf for PSF (default: False)
-      --date DATE [DATE ...]
-                            Date of observation to simulate: year month day hour minute second
-                            microsecond (default: None)
+      --boresight           radec specifies location of boresight, not center of WFI. (default: False)
+      --catalog CATALOG     input catalog (ecsv) (default: None)
+      --config CONFIG       input parameter override file (yaml) (default: None)
+      --date DATE           UTC Date and Time of observation to simulate in ISOT format. (default: None)
       --level LEVEL         1 or 2, for L1 or L2 output (default: 2)
       --ma_table_number MA_TABLE_NUMBER
-      --rng_seed SEED
       --nobj NOBJ
-      --boresight           radec specifies location of boresight, not center of WFI. (default: False)
-      --previous PREVIOUS   previous simulated file in chronological order used for persistence modeling.
-                            (default: None)
-    
-    EXAMPLE: romanisim-make-image output_image.fits
+      --previous PREVIOUS   previous simulated file in chronological order used for persistence modeling. (default:
+			    None)
+      --radec RADEC RADEC   ra and dec (deg) (default: None)
+      --rng_seed RNG_SEED
+      --roll ROLL           Position angle (North towards YIdl) measured at the V2Ref/V3Ref of the aperture used.
+			    (default: 0)
+      --sca SCA             SCA to simulate. Use -1 to generate images for all SCAs; include {} in filename for this
+			    mode to indicate where the SCA number should be filled, e.g. l1_wfi{}.asdf (default: 7)
+      --usecrds             Use CRDS for distortion map (default: False)
+      --webbpsf             Use webbpsf for PSF (default: False)
+      --truncate TRUNCATE   If set, truncate the MA table at given number of resultants. (default: None)
+      --pretend-spectral PRETEND_SPECTRAL
+			    Pretend the image is spectral. exposure.type and instrument.element are updated to be
+			    grism / prism. (default: None)
+      --drop-extra-dq       Do not store the optional simulated dq array. (default: False)
+
+    EXAMPLE: romanisim-make-image output_image.asdf
 
 Expected arguments controlling things like the input :doc:`here </romanisim/catalog>` to
 simulate, the right ascension and declination of the telescope

romanisim/image.py

@@ -259,7 +259,7 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf,
     if rng is None:
         rng = galsim.UniformDeviate(seed)
 
-    log.info('Adding sources to image...')
+    log.info(f'Adding {len(objlist)} sources to image...')
     nrender = 0
 
     chromatic = False
@@ -791,7 +791,7 @@ def simulate(metadata, objlist,
     counts, simcatobj = simulate_counts(
         image_mod.meta, objlist, rng=rng, usecrds=usecrds, darkrate=darkrate,
         webbpsf=webbpsf, flat=flat, psf_keywords=psf_keywords)
-    util.update_aperture_and_wcsinfo_metadata(image_mod.meta, counts.wcs)
+    util.update_pointing_and_wcsinfo_metadata(image_mod.meta, counts.wcs)
     if level == 0:
         im = dict(data=counts.array, meta=dict(image_mod.meta.items()))
     else:

romanisim/l3.py

@@ -534,6 +534,7 @@ def simulate_cps(image, filter_name, efftimes, objlist=None, psf=None,
         objlist = []
     if len(objlist) > 0 and xpos is None:
         if isinstance(objlist, table.Table):
+            objlist = romanisim.image.trim_objlist(objlist, image)
             coord = np.array([[o['ra'], o['dec']] for o in objlist])
         else:
             coord = np.array([[o.sky_pos.ra.deg, o.sky_pos.dec.deg]
@@ -730,9 +731,8 @@ def add_more_metadata(metadata, efftimes, filter_name, wcs, shape, nexposures):
     metadata['basic']['individual_image_meta'] = None
     metadata['model_type'] = 'WfiMosaic'
     metadata['photometry']['conversion_microjanskys'] = (
-        (1e12 * (u.rad / u.arcsec) ** 2).to(u.dimensionless_unscaled) *
-        u.uJy / u.arcsec ** 2)
-    metadata['photometry']['conversion_megajanskys'] = 1 * u.MJy / u.sr
+        (1e12 * (u.rad / u.arcsec) ** 2).to(u.dimensionless_unscaled))
+    metadata['photometry']['conversion_megajanskys'] = 1
 
     cenx, ceny = ((shape[1] - 1) / 2, (shape[0] - 1) / 2)
     c1 = wcs.pixel_to_world(cenx, ceny)
@@ -742,10 +742,8 @@ def add_more_metadata(metadata, efftimes, filter_name, wcs, shape, nexposures):
     metadata['photometry']['pixelarea_steradians'] = (pscale ** 2).to(u.sr)
     metadata['photometry']['pixelarea_arcsecsq'] = (
         pscale.to(u.arcsec) ** 2)
-    metadata['photometry']['conversion_microjanskys_uncertainty'] = (
-        0 * u.uJy / u.arcsec ** 2)
-    metadata['photometry']['conversion_megajanskys_uncertainty'] = (
-        0 * u.MJy / u.sr)
+    metadata['photometry']['conversion_microjanskys_uncertainty'] = 0
+    metadata['photometry']['conversion_megajanskys_uncertainty'] = 0
     metadata['resample']['pixel_scale_ratio'] = (
         pscale.to(u.arcsec).value / romanisim.parameters.pixel_scale)
     metadata['resample']['pixfrac'] = 0

romanisim/parameters.py

@@ -40,8 +40,9 @@
                  'ma_table_number': 1,
                  'read_pattern': read_pattern[1],
                  },
-    'pointing': {'ra_v1': 270.0,
-                 'dec_v1': 66.0,
+    'pointing': {'target_ra': 270.0,
+                 'target_dec': 66.0,
+                 'target_aperture': 'WFI_CEN',
                  },
     'wcsinfo': {'ra_ref': 270.0,
                 'dec_ref': 66.0,
@@ -53,9 +54,6 @@
                 # I don't know what vparity and v3yangle should really be,
                 # but they are always -1 and -60 in existing files.
                 },
-    'aperture': {'name': 'WFI_CEN',
-                 'position_angle': 0
-                 },
 }
 
 # Default metadata for level 3 mosaics

romanisim/ris_make_utils.py

@@ -202,9 +202,9 @@ def parse_filename(filename):
     match = pattern.match(filename)
     if match is None:
         return None
-    out = dict(obs_id=filename.replace('_', '')[1:],
+    out = dict(observation_id=filename.replace('_', '')[1:],
                visit_id=filename[1:20],
-               program=match.group(1),  # this one is a string
+               program=int(match.group(1)),  # this one is a string
                execution_plan=int(match.group(2)),
                # pass = int(match.group(3))
                segment=int(match.group(4)),
@@ -259,8 +259,8 @@ def simulate_image_file(args, metadata, cat, rng=None, persist=None):
             'WFI_' + args.pretend_spectral.upper())
         im['meta']['instrument']['optical_element'] = (
             args.pretend_spectral.upper())
-        im['meta']['guidestar']['gw_window_xsize'] = 170
-        im['meta']['guidestar']['gw_window_ysize'] = 24
+        im['meta']['guide_star']['window_xsize'] = 170
+        im['meta']['guide_star']['window_ysize'] = 24
 
     drop_extra_dq = getattr(args, 'drop_extra_dq', False)
     if drop_extra_dq and ('dq' in romanisimdict):

romanisim/tests/test_ris_make_utils.py

@@ -55,5 +55,5 @@ def test_parse_filename():
     obs = ris_make_utils.parse_filename(
         'r9999901001001001001_01101_0001_uncal.asdf')
     assert obs is not None
-    assert obs['program'] == '99999'
+    assert obs['program'] == 99999
     assert obs['pass'] == 1

romanisim/tests/test_wcs.py

@@ -98,10 +98,10 @@ def test_wcs():
     # should be close to the reference v2 & v3 offset.
     assert np.abs(cc3.separation(cc4).to(u.arcsec).value
                   - np.hypot(*parameters.v2v3_wficen)) < 1
-    metadata = dict(aperture=dict(), instrument=dict(), wcsinfo=dict())
+    metadata = dict(pointing=dict(), instrument=dict(), wcsinfo=dict())
     metadata['instrument']['detector'] = 'WFI01'
-    util.update_aperture_and_wcsinfo_metadata(metadata, wcs.GWCS(gwcs))
-    assert metadata['aperture']['name'] == 'WFI_01_FULL'
+    util.update_pointing_and_wcsinfo_metadata(metadata, wcs.GWCS(gwcs))
+    assert metadata['wcsinfo']['aperture_name'] == 'WFI01_FULL'
 
 
 def test_wcs_from_fits_header():

romanisim/util.py

@@ -189,8 +189,8 @@ def add_more_metadata(metadata):
 
     if 'exposure' not in metadata.keys():
         metadata['exposure'] = {}
-    if 'guidestar' not in metadata.keys():
-        metadata['guidestar'] = {}
+    if 'guide_star' not in metadata.keys():
+        metadata['guide_star'] = {}
     read_pattern = metadata['exposure'].get(
         'read_pattern',
         parameters.read_pattern[metadata['exposure']['ma_table_number']])
@@ -204,73 +204,32 @@ def add_more_metadata(metadata):
     for prefix, offset in offsets.items():
         metadata['exposure'][f'{prefix}_time'] = Time((
             starttime + offset).isot)
-        metadata['exposure'][f'{prefix}_time_mjd'] = (
-            starttime + offset).mjd
-        metadata['exposure'][f'{prefix}_time_tdb'] = (
-            starttime + offset).tdb.mjd
-    metadata['exposure']['ngroups'] = len(read_pattern)
-    metadata['exposure']['sca_number'] = (
-        int(metadata['instrument']['detector'][-2:]))
-    metadata['exposure']['integration_time'] = openshuttertime
-    metadata['exposure']['elapsed_exposure_time'] = openshuttertime
-    # ???
-    metadata['exposure']['groupgap'] = 0
+    metadata['exposure']['nresultants'] = len(read_pattern)
     metadata['exposure']['frame_time'] = parameters.read_time
     metadata['exposure']['exposure_time'] = openshuttertime
-    metadata['exposure']['effective_exposure_time'] = openshuttertime
-    metadata['exposure']['duration'] = openshuttertime
-    metadata['guidestar']['gw_window_xsize'] = 16
-    metadata['guidestar']['gw_window_ysize'] = 16
-    if 'gw_window_xstart' in metadata['guidestar']:
-        metadata['guidestar']['gw_window_xstop'] = (
-            metadata['guidestar']['gw_window_xstart'])
-        metadata['guidestar']['gw_window_ystop'] = (
-            metadata['guidestar']['gw_window_ystart'])
-    # integration_start?  integration_end?  nints = 1?  ...
-
-    if 'target' not in metadata.keys():
-        metadata['target'] = {}
-    target = metadata['target']
-    target['type'] = 'FIXED'
-    if 'wcsinfo' in metadata.keys():
-        target['ra'] = metadata['wcsinfo']['ra_ref']
-        target['dec'] = metadata['wcsinfo']['dec_ref']
-        target['proposer_ra'] = target['ra']
-        target['proposer_dec'] = target['dec']
-    target['ra_uncertainty'] = 0
-    target['dec_uncertainty'] = 0
-    target['proper_motion_ra'] = 0
-    target['proper_motion_dec'] = 0
-    target['proper_motion_epoch'] = 'J2000'
-    target['source_type'] = 'EXTENDED'
-
-    # there are a few metadata keywords that have problematic, too-long
-    # defaults in RDM.
-    # program.category
-    # ephemeris.ephemeris_reference_frame
-    # guidestar.gs_epoch
-    # this truncates these to the maximum allowed characters.  Alternative
-    # solutions would include doing things like:
-    #   making the roman_datamodels defaults archivable
-    #   making the roman_datamodels validation check lengths of strings
-    if 'program' in metadata:
-        metadata['program']['category'] = metadata['program']['category'][:6]
-    if 'ephemeris' in metadata:
-        metadata['ephemeris']['ephemeris_reference_frame'] = (
-            metadata['ephemeris']['ephemeris_reference_frame'][:10])
-    if 'guidestar' in metadata and 'gs_epoch' in metadata['guidestar']:
-        metadata['guidestar']['gs_epoch'] = (
-            metadata['guidestar']['gs_epoch'][:10])
-
-
-def update_aperture_and_wcsinfo_metadata(metadata, gwcs):
-    """Update aperture and wcsinfo keywords to use the aperture for this SCA.
+    effexptime = parameters.read_time * (
+        np.mean(read_pattern[-1]) - np.mean(read_pattern[0]))
+    metadata['exposure']['effective_exposure_time'] = effexptime
+    metadata['guide_star']['window_xsize'] = 16
+    metadata['guide_star']['window_ysize'] = 16
+    if 'window_xstart' in metadata['guide_star']:
+        metadata['guide_star']['window_xstop'] = (
+            metadata['guide_star']['window_xstart'])
+        metadata['guide_star']['window_ystop'] = (
+            metadata['guide_star']['window_ystart'])
+    if 'visit' not in metadata.keys():
+        metadata['visit'] = dict()
+    metadata['visit']['status'] = 'SUCCESSFUL'
+
+
+def update_pointing_and_wcsinfo_metadata(metadata, gwcs):
+    """Update pointing and wcsinfo keywords to use the aperture for this SCA.
 
     Updates metadata in place, setting v2/v3_ref to be equal to the V2 and V3 of
     the center of the detector, and ra/dec_ref accordingly.  Also updates the
-    aperture to refer to this SCA.
+    pointing to refer to this SCA and ra/dec_v1 to point along the boresight.
 
-    No updates are  performed if gwcs is not a gWCS object or if aperture and
+    No updates are  performed if gwcs is not a gWCS object or if pointing and
     wcsinfo are not present in metadata.
 
     Parameters
@@ -280,15 +239,14 @@ def update_aperture_and_wcsinfo_metadata(metadata, gwcs):
     gwcs : WCS object
         image WCS
     """
-    if 'aperture' not in metadata or 'wcsinfo' not in metadata:
+    if 'pointing' not in metadata or 'wcsinfo' not in metadata:
         return
     if isinstance(gwcs, wcs.GWCS):
         gwcs = gwcs.wcs
     if not isinstance(gwcs, gwcsmod.wcs.WCS):
         return
-    metadata['aperture']['name'] = (
-        metadata['instrument']['detector'][:3] + '_'
-        + metadata['instrument']['detector'][3:] + '_FULL')
+    metadata['wcsinfo']['aperture_name'] = (
+        metadata['instrument']['detector'] + '_FULL')
     distortion = gwcs.get_transform('detector', 'v2v3')
     center = (galsim.roman.n_pix / 2 - 0.5, galsim.roman.n_pix / 2 - 0.5)
     v2v3 = distortion(*center)
@@ -314,6 +272,16 @@ def update_aperture_and_wcsinfo_metadata(metadata, gwcs):
     metadata['wcsinfo']['v3_ref'] = v2v3[1]
     metadata['wcsinfo']['roll_ref'] = roll_ref
 
+    boresight = t2sky(0, 0)
+    metadata['pointing']['ra_v1'] = boresight[0]
+    metadata['pointing']['dec_v1'] = boresight[1]
+    boresightn = t2sky(0, 1)
+    pa_v3 = (
+        SkyCoord(boresight[0] * u.deg, boresight[1] * u.deg).position_angle(
+        SkyCoord(boresightn[0] * u.deg, boresightn[1] * u.deg)))
+    pa_v3 = pa_v3.to(u.deg).value
+    metadata['pointing']['pa_v3'] = pa_v3
+
 
 def king_profile(r, rc, rt):
     """Compute the King profile.
@@ -532,12 +500,10 @@ def update_photom_keywords(im, gain=None):
                  cc[0].position_angle(cc[2]))
         area = (cc[0].separation(cc[1]) * cc[0].separation(cc[2])
                 * np.sin(angle.to(u.rad).value))
-        im['meta']['photometry']['pixelarea_steradians'] = area.to(u.sr).value
-        im['meta']['photometry']['pixelarea_arcsecsq'] = (
-            area.to(u.arcsec ** 2)).value
+        im['meta']['photometry']['pixel_area'] = area.to(u.sr).value
         val = (gain * (3631 / bandpass.get_abflux(
-            im.meta['instrument']['optical_element']) /
-            10 ** 6 / im['meta']['photometry']['pixelarea_steradians']))
+             im.meta['instrument']['optical_element']) /
+             10 ** 6 / im['meta']['photometry']['pixel_area']))
         im['meta']['photometry']['conversion_megajanskys'] = val
         im['meta']['photometry']['conversion_microjanskys'] = (
             val * u.MJy / u.sr).to(u.uJy / u.arcsec ** 2).value

romanisim/wcs.py

@@ -59,20 +59,14 @@ def fill_in_parameters(parameters, coord, pa_aper=0, boresight=True):
 
     if 'pointing' not in parameters.keys():
         parameters['pointing'] = {}
-
-    parameters['pointing']['ra_v1'] = coord.ra.to(u.deg).value
-
     if 'wcsinfo' not in parameters.keys():
         parameters['wcsinfo'] = {}
-    if 'aperture' not in parameters.keys():
-        parameters['aperture'] = {}
 
-    parameters['wcsinfo']['ra_ref'] = (
-        parameters['pointing']['ra_v1'])
+    parameters['pointing']['target_ra'] = coord.ra.to(u.deg).value
+    parameters['wcsinfo']['ra_ref'] = parameters['pointing']['target_ra']
 
-    parameters['pointing']['dec_v1'] = coord.dec.to(u.deg).value
-    parameters['wcsinfo']['dec_ref'] = (
-        parameters['pointing']['dec_v1'])
+    parameters['pointing']['target_dec'] = coord.dec.to(u.deg).value
+    parameters['wcsinfo']['dec_ref'] = parameters['pointing']['target_dec']
 
     # Romanisim uses ROLL_REF = PA_APER - V3IdlYAngle
     parameters['wcsinfo']['roll_ref'] = (
@@ -81,14 +75,16 @@ def fill_in_parameters(parameters, coord, pa_aper=0, boresight=True):
     if boresight:
         parameters['wcsinfo']['v2_ref'] = 0
         parameters['wcsinfo']['v3_ref'] = 0
-        parameters['aperture']['name'] = 'BORESIGHT'
+        parameters['wcsinfo']['aperture_name'] = 'BORESIGHT'
+        parameters['pointing']['target_aperture'] = 'BORESIGHT'
     else:
         from .parameters import v2v3_wficen
         parameters['wcsinfo']['v2_ref'] = v2v3_wficen[0]
         parameters['wcsinfo']['v3_ref'] = v2v3_wficen[1]
-        parameters['aperture']['name'] = 'WFI_CEN'
+        parameters['wcsinfo']['aperture_name'] = 'WFI_CEN'
+        parameters['pointing']['target_aperture'] = 'WFI_CEN'
 
-    parameters['aperture']['position_angle'] = pa_aper
+    parameters['wcsinfo']['pa_aperture'] = pa_aper
 
 
 def get_wcs(image, usecrds=True, distortion=None):