Merge pull request #58 from PhangsTeam/phangs-jwst-v1p1

Update PHANGS Cy1 to v1.1
PhangsTeam · Mar 4, 2024 · b59226b · b59226b
2 parents d8303ec + ac6d633
commit b59226b
Show file tree

Hide file tree

Showing 14 changed files with 313 additions and 192 deletions.
diff --git a/CHANGES.rst b/CHANGES.rst
@@ -1,6 +1,16 @@
-1.0.5 (Unreleased)
-==================
+1.1 (Unreleased)
+================
 
+- If science lv1 data already exists, just copy it over for background data to save processing time
+- Include BMAJ in ``psf_matching_step``
+- Auto rotate the WCS in ``level_match_step``
+- Fix SIP approximation issues in ``apply_wcs_adjust_step``
+- ``astrometric_catalog_step`` can now use either DAOStarFinder or IRAFStarFinder. Observatory recommendations seem
+    to be IRAFStarFinder in general
+- Removed factor of 2 for FWHM in lv3 tweakreg-related steps to more closely match observatory recommendations
+- Add lower version pins to project requirements to avoid incompatibilities
+- Fixed bug in ``astrometric_align_step`` where file might already be closed before saving
+- Updated PHANGS Cycle 1 config to v1p1 (last couple of observations, and a few small changes to the config)
 - Fix bug with ``process_bgr_like_science`` in ``lv3_step``
 - Added 'filename' option for background association
 - Fix in ``parse_fits_to_table`` when there is not a defined observation label

diff --git a/config/2107/astronode.toml b/config/2107/astronode.toml
@@ -1,6 +1,6 @@
 crds_path = '/data/beegfs/astro-storage/groups/schinnerer/williams/crds/'
-crds_context = 'jwst_1147.pmap'
-raw_dir = '/data/beegfs/astro-storage/groups/schinnerer/williams/jwst_raw/2107/archive_20231015/'
+crds_context = 'jwst_1201.pmap'
+raw_dir = '/data/beegfs/astro-storage/groups/schinnerer/williams/jwst_raw/2107/archive_20240211/'
 reprocess_dir = '/data/beegfs/astro-storage/groups/schinnerer/williams/jwst_phangs_reprocessed/'
 alignment_dir = '/data/beegfs/astro-storage/groups/schinnerer/williams/pjpipe/config/2107/alignment/'
 webb_psf_data = '/data/beegfs/astro-storage/groups/schinnerer/williams/webbpsf-data'

