JP-2076 moved wrap_ra to assign_wcs

CHANGES.rst

@@ -9,6 +9,9 @@ ami_analyze
 assign_wcs
 ----------
 
+- Moved the routine wrap_ra from cube_build to assign_wcs.ulit. The correctly
+  determines for s_region for data that cross ra boundary. [#6026]
+
 - Changed evaluation of grism bounding box center from averaged extrema of
   transformed bounding box to transformed centroid of source_cat object [#5809]
 

jwst/assign_wcs/util.py

@@ -32,7 +32,7 @@
 
 __all__ = ["reproject", "wcs_from_footprints", "velocity_correction",
            "MSAFileError", "NoDataOnDetectorError", "compute_scale",
-           "calc_rotation_matrix"]
+           "calc_rotation_matrix", "wrap_ra"]
 
 
 class MSAFileError(Exception):
@@ -883,10 +883,22 @@ def compute_footprint_spectral(model):
 
     x, y = grid_from_bounding_box(bbox)
     ra, dec, lam = swcs(x, y)
-    footprint = np.array([[np.nanmin(ra), np.nanmin(dec)],
-                          [np.nanmax(ra), np.nanmin(dec)],
-                          [np.nanmax(ra), np.nanmax(dec)],
-                          [np.nanmin(ra), np.nanmax(dec)]])
+
+    # the wrapped ra values are forced to be on one side of ra-boarder
+    # the wrapped ra are used to determine the correct  min and max ra
+    ra = wrap_ra(ra)
+    min_ra = np.nanmin(ra)
+    max_ra = np.nanmax(ra)
+
+    # for the footprint we want the ra values to fall between 0 to 360
+    if(min_ra < 0):
+        min_ra = min_ra + 360.0
+    if(max_ra > 360.0):
+        max_ra = max_ra - 360.0
+    footprint = np.array([[min_ra, np.nanmin(dec)],
+                          [max_ra, np.nanmin(dec)],
+                          [max_ra, np.nanmax(dec)],
+                          [min_ra, np.nanmax(dec)]])
     lam_min = np.nanmin(lam)
     lam_max = np.nanmax(lam)
     return footprint, (lam_min, lam_max)
@@ -1002,8 +1014,17 @@ def compute_footprint_nrs_ifu(dmodel, mod):
         ra_total.extend(np.ravel(ra))
         dec_total.extend(np.ravel(dec))
         lam_total.extend(np.ravel(lam))
+    # the wrapped ra values are forced to be on one side of ra-boa`rder
+    # the wrapped ra are used to determine the correct  min and max ra
+    ra_total = wrap_ra(ra_total)
     ra_max = np.nanmax(ra_total)
     ra_min = np.nanmin(ra_total)
+    # for the footprint we want ra to be between 0 to 360
+    if(ra_min < 0):
+        ra_min = ra_min + 360.0
+    if(ra_max > 360.0):
+        ra_max = ra_max - 360.0
+
     dec_max = np.nanmax(dec_total)
     dec_min = np.nanmin(dec_total)
     lam_max = np.nanmax(lam_total)
@@ -1056,3 +1077,38 @@ def velocity_correction(velosys):
     model.inverse = astmodels.Identity(1) / astmodels.Const1D(correction, name="inv_vel_correciton")
 
     return model
+
+
+def wrap_ra(ravalues):
+    """Test for 0/360 wrapping in ra values.
+
+    If exists it makes it difficult to determine
+    ra range of a region on the sky. This problem is solved by putting them all
+    on "one side" of 0/360 border
+
+    Parameters
+    ----------
+    input : ravalues
+        RA values numpy.ndarray
+
+    Returns
+    ------
+    a numpy array of ra values all on "same side" of 0/360 border
+    """
+
+    index_good = np.where(np.isfinite(ravalues))
+    ravalues_wrap = ravalues[index_good].copy()
+    median_ra = np.nanmedian(ravalues_wrap)
+
+    # using median to test if there is any wrapping going on
+    wrap_index = np.where(np.fabs(ravalues_wrap - median_ra) > 180.0)
+    nwrap = wrap_index[0].size
+
+    # get all the ra on the same "side" of 0/360
+    if(nwrap != 0 and median_ra < 180):
+        ravalues_wrap[wrap_index] = ravalues_wrap[wrap_index] - 360.0
+
+    if(nwrap != 0 and median_ra > 180):
+        ravalues_wrap[wrap_index] = ravalues_wrap[wrap_index] + 360.0
+
+    return ravalues_wrap

jwst/cube_build/cube_build_wcs_util.py

@@ -2,6 +2,7 @@
 """
 import numpy as np
 from ..assign_wcs import nirspec
+from ..assign_wcs.util import wrap_ra
 from gwcs import wcstools
 import logging
 
@@ -91,6 +92,7 @@ def find_corners_MIRI(input, this_channel, instrument_info, coord_system):
     a_min = np.nanmin(coord1)
     a_max = np.nanmax(coord1)
 
+    # find index of min a value
     a1_index = np.argmin(coord1)
     a2_index = np.argmax(coord1)
 
@@ -100,6 +102,7 @@ def find_corners_MIRI(input, this_channel, instrument_info, coord_system):
     b_min = np.nanmin(coord2)
     b_max = np.nanmax(coord2)
 
+    # find index of min b value
     b1_index = np.argmin(coord2)
     b2_index = np.argmax(coord2)
     a1 = coord1[b1_index]
@@ -108,6 +111,22 @@ def find_corners_MIRI(input, this_channel, instrument_info, coord_system):
     lambda_min = np.nanmin(lam)
     lambda_max = np.nanmax(lam)
 
+    # before returning ra should be between 0 to 360
+    if a_min < 0:
+        a_min = a_min + 360
+    if a_max > 360.0:
+        a_max = a_max - 360.0
+
+    if a1 < 0:
+        a1 = a1 + 360
+    if a1 > 360.0:
+        a1 = a1 - 360.0
+
+    if a2 < 0:
+        a2 = a2 + 360
+    if a2 > 360.0:
+        a2 = a2 - 360.0
+
     return a_min, b1, a_max, b2, a1, b_min, a2, b_max, lambda_min, lambda_max
 # *****************************************************************************
 
@@ -224,38 +243,3 @@ def find_corners_NIRSPEC(input, this_channel, instrument_info, coord_system):
     lambda_max = max(lambda_slice)
     return a_min, b1, a_max, b2, a1, b_min, a2, b_max, lambda_min, lambda_max
 # ______________________________________________________________________________
-
-
-def wrap_ra(ravalues):
-    """Test for 0/360 wrapping in ra values.
-
-    If exists it makes it difficult to determine
-    ra range of IFU cube. So put them all on "one side" of 0/360 border
-
-    Parameters
-    ----------
-    input : ravalues
-        RA values numpy.ndarray
-
-    Returns
-    ------
-    a numpy array of ra values all on "same side" of 0/360 border
-    """
-
-    index_good = np.where(np.isfinite(ravalues))
-    ravalues_wrap = ravalues[index_good].copy()
-    median_ra = np.nanmedian(ravalues_wrap)
-
-    # using median to test if there is any wrapping going on
-    wrap_index = np.where(np.fabs(ravalues_wrap - median_ra) > 180.0)
-    nwrap = wrap_index[0].size
-
-    # get all the ra on the same "side" of 0/360
-    if(nwrap != 0 and median_ra < 180):
-        ravalues_wrap[wrap_index] = ravalues_wrap[wrap_index] - 360.0
-
-    if(nwrap != 0 and median_ra > 180):
-        ravalues_wrap[wrap_index] = ravalues_wrap[wrap_index] + 360.0
-
-    return ravalues_wrap
-# ______________________________________________________________________________

jwst/cube_build/ifu_cube.py

@@ -16,6 +16,7 @@
 from astropy import units as u
 from gwcs import wcstools
 from ..assign_wcs import nirspec
+from ..assign_wcs.util import wrap_ra
 from ..datamodels import dqflags
 from . import cube_build_wcs_util
 from . import cube_overlap
@@ -1324,13 +1325,15 @@ def setup_ifucube_wcs(self):
                 lambda_max.append(lmax)
     # ________________________________________________________________________________
     # done looping over files determine final size of cube
+        corner_a = np.array(corner_a)
+        corner_a = wrap_ra(corner_a)
 
         final_a_min = min(corner_a)
         final_a_max = max(corner_a)
         final_b_min = min(corner_b)
         final_b_max = max(corner_b)
 
-        log.debug(f'final a and b:{final_a_min, final_b_min, final_a_max, final_b_max}')
+        log.info(f'final a and b:{final_a_min, final_b_min, final_a_max, final_b_max}')
         log.debug(f'final wave:   {min(lambda_min), max(lambda_max)}')
 
         # for MIRI wavelength range read in from cube pars reference file

jwst/resample/resample_spec.py

@@ -10,7 +10,7 @@
 from astropy.modeling.fitting import LinearLSQFitter
 from gwcs import wcstools, WCS
 from gwcs import coordinate_frames as cf
-from ..cube_build.cube_build_wcs_util import wrap_ra
+from ..assign_wcs.util import wrap_ra
 
 from .. import datamodels
 from . import gwcs_drizzle