From 85ecd6d4c76b591a24990acc1f2f4e262306fd69 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Fri, 14 Feb 2025 08:17:42 -0700
Subject: [PATCH 01/21] updates to MIRI LRS s_region and resample WCS

---
 jwst/assign_wcs/assign_wcs.py            |  12 +-
 jwst/assign_wcs/assign_wcs_step.py       |   1 +
 jwst/assign_wcs/miri.py                  |  68 +++++-----
 jwst/assign_wcs/util.py                  |  45 ++++++-
 jwst/regtest/test_miri_lrs_slit_spec2.py |  11 +-
 jwst/resample/resample_spec.py           | 151 +++++++++++++++++++----
 jwst/resample/resample_spec_step.py      |   4 +-
 7 files changed, 222 insertions(+), 70 deletions(-)

diff --git a/jwst/assign_wcs/assign_wcs.py b/jwst/assign_wcs/assign_wcs.py
index ca61e6d188..776d50ae45 100644
--- a/jwst/assign_wcs/assign_wcs.py
+++ b/jwst/assign_wcs/assign_wcs.py
@@ -2,7 +2,8 @@
 import importlib
 from gwcs.wcs import WCS
 from .util import (update_s_region_spectral, update_s_region_imaging,
-                   update_s_region_nrs_ifu, update_s_region_mrs)
+                   update_s_region_nrs_ifu, update_s_region_mrs,
+                   update_s_region_lrs)
 from ..lib.exposure_types import IMAGING_TYPES, SPEC_TYPES, NRS_LAMP_MODE_SPEC_TYPES
 from ..lib.dispaxis import get_dispersion_direction
 from ..lib.wcs_utils import get_wavelengths
@@ -72,8 +73,13 @@ def load_wcs(input_model, reference_files={}, nrs_slit_y_range=None):
 
         if output_model.meta.exposure.type.lower() not in exclude_types:
             imaging_types = IMAGING_TYPES.copy()
-            imaging_types.update(['mir_lrs-fixedslit', 'mir_lrs-slitless'])
-            if output_model.meta.exposure.type.lower() in imaging_types:
+            imaging_lrs_types = ['mir_lrs-fixedslit'] # keep it separate
+            imaging_types.update(['mir_lrs-slitless'])
+            if output_model.meta.exposure.type.lower() in imaging_lrs_types:
+                # uses slits corners in V2, V3 that are read in from the
+                # lrs specwcs reference file
+                update_s_region_lrs(output_model, reference_files)
+            elif output_model.meta.exposure.type.lower() in IMAGING_TYPES:
                 try:
                     update_s_region_imaging(output_model)
                 except Exception as exc:
diff --git a/jwst/assign_wcs/assign_wcs_step.py b/jwst/assign_wcs/assign_wcs_step.py
index 8d5ba034d8..9c2a6c268f 100755
--- a/jwst/assign_wcs/assign_wcs_step.py
+++ b/jwst/assign_wcs/assign_wcs_step.py
@@ -95,6 +95,7 @@ def process(self, input, *args, **kwargs):
                     message = "MSA metadata file (MSAMETFL) is required for NRS_MSASPEC exposures."
                     log.error(message)
                     raise MSAFileError(message)
+
             slit_y_range = [self.slit_y_low, self.slit_y_high]
             result = load_wcs(input_model, reference_file_names, slit_y_range)
 
diff --git a/jwst/assign_wcs/miri.py b/jwst/assign_wcs/miri.py
index fd3df3a264..adf14a8dca 100644
--- a/jwst/assign_wcs/miri.py
+++ b/jwst/assign_wcs/miri.py
@@ -5,7 +5,6 @@
 from astropy.modeling import models
 from astropy import coordinates as coord
 from astropy import units as u
-from astropy.io import fits
 
 from scipy.interpolate import UnivariateSpline
 import gwcs.coordinate_frames as cf
@@ -13,7 +12,7 @@
 
 from stdatamodels.jwst.datamodels import (DistortionModel, FilteroffsetModel,
                                           DistortionMRSModel, WavelengthrangeModel,
-                                          RegionsModel, SpecwcsModel)
+                                          RegionsModel, SpecwcsModel, MiriLRSSpecwcsModel)
 from stdatamodels.jwst.transforms.models import (MIRI_AB2Slice, IdealToV2V3)
 
 from . import pointing
@@ -239,7 +238,6 @@ def lrs_xytoabl(input_model, reference_files):
     the "specwcs" and "distortion" reference files.
 
     """
-
     # subarray to full array transform
     subarray2full = subarray_transform(input_model)
 
@@ -253,19 +251,15 @@ def lrs_xytoabl(input_model, reference_files):
     else:
         subarray_dist = distortion
 
-    ref = fits.open(reference_files['specwcs'])
+    refmodel = MiriLRSSpecwcsModel(reference_files['specwcs'])
+    if input_model.meta.exposure.type.lower() == 'mir_lrs-fixedslit':
+        zero_point = refmodel.meta.x_ref - 1, refmodel.meta.y_ref - 1
+    elif input_model.meta.exposure.type.lower() == 'mir_lrs-slitless':
+        zero_point = refmodel.meta.x_ref_slitless - 1, \
+            refmodel.meta.y_ref_slitless - 1
+        # Transform to slitless subarray from full array
+        zero_point = subarray2full.inverse(zero_point[0], zero_point[1])
 
-    with ref:
-        lrsdata = np.array([d for d in ref[1].data])
-        # Get the zero point from the reference data.
-        # The zero_point is X, Y  (which should be COLUMN, ROW)
-        # These are 1-indexed in CDP-7 (i.e., SIAF convention) so must be converted to 0-indexed
-        if input_model.meta.exposure.type.lower() == 'mir_lrs-fixedslit':
-            zero_point = ref[0].header['imx'] - 1, ref[0].header['imy'] - 1
-        elif input_model.meta.exposure.type.lower() == 'mir_lrs-slitless':
-            zero_point = ref[0].header['imxsltl'] - 1, ref[0].header['imysltl'] - 1
-            # Transform to slitless subarray from full array
-            zero_point = subarray2full.inverse(zero_point[0], zero_point[1])
 
     # Figure out the typical along-slice pixel scale at the center of the slit
     v2_cen, v3_cen = subarray_dist(zero_point[0], zero_point[1])
@@ -276,14 +270,14 @@ def lrs_xytoabl(input_model, reference_files):
     # centroid trace along the detector in pixels relative to nominal location.
     # x0,y0(ul) x1,y1 (ur) x2,y2(lr) x3,y3(ll) define corners of the box within which the distortion
     # and wavelength calibration was derived
-    xcen = lrsdata[:, 0]
-    ycen = lrsdata[:, 1]
-    wavetab = lrsdata[:, 2]
-    x0 = lrsdata[:, 3]
-    y0 = lrsdata[:, 4]
-    x1 = lrsdata[:, 5]
-    y2 = lrsdata[:, 8]
-
+    xcen = refmodel.wavetable.x_center
+    ycen = refmodel.wavetable.y_center
+    wavetab = refmodel.wavetable.wavelength
+    x0 = refmodel.wavetable.x0
+    y0 = refmodel.wavetable.y0
+    x1 = refmodel.wavetable.x1
+    y2 = refmodel.wavetable.y2
+    refmodel.close()
     # If in fixed slit mode, define the bounding box using the corner locations provided in
     # the CDP reference file.
     if input_model.meta.exposure.type.lower() == 'mir_lrs-fixedslit':
@@ -313,9 +307,10 @@ def lrs_xytoabl(input_model, reference_files):
     # This function will give slit dX as a function of Y subarray pixel value
     dxmodel = models.Tabular1D(lookup_table=xshiftref, points=ycen_subarray, name='xshiftref',
                                  bounds_error=False, fill_value=np.nan)
+    # remove commented out how after testing
+    #if input_model.meta.exposure.type.lower() == 'mir_lrs-fixedslit':
+    #    bb_sub = (bb_sub[0], (dxmodel.points[0].min(), dxmodel.points[0].max()))
 
-    if input_model.meta.exposure.type.lower() == 'mir_lrs-fixedslit':
-        bb_sub = (bb_sub[0], (dxmodel.points[0].min(), dxmodel.points[0].max()))
     # Fit for the wavelength as a function of Y
     # Reverse the vectors so that yinv is increasing (needed for spline fitting function)
     # Spline fit with enforced smoothness
@@ -405,19 +400,16 @@ def lrs_abltov2v3l(input_model, reference_files):
     else:
         subarray_dist = distortion
 
-    ref = fits.open(reference_files['specwcs'])
-
-    with ref:
-        # Get the zero point from the reference data.
-        # The zero_point is X, Y  (which should be COLUMN, ROW)
-        # These are 1-indexed in CDP-7 (i.e., SIAF convention) so must be converted to 0-indexed
-        if input_model.meta.exposure.type.lower() == 'mir_lrs-fixedslit':
-            zero_point = ref[0].header['imx'] - 1, ref[0].header['imy'] - 1
-        elif input_model.meta.exposure.type.lower() == 'mir_lrs-slitless':
-            zero_point = ref[0].header['imxsltl'] - 1, ref[0].header['imysltl'] - 1
-            # Transform to slitless subarray from full array
-            zero_point = subarray2full.inverse(zero_point[0], zero_point[1])
-
+    refmodel = MiriLRSSpecwcsModel(reference_files['specwcs'])
+    if input_model.meta.exposure.type.lower() == 'mir_lrs-fixedslit':
+        zero_point = refmodel.meta.x_ref - 1, refmodel.meta.y_ref - 1
+    elif input_model.meta.exposure.type.lower() == 'mir_lrs-slitless':
+        zero_point = refmodel.meta.x_ref_slitless - 1, \
+            refmodel.meta.y_ref_slitless - 1
+        # Transform to slitless subarray from full array
+        zero_point = subarray2full.inverse(zero_point[0], zero_point[1])
+
+    refmodel.close()
     # Figure out the typical along-slice pixel scale at the center of the slit
     v2_cen, v3_cen = subarray_dist(zero_point[0], zero_point[1])
     v2_off, v3_off = subarray_dist(zero_point[0] + 1, zero_point[1])
diff --git a/jwst/assign_wcs/util.py b/jwst/assign_wcs/util.py
index 695f966502..f3714455e6 100644
--- a/jwst/assign_wcs/util.py
+++ b/jwst/assign_wcs/util.py
@@ -19,10 +19,10 @@
 from stpipe.exceptions import StpipeExitException
 from stcal.alignment.util import compute_s_region_keyword, compute_s_region_imaging
 
-from stdatamodels.jwst.datamodels import WavelengthrangeModel
+from stdatamodels.jwst.datamodels import WavelengthrangeModel, MiriLRSSpecwcsModel
 from stdatamodels.jwst.transforms.models import GrismObject
 
-from ..lib.catalog_utils import SkyObject
+from jwst.lib.catalog_utils import SkyObject
 
 
 log = logging.getLogger(__name__)
@@ -807,6 +807,45 @@ def update_s_region_imaging(model):
         model.meta.wcsinfo.s_region = s_region
 
 
+def update_s_region_lrs(model, reference_files):
+    """
+    Determine ``S_REGION`` using V2,V3 of the slit corners from reference file.
+    """
+
+    refmodel = MiriLRSSpecwcsModel(reference_files['specwcs'])
+
+    v2vert1 = refmodel.meta.v2_vert1
+    v2vert2 = refmodel.meta.v2_vert2
+    v2vert3 = refmodel.meta.v2_vert3
+    v2vert4 = refmodel.meta.v2_vert4
+
+    v3vert1 = refmodel.meta.v3_vert1
+    v3vert2 = refmodel.meta.v3_vert2
+    v3vert3 = refmodel.meta.v3_vert3
+    v3vert4 = refmodel.meta.v3_vert4
+
+    refmodel.close()
+    v2 = [v2vert1, v2vert2, v2vert3, v2vert4]
+    v3 = [v3vert1, v3vert2, v3vert3, v3vert4]
+
+    if (any(elem is None for elem in v2) or
+        any(elem is None for elem in v3)):
+        log.info("The V2,V3 coordinates of the MIRI LRS-Fixed slit contains NaN values.")
+        log.info("The s_region will not be updated")
+
+    lam = 7.0 # wavelength does not matter for s region assign a value in
+              # wavelength of MIRI LRS
+    s = model.meta.wcs.transform('v2v3', 'world', v2, v3, lam)
+    a = s[0]
+    b = s[1]
+    footprint = np.array([[a[0], b[0]],
+                          [a[1], b[1]],
+                          [a[2], b[2]],
+                          [a[3], b[3]]])
+
+    update_s_region_keyword(model, footprint)
+    print('in  assign_wcs util.py',model.meta.wcsinfo.s_region)
+
 def compute_footprint_spectral(model):
     """
     Determine spatial footprint for spectral observations using the instrument model.