diff --git a/config/2107/config.toml b/config/2107/config.toml
@@ -1,5 +1,6 @@
-# List of targets
+# List of targets. NGC7496 goes first, since IC5332 uses the background from those obs
 targets = [
+    'ngc7496',
     'ic5332',
     'ngc0628',
     'ngc1087',
@@ -18,11 +19,10 @@ targets = [
     'ngc4321',
     'ngc4535',
     'ngc5068',
-    'ngc7496',
 ]
 
 # Version for the reprocessing
-version = 'v1p0p1'
+version = 'v1p1'
 
 # Bands to consider
 bands = [
@@ -88,6 +88,12 @@ ngc1566 = 'jw02107007003_02105_00003*'
 # New IC5332 MIRI offs, which we need to skip
 ic5332.sci = 'jw02107042001_0?2*'
 
+# New NGC4303 MIRI offs, which we need to skip
+ngc4303.sci = 'jw02107044001_0?2*'
+
+# New NGC4321 MIRI offs, which we need to skip
+ngc4321.sci = 'jw02107043001_0?2*'
+
 # For IC5332, pull the backgrounds from NGC7496 but only for science images
 [parameters.move_raw_obs.extra_obs_to_include]
 ic5332.sci.ngc7496 = 'jw02107041001_0?2'
@@ -124,7 +130,7 @@ group_dithers = [
 
 [parameters.get_wcs_adjust.tweakreg_parameters]
 
-align_to_gaia = false
+starfinder = 'iraf'
 brightest = 500
 snr_threshold = 3
 expand_refcat = true
@@ -134,11 +140,11 @@ roundhi.nircam = 0.5
 
 # Parameters to get wcs_adjust shifts
 minobj = 3
+sigma = 2
 searchrad = 2
 separation.miri = 1
-tolerance.miri = 0.7
 separation.nircam = 2
-tolerance.nircam = 1
+tolerance = 0.3
 use2dhist = true
 fitgeometry = 'shift'
 nclip = 5
@@ -199,7 +205,7 @@ skip.F1000W = true
 skip.F1130W = true
 skip.F2100W = true
 
-align_to_gaia = false
+starfinder = 'iraf'
 brightest.nircam_short = 125
 brightest.nircam_long = 500
 brightest.miri = 500
@@ -215,6 +221,7 @@ roundhi.nircam = 0.5
 # close to the correct solution
 separation.miri = 1
 separation.nircam = 2
+sigma = 2
 tolerance.nircam_short = 0.3
 tolerance.nircam_long = 0.3
 tolerance.miri = 1
@@ -253,8 +260,9 @@ deblend = true
 
 [parameters.astrometric_catalog]
 snr = 10
+starfind_method = "iraf"
 
-[parameters.astrometric_catalog.dao_parameters]
+[parameters.astrometric_catalog.starfind_parameters]
 sharplo = 0.2
 sharphi = 1.0
 roundlo = -0.5
@@ -354,4 +362,4 @@ move_psf_matched = true
 move_diagnostic_plots = true
 
 [parameters.regress_against_previous]
-prev_version = 'v1p0'
+prev_version = 'v1p0p1'
diff --git a/config/2107/mac.toml b/config/2107/mac.toml
@@ -1,5 +1,5 @@
 crds_path = '/Users/thomaswilliams/Documents/phangs/jwst/crds/'
 raw_dir = '/Users/thomaswilliams/Documents/phangs/jwst/jwst_raw/'
 reprocess_dir = '/Users/thomaswilliams/Documents/phangs/jwst/jwst_phangs_reprocessed/'
-alignment_dir = '/Users/thomaswilliams/PycharmProjects/jwst_scripts/config/2107/alignment/2107/alignment/'
+alignment_dir = '/Users/thomaswilliams/PycharmProjects/pjpipe/config/2107/alignment/2107/alignment/'
 webb_psf_data = '/Users/thomaswilliams/Documents/phangs/jwst/webb_psf_data/'
diff --git a/pjpipe/apply_wcs_adjust/apply_wcs_adjust_step.py b/pjpipe/apply_wcs_adjust/apply_wcs_adjust_step.py
@@ -1,3 +1,4 @@
+import copy
 import gc
 import glob
 import logging
@@ -146,10 +147,18 @@ def parallel_wcs_adjust(
 
         # Set up the WCSCorrector per tweakreg
         with datamodels.open(file) as input_im:
-            ref_wcs = input_im.meta.wcs
-            ref_wcsinfo = input_im.meta.wcsinfo.instance
 
-            im = JWSTWCSCorrector(ref_wcs, ref_wcsinfo)
+            model_name = os.path.splitext(input_im.meta.filename)[0].strip('_- ')
+
+            refang = input_im.meta.wcsinfo.instance
+            im = JWSTWCSCorrector(
+                wcs=input_im.meta.wcs,
+                wcsinfo={'roll_ref': refang['roll_ref'],
+                         'v2_ref': refang['v2_ref'],
+                         'v3_ref': refang['v3_ref']},
+                meta={'image_model': input_im,
+                      'name': model_name},
+            )
 
             # Pull out the info we need to shift. If we have both
             # dithers ungrouped and grouped, prefer the ungrouped
@@ -179,27 +188,38 @@ def parallel_wcs_adjust(
                         visit_found = True
 
             if not visit_found:
-                log.info(f"No shifts found for {file_short}. Defaulting to no shift")
+                log.info(f"No shifts found for {file_short}. Will write out without shifting")
+
+            if visit_found:
+
+                im.set_correction(matrix=matrix, shift=shift)
+
+                image_model = im.meta["image_model"]
+                image_model.meta.wcs = im.wcs
 
-            im.set_correction(matrix=matrix, shift=shift)
+                with warnings.catch_warnings():
+                    warnings.simplefilter("ignore")
+                    try:
+                        update_fits_wcsinfo(
+                            image_model,
+                            degree=6,
+                            max_pix_error=0.01,
+                            npoints=32,
+                        )
+                    except (ValueError, RuntimeError) as e:
+                        log.warning(
+                            "Failed to update 'meta.wcsinfo' with FITS SIP "
+                            f"approximation. Reported error is:\n'{e.args[0]}'"
+                        )
 
-            input_im.meta.wcs = im.wcs
+            else:
 
-            with warnings.catch_warnings():
-                warnings.simplefilter("ignore")
-                try:
-                    update_fits_wcsinfo(
-                        input_im,
-                    )
-                except (ValueError, RuntimeError) as e:
-                    log.warning(
-                        "Failed to update 'meta.wcsinfo' with FITS SIP "
-                        f"approximation. Reported error is:\n'{e.args[0]}'"
-                    )
+                image_model = copy.deepcopy(im.meta["image_model"])
 
-            input_im.save(output_file)
+            image_model.save(output_file)
 
         del input_im
+        del image_model
         del im
         gc.collect()
 

diff --git a/pjpipe/astrometric_align/astrometric_align_step.py b/pjpipe/astrometric_align/astrometric_align_step.py
@@ -156,8 +156,9 @@ def parallel_tweakback(
         try:
             update_fits_wcsinfo(
                 input_im,
+                degree=6,
                 max_pix_error=0.01,
-                npoints=16,
+                npoints=32,
             )
         except (ValueError, RuntimeError) as e:
             logging.warning(
@@ -363,81 +364,82 @@ def align_to_aligned_image(
                 shift = shift.astype(np.ndarray).astype(float)
                 matrix = matrix.astype(np.ndarray).astype(float)
 
-            with datamodels.open(file) as target_im:
-                target_wcs = get_lv3_wcs(target_im)
-                target_wcs_corrector = FITSWCSCorrector(target_wcs)
-                target_wcs_corrector_orig = copy.deepcopy(target_wcs_corrector)
-                target_data = copy.deepcopy(target_im.data)
-                target_err = copy.deepcopy(target_im.err)
-                target_data[target_data == 0] = np.nan
-
-            if self.align_mapping_mode == "cross_corr":
-                with warnings.catch_warnings():
-                    warnings.simplefilter("ignore")
-                    ref_data = reproject_interp(
-                        (ref_data, ref_wcs),
-                        target_wcs,
-                        shape_out=target_data.shape,
-                        return_footprint=False,
-                    )
-
-                    ref_err = reproject_interp(
-                        (ref_err, ref_wcs),
-                        target_wcs,
-                        shape_out=target_data.shape,
-                        return_footprint=False,
-                    )
-
-                nan_idx = np.logical_or(np.isnan(ref_data), np.isnan(target_data))
-
-                ref_data[nan_idx] = np.nan
-                target_data[nan_idx] = np.nan
+                with datamodels.open(file) as target_im:
+                    target_wcs = get_lv3_wcs(target_im)
+                    target_wcs_corrector = FITSWCSCorrector(target_wcs)
+                    target_wcs_corrector_orig = copy.deepcopy(target_wcs_corrector)
+                    target_data = copy.deepcopy(target_im.data)
+                    target_err = copy.deepcopy(target_im.err)
+                    target_data[target_data == 0] = np.nan
+
+                    if self.align_mapping_mode == "cross_corr":
+                        with warnings.catch_warnings():
+                            warnings.simplefilter("ignore")
+                            ref_data = reproject_interp(
+                                (ref_data, ref_wcs),
+                                target_wcs,
+                                shape_out=target_data.shape,
+                                return_footprint=False,
+                            )
+
+                            ref_err = reproject_interp(
+                                (ref_err, ref_wcs),
+                                target_wcs,
+                                shape_out=target_data.shape,
+                                return_footprint=False,
+                            )
+
+                        nan_idx = np.logical_or(np.isnan(ref_data), np.isnan(target_data))
+
+                        ref_data[nan_idx] = np.nan
+                        target_data[nan_idx] = np.nan
+
+                        ref_err[nan_idx] = np.nan
+                        target_err[nan_idx] = np.nan
+
+                        # Make sure we're square, since apparently this causes weirdness
+                        data_size_min = min(target_data.shape)
+                        data_slice_i = slice(
+                            target_data.shape[0] // 2 - data_size_min // 2,
+                            target_data.shape[0] // 2 + data_size_min // 2,
+                        )
+                        data_slice_j = slice(
+                            target_data.shape[1] // 2 - data_size_min // 2,
+                            target_data.shape[1] // 2 + data_size_min // 2,
+                        )
 
-                ref_err[nan_idx] = np.nan
-                target_err[nan_idx] = np.nan
+                        x_off, y_off = cross_correlation_shifts(
+                            ref_data[data_slice_i, data_slice_j],
+                            target_data[data_slice_i, data_slice_j],
+                            errim1=ref_err[data_slice_i, data_slice_j],
+                            errim2=target_err[data_slice_i, data_slice_j],
+                        )
+                        shift = [-x_off, -y_off]
+                        matrix = [[1, 0], [0, 1]]
 
-                # Make sure we're square, since apparently this causes weirdness
-                data_size_min = min(target_data.shape)
-                data_slice_i = slice(
-                    target_data.shape[0] // 2 - data_size_min // 2,
-                    target_data.shape[0] // 2 + data_size_min // 2,
-                )
-                data_slice_j = slice(
-                    target_data.shape[1] // 2 - data_size_min // 2,
-                    target_data.shape[1] // 2 + data_size_min // 2,
-                )
+                        log.info(f"Found offset of {shift}")
 
-                x_off, y_off = cross_correlation_shifts(
-                    ref_data[data_slice_i, data_slice_j],
-                    target_data[data_slice_i, data_slice_j],
-                    errim1=ref_err[data_slice_i, data_slice_j],
-                    errim2=target_err[data_slice_i, data_slice_j],
-                )
-                shift = [-x_off, -y_off]
-                matrix = [[1, 0], [0, 1]]
+                    elif self.align_mapping_mode == "shift":
+                        # Add in shift metadata
+                        target_im.meta.abs_astro_alignment = {
+                            "shift": shift,
+                            "matrix": matrix,
+                        }
 
-                log.info(f"Found offset of {shift}")
-
-            elif self.align_mapping_mode == "shift":
-                # Add in shift metadata
-                target_im.meta.abs_astro_alignment = {
-                    "shift": shift,
-                    "matrix": matrix,
-                }
+                    # Apply correction
+                    target_wcs_corrector.set_correction(
+                        shift=shift,
+                        matrix=matrix,
+                        ref_tpwcs=target_wcs_corrector_orig,
+                    )
 
-            # Apply correction
-            target_wcs_corrector.set_correction(
-                shift=shift,
-                matrix=matrix,
-                ref_tpwcs=target_wcs_corrector_orig,
-            )
+                    target_hdr, new_gwcs = transform_wcs_gwcs(target_wcs_corrector.wcs)
+                    target_im.meta.wcs = new_gwcs
 
-            target_hdr, new_gwcs = transform_wcs_gwcs(target_wcs_corrector.wcs)
-            target_im.meta.wcs = new_gwcs
+                    # Update WCS info
+                    updated_im = lv3_update_fits_wcsinfo(im=target_im, hdr=target_hdr)
 
-            # Update WCS info
-            updated_im = lv3_update_fits_wcsinfo(im=target_im, hdr=target_hdr)
-            updated_im.write(aligned_file)
+                    updated_im.write(aligned_file)
 
             # Also apply this to each individual crf file
             crf_files = glob.glob(