JP-2500 Allow cube build to accept a tangent point, position angle an…

…d number of spaxels in x, y cube dimensions (#7882)

CHANGES.rst

@@ -85,6 +85,9 @@ cube_build
 
 - Fix a bug with memory allocation in C extensions when weighting=emsm. [#7847]
 
+- Add options to set ra,dec tangent projection center, position angle and size of cube
+  in tangent plane. [#7882]
+
 datamodels
 ----------
 

docs/jwst/cube_build/arguments.rst

@@ -76,6 +76,21 @@ The following arguments control the size and sampling characteristics of the out
 ``wavemax``
   The maximum wavelength, in microns, to use in constructing the IFU cube.
 
+``ra_center``
+  Right ascension center, in decimal degrees, of the IFU cube that defines the location of xi/eta tangent plane projection origin.
+
+``dec_center``
+  Declination center, in decimal degrees, of the IFU cube that defines the location of xi/eta tangent plane projection origin.
+
+``cube_pa``
+  The position angle of the IFU cube in decimal degrees (E from N).
+
+``nspax_x``
+  The odd integer number of spaxels to use in the x dimension of the tangent plane.
+
+``nspax_y``
+  The odd integer number of spaxels to use in the y dimension of the tangent plane.
+
 ``coord_system [string]``
   The default IFU cubes are built on the ra-dec coordinate system (``coord_system=skyalign``). In these cubes north is up 
   and east is left. There are two other coordinate systems an IFU cube can be built on:

jwst/cube_build/cube_build_step.py

@@ -47,6 +47,11 @@ class CubeBuildStep (Step):
          scalew = float(default=0.0) # cube sample size to use for axis 3, microns
          weighting = option('emsm','msm','drizzle',default = 'drizzle') # Type of weighting function
          coord_system = option('skyalign','world','internal_cal','ifualign',default='skyalign') # Output Coordinate system.
+         ra_center = float(default=None) # RA center of the IFU cube
+         dec_center = float(default=None) # Declination center of the IFU cube
+         cube_pa = float(default=None) # The position angle of the desired cube in decimal degrees E from N
+         nspax_x = integer(default=None) # The odd integer number of spaxels to use in the x dimension of cube tangent plane.
+         nspax_y = integer(default=None) # The odd integer number of spaxels to use in the y dimension of cube tangent plane.
          rois = float(default=0.0) # region of interest spatial size, arc seconds
          roiw = float(default=0.0) # region of interest wavelength size, microns
          weight_power = float(default=2.0) # Weighting option to use for Modified Shepard Method
@@ -108,6 +113,19 @@ def process(self, input):
         if self.roiw != 0.0:
             self.log.info(f'Input Wave ROI size {self.roiw}')
 
+        # check that if self.nspax_x or self.nspax_y is provided they must be odd numbers
+        if self.nspax_x is not None:
+            if self.nspax_x % 2 == 0:
+                self.log.info(f'Input nspax_x must be an odd number {self.nspax_x}')
+                self.nspax_x = self.nspax_x + 1
+                self.log.info(f'Updating nspa by 1. New value {self.nspax_x}')
+
+        if self.nspax_y is not None:
+            if self.nspax_y % 2 == 0:
+                self.log.info(f'Input nspax_y must be an odd number {self.nspax_y}')
+                self.nspax_y = self.nspax_y + 1
+                self.log.info(f'Updating nspax_y by 1. New value {self.nspax_y}')
+
         # valid coord_system:
         # 1. skyalign (ra dec) (aka world)
         # 2. ifualign (ifu cube aligned with slicer plane/ MRS local coord system)
@@ -214,9 +232,8 @@ def process(self, input):
             elif instrument == 'MIRI':
                 self.output_type = 'channel'
         self.pars_input['output_type'] = self.output_type
-
         self.log.info(f'Setting output type to: {self.output_type}')
-         
+
 # Read in Cube Parameter Reference file
 # identify what reference file has been associated with these input
 
@@ -247,6 +264,11 @@ def process(self, input):
             'weighting': self.weighting,
             'weight_power': self.weight_power,
             'coord_system': self.pars_input['coord_system'],
+            'ra_center': self.ra_center,
+            'dec_center': self.dec_center,
+            'cube_pa': self.cube_pa,
+            'nspax_x': self.nspax_x,
+            'nspax_y': self.nspax_y,
             'rois': self.rois,
             'roiw': self.roiw,
             'wavemin': self.wavemin,

jwst/cube_build/ifu_cube.py

@@ -64,6 +64,11 @@ def __init__(self,
 
         self.scalexy = pars_cube.get('scalexy')
         self.scalew = pars_cube.get('scalew')
+        self.ra_center = pars_cube.get('ra_center')
+        self.dec_center = pars_cube.get('dec_center')
+        self.cube_pa = pars_cube.get('cube_pa')
+        self.nspax_x = pars_cube.get('nspax_x')
+        self.nspax_y = pars_cube.get('nspax_y')
         self.rois = pars_cube.get('rois')
         self.roiw = pars_cube.get('roiw')
         self.debug_spaxel = pars_cube.get('debug_spaxel')
@@ -231,6 +236,7 @@ def set_geometry(self, corner_a, corner_b, lambda_min, lambda_max):
         """ Based on the ra,dec and wavelength footprint set up the size
         of the cube in the tangent plane projected coordinate system.
 
+
         Parameters
         ----------
         footprint: tuple
@@ -255,8 +261,16 @@ def set_geometry(self, corner_a, corner_b, lambda_min, lambda_max):
         # we have angles in degrees
         ra_ave = circmean(ravalues * u.deg).value
 
-        self.crval1 = ra_ave
-        self.crval2 = dec_ave
+        if self.ra_center is not None:
+            self.crval1 = self.ra_center
+        else:
+            self.crval1 = ra_ave
+
+        if self.dec_center is not None:
+            self.crval2 = self.dec_center
+        else:
+            self.crval2 = dec_ave
+
         rot_angle = self.rot_angle
 
         # find the 4 corners
@@ -269,7 +283,6 @@ def set_geometry(self, corner_a, corner_b, lambda_min, lambda_max):
                                       corner_a[i], corner_b[i], rot_angle)
             xi_corner.append(xi)
             eta_corner.append(eta)
-
         xi_min = min(xi_corner)
         xi_max = max(xi_corner)
         eta_min = min(eta_corner)
@@ -284,6 +297,15 @@ def set_geometry(self, corner_a, corner_b, lambda_min, lambda_max):
         na = math.ceil(xilimit / self.cdelt1) + 1
         nb = math.ceil(etalimit / self.cdelt2) + 1
 
+        # if the user set the nspax_x or nspax_y then redefine na, nb
+        # it is assumed that both values are ODD numbers
+        # We want the central pixel to be the tangent point with na/nb pixels on either
+        # side of central pixel. 
+        if self.nspax_x is not None:
+            na = int(self.nspax_x/2)
+        if self.nspax_y is not None:
+            nb = int(self.nspax_y/2)
+
         xi_min = 0.0 - (na * self.cdelt1) - (self.cdelt1 / 2.0)
         xi_max = (na * self.cdelt1) + (self.cdelt1 / 2.0)
 
@@ -550,7 +572,7 @@ def build_ifucube(self):
 
         if self.spaxel_z == -1 and self.spaxel_x == -1 and self.spaxel_y == -1:
             debug_cube_index = -1
-            
+
         elif(self.spaxel_z < 0 or self.spaxel_x < 0 or self.spaxel_y < 0):
             print('Incorrect input for Debug Spaxel values. Counting starts at 0')
             debug_cube_index = -1
@@ -1209,11 +1231,11 @@ def setup_ifucube_wcs(self):
 
         parameter1 = self.list_par1
         parameter2 = self.list_par2
-
 # ________________________________________________________________________________
-# Define the rotation angle between the ra-dec and alpha beta coord system using
-# the first input model. Use first file in first band to set up rotation angle
+# Define the rotation angle
 
+# If coord_system = ifualign then the angle is between the ra-dec and alpha beta
+# coord system using the first input model. Use first file in first band to set up rotation angle
 # Compute the rotation angle between local IFU system  and RA-DEC
 
         if self.coord_system == 'ifualign':
@@ -1257,6 +1279,12 @@ def setup_ifucube_wcs(self):
             self.rot_angle = 90 + np.arctan2(dra, ddec) * 180. / np.pi
             log.info(f'Rotation angle between ifu and sky: {self.rot_angle}')
 
+        # If coord_system = iskyalign and the user provided a position angle. Define the rotation angle
+        # to be the user provided value.
+
+        if self.coord_system == 'skyalign' and self.cube_pa is not None:
+            self.rot_angle = self.cube_pa
+            log.info(f'Setting rotation angle between ifu and sky: {self.rot_angle}')
 # ________________________________________________________________________________
 # now loop over data and find min and max ranges data covers
 
@@ -1285,8 +1313,9 @@ def setup_ifucube_wcs(self):
                 spectral_found = hasattr(input_model.meta.wcsinfo, 'spectral_region')
                 spatial_found = hasattr(input_model.meta.wcsinfo, 's_region')
                 world = False
-                if self.coord_system == 'skyalign':
+                if self.coord_system == 'skyalign' and self.cube_pa is None:
                     world = True
+
                 # Do not use the default spatial or spectral region found in the wcs if
                 # 1. instrument is MIRI and
                 # 2. Output type is not multi and (not default calspec2) and