@@ -844,9 +883,11 @@ def compute_footprint_spectral(model):
     return footprint, (lam_min, lam_max)
 
 
+
 def update_s_region_spectral(model):
     """ Update the S_REGION keyword.
     """
+
     footprint, spectral_region = compute_footprint_spectral(model)
     update_s_region_keyword(model, footprint)
     model.meta.wcsinfo.spectral_region = spectral_region
diff --git a/jwst/regtest/test_miri_lrs_slit_spec2.py b/jwst/regtest/test_miri_lrs_slit_spec2.py
index 65765fd86a..4d819fe46e 100644
--- a/jwst/regtest/test_miri_lrs_slit_spec2.py
+++ b/jwst/regtest/test_miri_lrs_slit_spec2.py
@@ -7,7 +7,7 @@
 
 from jwst.stpipe import Step
 from jwst.extract_1d import Extract1dStep
-
+from stcal.alignment import util
 
 @pytest.fixture(scope="module")
 def run_pipeline(rtdata_module):
@@ -87,3 +87,12 @@ def test_miri_lrs_slit_wcs(run_pipeline, rtdata_module, fitsdiff_default_kwargs)
         xtruth, ytruth = im_truth.meta.wcs.backward_transform(ratruth, dectruth, lamtruth)
         assert_allclose(xtest, xtruth)
         assert_allclose(ytest, ytruth)
+
+        # Test the s_region. S_region is formed by footprint which contains
+        # floats rather than a string. Test footprint
+        sregion = im.meta.wcsinfo.s_region
+        sregion_test = im_truth.meta.wcsinfo.s_region
+        footprint=util._sregion_to_footprint(sregion)
+        footprint_test = util._sregion_to_footprint(sregion_test)
+        assert_allclose(footprint, footprint_test)
+
diff --git a/jwst/resample/resample_spec.py b/jwst/resample/resample_spec.py
index 4f991d76ce..e0a216b7ee 100644
--- a/jwst/resample/resample_spec.py
+++ b/jwst/resample/resample_spec.py
@@ -10,17 +10,20 @@
 )
 from astropy.modeling.fitting import LinearLSQFitter
 from astropy.stats import sigma_clip
+from astropy.coordinates import SkyCoord
+
 from astropy.utils.exceptions import AstropyUserWarning
 from gwcs import wcstools, WCS
 from gwcs import coordinate_frames as cf
 
 from stdatamodels.jwst import datamodels
 
-from jwst.assign_wcs.util import compute_scale, wcs_bbox_from_shape, wrap_ra
+from jwst.assign_wcs.util import compute_scale, wcs_bbox_from_shape,\
+    wrap_ra
 from jwst.resample import resample_utils
 from jwst.resample.resample import ResampleImage
 from jwst.datamodels import ModelLibrary
-
+from stcal.alignment.util import compute_s_region_keyword
 
 log = logging.getLogger(__name__)
 log.setLevel(logging.DEBUG)
@@ -66,6 +69,7 @@ def __init__(self, input_models, pixfrac=1.0, kernel="square",
         pixel_scale = None
         pixel_area = None
         pixel_scale_ratio = 1.0
+        s_region = None
 
         if isinstance(output_wcs, dict):
             output_wcs_dict = {
@@ -163,7 +167,7 @@ def __init__(self, input_models, pixfrac=1.0, kernel="square",
             # Any other customizations (crpix, crval, rotation) are ignored.
             if resample_utils.is_sky_like(input_models[0].meta.wcs.output_frame):
                 if input_models[0].meta.instrument.name != "NIRSPEC":
-                    output_wcs = self.build_interpolated_output_wcs(
+                    output_wcs, s_region = self.build_interpolated_output_wcs(
                         input_models,
                         pixel_scale_ratio=pixel_scale_ratio
                     )
@@ -189,6 +193,7 @@ def __init__(self, input_models, pixfrac=1.0, kernel="square",
                 output_pix_area = None
 
         self._spec_output_pix_area = output_pix_area
+        self._s_region = s_region
 
         if pixel_scale is None:
             log.info(f'Specified output pixel scale ratio: {pixel_scale_ratio}.')
@@ -203,6 +208,7 @@ def __init__(self, input_models, pixfrac=1.0, kernel="square",
         output_wcs_dict["pixel_scale"] = pixel_scale
         output_wcs_dict["pixel_scale_ratio"] = pixel_scale_ratio
 
+
         library = ModelLibrary(input_models, on_disk=False)
 
         super().__init__(
@@ -242,7 +248,9 @@ def update_output_model(self, model, info_dict):
             model.meta.photometry.pixelarea_arcsecsq = (
                 self._spec_output_pix_area * np.rad2deg(3600)**2
             )
-
+        if self._s_region is not None:
+            model.meta.wcsinfo.s_region = self._s_region
+            log.info(f'Updating S_REGION: {self._s_region}.')
         # TODO: this is helpful info that should be stored in products.
         #       Not storing this at this time in order to reduce the number of
         #       failures in the regression tests.
@@ -576,13 +584,18 @@ def build_interpolated_output_wcs(self, input_models, pixel_scale_ratio=1.0):
         all_dec_slit = []
         xstop = 0
 
+        sregion_list = []
+        all_wcs = [m.meta.wcs for m in input_models]
         for im, model in enumerate(input_models):
             wcs = model.meta.wcs
             bbox = wcs.bounding_box
             grid = wcstools.grid_from_bounding_box(bbox)
             ra, dec, lam = np.array(wcs(*grid))
-            # Handle vertical (MIRI) or horizontal (NIRSpec) dispersion.  The
-            # following 2 variables are 0 or 1, i.e. zero-indexed in x,y WCS order
+
+            sregion_list.append(model.meta.wcsinfo.s_region)
+
+            # Handle vertical (MIRI).  The following 2 variables are
+            # 0 or 1, i.e. zero-indexed in x,y WCS order
             spectral_axis = find_dispersion_axis(model)
             spatial_axis = spectral_axis ^ 1
 
@@ -599,7 +612,7 @@ def build_interpolated_output_wcs(self, input_models, pixel_scale_ratio=1.0):
             # sampling.
 
             # Steps to do this for first input model:
-            # 1. find the middle of the spectrum in wavelength
+            # 1. Find the middle of the spectrum in wavelength
             # 2. Pull out the ra and dec at the center of the slit.
             # 3. Find the mean ra,dec and the center of the slit this will
             #    represent the tangent point
@@ -614,7 +627,7 @@ def build_interpolated_output_wcs(self, input_models, pixel_scale_ratio=1.0):
                 lam_center_index = int((bbox[spectral_axis][1] -
                                         bbox[spectral_axis][0]) / 2)
                 if spatial_axis == 0:
-                    # MIRI LRS, the WCS x axis is spatial
+                    # MIRI LRS spectral = 1, the spatial axis = 0
                     ra_slice = ra[lam_center_index, :]
                     dec_slice = dec[lam_center_index, :]
                 else:
@@ -637,10 +650,10 @@ def build_interpolated_output_wcs(self, input_models, pixel_scale_ratio=1.0):
                 x_tan, y_tan = undist2sky1.inverse(ra, dec)
 
                 # pull out data from center
-                if spectral_axis == 0:  # MIRI LRS, the WCS x axis is spatial
+                if spectral_axis == 0:
                     x_tan_array = x_tan.T[lam_center_index]
                     y_tan_array = y_tan.T[lam_center_index]
-                else:
+                else: # MIRI LRS Spectral Axis = 1, the WCS x axis is spatial
                     x_tan_array = x_tan[lam_center_index]
                     y_tan_array = y_tan[lam_center_index]
 
@@ -750,26 +763,31 @@ def build_interpolated_output_wcs(self, input_models, pixel_scale_ratio=1.0):
 
         native2celestial = RotateNative2Celestial(ra_center_final, dec_center_final, 180)
         undist2sky = tan | native2celestial
-        # find the spatial size of the output - same in x,y
-        if swap_xy:
-            _, x_tan_all = undist2sky.inverse(all_ra, all_dec)
-            pix_to_tan_slope = pix_to_ytan.slope
-        else:
-            x_tan_all, _ = undist2sky.inverse(all_ra, all_dec)
-            pix_to_tan_slope = pix_to_xtan.slope
 
-        x_min = np.amin(x_tan_all)
-        x_max = np.amax(x_tan_all)
-        x_size = int(np.ceil((x_max - x_min) / np.absolute(pix_to_tan_slope)))
+        ## Use all the wcs
+        min_tan_x, max_tan_x, min_tan_y, max_tan_y = self._max_spatial_extent(
+            all_wcs, undist2sky.inverse)
+        diff_y = np.abs(max_tan_y - min_tan_y)
+        diff_x = np.abs(max_tan_x - min_tan_x)
+        pix_to_tan_slope_y = np.abs(pix_to_ytan.slope)
+        slope_sign_y = np.sign(pix_to_ytan.slope)
+        pix_to_tan_slope_x = np.abs(pix_to_xtan.slope)
+        slope_sign_x = np.sign(pix_to_xtan.slope)
+
         if swap_xy:
-            pix_to_ytan.intercept = -0.5 * (x_size - 1) * pix_to_ytan.slope
+            ny = int(np.ceil(diff_y / pix_to_tan_slope_y)) + 1
         else:
-            pix_to_xtan.intercept = -0.5 * (x_size - 1) * pix_to_xtan.slope
+            ny = int(np.ceil(diff_x / pix_to_tan_slope_x)) + 1
+
+        offset_y  = (ny)/2 * pix_to_tan_slope_y
+        offset_x  = (ny)/2 * pix_to_tan_slope_x
+        pix_to_ytan.intercept =  - slope_sign_y * offset_y
+        pix_to_xtan.intercept =  - slope_sign_x * offset_x
 
         # single model: use size of x_tan_array
         # to be consistent with method before
         if len(input_models) == 1:
-            x_size = int(np.ceil(xstop))
+            ny = int(np.ceil(xstop))
 
         # define the output wcs
         transform = mapping | (pix_to_xtan & pix_to_ytan | undist2sky) & pix_to_wavelength
@@ -789,7 +807,7 @@ def build_interpolated_output_wcs(self, input_models, pixel_scale_ratio=1.0):
         # compute the output array size in WCS axes order, i.e. (x, y)
         output_array_size = [0, 0]
         output_array_size[spectral_axis] = int(np.ceil(len(wavelength_array)))
-        output_array_size[spatial_axis] = x_size
+        output_array_size[spatial_axis] = ny
 
         # turn the size into a numpy shape in (y, x) order
         output_wcs.array_shape = output_array_size[::-1]
@@ -797,7 +815,8 @@ def build_interpolated_output_wcs(self, input_models, pixel_scale_ratio=1.0):
         bounding_box = wcs_bbox_from_shape(output_array_size[::-1])
         output_wcs.bounding_box = bounding_box
 
-        return output_wcs
+        s_region = find_miri_lrs_sregion(sregion_list,output_wcs)
+        return output_wcs, s_region
 
     def build_nirspec_lamp_output_wcs(self, input_models, pixel_scale_ratio):
         """
@@ -1011,3 +1030,85 @@ def compute_spectral_pixel_scale(wcs, fiducial=None, disp_axis=1):
 
     pixel_scale = compute_scale(wcs, fiducial, disp_axis=disp_axis)
     return float(pixel_scale)
+
+def find_miri_lrs_sregion(sregion_list,wcs):
+    """ Find s region for MIRI LRS resampled data.
+
+    Parameters
+    ----------
+    sregion_list : list
+        List of s_regions.
+    wcs : gwcs.WCS
+        Spatial/spectral WCS.
+
+    Returns
+    -------
+    sregion : string
+        S_region for the resample data.
+    """
+    # use the first sregion to set the width of the slit
+
+    spatial_box = sregion_list[0]
+    s = spatial_box.split(' ')
+    a1 = float(s[3])
+    b1 = float(s[4])
+    a2 = float(s[5])
+    b2 = float(s[6])
+    a3 = float(s[7])
+    b3 = float(s[8])
+    a4 = float(s[9])
+    b4 = float(s[10])
+
+    # convert each corner to SkyCoord
+    coord1 = SkyCoord(a1, b1, unit='deg')
+    coord2 = SkyCoord(a2, b2, unit='deg')
+    coord3 = SkyCoord(a3, b3, unit='deg')
+    coord4 = SkyCoord(a4, b4, unit='deg')
+
+    # Find the distance between the corners
+    # corners are counter clockwize from 1,2,3,4
+    sep1 = coord1.separation(coord2)
+    sep2 = coord2.separation(coord3)
+    sep3 = coord3.separation(coord4)
+    sep4 = coord4.separation(coord1)
+
+    # use the separation values so we can find the min value later
+    sep = [sep1.value, sep2.value, sep3.value, sep4.value]
+
+    # the minimum separation is the slit width
+    min_sep = np.min(sep)
+    min_sep = min_sep* u.deg # set the units to degrees
+
+    log.info(f'Estimated MIRI LRS slit width: {min_sep*3600} arcsec.')
+    # now use the combined WCS to map all pixels to the slit center
+    bbox = wcs.bounding_box
+    grid = wcstools.grid_from_bounding_box(bbox)
+    ra, dec, _ = np.array(wcs(*grid))
+    ra = ra.flatten()
+    dec = dec.flatten()
+    # ra and dec are the values along the output resampled slit center
+    # using the first point and last point find the position angle
+    star1 = SkyCoord(ra[0]*u.deg, dec[0]*u.deg, frame='icrs')
+    star2 = SkyCoord(ra[-1]*u.deg, dec[-1]*u.deg, frame='icrs')
+    position_angle = star1.position_angle(star2).to(u.deg)
+
+    # 90 degrees to the position angle of the slit will define s_region
+    pos_angle = position_angle - 90.0*u.deg
+
+    star_c1 = star1.directional_offset_by(pos_angle, min_sep/2)
+    star_c2 = star1.directional_offset_by(pos_angle, -min_sep/2)
+    star_c3 = star2.directional_offset_by(pos_angle, min_sep/2)
+    star_c4 = star2.directional_offset_by(pos_angle, -min_sep/2)
+
+    # set  these values to footprint
+    # ra,dec corners are in counter-clockwise direction
+    footprint = [star_c1.ra.value, star_c1.dec.value,
+                 star_c3.ra.value, star_c3.dec.value,
+                 star_c4.ra.value, star_c4.dec.value,
+                 star_c2.ra.value, star_c2.dec.value]
+    footprint = np.array(footprint)
+    s_region = compute_s_region_keyword(footprint)
+    return s_region
+
+
+
diff --git a/jwst/resample/resample_spec_step.py b/jwst/resample/resample_spec_step.py
index 67a5f6e91a..7d5cc03e46 100755
--- a/jwst/resample/resample_spec_step.py
+++ b/jwst/resample/resample_spec_step.py
@@ -265,8 +265,10 @@ def _process_slit(self, input_models):
         else:
             result.meta.resample.pixel_scale_ratio = resamp.pixel_scale_ratio
         result.meta.resample.pixfrac = self.pixfrac
+
         self.update_slit_metadata(result)
-        update_s_region_spectral(result)
+        if result.meta.exposure.type.lower() != 'mir_lrs-fixedslit':
+            update_s_region_spectral(result)
 
         return result
 

From ed04b54f80f3806461bb4c2327b613ce03350a5e Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Fri, 14 Feb 2025 08:41:58 -0700
Subject: [PATCH 02/21] changed pyproject.toml to point to stdatamodels with
 new MIRI LRS specwcs model

---
 pyproject.toml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 87a88914b1..90f6616258 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -34,7 +34,8 @@ dependencies = [
     "scikit-image>=0.20.0",
     "scipy>=1.14.1",
     "spherical-geometry>=1.2.22",
-    "stdatamodels @ git+https://github.com/spacetelescope/stdatamodels.git@main",
+    #"stdatamodels @ git+https://github.com/spacetelescope/stdatamodels.git@main",
+    "stdatamodels @ stdatamodels-2.2.1.dev20+gf85885d",
     "stcal @ git+https://github.com/spacetelescope/stcal.git@main",
     "stpipe>=0.8.0,<0.9.0",
     "stsci.imagestats>=1.6.3",

From a5626eb21e71966b942f98c207de71840d2f215c Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Fri, 14 Feb 2025 08:53:22 -0700
Subject: [PATCH 03/21] added change log

---
 changes/9193.assign_wcs.rst | 1 +
 changes/9193.resample.rst   | 1 +
 2 files changed, 2 insertions(+)
 create mode 100644 changes/9193.assign_wcs.rst
 create mode 100644 changes/9193.resample.rst

diff --git a/changes/9193.assign_wcs.rst b/changes/9193.assign_wcs.rst
new file mode 100644
index 0000000000..5090d580f7
--- /dev/null
+++ b/changes/9193.assign_wcs.rst
@@ -0,0 +1 @@
+Fix MIRI LRS s_region in assign_wcs
diff --git a/changes/9193.resample.rst b/changes/9193.resample.rst
new file mode 100644
index 0000000000..f4e305c174
--- /dev/null
+++ b/changes/9193.resample.rst
@@ -0,0 +1 @@
+Fix MIRI LRS s_region and WCS  in resample_spec

From 3b79df2c6f69e1c193a1fa43776ce070181b3f14 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Fri, 14 Feb 2025 11:09:41 -0700
Subject: [PATCH 04/21] update location of stdatamodels to use

---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 90f6616258..f513e9534b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -35,7 +35,7 @@ dependencies = [
     "scipy>=1.14.1",
     "spherical-geometry>=1.2.22",
     #"stdatamodels @ git+https://github.com/spacetelescope/stdatamodels.git@main",
-    "stdatamodels @ stdatamodels-2.2.1.dev20+gf85885d",
+    "stdatamodels @ git+https://github.com/jemorrison/stdatamodels.git@MIRI_LRS_Specwcs",
     "stcal @ git+https://github.com/spacetelescope/stcal.git@main",
     "stpipe>=0.8.0,<0.9.0",
     "stsci.imagestats>=1.6.3",

From fbe54c13d01ffd2038e968241933fb13d3541e81 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Mon, 17 Feb 2025 08:14:05 -0700
Subject: [PATCH 05/21] fix ny of s2d region

---
 jwst/resample/resample_spec.py | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/jwst/resample/resample_spec.py b/jwst/resample/resample_spec.py
index e0a216b7ee..e726e41ae5 100644
--- a/jwst/resample/resample_spec.py
+++ b/jwst/resample/resample_spec.py
@@ -784,11 +784,6 @@ def build_interpolated_output_wcs(self, input_models, pixel_scale_ratio=1.0):
         pix_to_ytan.intercept =  - slope_sign_y * offset_y
         pix_to_xtan.intercept =  - slope_sign_x * offset_x
 
-        # single model: use size of x_tan_array
-        # to be consistent with method before
-        if len(input_models) == 1:
-            ny = int(np.ceil(xstop))
-
         # define the output wcs
         transform = mapping | (pix_to_xtan & pix_to_ytan | undist2sky) & pix_to_wavelength
 

From d6bcb7841c9351e955f4af6b75e36ce15f7d66f4 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Mon, 17 Feb 2025 12:28:17 -0700
Subject: [PATCH 06/21] update from review

---
 changes/9193.resample.rst     | 2 +-
 jwst/assign_wcs/assign_wcs.py | 2 +-
 jwst/assign_wcs/miri.py       | 3 ---
 jwst/assign_wcs/util.py       | 4 +---
 4 files changed, 3 insertions(+), 8 deletions(-)

diff --git a/changes/9193.resample.rst b/changes/9193.resample.rst
index f4e305c174..0836cebda5 100644
--- a/changes/9193.resample.rst
+++ b/changes/9193.resample.rst
@@ -1 +1 @@
-Fix MIRI LRS s_region and WCS  in resample_spec
+Fix MIRI LRS s_region and WCS in resample_spec
diff --git a/jwst/assign_wcs/assign_wcs.py b/jwst/assign_wcs/assign_wcs.py
index 776d50ae45..db9eeee00a 100644
--- a/jwst/assign_wcs/assign_wcs.py
+++ b/jwst/assign_wcs/assign_wcs.py
@@ -73,7 +73,7 @@ def load_wcs(input_model, reference_files={}, nrs_slit_y_range=None):
 
         if output_model.meta.exposure.type.lower() not in exclude_types:
             imaging_types = IMAGING_TYPES.copy()
-            imaging_lrs_types = ['mir_lrs-fixedslit'] # keep it separate
+            imaging_lrs_types = ['mir_lrs-fixedslit']
             imaging_types.update(['mir_lrs-slitless'])
             if output_model.meta.exposure.type.lower() in imaging_lrs_types:
                 # uses slits corners in V2, V3 that are read in from the
diff --git a/jwst/assign_wcs/miri.py b/jwst/assign_wcs/miri.py
index adf14a8dca..0b151f80ba 100644
--- a/jwst/assign_wcs/miri.py
+++ b/jwst/assign_wcs/miri.py
@@ -307,9 +307,6 @@ def lrs_xytoabl(input_model, reference_files):
     # This function will give slit dX as a function of Y subarray pixel value
     dxmodel = models.Tabular1D(lookup_table=xshiftref, points=ycen_subarray, name='xshiftref',
                                  bounds_error=False, fill_value=np.nan)
-    # remove commented out how after testing
-    #if input_model.meta.exposure.type.lower() == 'mir_lrs-fixedslit':
-    #    bb_sub = (bb_sub[0], (dxmodel.points[0].min(), dxmodel.points[0].max()))
 
     # Fit for the wavelength as a function of Y
     # Reverse the vectors so that yinv is increasing (needed for spline fitting function)
diff --git a/jwst/assign_wcs/util.py b/jwst/assign_wcs/util.py
index f3714455e6..824ecbfd27 100644
--- a/jwst/assign_wcs/util.py
+++ b/jwst/assign_wcs/util.py
@@ -833,8 +833,7 @@ def update_s_region_lrs(model, reference_files):
         log.info("The V2,V3 coordinates of the MIRI LRS-Fixed slit contains NaN values.")
         log.info("The s_region will not be updated")
 
-    lam = 7.0 # wavelength does not matter for s region assign a value in
-              # wavelength of MIRI LRS
+    lam = 7.0 # wavelength does not matter for s region assign a value in range of LRS
     s = model.meta.wcs.transform('v2v3', 'world', v2, v3, lam)
     a = s[0]
     b = s[1]
@@ -844,7 +843,6 @@ def update_s_region_lrs(model, reference_files):
                           [a[3], b[3]]])
 
     update_s_region_keyword(model, footprint)
-    print('in  assign_wcs util.py',model.meta.wcsinfo.s_region)
 
 def compute_footprint_spectral(model):
     """

From 9dc38cae0cd63729eb4bf9f775d1079490933340 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Tue, 18 Feb 2025 11:03:23 -0700
Subject: [PATCH 07/21] update to where resampled s_region determined

---
 jwst/assign_wcs/util.py             |  13 ++++
 jwst/resample/resample_spec.py      | 101 ++--------------------------
 jwst/resample/resample_spec_step.py |  10 ++-
 jwst/resample/resample_utils.py     |  87 ++++++++++++++++++++++++
 4 files changed, 112 insertions(+), 99 deletions(-)

diff --git a/jwst/assign_wcs/util.py b/jwst/assign_wcs/util.py
index 824ecbfd27..263b3987c8 100644
--- a/jwst/assign_wcs/util.py
+++ b/jwst/assign_wcs/util.py
@@ -810,6 +810,19 @@ def update_s_region_imaging(model):
 def update_s_region_lrs(model, reference_files):
     """
     Determine ``S_REGION`` using V2,V3 of the slit corners from reference file.
+
+    Parameters
+    ----------
+    model : DataModel
+        Input model
+    reference_files : list
+        List of reference files for assign_wcs.
+
+    Returns
+    -------
+    None
+        s_region for model is updated in place .
+
     """
 
     refmodel = MiriLRSSpecwcsModel(reference_files['specwcs'])
diff --git a/jwst/resample/resample_spec.py b/jwst/resample/resample_spec.py
index e726e41ae5..9b6cb08ef5 100644
--- a/jwst/resample/resample_spec.py
+++ b/jwst/resample/resample_spec.py
@@ -10,7 +10,6 @@
 )
 from astropy.modeling.fitting import LinearLSQFitter
 from astropy.stats import sigma_clip
-from astropy.coordinates import SkyCoord
 
 from astropy.utils.exceptions import AstropyUserWarning
 from gwcs import wcstools, WCS
@@ -23,7 +22,6 @@
 from jwst.resample import resample_utils
 from jwst.resample.resample import ResampleImage
 from jwst.datamodels import ModelLibrary
-from stcal.alignment.util import compute_s_region_keyword
 
 log = logging.getLogger(__name__)
 log.setLevel(logging.DEBUG)
@@ -69,7 +67,6 @@ def __init__(self, input_models, pixfrac=1.0, kernel="square",
         pixel_scale = None
         pixel_area = None
         pixel_scale_ratio = 1.0
-        s_region = None
 
         if isinstance(output_wcs, dict):
             output_wcs_dict = {
@@ -167,7 +164,8 @@ def __init__(self, input_models, pixfrac=1.0, kernel="square",
             # Any other customizations (crpix, crval, rotation) are ignored.
             if resample_utils.is_sky_like(input_models[0].meta.wcs.output_frame):
                 if input_models[0].meta.instrument.name != "NIRSPEC":
-                    output_wcs, s_region = self.build_interpolated_output_wcs(
+
+                    output_wcs = self.build_interpolated_output_wcs(
                         input_models,
                         pixel_scale_ratio=pixel_scale_ratio
                     )
@@ -193,7 +191,6 @@ def __init__(self, input_models, pixfrac=1.0, kernel="square",
                 output_pix_area = None
 
         self._spec_output_pix_area = output_pix_area
-        self._s_region = s_region
 
         if pixel_scale is None:
             log.info(f'Specified output pixel scale ratio: {pixel_scale_ratio}.')
@@ -208,7 +205,6 @@ def __init__(self, input_models, pixfrac=1.0, kernel="square",
         output_wcs_dict["pixel_scale"] = pixel_scale
         output_wcs_dict["pixel_scale_ratio"] = pixel_scale_ratio
 
-
         library = ModelLibrary(input_models, on_disk=False)
 
         super().__init__(
@@ -248,9 +244,7 @@ def update_output_model(self, model, info_dict):
             model.meta.photometry.pixelarea_arcsecsq = (
                 self._spec_output_pix_area * np.rad2deg(3600)**2
             )
-        if self._s_region is not None:
-            model.meta.wcsinfo.s_region = self._s_region
-            log.info(f'Updating S_REGION: {self._s_region}.')
+
         # TODO: this is helpful info that should be stored in products.
         #       Not storing this at this time in order to reduce the number of
         #       failures in the regression tests.
@@ -584,7 +578,7 @@ def build_interpolated_output_wcs(self, input_models, pixel_scale_ratio=1.0):
         all_dec_slit = []
         xstop = 0
 
-        sregion_list = []
+        #sregion_list = []
         all_wcs = [m.meta.wcs for m in input_models]
         for im, model in enumerate(input_models):
             wcs = model.meta.wcs
@@ -592,8 +586,6 @@ def build_interpolated_output_wcs(self, input_models, pixel_scale_ratio=1.0):
             grid = wcstools.grid_from_bounding_box(bbox)
             ra, dec, lam = np.array(wcs(*grid))
 
-            sregion_list.append(model.meta.wcsinfo.s_region)
-
             # Handle vertical (MIRI).  The following 2 variables are
             # 0 or 1, i.e. zero-indexed in x,y WCS order
             spectral_axis = find_dispersion_axis(model)
@@ -809,9 +801,7 @@ def build_interpolated_output_wcs(self, input_models, pixel_scale_ratio=1.0):
         output_wcs.pixel_shape = output_array_size
         bounding_box = wcs_bbox_from_shape(output_array_size[::-1])
         output_wcs.bounding_box = bounding_box
-
-        s_region = find_miri_lrs_sregion(sregion_list,output_wcs)
-        return output_wcs, s_region
+        return output_wcs
 
     def build_nirspec_lamp_output_wcs(self, input_models, pixel_scale_ratio):
         """
@@ -1026,84 +1016,3 @@ def compute_spectral_pixel_scale(wcs, fiducial=None, disp_axis=1):
     pixel_scale = compute_scale(wcs, fiducial, disp_axis=disp_axis)
     return float(pixel_scale)
 
-def find_miri_lrs_sregion(sregion_list,wcs):
-    """ Find s region for MIRI LRS resampled data.
-
-    Parameters
-    ----------
-    sregion_list : list
-        List of s_regions.
-    wcs : gwcs.WCS
-        Spatial/spectral WCS.
-
-    Returns
-    -------
-    sregion : string
-        S_region for the resample data.
-    """
-    # use the first sregion to set the width of the slit
-
-    spatial_box = sregion_list[0]
-    s = spatial_box.split(' ')
-    a1 = float(s[3])
-    b1 = float(s[4])
-    a2 = float(s[5])
-    b2 = float(s[6])
-    a3 = float(s[7])
-    b3 = float(s[8])
-    a4 = float(s[9])
-    b4 = float(s[10])
-
-    # convert each corner to SkyCoord
-    coord1 = SkyCoord(a1, b1, unit='deg')
-    coord2 = SkyCoord(a2, b2, unit='deg')
-    coord3 = SkyCoord(a3, b3, unit='deg')
-    coord4 = SkyCoord(a4, b4, unit='deg')
-
-    # Find the distance between the corners
-    # corners are counter clockwize from 1,2,3,4
-    sep1 = coord1.separation(coord2)
-    sep2 = coord2.separation(coord3)
-    sep3 = coord3.separation(coord4)
-    sep4 = coord4.separation(coord1)
-
-    # use the separation values so we can find the min value later
-    sep = [sep1.value, sep2.value, sep3.value, sep4.value]
-
-    # the minimum separation is the slit width
-    min_sep = np.min(sep)
-    min_sep = min_sep* u.deg # set the units to degrees
-
-    log.info(f'Estimated MIRI LRS slit width: {min_sep*3600} arcsec.')
-    # now use the combined WCS to map all pixels to the slit center
-    bbox = wcs.bounding_box
-    grid = wcstools.grid_from_bounding_box(bbox)
-    ra, dec, _ = np.array(wcs(*grid))
-    ra = ra.flatten()
-    dec = dec.flatten()
-    # ra and dec are the values along the output resampled slit center
-    # using the first point and last point find the position angle
-    star1 = SkyCoord(ra[0]*u.deg, dec[0]*u.deg, frame='icrs')
-    star2 = SkyCoord(ra[-1]*u.deg, dec[-1]*u.deg, frame='icrs')
-    position_angle = star1.position_angle(star2).to(u.deg)
-
-    # 90 degrees to the position angle of the slit will define s_region
-    pos_angle = position_angle - 90.0*u.deg
-
-    star_c1 = star1.directional_offset_by(pos_angle, min_sep/2)
-    star_c2 = star1.directional_offset_by(pos_angle, -min_sep/2)
-    star_c3 = star2.directional_offset_by(pos_angle, min_sep/2)
-    star_c4 = star2.directional_offset_by(pos_angle, -min_sep/2)
-
-    # set  these values to footprint
-    # ra,dec corners are in counter-clockwise direction
-    footprint = [star_c1.ra.value, star_c1.dec.value,
-                 star_c3.ra.value, star_c3.dec.value,
-                 star_c4.ra.value, star_c4.dec.value,
-                 star_c2.ra.value, star_c2.dec.value]
-    footprint = np.array(footprint)
-    s_region = compute_s_region_keyword(footprint)
-    return s_region
-
-
-
diff --git a/jwst/resample/resample_spec_step.py b/jwst/resample/resample_spec_step.py
index 7d5cc03e46..12a6d1cd1b 100755
--- a/jwst/resample/resample_spec_step.py
+++ b/jwst/resample/resample_spec_step.py
@@ -6,7 +6,7 @@
 from jwst.datamodels import ModelContainer, ModelLibrary
 from jwst.lib.pipe_utils import match_nans_and_flags
 from jwst.lib.wcs_utils import get_wavelengths
-from jwst.resample.resample_utils import load_custom_wcs
+from jwst.resample.resample_utils import load_custom_wcs, find_miri_lrs_sregion
 
 from . import resample_spec, ResampleStep
 from ..exp_to_source import multislit_to_container
@@ -267,9 +267,13 @@ def _process_slit(self, input_models):
         result.meta.resample.pixfrac = self.pixfrac
 
         self.update_slit_metadata(result)
-        if result.meta.exposure.type.lower() != 'mir_lrs-fixedslit':
+        if result.meta.exposure.type.lower() == 'mir_lrs-fixedslit':
+            s_region_model1 = input_models[0].meta.wcsinfo.s_region
+            s_region = find_miri_lrs_sregion(s_region_model1, result.meta.wcs)
+            result.meta.wcsinfo.s_region = s_region
+            self.log.info(f'Updating S_REGION: {s_region}.')
+        else:
             update_s_region_spectral(result)
-
         return result
 
     def update_slit_metadata(self, model):
diff --git a/jwst/resample/resample_utils.py b/jwst/resample/resample_utils.py
index b4ac1bac32..abff0d6184 100644
--- a/jwst/resample/resample_utils.py
+++ b/jwst/resample/resample_utils.py
@@ -9,11 +9,15 @@
 from drizzle.utils import decode_context as _drizzle_decode_context
 
 from stdatamodels.jwst.datamodels.dqflags import pixel
+from astropy.coordinates import SkyCoord
 
 from stcal.alignment.util import (
     compute_scale,
     wcs_from_sregions,
 )
+from gwcs import wcstools
+
+from stcal.alignment.util import compute_s_region_keyword
 from stcal.resample import UnsupportedWCSError
 from stcal.resample.utils import compute_mean_pixel_area
 from stcal.resample.utils import (
@@ -462,3 +466,86 @@ def load_custom_wcs(asdf_wcs_file, output_shape=None):
         "pixel_scale": user_pixel_scale,
     }
     return wcs_dict
+
+
+def find_miri_lrs_sregion(sregion_model1, wcs):
+    """ Find s region for MIRI LRS resampled data.
+
+    Parameters
+    ----------
+    sregion_model1 : string
+        s_regions of the first input model
+    wcs : gwcs.WCS
+        Spatial/spectral WCS.
+
+    Returns
+    -------
+    sregion : string
+        s_region for the resample data.
+    """
+    # use the first sregion to set the width of the slit
+
+    spatial_box = sregion_model1
+    s = spatial_box.split(' ')
+    a1 = float(s[3])
+    b1 = float(s[4])
+    a2 = float(s[5])
+    b2 = float(s[6])
+    a3 = float(s[7])
+    b3 = float(s[8])
+    a4 = float(s[9])
+    b4 = float(s[10])
+
+    # convert each corner to SkyCoord
+    coord1 = SkyCoord(a1, b1, unit='deg')
+    coord2 = SkyCoord(a2, b2, unit='deg')
+    coord3 = SkyCoord(a3, b3, unit='deg')
+    coord4 = SkyCoord(a4, b4, unit='deg')
+
+    # Find the distance between the corners
+    # corners are counter clockwize from 1,2,3,4
+    sep1 = coord1.separation(coord2)
+    sep2 = coord2.separation(coord3)
+    sep3 = coord3.separation(coord4)
+    sep4 = coord4.separation(coord1)
+
+    # use the separation values so we can find the min value later
+    sep = [sep1.value, sep2.value, sep3.value, sep4.value]
+
+    # the minimum separation is the slit width
+    min_sep = np.min(sep)
+    min_sep = min_sep* u.deg # set the units to degrees
+
+    log.info(f'Estimated MIRI LRS slit width: {min_sep*3600} arcsec.')
+    # now use the combined WCS to map all pixels to the slit center
+    bbox = wcs.bounding_box
+    grid = wcstools.grid_from_bounding_box(bbox)
+    ra, dec, _ = np.array(wcs(*grid))
+    ra = ra.flatten()
+    dec = dec.flatten()
+    # ra and dec are the values along the output resampled slit center
+    # using the first point and last point find the position angle
+    star1 = SkyCoord(ra[0]*u.deg, dec[0]*u.deg, frame='icrs')
+    star2 = SkyCoord(ra[-1]*u.deg, dec[-1]*u.deg, frame='icrs')
+    position_angle = star1.position_angle(star2).to(u.deg)
+
+    # 90 degrees to the position angle of the slit will define s_region
+    pos_angle = position_angle - 90.0*u.deg
+
+    star_c1 = star1.directional_offset_by(pos_angle, min_sep/2)
+    star_c2 = star1.directional_offset_by(pos_angle, -min_sep/2)
+    star_c3 = star2.directional_offset_by(pos_angle, min_sep/2)
+    star_c4 = star2.directional_offset_by(pos_angle, -min_sep/2)
+
+    # set  these values to footprint
+    # ra,dec corners are in counter-clockwise direction
+    footprint = [star_c1.ra.value, star_c1.dec.value,
+                 star_c3.ra.value, star_c3.dec.value,
+                 star_c4.ra.value, star_c4.dec.value,
+                 star_c2.ra.value, star_c2.dec.value]
+    footprint = np.array(footprint)
+    s_region = compute_s_region_keyword(footprint)
+    return s_region
+
+
+

From 4dc9a44bcee1699852e62f0592da29817c17a73b Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Thu, 20 Feb 2025 09:14:37 -0700
Subject: [PATCH 08/21] working on unit tests

---
 jwst/resample/tests/test_resample_step.py | 22 +++++++++++++++-------
 1 file changed, 15 insertions(+), 7 deletions(-)

diff --git a/jwst/resample/tests/test_resample_step.py b/jwst/resample/tests/test_resample_step.py
index 558e6f0f16..3b4e4ffd05 100644
--- a/jwst/resample/tests/test_resample_step.py
+++ b/jwst/resample/tests/test_resample_step.py
@@ -44,6 +44,9 @@ def _set_photom_kwd(im):
 def miri_rate_model():
     xsize = 72
     ysize = 416
+    sregion  = 'POLYGON ICRS  10.355323877 -22.353560934 10.355437846 -22.353464295 ' + \
+    '10.354477543 -22.352498313 10.354363599 -22.352595345'
+
     shape = (ysize, xsize)
     im = ImageModel(shape)
     im.data += 5
@@ -55,7 +58,8 @@ def miri_rate_model():
         'v2_ref': -453.5134,
         'v3_ref': -373.4826,
         'v3yangle': 0.0,
-        'vparity': -1}
+        'vparity': -1,
+        's_region': sregion}
     im.meta.instrument = {
         'detector': 'MIRIMAGE',
         'filter': 'P750L',
@@ -460,7 +464,8 @@ def test_pixel_scale_ratio_imaging(nircam_rate, ratio):
 
 
 @pytest.mark.parametrize("units", ["MJy", "MJy/sr"])
-@pytest.mark.parametrize("ratio", [0.7, 1.0, 1.3])
+#@pytest.mark.parametrize("ratio", [0.7, 1.0, 1.3])
+@pytest.mark.parametrize("ratio", [0.7])
 def test_pixel_scale_ratio_spec_miri(miri_cal, ratio, units):
     miri_cal.meta.bunit_data = units
 
@@ -468,10 +473,11 @@ def test_pixel_scale_ratio_spec_miri(miri_cal, ratio, units):
     input_scale = compute_spectral_pixel_scale(miri_cal.meta.wcs, disp_axis=2)
     pscale = 3600.0 * input_scale / ratio
 
-    result1 = ResampleSpecStep.call(miri_cal)
-    result2 = ResampleSpecStep.call(miri_cal, pixel_scale_ratio=ratio)
-    result3 = ResampleSpecStep.call(miri_cal, pixel_scale=pscale)
+    result1 = ResampleSpecStep.call(miri_cal, output_file='test1.fits')
+    result2 = ResampleSpecStep.call(miri_cal, pixel_scale_ratio=ratio, output_file='test2.fits')
+    result3 = ResampleSpecStep.call(miri_cal, pixel_scale=pscale,output_file='test_3.fits')
 
+    print(miri_cal.meta.wcsinfo.s_region)
     # pixel_scale and pixel_scale_ratio should be equivalent
     nn = np.isnan(result2.data) | np.isnan(result3.data)
     assert np.allclose(result2.data[~nn], result3.data[~nn])
@@ -482,6 +488,7 @@ def test_pixel_scale_ratio_spec_miri(miri_cal, ratio, units):
     assert result1.data.shape[0] == result2.data.shape[0]
 
     # spatial dimension is scaled
+    print(result1.data.shape[1], result2.data.shape[1])
     assert np.isclose(result1.data.shape[1], result2.data.shape[1] / ratio, atol=1)
 
     # data is non-trivial
@@ -518,8 +525,9 @@ def test_pixel_scale_ratio_spec_miri(miri_cal, ratio, units):
 
 
 @pytest.mark.parametrize("units", ["MJy", "MJy/sr"])
-@pytest.mark.parametrize("ratio", [0.7, 1.0, 1.3])
-def test_pixel_scale_ratio_spec_miri_pair(miri_rate_pair, ratio, units):
+#@pytest.mark.parametrize("ratio", [0.7, 1.0, 1.3])
+@pytest.mark.parametrize("ratio", [0.7])
+def test_pixel_scale_ratio_1spec_miri_pair(miri_rate_pair, ratio, units):
     im1, im2 = miri_rate_pair
     _set_photom_kwd(im1)
     _set_photom_kwd(im2)

From ab57dd159b92d6597d23fd8d8b362535ee61ae5c Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Thu, 20 Feb 2025 09:38:15 -0700
Subject: [PATCH 09/21] fix unit tests

---
 jwst/resample/resample_spec.py            |  5 ++---
 jwst/resample/tests/test_resample_step.py | 18 +++++++++---------
 2 files changed, 11 insertions(+), 12 deletions(-)

diff --git a/jwst/resample/resample_spec.py b/jwst/resample/resample_spec.py
index 9b6cb08ef5..f2f59f0614 100644
--- a/jwst/resample/resample_spec.py
+++ b/jwst/resample/resample_spec.py
@@ -578,7 +578,6 @@ def build_interpolated_output_wcs(self, input_models, pixel_scale_ratio=1.0):
         all_dec_slit = []
         xstop = 0
 
-        #sregion_list = []
         all_wcs = [m.meta.wcs for m in input_models]
         for im, model in enumerate(input_models):
             wcs = model.meta.wcs
@@ -767,9 +766,9 @@ def build_interpolated_output_wcs(self, input_models, pixel_scale_ratio=1.0):
         slope_sign_x = np.sign(pix_to_xtan.slope)
 
         if swap_xy:
-            ny = int(np.ceil(diff_y / pix_to_tan_slope_y)) + 1
+            ny = int(np.ceil(diff_y / pix_to_tan_slope_y))# + 1
         else:
-            ny = int(np.ceil(diff_x / pix_to_tan_slope_x)) + 1
+            ny = int(np.ceil(diff_x / pix_to_tan_slope_x))# + 1
 
         offset_y  = (ny)/2 * pix_to_tan_slope_y
         offset_x  = (ny)/2 * pix_to_tan_slope_x
diff --git a/jwst/resample/tests/test_resample_step.py b/jwst/resample/tests/test_resample_step.py
index 3b4e4ffd05..f0cade0784 100644
--- a/jwst/resample/tests/test_resample_step.py
+++ b/jwst/resample/tests/test_resample_step.py
@@ -119,6 +119,8 @@ def miri_cal(miri_rate):
 def miri_rate_zero_crossing():
     xsize = 1032
     ysize = 1024
+    sregion  = 'POLYGON ICRS  10.355323877 -22.353560934 10.355437846 -22.353464295 ' + \
+        '10.354477543 -22.352498313 10.354363599 -22.352595345'
     shape = (ysize, xsize)
     im = ImageModel(shape)
     im.var_rnoise = np.random.random(shape)
@@ -129,7 +131,8 @@ def miri_rate_zero_crossing():
         'v2_ref': -415.0690466121227,
         'v3_ref': -400.575920398547,
         'v3yangle': 0.0,
-        'vparity': -1}
+        'vparity': -1,
+        's_region': sregion}
     im.meta.instrument = {
         'detector': 'MIRIMAGE',
         'filter': 'P750L',
@@ -464,8 +467,7 @@ def test_pixel_scale_ratio_imaging(nircam_rate, ratio):
 
 
 @pytest.mark.parametrize("units", ["MJy", "MJy/sr"])
-#@pytest.mark.parametrize("ratio", [0.7, 1.0, 1.3])
-@pytest.mark.parametrize("ratio", [0.7])
+@pytest.mark.parametrize("ratio", [0.7, 1.0, 1.3])
 def test_pixel_scale_ratio_spec_miri(miri_cal, ratio, units):
     miri_cal.meta.bunit_data = units
 
@@ -473,11 +475,10 @@ def test_pixel_scale_ratio_spec_miri(miri_cal, ratio, units):
     input_scale = compute_spectral_pixel_scale(miri_cal.meta.wcs, disp_axis=2)
     pscale = 3600.0 * input_scale / ratio
 
-    result1 = ResampleSpecStep.call(miri_cal, output_file='test1.fits')
-    result2 = ResampleSpecStep.call(miri_cal, pixel_scale_ratio=ratio, output_file='test2.fits')
-    result3 = ResampleSpecStep.call(miri_cal, pixel_scale=pscale,output_file='test_3.fits')
+    result1 = ResampleSpecStep.call(miri_cal)
+    result2 = ResampleSpecStep.call(miri_cal, pixel_scale_ratio=ratio)
+    result3 = ResampleSpecStep.call(miri_cal, pixel_scale=pscale)
 
-    print(miri_cal.meta.wcsinfo.s_region)
     # pixel_scale and pixel_scale_ratio should be equivalent
     nn = np.isnan(result2.data) | np.isnan(result3.data)
     assert np.allclose(result2.data[~nn], result3.data[~nn])
@@ -525,8 +526,7 @@ def test_pixel_scale_ratio_spec_miri(miri_cal, ratio, units):
 
 
 @pytest.mark.parametrize("units", ["MJy", "MJy/sr"])
-#@pytest.mark.parametrize("ratio", [0.7, 1.0, 1.3])
-@pytest.mark.parametrize("ratio", [0.7])
+@pytest.mark.parametrize("ratio", [0.7, 1.0, 1.3])
 def test_pixel_scale_ratio_1spec_miri_pair(miri_rate_pair, ratio, units):
     im1, im2 = miri_rate_pair
     _set_photom_kwd(im1)

From 5fde36a5e5ba622e47df8c1f83a4e816550e08c4 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Thu, 20 Feb 2025 09:43:34 -0700
Subject: [PATCH 10/21] remove comment

---
 jwst/resample/resample_spec.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/jwst/resample/resample_spec.py b/jwst/resample/resample_spec.py
index f2f59f0614..b19ed16ad5 100644
--- a/jwst/resample/resample_spec.py
+++ b/jwst/resample/resample_spec.py
@@ -766,9 +766,9 @@ def build_interpolated_output_wcs(self, input_models, pixel_scale_ratio=1.0):
         slope_sign_x = np.sign(pix_to_xtan.slope)
 
         if swap_xy:
-            ny = int(np.ceil(diff_y / pix_to_tan_slope_y))# + 1
+            ny = int(np.ceil(diff_y / pix_to_tan_slope_y))
         else:
-            ny = int(np.ceil(diff_x / pix_to_tan_slope_x))# + 1
+            ny = int(np.ceil(diff_x / pix_to_tan_slope_x))
 
         offset_y  = (ny)/2 * pix_to_tan_slope_y
         offset_x  = (ny)/2 * pix_to_tan_slope_x

From cb8b93c94a8e6ff9971bb2b8f61c034c62aab927 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Wed, 26 Feb 2025 09:16:34 -0700
Subject: [PATCH 11/21] updates from review

---
 jwst/assign_wcs/miri.py                   | 3 +--
 jwst/assign_wcs/util.py                   | 7 ++++---
 jwst/resample/resample_utils.py           | 1 -
 jwst/resample/tests/test_resample_step.py | 1 -
 4 files changed, 5 insertions(+), 7 deletions(-)

diff --git a/jwst/assign_wcs/miri.py b/jwst/assign_wcs/miri.py
index 0b151f80ba..93b152d1c4 100644
--- a/jwst/assign_wcs/miri.py
+++ b/jwst/assign_wcs/miri.py
@@ -255,8 +255,7 @@ def lrs_xytoabl(input_model, reference_files):
     if input_model.meta.exposure.type.lower() == 'mir_lrs-fixedslit':
         zero_point = refmodel.meta.x_ref - 1, refmodel.meta.y_ref - 1
     elif input_model.meta.exposure.type.lower() == 'mir_lrs-slitless':
-        zero_point = refmodel.meta.x_ref_slitless - 1, \
-            refmodel.meta.y_ref_slitless - 1
+        zero_point = refmodel.meta.x_ref_slitless - 1, refmodel.meta.y_ref_slitless - 1
         # Transform to slitless subarray from full array
         zero_point = subarray2full.inverse(zero_point[0], zero_point[1])
 
diff --git a/jwst/assign_wcs/util.py b/jwst/assign_wcs/util.py
index 263b3987c8..f7fbfcd9f2 100644
--- a/jwst/assign_wcs/util.py
+++ b/jwst/assign_wcs/util.py
@@ -809,7 +809,7 @@ def update_s_region_imaging(model):
 
 def update_s_region_lrs(model, reference_files):
     """
-    Determine ``S_REGION`` using V2,V3 of the slit corners from reference file.
+    Update ``S_REGION`` using V2,V3 of the slit corners from reference file.
 
     Parameters
     ----------
@@ -821,7 +821,8 @@ def update_s_region_lrs(model, reference_files):
     Returns
     -------
     None
-        s_region for model is updated in place .
+
+        s_region for model is updated in place.
 
     """
 
@@ -846,7 +847,7 @@ def update_s_region_lrs(model, reference_files):
         log.info("The V2,V3 coordinates of the MIRI LRS-Fixed slit contains NaN values.")
         log.info("The s_region will not be updated")
 
-    lam = 7.0 # wavelength does not matter for s region assign a value in range of LRS
+    lam = 7.0 # wavelength does not matter for s region so jwst assign a value in range of LRS
     s = model.meta.wcs.transform('v2v3', 'world', v2, v3, lam)
     a = s[0]
     b = s[1]
diff --git a/jwst/resample/resample_utils.py b/jwst/resample/resample_utils.py
index abff0d6184..272ef314a0 100644
--- a/jwst/resample/resample_utils.py
+++ b/jwst/resample/resample_utils.py
@@ -484,7 +484,6 @@ def find_miri_lrs_sregion(sregion_model1, wcs):
         s_region for the resample data.
     """
     # use the first sregion to set the width of the slit
-
     spatial_box = sregion_model1
     s = spatial_box.split(' ')
     a1 = float(s[3])
diff --git a/jwst/resample/tests/test_resample_step.py b/jwst/resample/tests/test_resample_step.py
index f0cade0784..bc23e8f920 100644
--- a/jwst/resample/tests/test_resample_step.py
+++ b/jwst/resample/tests/test_resample_step.py
@@ -489,7 +489,6 @@ def test_pixel_scale_ratio_spec_miri(miri_cal, ratio, units):
     assert result1.data.shape[0] == result2.data.shape[0]
 
     # spatial dimension is scaled
-    print(result1.data.shape[1], result2.data.shape[1])
     assert np.isclose(result1.data.shape[1], result2.data.shape[1] / ratio, atol=1)
 
     # data is non-trivial

From 48c5c84c181882002024522fb67857624d2091a3 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Mon, 3 Mar 2025 10:39:31 -0700
Subject: [PATCH 12/21] pyproject pointing to jwst/main

---
 pyproject.toml | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index f513e9534b..87a88914b1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -34,8 +34,7 @@ dependencies = [
     "scikit-image>=0.20.0",
     "scipy>=1.14.1",
     "spherical-geometry>=1.2.22",
-    #"stdatamodels @ git+https://github.com/spacetelescope/stdatamodels.git@main",
-    "stdatamodels @ git+https://github.com/jemorrison/stdatamodels.git@MIRI_LRS_Specwcs",
+    "stdatamodels @ git+https://github.com/spacetelescope/stdatamodels.git@main",
     "stcal @ git+https://github.com/spacetelescope/stcal.git@main",
     "stpipe>=0.8.0,<0.9.0",
     "stsci.imagestats>=1.6.3",

From 550a1e173fd216ce8d3e75af2e5144ea20139276 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Mon, 3 Mar 2025 12:33:50 -0700
Subject: [PATCH 13/21] update doc

---
 jwst/assign_wcs/util.py | 10 ++--------
 1 file changed, 2 insertions(+), 8 deletions(-)

diff --git a/jwst/assign_wcs/util.py b/jwst/assign_wcs/util.py
index f7fbfcd9f2..4048147d6e 100644
--- a/jwst/assign_wcs/util.py
+++ b/jwst/assign_wcs/util.py
@@ -811,21 +811,15 @@ def update_s_region_lrs(model, reference_files):
     """
     Update ``S_REGION`` using V2,V3 of the slit corners from reference file.
 
+    s_region for model is updated in place.
+
     Parameters
     ----------
     model : DataModel
         Input model
     reference_files : list
         List of reference files for assign_wcs.
-
-    Returns
-    -------
-    None
-
-        s_region for model is updated in place.
-
     """
-
     refmodel = MiriLRSSpecwcsModel(reference_files['specwcs'])
 
     v2vert1 = refmodel.meta.v2_vert1

From fba36c909a5516cc79047b1890beafcd7fc5338e Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Mon, 3 Mar 2025 16:06:43 -0700
Subject: [PATCH 14/21] remove blank line

---
 jwst/assign_wcs/assign_wcs_step.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/jwst/assign_wcs/assign_wcs_step.py b/jwst/assign_wcs/assign_wcs_step.py
index 9c2a6c268f..8d5ba034d8 100755
--- a/jwst/assign_wcs/assign_wcs_step.py
+++ b/jwst/assign_wcs/assign_wcs_step.py
@@ -95,7 +95,6 @@ def process(self, input, *args, **kwargs):
                     message = "MSA metadata file (MSAMETFL) is required for NRS_MSASPEC exposures."
                     log.error(message)
                     raise MSAFileError(message)
-
             slit_y_range = [self.slit_y_low, self.slit_y_high]
             result = load_wcs(input_model, reference_file_names, slit_y_range)
 

From 1bf6725b7743eaee904aec82af12570df109e2bd Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Mon, 3 Mar 2025 16:24:32 -0700
Subject: [PATCH 15/21] fix error with miri slitless

---
 jwst/assign_wcs/assign_wcs.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/jwst/assign_wcs/assign_wcs.py b/jwst/assign_wcs/assign_wcs.py
index db9eeee00a..7d91fee238 100644
--- a/jwst/assign_wcs/assign_wcs.py
+++ b/jwst/assign_wcs/assign_wcs.py
@@ -73,13 +73,14 @@ def load_wcs(input_model, reference_files={}, nrs_slit_y_range=None):
 
         if output_model.meta.exposure.type.lower() not in exclude_types:
             imaging_types = IMAGING_TYPES.copy()
-            imaging_lrs_types = ['mir_lrs-fixedslit']
             imaging_types.update(['mir_lrs-slitless'])
+            imaging_lrs_types = ['mir_lrs-fixedslit']
             if output_model.meta.exposure.type.lower() in imaging_lrs_types:
                 # uses slits corners in V2, V3 that are read in from the
                 # lrs specwcs reference file
                 update_s_region_lrs(output_model, reference_files)
-            elif output_model.meta.exposure.type.lower() in IMAGING_TYPES:
+            elif output_model.meta.exposure.type.lower() in imaging_types:
+                print(output_model.meta.exposure.type.lower())
                 try:
                     update_s_region_imaging(output_model)
                 except Exception as exc:

From 69b345aa614cf8dd8b0183ab03026b85ab04b8e8 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Tue, 4 Mar 2025 09:52:02 -0700
Subject: [PATCH 16/21] testing reading reference file

---
 jwst/assign_wcs/miri.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/jwst/assign_wcs/miri.py b/jwst/assign_wcs/miri.py
index 93b152d1c4..3284c4f35d 100644
--- a/jwst/assign_wcs/miri.py
+++ b/jwst/assign_wcs/miri.py
@@ -258,7 +258,7 @@ def lrs_xytoabl(input_model, reference_files):
         zero_point = refmodel.meta.x_ref_slitless - 1, refmodel.meta.y_ref_slitless - 1
         # Transform to slitless subarray from full array
         zero_point = subarray2full.inverse(zero_point[0], zero_point[1])
-
+    print('zero point', zero_point)
 
     # Figure out the typical along-slice pixel scale at the center of the slit
     v2_cen, v3_cen = subarray_dist(zero_point[0], zero_point[1])
@@ -276,6 +276,10 @@ def lrs_xytoabl(input_model, reference_files):
     y0 = refmodel.wavetable.y0
     x1 = refmodel.wavetable.x1
     y2 = refmodel.wavetable.y2
+    print('xcen', xcen)
+    print('ycen', ycen)
+    print('wavetab', wavetab)
+
     refmodel.close()
     # If in fixed slit mode, define the bounding box using the corner locations provided in
     # the CDP reference file.
@@ -291,6 +295,8 @@ def lrs_xytoabl(input_model, reference_files):
         bb_sub = ((input_model.meta.subarray.xstart - 1 + 4 - 0.5, input_model.meta.subarray.xsize - 1 + 0.5),
                   (np.floor(y2.min() + zero_point[1]) - 0.5, np.ceil(y0.max() + zero_point[1]) + 0.5))
 
+    print('bb_sub', bb_sub)
+
     # Now deal with the fact that the spectral trace isn't perfectly up and down along detector.
     # This information is contained in the xcenter/ycenter values in the CDP table, but we'll handle it
     # as a simple x shift using a linear fit to this relation provided by the CDP.
@@ -405,6 +411,7 @@ def lrs_abltov2v3l(input_model, reference_files):
         # Transform to slitless subarray from full array
         zero_point = subarray2full.inverse(zero_point[0], zero_point[1])
 
+    print('in lrs_abltov2v3l zero pt',zero_point)
     refmodel.close()
     # Figure out the typical along-slice pixel scale at the center of the slit
     v2_cen, v3_cen = subarray_dist(zero_point[0], zero_point[1])

From d17903a040fbe1b222c90818739560701c33eae6 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Tue, 4 Mar 2025 15:33:36 -0700
Subject: [PATCH 17/21] added bb from before

---
 jwst/assign_wcs/miri.py | 18 +++++++-----------
 1 file changed, 7 insertions(+), 11 deletions(-)

diff --git a/jwst/assign_wcs/miri.py b/jwst/assign_wcs/miri.py
index 3284c4f35d..e5ab0ea11b 100644
--- a/jwst/assign_wcs/miri.py
+++ b/jwst/assign_wcs/miri.py
@@ -5,7 +5,6 @@
 from astropy.modeling import models
 from astropy import coordinates as coord
 from astropy import units as u
-
 from scipy.interpolate import UnivariateSpline
 import gwcs.coordinate_frames as cf
 from gwcs import selector
@@ -258,7 +257,6 @@ def lrs_xytoabl(input_model, reference_files):
         zero_point = refmodel.meta.x_ref_slitless - 1, refmodel.meta.y_ref_slitless - 1
         # Transform to slitless subarray from full array
         zero_point = subarray2full.inverse(zero_point[0], zero_point[1])
-    print('zero point', zero_point)
 
     # Figure out the typical along-slice pixel scale at the center of the slit
     v2_cen, v3_cen = subarray_dist(zero_point[0], zero_point[1])
@@ -276,10 +274,6 @@ def lrs_xytoabl(input_model, reference_files):
     y0 = refmodel.wavetable.y0
     x1 = refmodel.wavetable.x1
     y2 = refmodel.wavetable.y2
-    print('xcen', xcen)
-    print('ycen', ycen)
-    print('wavetab', wavetab)
-
     refmodel.close()
     # If in fixed slit mode, define the bounding box using the corner locations provided in
     # the CDP reference file.
@@ -295,8 +289,6 @@ def lrs_xytoabl(input_model, reference_files):
         bb_sub = ((input_model.meta.subarray.xstart - 1 + 4 - 0.5, input_model.meta.subarray.xsize - 1 + 0.5),
                   (np.floor(y2.min() + zero_point[1]) - 0.5, np.ceil(y0.max() + zero_point[1]) + 0.5))
 
-    print('bb_sub', bb_sub)
-
     # Now deal with the fact that the spectral trace isn't perfectly up and down along detector.
     # This information is contained in the xcenter/ycenter values in the CDP table, but we'll handle it
     # as a simple x shift using a linear fit to this relation provided by the CDP.
@@ -312,7 +304,11 @@ def lrs_xytoabl(input_model, reference_files):
     # This function will give slit dX as a function of Y subarray pixel value
     dxmodel = models.Tabular1D(lookup_table=xshiftref, points=ycen_subarray, name='xshiftref',
                                  bounds_error=False, fill_value=np.nan)
-
+    print('bb_sub', bb_sub)
+    #### CHECK WITH NADIA
+    if input_model.meta.exposure.type.lower() == 'mir_lrs-fixedslit':
+        bb_sub = (bb_sub[0], (dxmodel.points[0].min(), dxmodel.points[0].max()))
+    print('second bb_sub', bb_sub)
     # Fit for the wavelength as a function of Y
     # Reverse the vectors so that yinv is increasing (needed for spline fitting function)
     # Spline fit with enforced smoothness
@@ -322,7 +318,6 @@ def lrs_xytoabl(input_model, reference_files):
     # This model will now give the wavelength corresponding to a given Y subarray pixel value
     wavemodel = models.Tabular1D(lookup_table=wavereference, points=ycen_subarray, name='waveref',
                                  bounds_error=False, fill_value=np.nan)
-
     # Wavelength barycentric correction
     try:
         velosys = input_model.meta.wcsinfo.velosys
@@ -380,6 +375,7 @@ def lrs_xytoabl(input_model, reference_files):
 
     return dettoabl
 
+
 def lrs_abltov2v3l(input_model, reference_files):
     """
     The second part of LRS-FIXEDSLIT and LRS-SLITLESS WCS pipeline.
@@ -411,7 +407,6 @@ def lrs_abltov2v3l(input_model, reference_files):
         # Transform to slitless subarray from full array
         zero_point = subarray2full.inverse(zero_point[0], zero_point[1])
 
-    print('in lrs_abltov2v3l zero pt',zero_point)
     refmodel.close()
     # Figure out the typical along-slice pixel scale at the center of the slit
     v2_cen, v3_cen = subarray_dist(zero_point[0], zero_point[1])
@@ -442,6 +437,7 @@ def lrs_abltov2v3l(input_model, reference_files):
 
     return abl_to_v2v3l
 
+
 def ifu(input_model, reference_files):
     """
     The MIRI MRS WCS pipeline.

From 7ac85e6b135b4eb9335fe6e2b4ff61f3420b1516 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Wed, 5 Mar 2025 09:33:49 -0700
Subject: [PATCH 18/21] added bbox using tabular model

---
 jwst/assign_wcs/miri.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/jwst/assign_wcs/miri.py b/jwst/assign_wcs/miri.py
index e5ab0ea11b..9cd1bb7ea8 100644
--- a/jwst/assign_wcs/miri.py
+++ b/jwst/assign_wcs/miri.py
@@ -304,11 +304,9 @@ def lrs_xytoabl(input_model, reference_files):
     # This function will give slit dX as a function of Y subarray pixel value
     dxmodel = models.Tabular1D(lookup_table=xshiftref, points=ycen_subarray, name='xshiftref',
                                  bounds_error=False, fill_value=np.nan)
-    print('bb_sub', bb_sub)
-    #### CHECK WITH NADIA
     if input_model.meta.exposure.type.lower() == 'mir_lrs-fixedslit':
         bb_sub = (bb_sub[0], (dxmodel.points[0].min(), dxmodel.points[0].max()))
-    print('second bb_sub', bb_sub)
+
     # Fit for the wavelength as a function of Y
     # Reverse the vectors so that yinv is increasing (needed for spline fitting function)
     # Spline fit with enforced smoothness

From 491ab17721cc29c05bcae29fb38216b779ecf1f9 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Thu, 6 Mar 2025 11:07:59 -0700
Subject: [PATCH 19/21] update test

---
 jwst/assign_wcs/miri.py                  | 1 -
 jwst/regtest/test_miri_lrs_slit_spec2.py | 4 ++--
 2 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/jwst/assign_wcs/miri.py b/jwst/assign_wcs/miri.py
index 9cd1bb7ea8..ed782f291a 100644
--- a/jwst/assign_wcs/miri.py
+++ b/jwst/assign_wcs/miri.py
@@ -306,7 +306,6 @@ def lrs_xytoabl(input_model, reference_files):
                                  bounds_error=False, fill_value=np.nan)
     if input_model.meta.exposure.type.lower() == 'mir_lrs-fixedslit':
         bb_sub = (bb_sub[0], (dxmodel.points[0].min(), dxmodel.points[0].max()))
-
     # Fit for the wavelength as a function of Y
     # Reverse the vectors so that yinv is increasing (needed for spline fitting function)
     # Spline fit with enforced smoothness
diff --git a/jwst/regtest/test_miri_lrs_slit_spec2.py b/jwst/regtest/test_miri_lrs_slit_spec2.py
index 4d819fe46e..034be6a144 100644
--- a/jwst/regtest/test_miri_lrs_slit_spec2.py
+++ b/jwst/regtest/test_miri_lrs_slit_spec2.py
@@ -92,7 +92,7 @@ def test_miri_lrs_slit_wcs(run_pipeline, rtdata_module, fitsdiff_default_kwargs)
         # floats rather than a string. Test footprint
         sregion = im.meta.wcsinfo.s_region
         sregion_test = im_truth.meta.wcsinfo.s_region
-        footprint=util._sregion_to_footprint(sregion)
-        footprint_test = util._sregion_to_footprint(sregion_test)
+        footprint=util.sregion_to_footprint(sregion)
+        footprint_test = util.sregion_to_footprint(sregion_test)
         assert_allclose(footprint, footprint_test)
 

From 0ef2da052607ee218d287b8ceec2845232d47e02 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Thu, 6 Mar 2025 14:31:40 -0700
Subject: [PATCH 20/21] update test

---
 jwst/regtest/test_miri_lrs_slit_spec2.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/jwst/regtest/test_miri_lrs_slit_spec2.py b/jwst/regtest/test_miri_lrs_slit_spec2.py
index 034be6a144..1c547e12e1 100644
--- a/jwst/regtest/test_miri_lrs_slit_spec2.py
+++ b/jwst/regtest/test_miri_lrs_slit_spec2.py
@@ -67,7 +67,6 @@ def test_miri_lrs_extract1d_from_cal(run_pipeline, rtdata_module, fitsdiff_defau
 @pytest.mark.bigdata
 def test_miri_lrs_slit_wcs(run_pipeline, rtdata_module, fitsdiff_default_kwargs):
     rtdata = rtdata_module
-
     # get input assign_wcs and truth file
     output = "jw01530005001_03103_00001_mirimage_assign_wcs.fits"
     rtdata.output = output

From 8a4396847f944e9b637dc3cd2bbfd6261d292eb0 Mon Sep 17 00:00:00 2001
From: jemorrison <morrison@as.arizona.edu>
Date: Fri, 7 Mar 2025 12:19:26 -0700
Subject: [PATCH 21/21] updates from review

---
 jwst/assign_wcs/assign_wcs.py   | 1 -
 jwst/assign_wcs/util.py         | 2 +-
 jwst/resample/resample_utils.py | 2 +-
 3 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/jwst/assign_wcs/assign_wcs.py b/jwst/assign_wcs/assign_wcs.py
index 7d91fee238..52354d29aa 100644
--- a/jwst/assign_wcs/assign_wcs.py
+++ b/jwst/assign_wcs/assign_wcs.py
@@ -80,7 +80,6 @@ def load_wcs(input_model, reference_files={}, nrs_slit_y_range=None):
                 # lrs specwcs reference file
                 update_s_region_lrs(output_model, reference_files)
             elif output_model.meta.exposure.type.lower() in imaging_types:
-                print(output_model.meta.exposure.type.lower())
                 try:
                     update_s_region_imaging(output_model)
                 except Exception as exc:
diff --git a/jwst/assign_wcs/util.py b/jwst/assign_wcs/util.py
index 4048147d6e..40dadd1228 100644
--- a/jwst/assign_wcs/util.py
+++ b/jwst/assign_wcs/util.py
@@ -841,7 +841,7 @@ def update_s_region_lrs(model, reference_files):
         log.info("The V2,V3 coordinates of the MIRI LRS-Fixed slit contains NaN values.")
         log.info("The s_region will not be updated")
 
-    lam = 7.0 # wavelength does not matter for s region so jwst assign a value in range of LRS
+    lam = 7.0 # wavelength does not matter for s region so just assign a value in range of LRS
     s = model.meta.wcs.transform('v2v3', 'world', v2, v3, lam)
     a = s[0]
     b = s[1]
diff --git a/jwst/resample/resample_utils.py b/jwst/resample/resample_utils.py
index 272ef314a0..b6e95ff7ca 100644
--- a/jwst/resample/resample_utils.py
+++ b/jwst/resample/resample_utils.py
@@ -502,7 +502,7 @@ def find_miri_lrs_sregion(sregion_model1, wcs):
     coord4 = SkyCoord(a4, b4, unit='deg')
 
     # Find the distance between the corners
-    # corners are counter clockwize from 1,2,3,4
+    # corners are counterclockwize from 1,2,3,4
     sep1 = coord1.separation(coord2)
     sep2 = coord2.separation(coord3)
     sep3 = coord3.separation(coord4)