diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index ba3f12949b..6e63f5e5af 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -39,7 +39,6 @@ repos: (?x)^( jwst/associations/.* | jwst/coron/.* | - jwst/cube_build/.* | jwst/cube_skymatch/.* | jwst/dark_current/.* | jwst/flatfield/.* | diff --git a/.ruff.toml b/.ruff.toml index 8905406e6c..9d00841a62 100644 --- a/.ruff.toml +++ b/.ruff.toml @@ -21,7 +21,6 @@ docstring-code-format = true exclude = [ "jwst/associations/**.py", "jwst/coron/**.py", - "jwst/cube_build/**.py", "jwst/cube_skymatch/**.py", "jwst/dark_current/**.py", "jwst/flatfield/**.py", @@ -115,7 +114,6 @@ ignore-fully-untyped = true # Turn off annotation checking for fully untyped co ] "jwst/associations/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"] "jwst/coron/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"] -"jwst/cube_build/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"] "jwst/cube_skymatch/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"] "jwst/dark_current/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"] "jwst/flatfield/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"] diff --git a/changes/9347.cube_build.rst b/changes/9347.cube_build.rst new file mode 100644 index 0000000000..85da90c8b1 --- /dev/null +++ b/changes/9347.cube_build.rst @@ -0,0 +1 @@ +Updated cube_build function names to adhere to Python standards and updated some of the exception names. diff --git a/jwst/cube_build/__init__.py b/jwst/cube_build/__init__.py index da6ed8e862..0897ad2d72 100644 --- a/jwst/cube_build/__init__.py +++ b/jwst/cube_build/__init__.py @@ -1,3 +1,5 @@ +"""Spectral Cube building program for JWST IFU data.""" + from .cube_build_step import CubeBuildStep -__all__ = ['CubeBuildStep'] +__all__ = ["CubeBuildStep"] diff --git a/jwst/cube_build/blot_cube_build.py b/jwst/cube_build/blot_cube_build.py index c3e2ba043c..ee487e7bf4 100644 --- a/jwst/cube_build/blot_cube_build.py +++ b/jwst/cube_build/blot_cube_build.py @@ -1,23 +1,23 @@ -""" Main module for blotting sky cube back to detector space -""" import numpy as np import logging from jwst.datamodels import ModelContainer - from ..assign_wcs import nirspec from gwcs import wcstools from jwst.assign_wcs.util import in_ifu_slice from . import instrument_defaults from .blot_median import blot_wrapper # c extension + log = logging.getLogger(__name__) log.setLevel(logging.DEBUG) -class CubeBlot(): +class CubeBlot: + """Main module for blotting a sky cube back to detector space.""" def __init__(self, median_model, input_models): - """Class Blot holds the main variables for blotting sky cube to detector + """ + Initialize main variables for blotting a sky cube to detector space. Information is pulled out of the median sky cube created by a previous run of cube_build in single mode and stored in the ClassBlot.These @@ -26,19 +26,14 @@ def __init__(self, median_model, input_models): data (instrument, channel, band, grating or filter). Parameters - --------- - median_model: ifucube model + ---------- + median_model : IFUCubeModel The median input sky cube is created from a median stack of all the individual input_models mapped to the full IFU cube imprint on the sky. - input_models: data model + input_models : ModelContainer The input models used to create the median sky cube. - - Returns - ------- - CubeBlot class initialized """ - # Pull out the needed information from the Median IFUCube self.median_skycube = median_model self.instrument = median_model.meta.instrument.name @@ -51,13 +46,13 @@ def __init__(self, median_model, input_models): self.par_median_select1 = None self.par_median_select2 = None - if self.instrument == 'MIRI': + if self.instrument == "MIRI": self.channel = median_model.meta.instrument.channel self.subchannel = median_model.meta.instrument.band.lower() self.par_median_select1 = self.channel self.par_median_select2 = self.subchannel - elif self.instrument == 'NIRSPEC': + elif self.instrument == "NIRSPEC": self.grating = median_model.meta.instrument.grating self.filter = median_model.meta.instrument.filter self.par_median_select1 = self.grating @@ -65,15 +60,16 @@ def __init__(self, median_model, input_models): # set up x,y,z of Median Cube # Median cube should have linear wavelength xcube, ycube, zcube = wcstools.grid_from_bounding_box( - self.median_skycube.meta.wcs.bounding_box, - step=(1, 1, 1)) + self.median_skycube.meta.wcs.bounding_box, step=(1, 1, 1) + ) - # using wcs of ifu cube determine ra,dec,lambda - self.cube_ra, self.cube_dec, self.cube_wave = \ - self.median_skycube.meta.wcs(xcube + 1, ycube + 1, zcube + 1) + # using wcs of ifu cube determine ra, dec, lambda + self.cube_ra, self.cube_dec, self.cube_wave = self.median_skycube.meta.wcs( + xcube + 1, ycube + 1, zcube + 1 + ) # pull out flux from the median sky cube that matches with - # cube_ra,dec,wave + # cube_ra, dec, wave self.cube_flux = self.median_skycube.data # remove all the nan values - just in case @@ -93,11 +89,11 @@ def __init__(self, median_model, input_models): self.input_list_number = [] for icount, model in enumerate(input_models): - if self.instrument == 'MIRI': + if self.instrument == "MIRI": par1 = model.meta.instrument.channel par2 = model.meta.instrument.band.lower() found2 = par2.find(self.par_median_select2) - if self.instrument == 'NIRSPEC': + if self.instrument == "NIRSPEC": par1 = model.meta.instrument.grating par2 = model.meta.instrument.grating found2 = 1 @@ -106,43 +102,57 @@ def __init__(self, median_model, input_models): if found1 > -1 and found2 > -1: self.input_models.append(model) self.input_list_number.append(icount) + # ********************************************************************** def blot_info(self): - """ Prints the basic parameters of the blot image and median sky cube - """ - log.info('Information on Blotting') - log.info('Working with instrument %s', self.instrument) - log.info('Shape of sky cube %f %f %f', - self.median_skycube.data.shape[2], - self.median_skycube.data.shape[1], - self.median_skycube.data.shape[0]) - - if self.instrument == 'MIRI': - log.info('Channel %s', self.channel) - log.info('Sub-channel %s', self.subchannel) - - elif self.instrument == 'NIRSPEC': - log.info('Grating %s', self.grating) - log.info('Filter %s', self.filter) + """Print the basic parameters of the blot image and median sky cube.""" + log.info("Information on Blotting") + log.info("Working with instrument %s", self.instrument) + log.info( + "Shape of sky cube %f %f %f", + self.median_skycube.data.shape[2], + self.median_skycube.data.shape[1], + self.median_skycube.data.shape[0], + ) + + if self.instrument == "MIRI": + log.info("Channel %s", self.channel) + log.info("Sub-channel %s", self.subchannel) + + elif self.instrument == "NIRSPEC": + log.info("Grating %s", self.grating) + log.info("Filter %s", self.filter) # *********************************************************************** def blot_images(self): - if self.instrument == 'MIRI': + """ + Call the instrument specific blotting code. + + Returns + ------- + blotmodels : ModelContainer + Blotted IFU image models + input_list_number : list of int + List containing index of blot model in input_models + """ + if self.instrument == "MIRI": blotmodels = self.blot_images_miri() - elif self.instrument == 'NIRSPEC': + elif self.instrument == "NIRSPEC": blotmodels = self.blot_images_nirspec() return blotmodels, self.input_list_number + # ************************************************************************ def blot_images_miri(self): - """ Core blotting routine for MIRI + """ + Core blotting routine for MIRI. This is the main routine for blotting the MIRI median sky cube back to the detector space and creating a blotting image for each input model - 1. Loop over every data model to be blotted and find ra,dec,wavelength - for every pixel in a valid slice on the detector. + 1. Loop over every data model to be blotted and find ra, dec and + wavelength for every pixel in a valid slice on the detector. 2. Loop over every input model and using the inverse (backwards) transform convert the median sky cube values ra, dec, lambda to the blotted x, y detector value (x_cube, y_cube). @@ -150,6 +160,11 @@ def blot_images_miri(self): the x, y detector values that fall within the roi region. The blotted flux = the weighted flux, where the weight is based on distance between the center of the blotted pixel and the detector pixel. + + Returns + ------- + blot_models : ModelContainer + Container of blotted IFUImage models """ blot_models = ModelContainer() instrument_info = instrument_defaults.InstrumentInfo() @@ -169,8 +184,7 @@ def blot_images_miri(self): ydet, xdet = np.mgrid[:ysize, :xsize] ydet = ydet.flatten() xdet = xdet.flatten() - self.ycenter_grid, self.xcenter_grid = np.mgrid[0:ysize, - 0:xsize] + self.ycenter_grid, self.xcenter_grid = np.mgrid[0:ysize, 0:xsize] xsize2 = xend - xstart + 1 xcenter = np.arange(xsize2) + xstart @@ -180,23 +194,23 @@ def blot_images_miri(self): xdet = xdet[valid_channel] ydet = ydet[valid_channel] # cube spaxel ra,dec values --> x, y on detector - x_cube, y_cube = model.meta.wcs.backward_transform(self.cube_ra, - self.cube_dec, - self.cube_wave) + x_cube, y_cube = model.meta.wcs.backward_transform( + self.cube_ra, self.cube_dec, self.cube_wave + ) x_cube = np.ndarray.flatten(x_cube) y_cube = np.ndarray.flatten(y_cube) flux_cube = np.ndarray.flatten(self.cube_flux) valid = ~np.isnan(y_cube) valid_channel = np.logical_and(x_cube >= xstart, x_cube <= xend) - valid_flux = (flux_cube != 0) + valid_flux = flux_cube != 0 fuse = np.where(valid & valid_channel & valid_flux) x_cube = x_cube[fuse] y_cube = y_cube[fuse] flux_cube = flux_cube[fuse] - log.info('Blotting back to %s', model.meta.filename) + log.info("Blotting back to %s", model.meta.filename) # ______________________________________________________________________________ # blot_wrapper is a c extension that finds: # the overlapping median cube spaxels with the detector pixels @@ -207,9 +221,9 @@ def blot_images_miri(self): roi_det = 1.0 # Just large enough that we don't get holes # set up c wrapper for blotting - result = blot_wrapper(roi_det, xsize, ysize, xstart, xsize2, - xcenter, ycenter, - x_cube, y_cube, flux_cube) + result = blot_wrapper( + roi_det, xsize, ysize, xstart, xsize2, xcenter, ycenter, x_cube, y_cube, flux_cube + ) blot_flux, blot_weight = result igood = np.where(blot_weight > 0) blot_flux[igood] = blot_flux[igood] / blot_weight[igood] @@ -222,7 +236,8 @@ def blot_images_miri(self): # ************************************************************************ def blot_images_nirspec(self): - """ Core blotting routine for NIRSPEC + """ + Core blotting routine for NIRSPEC. This is the main routine for blotting the NIRSPEC median sky cube back to the detector space and creating a blotting image for each input model. @@ -243,6 +258,10 @@ def blot_images_nirspec(self): determines the overlap of the blotted x,y values with a regular grid setup in the detector plane which is the blotted image. + Returns + ------- + blot_models : ModelContainer + Container of blotted IFUImage models """ blot_models = ModelContainer() @@ -261,18 +280,19 @@ def blot_images_nirspec(self): # for NIRSPEC wcs information accessed separately for each slice nslices = 30 - log.info('Blotting 30 slices on NIRSPEC detector') + log.info("Blotting 30 slices on NIRSPEC detector") roi_det = 1.0 # Just large enough that we don't get holes - wcsobj, tr1, tr2, tr3 = nirspec._get_transforms(model, np.arange(nslices)) + wcsobj, tr1, tr2, tr3 = nirspec._get_transforms(model, np.arange(nslices)) # noqa: SLF001 for ii in range(nslices): # for each slice pull out the blotted values that actually fall on the slice region # use the bounding box of each slice to determine the slice limits - slice_wcs = nirspec._nrs_wcs_set_input_lite(model, wcsobj, ii, - [tr1, tr2[ii], tr3[ii]]) - slicer2world = slice_wcs.get_transform('slicer','world') - detector2slicer = slice_wcs.get_transform('detector','slicer') + slice_wcs = nirspec._nrs_wcs_set_input_lite( # noqa: SLF001 + model, wcsobj, ii, [tr1, tr2[ii], tr3[ii]] + ) + slicer2world = slice_wcs.get_transform("slicer", "world") + detector2slicer = slice_wcs.get_transform("detector", "slicer") # find some rough limits on ra,dec, lambda using the x,y -> ra,dec,lambda x, y = wcstools.grid_from_bounding_box(slice_wcs.bounding_box) @@ -297,7 +317,7 @@ def blot_images_nirspec(self): use1 = np.logical_and(self.cube_ra >= ramin, self.cube_ra <= ramax) use2 = np.logical_and(self.cube_dec >= decmin, self.cube_dec <= decmax) use3 = np.logical_and(self.cube_wave >= lam_min, self.cube_wave <= lam_max) - use = np.logical_and(np.logical_and(use1, use2),use3) + use = np.logical_and(np.logical_and(use1, use2), use3) ra_use = self.cube_ra[use] dec_use = self.cube_dec[use] @@ -305,10 +325,11 @@ def blot_images_nirspec(self): flux_use = self.cube_flux[use] # get the indices of elements on the slice - onslice_ind = in_ifu_slice(slice_wcs,ra_use,dec_use,wave_use) + onslice_ind = in_ifu_slice(slice_wcs, ra_use, dec_use, wave_use) slx, sly, sllam = slicer2world.inverse(ra_use, dec_use, wave_use) - xslice,yslice = detector2slicer.inverse(slx[onslice_ind], sly[onslice_ind], - sllam[onslice_ind]) + xslice, yslice = detector2slicer.inverse( + slx[onslice_ind], sly[onslice_ind], sllam[onslice_ind] + ) # pull out region for slice fluxslice = flux_use[onslice_ind] @@ -317,7 +338,7 @@ def blot_images_nirspec(self): xlimit, ylimit = slice_wcs.bounding_box xuse = np.logical_and(xslice >= xlimit[0], xslice <= xlimit[1]) yuse = np.logical_and(yslice >= ylimit[0], yslice <= ylimit[1]) - use = np.logical_and(xuse,yuse) + use = np.logical_and(xuse, yuse) xuse = xslice[use] yuse = yslice[use] flux_use = fluxslice[use] @@ -336,9 +357,18 @@ def blot_images_nirspec(self): xstart = 0 xsize2 = blot_xsize - result = blot_wrapper(roi_det, blot_xsize, blot_ysize, xstart, xsize2, - xcenter, ycenter, - x_total, y_total, flux_total) + result = blot_wrapper( + roi_det, + blot_xsize, + blot_ysize, + xstart, + xsize2, + xcenter, + ycenter, + x_total, + y_total, + flux_total, + ) blot_flux_slice, blot_weight_slice = result blot_flux = blot_flux + blot_flux_slice blot_weight = blot_weight + blot_weight_slice diff --git a/jwst/cube_build/coord.py b/jwst/cube_build/coord.py index c1fa35d617..a12a9a0a2f 100644 --- a/jwst/cube_build/coord.py +++ b/jwst/cube_build/coord.py @@ -1,34 +1,34 @@ -""" A set of routines to assist in the WCS transforms used in the -cube_build step -""" +"""A set of routines to assist in the WCS transforms used in the cube_build step.""" + import numpy as np import math # _______________________________________________________________________ def radec2std(crval1, crval2, ra, dec, rot_angle=None): - """ Compute the tangent projection coordinates (xi,eta) from ra,dec - using crval1 and crval2 (the tangent point). + """ + Compute the tangent projection coordinates (xi,eta) from ra,dec using crval1 and crval2. Parameters - ____________ + ---------- crval1 : float RA value of tangent point crval2 : float - DEC value of tangent point - rot_angle: float - rotation angle given in degrees + Dec value of tangent point ra : numpy.ndarray or float - A list (or single value) of ra points to convert + A list (or single value) of RA points to convert dec : numpy.ndarray or float - A list (or single value) of ra points to convert - - Return Values - _____________ - xi, eta - rectangular coordinates of tangent plane projected ra,dec + A list (or single value) of Dec points to convert + rot_angle : float or None, optional + Rotation angle given in degrees + Returns + ------- + xi : float + X-axis tangent plane coordinate of ra,dec + eta : float + Y-axis tangent plane coordinate of ra,dec """ - if np.isscalar(ra): ra = np.asarray([ra]) dec = np.asarray([dec]) @@ -44,8 +44,7 @@ def radec2std(crval1, crval2, ra, dec, rot_angle=None): h = np.sin(decr) * math.sin(dec0) + np.cos(decr) * math.cos(dec0) * np.cos(radiff) xi = np.cos(decr) * np.sin(radiff) / h - eta = (np.sin(decr) * math.cos(dec0) - - np.cos(decr) * math.sin(dec0) * np.cos(radiff)) / h + eta = (np.sin(decr) * math.cos(dec0) - np.cos(decr) * math.sin(dec0) * np.cos(radiff)) / h xi = -xi xi = xi * rad2arcsec @@ -61,36 +60,38 @@ def radec2std(crval1, crval2, ra, dec, rot_angle=None): eta = temp2 return xi, eta + + # ________________________________________________________________________________ def std2radec(crval1, crval2, xi, eta): - """ Compute ra,dec from the tangent plane rectangular coordinates + """ + Compute ra,dec from the tangent plane rectangular coordinates. Compute the ra,dec values of tangent plane rectangular coordinates using crval1, crval2(the tangent point). This routine takes the rectangular - plane and projects it to the spherical plane using crval1, crval2 as + plane and projects it onto the spherical plane using crval1, crval2 as the tangent plane. Parameters - ____________ + ---------- crval1 : float RA value of tangent point crval2 : float - DEC value of tangent point + Dec value of tangent point xi : float - xi rectangular coordinates of tangent plane projected ra,dec - eta : float - eta rectangular coordinates of tangent plane projected ra,dec + Xi rectangular coordinate of tangent plane projected ra,dec + eta : float + Eta rectangular coordinate of tangent plane projected ra,dec - Return Values - _____________ + Returns + ------- ra : float - list (or single value) of ra computed values + List (or single value) of ra computed values dec : float - list (or single value) of dec computed values + List (or single value) of dec computed values """ - if np.isscalar(xi): eta = np.asarray([eta]) xi = np.asarray([xi]) @@ -128,11 +129,13 @@ def std2radec(crval1, crval2, xi, eta): ra[mask] -= 360.0 return ra, dec + # _______________________________________________________________________ def v2v32radec_estimate(ra_ref, dec_ref, roll_ref, v2_ref, v3_ref, v2, v3): - """ Estimation of ra and dec from the v2, v3 coordinates + """ + Estimation of ra and dec from the v2, v3 coordinates. This routine is used for debugging purposes. It is not actually used in the cube_build step for routine IFU cube building. @@ -142,33 +145,35 @@ def v2v32radec_estimate(ra_ref, dec_ref, roll_ref, v2_ref, v3_ref, v2, v3): Parameters ---------- ra_ref : float - ra of reference point given in arc seconds + RA of reference point given in arc seconds dec_ref : float - dec of reference point given in arc seconds + Dec of reference point given in arc seconds roll_ref : float - roll angle given in degrees + Roll angle given in degrees v2_ref : float - v2 coordinate of reference point given in arc seconds + V2 coordinate of reference point given in arc seconds v3_ref : float - v3 coordinate of reference point given in arc seconds + V3 coordinate of reference point given in arc seconds v2 : float - v2 coordinate given in arc seconds + V2 coordinate given in arc seconds v3 : float - v3 coordinate given in arc seconds + V3 coordinate given in arc seconds Returns ------- ra : float + Calculated ra from v2, v3 dec : float + Calculate dec from v2, v3 Notes - ---- - it is assumed that the v2,v3 coordinates have the effects of dithering included + ----- + It is assumed that the v2,v3 coordinates have the effects of dithering included. """ d2r = math.pi / 180.0 - v2deg = v2.copy() / 3600.0 # convert to degrees - v3deg = v3.copy() / 3600.0 # convert to degrees + v2deg = v2.copy() / 3600.0 # convert to degrees + v3deg = v3.copy() / 3600.0 # convert to degrees v2_ref = v2_ref / 3600.0 # convert to degrees v3_ref = v3_ref / 3600.0 # convert to degrees @@ -176,7 +181,7 @@ def v2v32radec_estimate(ra_ref, dec_ref, roll_ref, v2_ref, v3_ref, v2, v3): roll_ref_rad = roll_ref * d2r delta_v2 = (v2deg - v2_ref) * math.cos(v3_ref_rad) - delta_v3 = (v3deg - v3_ref) + delta_v3 = v3deg - v3_ref delta_ra = delta_v2 * math.cos(roll_ref_rad) + delta_v3 * math.sin(roll_ref_rad) delta_dec = -delta_v2 * math.sin(roll_ref_rad) + delta_v3 * math.cos(roll_ref_rad) @@ -184,42 +189,44 @@ def v2v32radec_estimate(ra_ref, dec_ref, roll_ref, v2_ref, v3_ref, v2, v3): dec = delta_dec + dec_ref return ra, dec + + # _______________________________________________________________________ def radec2v2v3_estimate(ra_ref, dec_ref, roll_ref, v2_ref, v3_ref, ra, dec): - """ Convert ra,dec to v2, v3 + """ + Convert ra,dec to v2, v3. This routine is used for debugging purposes. It is not actually used in the cube_build step for routine IFU cube building. - The conversion from Ra,Dec to V2,V3 is handled more accurately by + The conversion from RA,Dec to V2,V3 is handled more accurately by the transforms provided by assign_wcs. Parameters ---------- ra_ref : float - ra of reference point given in degrees + RA of reference point given in degrees dec_ref : float - dec of reference point given in degrees + Dec of reference point given in degrees roll_ref : float - roll angle given in degrees + Roll angle given in degrees v2_ref : float - v2 coordinate of reference point given in arc seconds + V2 coordinate of reference point given in arc seconds v3_ref : float - v3 coordinate of reference point given in arc seconds + V3 coordinate of reference point given in arc seconds ra : float - ra coordinate given in degrees + RA coordinate given in degrees dec : float - dec coordinate given in degrees + Dec coordinate given in degrees Returns ------- v2 : float - v2 coordinate in arc seconds + V2 coordinate in arc seconds v3 : float - v2 coordinate in arc seconds - """ - + V3 coordinate in arc seconds + """ d2r = math.pi / 180.0 r2d = 180.0 / math.pi @@ -234,7 +241,7 @@ def radec2v2v3_estimate(ra_ref, dec_ref, roll_ref, v2_ref, v3_ref, ra, dec): this_dec = dec * d2r delta_ra = (this_ra - ra_ref_rad) * math.cos(dec_ref_rad) - delta_dec = (this_dec - dec_ref_rad) + delta_dec = this_dec - dec_ref_rad dv2 = delta_ra * math.cos(roll_ref_rad) - delta_dec * math.sin(roll_ref_rad) dv3 = delta_ra * math.sin(roll_ref_rad) + delta_dec * math.cos(roll_ref_rad) diff --git a/jwst/cube_build/cube_build.py b/jwst/cube_build/cube_build.py index 845bf8888a..716f7b7b2f 100644 --- a/jwst/cube_build/cube_build.py +++ b/jwst/cube_build/cube_build.py @@ -1,5 +1,5 @@ -""" Basic routines used to set up IFU cubes -""" +"""Basic routines used to set up IFU cubes.""" + import logging from . import cube_build_io_util from . import file_table @@ -9,16 +9,12 @@ log.setLevel(logging.DEBUG) -class CubeData(): - """ class CubeData holds top level information on the ifucube - - """ +class CubeData: + """Hold top level information on the ifucube.""" - def __init__(self, - input_models, - par_filename, - **pars): - """ Initialize the high level of information for the ifu cube + def __init__(self, input_models, par_filename, **pars): + """ + Initialize high level of information for the ifu cube. The class CubeData holds information on what type of cube is being created and how the cube is to be constructed. This information includes: @@ -29,21 +25,22 @@ def __init__(self, Parameters ---------- - input_models : list of data models - par_filename: str - cube parameter reference filename - pars : dictionary holding top level cube parameters + input_models : list + List of data models + par_filename : str + Cube parameter reference filename + **pars : dict + Holds top level cube parameters """ - self.input_models = input_models self.par_filename = par_filename - self.single = pars.get('single') - self.channel = pars.get('channel') - self.subchannel = pars.get('subchannel') - self.grating = pars.get('grating') - self.filter = pars.get('filter') - self.weighting = pars.get('weighting') - self.output_type = pars.get('output_type') + self.single = pars.get("single") + self.channel = pars.get("channel") + self.subchannel = pars.get("subchannel") + self.grating = pars.get("grating") + self.filter = pars.get("filter") + self.weighting = pars.get("weighting") + self.output_type = pars.get("output_type") self.instrument = None self.all_channel = [] @@ -51,77 +48,84 @@ def __init__(self, self.all_grating = [] self.all_filter = [] - self.output_name = '' -# _____________________________________________________________________________ + self.output_name = "" + + # _____________________________________________________________________________ def setup(self): - """ Set up IFU Cube. Determine band coverage and read in reference files + """ + Set up IFU Cube, determine band coverage, and read in reference files. Read in the input_models and fill in the dictionary, master_table, that stores the data for each channel/subchannel or grating/filter. - If the channel/subchannel or grating/filter are not set by the user, then determine which ones are found in the data - Read in necessary reference data: - * cube parameter reference file - This routine fills in the instrument_info dictionary, which holds the default spatial and spectral size of the output cube, as well as, the region of influence size in the spatial and spectral dimensions. The dictionary master_table is also filled in. This dictionary contains list of datamodels for each band (channel/subchannel or grating/filter combination). + + Returns + ------- + instrument : str + Instrument name, either "MIRI" or "NIRSpec" + instrument_info : dict + Dictionary containing information on bands + master_table : dict + Dictionary containing files organized by type of data: + instrument, channel/band or grating/filter """ -# _____________________________________________________________________________ -# Read in the input data (association table or single file) -# Fill in MasterTable based on Channel/Subchannel or filter/grating -# ______________________________________________________________________________ + # Read in the input data (association table or single file) + # Fill in MasterTable based on Channel/Subchannel or filter/grating + master_table = file_table.FileTable() instrument = master_table.set_file_table(self.input_models) -# _______________________________________________________________________________ -# find out how many files are in the association table or if it is an single -# file store the input_filenames and input_models + # find out how many files are in the association table or if it is an single + # file store the input_filenames and input_models self.instrument = instrument -# _______________________________________________________________________________ -# Determine which channels/subchannels or filter/grating will be covered by the -# spectral cube. -# fills in band_channel, band_subchannel, band_grating, band_filer -# _______________________________________________________________________________ + + # Determine which channels/subchannels or filter/grating will be covered by the + # spectral cube. + # fills in band_channel, band_subchannel, band_grating, band_filter self.determine_band_coverage(master_table) -# _______________________________________________________________________________ -# instrument_defaults is an dictionary class that holds default parameters for -# each band for the different instruments -# _______________________________________________________________________________ + + # instrument_defaults is a dictionary class that holds default parameters for + # each band for the different instruments instrument_info = instrument_defaults.InstrumentInfo() -# ------------------------------------------------------------------------------- -# Read the cube pars reference file - log.info('Reading cube parameter file %s', self.par_filename) - cube_build_io_util.read_cubepars(self.par_filename, - self.instrument, - self.weighting, - self.all_channel, - self.all_subchannel, - self.all_grating, - self.all_filter, - instrument_info) -# ------------------------------------------------------------------------------- + + # Read the cube pars reference file + log.info("Reading cube parameter file %s", self.par_filename) + cube_build_io_util.read_cubepars( + self.par_filename, + self.instrument, + self.weighting, + self.all_channel, + self.all_subchannel, + self.all_grating, + self.all_filter, + instrument_info, + ) + self.instrument_info = instrument_info -# _______________________________________________________________________________ -# Set up values to return and access for other parts of cube_build + # Set up values to return and access for other parts of cube_build self.master_table = master_table - return {'instrument': self.instrument, - 'instrument_info': self.instrument_info, - 'master_table': self.master_table} + return { + "instrument": self.instrument, + "instrument_info": self.instrument_info, + "master_table": self.master_table, + } -# ******************************************************************************** + # ******************************************************************************** def determine_band_coverage(self, master_table): - """Either by user parameters or the data itself, determine which bands - cover by cube + """ + Determine which bands are covered by the cube. To determine which bands are covered by the output cube: 1. If MIRI data then check if the user has set either the channel or @@ -138,30 +142,39 @@ def determine_band_coverage(self, master_table): Parameters ---------- - master_table: dictionary + master_table : dict A dictionary for each band that contains of list of datamodels covering that band. Raises ------ - ErrorNoChannels + NoChannelsError The user selected channels are not in the data - ErrorNoSubChannels + NoSubChannelsError The user selected subchannels are not in the data - ErrorNoGratings + NoGratingsError The user selected gratings are not in the data - ErrorNoFilters + NoFiltersError The user selected filters are not in the data - ErrorMissingParameter + MissingParameterError The user selected grating but not filter or vice versa """ -# ________________________________________________________________________________ -# IF INSTRUMENT = MIRI -# loop over the file names - if self.instrument == 'MIRI': - valid_channel = ['1', '2', '3', '4'] - valid_subchannel = ['short', 'medium', 'long', 'short-medium', 'short-long', - 'medium-short', 'medium-long', 'long-short', 'long-medium'] + # ________________________________________________________________________________ + # IF INSTRUMENT = MIRI + # loop over the file names + if self.instrument == "MIRI": + valid_channel = ["1", "2", "3", "4"] + valid_subchannel = [ + "short", + "medium", + "long", + "short-medium", + "short-long", + "medium-short", + "medium-long", + "long-short", + "long-medium", + ] nchannels = len(valid_channel) nsubchannels = len(valid_subchannel) @@ -174,18 +187,18 @@ def determine_band_coverage(self, master_table): user_clen = len(self.channel) user_slen = len(self.subchannel) - if user_clen == 1 and self.channel[0] == 'all': + if user_clen == 1 and self.channel[0] == "all": user_clen = 0 - if user_slen == 1 and self.subchannel[0] == 'all': + if user_slen == 1 and self.subchannel[0] == "all": user_slen = 0 - # _______________________________________________________________________________ for i in range(nchannels): for j in range(nsubchannels): - nfiles = len(master_table.FileMap['MIRI'][valid_channel[i]][valid_subchannel[j]]) + nfiles = len( + master_table.FileMap["MIRI"][valid_channel[i]][valid_subchannel[j]] + ) if nfiles > 0: - # ______________________________________________________ # neither parameters are set if user_clen == 0 and user_slen == 0: self.all_channel.append(valid_channel[i]) @@ -202,68 +215,99 @@ def determine_band_coverage(self, master_table): self.all_subchannel.append(valid_subchannel[j]) # both parameters set else: - if (valid_channel[i] in self.channel and - valid_subchannel[j] in self.subchannel): + if ( + valid_channel[i] in self.channel + and valid_subchannel[j] in self.subchannel + ): self.all_channel.append(valid_channel[i]) self.all_subchannel.append(valid_subchannel[j]) - log.info('The desired cubes cover the MIRI Channels: %s', - self.all_channel) - log.info('The desired cubes cover the MIRI subchannels: %s', - self.all_subchannel) + log.info("The desired cubes cover the MIRI Channels: %s", self.all_channel) + log.info("The desired cubes cover the MIRI subchannels: %s", self.all_subchannel) number_channels = len(self.all_channel) number_subchannels = len(self.all_subchannel) if number_channels == 0: - raise ErrorNoChannels( - "The cube does not cover any channels, change channel parameter") + raise NoChannelsError( + "The cube does not cover any channels, change channel parameter" + ) if number_subchannels == 0: - raise ErrorNoSubchannels( - "The cube does not cover any subchannels, change band parameter") + raise NoSubchannelsError( + "The cube does not cover any subchannels, change band parameter" + ) # _______________________________________________________________ - if self.instrument == 'NIRSPEC': + if self.instrument == "NIRSPEC": # 1 to 1 mapping valid_gwa[i] -> valid_fwa[i] - valid_gwa = ['g140m', 'g140h', 'g140m', 'g140h', 'g235m', - 'g235h', 'g395m', 'g395h', 'prism', - 'prism', 'g140m', 'g140h', 'g235m', 'g235h', 'g395m', 'g395h'] - valid_fwa = ['f070lp', 'f070lp', 'f100lp', 'f100lp', 'f170lp', - 'f170lp', 'f290lp', 'f290lp', 'clear', - 'opaque', 'opaque', 'opaque', 'opaque', 'opaque', 'opaque', 'opaque'] + valid_gwa = [ + "g140m", + "g140h", + "g140m", + "g140h", + "g235m", + "g235h", + "g395m", + "g395h", + "prism", + "prism", + "g140m", + "g140h", + "g235m", + "g235h", + "g395m", + "g395h", + ] + valid_fwa = [ + "f070lp", + "f070lp", + "f100lp", + "f100lp", + "f170lp", + "f170lp", + "f290lp", + "f290lp", + "clear", + "opaque", + "opaque", + "opaque", + "opaque", + "opaque", + "opaque", + "opaque", + ] nbands = len(valid_fwa) user_glen = len(self.grating) user_flen = len(self.filter) - if user_glen == 1 and self.grating[0] == 'all': + if user_glen == 1 and self.grating[0] == "all": user_glen = 0 - if user_flen == 1 and self.filter[0] == 'all': + if user_flen == 1 and self.filter[0] == "all": user_flen = 0 # check if input filter or grating has been set if user_glen == 0 and user_flen != 0: - raise ErrorMissingParameter("Filter specified, but Grating was not") + raise MissingParameterError("Filter specified, but Grating was not") + elif user_glen != 0 and user_flen == 0: + raise MissingParameterError("Grating specified, but Filter was not") for i in range(nbands): - nfiles = len(master_table.FileMap['NIRSPEC'][valid_gwa[i]][valid_fwa[i]]) + nfiles = len(master_table.FileMap["NIRSPEC"][valid_gwa[i]][valid_fwa[i]]) if nfiles > 0: - # _________________________________________________ # neither parameters are set if user_glen == 0 and user_flen == 0: self.all_grating.append(valid_gwa[i]) self.all_filter.append(valid_fwa[i]) - # __________________________________________________ + # grating was set by user but not filter elif user_glen != 0 and user_flen == 0: if valid_gwa[i] in self.grating: self.all_grating.append(valid_gwa[i]) self.all_filter.append(valid_fwa[i]) - # __________________________________________________ # both parameters set else: - if (valid_fwa[i] in self.filter and - valid_gwa[i] in self.grating): + if valid_fwa[i] in self.filter and valid_gwa[i] in self.grating: self.all_grating.append(valid_gwa[i]) self.all_filter.append(valid_fwa[i]) @@ -271,63 +315,84 @@ def determine_band_coverage(self, master_table): number_gratings = len(self.all_grating) if number_filters == 0: - raise ErrorNoFilters("The cube does not cover any filters") + raise NoFiltersError("The cube does not cover any filters") if number_gratings == 0: - raise ErrorNoGratings("The cube does not cover any gratings") -# ______________________________________________________________________ + raise NoGratingsError("The cube does not cover any gratings") + + # ______________________________________________________________________ def number_cubes(self): - """Determine the number of IFUcubes to created based on: - Type of cube (single band, multiple bands, or Single mode) + """ + Determine the number of IFUcubes to create. + + The number of cubes depends if the cube is created from a single + band or multiple bands. In calspec2, 1 cube is created. If the data + is MIRI MRS then 2 channels worth of data is combined into one IFUCube. + For calspec3 data, the default is to create multiple cubes containing + single band data. However, the user can override this option and combine + multiple bands in a single IFUcube. + + Returns + ------- + num_cubes : int + Number of IFU cubes being created + cube_pars : list + For each IFUcube, cube_pars holds 2 parameters. For MIRI data the + parameters hold channel type and band type, while for NIRSpec data + the pamaters hold grating type and filter type. """ num_cubes = 0 cube_pars = {} -# ______________________________________________________________________ -# MIRI -# ______________________________________________________________________ - if self.instrument == 'MIRI': + # ______________________________________________________________________ + # MIRI + # ______________________________________________________________________ + if self.instrument == "MIRI": band_channel = self.all_channel band_subchannel = self.all_subchannel -# user, single, or multi - if (self.output_type == 'user' or self.output_type == 'single' or - self.output_type == 'multi'): - - if self.output_type == 'multi': - log.info('Output IFUcube are constructed from all the data ') + # user, single, or multi + if ( + self.output_type == "user" + or self.output_type == "single" + or self.output_type == "multi" + ): + if self.output_type == "multi": + log.info("Output IFUcube are constructed from all the data ") if self.single: - log.info('Single = true, creating a set of single exposures mapped' + - ' to output IFUCube coordinate system') - if self.output_type == 'user': - log.info('The user has selected the type of IFU cube to make') + log.info( + "Single = true, creating a set of single exposures mapped" + " to output IFUCube coordinate system" + ) + if self.output_type == "user": + log.info("The user has selected the type of IFU cube to make") num_cubes = 1 - cube_pars['1'] = {} - cube_pars['1']['par1'] = {} - cube_pars['1']['par2'] = {} - cube_pars['1']['par1'] = self.all_channel - cube_pars['1']['par2'] = self.all_subchannel + cube_pars["1"] = {} + cube_pars["1"]["par1"] = {} + cube_pars["1"]["par2"] = {} + cube_pars["1"]["par1"] = self.all_channel + cube_pars["1"]["par2"] = self.all_subchannel -# default band cubes - if self.output_type == 'band': - log.info('Output Cubes are single channel, single sub-channel IFU Cubes') + # default band cubes + if self.output_type == "band": + log.info("Output Cubes are single channel, single sub-channel IFU Cubes") for i in range(len(band_channel)): num_cubes = num_cubes + 1 cube_no = str(num_cubes) cube_pars[cube_no] = {} - cube_pars[cube_no]['pars1'] = {} - cube_pars[cube_no]['pars2'] = {} + cube_pars[cube_no]["pars1"] = {} + cube_pars[cube_no]["pars2"] = {} this_channel = [] this_subchannel = [] this_channel.append(band_channel[i]) this_subchannel.append(band_subchannel[i]) - cube_pars[cube_no]['par1'] = this_channel - cube_pars[cube_no]['par2'] = this_subchannel + cube_pars[cube_no]["par1"] = this_channel + cube_pars[cube_no]["par2"] = this_subchannel -# default channel cubes - if self.output_type == 'channel': - log.info('Output cubes are single channel and all subchannels in data') + # default channel cubes + if self.output_type == "channel": + log.info("Output cubes are single channel and all subchannels in data") num_cubes = 0 channel_no_repeat = list(set(band_channel)) @@ -335,104 +400,107 @@ def number_cubes(self): num_cubes = num_cubes + 1 cube_no = str(num_cubes) cube_pars[cube_no] = {} - cube_pars[cube_no]['pars1'] = {} - cube_pars[cube_no]['pars2'] = {} + cube_pars[cube_no]["pars1"] = {} + cube_pars[cube_no]["pars2"] = {} this_channel = [] this_subchannel = [] for k, j in enumerate(band_channel): if j == i: this_subchannel.append(band_subchannel[k]) this_channel.append(i) - cube_pars[cube_no]['par1'] = this_channel - cube_pars[cube_no]['par2'] = this_subchannel - -# ______________________________________________________________________ -# NIRSPEC -# ______________________________________________________________________ - if self.instrument == 'NIRSPEC': + cube_pars[cube_no]["par1"] = this_channel + cube_pars[cube_no]["par2"] = this_subchannel + # ______________________________________________________________________ + # NIRSPEC + # ______________________________________________________________________ + if self.instrument == "NIRSPEC": band_grating = self.all_grating band_filter = self.all_filter - if (self.output_type == 'user' or self.output_type == 'single' or - self.output_type == 'multi'): - if self.output_type == 'multi': - log.info('Output IFUcube are constructed from all the data ') + if ( + self.output_type == "user" + or self.output_type == "single" + or self.output_type == "multi" + ): + if self.output_type == "multi": + log.info("Output IFUcube are constructed from all the data ") if self.single: - log.info('Single = true, creating a set of single exposures' + - ' mapped to output IFUCube coordinate system') - if self.output_type == 'user': - log.info('The user has selected the type of IFU cube to make') + log.info( + "Single = true, creating a set of single exposures" + " mapped to output IFUCube coordinate system" + ) + if self.output_type == "user": + log.info("The user has selected the type of IFU cube to make") num_cubes = 1 - cube_pars['1'] = {} - cube_pars['1']['par1'] = {} - cube_pars['1']['par2'] = {} - cube_pars['1']['par1'] = self.all_grating - cube_pars['1']['par2'] = self.all_filter - -# default band cubes - if self.output_type == 'band': - log.info('Output Cubes are single grating, single filter IFU Cubes') + cube_pars["1"] = {} + cube_pars["1"]["par1"] = {} + cube_pars["1"]["par2"] = {} + cube_pars["1"]["par1"] = self.all_grating + cube_pars["1"]["par2"] = self.all_filter + + # default band cubes + if self.output_type == "band": + log.info("Output Cubes are single grating, single filter IFU Cubes") for i in range(len(band_grating)): num_cubes = num_cubes + 1 cube_no = str(num_cubes) cube_pars[cube_no] = {} - cube_pars[cube_no]['pars1'] = {} - cube_pars[cube_no]['pars2'] = {} + cube_pars[cube_no]["pars1"] = {} + cube_pars[cube_no]["pars2"] = {} this_grating = [] this_filter = [] this_grating.append(band_grating[i]) this_filter.append(band_filter[i]) - cube_pars[cube_no]['par1'] = this_grating - cube_pars[cube_no]['par2'] = this_filter -# default grating cubes - if self.output_type == 'grating': - log.info('Output cubes are single grating & all filters in data') + cube_pars[cube_no]["par1"] = this_grating + cube_pars[cube_no]["par2"] = this_filter + # default grating cubes + if self.output_type == "grating": + log.info("Output cubes are single grating & all filters in data") num_cubes = 0 for i in range(len(band_grating)): num_cubes = num_cubes + 1 cube_no = str(num_cubes) cube_pars[cube_no] = {} - cube_pars[cube_no]['pars1'] = {} - cube_pars[cube_no]['pars2'] = {} + cube_pars[cube_no]["pars1"] = {} + cube_pars[cube_no]["pars2"] = {} this_grating = [] this_filter = band_subchannel this_grating.append(i) - cube_pars[cube_no]['par1'] = this_grating - cube_pars[cube_no]['par2'] = this_filter + cube_pars[cube_no]["par1"] = this_grating + cube_pars[cube_no]["par2"] = this_filter self.num_cubes = num_cubes self.cube_pars = cube_pars return self.num_cubes, self.cube_pars -# _____________________________________________________________________________ -class ErrorNoChannels(Exception): - """ Raises Exception if the user selected channels are not in the data - """ +class NoChannelsError(Exception): + """Raises Exception if the user-selected channels are not in the data.""" + pass -class ErrorNoSubchannels(Exception): - """ Raises Exception if the user selected subchannels are not in the data - """ +class NoSubchannelsError(Exception): + """Raises Exception if the user-selected subchannels are not in the data.""" + pass -class ErrorNoFilters(Exception): - """ Raises Exception if the user selected filters are not in the data - """ +class NoFiltersError(Exception): + """Raises Exception if the user-selected filters are not in the data.""" + pass -class ErrorNoGratings(Exception): - """ Raises Exception if the user selected gratings are not in the data - """ +class NoGratingsError(Exception): + """Raises Exception if the user-selected gratings are not in the data.""" + pass -class ErrorMissingParameter(Exception): - """ Raises Exception if provided grating but not filter or vice versa - """ +class MissingParameterError(Exception): + """Raises Exception if provided grating but not filter or vice versa.""" + pass diff --git a/jwst/cube_build/cube_build_io_util.py b/jwst/cube_build/cube_build_io_util.py index cc5bf0a03a..c7764d746a 100644 --- a/jwst/cube_build/cube_build_io_util.py +++ b/jwst/cube_build/cube_build_io_util.py @@ -1,21 +1,25 @@ -""" Read in reference files for the cube_build setp -""" +"""Read in reference files for the cube_build step.""" + from stdatamodels.jwst import datamodels import logging + log = logging.getLogger(__name__) log.setLevel(logging.DEBUG) -def read_cubepars(par_filename, - instrument, - weighting, - all_channel, - all_subchannel, - all_grating, - all_filter, - instrument_info): - """ Read in cube parameter reference file +def read_cubepars( + par_filename, + instrument, + weighting, + all_channel, + all_subchannel, + all_grating, + all_filter, + instrument_info, +): + """ + Read in cube parameter reference file. Based on the instrument and channel/subchannels (MIRI) or grating/filter(NIRSPEC), read in the appropriate columns in the @@ -25,31 +29,25 @@ def read_cubepars(par_filename, Parameters ---------- par_filename : str - cube parameter reference filename + Cube parameter reference filename instrument : str Either MIRI or NIRSPEC weighting : str Type of weighting, msm, emem or drizzle all_channel : list - all the channels contained in input data + All the channels contained in input data all_subchannel : list - all subchannels contained in input data - all_grating: list - all the gratings contained in the input data - all_filter: list - all the filters contained in the input data + All subchannels contained in input data + all_grating : list + All the gratings contained in the input data + all_filter : list + All the filters contained in the input data instrument_info : dictionary - holds the defaults spatial scales, spectral scales, roi size, + Holds the default spatial scales, spectral scales, roi size, weighting parameters, and min and max wavelengths for each for each band - - Returns - ------- - The dictionary, instrument_info, is filled in for each band covered - by the input data - """ - if instrument == 'MIRI': + if instrument == "MIRI": with datamodels.MiriIFUCubeParsModel(par_filename) as ptab: number_bands = len(all_channel) # pull out the channels and subchannels that cover the cube @@ -59,187 +57,203 @@ def read_cubepars(par_filename, this_sub = all_subchannel[i] # find the table entries for this combination for tabdata in ptab.ifucubepars_table: - table_channel = tabdata['channel'] - table_band = tabdata['band'].lower() - table_spaxelsize = tabdata['SPAXELSIZE'] - table_spectralstep = tabdata['SPECTRALSTEP'] - table_wavemin = tabdata['WAVEMIN'] - table_wavemax = tabdata['WAVEMAX'] + table_channel = tabdata["channel"] + table_band = tabdata["band"].lower() + table_spaxelsize = tabdata["SPAXELSIZE"] + table_spectralstep = tabdata["SPECTRALSTEP"] + table_wavemin = tabdata["WAVEMIN"] + table_wavemax = tabdata["WAVEMAX"] # match on this_channel and this_sub if this_channel == table_channel and this_sub == table_band: - instrument_info.SetSpatialSize(table_spaxelsize, this_channel, this_sub) - instrument_info.SetSpectralStep(table_spectralstep, this_channel, this_sub) - instrument_info.SetWaveMin(table_wavemin, this_channel, this_sub) - instrument_info.SetWaveMax(table_wavemax, this_channel, this_sub) + instrument_info.set_spatial_size(table_spaxelsize, this_channel, this_sub) + instrument_info.set_spectral_step( + table_spectralstep, this_channel, this_sub + ) + instrument_info.set_wave_min(table_wavemin, this_channel, this_sub) + instrument_info.set_wave_max(table_wavemax, this_channel, this_sub) # modified Shepard method 1/r weighting - if weighting == 'msm': + if weighting == "msm": for tabdata in ptab.ifucubepars_msm_table: - table_channel = tabdata['channel'] - table_band = tabdata['band'].lower() - table_sroi = tabdata['ROISPATIAL'] - table_wroi = tabdata['ROISPECTRAL'] - table_power = tabdata['POWER'] - table_softrad = tabdata['SOFTRAD'] - # match on this_channel and this_sub + table_channel = tabdata["channel"] + table_band = tabdata["band"].lower() + table_sroi = tabdata["ROISPATIAL"] + table_wroi = tabdata["ROISPECTRAL"] + table_power = tabdata["POWER"] + table_softrad = tabdata["SOFTRAD"] + # match on this_channel and this_sub if this_channel == table_channel and this_sub == table_band: - instrument_info.SetMSM(this_channel, this_sub, - table_sroi, table_wroi, - table_power, table_softrad) + instrument_info.set_msm( + this_channel, + this_sub, + table_sroi, + table_wroi, + table_power, + table_softrad, + ) # modified Shepard method e^-r weighting - elif weighting == 'emsm': + elif weighting == "emsm": for tabdata in ptab.ifucubepars_emsm_table: - table_channel = tabdata['channel'] - table_band = tabdata['band'].lower() - table_sroi = tabdata['ROISPATIAL'] - table_wroi = tabdata['ROISPECTRAL'] - table_scalerad = tabdata['SCALERAD'] - # match on this_channel and this_sub + table_channel = tabdata["channel"] + table_band = tabdata["band"].lower() + table_sroi = tabdata["ROISPATIAL"] + table_wroi = tabdata["ROISPECTRAL"] + table_scalerad = tabdata["SCALERAD"] + # match on this_channel and this_sub if this_channel == table_channel and this_sub == table_band: - instrument_info.SetEMSM(this_channel, this_sub, - table_sroi, table_wroi, - table_scalerad) + instrument_info.set_emsm( + this_channel, this_sub, table_sroi, table_wroi, table_scalerad + ) # read in wavelength table for multi-band data - if weighting == 'msm': + if weighting == "msm": for tabdata in ptab.ifucubepars_multichannel_msm_wavetable: - table_wave = tabdata['WAVELENGTH'] - table_sroi = tabdata['ROISPATIAL'] - table_wroi = tabdata['ROISPECTRAL'] - table_power = tabdata['POWER'] - table_softrad = tabdata['SOFTRAD'] - instrument_info.SetMultiChannelTable(table_wave, table_sroi, - table_wroi, table_power, - table_softrad) + table_wave = tabdata["WAVELENGTH"] + table_sroi = tabdata["ROISPATIAL"] + table_wroi = tabdata["ROISPECTRAL"] + table_power = tabdata["POWER"] + table_softrad = tabdata["SOFTRAD"] + instrument_info.set_multi_channel_table( + table_wave, table_sroi, table_wroi, table_power, table_softrad + ) # read in wavelength table for modified Shepard method 1/r weighting - elif weighting == 'emsm': + elif weighting == "emsm": for tabdata in ptab.ifucubepars_multichannel_emsm_wavetable: - table_wave = tabdata['WAVELENGTH'] - table_sroi = tabdata['ROISPATIAL'] - table_wroi = tabdata['ROISPECTRAL'] - table_scalerad = tabdata['SCALERAD'] - instrument_info.SetMultiChannelEMSMTable(table_wave, table_sroi, - table_wroi, table_scalerad) - elif weighting == 'drizzle': + table_wave = tabdata["WAVELENGTH"] + table_sroi = tabdata["ROISPATIAL"] + table_wroi = tabdata["ROISPECTRAL"] + table_scalerad = tabdata["SCALERAD"] + instrument_info.set_multi_channel_emsm_table( + table_wave, table_sroi, table_wroi, table_scalerad + ) + elif weighting == "drizzle": for tabdata in ptab.ifucubepars_multichannel_driz_wavetable: - table_wave = tabdata['WAVELENGTH'] - instrument_info.SetMultiChannelDrizTable(table_wave) + table_wave = tabdata["WAVELENGTH"] + instrument_info.set_multi_channel_driz_table(table_wave) ptab.close() del ptab # Read in NIRSPEC Values - elif instrument == 'NIRSPEC': + elif instrument == "NIRSPEC": with datamodels.NirspecIFUCubeParsModel(par_filename) as ptab: number_gratings = len(all_grating) for i in range(number_gratings): this_gwa = all_grating[i] this_filter = all_filter[i] for tabdata in ptab.ifucubepars_table: - table_grating = tabdata['DISPERSER'].lower() - table_filter = tabdata['FILTER'].lower() - table_spaxelsize = tabdata['SPAXELSIZE'] - table_spectralstep = tabdata['SPECTRALSTEP'] - table_wavemin = tabdata['WAVEMIN'] - table_wavemax = tabdata['WAVEMAX'] + table_grating = tabdata["DISPERSER"].lower() + table_filter = tabdata["FILTER"].lower() + table_spaxelsize = tabdata["SPAXELSIZE"] + table_spectralstep = tabdata["SPECTRALSTEP"] + table_wavemin = tabdata["WAVEMIN"] + table_wavemax = tabdata["WAVEMAX"] if this_gwa == table_grating and this_filter == table_filter: - instrument_info.SetSpatialSize(table_spaxelsize, this_gwa, this_filter) - instrument_info.SetSpectralStep(table_spectralstep, this_gwa, this_filter) - instrument_info.SetWaveMin(table_wavemin, this_gwa, this_filter) - instrument_info.SetWaveMax(table_wavemax, this_gwa, this_filter) + instrument_info.set_spatial_size(table_spaxelsize, this_gwa, this_filter) + instrument_info.set_spectral_step(table_spectralstep, this_gwa, this_filter) + instrument_info.set_wave_min(table_wavemin, this_gwa, this_filter) + instrument_info.set_wave_max(table_wavemax, this_gwa, this_filter) # modified Shepard method 1/r weighting - if weighting == 'msm': + if weighting == "msm": for tabdata in ptab.ifucubepars_msm_table: - table_grating = tabdata['DISPERSER'].lower() - table_filter = tabdata['FILTER'].lower() - table_sroi = tabdata['ROISPATIAL'] - table_wroi = tabdata['ROISPECTRAL'] - table_power = tabdata['POWER'] - table_softrad = tabdata['SOFTRAD'] + table_grating = tabdata["DISPERSER"].lower() + table_filter = tabdata["FILTER"].lower() + table_sroi = tabdata["ROISPATIAL"] + table_wroi = tabdata["ROISPECTRAL"] + table_power = tabdata["POWER"] + table_softrad = tabdata["SOFTRAD"] if this_gwa == table_grating and this_filter == table_filter: - instrument_info.SetMSM(this_gwa, this_filter, - table_sroi, table_wroi, - table_power, table_softrad) + instrument_info.set_msm( + this_gwa, + this_filter, + table_sroi, + table_wroi, + table_power, + table_softrad, + ) # modified Shepard method e^-r weighting - elif weighting == 'emsm': + elif weighting == "emsm": for tabdata in ptab.ifucubepars_emsm_table: - table_grating = tabdata['DISPERSER'].lower() - table_filter = tabdata['FILTER'].lower() - table_sroi = tabdata['ROISPATIAL'] - table_wroi = tabdata['ROISPECTRAL'] - table_scalerad = tabdata['SCALERAD'] + table_grating = tabdata["DISPERSER"].lower() + table_filter = tabdata["FILTER"].lower() + table_sroi = tabdata["ROISPATIAL"] + table_wroi = tabdata["ROISPECTRAL"] + table_scalerad = tabdata["SCALERAD"] if this_gwa == table_grating and this_filter == table_filter: - instrument_info.SetEMSM(this_gwa, this_filter, - table_sroi, table_wroi, - table_scalerad) + instrument_info.set_emsm( + this_gwa, this_filter, table_sroi, table_wroi, table_scalerad + ) # read in wavelength tables - if weighting == 'msm': + if weighting == "msm": for tabdata in ptab.ifucubepars_prism_msm_wavetable: - table_wave = tabdata['WAVELENGTH'] - table_sroi = tabdata['ROISPATIAL'] - table_wroi = tabdata['ROISPECTRAL'] - table_power = tabdata['POWER'] - table_softrad = tabdata['SOFTRAD'] - instrument_info.SetPrismTable(table_wave, table_sroi, - table_wroi, table_power, - table_softrad) + table_wave = tabdata["WAVELENGTH"] + table_sroi = tabdata["ROISPATIAL"] + table_wroi = tabdata["ROISPECTRAL"] + table_power = tabdata["POWER"] + table_softrad = tabdata["SOFTRAD"] + instrument_info.set_prism_table( + table_wave, table_sroi, table_wroi, table_power, table_softrad + ) for tabdata in ptab.ifucubepars_med_msm_wavetable: - table_wave = tabdata['WAVELENGTH'] - table_sroi = tabdata['ROISPATIAL'] - table_wroi = tabdata['ROISPECTRAL'] - table_power = tabdata['POWER'] - table_softrad = tabdata['SOFTRAD'] - instrument_info.SetMedTable(table_wave, table_sroi, - table_wroi, table_power, - table_softrad) + table_wave = tabdata["WAVELENGTH"] + table_sroi = tabdata["ROISPATIAL"] + table_wroi = tabdata["ROISPECTRAL"] + table_power = tabdata["POWER"] + table_softrad = tabdata["SOFTRAD"] + instrument_info.set_med_table( + table_wave, table_sroi, table_wroi, table_power, table_softrad + ) for tabdata in ptab.ifucubepars_high_msm_wavetable: - table_wave = tabdata['WAVELENGTH'] - table_sroi = tabdata['ROISPATIAL'] - table_wroi = tabdata['ROISPECTRAL'] - table_power = tabdata['POWER'] - table_softrad = tabdata['SOFTRAD'] - instrument_info.SetHighTable(table_wave, table_sroi, - table_wroi, table_power, - table_softrad) - - elif weighting == 'emsm': + table_wave = tabdata["WAVELENGTH"] + table_sroi = tabdata["ROISPATIAL"] + table_wroi = tabdata["ROISPECTRAL"] + table_power = tabdata["POWER"] + table_softrad = tabdata["SOFTRAD"] + instrument_info.set_high_table( + table_wave, table_sroi, table_wroi, table_power, table_softrad + ) + + elif weighting == "emsm": for tabdata in ptab.ifucubepars_prism_emsm_wavetable: - table_wave = tabdata['WAVELENGTH'] - table_sroi = tabdata['ROISPATIAL'] - table_wroi = tabdata['ROISPECTRAL'] - table_scalerad = tabdata['SCALERAD'] - instrument_info.SetPrismEMSMTable(table_wave, table_sroi, - table_wroi, table_scalerad) + table_wave = tabdata["WAVELENGTH"] + table_sroi = tabdata["ROISPATIAL"] + table_wroi = tabdata["ROISPECTRAL"] + table_scalerad = tabdata["SCALERAD"] + instrument_info.set_prism_emsm_table( + table_wave, table_sroi, table_wroi, table_scalerad + ) for tabdata in ptab.ifucubepars_med_emsm_wavetable: - table_wave = tabdata['WAVELENGTH'] - table_sroi = tabdata['ROISPATIAL'] - table_wroi = tabdata['ROISPECTRAL'] - table_scalerad = tabdata['SCALERAD'] - instrument_info.SetMedEMSMTable(table_wave, table_sroi, - table_wroi, table_scalerad) + table_wave = tabdata["WAVELENGTH"] + table_sroi = tabdata["ROISPATIAL"] + table_wroi = tabdata["ROISPECTRAL"] + table_scalerad = tabdata["SCALERAD"] + instrument_info.set_med_emsm_table( + table_wave, table_sroi, table_wroi, table_scalerad + ) for tabdata in ptab.ifucubepars_high_emsm_wavetable: - table_wave = tabdata['WAVELENGTH'] - table_sroi = tabdata['ROISPATIAL'] - table_wroi = tabdata['ROISPECTRAL'] - table_scalerad = tabdata['SCALERAD'] - instrument_info.SetHighEMSMTable(table_wave, table_sroi, - table_wroi, table_scalerad) - - elif weighting == 'drizzle': + table_wave = tabdata["WAVELENGTH"] + table_sroi = tabdata["ROISPATIAL"] + table_wroi = tabdata["ROISPECTRAL"] + table_scalerad = tabdata["SCALERAD"] + instrument_info.set_high_emsm_table( + table_wave, table_sroi, table_wroi, table_scalerad + ) + + elif weighting == "drizzle": for tabdata in ptab.ifucubepars_prism_driz_wavetable: - table_wave = tabdata['WAVELENGTH'] - instrument_info.SetPrismDrizTable(table_wave) + table_wave = tabdata["WAVELENGTH"] + instrument_info.set_prism_driz_table(table_wave) for tabdata in ptab.ifucubepars_med_driz_wavetable: - table_wave = tabdata['WAVELENGTH'] - instrument_info.SetMedDrizTable(table_wave) + table_wave = tabdata["WAVELENGTH"] + instrument_info.set_med_driz_table(table_wave) for tabdata in ptab.ifucubepars_high_driz_wavetable: - table_wave = tabdata['WAVELENGTH'] - instrument_info.SetHighDrizTable(table_wave) + table_wave = tabdata["WAVELENGTH"] + instrument_info.set_high_driz_table(table_wave) ptab.close() del ptab diff --git a/jwst/cube_build/cube_build_step.py b/jwst/cube_build/cube_build_step.py index daa8e0d5fd..3bd220e421 100755 --- a/jwst/cube_build/cube_build_step.py +++ b/jwst/cube_build/cube_build_step.py @@ -1,6 +1,3 @@ -""" This is the main ifu spectral cube building routine. -""" - import time from jwst.datamodels import ModelContainer from jwst.lib.pipe_utils import match_nans_and_flags @@ -8,21 +5,20 @@ from . import ifu_cube from . import data_types import asdf -from asdf.exceptions import ValidationError from ..assign_wcs.util import update_s_region_keyword from ..stpipe import Step, record_step_status from pathlib import Path from astropy import units + __all__ = ["CubeBuildStep"] -class CubeBuildStep (Step): - """CubeBuildStep: Creates a 3-D spectral cube +class CubeBuildStep(Step): + """ + Create a 3-D spectral cube from IFU data. - Notes - ----- This is the controlling routine for building IFU Spectral Cubes. - It loads and sets the various input data and parameters need by + It loads and sets the various input data and parameters needed by the cube_build_step. This routine does the following operations: @@ -30,11 +26,10 @@ class CubeBuildStep (Step): 1. Extracts the input parameters from the cubepars reference file and merges them with any user-provided values. 2. Creates the output WCS from the input images and defines the mapping - between all the input arrays and the output array + between all the input arrays and the output array. 3. Passes the input data to the function to map all their input data to the output array. - 4. Updates the output data model with correct meta data - + 4. Updates the output data model with correct meta data. """ class_alias = "cube_build" @@ -67,27 +62,32 @@ class CubeBuildStep (Step): suffix = string(default='s3d') offset_file = string(default=None) # Filename containing a list of Ra and Dec offsets to apply to files. debug_spaxel = string(default='-1 -1 -1') # Default not used - """ # noqa: E501 + """ # noqa: E501 - reference_file_types = ['cubepar'] + reference_file_types = ["cubepar"] -# ________________________________________________________________________________ - def process(self, input): - """This is the main routine for IFU spectral cube building. + # ________________________________________________________________________________ + def process(self, input_data): + """ + Build an IFUCube from overlapping IFUImage data. Parameters ---------- - input : list of objects or str - list of datamodels or string name of input fits file or association. - """ + input_data : list of DataModel or str + List of datamodels or string of input fits filenames or association name. - self.log.info('Starting IFU Cube Building Step') + Returns + ------- + cube_container : ModelContainer + Container (list) of IFUCube models + """ + self.log.info("Starting IFU Cube Building Step") t0 = time.time() -# ________________________________________________________________________________ -# For all parameters convert to a standard format -# Report read in values to screen -# ________________________________________________________________________________ + # ________________________________________________________________________________ + # For all parameters convert to a standard format + # Report read in values to screen + # ________________________________________________________________________________ self.subchannel = self.band if not self.subchannel.islower(): @@ -103,240 +103,239 @@ def process(self, input): self.weighting = self.weighting.lower() if self.scalexy != 0.0: - self.log.info(f'Input Scale of axis 1 and 2 {self.scalexy}') + self.log.info(f"Input Scale of axis 1 and 2 {self.scalexy}") if self.scalew != 0.0: - self.log.info(f'Input wavelength scale {self.scalew}') + self.log.info(f"Input wavelength scale {self.scalew}") if self.wavemin is not None: - self.log.info(f'Setting minimum wavelength of spectral cube to: {self.wavemin}') + self.log.info(f"Setting minimum wavelength of spectral cube to: {self.wavemin}") if self.wavemax is not None: - self.log.info(f'Setting maximum wavelength of spectral cube to: {self.wavemax}') + self.log.info(f"Setting maximum wavelength of spectral cube to: {self.wavemax}") if self.rois != 0.0: - self.log.info(f'Input Spatial ROI size {self.rois}') + self.log.info(f"Input Spatial ROI size {self.rois}") if self.roiw != 0.0: - self.log.info(f'Input Wave ROI size {self.roiw}') + 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.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}') + 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.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}') + 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) # 3. internal_cal (local IFU - ifu cubes built in local IFU system) - if self.coord_system == 'world': - self.coord_system = 'skyalign' + if self.coord_system == "world": + self.coord_system = "skyalign" - self.interpolation = 'pointcloud' # initialize + self.interpolation = "pointcloud" # initialize # coord system = internal_cal only option for weighting = area ONLY USED FOR NIRSPEC - if self.coord_system == 'internal_cal': - self.interpolation = 'area' - self.weighting = 'emsm' + if self.coord_system == "internal_cal": + self.interpolation = "area" + self.weighting = "emsm" # if interpolation is point cloud then weighting can be # 1. MSM: modified Shepard method # 2. EMSM - if self.coord_system == 'skyalign': - self.interpolation = 'pointcloud' + if self.coord_system == "skyalign": + self.interpolation = "pointcloud" - if self.coord_system == 'ifualign': - self.interpolation = 'pointcloud' + if self.coord_system == "ifualign": + self.interpolation = "pointcloud" - if self.weighting == 'drizzle': - self.interpolation = 'drizzle' + if self.weighting == "drizzle": + self.interpolation = "drizzle" - self.log.info(f'Input interpolation: {self.interpolation}') - self.log.info(f'Coordinate system to use: {self.coord_system}') - if self.interpolation == 'pointcloud': - self.log.info(f'Weighting method for point cloud: {self.weighting}') + self.log.info(f"Input interpolation: {self.interpolation}") + self.log.info(f"Coordinate system to use: {self.coord_system}") + if self.interpolation == "pointcloud": + self.log.info(f"Weighting method for point cloud: {self.weighting}") if self.weight_power != 0: - self.log.info(f'Power weighting distance: {self.weight_power}') + self.log.info(f"Power weighting distance: {self.weight_power}") -# ________________________________________________________________________________ -# read input parameters - Channel, Band (Subchannel), Grating, Filter -# ________________________________________________________________________________ + # ________________________________________________________________________________ + # read input parameters - Channel, Band (Subchannel), Grating, Filter + # ________________________________________________________________________________ self.pars_input = {} -# the following parameters are set either by the an input parameter -# or -# if not set on the command line then from reading in the data. - self.pars_input['channel'] = [] - self.pars_input['subchannel'] = [] + # the following parameters are set either by the an input parameter + # or + # if not set on the command line then from reading in the data. + self.pars_input["channel"] = [] + self.pars_input["subchannel"] = [] - self.pars_input['filter'] = [] - self.pars_input['grating'] = [] + self.pars_input["filter"] = [] + self.pars_input["grating"] = [] # including values in pars_input that could get updated in cube_build_step.py - self.pars_input['coord_system'] = self.coord_system + self.pars_input["coord_system"] = self.coord_system if self.single: - self.pars_input['output_type'] = 'single' - self.log.info('Cube Type: Single cubes') - self.pars_input['coord_system'] = 'skyalign' + self.pars_input["output_type"] = "single" + self.log.info("Cube Type: Single cubes") + self.pars_input["coord_system"] = "skyalign" # Don't allow anything but drizzle, msm, or emsm weightings - if self.weighting not in ['msm', 'emsm', 'drizzle']: - self.weighting = 'drizzle' + if self.weighting not in ["msm", "emsm", "drizzle"]: + self.weighting = "drizzle" - if self.weighting == 'drizzle': - self.interpolation = 'drizzle' + if self.weighting == "drizzle": + self.interpolation = "drizzle" - if self.weighting == 'msm': - self.interpolation = 'pointcloud' + if self.weighting == "msm": + self.interpolation = "pointcloud" - if self.weighting == 'emsm': - self.interpolation = 'pointcloud' + if self.weighting == "emsm": + self.interpolation = "pointcloud" # read_user_input: - # see if options channel, band,grating filter are set on the command lines - # if they are then self.pars_input['output_type'] = 'user' and fill in par_input with values + # if options channel, band, grating, or filter are set on the command lines + # then set self.pars_input['output_type'] = 'user' and fill in par_input with values self.read_user_input() -# ________________________________________________________________________________ -# DataTypes -# Read in the input data - 4 formats are allowed: -# 1. filename -# 2. single model -# 3. ASN table -# 4. model container -# figure out what type of data we have. Fill in the input_table.input_models. -# input_table.input_models is used in the rest of IFU Cube Building -# We need to do this in cube_build_step because we need to pass the data_model -# to CRDS to figure out what type of reference files to grab (MIRI or NIRSPEC) -# if the user has provided the filename - strip out .fits and pull out the -# base name. The cube_build software will attached the needed information: -# channel, sub-channel grating or filter to filename -# ________________________________________________________________________________ - input_table = data_types.DataTypes(input, self.single, - self.output_file, - self.output_dir) + # ________________________________________________________________________________ + # DataTypes + # Read in the input data - 4 formats are allowed: + # 1. filename + # 2. single model + # 3. ASN table + # 4. model container + # figure out what type of data we have. Fill in the input_table.input_models. + # input_table.input_models is used in the rest of IFU Cube Building + # We need to do this in cube_build_step because we need to pass the data_model + # to CRDS to figure out what type of reference files to grab (MIRI or NIRSPEC) + # if the user has provided the filename - strip out .fits and pull out the + # base name. The cube_build software will attached the needed information: + # channel, sub-channel grating or filter to filename + # ________________________________________________________________________________ + input_table = data_types.DataTypes( + input_data, self.single, self.output_file, self.output_dir + ) self.input_models = input_table.input_models self.output_name_base = input_table.output_name -# ________________________________________________________________________________ + # ________________________________________________________________________________ -# Read in the first input model to determine with instrument we have -# output type is by default 'Channel' for MIRI and 'Band' for NIRSpec + # Read in the first input model to determine with instrument we have + # output type is by default 'Channel' for MIRI and 'Band' for NIRSpec instrument = self.input_models[0].meta.instrument.name.upper() if self.output_type is None: - if instrument == 'NIRSPEC': - self.output_type = 'band' - - 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}') -# ________________________________________________________________________________ -# If an offset file is provided do some basic checks on the file and its contents. -# The offset list contains a matching list to the files in the association -# used in calspec3 (for offline cube building). -# The offset list is an asdf file. + if instrument == "NIRSPEC": + self.output_type = "band" + + 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}") + # ________________________________________________________________________________ + # If an offset file is provided do some basic checks on the file and its contents. + # The offset list contains a matching list to the files in the association + # used in calspec3 (for offline cube building). + # The offset list is an asdf file. self.offsets = None if self.offset_file is not None: offsets = self.check_offset_file() if offsets is not None: self.offsets = offsets -# ________________________________________________________________________________ -# Read in Cube Parameter Reference file -# identify what reference file has been associated with these input + # ________________________________________________________________________________ + # Read in Cube Parameter Reference file + # identify what reference file has been associated with these input - par_filename = self.get_reference_file(self.input_models[0], 'cubepar') + par_filename = self.get_reference_file(self.input_models[0], "cubepar") # Check for a valid reference file - if par_filename == 'N/A': - self.log.warning('No default cube parameters reference file found') + if par_filename == "N/A": + self.log.warning("No default cube parameters reference file found") input_table.close() return -# ________________________________________________________________________________ -# shove the input parameters in to pars to pull out in general cube_build.py + # ________________________________________________________________________________ + # shove the input parameters in to pars to pull out in general cube_build.py pars = { - 'channel': self.pars_input['channel'], - 'subchannel': self.pars_input['subchannel'], - 'grating': self.pars_input['grating'], - 'filter': self.pars_input['filter'], - 'weighting': self.weighting, - 'single': self.single, - 'output_type': self.pars_input['output_type']} + "channel": self.pars_input["channel"], + "subchannel": self.pars_input["subchannel"], + "grating": self.pars_input["grating"], + "filter": self.pars_input["filter"], + "weighting": self.weighting, + "single": self.single, + "output_type": self.pars_input["output_type"], + } # shove the input parameters in to pars_cube to pull out ifu_cube.py # these parameters are related to the building a single ifucube_model pars_cube = { - 'scalexy': self.scalexy, - 'scalew': self.scalew, - 'interpolation': self.interpolation, - '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, - 'wavemax': self.wavemax, - 'offsets': self.offsets, - 'skip_dqflagging': self.skip_dqflagging, - 'suffix': self.suffix, - 'debug_spaxel': self.debug_spaxel} - -# ________________________________________________________________________________ -# create an instance of class CubeData + "scalexy": self.scalexy, + "scalew": self.scalew, + "interpolation": self.interpolation, + "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, + "wavemax": self.wavemax, + "offsets": self.offsets, + "skip_dqflagging": self.skip_dqflagging, + "suffix": self.suffix, + "debug_spaxel": self.debug_spaxel, + } + + # ________________________________________________________________________________ + # create an instance of class CubeData # Make sure all input models have consistent NaN and DO_NOT_USE values for model in self.input_models: match_nans_and_flags(model) - cubeinfo = cube_build.CubeData( - self.input_models, - par_filename, - **pars) -# ________________________________________________________________________________ -# cubeinfo.setup: -# read in all the input files, information from cube_pars, read in input data -# and fill in master_table holding what files are associated with each -# ch/sub-ch or grating/filter. -# Fill in all_channel, all_subchannel,all_filter, all_grating and instrument + cubeinfo = cube_build.CubeData(self.input_models, par_filename, **pars) + # ________________________________________________________________________________ + # cubeinfo.setup: + # read in all the input files, information from cube_pars, read in input data + # and fill in master_table holding what files are associated with each + # ch/sub-ch or grating/filter. + # Fill in all_channel, all_subchannel,all_filter, all_grating and instrument result = cubeinfo.setup() - instrument = result['instrument'] - instrument_info = result['instrument_info'] - master_table = result['master_table'] + instrument = result["instrument"] + instrument_info = result["instrument_info"] + master_table = result["master_table"] - if instrument == 'MIRI' and self.coord_system == 'internal_cal': - self.log.warning('The output coordinate system of internal_cal is not valid for MIRI') - self.log.warning('use output_coord = ifualign instead') + if instrument == "MIRI" and self.coord_system == "internal_cal": + self.log.warning("The output coordinate system of internal_cal is not valid for MIRI") + self.log.warning("use output_coord = ifualign instead") input_table.close() return - filenames = master_table.FileMap['filename'] + filenames = master_table.FileMap["filename"] self.pipeline = 3 - if self.output_type == 'multi' and len(filenames) == 1: + if self.output_type == "multi" and len(filenames) == 1: self.pipeline = 2 -# ________________________________________________________________________________ -# How many and what type of cubes will be made. -# send self.pars_input['output_type'], all_channel, all_subchannel, all_grating, all_filter -# return number of cubes and for each cube the fill in -# list_pars1 (valid channel or grating) and -# list_pars2 (value subchannel or filter) + # ________________________________________________________________________________ + # How many and what type of cubes will be made. + # send self.pars_input['output_type'], all_channel, all_subchannel, all_grating, all_filter + # return number of cubes and for each cube the fill in + # list_pars1 (valid channel or grating) and + # list_pars2 (value subchannel or filter) num_cubes, cube_pars = cubeinfo.number_cubes() if not self.single: - self.log.info(f'Number of IFU cubes produced by this run = {num_cubes}') + self.log.info(f"Number of IFU cubes produced by this run = {num_cubes}") # ModelContainer of ifucubes cube_container = ModelContainer() @@ -346,51 +345,52 @@ def process(self, input): # bands is done in outlier detection. for i in range(num_cubes): icube = str(i + 1) - list_par1 = cube_pars[icube]['par1'] - list_par2 = cube_pars[icube]['par2'] + list_par1 = cube_pars[icube]["par1"] + list_par2 = cube_pars[icube]["par2"] thiscube = ifu_cube.IFUCubeData( self.pipeline, self.input_models, self.output_name_base, - self.pars_input['output_type'], + self.pars_input["output_type"], instrument, list_par1, list_par2, instrument_info, master_table, - **pars_cube) -# ________________________________________________________________________________ + **pars_cube, + ) + # ________________________________________________________________________________ thiscube.check_ifucube() # basic checks -# Based on channel/subchannel or grating/prism find: -# spatial spaxel size, min wave, max wave -# Set linear_wavelength to true or false depending on which type of IFUcube is -# being created. -# If linear wavelength (single band) single values are assigned to: -# rois, roiw, weight_power, softrad -# If not linear wavelength (multi bands) an arrays are created for: -# rois, roiw, weight_power, softrad + # Based on channel/subchannel or grating/prism find: + # spatial spaxel size, min wave, max wave + # Set linear_wavelength to true or false depending on which type of IFUcube is + # being created. + # If linear wavelength (single band) single values are assigned to: + # rois, roiw, weight_power, softrad + # If not linear wavelength (multi bands) an arrays are created for: + # rois, roiw, weight_power, softrad - if self.coord_system == 'internal_cal': + if self.coord_system == "internal_cal": thiscube.determine_cube_parameters_internal() else: thiscube.determine_cube_parameters() thiscube.setup_ifucube_wcs() -# _______________________________________________________________________________ -# build the IFU Cube + # _______________________________________________________________________________ + # build the IFU Cube status = 0 -# If single = True: map each file to output grid and return single mapped file -# to output grid. # This option is used for background matching and outlier rejection + # If single = True: map each file to output grid and return single mapped file + # to output grid. # This option is used for background matching and outlier rejection if self.single: self.output_file = None cube_container = thiscube.build_ifucube_single() - self.log.info("Number of Single IFUCube models returned %i ", - len(cube_container)) + self.log.info("Number of Single IFUCube models returned %i ", len(cube_container)) -# Else standard IFU cube building - the result returned from build_ifucube will be 1 IFU CUBR + # Else standard IFU cube building + # the result returned from build_ifucube will be 1 IFU CUBE else: result, status = thiscube.build_ifucube() @@ -404,8 +404,7 @@ def process(self, input): del thiscube # irrelevant WCS keywords we will remove from final product - rm_keys = ['v2_ref', 'v3_ref', 'ra_ref', 'dec_ref', 'roll_ref', - 'v3yangle', 'vparity'] + rm_keys = ["v2_ref", "v3_ref", "ra_ref", "dec_ref", "roll_ref", "v3yangle", "vparity"] for cube in cube_container: footprint = cube.meta.wcs.footprint(axis_type="spatial") @@ -421,161 +420,177 @@ def process(self, input): record_step_status(cube_container, "cube_build", success=True) t1 = time.time() - self.log.debug(f'Time to build all cubes {t1-t0}') + self.log.debug(f"Time to build all cubes {t1 - t0}") input_table.close() return cube_container -# ****************************************************************************** + + # ****************************************************************************** def read_user_input(self): - """Read user input options for channel, subchannel, filter, or grating""" - - # Determine if any of the input parameters channel, band, filter or - # grating have been set. - - # This routine updates the dictionary self.pars_input with any user - # provided inputs. In particular it sets pars_input['channel'], - # pars_input['sub_channel'], pars_input['grating'], and - # pars_input['filter'] with user provided values. - - valid_channel = ['1', '2', '3', '4', 'all'] - valid_subchannel = ['short', 'medium', 'long', 'all', 'short-medium', 'short-long', - 'medium-short', 'medium-long', 'long-short', 'long-medium'] - - valid_fwa = ['f070lp', 'f100lp', - 'g170lp', 'f290lp', 'clear', 'opaque', 'all'] - valid_gwa = ['g140m', 'g140h', 'g235m', 'g235h', - 'g395m', 'g395h', 'prism', 'all'] -# ________________________________________________________________________________ -# For MIRI we can set the channel. -# If channel is set to 'all' then let the determine_band_coverage figure out -# which channels are covered by the data. - - if self.channel == 'all': - self.pars_input['channel'].append('all') + """ + Read user input options for channel, subchannel, filter, or grating. + + Determine if any of the input parameters channel, band, filter or + grating have been set by the user. + This routine updates the dictionary self.pars_input with any user + provided inputs. In particular it sets pars_input['channel'], + pars_input['sub_channel'], pars_input['grating'], and + pars_input['filter'] with user provided values. + """ + valid_channel = ["1", "2", "3", "4", "all"] + valid_subchannel = [ + "short", + "medium", + "long", + "all", + "short-medium", + "short-long", + "medium-short", + "medium-long", + "long-short", + "long-medium", + ] + + valid_fwa = ["f070lp", "f100lp", "g170lp", "f290lp", "clear", "opaque", "all"] + valid_gwa = ["g140m", "g140h", "g235m", "g235h", "g395m", "g395h", "prism", "all"] + # ________________________________________________________________________________ + # For MIRI we can set the channel. + # If channel is set to 'all' then let the determine_band_coverage figure out + # which channels are covered by the data. + + if self.channel == "all": + self.pars_input["channel"].append("all") else: # user has set value if not self.single: - self.pars_input['output_type'] = 'user' - channellist = self.channel.split(',') + self.pars_input["output_type"] = "user" + channellist = self.channel.split(",") user_clen = len(channellist) for j in range(user_clen): ch = channellist[j] if user_clen > 1: - ch = ch.strip('[') - ch = ch.strip(']') - ch = ch.strip(' ') + ch = ch.strip("[") + ch = ch.strip("]") + ch = ch.strip(" ") ch = ch[1:-1] ch = str(ch) if ch in valid_channel: - self.pars_input['channel'].append(ch) -# remove duplicates if needed - self.pars_input['channel'] = list(set(self.pars_input['channel'])) -# ________________________________________________________________________________ -# For MIRI we can set the subchannel -# if set to all then let the determine_band_coverage figure out what subchannels -# are covered by the data - - if self.subchannel == 'all': - self.pars_input['subchannel'].append('all') + self.pars_input["channel"].append(ch) + # remove duplicates if needed + self.pars_input["channel"] = list(set(self.pars_input["channel"])) + # ________________________________________________________________________________ + # For MIRI we can set the subchannel + # if set to all then let the determine_band_coverage figure out what subchannels + # are covered by the data + + if self.subchannel == "all": + self.pars_input["subchannel"].append("all") else: # user has set value if not self.single: - self.pars_input['output_type'] = 'user' - subchannellist = self.subchannel.split(',') + self.pars_input["output_type"] = "user" + subchannellist = self.subchannel.split(",") user_blen = len(subchannellist) for j in range(user_blen): b = subchannellist[j] if user_blen > 1: - b = b.strip('[') - b = b.strip(']') - b = b.strip(' ') + b = b.strip("[") + b = b.strip("]") + b = b.strip(" ") b = b[1:-1] b = str(b) if b in valid_subchannel: - self.pars_input['subchannel'].append(b) -# remove duplicates if needed - self.pars_input['subchannel'] = list(set(self.pars_input['subchannel'])) -# ________________________________________________________________________________ -# For NIRSPEC we can set the filter -# If set to all then let the determine_band_coverage figure out what filters are -# covered by the data. - - if self.filter == 'all': - self.pars_input['filter'].append('all') - else: # User has set value + self.pars_input["subchannel"].append(b) + # remove duplicates if needed + self.pars_input["subchannel"] = list(set(self.pars_input["subchannel"])) + # ________________________________________________________________________________ + # For NIRSPEC we can set the filter + # If set to all then let the determine_band_coverage figure out what filters are + # covered by the data. + + if self.filter == "all": + self.pars_input["filter"].append("all") + else: # User has set value if not self.single: - self.pars_input['output_type'] = 'user' - filterlist = self.filter.split(',') + self.pars_input["output_type"] = "user" + filterlist = self.filter.split(",") user_flen = len(filterlist) for j in range(user_flen): f = filterlist[j] if user_flen > 1: - f = f.strip('[') - f = f.strip(']') - f = f.strip(' ') + f = f.strip("[") + f = f.strip("]") + f = f.strip(" ") f = f[1:-1] f = str(f) if f in valid_fwa: - self.pars_input['filter'].append(f) -# remove duplicates if needed - self.pars_input['filter'] = list(set(self.pars_input['filter'])) -# ________________________________________________________________________________ -# For NIRSPEC we can set the grating -# If set to all then let the determine_band_coverage figure out what gratings are -# covered by the data - if self.grating == 'all': - self.pars_input['grating'].append('all') - else: # user has set value + self.pars_input["filter"].append(f) + # remove duplicates if needed + self.pars_input["filter"] = list(set(self.pars_input["filter"])) + # ________________________________________________________________________________ + # For NIRSPEC we can set the grating + # If set to all then let the determine_band_coverage figure out what gratings are + # covered by the data + if self.grating == "all": + self.pars_input["grating"].append("all") + else: # user has set value if not self.single: - self.pars_input['output_type'] = 'user' - gratinglist = self.grating.split(',') + self.pars_input["output_type"] = "user" + gratinglist = self.grating.split(",") user_glen = len(gratinglist) for j in range(user_glen): - g = gratinglist[j] if user_glen > 1: - g = g.strip('[') - g = g.strip(']') - g = g.strip(' ') + g = g.strip("[") + g = g.strip("]") + g = g.strip(" ") g = g[1:-1] g = str(g) if g in valid_gwa: - self.pars_input['grating'].append(g) -# remove duplicates if needed - self.pars_input['grating'] = list(set(self.pars_input['grating'])) -# ________________________________________________________________________________ + self.pars_input["grating"].append(g) + # remove duplicates if needed + self.pars_input["grating"] = list(set(self.pars_input["grating"])) + + # ________________________________________________________________________________ def check_offset_file(self): - """Read in an optional ra and dec offset for each file. + """ + Read in an optional ra and dec offset for each file. - Summary - ---------- - Check that is file is asdf file. - Check the file has the correct format using an local schema file. - The schema file, ifuoffset.schema.yaml, is located in the jwst/cube_build directory. - For each file in the input association check that there is a corresponding + Check that the offset file is an asdf file. + Check that the file has the correct format using an local schema file. + For each file in the input association check that there is a corresponding file in the offset file. - """ - + Returns + ------- + offsets : dict + Dictionary containing offset information. + """ # validate the offset file using the schema file - DATA_PATH = Path(__file__).parent + data_path = Path(__file__).parent try: - af = asdf.open(self.offset_file, custom_schema=DATA_PATH/'ifuoffset.schema.yaml') - except ValidationError: - schema_message = ('Validation Error for offset file. Fix the offset file. \n' + \ - 'The offset file needs to have the same number of elements ' + \ - 'in the three lists: filename, raoffset and decoffset.\n' +\ - 'The units need to provided and only arcsec is allowed.') - - raise Exception(schema_message) - - offset_filename = af['filename'] - offset_ra = af['raoffset'] - offset_dec = af['decoffset'] + af = asdf.open(self.offset_file, custom_schema=data_path / "ifuoffset.schema.yaml") + except ValueError: + self.log.error( + "Validation Error for offset file. Fix the offset file. " + " The offset file needs to have the same number of elements " + " in the three lists: filename, raoffset and decoffset." + " The units need to provided and only arcsec is allowed." + ) + + raise ValueError( + "Offset file is not correct. Offset file needs to have three lists: " + "filename, raoffset and decoffset all of the same length." + ) from None + + offset_filename = af["filename"] + offset_ra = af["raoffset"] + offset_dec = af["decoffset"] # Note: - # af['units'] is checked by the schema validation. It must be arcsec or a validation error occurs. + # af['units'] is checked by the schema validation. + # It must be arcsec or a validation error occurs. # check that all the file names in input_model are in the offset filename for model in self.input_models: @@ -584,24 +599,30 @@ def check_offset_file(self): continue else: af.close() - raise ValueError('Error in offset file. A file in the association is not found in offset list %s', file_check) + raise ValueError( + "Error in offset file. A file in the association is not found in offset " + f" list {file_check}" + ) # check that all the lists have the same length - len_file = len(offset_filename) + len_file = len(offset_filename) len_ra = len(offset_ra) len_dec = len(offset_dec) - if (len_file != len_ra or len_ra != len_dec or len_file != len_dec): + if len_file != len_ra or len_ra != len_dec or len_file != len_dec: af.close() - raise ValueError('The offset file does not have the same number of values for filename, raoffset, decoffset') + raise ValueError( + "The offset file does not have the same number of values for " + " filename, raoffset, decoffset" + ) - offset_ra = offset_ra* units.arcsec - offset_dec = offset_dec* units.arcsec + offset_ra = offset_ra * units.arcsec + offset_dec = offset_dec * units.arcsec # The offset file has passed tests so set the offset dictionary offsets = {} - offsets['filename'] = offset_filename - offsets['raoffset'] = offset_ra - offsets['decoffset'] = offset_dec + offsets["filename"] = offset_filename + offsets["raoffset"] = offset_ra + offsets["decoffset"] = offset_dec af.close() return offsets diff --git a/jwst/cube_build/cube_build_wcs_util.py b/jwst/cube_build/cube_build_wcs_util.py index 86a1177ea0..d9b631e3a6 100644 --- a/jwst/cube_build/cube_build_wcs_util.py +++ b/jwst/cube_build/cube_build_wcs_util.py @@ -1,5 +1,5 @@ -""" Routines related to WCS procedures of cube_build -""" +"""Routines related to WCS procedures of cube_build.""" + import numpy as np from ..assign_wcs import nirspec from ..assign_wcs.util import wrap_ra @@ -11,8 +11,9 @@ # ****************************************************************************** -def find_corners_MIRI(input, this_channel, instrument_info, coord_system): - """ For MIRI channel data find the footprint of this data on the sky +def find_corners_miri(input_data, this_channel, instrument_info, coord_system): + """ + For MIRI channel data find the footprint of this data on the sky. For a specific channel on an exposure find the min and max of the spatial coordinates, either in alpha,beta or ra,dec depending @@ -21,34 +22,50 @@ def find_corners_MIRI(input, this_channel, instrument_info, coord_system): Parameters ---------- - input : data model - input model (or file) + input_data : IFUImage model + Input model (or file) this_channel : str - channel working with - instrument_info : dictionary - dictionary holding x pixel min and max values for each channel + Channel working with + instrument_info : dict + Dictionary holding x pixel min and max values for each channel coord_system : str - coordinate system of output cube: skyalign, ifualign, internal_cal + Coordinate system of output cube: skyalign, ifualign, internal_cal Returns ------- - corners of the spatial system determined from min and max of spatial values - min and max of wavelength - + a_min : float + Minimum value of coord 1 - along axis 1 of the IFUcube + b1 : float + Coord 2 value corresponding to a_min + a_max : float + Maximum value of coord 1 - along the axis 1 of the IFUcube + b2 : float + Coord 2 value corresponding to a_min + a1 : float + Coord 1 value coorsponding to b_min + b_min : float + Minimum value of coord 2 - along the axis 2 of the IFU cube + a2 : float + Coord 1 value coorsponding to b_max + b_max : float + Maximum value of coord 2 - along the axis 2 of the IFU cube + lambda_min : float + Minimum wavelength + lambda_max : float + Maximum wavelength """ # x,y values for channel - convert to output coordinate system # return the min & max of spatial coords and wavelength - xstart, xend = instrument_info.GetMIRISliceEndPts(this_channel) - ysize = input.data.shape[0] + xstart, xend = instrument_info.get_miri_slice_endpts(this_channel) + ysize = input_data.data.shape[0] y, x = np.mgrid[:ysize, xstart:xend] - if coord_system == 'internal_cal': + if coord_system == "internal_cal": # coord1 = along slice # coord2 = across slice - detector2alpha_beta = input.meta.wcs.get_transform('detector', - 'alpha_beta') + detector2alpha_beta = input_data.meta.wcs.get_transform("detector", "alpha_beta") coord1, coord2, lam = detector2alpha_beta(x, y) valid = np.logical_and(np.isfinite(coord1), np.isfinite(coord2)) @@ -78,7 +95,7 @@ def find_corners_MIRI(input, this_channel, instrument_info, coord_system): else: # skyalign or ifualign # coord1 = ra # coord2 = dec - coord1, coord2, lam = input.meta.wcs(x, y) + coord1, coord2, lam = input_data.meta.wcs(x, y) valid = np.logical_and(np.isfinite(coord1), np.isfinite(coord2)) coord1 = coord1[valid] coord2 = coord2[valid] @@ -111,7 +128,7 @@ def find_corners_MIRI(input, this_channel, instrument_info, coord_system): lambda_min = np.nanmin(lam) lambda_max = np.nanmax(lam) - if coord_system != 'internal_cal': + if coord_system != "internal_cal": # before returning, ra should be between 0 to 360 a_min %= 360 a_max %= 360 @@ -120,11 +137,14 @@ def find_corners_MIRI(input, this_channel, instrument_info, coord_system): a2 %= 360 return a_min, b1, a_max, b2, a1, b_min, a2, b_max, lambda_min, lambda_max + + # ***************************************************************************** -def find_corners_NIRSPEC(input, instrument_info, coord_system): - """Find the sky footprint of a slice of a NIRSpec exposure +def find_corners_nirspec(input_data, coord_system): + """ + Find the sky footprint of a slice of a NIRSpec exposure. For each slice find: a. the min and max spatial coordinates (along slice, across slice) or @@ -133,17 +153,33 @@ def find_corners_NIRSPEC(input, instrument_info, coord_system): Parameters ---------- - input: data model - input model (or file) + input_data : IFUImageModel + Input calibrated model (or file) coord_system : str - coordinate system of output cube: skyalign, ifualign, internal_cal + Coordinate system of output cube: skyalign, ifualign, internal_cal - Notes - ----- Returns ------- - min and max spatial coordinates and wavelength for slice. - + a_min : float + Minimum value of coord 1 - along axis 1 of the IFUcube + b1 : float + Coord 2 value corresponding to a_min + a_max : float + Maximum value of coord 1 - along the axis 1 of the IFUcube + b2 : float + Coord 2 value corresponding to a_min + a1 : float + Coord 1 value coorsponding to b_min + b_min : float + Minimum value of coord 2 - along the axis 2 of the IFU cube + a2 : float + Coord 1 value coorsponding to b_max + b_max : float + Maximum value of coord 2 - along the axis 2 of the IFU cube + lambda_min : float + Minimum wavelength + lambda_max : float + Maximum wavelength """ nslices = 30 a_slice = np.zeros(nslices * 2) @@ -153,19 +189,17 @@ def find_corners_NIRSPEC(input, instrument_info, coord_system): lambda_slice = np.zeros(nslices * 2) k = 0 # for NIRSPEC there are 30 regions - log.info('Looping over slices to determine cube size') + log.info("Looping over slices to determine cube size") - wcsobj, tr1, tr2, tr3 = nirspec._get_transforms(input, np.arange(nslices)) + wcsobj, tr1, tr2, tr3 = nirspec._get_transforms(input_data, np.arange(nslices)) # noqa: SLF001 for i in range(nslices): - slice_wcs = nirspec._nrs_wcs_set_input_lite(input, wcsobj, i, - [tr1, tr2[i], tr3[i]]) - x, y = wcstools.grid_from_bounding_box(slice_wcs.bounding_box, - step=(1, 1), center=True) - if coord_system == 'internal_cal': + slice_wcs = nirspec._nrs_wcs_set_input_lite(input_data, wcsobj, i, [tr1, tr2[i], tr3[i]]) # noqa: SLF001 + x, y = wcstools.grid_from_bounding_box(slice_wcs.bounding_box, step=(1, 1), center=True) + if coord_system == "internal_cal": # coord1 = along slice # coord2 = across slice - detector2slicer = slice_wcs.get_transform('detector', 'slicer') + detector2slicer = slice_wcs.get_transform("detector", "slicer") coord2, coord1, lam = detector2slicer(x, y) # lam ~0 for this transform valid = np.logical_and(np.isfinite(coord1), np.isfinite(coord2)) coord1 = coord1[valid] @@ -189,7 +223,7 @@ def find_corners_NIRSPEC(input, instrument_info, coord_system): # ra range coord1_wrap = wrap_ra(coord1) coord1 = coord1_wrap -# ________________________________________________________________________________ + # ________________________________________________________________________________ coord1 = coord1.flatten() coord2 = coord2.flatten() lam = lam.flatten() @@ -217,8 +251,8 @@ def find_corners_NIRSPEC(input, instrument_info, coord_system): lambda_slice[k + 1] = np.nanmax(lam) k = k + 2 -# ________________________________________________________________________________ -# now test the ra slices for consistency. Adjust if needed. + # ________________________________________________________________________________ + # now test the ra slices for consistency. Adjust if needed. a_min = np.nanmin(a_slice) a_max = np.nanmax(a_slice) diff --git a/jwst/cube_build/cube_internal_cal.py b/jwst/cube_build/cube_internal_cal.py index 105768425a..73f48336f6 100644 --- a/jwst/cube_build/cube_internal_cal.py +++ b/jwst/cube_build/cube_internal_cal.py @@ -1,6 +1,5 @@ -""" Routines for creating single band, single exposure IFU Cubes with -the interpolation method = area , coord_system = internal_cal -""" +"""Routines for creating IFUCubes with interpolation method = area , coord_system = internal_cal.""" + import numpy as np from stdatamodels.jwst.datamodels import dqflags from stdatamodels.jwst.transforms.models import _toindex @@ -8,51 +7,88 @@ from .cube_match_internal import cube_wrapper_internal # c extension -def match_det2cube(instrument, - x, y, sliceno, - input_model, - transform, - acoord, zcoord, - crval_along, crval3, - cdelt_along, cdelt3, - naxis1, naxis2): - """ Match detector pixels to output plane in local IFU coordinate system +def match_det2cube( + instrument, + x, + y, + sliceno, + input_model, + transform, + acoord, + zcoord, + crval_along, + crval3, + cdelt_along, + cdelt3, + naxis1, + naxis2, +): + """ + Match detector pixels to output plane in local IFU coordinate system. This routine assumes a 1-1 mapping in across slice to slice no. This routine assumes the output coordinate systems is local IFU plane. The user can not change scaling in across slice dimension - Map the corners of the x,y detector values to a cube defined by local IFU plane. + Map the corners of the x, y detector values to a cube defined by local IFU plane. In the along slice, lambda plane find the % area of the detector pixel which it overlaps with in the cube. For each spaxel record the detector pixels that overlap with it - store flux, % overlap, beta_distance. Parameters ---------- - x : numpy.ndarray - x values of pixels in slice - y : numpy.ndarray - y values of pixels in slice + x : ndarray + Detector x pixel values in the slice + y : ndarray + Detector y pixel values in the slice sliceno : int - slice number - input_model : datamodel - input slope model or file + Slice number + input_model : IFUImage model + Input calibrated model or file transform : transform - wcs transform to transform x,y to alpha,beta, lambda - spaxel : list - list of spaxels holding information on each cube pixel. + WCS transform to transform x, y to alpha, beta, lambda + acoord : ndarray + Array of along slice valuee defining the along slice spatial dimension the of IFU cube + zcoord : ndarray + Array of wavelength values defining the wavelength dimension of the IFU cube + crval_along : float + Along slixe reference value in the IFU cube + crval3 : float + Wavelength reference value in the IFU cube + cdelt_along : float + Along slice reference value in the IFU cube + cdelt3 : float + Wavelength sampling in the IFU cube + naxis1 : int + Size of axis 1 of the IFU cube + naxis2 : int + Size of axis 2 of the IFU cube Returns ------- - spaxel filled in with needed information on overlapping detector pixels + result : tuple + Results from running c code to determine internal coordinate system IFU Cubes + Values contained in result: + instrument_no: integer. NIRSpec = 1, MIRI = 0 + naxis1 & naxis2: output axis of IFU cube + crval_along : reference value in along slice dimension + cdelt_along : sampling in along slice dimension + crval3 : reference value in wavelength dimension + cdelt3 : wavelength sampling + a1, a2, a3, a4: Array of corners of pixels holding along slice coordinates + lam1, lam2, lam3, lam4: Array of corners of pixels holding wavelength coordinates + acoord: array holding slice coordinates of the IFU internal cube + zcoord: array holding wavelength coordinates of the IFU internal cube + ss: array holding the slice number of each pixel + pixel_flux: array of pixel fluxes + pixel_err: array of pixel errors """ - x = _toindex(x) y = _toindex(y) pixel_dq = input_model.dq[y, x] - all_flags = (dqflags.pixel['DO_NOT_USE'] + dqflags.pixel['NON_SCIENCE']) + all_flags = dqflags.pixel["DO_NOT_USE"] + dqflags.pixel["NON_SCIENCE"] # find the location of all the values to reject in cube building - good_data = np.where((np.bitwise_and(pixel_dq, all_flags) == 0)) + good_data = np.where(np.bitwise_and(pixel_dq, all_flags) == 0) # good data holds the location of pixels we want to map to cube x = x[good_data] @@ -67,7 +103,7 @@ def match_det2cube(instrument, xx_left = x xx_right = x + 1 # along slice dimension is second coordinate returned from transform - if instrument == 'NIRSPEC': + if instrument == "NIRSPEC": b1, a1, lam1 = transform(xx_left, yy_bot) b2, a2, lam2 = transform(xx_right, yy_bot) b3, a3, lam3 = transform(xx_right, yy_top) @@ -81,7 +117,7 @@ def match_det2cube(instrument, lam3 = lam3 * 1.0e6 lam4 = lam4 * 1.0e6 - if instrument == 'MIRI': + if instrument == "MIRI": a1, b1, lam1 = transform(xx_left, yy_bot) a2, b2, lam2 = transform(xx_right, yy_bot) a3, b3, lam3 = transform(xx_right, yy_top) @@ -116,16 +152,33 @@ def match_det2cube(instrument, pixel_err = input_model.err[y, x] # 1-1 mapping in across slice direction (x for NIRSPEC, y for MIRI) - if instrument == 'NIRSPEC': + if instrument == "NIRSPEC": ss = sliceno instrument_no = 1 - if instrument == 'MIRI': + if instrument == "MIRI": ss = sliceno instrument_no = 0 - result = cube_wrapper_internal(instrument_no, naxis1, naxis2, - crval_along, cdelt_along, crval3, cdelt3, - a1, a2, a3, a4, lam1, lam2, lam3, lam4, - acoord, zcoord, ss, - pixel_flux, pixel_err) + result = cube_wrapper_internal( + instrument_no, + naxis1, + naxis2, + crval_along, + cdelt_along, + crval3, + cdelt3, + a1, + a2, + a3, + a4, + lam1, + lam2, + lam3, + lam4, + acoord, + zcoord, + ss, + pixel_flux, + pixel_err, + ) return result diff --git a/jwst/cube_build/data_types.py b/jwst/cube_build/data_types.py index 123f214ed2..54ddc4eb45 100644 --- a/jwst/cube_build/data_types.py +++ b/jwst/cube_build/data_types.py @@ -1,51 +1,39 @@ -""" Class Data Type is used to read in the input data. It also determines -if the input data is a single science exposure, an association table, a -single datamodel or several data models stored in a ModelContainer. -""" -import os - from stdatamodels.jwst import datamodels - from jwst.datamodels import ModelContainer - import logging +from pathlib import Path + log = logging.getLogger(__name__) log.setLevel(logging.DEBUG) -# ****************************************************************************** -class DataTypes(): +class DataTypes: + """Class to handle reading input data to cube_build.""" - """ Class to handle reading input data to cube_build. - """ + template = { + "asn_rule": "", + "target": "", + "asn_pool": "", + "asn_type": "", + "products": [{"name": "", "members": [{"exptype": "", "expname": ""}]}], + } - template = {"asn_rule": "", - "target": "", - "asn_pool": "", - "asn_type": "", - "products": [ - {"name": "", - "members": [ - {"exptype": "", - "expname": ""} - ] - } - ] - } - - def __init__(self, input, single, output_file, output_dir): - """ Read in input data and determine what type of input data. + def __init__(self, input_data, single, output_file, output_dir): + """ + Assemble input data for processing. Open the input data using datamodels and determine if data is a single input model, an association, or a set of input models - contained in a ModelContainer. + contained in a ModelContainer. The method populates the self.input_models + which is a list of input models. An initial base name for the output file + is constructed. Parameters ---------- - input : datamodel or ModelContainer + input_data : str, IFUImageModel or ModelContainer Input data to cube_build either a filename, single model, association table, or a ModelContainer - single : boolean + single : bool If True then creating single mode IFUCubes for outlier detection or mrs_matching. If false then creating standard IFUcubes output_file : str @@ -53,67 +41,58 @@ def __init__(self, input, single, output_file, output_dir): output_dir : str Optional user provided output directory name - Notes - ----- - The method populates the self.input_models which is list of input models. - An initial base name for the output file is constructed. - Raises ------ - NotIFUImageModel + NotIFUImageModelError IFU data was not the input data TypeError Input data was not of correct form - """ - self.input_models = [] self.output_name = None - # open the input with datamodels - # if input is filename or model when it is opened it is a model - # if input if an association name or ModelContainer then it is opened as a container + # open the input_data with datamodels + # if input_data is filename or model when it is opened it is a model + # if input_data is an association name or ModelContainer then it is opened as a container - input_models = datamodels.open(input) - # if input is a filename, we will need to close the opened file + input_models = datamodels.open(input_data) + # if input_data is a filename, we will need to close the opened file self._opened = [input_models] if isinstance(input_models, datamodels.IFUImageModel): # It's a single image that's been passed in as a model - # input is a model + # input_data is a model filename = input_models.meta.filename self.input_models.append(input_models) self.output_name = self.build_product_name(filename) elif isinstance(input_models, ModelContainer): - self.output_name = 'Temp' + self.output_name = "Temp" self.input_models = input_models if not single: # find the name of the output file from the association self.output_name = input_models.asn_table["products"][0]["name"] else: # close files opened above self.close() - raise TypeError("Failed to process file type {}".format(type(input_models))) + raise TypeError(f"Failed to process file type {type(input_models)}") # If the user has set the output name, strip off *.fits. # Suffixes will be added to this name later, to designate the # channel+subchannel (MIRI MRS) or grating+filter (NRS IFU) the output cube covers. if output_file is not None: - basename, ext = os.path.splitext(os.path.basename(output_file)) - self.output_name = basename + self.output_name = Path(output_file).stem if output_dir is not None: - self.output_name = output_dir + '/' + self.output_name + self.output_name = output_dir + "/" + self.output_name def close(self): - """ - Close any files opened by this instance - """ + """Close any files opened by this instance.""" [f.close() for f in self._opened] -# _______________________________________________________________________________ + # _______________________________________________________________________________ def build_product_name(self, filename): - """ Determine the base of output name if an input data is a fits filename. + """ + Determine the base of output name if an input data is a FITS filename. Parameters ---------- @@ -126,10 +105,9 @@ def build_product_name(self, filename): single_product : str Output base filename. """ - - indx = filename.rfind('.fits') - indx_try = filename.rfind('_rate.fits') - indx_try2 = filename.rfind('_cal.fits') + indx = filename.rfind(".fits") + indx_try = filename.rfind("_rate.fits") + indx_try2 = filename.rfind("_cal.fits") if indx_try > 0: single_product = filename[:indx_try] @@ -139,10 +117,11 @@ def build_product_name(self, filename): single_product = filename[:indx] return single_product + # _______________________________________________________________________________ -class NotIFUImageModel(Exception): - """ Raise Exception if data is not of type IFUImageModel - """ +class NotIFUImageModelError(Exception): + """Raise Exception if data is not of type IFUImageModel.""" + pass diff --git a/jwst/cube_build/file_table.py b/jwst/cube_build/file_table.py index 73863a83d5..8917027883 100644 --- a/jwst/cube_build/file_table.py +++ b/jwst/cube_build/file_table.py @@ -1,155 +1,158 @@ -""" Dictionary holding defaults for cube_build -""" from stdatamodels.jwst import datamodels import logging + log = logging.getLogger(__name__) log.setLevel(logging.DEBUG) -class FileTable(): - """ Dictionary contains defaults for MIRI and NIRSPEC data - """ +class FileTable: + """Dictionary contains defaults for MIRI and NIRSPEC data.""" def __init__(self): - self.FileMap = {} - self.FileMap['filename'] = [] - self.FileMap['MIRI'] = {} - self.FileMap['MIRI']['1'] = {} - self.FileMap['MIRI']['1']['short'] = [] - self.FileMap['MIRI']['1']['medium'] = [] - self.FileMap['MIRI']['1']['long'] = [] - self.FileMap['MIRI']['1']['short-medium'] = [] - self.FileMap['MIRI']['1']['short-long'] = [] - self.FileMap['MIRI']['1']['medium-short'] = [] - self.FileMap['MIRI']['1']['medium-long'] = [] - self.FileMap['MIRI']['1']['long-short'] = [] - self.FileMap['MIRI']['1']['long-medium'] = [] - - self.FileMap['MIRI']['2'] = {} - self.FileMap['MIRI']['2']['short'] = [] - self.FileMap['MIRI']['2']['medium'] = [] - self.FileMap['MIRI']['2']['long'] = [] - self.FileMap['MIRI']['2']['short-medium'] = [] - self.FileMap['MIRI']['2']['short-long'] = [] - self.FileMap['MIRI']['2']['medium-short'] = [] - self.FileMap['MIRI']['2']['medium-long'] = [] - self.FileMap['MIRI']['2']['long-short'] = [] - self.FileMap['MIRI']['2']['long-medium'] = [] - - self.FileMap['MIRI']['3'] = {} - self.FileMap['MIRI']['3']['short'] = [] - self.FileMap['MIRI']['3']['medium'] = [] - self.FileMap['MIRI']['3']['long'] = [] - self.FileMap['MIRI']['3']['short-medium'] = [] - self.FileMap['MIRI']['3']['short-long'] = [] - self.FileMap['MIRI']['3']['medium-short'] = [] - self.FileMap['MIRI']['3']['medium-long'] = [] - self.FileMap['MIRI']['3']['long-short'] = [] - self.FileMap['MIRI']['3']['long-medium'] = [] - - self.FileMap['MIRI']['4'] = {} - self.FileMap['MIRI']['4']['short'] = [] - self.FileMap['MIRI']['4']['medium'] = [] - self.FileMap['MIRI']['4']['long'] = [] - self.FileMap['MIRI']['4']['short-medium'] = [] - self.FileMap['MIRI']['4']['short-long'] = [] - self.FileMap['MIRI']['4']['medium-short'] = [] - self.FileMap['MIRI']['4']['medium-long'] = [] - self.FileMap['MIRI']['4']['long-short'] = [] - self.FileMap['MIRI']['4']['long-medium'] = [] - - self.FileMap['NIRSPEC'] = {} - self.FileMap['NIRSPEC']['prism'] = {} - self.FileMap['NIRSPEC']['prism']['clear'] = [] - self.FileMap['NIRSPEC']['prism']['opaque'] = [] - - self.FileMap['NIRSPEC']['g140m'] = {} - self.FileMap['NIRSPEC']['g140m']['f070lp'] = [] - self.FileMap['NIRSPEC']['g140m']['f100lp'] = [] - self.FileMap['NIRSPEC']['g140m']['opaque'] = [] - - self.FileMap['NIRSPEC']['g140h'] = {} - self.FileMap['NIRSPEC']['g140h']['f070lp'] = [] - self.FileMap['NIRSPEC']['g140h']['f100lp'] = [] - self.FileMap['NIRSPEC']['g140h']['opaque'] = [] - - self.FileMap['NIRSPEC']['g235m'] = {} - self.FileMap['NIRSPEC']['g235m']['f170lp'] = [] - self.FileMap['NIRSPEC']['g235m']['opaque'] = [] - - self.FileMap['NIRSPEC']['g235h'] = {} - self.FileMap['NIRSPEC']['g235h']['f170lp'] = [] - self.FileMap['NIRSPEC']['g235h']['opaque'] = [] - - self.FileMap['NIRSPEC']['g395m'] = {} - self.FileMap['NIRSPEC']['g395m']['f290lp'] = [] - self.FileMap['NIRSPEC']['g395m']['opaque'] = [] - - self.FileMap['NIRSPEC']['g395h'] = {} - self.FileMap['NIRSPEC']['g395h']['f290lp'] = [] - self.FileMap['NIRSPEC']['g395h']['opaque'] = [] - -# ******************************************************************************** - - def set_file_table(self, - input_models): + self.FileMap["filename"] = [] + self.FileMap["MIRI"] = {} + self.FileMap["MIRI"]["1"] = {} + self.FileMap["MIRI"]["1"]["short"] = [] + self.FileMap["MIRI"]["1"]["medium"] = [] + self.FileMap["MIRI"]["1"]["long"] = [] + self.FileMap["MIRI"]["1"]["short-medium"] = [] + self.FileMap["MIRI"]["1"]["short-long"] = [] + self.FileMap["MIRI"]["1"]["medium-short"] = [] + self.FileMap["MIRI"]["1"]["medium-long"] = [] + self.FileMap["MIRI"]["1"]["long-short"] = [] + self.FileMap["MIRI"]["1"]["long-medium"] = [] + + self.FileMap["MIRI"]["2"] = {} + self.FileMap["MIRI"]["2"]["short"] = [] + self.FileMap["MIRI"]["2"]["medium"] = [] + self.FileMap["MIRI"]["2"]["long"] = [] + self.FileMap["MIRI"]["2"]["short-medium"] = [] + self.FileMap["MIRI"]["2"]["short-long"] = [] + self.FileMap["MIRI"]["2"]["medium-short"] = [] + self.FileMap["MIRI"]["2"]["medium-long"] = [] + self.FileMap["MIRI"]["2"]["long-short"] = [] + self.FileMap["MIRI"]["2"]["long-medium"] = [] + + self.FileMap["MIRI"]["3"] = {} + self.FileMap["MIRI"]["3"]["short"] = [] + self.FileMap["MIRI"]["3"]["medium"] = [] + self.FileMap["MIRI"]["3"]["long"] = [] + self.FileMap["MIRI"]["3"]["short-medium"] = [] + self.FileMap["MIRI"]["3"]["short-long"] = [] + self.FileMap["MIRI"]["3"]["medium-short"] = [] + self.FileMap["MIRI"]["3"]["medium-long"] = [] + self.FileMap["MIRI"]["3"]["long-short"] = [] + self.FileMap["MIRI"]["3"]["long-medium"] = [] + + self.FileMap["MIRI"]["4"] = {} + self.FileMap["MIRI"]["4"]["short"] = [] + self.FileMap["MIRI"]["4"]["medium"] = [] + self.FileMap["MIRI"]["4"]["long"] = [] + self.FileMap["MIRI"]["4"]["short-medium"] = [] + self.FileMap["MIRI"]["4"]["short-long"] = [] + self.FileMap["MIRI"]["4"]["medium-short"] = [] + self.FileMap["MIRI"]["4"]["medium-long"] = [] + self.FileMap["MIRI"]["4"]["long-short"] = [] + self.FileMap["MIRI"]["4"]["long-medium"] = [] + + self.FileMap["NIRSPEC"] = {} + self.FileMap["NIRSPEC"]["prism"] = {} + self.FileMap["NIRSPEC"]["prism"]["clear"] = [] + self.FileMap["NIRSPEC"]["prism"]["opaque"] = [] + + self.FileMap["NIRSPEC"]["g140m"] = {} + self.FileMap["NIRSPEC"]["g140m"]["f070lp"] = [] + self.FileMap["NIRSPEC"]["g140m"]["f100lp"] = [] + self.FileMap["NIRSPEC"]["g140m"]["opaque"] = [] + + self.FileMap["NIRSPEC"]["g140h"] = {} + self.FileMap["NIRSPEC"]["g140h"]["f070lp"] = [] + self.FileMap["NIRSPEC"]["g140h"]["f100lp"] = [] + self.FileMap["NIRSPEC"]["g140h"]["opaque"] = [] + + self.FileMap["NIRSPEC"]["g235m"] = {} + self.FileMap["NIRSPEC"]["g235m"]["f170lp"] = [] + self.FileMap["NIRSPEC"]["g235m"]["opaque"] = [] + + self.FileMap["NIRSPEC"]["g235h"] = {} + self.FileMap["NIRSPEC"]["g235h"]["f170lp"] = [] + self.FileMap["NIRSPEC"]["g235h"]["opaque"] = [] + + self.FileMap["NIRSPEC"]["g395m"] = {} + self.FileMap["NIRSPEC"]["g395m"]["f290lp"] = [] + self.FileMap["NIRSPEC"]["g395m"]["opaque"] = [] + + self.FileMap["NIRSPEC"]["g395h"] = {} + self.FileMap["NIRSPEC"]["g395h"]["f290lp"] = [] + self.FileMap["NIRSPEC"]["g395h"]["opaque"] = [] + + # ******************************************************************************** + + def set_file_table(self, input_models): """ - Short Summary - ------------- - Fill in the MasterTable which holds the files that the cube will be constructed - from. Since MIRI has 2 channels per image this MASTERTable helps to figure out + Set up the master_table dictionary. + + Fill in the master_table which holds the files that the cube will be constructed + from. Since MIRI has 2 channels per image this MASTERTable dictionary helps to figure out which data needs to be use. - THe MasterTable for MIRI is broken down by channel and subchannel. + The master_table for MIRI is broken down by channel and subchannel. For each channel/subchannel combination - a file is listed that covers those options - For NIRSPEC the table contains the Grating and Filter for each file. + For NIRSPEC the table contains the grating and filter for each file. + Parameters + ---------- + input_models : IFUImageModel + The input datamodels used the set up the class FileTable Returns ------- - MasterTable filled in with files needed - instrument name + master_table : dict + Dictionary containing the filename/model for each channel/band or grating/filter """ -# ________________________________________________________________________________ -# Loop over input list of files and assign fill in the MasterTable with filename -# for the correct (channel-subchannel) or (grating-subchannel) + # ________________________________________________________________________________ + # Loop over input list of files and fill in the master_table with filename + # for the correct (channel-subchannel) or (grating-subchannel) for model in input_models: instrument = model.meta.instrument.name.upper() assign_wcs = model.meta.cal_step.assign_wcs filename = model.meta.filename - self.FileMap['filename'].append(filename) + self.FileMap["filename"].append(filename) if not isinstance(model, datamodels.IFUImageModel): - raise NotIFUImageModel(f"Input data is not a IFUImageModel, instead it is {model}") - if assign_wcs != 'COMPLETE': - raise ErrorNoAssignWCS("Assign WCS has not been run on file %s", - model.meta.filename) + raise NotIFUImageModelError( + f"Input data is not a IFUImageModel, instead it is {model}" + ) + if assign_wcs != "COMPLETE": + raise NoAssignWCSError(f"Assign WCS has not been run on file {model.meta.filename}") # _____________________________________________________________________ # MIRI instrument - if instrument == 'MIRI': + if instrument == "MIRI": channel = model.meta.instrument.channel subchannel = model.meta.instrument.band.lower() clenf = len(channel) for k in range(clenf): - self.FileMap['MIRI'][channel[k]][subchannel].append(model) + self.FileMap["MIRI"][channel[k]][subchannel].append(model) # _____________________________________________________________________ # NIRSPEC instrument - elif instrument == 'NIRSPEC': + elif instrument == "NIRSPEC": fwa = model.meta.instrument.filter.lower() gwa = model.meta.instrument.grating.lower() - self.FileMap['NIRSPEC'][gwa][fwa].append(model) + self.FileMap["NIRSPEC"][gwa][fwa].append(model) else: - log.info('Instrument not valid for cube') + log.info("Instrument not valid for cube") pass return instrument -class ErrorNoAssignWCS(Exception): +class NoAssignWCSError(Exception): + """Raise Exception if assign_wcs has not been run.""" + pass -class NotIFUImageModel(Exception): - """ Raise Exception if data is not of type IFUImageModel - """ +class NotIFUImageModelError(Exception): + """Raise Exception if data is not of type IFUImageModel.""" + pass diff --git a/jwst/cube_build/ifu_cube.py b/jwst/cube_build/ifu_cube.py index 7f82444227..fa166f965d 100644 --- a/jwst/cube_build/ifu_cube.py +++ b/jwst/cube_build/ifu_cube.py @@ -1,7 +1,4 @@ -""" Work horse routines used for building ifu spectra cubes -(including the main loop over files and the construction of -final spaxel fluxes) -""" +"""Work horse routines used for building ifu spectra cubes.""" import numpy as np import logging @@ -32,20 +29,59 @@ log.setLevel(logging.DEBUG) -class IFUCubeData(): - - def __init__(self, - pipeline, - input_models, - output_name_base, - output_type, - instrument, - list_par1, - list_par2, - instrument_info, - master_table, - **pars_cube): - """ Class IFUCube holds the high level data for each IFU Cube +class IFUCubeData: + """Combine IFU data onto a regular grid.""" + + def __init__( + self, + pipeline, + input_models, + output_name_base, + output_type, + instrument, + list_par1, + list_par2, + instrument_info, + master_table, + **pars_cube, + ): + """ + Initialize the IFUCube. + + Parameters + ---------- + pipeline : int + Integer that indicates which pipeline is being run. Pipeline = 2 for + calwebb_spec and pipeline = 3 for calwebb_spec3. + input_models : ModelContainer + Container of IFUIMageModels + output_name_base : str + String defining the name of the output product. Unlike other steps, + cube_build determines the name of the output file because the bands + of the data are determined after reading in the data. In addition, + the name of the output also depend if calwebb_spec2 or calwebb_spec3 + is being run and if the user has selected particular band to use. + output_type : str + The type of cube to created. The possibilities are 'band','channel', + 'grating','multi'. The default for calwebb_spec2 is to create 'multi' + band cubes. In the case of MIRI, the "multi' band cubes from + calwebb_spec2 contain the two channels from single calibration fits file. + The default type of cube for calwebb_spec3 is 'band'. These cubes will + contain a single band of data. + instrument : str + Instrument name, either "MIRI" or "NIRSpec" + list_par1 : list + For MIRI the list contains the channels used to build the IFU, while + for NIRSpec the list contains the gratinging used. + list_par2 : list + For MIRI the list contains the sub-channels uses to build the IFU, while + for NIRSpec the list contains the filters used. + instrument_info : dict + Dictionary containing information on the basic instrument parameters. + master_table : dict + Dictionary of containing the files covering each band. + **pars_cube : dict + Dictionary of parameters controlling how the cube is built. """ self.input_models_this_cube = [] # list of files use to make cube working on @@ -64,31 +100,33 @@ def __init__(self, self.master_table = master_table self.output_type = output_type - 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.offsets = pars_cube.get('offsets') - self.rois = pars_cube.get('rois') - self.roiw = pars_cube.get('roiw') - self.debug_spaxel = pars_cube.get('debug_spaxel') - - self.spaxel_x, self.spaxel_y, self.spaxel_z = [int(val) for val in self.debug_spaxel.split()] + 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.offsets = pars_cube.get("offsets") + self.rois = pars_cube.get("rois") + self.roiw = pars_cube.get("roiw") + self.debug_spaxel = pars_cube.get("debug_spaxel") + + self.spaxel_x, self.spaxel_y, self.spaxel_z = [ + int(val) for val in self.debug_spaxel.split() + ] self.spatial_size = None self.spectral_size = None - self.interpolation = pars_cube.get('interpolation') - self.coord_system = pars_cube.get('coord_system') - self.wavemin = pars_cube.get('wavemin') - self.wavemax = pars_cube.get('wavemax') - self.weighting = pars_cube.get('weighting') - self.weight_power = pars_cube.get('weight_power') - self.skip_dqflagging = pars_cube.get('skip_dqflagging') - self.suffix = pars_cube.get('suffix') + self.interpolation = pars_cube.get("interpolation") + self.coord_system = pars_cube.get("coord_system") + self.wavemin = pars_cube.get("wavemin") + self.wavemax = pars_cube.get("wavemax") + self.weighting = pars_cube.get("weighting") + self.weight_power = pars_cube.get("weight_power") + self.skip_dqflagging = pars_cube.get("skip_dqflagging") + self.suffix = pars_cube.get("suffix") self.num_bands = 0 - self.output_name = '' + self.output_name = "" self.wavemin_user = False # Check for NIRSpec if user has set wavelength limits self.wavemax_user = False @@ -129,21 +167,20 @@ def __init__(self, self.tolerance_dq_overlap = 0.05 # spaxel has to have 5% overlap to flag in FOV self.overlap_partial = 4 # intermediate flag - self.overlap_full = 2 # intermediate flag - self.overlap_hole = dqflags.pixel['DO_NOT_USE'] - self.overlap_no_coverage = dqflags.pixel['NON_SCIENCE'] + self.overlap_full = 2 # intermediate flag + self.overlap_hole = dqflags.pixel["DO_NOT_USE"] + self.overlap_no_coverage = dqflags.pixel["NON_SCIENCE"] - # ************************************************************** + # ________________________________________________________________________________ def check_ifucube(self): - """ Perform some quick checks that the type of cube to be produced - conforms to rules + """ + Perform some quick checks that the type of cube to be produced conforms to the rules. Raises ------ - IncorrectInput + IncorrectInputError Interpolation = area was selected for when input data is more than one file or model - """ num1 = len(self.list_par1) num_files = 0 @@ -157,30 +194,46 @@ def check_ifucube(self): # do some basic checks on the cubes if self.coord_system == "internal_cal": if num_files > 1: - raise IncorrectInput("Cubes built in internal_cal coordinate system" + - " are built from a single file, not multiple exposures") + raise IncorrectInputError( + "Cubes built in internal_cal coordinate system" + " are built from a single file, not multiple exposures" + ) if len(self.list_par1) > 1: - raise IncorrectInput("Only a single channel or grating " + - " can be used to create cubes in internal_cal coordinate system." + - " Use --output_type=band") + raise IncorrectInputError( + "Only a single channel or grating " + " can be used to create cubes in internal_cal coordinate system." + " Use --output_type=band" + ) # ________________________________________________________________________________ def define_cubename(self): - """ Define the base output name + """ + Define the base output name. + + Usually the output name is defined by the association table. However in the case + of cube_build several cubes can be created from a single call of cube_build. The + user can override the type of data to combine to make a cube. It is left to cube_build + to determine which channels, bands, gratings or filters are used to make the IFUCube. + The final name includes the channel/band (MIRI) or grating/filter (NIRSpec). + + Returns + ------- + newname : str + Output name of the IFU cube. """ if self.pipeline == 2: - newname = self.output_name_base + '_' + self.suffix + '.fits' + newname = self.output_name_base + "_" + self.suffix + ".fits" else: - if self.instrument == 'MIRI': - + if self.instrument == "MIRI": # Check to see if the output base name already contains the # field "clear", which sometimes shows up in IFU product # names created by the ASN rules. If so, strip it off, so # that the remaining suffixes created below form the entire # list of optical elements in the final output name. - suffix = self.output_name_base[self.output_name_base.rfind('_') + 1:] - if suffix in ['clear']: - self.output_name_base = self.output_name_base[:self.output_name_base.rfind('_')] + basename = self.output_name_base + suffix = basename[basename.rfind("_") + 1 :] + if suffix in ["clear"]: + self.output_name_base = basename[: basename.rfind("_")] # Now compose the appropriate list of optical element suffix names # based on MRS channel and sub-channel @@ -189,64 +242,67 @@ def define_cubename(self): if ch not in channels: channels.append(ch) number_channels = len(channels) - ch_name = '_ch' + ch_name = "_ch" for i in range(number_channels): ch_name = ch_name + channels[i] if i < number_channels - 1: - ch_name = ch_name + '-' + ch_name = ch_name + "-" # Sort by inverse alphabetical, e.g. short -> medium -> long - subchannels = sorted(list(set(self.list_par2)))[::-1] + subchannels = sorted(set(self.list_par2))[::-1] log.info(f"Subchannel listing: {subchannels}") number_subchannels = len(subchannels) - b_name = '' + b_name = "" for i in range(number_subchannels): b_name = b_name + subchannels[i] b_name = b_name.lower() - newname = self.output_name_base + ch_name + '-' + b_name - if self.coord_system == 'internal_cal': - newname = self.output_name_base + ch_name + '-' + b_name + '_internal' - if self.output_type == 'single': - newname = self.output_name_base + ch_name + '-' + b_name + '_single' + newname = self.output_name_base + ch_name + "-" + b_name + if self.coord_system == "internal_cal": + newname = self.output_name_base + ch_name + "-" + b_name + "_internal" + if self.output_type == "single": + newname = self.output_name_base + ch_name + "-" + b_name + "_single" # ________________________________________________________________________________ - elif self.instrument == 'NIRSPEC': - + elif self.instrument == "NIRSPEC": # Check to see if the output base name already has a grating/prism # suffix attached. If so, strip it off, and let the following logic # add all necessary grating and filter suffixes. - suffix = self.output_name_base[self.output_name_base.rfind('_') + 1:] - if suffix in ['g140m', 'g235m', 'g395m', 'g140h', 'g235h', 'g395h', 'prism']: - self.output_name_base = self.output_name_base[:self.output_name_base.rfind('_')] + basename = self.output_name_base + suffix = basename[basename.rfind("_") + 1 :] + if suffix in ["g140m", "g235m", "g395m", "g140h", "g235h", "g395h", "prism"]: + self.output_name_base = basename[: basename.rfind("_")] - fg_name = '_' + fg_name = "_" for i in range(len(self.list_par1)): - fg_name = fg_name + self.list_par1[i] + '-' + self.list_par2[i] + fg_name = fg_name + self.list_par1[i] + "-" + self.list_par2[i] if i < self.num_bands - 1: - fg_name = fg_name + '-' + fg_name = fg_name + "-" fg_name = fg_name.lower() newname = self.output_name_base + fg_name - if self.output_type == 'single': - newname = self.output_name_base + fg_name + '_single' - if self.coord_system == 'internal_cal': - newname = self.output_name_base + fg_name + '_internal' - # ______________________________________________________________________________ - if self.output_type != 'single': - log.info(f'Output Name: {newname}') + if self.output_type == "single": + newname = self.output_name_base + fg_name + "_single" + if self.coord_system == "internal_cal": + newname = self.output_name_base + fg_name + "_internal" + + if self.output_type != "single": + log.info(f"Output Name: {newname}") return newname - # _______________________________________________________________________ + # _______________________________________________________________________ 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. - + """ + Set up the WCS of the cube in the tangent plane. Parameters ---------- - footprint: tuple - holds min and max or ra,dec, and wavelength for the cube - footprint + corner_a : ndarray + Array of RA corners of the footprint of all input data. + corner_b : ndarray + Array of Dec corners of the footprint of all input data. + lambda_min : float + Minimum wavelength value of the data. + lambda_max : float + Maximum wavelength value of the data. """ - ra_min = np.min(corner_a) ra_max = np.max(corner_a) @@ -282,15 +338,14 @@ def set_geometry(self, corner_a, corner_b, lambda_min, lambda_max): num = len(corner_a) for i in range(num): - xi, eta = coord.radec2std(self.crval1, self.crval2, - corner_a[i], corner_b[i], rot_angle) + xi, eta = coord.radec2std(self.crval1, self.crval2, 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) eta_max = max(eta_corner) - # ________________________________________________________________________________ + # find the CRPIX1 CRPIX2 - xi and eta centered at 0,0 # to find location of center abs of min values is how many pixels # we want a symmetric cube centered on xi,eta = 0 @@ -305,9 +360,9 @@ def set_geometry(self, corner_a, corner_b, lambda_min, lambda_max): # 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) + na = int(self.nspax_x / 2) if self.nspax_y is not None: - nb = int(self.nspax_y/2) + 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) @@ -328,20 +383,24 @@ def set_geometry(self, corner_a, corner_b, lambda_min, lambda_max): # center of spaxels self.xcoord = np.zeros(self.naxis1) xstart = xi_min + self.cdelt1 / 2.0 - self.xcoord = np.arange(start=xstart, stop=xstart + self.naxis1 * self.cdelt1, step=self.cdelt1) + self.xcoord = np.arange( + start=xstart, stop=xstart + self.naxis1 * self.cdelt1, step=self.cdelt1 + ) self.ycoord = np.zeros(self.naxis2) ystart = eta_min + self.cdelt2 / 2.0 - self.ycoord = np.arange(start=ystart, stop=ystart + self.naxis2 * self.cdelt2, step=self.cdelt2) + self.ycoord = np.arange( + start=ystart, stop=ystart + self.naxis2 * self.cdelt2, step=self.cdelt2 + ) # depending on the naxis and cdelt values the x,ycoord can have 1 more element than naxis. # Clean up arrays dropping extra values at the end. - self.xcoord = self.xcoord[0:self.naxis1] - self.ycoord = self.ycoord[0:self.naxis2] + self.xcoord = self.xcoord[0 : self.naxis1] + self.ycoord = self.ycoord[0 : self.naxis2] xv, yv = np.meshgrid(self.xcoord, self.ycoord) self.xcenters = xv.flatten() self.ycenters = yv.flatten() - # _______________________________________________________________________ + # set up the lambda (z) coordinate of the cube self.cdelt3_normal = None if self.linear_wavelength: @@ -361,7 +420,7 @@ def set_geometry(self, corner_a, corner_b, lambda_min, lambda_max): self.crpix3 = 1.0 zstart = self.lambda_min + self.cdelt3 / 2.0 self.zcoord = np.arange(start=zstart, stop=self.lambda_max, step=self.cdelt3) - self.zcoord = self.zcoord[0:self.naxis3] + self.zcoord = self.zcoord[0 : self.naxis3] else: self.naxis3 = len(self.wavelength_table) self.zcoord = np.asarray(self.wavelength_table) @@ -374,22 +433,26 @@ def set_geometry(self, corner_a, corner_b, lambda_min, lambda_max): cdelt3_normal[self.naxis3 - 1] = cdelt3_normal[self.naxis3 - 2] self.cdelt3_normal = cdelt3_normal - # _______________________________________________________________________ - def set_geometryAB(self, corner_a, corner_b, lambda_min, lambda_max): - """Based on the along slice, across slice and wavelength footprint set up the - size of the cube in internal IFU plane. + # _______________________________________________________________________ + def set_geometry_slicer(self, corner_a, corner_b, lambda_min, lambda_max): + """ + Set up the size of the cube in the internal IFU plane. - This will be a single exposure cube - small FOV assume - rectangular coord system. + This will be a single exposure cube. The internal IFU Cube is in the + slicer plane and is defined by along slice and across slice coordinates. Parameters ---------- - footprint : tuple - Holds the min and max alpha, beta and wavelength values of - cube on sky + corner_a : ndarray + Array of along slice corners of the footprint of all input data. + corner_b : ndarray + Array of across slice corners of the footprint of all input data. + lambda_min : float + Minimum wavelength value of the data. + lambda_max : float + Maximum wavelength value of the data. """ - self.a_min = np.min(corner_a) self.a_max = np.max(corner_a) @@ -403,8 +466,7 @@ def set_geometryAB(self, corner_a, corner_b, lambda_min, lambda_max): alimit = max(np.abs(self.a_min), np.abs(self.a_max)) range_b = self.b_max - self.b_min - if self.instrument == 'MIRI': - # self.cdelt1 = self.cdelt2 # make cubes same scaling. MIRI EC team requested this removed (2/16/21) + if self.instrument == "MIRI": along_cdelt = self.cdelt1 n1a = math.ceil(alimit / along_cdelt) @@ -419,7 +481,7 @@ def set_geometryAB(self, corner_a, corner_b, lambda_min, lambda_max): across_naxis = self.naxis2 across_cdelt = self.cdelt2 - if self.instrument == 'NIRSPEC': + if self.instrument == "NIRSPEC": along_cdelt = self.cdelt2 n1a = math.ceil(alimit / along_cdelt) @@ -452,19 +514,19 @@ def set_geometryAB(self, corner_a, corner_b, lambda_min, lambda_max): across_coord[i] = start start = start + across_cdelt - if self.instrument == 'MIRI': + if self.instrument == "MIRI": self.crval1 = self.a_min self.xcoord = acoord self.ycoord = across_coord self.crval2 = self.b_min - if self.instrument == 'NIRSPEC': + if self.instrument == "NIRSPEC": self.crval2 = self.a_min self.ycoord = acoord self.xcoord = across_coord self.crval1 = self.b_min -# _______________________________________________________________________ -# common to both MIRI and NIRSPEC + # _______________________________________________________________________ + # common to both MIRI and NIRSPEC self.crpix1 = 0.5 self.crpix2 = 0.5 @@ -486,57 +548,82 @@ def set_geometryAB(self, corner_a, corner_b, lambda_min, lambda_max): for i in range(self.naxis3): self.zcoord[i] = zstart zstart = zstart + self.cdelt3 -# _______________________________________________________________________ + # _______________________________________________________________________ def print_cube_geometry(self): - """Print out the general properties of the size of the IFU Cube - """ - - log.info('Cube Geometry:') - if self.coord_system == 'internal_cal': - log.info('axis# Naxis CRPIX CRVAL CDELT(arcsec) Min & Max (along slice, across slice)') + """Print out the general properties of the size of the IFU Cube.""" + log.info("Cube Geometry:") + if self.coord_system == "internal_cal": + log.info( + "axis# Naxis CRPIX CRVAL CDELT(arcsec) " + " Min & Max (along slice, across slice)" + ) else: - log.info('axis# Naxis CRPIX CRVAL CDELT(arcsec) Min & Max (xi, eta arcsec)') - log.info('Axis 1 %5d %5.2f %12.8f %12.8f %12.8f %12.8f', - self.naxis1, self.crpix1, self.crval1, self.cdelt1, - self.a_min, self.a_max) - log.info('Axis 2 %5d %5.2f %12.8f %12.8f %12.8f %12.8f', - self.naxis2, self.crpix2, self.crval2, self.cdelt2, - self.b_min, self.b_max) + log.info("axis# Naxis CRPIX CRVAL CDELT(arcsec) Min & Max (xi, eta arcsec)") + log.info( + "Axis 1 %5d %5.2f %12.8f %12.8f %12.8f %12.8f", + self.naxis1, + self.crpix1, + self.crval1, + self.cdelt1, + self.a_min, + self.a_max, + ) + log.info( + "Axis 2 %5d %5.2f %12.8f %12.8f %12.8f %12.8f", + self.naxis2, + self.crpix2, + self.crval2, + self.cdelt2, + self.b_min, + self.b_max, + ) if self.linear_wavelength: - log.info('axis# Naxis CRPIX CRVAL CDELT(microns) Min & Max (microns)') - log.info('Axis 3 %5d %5.2f %12.8f %12.8f %12.8f %12.8f', - self.naxis3, self.crpix3, self.crval3, self.cdelt3, - self.lambda_min, self.lambda_max) + log.info("axis# Naxis CRPIX CRVAL CDELT(microns) Min & Max (microns)") + log.info( + "Axis 3 %5d %5.2f %12.8f %12.8f %12.8f %12.8f", + self.naxis3, + self.crpix3, + self.crval3, + self.cdelt3, + self.lambda_min, + self.lambda_max, + ) if not self.linear_wavelength: - log.info('Non-linear wavelength dimension; CDELT3 variable') - log.info('axis# Naxis CRPIX CRVAL Min & Max (microns)') - log.info('Axis 3 %5d %5.2f %12.8f %12.8f %12.8f', - self.naxis3, self.crpix3, self.crval3, - self.wavelength_table[0], self.wavelength_table[self.naxis3 - 1]) + log.info("Non-linear wavelength dimension; CDELT3 variable") + log.info("axis# Naxis CRPIX CRVAL Min & Max (microns)") + log.info( + "Axis 3 %5d %5.2f %12.8f %12.8f %12.8f", + self.naxis3, + self.crpix3, + self.crval3, + self.wavelength_table[0], + self.wavelength_table[self.naxis3 - 1], + ) if self.rot_angle is not None: - log.info('Rotation angle between Ra-Dec and Slicer-Plane %12.8f', self.rot_angle) + log.info("Rotation angle between RA-Dec and Slicer-Plane %12.8f", self.rot_angle) - if self.instrument == 'MIRI': + if self.instrument == "MIRI": # length of channel and subchannel are the same number_bands = len(self.list_par1) for i in range(number_bands): this_channel = self.list_par1[i] this_subchannel = self.list_par2[i] - log.info(f'Cube covers channel, subchannel: {this_channel}, {this_subchannel}') - elif self.instrument == 'NIRSPEC': + log.info(f"Cube covers channel, subchannel: {this_channel}, {this_subchannel}") + elif self.instrument == "NIRSPEC": # number of filters and gratings are the same number_bands = len(self.list_par1) for i in range(number_bands): this_fwa = self.list_par2[i] this_gwa = self.list_par1[i] - log.info(f'Cube covers grating, filter: {this_gwa}, {this_fwa}') -# ________________________________________________________________________________ + log.info(f"Cube covers grating, filter: {this_gwa}, {this_fwa}") + # ________________________________________________________________________________ def build_ifucube(self): - """ Create the IFU cube + """ + Create an IFU cube. 1. Loop over every band contained in the IFU cube and read in the data associated with the band @@ -553,14 +640,12 @@ def build_ifucube(self): to determine the overlap in the along slice-wavelength plane. 4. find_spaxel_flux: find the final flux associated with each spaxel 5. setup_final_ifucube_model - 6. output_ifucube Returns ------- - Returns an ifu cube - + result : IFUCubeModel + An IFU cube of combined IFU image data. """ - self.output_name = self.define_cubename() total_num = self.naxis1 * self.naxis2 * self.naxis3 @@ -569,25 +654,25 @@ def build_ifucube(self): self.spaxel_var = np.zeros(total_num, dtype=np.float64) self.spaxel_iflux = np.zeros(total_num, dtype=np.float64) self.spaxel_dq = np.zeros(total_num, dtype=np.uint32) - # ______________________________________________________________________________ nxyplane = self.naxis1 * self.naxis2 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') + elif self.spaxel_z < 0 or self.spaxel_x < 0 or self.spaxel_y < 0: + log.info("Incorrect input for Debug Spaxel values. Counting starts at 0") debug_cube_index = -1 - print(self.spaxel_z, self.spaxel_x, self.spaxel_y) + log.info(f"{self.spaxel_z} {self.spaxel_x} {self.spaxel_y}") else: spaxel_z = self.spaxel_z spaxel_x = self.spaxel_x spaxel_y = self.spaxel_y debug_cube_index = spaxel_z * (nxyplane) + spaxel_y * self.naxis1 + spaxel_x - log.info(f"Printing debug information for cube spaxel: {spaxel_x} {spaxel_y} {spaxel_z}") + log.info( + f"Printing debug information for cube spaxel: {spaxel_x} {spaxel_y} {spaxel_z}" + ) - # ______________________________________________________________________________ subtract_background = True # loop over every file that covers this channel/subchannel (MIRI) or @@ -600,7 +685,6 @@ def build_ifucube(self): this_par1 = self.list_par1[ib] this_par2 = self.list_par2[ib] for input_model in self.master_table.FileMap[self.instrument][this_par1][this_par2]: - # ________________________________________________________________________________ # loop over the files that cover the spectral range the cube is for self.input_models_this_cube.append(input_model.copy()) @@ -610,17 +694,29 @@ def build_ifucube(self): input_model_ref = input_model log.debug(f"Working on Band defined by: {this_par1} {this_par2}") - # -------------------------------------------------------------------------------- - # POINTCLOUD used for skyalign and IFUalign - # -------------------------------------------------------------------------------- - if self.interpolation in ['pointcloud', 'drizzle']: - pixelresult = self.map_detector_to_outputframe(this_par1, - subtract_background, - input_model) - - coord1, coord2, corner_coord, wave, dwave, flux, err, slice_no, rois_pixel, \ - roiw_pixel, weight_pixel, softrad_pixel, scalerad_pixel, \ - x_det, y_det = pixelresult + + if self.interpolation in ["pointcloud", "drizzle"]: + pixelresult = self.map_detector_to_outputframe( + this_par1, subtract_background, input_model + ) + + ( + coord1, + coord2, + corner_coord, + wave, + dwave, + flux, + err, + slice_no, + rois_pixel, + roiw_pixel, + weight_pixel, + softrad_pixel, + scalerad_pixel, + x_det, + y_det, + ) = pixelresult # by default flag the dq plane based on the FOV of the detector projected to sky flag_dq_plane = 1 @@ -631,42 +727,62 @@ def build_ifucube(self): # If all the data is flagged as DO_NOT_USE - not common- then log warning build_cube = True if wave is None: - log.warning(f'No valid data found on file {input_model.meta.filename}') + log.warning(f"No valid data found on file {input_model.meta.filename}") flag_dq_plane = 0 build_cube = False - # ______________________________________________________________________ + # C extension setup - # ______________________________________________________________________ start_region = 0 end_region = 0 - if self.instrument == 'MIRI': + if self.instrument == "MIRI": instrument = 0 - start_region = self.instrument_info.GetStartSlice(this_par1) - end_region = self.instrument_info.GetEndSlice(this_par1) + start_region = self.instrument_info.get_start_slice(this_par1) + end_region = self.instrument_info.get_end_slice(this_par1) else: # NIRSPEC instrument = 1 result = None weight_type = 0 # default to emsm instead of msm - if self.weighting == 'msm': + if self.weighting == "msm": weight_type = 1 - if self.interpolation == 'pointcloud' and build_cube: + if self.interpolation == "pointcloud" and build_cube: roiw_ave = np.mean(roiw_pixel) - result = cube_wrapper(instrument, flag_dq_plane, weight_type, start_region, end_region, - self.overlap_partial, self.overlap_full, - self.xcoord, self.ycoord, self.zcoord, - coord1, coord2, wave, flux, err, slice_no, - rois_pixel, roiw_pixel, scalerad_pixel, - weight_pixel, softrad_pixel, - self.cdelt3_normal, - roiw_ave, self.cdelt1, self.cdelt2) + result = cube_wrapper( + instrument, + flag_dq_plane, + weight_type, + start_region, + end_region, + self.overlap_partial, + self.overlap_full, + self.xcoord, + self.ycoord, + self.zcoord, + coord1, + coord2, + wave, + flux, + err, + slice_no, + rois_pixel, + roiw_pixel, + scalerad_pixel, + weight_pixel, + softrad_pixel, + self.cdelt3_normal, + roiw_ave, + self.cdelt1, + self.cdelt2, + ) spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux, spaxel_dq = result self.spaxel_flux = self.spaxel_flux + np.asarray(spaxel_flux, np.float64) - self.spaxel_weight = self.spaxel_weight + np.asarray(spaxel_weight, np.float64) + self.spaxel_weight = self.spaxel_weight + np.asarray( + spaxel_weight, np.float64 + ) self.spaxel_var = self.spaxel_var + np.asarray(spaxel_var, np.float64) self.spaxel_iflux = self.spaxel_iflux + np.asarray(spaxel_iflux, np.float64) spaxel_dq.astype(np.uint) @@ -674,7 +790,7 @@ def build_ifucube(self): result = None del result del spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux, spaxel_dq - if self.weighting == 'drizzle' and build_cube: + if self.weighting == "drizzle" and build_cube: cdelt3_mean = np.nanmean(self.cdelt3_normal) xi1, eta1, xi2, eta2, xi3, eta3, xi4, eta4 = corner_coord linear = 0 @@ -682,20 +798,46 @@ def build_ifucube(self): linear = 1 if debug_cube_index >= 0: log.info(f"Input filename: {input_model.meta.filename}") - result = cube_wrapper_driz(instrument, flag_dq_plane, - start_region, end_region, - self.overlap_partial, self.overlap_full, - self.xcoord, self.ycoord, self.zcoord, - coord1, coord2, wave, flux, err, slice_no, - xi1, eta1, xi2, eta2, xi3, eta3, xi4, eta4, - dwave, - self.cdelt3_normal, - self.cdelt1, self.cdelt2, cdelt3_mean, linear, - x_det, y_det, debug_cube_index) + result = cube_wrapper_driz( + instrument, + flag_dq_plane, + start_region, + end_region, + self.overlap_partial, + self.overlap_full, + self.xcoord, + self.ycoord, + self.zcoord, + coord1, + coord2, + wave, + flux, + err, + slice_no, + xi1, + eta1, + xi2, + eta2, + xi3, + eta3, + xi4, + eta4, + dwave, + self.cdelt3_normal, + self.cdelt1, + self.cdelt2, + cdelt3_mean, + linear, + x_det, + y_det, + debug_cube_index, + ) spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux, spaxel_dq = result self.spaxel_flux = self.spaxel_flux + np.asarray(spaxel_flux, np.float64) - self.spaxel_weight = self.spaxel_weight + np.asarray(spaxel_weight, np.float64) + self.spaxel_weight = self.spaxel_weight + np.asarray( + spaxel_weight, np.float64 + ) self.spaxel_var = self.spaxel_var + np.asarray(spaxel_var, np.float64) self.spaxel_iflux = self.spaxel_iflux + np.asarray(spaxel_iflux, np.float64) spaxel_dq.astype(np.uint) @@ -703,80 +845,110 @@ def build_ifucube(self): result = None del result del spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux, spaxel_dq - # -------------------------------------------------------------------------------- - # # AREA - 2d method only works for single files local slicer plane (internal_cal) - # -------------------------------------------------------------------------------- - elif self.interpolation == 'area': - # -------------------------------------------------------------------------------- + + # AREA - 2d method only works for single files local slicer plane (internal_cal) + elif self.interpolation == "area": + # ---------------------------------------------------------------------------- # MIRI - # -------------------------------------------------------------------------------- - if self.instrument == 'MIRI': - det2ab_transform = input_model.meta.wcs.get_transform('detector', - 'alpha_beta') - start_region = self.instrument_info.GetStartSlice(this_par1) - end_region = self.instrument_info.GetEndSlice(this_par1) + # ---------------------------------------------------------------------------- + if self.instrument == "MIRI": + det2ab_transform = input_model.meta.wcs.get_transform( + "detector", "alpha_beta" + ) + start_region = self.instrument_info.get_start_slice(this_par1) + end_region = self.instrument_info.get_end_slice(this_par1) regions = list(range(start_region, end_region + 1)) for i in regions: - log.info('Working on Slice # %d', i) + log.info("Working on Slice # %d", i) y, x = (det2ab_transform.label_mapper.mapper == i).nonzero() # getting pixel corner - ytop = y + 1 (routine fails for y = 1024) index = np.where(y < 1023) y = y[index] x = x[index] - slice = i - start_region - result = cube_internal_cal.match_det2cube(self.instrument, - x, y, slice, - input_model, - det2ab_transform, - self.xcoord, self.zcoord, - self.crval1, self.crval3, - self.cdelt1, self.cdelt3, - self.naxis1, self.naxis2) + slice_ifu = i - start_region + result = cube_internal_cal.match_det2cube( + self.instrument, + x, + y, + slice_ifu, + input_model, + det2ab_transform, + self.xcoord, + self.zcoord, + self.crval1, + self.crval3, + self.cdelt1, + self.cdelt3, + self.naxis1, + self.naxis2, + ) spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux = result - self.spaxel_flux = self.spaxel_flux + np.asarray(spaxel_flux, np.float64) - self.spaxel_weight = self.spaxel_weight + np.asarray(spaxel_weight, np.float64) + self.spaxel_flux = self.spaxel_flux + np.asarray( + spaxel_flux, np.float64 + ) + self.spaxel_weight = self.spaxel_weight + np.asarray( + spaxel_weight, np.float64 + ) self.spaxel_var = self.spaxel_var + np.asarray(spaxel_var, np.float64) - self.spaxel_iflux = self.spaxel_iflux + np.asarray(spaxel_iflux, np.float64) + self.spaxel_iflux = self.spaxel_iflux + np.asarray( + spaxel_iflux, np.float64 + ) result = None del spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux, result - # -------------------------------------------------------------------------------- + # ---------------------------------------------------------------------------- # NIRSPEC - # -------------------------------------------------------------------------------- - if self.instrument == 'NIRSPEC': + # ---------------------------------------------------------------------------- + if self.instrument == "NIRSPEC": nslices = 30 - slicemap = [15, 14, 16, 13, 17, 12, 18, 11, 19, 10, - 20, 9, 21, 8, 22, 7, 23, 6, 24, 5, 25, - 4, 26, 3, 27, 2, 28, 1, 29, 0] + slicemap = [15, 14, 16, 13, 17, 12, 18, 11, 19, 10, 20, 9, 21, 8, 22, 7, + 23, 6, 24, 5, 25, 4, 26, 3, 27, 2, 28, 1, 29, 0] # fmt: skip - wcsobj, tr1, tr2, tr3 = nirspec._get_transforms(input_model, np.arange(nslices)) + wcsobj, tr1, tr2, tr3 = nirspec._get_transforms( # noqa: SLF001 + input_model, np.arange(nslices) + ) for i in range(nslices): - - slice_wcs = nirspec._nrs_wcs_set_input_lite(input_model, wcsobj, i, - [tr1, tr2[i], tr3[i]]) - x, y = wcstools.grid_from_bounding_box(slice_wcs.bounding_box, step=(1, 1), center=True) - detector2slicer = slice_wcs.get_transform('detector', 'slicer') - - result = cube_internal_cal.match_det2cube(self.instrument, - x, y, slicemap[i], - input_model, - detector2slicer, - self.ycoord, self.zcoord, - self.crval2, self.crval3, - self.cdelt2, self.cdelt3, - self.naxis1, self.naxis2) + slice_wcs = nirspec._nrs_wcs_set_input_lite( # noqa: SLF001 + input_model, wcsobj, i, [tr1, tr2[i], tr3[i]] + ) + x, y = wcstools.grid_from_bounding_box( + slice_wcs.bounding_box, step=(1, 1), center=True + ) + detector2slicer = slice_wcs.get_transform("detector", "slicer") + + result = cube_internal_cal.match_det2cube( + self.instrument, + x, + y, + slicemap[i], + input_model, + detector2slicer, + self.ycoord, + self.zcoord, + self.crval2, + self.crval3, + self.cdelt2, + self.cdelt3, + self.naxis1, + self.naxis2, + ) spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux = result - self.spaxel_flux = self.spaxel_flux + np.asarray(spaxel_flux, np.float64) - self.spaxel_weight = self.spaxel_weight + np.asarray(spaxel_weight, np.float64) + self.spaxel_flux = self.spaxel_flux + np.asarray( + spaxel_flux, np.float64 + ) + self.spaxel_weight = self.spaxel_weight + np.asarray( + spaxel_weight, np.float64 + ) self.spaxel_var = self.spaxel_var + np.asarray(spaxel_var, np.float64) - self.spaxel_iflux = self.spaxel_iflux + np.asarray(spaxel_iflux, np.float64) + self.spaxel_iflux = self.spaxel_iflux + np.asarray( + spaxel_iflux, np.float64 + ) result = None del spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux, result k = k + 1 - # _______________________________________________________________________ # done looping over files self.find_spaxel_flux() @@ -786,21 +958,25 @@ def build_ifucube(self): result = self.setup_final_ifucube_model(input_model_ref) return result - # ******************************************************************************** + # ________________________________________________________________________________ def build_ifucube_single(self): - """ Build a set of single mode IFU cubes used for outlier detection - and background matching + """ + Build a set of single mode IFU cubes. Loop over every band contained in the IFU cube and read in the data associated with the band. Map each band to the output cube coordinate - system + system. + Returns + ------- + single_ifucube_container : IFUCubeModel + A single type IFU cube datamodel """ # loop over input models single_ifucube_container = ModelContainer() weight_type = 0 # default to emsm instead of msm - if self.weighting == 'msm': + if self.weighting == "msm": weight_type = 1 number_bands = len(self.list_par1) this_par1 = self.list_par1[0] # single IFUcube only have a single channel @@ -808,12 +984,12 @@ def build_ifucube_single(self): for i in range(number_bands): this_par2 = self.list_par2[i] nfiles = len(self.master_table.FileMap[self.instrument][this_par1][this_par2]) - # ________________________________________________________________________________ + # loop over the files that cover the spectral range the cube is for for k in range(nfiles): input_model = self.master_table.FileMap[self.instrument][this_par1][this_par2][k] self.input_models_this_cube.append(input_model) - log.debug("Working on next Single IFU Cube = %i" % (j + 1)) + log.debug(f"Working on next Single IFU Cube = {j + 1}") # for each new data model create a new spaxel total_num = self.naxis1 * self.naxis2 * self.naxis3 @@ -825,41 +1001,75 @@ def build_ifucube_single(self): subtract_background = False - pixelresult = self.map_detector_to_outputframe(this_par1, - subtract_background, - input_model) - - coord1, coord2, corner_coord, wave, dwave, flux, err, slice_no, \ - rois_pixel, roiw_pixel, weight_pixel, \ - softrad_pixel, scalerad_pixel, _, _ = pixelresult + pixelresult = self.map_detector_to_outputframe( + this_par1, subtract_background, input_model + ) + + ( + coord1, + coord2, + corner_coord, + wave, + dwave, + flux, + err, + slice_no, + rois_pixel, + roiw_pixel, + weight_pixel, + softrad_pixel, + scalerad_pixel, + _, + _, + ) = pixelresult build_cube = True - if wave is None: # there is no valid data on the detector. Pixels are flagged as DO_NOT_USE. + # there is no valid data on the detector. Pixels are flagged as DO_NOT_USE. + if wave is None: build_cube = False - # the following values are not needed in cube_wrapper because the DQ plane is not being - # filled in + # The following values are not needed in cube_wrapper because the DQ plane + # is not being filled in. flag_dq_plane = 0 start_region = 0 end_region = 0 roiw_ave = 0 - if self.instrument == 'MIRI': + if self.instrument == "MIRI": instrument = 0 else: # NIRSPEC instrument = 1 result = None - if self.interpolation == 'pointcloud' and build_cube: + if self.interpolation == "pointcloud" and build_cube: roiw_ave = np.mean(roiw_pixel) - result = cube_wrapper(instrument, flag_dq_plane, weight_type, start_region, end_region, - self.overlap_partial, self.overlap_full, - self.xcoord, self.ycoord, self.zcoord, - coord1, coord2, wave, flux, err, slice_no, - rois_pixel, roiw_pixel, scalerad_pixel, - weight_pixel, softrad_pixel, - self.cdelt3_normal, - roiw_ave, self.cdelt1, self.cdelt2) + result = cube_wrapper( + instrument, + flag_dq_plane, + weight_type, + start_region, + end_region, + self.overlap_partial, + self.overlap_full, + self.xcoord, + self.ycoord, + self.zcoord, + coord1, + coord2, + wave, + flux, + err, + slice_no, + rois_pixel, + roiw_pixel, + scalerad_pixel, + weight_pixel, + softrad_pixel, + self.cdelt3_normal, + roiw_ave, + self.cdelt1, + self.cdelt2, + ) spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux, _ = result self.spaxel_flux = self.spaxel_flux + np.asarray(spaxel_flux, np.float64) @@ -869,7 +1079,7 @@ def build_ifucube_single(self): result = None del result, spaxel_flux, spaxel_var, spaxel_iflux - if self.weighting == 'drizzle' and build_cube: + if self.weighting == "drizzle" and build_cube: cdelt3_mean = np.nanmean(self.cdelt3_normal) xi1, eta1, xi2, eta2, xi3, eta3, xi4, eta4 = corner_coord linear = 0 @@ -878,16 +1088,40 @@ def build_ifucube_single(self): x_det = None y_det = None debug_cube = -1 - result = cube_wrapper_driz(instrument, flag_dq_plane, - start_region, end_region, - self.overlap_partial, self.overlap_full, - self.xcoord, self.ycoord, self.zcoord, - coord1, coord2, wave, flux, err, slice_no, - xi1, eta1, xi2, eta2, xi3, eta3, xi4, eta4, - dwave, - self.cdelt3_normal, - self.cdelt1, self.cdelt2, cdelt3_mean, linear, - x_det, y_det, debug_cube) + result = cube_wrapper_driz( + instrument, + flag_dq_plane, + start_region, + end_region, + self.overlap_partial, + self.overlap_full, + self.xcoord, + self.ycoord, + self.zcoord, + coord1, + coord2, + wave, + flux, + err, + slice_no, + xi1, + eta1, + xi2, + eta2, + xi3, + eta3, + xi4, + eta4, + dwave, + self.cdelt3_normal, + self.cdelt1, + self.cdelt2, + cdelt3_mean, + linear, + x_det, + y_det, + debug_cube, + ) spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux, _ = result self.spaxel_flux = self.spaxel_flux + np.asarray(spaxel_flux, np.float64) @@ -896,7 +1130,7 @@ def build_ifucube_single(self): self.spaxel_iflux = self.spaxel_iflux + np.asarray(spaxel_iflux, np.float64) result = None del result, spaxel_flux, spaxel_var, spaxel_iflux - # ______________________________________________________________________ + # shove Flux and iflux in the final ifucube self.find_spaxel_flux() @@ -911,26 +1145,21 @@ def build_ifucube_single(self): j = j + 1 return single_ifucube_container - # ************************************************************************** + # ________________________________________________________________________________ def determine_cube_parameters_internal(self): - """Determine the spatial and spectral ifu size for coord_system = internal_cal - - """ - - # ____________________________________________________________ + """Determine the spatial and spectral IFU size for coord_system = internal_cal.""" # internal_cal is for only 1 file and weighting= area # no msm or emsm information is needed par1 = self.list_par1[0] par2 = self.list_par2[0] - a_scale, b_scale, w_scale = self.instrument_info.GetScale(par1, - par2) + a_scale, b_scale, w_scale = self.instrument_info.get_scale(par1, par2) self.spatial_size = a_scale if self.scalexy != 0: self.spatial_size = self.scalexy - min_wave = self.instrument_info.GetWaveMin(par1, par2) - max_wave = self.instrument_info.GetWaveMax(par1, par2) + min_wave = self.instrument_info.get_wave_min(par1, par2) + max_wave = self.instrument_info.get_wave_max(par1, par2) if self.wavemin is None: self.wavemin = min_wave else: @@ -950,20 +1179,13 @@ def determine_cube_parameters_internal(self): self.spectral_size = w_scale self.linear_wavelength = True -# ************************************************************************** - + # ________________________________________________________________________________ def determine_cube_parameters(self): - """Determine the spatial and wavelength roi size to use for - selecting point cloud elements around the spaxel centers. - - If the IFU cube covers more than 1 band - then use the rules to - define the Spatial and Wavelength roi size to use for the cube - Current Rule: using the minimum - - Returns - ------- - roi size for spatial and wavelength + """ + Determine the spatial and wavelength roi size if IFU covers more than 1 band of data. + If the IFU cube covers more than 1 band, then use the rules to + define the spatial and wavelength roi size to use for the cube. """ # initialize wave_roi = None @@ -980,29 +1202,27 @@ def determine_cube_parameters(self): minwave = np.zeros(number_bands) maxwave = np.zeros(number_bands) - # ____________________________________________________________ for i in range(number_bands): - if self.instrument == 'MIRI': + if self.instrument == "MIRI": par1 = self.list_par1[i] par2 = self.list_par2[i] - elif self.instrument == 'NIRSPEC': + elif self.instrument == "NIRSPEC": par1 = self.list_par1[i] par2 = self.list_par2[i] - a_scale, b_scale, w_scale = self.instrument_info.GetScale(par1, - par2) + a_scale, b_scale, w_scale = self.instrument_info.get_scale(par1, par2) spaxelsize[i] = a_scale spectralsize[i] = w_scale - minwave[i] = self.instrument_info.GetWaveMin(par1, par2) - maxwave[i] = self.instrument_info.GetWaveMax(par1, par2) + minwave[i] = self.instrument_info.get_wave_min(par1, par2) + maxwave[i] = self.instrument_info.get_wave_max(par1, par2) # pull out the values from the cube pars reference file - roiw[i] = self.instrument_info.GetWaveRoi(par1, par2) - rois[i] = self.instrument_info.GetSpatialRoi(par1, par2) - power[i] = self.instrument_info.GetMSMPower(par1, par2) - softrad[i] = self.instrument_info.GetSoftRad(par1, par2) - scalerad[i] = self.instrument_info.GetScaleRad(par1, par2) + roiw[i] = self.instrument_info.get_wave_roi(par1, par2) + rois[i] = self.instrument_info.get_spatial_roi(par1, par2) + power[i] = self.instrument_info.get_msm_power(par1, par2) + softrad[i] = self.instrument_info.get_soft_rad(par1, par2) + scalerad[i] = self.instrument_info.get_scale_rad(par1, par2) # Check the spatial size. If it is the same for the array set up the parameters all_same = np.all(spaxelsize == spaxelsize[0]) @@ -1056,31 +1276,40 @@ def determine_cube_parameters(self): self.scalerad = scalerad[0] else: self.linear_wavelength = False - if self.instrument == 'MIRI': - - table = self.instrument_info.Get_multichannel_table(self.weighting) - (table_wavelength, table_sroi, - table_wroi, table_power, - table_softrad, table_scalerad) = table + if self.instrument == "MIRI": + table = self.instrument_info.get_multichannel_table() + ( + table_wavelength, + table_sroi, + table_wroi, + table_power, + table_softrad, + table_scalerad, + ) = table # getting NIRSPEC Table Values - elif self.instrument == 'NIRSPEC': + elif self.instrument == "NIRSPEC": # determine if have Prism, Medium or High resolution - med = ['g140m', 'g235m', 'g395m'] - high = ['g140h', 'g235h', 'g395h'] - prism = ['prism'] + med = ["g140m", "g235m", "g395m"] + high = ["g140h", "g235h", "g395h"] + prism = ["prism"] for i in range(number_bands): par1 = self.list_par1[i] if par1 in prism: - table = self.instrument_info.Get_prism_table() + table = self.instrument_info.get_prism_table() if par1 in med: - table = self.instrument_info.Get_med_table() + table = self.instrument_info.get_med_table() if par1 in high: - table = self.instrument_info.Get_high_table() - (table_wavelength, table_sroi, - table_wroi, table_power, - table_softrad, table_scalerad) = table + table = self.instrument_info.get_high_table() + ( + table_wavelength, + table_sroi, + table_wroi, + table_power, + table_softrad, + table_scalerad, + ) = table # based on Min and Max wavelength - pull out the tables values that fall in this range # find the closest table entries to the self.wavemin and self.wavemax limits imin = (np.abs(table_wavelength - self.wavemin)).argmin() @@ -1088,30 +1317,29 @@ def determine_cube_parameters(self): if imin > 1 and table_wavelength[imin] > self.wavemin: imin = imin - 1 - if (imax < (len(table_wavelength) - 1) and - self.wavemax > table_wavelength[imax]): + if imax < (len(table_wavelength) - 1) and self.wavemax > table_wavelength[imax]: imax = imax + 1 num_table = imax - imin + 1 self.scalerad_table = np.zeros(num_table) - self.scalerad_table[:] = table_scalerad[imin:imax + 1] + self.scalerad_table[:] = table_scalerad[imin : imax + 1] self.softrad_table = np.zeros(num_table) - self.softrad_table[:] = table_softrad[imin:imax + 1] + self.softrad_table[:] = table_softrad[imin : imax + 1] self.roiw_table = np.zeros(num_table) - self.roiw_table[:] = table_wroi[imin:imax + 1] + self.roiw_table[:] = table_wroi[imin : imax + 1] self.rois_table = np.zeros(num_table) - self.rois_table[:] = table_sroi[imin:imax + 1] + self.rois_table[:] = table_sroi[imin : imax + 1] if self.num_files < 4: self.rois_table = self.rois_table * 1.5 self.weight_power_table = np.zeros(num_table) - self.weight_power_table[:] = table_power[imin:imax + 1] + self.weight_power_table[:] = table_power[imin : imax + 1] self.wavelength_table = np.zeros(num_table) - self.wavelength_table[:] = table_wavelength[imin:imax + 1] + self.wavelength_table[:] = table_wavelength[imin : imax + 1] # check if using default values from the table (not user set) @@ -1121,12 +1349,14 @@ def determine_cube_parameters(self): # default rois in tables is designed with a 4 dither pattern # increase rois if less than 4 file - if self.output_type == 'single' or self.num_files < 4: + if self.output_type == "single" or self.num_files < 4: # We don't need to increase it if using 'emsm' weighting - if self.weighting.lower() != 'emsm': + if self.weighting.lower() != "emsm": self.rois = self.rois * 1.5 - log.info('Increasing spatial region of interest ' - f'default value set for 4 dithers {self.rois}') + log.info( + "Increasing spatial region of interest " + f"default value set for 4 dithers {self.rois}" + ) # set wave_roi and weight_power to same values if they are in list if self.roiw == 0: @@ -1141,51 +1371,61 @@ def determine_cube_parameters(self): found_error = False if self.linear_wavelength: # check we have valid data for key values - if self.interpolation == 'pointcloud': + if self.interpolation == "pointcloud": if np.isnan(self.rois): - log.error('Spatial roi is nan, possible reference file value error') + log.error("Spatial roi is nan, possible reference file value error") found_error = True if np.isnan(self.roiw): - log.error('Spectral roi is nan, possible reference file value error') + log.error("Spectral roi is nan, possible reference file value error") found_error = True - if self.weighting == 'msm': + if self.weighting == "msm": if np.isnan(self.weight_power): - log.error('Weight power is nan, possible reference file value error') + log.error("Weight power is nan, possible reference file value error") found_error = True if np.isnan(self.soft_rad): - log.error('Soft rad is nan, possible reference file value error') + log.error("Soft rad is nan, possible reference file value error") found_error = True - if self.weighting == 'emsm': + if self.weighting == "emsm": if np.isnan(self.scalerad): - log.error('Scalerad is nan, possible reference file value error') + log.error("Scalerad is nan, possible reference file value error") found_error = True else: if np.isnan(self.wavelength_table).all(): - log.error('Wavelength table contains all nans, possible reference file value error') + log.error("Wavelength table contains all nans, possible reference file value error") found_error = True - if self.interpolation == 'pointcloud': + if self.interpolation == "pointcloud": if np.isnan(self.rois_table).all(): - log.error('Spatial roi table contains all nans, possible reference file value error') + log.error( + "Spatial roi table contains all nans, possible reference file value error" + ) found_error = True if np.isnan(self.roiw_table).all(): - log.error('Spectral roi table contains all nans, possible reference file value error') + log.error( + "Spectral roi table contains all nans, possible reference file value error" + ) found_error = True - if self.weighting == 'msm': + if self.weighting == "msm": if np.isnan(self.softrad_table).all(): - log.error('Soft rad table contains all nans, possible reference file value error') + log.error( + "Soft rad table contains all nans, possible reference file value error" + ) found_error = True if np.isnan(self.weight_power_table).all(): - log.error('Weight power table contains all nans, possible reference file value error') + log.error( + "Weight power table contains all nans, possible reference file value error" + ) found_error = True - if self.weighting == 'emsm': + if self.weighting == "emsm": if np.isnan(self.scalerad_table).all(): - log.error('Scalerad table contains all nans, possible reference file value error') + log.error( + "Scalerad table contains all nans, possible reference file value error" + ) found_error = True if found_error: - raise IncorrectParameter("An essential parameter is = nan, refer to apply error message") + raise IncorrectParameterError("An essential parameter is = nan, refer to error message") # catch where self.weight_power = nan weighting = msm written to header # TODO update writing to header scalerad if weighting = emsm @@ -1194,32 +1434,27 @@ def determine_cube_parameters(self): if np.isnan(self.weight_power): self.weight_power = None - log.debug(f'spatial size {self.spatial_size}') - log.debug(f'spectral size {self.spectral_size}') - log.debug(f'spatial roi {self.rois}') - log.debug(f'wave min and max {self.wavemin} {self.wavemax}') - log.debug(f'linear wavelength {self.linear_wavelength}') - log.debug(f'roiw {self.roiw}') - log.debug(f'output_type {self.output_type}') - if self.weighting == 'msm': - log.debug(f'weight_power {self.weight_power}') - log.debug(f'softrad {self.soft_rad}') - if self.weighting == 'emsm': - log.debug(f'scalerad {self.scalerad}') - -# ****************************************************************************** + log.debug(f"spatial size {self.spatial_size}") + log.debug(f"spectral size {self.spectral_size}") + log.debug(f"spatial roi {self.rois}") + log.debug(f"wave min and max {self.wavemin} {self.wavemax}") + log.debug(f"linear wavelength {self.linear_wavelength}") + log.debug(f"roiw {self.roiw}") + log.debug(f"output_type {self.output_type}") + if self.weighting == "msm": + log.debug(f"weight_power {self.weight_power}") + log.debug(f"softrad {self.soft_rad}") + if self.weighting == "emsm": + log.debug(f"scalerad {self.scalerad}") + # ________________________________________________________________________________ def setup_ifucube_wcs(self): - """Function to determine the min and max coordinates of the spectral - cube + """ + Set up the WCS of the IFU cube. Loop over every datamodel contained in the cube and find the WCS of the output cube that contains all the data. - Returns - ------- - Footprint of cube: min and max of coordinates of cube. - Notes ----- If the coordinate system is internal_cal then min and max @@ -1229,10 +1464,8 @@ def setup_ifucube_wcs(self): For NIRSPEC the units along/across slice dimension are meters If the coordinate system is skyalign/ifualign then the min and max of - ra(degrees), dec (degrees) and lambda (microns) is returned. - + RA(degrees), dec (degrees) and lambda (microns) is returned. """ -# _____________________________________________________________________________ self.cdelt1 = self.spatial_size self.cdelt2 = self.spatial_size if self.linear_wavelength: @@ -1240,105 +1473,109 @@ def setup_ifucube_wcs(self): parameter1 = self.list_par1 parameter2 = self.list_par2 -# ________________________________________________________________________________ -# 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 + # Define the rotation angle - if self.coord_system == 'ifualign': + # 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": this_a = parameter1[0] # 0 is first band - this_a is channel this_b = parameter2[0] # 0 is first band - this_b is sub-channel - log.info(f'Defining rotation between ra-dec and IFU plane using {this_a}, {this_b}') + log.info(f"Defining rotation between ra-dec and IFU plane using {this_a}, {this_b}") # first file for this band input_model = self.master_table.FileMap[self.instrument][this_a][this_b][0] - if self.instrument == 'MIRI': - xstart, xend = self.instrument_info.GetMIRISliceEndPts(this_a) + if self.instrument == "MIRI": + xstart, xend = self.instrument_info.get_miri_slice_endpts(this_a) ysize = input_model.data.shape[0] y, x = np.mgrid[:ysize, xstart:xend] - detector2alpha_beta = input_model.meta.wcs.get_transform('detector', - 'alpha_beta') + detector2alpha_beta = input_model.meta.wcs.get_transform("detector", "alpha_beta") alpha, beta, lam = detector2alpha_beta(x, y) lam_med = np.nanmedian(lam) # pick two alpha, beta values to determine rotation angle # values in arc seconds - alpha_beta2world = input_model.meta.wcs.get_transform('alpha_beta', - input_model.meta.wcs.output_frame.name) + alpha_beta2world = input_model.meta.wcs.get_transform( + "alpha_beta", input_model.meta.wcs.output_frame.name + ) temp_ra1, temp_dec1, lam_temp = alpha_beta2world(0, 0, lam_med) temp_ra2, temp_dec2, lam_temp = alpha_beta2world(2, 0, lam_med) - elif self.instrument == 'NIRSPEC': + elif self.instrument == "NIRSPEC": slice_wcs = nirspec.nrs_wcs_set_input(input_model, 0) - x, y = wcstools.grid_from_bounding_box(slice_wcs.bounding_box, step=(1, 1), center=True) - detector2slicer = slice_wcs.get_transform('detector', 'slicer') + x, y = wcstools.grid_from_bounding_box( + slice_wcs.bounding_box, step=(1, 1), center=True + ) + detector2slicer = slice_wcs.get_transform("detector", "slicer") across, along, lam = detector2slicer(x, y) # lam ~0 for this transform lam_med = np.nanmedian(lam) # pick two along slice, across slice values to determine rotation angle # values in meters - slicer2world = slice_wcs.get_transform('slicer', slice_wcs.output_frame.name) + slicer2world = slice_wcs.get_transform("slicer", slice_wcs.output_frame.name) temp_ra1, temp_dec1, lam_temp = slicer2world(0, 0, lam_med) temp_ra2, temp_dec2, lam_temp = slicer2world(0, 0.005, lam_med) - # ________________________________________________________________________________ + # temp_dec1 is in degrees - dra, ddec = (temp_ra2 - temp_ra1) * np.cos(temp_dec1 * np.pi / 180.0), (temp_dec2 - temp_dec1) - self.rot_angle = 90 + np.arctan2(dra, ddec) * 180. / np.pi - log.info(f'Rotation angle between ifu and sky: {self.rot_angle}') + dra, ddec = ( + (temp_ra2 - temp_ra1) * np.cos(temp_dec1 * np.pi / 180.0), + (temp_dec2 - temp_dec1), + ) + self.rot_angle = 90 + np.arctan2(dra, ddec) * 180.0 / 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 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: + 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 + 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 corner_a = [] corner_b = [] lambda_min = [] lambda_max = [] self.num_bands = len(self.list_par1) - log.debug('Number of bands in cube: %i', self.num_bands) + log.debug(f"Number of bands in cube: {self.num_bands}") for i in range(self.num_bands): this_a = parameter1[i] this_b = parameter2[i] - log.debug(f'Working on data from {this_a}, {this_b}') + log.debug(f"Working on data from {this_a}, {this_b}") n = len(self.master_table.FileMap[self.instrument][this_a][this_b]) - log.debug('number of files %d', n) + log.debug(f"number of files, {n}") for k in range(n): lmin = 0.0 lmax = 0.0 input_model = self.master_table.FileMap[self.instrument][this_a][this_b][k] -# ________________________________________________________________________________ # If offsets are provided. Pull in ra and dec offsets. raoffset = 0.0 decoffset = 0.0 # pull out ra dec offset if it exists if self.offsets is not None: raoffset, decoffset = self.find_ra_dec_offset(input_model.meta.filename) -# ________________________________________________________________________________ + # Find the footprint of the image - spectral_found = hasattr(input_model.meta.wcsinfo, 'spectral_region') - spatial_found = hasattr(input_model.meta.wcsinfo, 's_region') + 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' and self.cube_pa is None: + 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 # 3. Channel is 1 or 3 - channel with smaller FOV on detector - if self.instrument == 'MIRI' and self.output_type != 'multi': - ch1 = '1' - ch3 = '3' + if self.instrument == "MIRI" and self.output_type != "multi": + ch1 = "1" + ch3 = "3" if ch1 in self.list_par1 or ch3 in self.list_par1: spatial_found = False spectral_found = False @@ -1347,12 +1584,12 @@ def setup_ifucube_wcs(self): # The S_REGION values contain the footprint # on the sky of the original WCS. target_type = input_model.meta.target.type - if target_type == 'MOVING': + if target_type == "MOVING": spatial_found = False if spectral_found and spatial_found and world: [lmin, lmax] = input_model.meta.wcsinfo.spectral_region spatial_box = input_model.meta.wcsinfo.s_region - s = spatial_box.split(' ') + s = spatial_box.split(" ") ca1 = float(s[3]) cb1 = float(s[4]) ca2 = float(s[5]) @@ -1362,20 +1599,17 @@ def setup_ifucube_wcs(self): ca4 = float(s[9]) cb4 = float(s[10]) else: - log.info('Mapping all pixels to output to determine IFU foot print') + log.info("Mapping all pixels to output to determine IFU foot print") - if self.instrument == 'NIRSPEC': - ch_corners = cube_build_wcs_util.find_corners_NIRSPEC( - input_model, - self.instrument_info, - self.coord_system) + if self.instrument == "NIRSPEC": + ch_corners = cube_build_wcs_util.find_corners_nirspec( + input_model, self.coord_system + ) ca1, cb1, ca2, cb2, ca3, cb3, ca4, cb4, lmin, lmax = ch_corners - if self.instrument == 'MIRI': - ch_corners = cube_build_wcs_util.find_corners_MIRI( - input_model, - this_a, - self.instrument_info, - self.coord_system) + if self.instrument == "MIRI": + ch_corners = cube_build_wcs_util.find_corners_miri( + input_model, this_a, self.instrument_info, self.coord_system + ) ca1, cb1, ca2, cb2, ca3, cb3, ca4, cb4, lmin, lmax = ch_corners @@ -1398,7 +1632,7 @@ def setup_ifucube_wcs(self): lambda_min.append(lmin) 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) @@ -1408,64 +1642,66 @@ def setup_ifucube_wcs(self): 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.debug(f' min and max wavelengths: {min(lambda_min), max(lambda_max)}') - # ______________________________________________________________________ + log.debug(f"final a and b:{final_a_min, final_b_min, final_a_max, final_b_max}") + log.debug(f" min and max wavelengths: {min(lambda_min), max(lambda_max)}") + # the wavelength limits of cube are determined from 1. User or 2. cubepars # reference file (in the priority) final_lambda_min = self.wavemin final_lambda_max = self.wavemax - log.debug(f' final min and max used in IFUcube: {final_lambda_min, final_lambda_max}') - - if self.instrument == 'MIRI' and self.coord_system == 'internal_cal': + log.debug(f" final min and max used in IFUcube: {final_lambda_min, final_lambda_max}") + # ______________________________________________________________________ + if self.instrument == "MIRI" and self.coord_system == "internal_cal": # we have a 1 to 1 mapping in y across slice dimension - nslice = self.instrument_info.GetNSlice(parameter1[0]) - log.info(f'Across slice scale {self.cdelt2}') + nslice = self.instrument_info.get_nslice(parameter1[0]) + log.info(f"Across slice scale {self.cdelt2}") self.cdelt2 = (final_b_max - final_b_min) / nslice final_b_max = final_b_min + (nslice) * self.cdelt2 - log.info('Changed the across slice scale dimension so we have 1-1 mapping between b and slice #') - log.info(f'New across slice scale {self.cdelt2}') + log.info( + "Changed the across slice scale dimension so we have 1-1" + " mapping between b and slice #" + ) + log.info(f"New across slice scale {self.cdelt2}") # ______________________________________________________________________ - if self.instrument == 'NIRSPEC' and self.coord_system == 'internal_cal': + if self.instrument == "NIRSPEC" and self.coord_system == "internal_cal": # we have a 1 to 1 mapping in x - across slice dimension. nslice = 30 - log.info(f'Across slice scale {self.cdelt1}') + log.info(f"Across slice scale {self.cdelt1}") self.cdelt1 = (final_b_max - final_b_min) / nslice final_b_max = final_b_min + (nslice) * self.cdelt1 - log.info('Changed the across slice scale dimension so we have 1-1 mapping between b and slice #') - log.info(f'New across slice Scale {self.cdelt1}') + log.info( + "Changed the across slice scale dimension so we have 1-1" + " mapping between b and slice #" + ) + log.info(f"New across slice Scale {self.cdelt1}") self.cdelt2 = self.cdelt1 / 2.0 - # ________________________________________________________________________________ # Test that we have data (NIRSPEC NRS2 only has IFU data for 3 configurations) test_a = final_a_max - final_a_min test_b = final_b_max - final_b_min tolerance1 = 0.00001 if test_a < tolerance1 or test_b < tolerance1: - log.info(f'No Valid IFU slice data found {test_a} {test_b}') - # ________________________________________________________________________________ + log.info(f"No Valid IFU slice data found {test_a} {test_b}") + # Based on Scaling and Min and Max values determine naxis1, naxis2, naxis3 # set cube CRVALs, CRPIXs - if self.coord_system == 'skyalign' or self.coord_system == 'ifualign': + if self.coord_system == "skyalign" or self.coord_system == "ifualign": self.set_geometry(corner_a, corner_b, final_lambda_min, final_lambda_max) else: - self.set_geometryAB(corner_a, corner_b, final_lambda_min, final_lambda_max) + self.set_geometry_slicer(corner_a, corner_b, final_lambda_min, final_lambda_max) self.print_cube_geometry() - # ************************************************************************** - def map_detector_to_outputframe(self, this_par1, - subtract_background, - input_model): - - """Loop over a file and map the detector pixels to the output cube + # ________________________________________________________________________________ + def map_detector_to_outputframe(self, this_par1, subtract_background, input_model): + """ + Loop over a file and map the detector pixels to the output cube. The output frame is on the SKY (ra-dec) - Return the coordinates of all the detector pixel in the output frame. In addition, an array of pixel fluxes and weighing parameters are determined. The pixel flux and weighing parameters are used later in @@ -1476,34 +1712,34 @@ def map_detector_to_outputframe(self, this_par1, Parameters ---------- this_par1 : str - for MIRI this is the channel # for NIRSPEC this is the grating name + For MIRI this is the channel number (1,2,3 or 4). For NIRSPEC this is the grating name. only need for MIRI to distinguish which channel on the detector we have - subtract_background : boolean - if TRUE then subtract the background found in the mrs_imatch step only + subtract_background : bool + If TRUE then subtract the background found in the mrs_imatch step. Only needed for MIRI data - input: datamodel - input data model + input_model : IFUImageModel + Input IFU image model to combine Returns ------- - coord1 : numpy.ndarray - coordinate for axis1 in output cube for mapped pixel + coord1 : ndarray + Coordinate for axis1 in output cube for mapped pixel coord2: numpy.ndarray - coordinate for axis2 in output cube for mapped pixel - wave: numpy.ndarray - wavelength associated with coord1,coord2 - flux: numpy.ndarray - flux associated with coord1, coord2 - err: numpy.ndarray - err associated with coord1, coord2 + Coordinate for axis2 in output cube for mapped pixel + wave: ndarray + Wavelength associated with coord1,coord2 + flux: ndarray + Flux associated with coord1, coord2 + err: ndarray + Err associated with coord1, coord2 rois_det: float - spatial roi size to use - roiw_det: numpy.ndarray - spectral roi size associated with coord1,coord2 - weight_det : numpy.ndarray - weighting parameter association with coord1,coord2 - softrad_det : numpy.ndarray - weighting parameter association with coord1,coord2 + Spatial roi size to use + roiw_det: ndarray + Spectral roi size associated with coord1,coord2 + weight_det : ndarray + Weighting parameter associated with coord1,coord2 + softrad_det : ndarray + Weighting parameter associated with coord1,coord2 """ # initialize alpha_det and beta_det to None. These are filled in # if the instrument is MIRI and the weighting is miripsf @@ -1531,12 +1767,13 @@ def map_detector_to_outputframe(self, this_par1, y_det = None offsets = self.offsets - if self.instrument == 'MIRI': - sky_result = self.map_miri_pixel_to_sky(input_model, this_par1, subtract_background, - offsets) + if self.instrument == "MIRI": + sky_result = self.map_miri_pixel_to_sky( + input_model, this_par1, subtract_background, offsets + ) (x, y, ra, dec, wave_all, slice_no_all, dwave_all, corner_coord_all) = sky_result - elif self.instrument == 'NIRSPEC': + elif self.instrument == "NIRSPEC": sky_result = self.map_nirspec_pixel_to_sky(input_model, offsets) (x, y, ra, dec, wave_all, slice_no_all, dwave_all, corner_coord_all) = sky_result @@ -1562,7 +1799,7 @@ def map_detector_to_outputframe(self, this_par1, else: min_wave_tolerance = self.zcoord[0] - np.max(self.roiw_table) max_wave_tolerance = self.zcoord[-1] + np.max(self.roiw_table) - if self.interpolation == 'drizzle': + if self.interpolation == "drizzle": dmax = np.nanmax(dwave_all) if self.linear_wavelength: @@ -1576,32 +1813,48 @@ def map_detector_to_outputframe(self, this_par1, valid_max = np.where(wave_all <= max_wave_tolerance) not_mapped_high = wave_all.size - len(valid_max[0]) if not_mapped_low > 0: - log.info('# of detector pixels not mapped to output plane: ' - f'{not_mapped_low} with wavelength below {min_wave_tolerance}') + log.info( + "# of detector pixels not mapped to output plane: " + f"{not_mapped_low} with wavelength below {min_wave_tolerance}" + ) if not_mapped_high > 0: - log.info('# of detector pixels not mapped to output plane: ' - f'{not_mapped_high} with wavelength above {max_wave_tolerance}') + log.info( + "# of detector pixels not mapped to output plane: " + f"{not_mapped_high} with wavelength above {max_wave_tolerance}" + ) - # ______________________________________________________________________________ # using the DQFlags from the input_image find pixels that should be excluded # from the cube mapping - valid3 = np.logical_and(wave_all >= min_wave_tolerance, - wave_all <= max_wave_tolerance) + valid3 = np.logical_and(wave_all >= min_wave_tolerance, wave_all <= max_wave_tolerance) # find the location of good data - bad1 = np.bitwise_and(dq_all, dqflags.pixel['DO_NOT_USE']).astype(bool) - bad2 = np.bitwise_and(dq_all, dqflags.pixel['NON_SCIENCE']).astype(bool) + bad1 = np.bitwise_and(dq_all, dqflags.pixel["DO_NOT_USE"]).astype(bool) + bad2 = np.bitwise_and(dq_all, dqflags.pixel["NON_SCIENCE"]).astype(bool) good_data = np.where(~bad1 & ~bad2 & valid2 & valid3) num_good = len(good_data[0]) if num_good == 0: # This can occur if all the pixels on the detector are marked DO_NOT_USE. - log.warning(f'No valid pixels found on detector {input_model.meta.filename}') - return coord1, coord2, corner_coord, wave, dwave, flux, err, \ - slice_no, rois_det, roiw_det, weight_det, \ - softrad_det, scalerad_det, x_det, y_det + log.warning(f"No valid pixels found on detector {input_model.meta.filename}") + return ( + coord1, + coord2, + corner_coord, + wave, + dwave, + flux, + err, + slice_no, + rois_det, + roiw_det, + weight_det, + softrad_det, + scalerad_det, + x_det, + y_det, + ) # good data holds the location of pixels we want to map to cube # define variables as numpy arrays (numba needs this defined) @@ -1621,7 +1874,7 @@ def map_detector_to_outputframe(self, this_par1, x_det = x_all[good_data] y_det = y_all[good_data] - log.debug(f'After removing pixels based on criteria min and max wave: {np.min(wave)}, {np.max(wave)}') + log.debug(f"After removing pixels min and max wave: {np.min(wave)} {np.max(wave)}") # based on the wavelength define the sroi, wroi, weight_power and # softrad to use in matching detector to spaxel values @@ -1633,7 +1886,7 @@ def map_detector_to_outputframe(self, this_par1, # ________________________________________________________________________ # if working with msm or emsm need roi sizes and other parameters defined: - if self.weighting == 'msm' or self.weighting == 'emsm': + if self.weighting == "msm" or self.weighting == "emsm": if self.linear_wavelength: rois_det[:] = self.rois roiw_det[:] = self.roiw @@ -1642,29 +1895,28 @@ def map_detector_to_outputframe(self, this_par1, scalerad_det[:] = self.scalerad else: # for each wavelength find the closest point in the self.wavelength_table - for iw, w in enumerate(wave): - self.find_closest_wave(iw, w, - self.wavelength_table, - self.rois_table, - self.roiw_table, - self.softrad_table, - self.weight_power_table, - self.scalerad_table, - rois_det, - roiw_det, - softrad_det, - weight_det, - scalerad_det) - # ________________________________________________________________________ + self.find_closest_wave( + iw, + w, + self.wavelength_table, + self.rois_table, + self.roiw_table, + self.softrad_table, + self.weight_power_table, + self.scalerad_table, + rois_det, + roiw_det, + softrad_det, + weight_det, + scalerad_det, + ) + ra_use = ra[good_data] dec_use = dec[good_data] - coord1, coord2 = coord.radec2std(self.crval1, - self.crval2, - ra_use, dec_use, - self.rot_angle) + coord1, coord2 = coord.radec2std(self.crval1, self.crval2, ra_use, dec_use, self.rot_angle) - if self.interpolation == 'drizzle': + if self.interpolation == "drizzle": dwave = np.zeros(good_shape) dwave[:] = dwave_all[good_data] ra1 = corner_coord_all[0] @@ -1684,54 +1936,62 @@ def map_detector_to_outputframe(self, this_par1, ra4 = ra4[good_data] dec4 = dec4[good_data] - xi1, eta1 = coord.radec2std(self.crval1, - self.crval2, - ra1, dec1, - self.rot_angle) - - xi2, eta2 = coord.radec2std(self.crval1, - self.crval2, - ra2, dec2, - self.rot_angle) - xi3, eta3 = coord.radec2std(self.crval1, - self.crval2, - ra3, dec3, - self.rot_angle) - xi4, eta4 = coord.radec2std(self.crval1, - self.crval2, - ra4, dec4, - self.rot_angle) + xi1, eta1 = coord.radec2std(self.crval1, self.crval2, ra1, dec1, self.rot_angle) + + xi2, eta2 = coord.radec2std(self.crval1, self.crval2, ra2, dec2, self.rot_angle) + xi3, eta3 = coord.radec2std(self.crval1, self.crval2, ra3, dec3, self.rot_angle) + xi4, eta4 = coord.radec2std(self.crval1, self.crval2, ra4, dec4, self.rot_angle) corner_coord = [xi1, eta1, xi2, eta2, xi3, eta3, xi4, eta4] - return coord1, coord2, corner_coord, wave, dwave, flux, err, \ - slice_no, rois_det, roiw_det, weight_det, \ - softrad_det, scalerad_det, x_det, y_det + return ( + coord1, + coord2, + corner_coord, + wave, + dwave, + flux, + err, + slice_no, + rois_det, + roiw_det, + weight_det, + softrad_det, + scalerad_det, + x_det, + y_det, + ) + # ______________________________________________________________________ + def map_miri_pixel_to_sky(self, input_model, this_par1, subtract_background, offsets): + """ + Loop over a MIRI model and map the detector pixels to the output cube. - def map_miri_pixel_to_sky(self, input_model, this_par1, subtract_background, - offsets): - """Loop over a file and map the detector pixels to the output cube The output frame is on the SKY (ra-dec) - Return the coordinates of all the detector pixel in the output frame. Parameters ---------- + input_model : IFUImageModel + Input IFU image model to combine this_par1 : str - for MIRI this is the channel # for NIRSPEC this is the grating name + For MIRI this is the channel # for NIRSPEC this is the grating name only need for MIRI to distinguish which channel on the detector we have - subtract_background : boolean - if TRUE then subtract the background found in the mrs_imatch step only + subtract_background : bool + If TRUE then subtract the background found in the mrs_imatch step only needed for MIRI data - input: datamodel - input data model - offsets: dictionary - optional dictionary of ra and dec offsets to apply + offsets : dict + Optional dictionary of ra and dec offsets to apply Returns ------- - x, y, ra, dec, lambda, slice_no of valid slice pixels - + sky_result : tuple + For every valid input pixel from the IFU image model it contains + x,y: the pixel values on the detector + ra, dec: detector values mapped to sky + wave: wavelength corresponding to pixel + slice_no: slice_no of the pixel + dwave: delta wavelength covered by pixel + corner_coord: the corners of the pixel mapped to ra,dec """ wave = None slice_no = None # Slice number @@ -1743,14 +2003,22 @@ def map_miri_pixel_to_sky(self, input_model, this_par1, subtract_background, # pull out ra dec offset if it exists if offsets is not None: raoffset, decoffset = self.find_ra_dec_offset(input_model.meta.filename) - log.info("Ra and Dec offset (arc seconds) applied to file :%8.6f, %8.6f, %s", - raoffset.value, decoffset.value, input_model.meta.filename) + log.info( + "RA and Dec offset (arc seconds) applied to file :%8.6f, %8.6f, %s", + raoffset.value, + decoffset.value, + input_model.meta.filename, + ) # check if background sky matching as been done in mrs_imatch step # If it has not been subtracted and the background has not been # subtracted - subtract it. num_ch_bgk = len(input_model.meta.background.polynomial_info) - if num_ch_bgk > 0 and subtract_background and input_model.meta.background.subtracted is False: + if ( + num_ch_bgk > 0 + and subtract_background + and input_model.meta.background.subtracted is False + ): for ich_num in range(num_ch_bgk): poly = input_model.meta.background.polynomial_info[ich_num] poly_ch = poly.channel @@ -1760,10 +2028,9 @@ def map_miri_pixel_to_sky(self, input_model, this_par1, subtract_background, # find the slice number of each pixel and fill in slice_det ysize, xsize = input_model.data.shape slice_det = np.zeros((ysize, xsize), dtype=int) - det2ab_transform = input_model.meta.wcs.get_transform('detector', - 'alpha_beta') - start_region = self.instrument_info.GetStartSlice(this_par1) - end_region = self.instrument_info.GetEndSlice(this_par1) + det2ab_transform = input_model.meta.wcs.get_transform("detector", "alpha_beta") + start_region = self.instrument_info.get_start_slice(this_par1) + end_region = self.instrument_info.get_end_slice(this_par1) regions = list(range(start_region, end_region + 1)) for i in regions: ys, xs = (det2ab_transform.label_mapper.mapper == i).nonzero() @@ -1774,7 +2041,7 @@ def map_miri_pixel_to_sky(self, input_model, this_par1, subtract_background, slice_det[yind, xind] = i # define the x,y detector values of channel to be mapped to desired coordinate system - xstart, xend = self.instrument_info.GetMIRISliceEndPts(this_par1) + xstart, xend = self.instrument_info.get_miri_slice_endpts(this_par1) y, x = np.mgrid[:ysize, xstart:xend] y = np.reshape(y, y.size) x = np.reshape(x, x.size) @@ -1799,7 +2066,7 @@ def map_miri_pixel_to_sky(self, input_model, this_par1, subtract_background, yind = np.ndarray.flatten(yind) slice_no = slice_det[yind, xind] - if self.interpolation == 'drizzle': + if self.interpolation == "drizzle": # Delta wavelengths _, _, wave1 = input_model.meta.wcs(x, y - 0.4999) _, _, wave2 = input_model.meta.wcs(x, y + 0.4999) @@ -1807,23 +2074,43 @@ def map_miri_pixel_to_sky(self, input_model, this_par1, subtract_background, # Pixel corners pixfrac = 1.0 - alpha1, beta, _ = input_model.meta.wcs.transform('detector', 'alpha_beta', x - 0.4999 * pixfrac, y) - alpha2, _, _ = input_model.meta.wcs.transform('detector', 'alpha_beta', x + 0.4999 * pixfrac, y) + alpha1, beta, _ = input_model.meta.wcs.transform( + "detector", "alpha_beta", x - 0.4999 * pixfrac, y + ) + alpha2, _, _ = input_model.meta.wcs.transform( + "detector", "alpha_beta", x + 0.4999 * pixfrac, y + ) # Find slice width allbetaval = np.unique(beta) dbeta = np.abs(allbetaval[1] - allbetaval[0]) - ra1, dec1, _ = input_model.meta.wcs.transform('alpha_beta', - input_model.meta.wcs.output_frame, alpha1, - beta - dbeta * pixfrac / 2., wave) - ra2, dec2, _ = input_model.meta.wcs.transform('alpha_beta', - input_model.meta.wcs.output_frame, alpha1, - beta + dbeta * pixfrac / 2., wave) - ra3, dec3, _ = input_model.meta.wcs.transform('alpha_beta', - input_model.meta.wcs.output_frame, alpha2, - beta + dbeta * pixfrac / 2., wave) - ra4, dec4, _ = input_model.meta.wcs.transform('alpha_beta', - input_model.meta.wcs.output_frame, alpha2, - beta - dbeta * pixfrac / 2., wave) + ra1, dec1, _ = input_model.meta.wcs.transform( + "alpha_beta", + input_model.meta.wcs.output_frame, + alpha1, + beta - dbeta * pixfrac / 2.0, + wave, + ) + ra2, dec2, _ = input_model.meta.wcs.transform( + "alpha_beta", + input_model.meta.wcs.output_frame, + alpha1, + beta + dbeta * pixfrac / 2.0, + wave, + ) + ra3, dec3, _ = input_model.meta.wcs.transform( + "alpha_beta", + input_model.meta.wcs.output_frame, + alpha2, + beta + dbeta * pixfrac / 2.0, + wave, + ) + ra4, dec4, _ = input_model.meta.wcs.transform( + "alpha_beta", + input_model.meta.wcs.output_frame, + alpha2, + beta - dbeta * pixfrac / 2.0, + wave, + ) # now offset the pixel corners if offsets is not None: @@ -1839,31 +2126,42 @@ def map_miri_pixel_to_sky(self, input_model, this_par1, subtract_background, # ______________________________________________________________________ def map_nirspec_pixel_to_sky(self, input_model, offsets): - - """Loop over a file and map the detector pixels to the output cube + """ + Loop over a NIRSpec model and map the detector pixels to the output cube. The output frame is on the SKY (ra-dec) Return the coordinates of all the detector pixel in the output frame. Parameters ---------- - input: datamodel - input data model + input_model : IFUImageModel + Input IFU image model to combine + offsets : numpy array of floats + RA and Dec offsets to apply to each file Returns ------- - x, y, ra, dec, lambda, slice_no - + sky_result : tuple + For every valid input pixel from the IFU image model it contains + x,y: the pixel values on the detector + ra, dec: detector values mapped to sky + wave: wavelength corresponding to pixel + slice_no: slice_no of the pixel + dwave: delta wavelength covered by pixel + corner_coord: the corners of the pixel mapped to ra,dec """ - # check if we have an ra and dec offset file raoffset = 0.0 decoffset = 0.0 # pull out ra dec offset if it exists if offsets is not None: raoffset, decoffset = self.find_ra_dec_offset(input_model.meta.filename) - log.info("Ra and Dec offset (arc seconds) applied to file :%8.6f, %8.6f, %s", - raoffset.value, decoffset.value, input_model.meta.filename) + log.info( + "RA and Dec offset (arc seconds) applied to file :%8.6f, %8.6f, %s", + raoffset.value, + decoffset.value, + input_model.meta.filename, + ) # initialize the ra,dec, and wavelength arrays # we will loop over slice_nos and fill in values @@ -1890,14 +2188,14 @@ def map_nirspec_pixel_to_sky(self, input_model, offsets): # determine the slice width using slice 1 and 3 slice_wcs1 = nirspec.nrs_wcs_set_input(input_model, 0) - detector2slicer = slice_wcs1.get_transform('detector', 'slicer') + detector2slicer = slice_wcs1.get_transform("detector", "slicer") x, y = wcstools.grid_from_bounding_box(slice_wcs1.bounding_box) across1, along1, _ = detector2slicer(x, y - 0.4999 * pixfrac) across1 = across1[~np.isnan(across1)] slice_loc1 = np.unique(across1) slice_wcs3 = nirspec.nrs_wcs_set_input(input_model, 2) - detector2slicer = slice_wcs3.get_transform('detector', 'slicer') + detector2slicer = slice_wcs3.get_transform("detector", "slicer") x, y = wcstools.grid_from_bounding_box(slice_wcs3.bounding_box) across3, along3, _ = detector2slicer(x, y - 0.4999 * pixfrac) across3 = across3[~np.isnan(across3)] @@ -1909,11 +2207,12 @@ def map_nirspec_pixel_to_sky(self, input_model, offsets): nslices = 30 log.info("Mapping each NIRSpec slice to sky for input file: %s", input_model.meta.filename) - wcsobj, tr1, tr2, tr3 = nirspec._get_transforms(input_model, np.arange(nslices)) + wcsobj, tr1, tr2, tr3 = nirspec._get_transforms(input_model, np.arange(nslices)) # noqa: SLF001 for ii in range(nslices): - slice_wcs = nirspec._nrs_wcs_set_input_lite(input_model, wcsobj, ii, - [tr1, tr2[ii], tr3[ii]]) + slice_wcs = nirspec._nrs_wcs_set_input_lite( # noqa: SLF001 + input_model, wcsobj, ii, [tr1, tr2[ii], tr3[ii]] + ) x, y = wcstools.grid_from_bounding_box(slice_wcs.bounding_box) ra, dec, lam = slice_wcs(x, y) @@ -1926,7 +2225,7 @@ def map_nirspec_pixel_to_sky(self, input_model, offsets): dec = dec[valid] lam = lam[valid] - if self.interpolation == 'drizzle': + if self.interpolation == "drizzle": # Delta wavelengths _, _, wave1 = slice_wcs(x - 0.4999, y) _, _, wave2 = slice_wcs(x + 0.4999, y) @@ -1934,8 +2233,8 @@ def map_nirspec_pixel_to_sky(self, input_model, offsets): # Pixel corners pixfrac = 1.0 - detector2slicer = slice_wcs.get_transform('detector', 'slicer') - slicer2world = slice_wcs.get_transform('slicer', slice_wcs.output_frame) + detector2slicer = slice_wcs.get_transform("detector", "slicer") + slicer2world = slice_wcs.get_transform("slicer", slice_wcs.output_frame) across1, along1, lam1 = detector2slicer(x, y - 0.49 * pixfrac) across2, along2, lam2 = detector2slicer(x, y + 0.49 * pixfrac) @@ -1993,7 +2292,7 @@ def map_nirspec_pixel_to_sky(self, input_model, offsets): dec3_det[yind, xind] = dec3 dec4_det[yind, xind] = dec4 - else: # not drizzling + else: # not drizzling xind = _toindex(x) yind = _toindex(y) xind = np.ndarray.flatten(xind) @@ -2012,7 +2311,6 @@ def map_nirspec_pixel_to_sky(self, input_model, offsets): valid_data = np.where(flag_det == 1) y, x = valid_data - wave = lam_det[valid_data] slice_no = slice_det[valid_data] dwave = dwave_det[valid_data] @@ -2042,22 +2340,54 @@ def map_nirspec_pixel_to_sky(self, input_model, offsets): sky_result = (x, y, ra, dec, wave, slice_no, dwave, corner_coord) return sky_result - # ******************************************************************************** - def find_closest_wave(self, iw, w, - wavelength_table, - rois_table, - roiw_table, - softrad_table, - weight_power_table, - scalerad_table, - rois_det, - roiw_det, - softrad_det, - weight_det, - scalerad_det): - - """ Given a specific wavelength, find the closest value in the wavelength_table + # ________________________________________________________________________________ + def find_closest_wave( + self, + iw, + w, + wavelength_table, + rois_table, + roiw_table, + softrad_table, + weight_power_table, + scalerad_table, + rois_det, + roiw_det, + softrad_det, + weight_det, + scalerad_det, + ): + """ + Given a specific wavelength, find the closest value in the wavelength_table. + Parameters + ---------- + iw : int + Index of wavelength array. + w : float + Wavelength array of data. + wavelength_table : ndarray + Wavelength array read from cubepars reference file. + rois_table : ndarray + Rois array read from cubepars reference file. + roiw_table : ndarray + Roiw array read from cubepars reference file. + softrad_table : ndarray + Softrad array read from cubepars reference file. + weight_power_table : ndarray + Weight power array read from cubepars reference file. + scalerad_table : ndarray + Scalerad array read from cubepars reference file. + rois_det : ndarray + Rois array of detector pixel for the associated wavelength of the pixel. + roiw_det : ndarray + Roiw array of detector pixel for the associated wavelength of the pixel. + softrad_det : ndarray + Softrad array of detector pixel for the associated wavelength of the pixel. + weight_det : ndarray + Weight array of detector pixel for the associated wavelength of the pixel. + scalerad_det : ndarray + Scalerad array of detector pixel for the associated wavelength of the pixel. """ ifound = (np.abs(wavelength_table - w)).argmin() rois_det[iw] = rois_table[ifound] @@ -2066,32 +2396,32 @@ def find_closest_wave(self, iw, w, weight_det[iw] = weight_power_table[ifound] scalerad_det[iw] = scalerad_table[ifound] - # ******************************************************************************** + # ________________________________________________________________________________ def find_spaxel_flux(self): - - """Depending on the interpolation method, find the flux for each spaxel value - """ + """Depending on the interpolation method, find the flux for each spaxel value.""" # currently these are the same but in the future there could be a difference in # how the spaxel flux is determined according to self.interpolation. - if self.interpolation == 'area': + if self.interpolation == "area": good = self.spaxel_iflux > 0 self.spaxel_flux[good] = self.spaxel_flux[good] / self.spaxel_weight[good] - self.spaxel_var[good] = self.spaxel_var[good] / (self.spaxel_weight[good] * self.spaxel_weight[good]) - elif self.interpolation == 'pointcloud' or self.interpolation == 'drizzle': - # Don't apply any normalization if no points contributed to a spaxel (i.e., don't divide by zero) + self.spaxel_var[good] = self.spaxel_var[good] / ( + self.spaxel_weight[good] * self.spaxel_weight[good] + ) + elif self.interpolation == "pointcloud" or self.interpolation == "drizzle": + # Don't apply any normalization if no points contributed to a spaxel + # (i.e., don't divide by zero) good = self.spaxel_iflux > 0 # Normalize the weighted sum of pixel fluxes by the sum of the weights self.spaxel_flux[good] = self.spaxel_flux[good] / self.spaxel_weight[good] # Normalize the variance by the square of the weights - self.spaxel_var[good] = self.spaxel_var[good] / (self.spaxel_weight[good] * self.spaxel_weight[good]) + self.spaxel_var[good] = self.spaxel_var[good] / ( + self.spaxel_weight[good] * self.spaxel_weight[good] + ) - # ******************************************************************************** + # ________________________________________________________________________________ def set_final_dq_flags(self): - - """ Set up the final dq flags, Good data(0) , NON_SCIENCE or DO_NOT_USE - """ - + """Set up the final dq flags, Good data(0) , NON_SCIENCE or DO_NOT_USE.""" # An initial set of dq flags was set in overlap_fov_with_spaxel or # overlap_slice_with_spaxel. The initial dq dlags are defined in ifu_cube # class: @@ -2115,8 +2445,9 @@ def set_final_dq_flags(self): # convert all remaining spaxel_dq of 0 to NON_SCIENCE + DO_NOT_USE # these pixel should have no overlap with the data non_science = self.spaxel_dq == 0 - self.spaxel_dq[non_science] = np.bitwise_or(self.overlap_no_coverage, - dqflags.pixel['DO_NOT_USE']) + self.spaxel_dq[non_science] = np.bitwise_or( + self.overlap_no_coverage, dqflags.pixel["DO_NOT_USE"] + ) # refine where good data should be ind_full = np.where(np.bitwise_and(self.spaxel_dq, self.overlap_full)) @@ -2143,10 +2474,10 @@ def set_final_dq_flags(self): found = 0 ij = 0 # do not allow holes to occur at the edge of IFU cube - if (yrem == 0 or yrem == (self.naxis2 - 1) or - xrem == 0 or xrem == (self.naxis1 - 1)): - spaxel_dq_temp[index[0][i]] = np.bitwise_or(self.overlap_no_coverage, - dqflags.pixel['DO_NOT_USE']) + if yrem == 0 or yrem == (self.naxis2 - 1) or xrem == 0 or xrem == (self.naxis1 - 1): + spaxel_dq_temp[index[0][i]] = np.bitwise_or( + self.overlap_no_coverage, dqflags.pixel["DO_NOT_USE"] + ) found = 1 # flag as NON_SCIENCE instead of hole if left, right, top, bottom pixel # is NON_SCIENCE @@ -2165,18 +2496,25 @@ def set_final_dq_flags(self): xcheck[3] = xrem ycheck[3] = yrem + 1 - while ((ij < 4) and (found == 0)): - if (xcheck[ij] > 0 and xcheck[ij] < self.naxis1 and - ycheck[ij] > 0 and ycheck[ij] < self.naxis2): + while (ij < 4) and (found == 0): + if ( + xcheck[ij] > 0 + and xcheck[ij] < self.naxis1 + and ycheck[ij] > 0 + and ycheck[ij] < self.naxis2 + ): index_check = iwave * nxy + ycheck[ij] * self.naxis1 + xcheck[ij] # If the nearby spaxel_dq contains overlap_no_coverage # then unmark dq flag as hole. A hole has to have nearby # pixels all in FOV. - check = (np.bitwise_and(self.spaxel_dq[index_check], - self.overlap_no_coverage) == self.overlap_no_coverage) + check = ( + np.bitwise_and(self.spaxel_dq[index_check], self.overlap_no_coverage) + == self.overlap_no_coverage + ) if check: - spaxel_dq_temp[index[0][i]] = np.bitwise_or(self.overlap_no_coverage, - dqflags.pixel['DO_NOT_USE']) + spaxel_dq_temp[index[0][i]] = np.bitwise_or( + self.overlap_no_coverage, dqflags.pixel["DO_NOT_USE"] + ) found = 1 ij = ij + 1 @@ -2185,17 +2523,25 @@ def set_final_dq_flags(self): ave_holes = len(location_holes[0]) / self.naxis3 if ave_holes < 1: - log.info('Average # of holes/wavelength plane is < 1') + log.info("Average # of holes/wavelength plane is < 1") else: - log.info('Average # of holes/wavelength plane: %i', ave_holes) - log.info('Total # of holes for IFU cube is : %i', len(location_holes[0])) + log.info("Average # of holes/wavelength plane: %i", ave_holes) + log.info("Total # of holes for IFU cube is : %i", len(location_holes[0])) - # ******************************************************************************** + # ________________________________________________________________________________ def setup_final_ifucube_model(self, model_ref): - """ Set up the final meta WCS info of IFUCube along with other fits keywords + """ + Set up the final meta WCS info of IFUCube along with other fits keywords. - return IFUCube model + Parameters + ---------- + model_ref : IFUImageModel + The first IFUImage model to use to fill in basic header values + Returns + ------- + result : IFUCubeModel + IFU cube datamodel with data arrays filled in. """ status = 0 # loop over the wavelength planes to confirm each plane has some data @@ -2203,22 +2549,18 @@ def setup_final_ifucube_model(self, model_ref): # from the IFUcube # Rearrange values from 1d vectors into 3d cubes - flux = self.spaxel_flux.reshape((self.naxis3, - self.naxis2, self.naxis1)) - wmap = self.spaxel_iflux.reshape((self.naxis3, - self.naxis2, self.naxis1)) - - var = self.spaxel_var.reshape((self.naxis3, - self.naxis2, self.naxis1)) - dq = self.spaxel_dq.reshape((self.naxis3, - self.naxis2, self.naxis1)) - - # For MIRI MRS, apply a quality cut to help fix spectral tearing at the ends of each band. - # This is largely taken care of by the WCS regions file, but there will still be 1-2 possibly - # problematic planes at the end of each band in multi-band cubes. - # Do this by looking for how many good spaxels there are at each wavelength and finding outliers - # from the trend. - if self.instrument == 'MIRI': + flux = self.spaxel_flux.reshape((self.naxis3, self.naxis2, self.naxis1)) + wmap = self.spaxel_iflux.reshape((self.naxis3, self.naxis2, self.naxis1)) + + var = self.spaxel_var.reshape((self.naxis3, self.naxis2, self.naxis1)) + dq = self.spaxel_dq.reshape((self.naxis3, self.naxis2, self.naxis1)) + + # For MIRI MRS, apply a quality cut to help fix spectral tearing at the ends of + # each band. This is largely taken care of by the WCS regions file, but there + # will still be 1-2 possibly problematic planes at the end of each band in + # multi-band cubes. Do this by looking for how many good spaxels there are at + # each wavelength and finding outliers from the trend. + if self.instrument == "MIRI": nz = flux.shape[0] # Create a vector of the number of good spaxels at each wavelength ngood = np.zeros(nz) @@ -2234,35 +2576,31 @@ def setup_final_ifucube_model(self, model_ref): # and zero out those arrays. lowcov = (np.where((ngood > 0) & (ngood < 0.75 * pctile)))[0] nlowcov = len(lowcov) - log.info('Number of spectral tear planes adjusted: %i', nlowcov) + log.info("Number of spectral tear planes adjusted: %i", nlowcov) for zz in range(0, nlowcov): flux[lowcov[zz], :, :] = 0 wmap[lowcov[zz], :, :] = 0 var[lowcov[zz], :, :] = 0 - dq[lowcov[zz], :, :] = dqflags.pixel['DO_NOT_USE'] + dqflags.pixel['NON_SCIENCE'] + dq[lowcov[zz], :, :] = ( + dqflags.pixel["DO_NOT_USE"] + dqflags.pixel["NON_SCIENCE"] + ) # Set np.nan values wherever the DO_NOT_USE flag is set - dnu = np.where((dq & dqflags.pixel['DO_NOT_USE']) != 0) + dnu = np.where((dq & dqflags.pixel["DO_NOT_USE"]) != 0) flux[dnu] = np.nan var[dnu] = np.nan var = np.sqrt(var) if self.linear_wavelength: - ifucube_model = datamodels.IFUCubeModel(data=flux, dq=dq, - err=var, - weightmap=wmap) + ifucube_model = datamodels.IFUCubeModel(data=flux, dq=dq, err=var, weightmap=wmap) else: wave = np.asarray(self.wavelength_table, dtype=np.float32) num = len(wave) - alldata = np.array( - [(wave[None].T, )], - dtype=[('wavelength', ' #include @@ -13,11 +21,23 @@ extern double sh_find_overlap(const double xcenter, const double ycenter, double xPixelCorner[], double yPixelCorner[]); -//________________________________________________________________________________ -// allocate the memory for the spaxel DQ array - int mem_alloc_dq(int nelem, int **idqv) { - + /* + Allocate the memory of the spaxel DQ array. + + Parameters + --------- + nelem : int + Number of elements to allocate + idqv : pointer + Pointer to memory allocated, set by this routine. + + Returns + ------- + status : int + 1 = error + 0 = success + */ const char *msg = "Couldn't allocate memory for output arrays."; if (!(*idqv = (int*)calloc(nelem, sizeof(int)))) { @@ -29,8 +49,6 @@ int mem_alloc_dq(int nelem, int **idqv) { } -//________________________________________________________________________________ -// Routine for MIRI DQ plane assignment int corner_wave_plane_miri(int w, int start_region, int end_region, double roiw_ave, @@ -40,26 +58,46 @@ int corner_wave_plane_miri(int w, int start_region, int end_region, long ncube, long npt, double *corner1, double *corner2, double *corner3, double *corner4) { /* - For wavelength plane determine the corners (in xi,eta) of the FOV for MIRI - Use the 2 extreme slices set by start_region and end_region to define the FOV of the wavelength plane FOV - Using the min and max coordinates for the extent of the slice on the sky for these two slice - set the + + Determine the MIRI DQ plane of IFU cubes. + + For each wavelength plane determine the corners (in xi,eta) of the FOV. + Use the 2 extreme slices set by start_region and end_region to define the FOV of the wavelength plane. + Using the min and max coordinates for the extent of the slice on the sky for these two slice, set the corners of the FOV. - w : wavelength plane - start_region : starting slice # for channel (slice # in nirspec) - end_region : ending slice # for channel (slice # in nirspec) - roiw_ave : average roiw for all wavelengths - zc : array of wavelengths - coord1 : point cloud xi values - coord2 : point cloud eta values - wave : point cloud wavelength values - sliceno : point cloud slice no - ncube : number of cube values - npt: number of point cloud elements - corner1 : xi, eta of corner 1 - corner2 : xi, eta of corner 2 - corner3 : xi, eta of corner 3 - corner4 : xi, eta of corner 4 + Parameters + ---------- + w : int + Wavelength plane + start_region : int + Starting slice # for channel (slice # in nirspec) + end_region : int + Ending slice # for channel (slice # in nirspec) + roiw_ave : double + Average roiw for all wavelengths + zc : ndarray + Array of wavelength + coord1 : ndarray + Detector xi values for the center of the pixel + coord2 : ndarray + Detector eta values for the center of the pixel + wave : ndarray + Wavelength values for the center of pixel + sliceno : ndarray + Slice no of the pixels + ncube: long + Number of cube elements + npt: int + Number of detector pxiels + corner1 : ndarray + xi, eta of corner 1 + corner2 : ndarray + xi, eta of corner 2 + corner3 : ndarray + xi, eta of corner 3 + corner4 : ndarray + xi, eta of corner 4 */ int slice, c1_use; @@ -210,8 +248,6 @@ int corner_wave_plane_miri(int w, int start_region, int end_region, } } -//________________________________________________________________________________ -// MIRI DQ routine. Find the overlap of the FOV for the wavelength slice in IFU cube int overlap_fov_with_spaxels(int overlap_partial, int overlap_full, double cdelt1, double cdelt2, @@ -220,40 +256,50 @@ int overlap_fov_with_spaxels(int overlap_partial, int overlap_full, double xi_corner[], double eta_corner[], int wave_slice_dq[]) { - /* find the amount of overlap of FOV each spaxel + /* MIRI routine to find the overlap of the FOV for a wavelength slice in the IFU cube. Given the corners of the FOV find the spaxels that - overlap with this FOV. Set the intermediate spaxel to + overlap with this FOV. Set the intermediate spaxel to a value based on the overlap between the FOV for each exposure and the spaxel area. The values assigned are: a. overlap_partial b overlap_full - bit_wise combination of these values is allowed to account for + Bitwise combination of these values is allowed to account for dithered FOVs. Parameters ---------- overlap_partial : int + DQ flag indicating there is a partial overlap of the cube spaxel for this + wavelength slice overlap_full : int + DQ flag indicated there is full overlap of the cube spaxel for this + wavelength slice. cdelt1 : double - IFU cube naxis 1 spatial scale + IFU cube naxis 1 spatial scale cdelt2 : double - IFU cube naxis 1 spatial scale + IFU cube naxis 1 spatial scale naxis1 : int - IFU cube naxis 1 size + IFU cube naxis 1 size naxis2 : int - IFU cube naxis 1 size + IFU cube naxis 1 size xcenters : double array - IFU naxis 1 array of xi values + IFU naxis 1 array of xi values ycenters : double array - IFU naxis 2 array of eta values - xi_corner: xi coordinates of the 4 corners of the FOV on the wavelength plane - eta_corner: eta coordinates of the 4 corners of the FOV on the wavelength plane + IFU naxis 2 array of eta values + xi_corner: ndarray + xi coordinates of the 4 corners of the FOV on the wavelength plane + eta_corner: ndarray + eta coordinates of the 4 corners of the FOV on the wavelength plane Sets - ------- + ---- wave_slice_dq: array containing intermediate dq flag + Returns + ------- + 0 = success + */ // loop over spaxels in the wavelength plane and set slice_dq @@ -307,8 +353,6 @@ int overlap_fov_with_spaxels(int overlap_partial, int overlap_full, return 0 ; } -//________________________________________________________________________________ -// Routine to setting NIRSpec dq plane for each wavelength plane long match_wave_plane_nirspec(double wave_plane, double roiw_ave, @@ -321,21 +365,42 @@ long match_wave_plane_nirspec(double wave_plane, double *c1_max, double *c2_max, int *match_slice){ - /* - NIRSpec dq plane is set by mapping each slice to IFU wavelength plane + /* Determine the NIRSpec DQ plane for each wavelength plane of an IFU cube. + + NIRSpec dq plane is set by mapping each slice to the IFU wavelength plane. This routine maps each slice to sky and finds the min and max coordinates on the sky of the slice. - wave_plane : wavelength of current plane - slicevalue : slice # 1 to 30 - roiw_ave : average roiw for all wavelengths - coord1, coord2: tangent project coordinate of pt cloud - wave : point cloud wavelength values - sliceno: slice value of point cloud. - npt: number of point cloud elements + Parameters: + ---------- + wave_plane : double + Wavelength of current plane + roiw_ave : double + Average roiw for all wavelengths + + coord1 : ndarray + Detector xi values for the center of the pixel + coord2 : ndarray + Detector eta values for the center of the pixel + wave : ndarray + Wavelength values for the center of pixel + sliceno : ndarray + Slice no of the pixels + npt: int + Number of detector pxiels + + Sets: + c1_min, c2_min, c1_max, c2_max : array of size 30 + C1_min and c1_max are the min and max xi coordinates of a slice. + C2_min and c2_max are the min and max eta coordinates of a slice. + + match_slice : + Integer flag that marks slices that the min and max coordinates have been determined. - return: - c1_min, c2_min, c1_max, c2_max, match_slice + Returns + ------- + ii : long + Number of pixels with DQ set. */ long ipt =0; @@ -404,8 +469,7 @@ long match_wave_plane_nirspec(double wave_plane, } -//________________________________________________________________________________ -// NIRSpec Find the overlap of the slices for the wavelength slice with sky + int overlap_slice_with_spaxels(int overlap_partial, double cdelt1, double cdelt2, int naxis1, int naxis2, @@ -415,29 +479,35 @@ int overlap_slice_with_spaxels(int overlap_partial, int wave_slice_dq[]) { /* - Set the initial dq plane indicating if the input data falls on a spaxel + Set the NIRSpec initial dq plane indicating if the input data falls on a spaxel. This algorithm assumes the input data falls on a line in the IFU cube, which is - the case for NIRSpec slices. The NIRSpec slice's endpoints are used to determine + the case for NIRSpec slices. The NIRSpec slice endpoints are used to determine which IFU spaxels the slice falls on to set an initial dq flag. - + The Bresenham's Line Algorithm is used to find the spaxels that intersects with a line. + Parameters ---------- overlap_partial: intermediate dq flag - + cdelt1 : spatial sampling size for naxis1 + cdelt2 : spatial sampling size for naxis2 + naxis1 : size of axis 1 + naxis2 : size of axis 2 + xstart : x start of slice + ystart : y start of slixe + xi_min : xi min of cube + eta_min : eta min of cube + xi_max : xi max of cube + eta_max : eta max of cube wavelength: the wavelength bin of the IFU cube working with Sets ---- wave_slice_dq : array containing intermediate dq flag - Bresenham's Line Algorithm to find points a line intersects with grid. - - Given the endpoints of a line find the spaxels this line intersects. - Returns ------- - wave_slice_dq filled in + 0 */ @@ -508,12 +578,27 @@ int overlap_slice_with_spaxels(int overlap_partial, } //________________________________________________________________________________ -// Set the spaxel dq = 0. This is used when not determining the FOV on the sky for -// setting the DQ plane. This is case for internalCal type cubes - int set_dqplane_to_zero(int ncube, int **spaxel_dq){ + /* + Set the spaxel dq = 0. + + This is used when not determining the FOV on the sky for + setting the DQ plane. This is case for internal cal type cubes. + + Parameters + ---------- + ncube : nint + Number of elements in cube + spaxel_dq : pointer + Pointer to memory + + Returns + ------- + 0 success + 1 failure + */ int *idqv; // int vector for spaxel if (mem_alloc_dq(ncube, &idqv)) return 1; *spaxel_dq = idqv; @@ -521,7 +606,6 @@ int set_dqplane_to_zero(int ncube, int **spaxel_dq){ } //________________________________________________________________________________ -// Main MIRI routine to set DQ plane int dq_miri(int start_region, int end_region, int overlap_partial, int overlap_full, int nx, int ny, int nz, @@ -532,6 +616,36 @@ int dq_miri(int start_region, int end_region, int overlap_partial, int overlap_f long ncube, long npt, int **spaxel_dq) { + /* + + Main MIRI routine to set DQ plane. + + Parameters + ---------- + start_region : x pixel for start of channel + end_region : x pixel for end of channel + overlap_partial : intermediate dq flag + overlap_full : dq flag + nx, ny, nz : size of cube + cdelt1 : spatial sampling size for naxis1 + cdelt2 : spatial sampling size for naxis2 + roiw_ave : average region of interest for wavelength + xc , xy , zc : vectors of cube centers for naxis1, 2, 3 + coord1, coord2 : xi, eta coordinates of the cube + wave : wavelength vector of the cube + sliceno : sliceno of pixel + ncube : number of elements in cube + npt : number detector pixels mapping + + Sets + ---- + spaxel_dq : array containing dq flag + + Returns + ------- + 0 : success + 1 : failure + */ int status, status_wave, w, nxy, i, istart, iend, in, ii; double xi_corner[4], eta_corner[4]; @@ -599,9 +713,7 @@ int dq_miri(int start_region, int end_region, int overlap_partial, int overlap_f return 0; } -//________________________________________________________________________________ -//Main NIRSpec routine to set up dq plane - + int dq_nirspec(int overlap_partial, int nx, int ny, int nz, double cdelt1, double cdelt2, double roiw_ave, @@ -612,7 +724,7 @@ int dq_nirspec(int overlap_partial, int **spaxel_dq) { /* - Set an initial DQ flag for the NIRSPEC IFU cube based on FOV of input data. + Set a DQ flag for the NIRSPEC IFU cube based on FOV of input data. Map the FOV of each NIRSpec slice to the DQ plane and set an initial DQ flagging. For NIRSpec the 30 different slices map to different FOV on the @@ -624,13 +736,29 @@ int dq_nirspec(int overlap_partial, Parameter --------- overlap_partial - intermediate DQ flag + nx, ny, nz : size of cube + cdelt1 : spatial sampling size for naxis1 + cdelt2 : spatial sampling size for naxis2 + roiw_ave : average region of interest for wavelength + xc , xy , zc : vectors of cube centers for naxis1, 2, 3 coord1: xi coordinates of input data (~x coordinate in IFU space) coord2: eta coordinates of input data (~y coordinate in IFU space) wave: wavelength of input data roiw_ave: average spectral roi used to determine which wavelength bins - the input values would be mapped to + the input values would be mapped to slice_no: integer slice value of input data (used in MIRI case to find - the points of the edge slices.) + the points of the edge slices.) + + Sets + ---- + spaxel_dq : ndarray + Spaxel dq plane of IFU cube. + + Returns + ------- + status : int + 1 = error + 0 = success */ int w, islice, status, status_wave, nxy, j; diff --git a/jwst/cube_build/src/cube_match_internal.c b/jwst/cube_build/src/cube_match_internal.c index 036143a811..c68e326455 100644 --- a/jwst/cube_build/src/cube_match_internal.c +++ b/jwst/cube_build/src/cube_match_internal.c @@ -1,5 +1,5 @@ /* -This module forms the IFU cube in the "detector plane". This method of creating cubes is +This module creates the IFU cube in the "slicer plane". This method of creating cubes is for an engineering mode. The detector pixels are mapped to the along slice coordinate and wavelength. The slice number represents the across slice coordinate. The pixels from each slice are assumed to not overlap on the output plane with pixels from other slices. This @@ -25,50 +25,54 @@ example output Parameters ---------- instrument: int - 0 = MIRI, 1 = NIRSPEC. Used for set the dq plane + 0 = MIRI, 1 = NIRSPEC. Used for set the dq plane naxis1 : int - axis 1 of IFU cube -nasix2 : int - axis 2 of IFU cube + Size of axis 1 of the IFU cube. +naxis2 : int + Size of axis 2 of the IFU cube. crval_along : double - along slice reference value in IFU cube + Along slice reference value in IFU cube cdelt_along : double - along slice sampling in IFU cube + Along slice sampling size in IFU cube crval3 : double wavelength reference value in IFU cube cdetl3 : double - wavelength sampling in IFU cube + wavelength sampling size in IFU cube a1, a2, a3, a4: double array - array of corners of pixels holding along slice coordinates + Array of corners of pixels holding along slice coordinates. a1 holds corner 1 a2 holds corner 2 a3 holds corner 3 a4 holds corner 4 lam1, lam2, lam3, lam4: double array - array of corners of pixels holding wavelength coordinates + Array of corners of pixels holding wavelength coordinates. lam1 holds corner 1 lam2 holds corner 2 lam3 holds corner 3 lam4 holds corner 4 acoord : double array - Array holds along slice coordinates in IFU cube + Array holding along slice coordinates in the IFU cube. zcoord : double array - Array holds the wavelength coordinates in the IFU cube + Array holding the wavelength coordinates in the IFU cube. ss : double array - Array holds the slice number + Array holding the slice number. pixel_flux : double array - Array that holds detector pixel flux + Array that holding detector pixel flux. pixel_error : double array - Array that holds detector pixel flux error + Array that holding detector pixel flux error. Returns ------- - -spaxel_flux : array -spaxel_weight : array -spaxel_iflux : array -spaxel_var : array -spaxel_dq : array +spaxel_flux : double array + Combined flux from overlapping detector pixels. +spaxel_weight : double array + Combined weighting used from overlapping detector pixels. +spaxel_iflux : double array + Combined counter from overlapping detector pixels. +spaxel_var : double array + Combined variance from overlapping detector pixels. +spaxel_dq : double array + Spaxel dq flux indicating if spaxel is overlapped with any detector pixels */ @@ -96,10 +100,6 @@ extern double find_area_quad(double MinX, double MinY, double Xcorner[], double extern int alloc_flux_arrays(int nelem, double **fluxv, double **weightv, double **varv, double **ifluxv); -//_______________________________________________________________________ -// based on a 2-D drizzling type of algorithm find the contribution of -// each detector flux to the IFU cube -//_______________________________________________________________________ int match_detector_cube(int instrument, int naxis1, int naxis2, int nz, int npt, int ncube, int na, double crval_along, double cdelt_along, double crval3, double cdelt3, @@ -109,6 +109,14 @@ int match_detector_cube(int instrument, int naxis1, int naxis2, int nz, int npt, double *pixel_flux, double *pixel_err, double **spaxel_flux, double **spaxel_weight, double **spaxel_var,double **spaxel_iflux){ + /* Find the overlap between the detector pixels and spaxels of the IFU cube when the IFU cube is of internal cal type. + + The overlap is based on a 2-D drizzling type of algorithm which finds the contribution of + each detector flux to the IFU cube spaxel + + A description of each variable is given at the top of this c routine. + */ + double *fluxv=NULL, *weightv=NULL, *varv=NULL, *ifluxv=NULL; // vectors for spaxel double along_corner[4], wave_corner[4], along_min, wave_min, along_max, wave_max, Area, MinW, MaxW, zcenter, acenter, area_overlap, AreaRatio, err; @@ -204,7 +212,8 @@ int match_detector_cube(int instrument, int naxis1, int naxis2, int nz, int npt, // C extension SETUP - +// This ensures that all the numpy arrays passed to the C routines +// follow C array rules. PyArrayObject * ensure_array(PyObject *obj, int *is_copy) { if (PyArray_CheckExact(obj) && PyArray_IS_C_CONTIGUOUS((PyArrayObject *) obj) && @@ -221,6 +230,8 @@ PyArrayObject * ensure_array(PyObject *obj, int *is_copy) { } +// Wrapper code that is called from python code and sets up interface with C code. + static PyObject *cube_wrapper_internal(PyObject *module, PyObject *args) { PyObject *result = NULL, *a1o, *a2o, *a3o, *a4o, *lam1o, *lam2o, *lam3o, *lam4o, diff --git a/jwst/cube_build/src/cube_match_sky_driz.c b/jwst/cube_build/src/cube_match_sky_driz.c index fe6e60c35f..1bd4ef2985 100644 --- a/jwst/cube_build/src/cube_match_sky_driz.c +++ b/jwst/cube_build/src/cube_match_sky_driz.c @@ -1,23 +1,39 @@ /* -The detector pixels are represented by a 'point cloud' on the sky. The IFU cube is -represented by a 3-D regular grid. This module finds the point cloud members contained -in a region centered on the center of the cube spaxel. The size of the spaxel is spatial -coordinates is cdetl1 and cdelt2, while the wavelength size is cdelt3. +An IFU spectral cube is built using the overlap between the detector pixels mapped +to the 3-D spectral grid. This method is similar to 2D drizzle type mapping techniques. +This is the top level c routine that interfaces with the python wrapper. Main function for Python: cube_wrapper_driz -Python signature: result = cube_wrapper_driz(instrument, flag_dq_plane, start_region, end_region, - overlap_partial, overlap_full, - xcoord, ycoord, zcoord, - coord1, coord2, wave, flux, err, slice_no, - rois_pixel, roiw_pixel, scalerad_pixel - weight_pixel, softrad_pixel,cdelt3_normal, - roiw_ave, cdelt1, cdelt2, x_det, y_det, debug_cube_index) -provide more details +Python signature: result = cube_wrapper_driz( + instrument, + flag_dq_plane, + start_region, end_region, + overlap_partial, + overlap_full, + xcoord, ycoord, zcoord, + coord1, coord2, + wave, + flux, + err, + slice_no, + xi1, eta1, + xi2, eta2, + xi3, eta3, + xi4, eta4, + dwave, + cdelt3_normal, + cdelt1, + cdelt2, + cdelt3_mean, + linear, + x_det, + y_det, + debug_cube_index + ) The output of this function is a tuple of 5 arrays:(spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux, spaxel_dq) -example output Parameters ---------- @@ -26,7 +42,6 @@ instrument : int flag_dq_plane : int 0 do set the DQ plane based on FOV, but set all values =0 1 set the DQ plane based on the FOV - start_region : int starting slice number for detector region used in dq flagging end_region: int @@ -35,33 +50,48 @@ overlap_partial : int a dq flag indicating that only a portion of the spaxel is overlapped by a mapped detector pixel overlap_full : int a dq flag indicating that the entire spaxel is overlapped by the mapped detector pixel -xcoord : double array +xcoord : ndarray of float size of naxis1. This array holds the center x axis values of the ifu cube -ycoord : double array +ycoord : ndarray of float size of naxis2. This array holds the center y axis values of the ifu cube -zcoord : double array - size of naxis3. This array holds the center x axis values of the ifu cube -flux : double array +zcoord : ndarray of float + size of naxis3. This array holds the center z axis /wavelength values of the ifu cube +coord1 : ndarray of float + The tangent projected xi values of the center of detector pixel. +coord2 : ndarray of float + The tangent projected wave values of the center of detector pixel.. +wave : ndarray of float + The wavelength corresponding to the center of the detector pixel. +flux : ndarray of float size: point cloud elements. Flux of each point cloud member -err : double array +err : ndarray of float size: point cloud elements. err of each point cloud member -slice_no: int +slice_no : int slice number of point cloud member to be in dq flagging -coord1 : double array - size: point cloud elements. Naxis 1 coordinate of point cloud member (xi) -coord2 : double array - size: point cloud elements. Naxis 2 coordinate of point cloud member (eta) -wave : double array - size: point cloud elements. Wavelength of each point cloud member -cdelt3: double array - size: point cloud elements. Spectral scale to use at wavelength of point cloud member +xi1, eta1 : ndarray of floats + xi, eta coordinates (tangent projection) of the detector pixel corner 1. +xi2, eta2 : ndarray of floats + xi, eta coordinates (tangent projection) of the detector pixel corner 2. +xi3, eta3 : ndarray of floats + xi, eta coordinates (tangent projection) of the detector pixel corner 3. +xi4, eta4 : ndarray of floats + xi, eta coordinates (tangent projection) of the detector pixel corner 4. +dwave : ndarray of float + Delta wavelength ranges for pixel +cdelt3_normal : ndarray of float + For linear wavelength range = 0 + For non-linear wavelength range = delta wavelength of neighboring wavelength planes cdelt1 : double Naxis 1 scale for cube cdelt2 : double Naxis 2 scale for cube -x_det : double array +cdelt3_mean: ndarray of float + Average of cdelt3_normal +linear : bool + Type of wavelength grid for IFU cube: linear or non-linear +x_det : ndarray of float size: point cloud elements. X detector value of each point cloud member -y_det : double array +y_det : ndarray of float size: point cloud elements. Y detector value of each point cloud member debug_cube_index : int if > 0, value of cube index to print information on @@ -122,9 +152,7 @@ extern double sh_find_overlap(double xcenter, double ycenter, // extern double find_area_quad(double MinX, double MinY, double Xcorner[], double Ycorner[]); - -// return values: spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux - +// C routine that does that does the drizzling int match_driz(double *xc, double *yc, double *zc, double *wave, double *flux, double *err, @@ -138,6 +166,31 @@ int match_driz(double *xc, double *yc, double *zc, double **spaxel_flux, double **spaxel_weight, double **spaxel_var, double **spaxel_iflux) { + // xc : IFU grid point values along x-axis + // yc : IFU grid point values along y-axis + // zc : IFU grid point values along the z-axis + // wave : wavelength of pixels + // flux : flux values of pixels + // err : err values of the pixels + // xi1, eta1 : xi, eta coordinates of a corner 1 of a pixel + // xi2, eta2 : xi, eta coordinates of a corner 2 of a pixel + // xi3, eta3 : xi, eta coordinates of a corner 3 of a pixel + // xi4, eta4 : xi, eta coordinates of a corner 4 of a pixel + // dwave : delta wavelength of pixel + // cdelt3: wavelength sampling of IFU cube + // x_det, y_det :x, y detector values + // cdelt1 : spatial size along x-axis + // cdelt2 : spatial size along y-axis + // nx, ny, nwave : size of IFU cube in x,y,z dimensions + // ncube : total numober cube elements + // npt : number of detector pixels + // linear: If T then, the IFU cube has a linear wavelength grid. + // If F then, the IFU cube has a non-linear wavelength grid. + // debug_cube_index: cube index of voxel to report information on + // spaxel_flux : return value of weighted combined flux + // spaxel_weight : return value of combined weights + // spaxel_var : return value of weighted combined variance + // spaxel_iflux :return value of weighted combined iweighting map double *fluxv=NULL, *weightv=NULL, *varv=NULL, *ifluxv=NULL; // vector for spaxel @@ -310,7 +363,9 @@ int match_driz(double *xc, double *yc, double *zc, } - +// This is a C structure that represents a NumPy array in the C API. +// This ensures that all the numpy arrays passed to the C routines +// follow C array rules. PyArrayObject * ensure_array(PyObject *obj, int *is_copy) { if (PyArray_CheckExact(obj) && PyArray_IS_C_CONTIGUOUS((PyArrayObject *) obj) && @@ -327,6 +382,9 @@ PyArrayObject * ensure_array(PyObject *obj, int *is_copy) { } +// Wrapper code that is called from python code and sets up interface with C code. +// cube_wrapper_driz. + static PyObject *cube_wrapper_driz(PyObject *module, PyObject *args) { PyObject *result = NULL, *xco, *yco, *zco, *fluxo, *erro, *coord1o, *coord2o, *waveo, *slicenoo; // codespell:ignore erro PyObject *cdelt3o; @@ -445,9 +503,8 @@ static PyObject *cube_wrapper_driz(PyObject *module, PyObject *args) { } - //______________________________________________________________________ - // if flag_dq_plane = 1, Set up the dq plane - //______________________________________________________________________ + // if flag_dq_plane = 1, Set up the dq plane of the IFU cube. + int status1 = 0; if(flag_dq_plane){ diff --git a/jwst/cube_build/src/cube_match_sky_pointcloud.c b/jwst/cube_build/src/cube_match_sky_pointcloud.c index c62b3605e0..adc2c0abbf 100644 --- a/jwst/cube_build/src/cube_match_sky_pointcloud.c +++ b/jwst/cube_build/src/cube_match_sky_pointcloud.c @@ -1,9 +1,12 @@ /* -The detector pixels are represented by a 'point cloud' on the sky. The IFU cube is -represented by a 3-D regular grid. This module finds the point cloud members contained -in a region centered on the center of the cube spaxel. The size of the spaxel is spatial +If the weighting method for IFU cube building is set to 'msm' or 'emsm' this c method +is used. +The detector pixels are represented by a 'point cloud' on the sky. Each point cloud member +contains the coordinates for the pixel center, the flux, error and variance of that pixel. +The IFU cube is represented by a 3-D regular grid. This module finds the point cloud members contained +in a region centered on the center of the cube spaxel. The size of the spaxel in spatial coordinates is cdetl1 and cdelt2, while the wavelength size is zcdelt3. -This module uses the modified shephard weighting method (emsm if weight_type =0 or msm if weight_type =1) +This module uses the modified Shepard weighting method (emsm if weight_type =0 or msm if weight_type =1) to determine how to weight each point cloud member in the spaxel. Main function for Python: cube_wrapper @@ -31,39 +34,39 @@ weight_type : int 0: use emsm weighting 1: use msm weighting start_region : int - starting slice number for detector region used in dq flagging + Starting slice number for detector region used in dq flagging end_region: int - ending slice number for detector region used in dq flagging + Ending slice number for detector region used in dq flagging overlap_partial : int - a dq flag indicating that only a portion of the spaxel is overlapped by a mapped detector pixel + A dq flag indicating that only a portion of the spaxel is overlapped by a mapped detector pixel overlap_full : int - a dq flag indicating that the entire spaxel is overlapped by the mapped detector pixel + A dq flag indicating that the entire spaxel is overlapped by the mapped detector pixel xcoord : double array - size of naxis1. This array holds the center x axis values of the ifu cube + This array holds the center x axis values of the ifu cube. It has a size of naxis1. ycoord : double array - size of naxis2. This array holds the center y axis values of the ifu cube + This array holds the center y axis values of the ifu cube. It has a size of naxis2. zcoord : double array - size of naxis3. This array holds the center x axis values of the ifu cube + This array holds the center x axis values of the ifu cube. It has a size of naxis3. flux : double array - size: point cloud elements. Flux of each point cloud member + Flux of each point cloud member. err : double array - size: point cloud elements. err of each point cloud member + Error of each point cloud member slice_no: int - slice number of point cloud member to be in dq flagging + Slice number of point cloud member to be in dq flagging coord1 : double array - size: point cloud elements. Naxis 1 coordinate of point cloud member (xi) + Naxis 1 coordinate of point cloud member (xi). coord2 : double array - size: point cloud elements. Naxis 2 coordinate of point cloud member (eta) + Naxis 2 coordinate of point cloud member (eta). wave : double array - size: point cloud elements. Wavelength of each point cloud member + Wavelength of each point cloud member. rois_pixel : double array - size: point cloud elements. Roi in spatial dimension to use for point cloud member + Roi in spatial dimension to use for point cloud member. roiw_pixel : double array - size: point cloud elements. Roi in wavelength dimension to use point cloud member + Roi in wavelength dimension to use point cloud member. scalerad_pixel : double array - size: point cloud elements. MSM weight parameter to use for point cloud member + MSM weight parameter to use for point cloud member zcdelt3: double array - size: point cloud elements. Spectral scale to use at wavelength of point cloud member + Spectral scale to use at wavelength of point cloud member roiw_ave : double Average roiw for all the wavelength planes. Used in dq flagging cdelt1 : double @@ -72,17 +75,17 @@ cdelt2 : double Naxis 2 scale for cube -Returns -------- -spaxel_flux : numpy.ndarray +Sets +---- +spaxel_flux : ndarray IFU spaxel cflux -spaxel_weight : numpy.ndarray +spaxel_weight : ndarray IFU spaxel weight -spaxel_iflux : numpy.ndarray +spaxel_iflux : ndarray IFU spaxel weight map (number of overlaps) -spaxel_var : numpy.ndarray +spaxel_var : ndarray IFU spaxel error -spaxel_dq : numpy.ndarray +spaxel_dq : ndarray IFU spaxel dq */ @@ -123,9 +126,11 @@ extern int dq_nirspec(int overlap_partial, extern int set_dqplane_to_zero(int ncube, int **spaxel_dq); // Match point cloud to sky and determine the weighting to assign to each point cloud member -// to matched spaxel based on ROI - weighting type - emsm +// to matched spaxel based on ROI - weighting type is emsm. -// return values: spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux +// sets values: spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux + +// returns : 0 success, 1 failure int match_point_emsm(double *xc, double *yc, double *zc, double *coord1, double *coord2, double *wave, @@ -283,8 +288,9 @@ int match_point_emsm(double *xc, double *yc, double *zc, // Match point cloud to sky and determine the weighting to assign to each point cloud member -// to matched spaxel based on ROI - weighting type - msm -//return values: spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux +// to matched spaxel based on ROI - weighting type = msm. +// set values: spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux. +// returns: 0 = success, 1 = failure. int match_point_msm(double *xc, double *yc, double *zc, double *coord1, double *coord2, double *wave, @@ -462,7 +468,7 @@ PyArrayObject * ensure_array(PyObject *obj, int *is_copy) { } } - +// Wrapper code that is called from python code and sets up interface with C code. static PyObject *cube_wrapper(PyObject *module, PyObject *args) { PyObject *result = NULL, *xco, *yco, *zco, *fluxo, *erro, *coord1o, *coord2o, *waveo, *slicenoo; // codespell:ignore erro PyObject *rois_pixelo, *roiw_pixelo, *scalerad_pixelo, *zcdelt3o, *softrad_pixelo, *weight_pixelo; diff --git a/jwst/cube_build/src/cube_utils.c b/jwst/cube_build/src/cube_utils.c index 92c9a80cfe..21a63b2611 100644 --- a/jwst/cube_build/src/cube_utils.c +++ b/jwst/cube_build/src/cube_utils.c @@ -1,4 +1,10 @@ -// This library contains c functions to support building IFU cubes +/* This library contains c functions to support building IFU cubes. +Many of these routines are used in the Sutherland-Hodgman algorithm. This +algorithm finds the clipped polygon that falls inside the cube spaxel. +A clipped polygon is the overlapping polygon of the detector pixel and +a spaxel. We are only dealing with the spatial dimensions in this routine. +*/ + #include #include @@ -12,12 +18,23 @@ #define CP_TOP 3 -//_______________________________________________________________________ -// Allocate the memory to for the spaxel arrays in the cube -//_______________________________________________________________________ - int alloc_flux_arrays(int nelem, double **fluxv, double **weightv, double **varv, double **ifluxv) { + /* + Allocate memory for the spaxel output vectors to be of size nelem. + + nelem : int + Number of elements to allocate memory + fluxv : double ndarray + Flux vector + weightv : double ndarray + Weight vector + varv : double ndarray + Variance vector + iflux : double ndarray + Counter index vector + */ + const char *msg = "Couldn't allocate memory for output arrays."; // flux: @@ -57,16 +74,49 @@ int alloc_flux_arrays(int nelem, double **fluxv, double **weightv, double **varv return 1; } -// support function for sh_find_overlap void addpoint (double x, double y, double xnew[], double ynew[], int *nVertices2){ + + /* + A support function for the sh_find_overlap (Sutherland-Hodgman algorithm). + Adds a new x and y point to the vertices describing the clipped polygon. + + x : double + X coordinate to add to clipped polygon that falls inside the cube spaxel + y : double + Y coordinate to add to clipped polygon that falls inside the cube spaxel + xnew : double array + Array containing the x vertices of clipped polygon + ynew : double array + Array containing the y vertices of clipped polygon + */ xnew[*nVertices2] = x; ynew[*nVertices2] = y; *nVertices2 = *nVertices2 + 1; } -// support function for sh_find_overlap + int insideWindow(int edge, double x, double y, double left,double right, double top, double bottom){ + /* + Support function for sh_find_overlap. Determine where a point is in relationship to + the clipped polygon. + + edge : int + Integer value which defines which edge of the polygon we are working with. + 0 = left, 1 = right, 2 = top, 3 = bottom. + x : double + X coordinate that is being tested if it falls inside the polygon. + y : double + Y coordinate that is being tested if it falls inside the polygon. + left : double + Value defining the left side of the clipped polygon + right : double + Value defining the right side of the clipped polygon + top : double + Value defining the top edge of the clipped polygon + bottom : double + Value defining the bottom edge of the clipped polygon + */ switch(edge) { case CP_LEFT: @@ -81,9 +131,36 @@ int insideWindow(int edge, double x, double y, return 0; } -// support function for sh_find_overlap + int calcCondition(int edge, double x1, double y1, double x2, double y2, double left, double right, double top, double bottom) { + + /* A support function for sh_find_overlap. Is a point in the polygon ? + + Calls the insideWindow routine to determine if x1,y1 from pixel j and/or + x2,y2 from pixel j+1 is inside the clipped polygon. Based on this information + the vertices of the clipped polygon are updated. + + edge : int + Integer value which defines which edge of the polygon we are working with. + 0 = left, 1 = right, 2 = top, 3 = bottom. + x1 : double + X value of a pixel[j] + y1 : double + Y value of a pixel[j] + x2 : double + X value of a pixel[j+1] + y2 : double + Y value of a pixel[j+1] + left : double + Value defining the left side of the clipped polygon + right : double + Value defining the right side of the clipped polygon + top : double + Value defining the top edge of the clipped polygon + bottom : double + Value defining the bottom edge of the clipped polygon + */ int stat1 = insideWindow(edge,x1,y1,left,right,top,bottom); int stat2 = insideWindow(edge,x2,y2,left,right,top,bottom); if(!stat1 && stat2) return 1; @@ -94,10 +171,40 @@ int calcCondition(int edge, double x1, double y1, double x2, double y2, } -// support function for sh_find_overlap -void solveIntersection(int edge ,double x1,double y1,double x2,double y2, - double *x,double *y, + +void solveIntersection(int edge, double x1, double y1, double x2, double y2, + double *x, double *y, double left, double right, double top, double bottom){ + + /* + A support function for sh_find_overlap. Find the intersection of a polygon and + a cube spaxel. The cube spaxel is defined on an evenly spaced regular grid. + + edge : int + Integer value that defines which edge of the polygon we are working with. + 0 = left, 1 = right, 2 = top, 3 = bottom. + x1 : double + X value of a pixel[j] + y1 : double + Y value of a pixel[j] + x2 : double + X value of a pixel[j+1] + y2 : double + Y value of a pixel[j+1] + x : double + The return x value of the x vertice of the clipped polygon + y : double + The return y value of the x vertice of the clipped polygon + left : double + Value defining the left side of the clipped polygon + right : double + Value defining the right side of the clipped polygon + top : double + Value defining the top edge of the clipped polygon + bottom : double + Value defining the bottom edge of the clipped polygon + + */ float m = 0; if(x2 != x1) m = ((double)(y2-y1)/(double)(x2-x1)); switch(edge) @@ -128,7 +235,7 @@ void solveIntersection(int edge ,double x1,double y1,double x2,double y2, } double find_area_quad(double MinX, double MinY, double Xcorner[], double Ycorner[]){ - /* Find the area of an quadrilateral + /* Find the area of an quadrilateral between clipped polygon and cube spaxel. Parameters ---------- @@ -136,14 +243,15 @@ double find_area_quad(double MinX, double MinY, double Xcorner[], double Ycorner Minimum X value MinY : float Minimum Y value - Xcorners : numpy.ndarray - x corner values (use first 4 corners) - YCorners : numpy.ndarray - y corner values (use first 4 corners) + Xcorners : ndarray + X corner values of a cube spaxel + YCorners : ndarray + Y corner values of a cube spaxe Returns ------- - Area + Area : double + Area of the overlap */ double PX[5]; @@ -171,9 +279,23 @@ double find_area_quad(double MinX, double MinY, double Xcorner[], double Ycorner } -// find the area of a closed polygon +// double find_area_poly(int nVertices,double xPixel[],double yPixel[]){ + /* + Find the area of a closed clipped polygon. + + nVertices : int + Number vertices of the clipped polygon. + xPixel : ndarray + X vertices of the clipped polygon + yPixel : ndarray + Y vertices of the clipped polygon + + Returns + areaPoly : double + Area of the clipped closed polygon + */ double area = 0; int i; double areaPoly = 0.0; @@ -199,9 +321,26 @@ double sh_find_overlap(const double xcenter, const double ycenter, const double xlength, const double ylength, double xPixelCorner[],double yPixelCorner[]) { - // Sutherland_hedgeman Polygon Clipping Algorithm to solve the overlap region - // between a 2-D detector mapped to IFU space with cube spaxel - + /*Sutherland-Hodgman Polygon Clipping Algorithm to solve the overlap region + between a 2-D detector mapped and the spatial plane of the IFU spaxel. + + xcenter : double + X center of the cube spaxel + ycenter : double + Y center of the cube spaxel + xlength : double + X length of the cube spaxel + ylength : double + Y length of the cube spaxel + xPixelCorner : ndarray + X corners of the polygon that is the formed by overlap between detector pixel and cube spaxel + yPixelCorner : ndarray + Y corners of the polygon that is the formed by overlap between detector pixel and cube spaxel + + Returns: + areaClipped : double + Area over the overlap clipped polygon + */ int i,j,k; int nVertices2 = 0; int condition; diff --git a/jwst/cube_build/tests/test_configuration.py b/jwst/cube_build/tests/test_configuration.py index b98fa754fb..931b44100b 100644 --- a/jwst/cube_build/tests/test_configuration.py +++ b/jwst/cube_build/tests/test_configuration.py @@ -7,7 +7,7 @@ from jwst.cube_build import cube_build from jwst.cube_build import file_table - +from jwst.cube_build.cube_build import NoFiltersError wcsinfo = { 'dec_ref': -0.00244536159612126, @@ -208,7 +208,7 @@ def test_calspec2_config(tmp_cwd, miri_ifushort_short): pars_input['subchannel'] = [] pars_input['filter'] = [] pars_input['grating'] = [] - weighting = 'msm' + weighting = 'drizzle' output_type = 'multi' # calspec 2 setup. Only 1 cube create from 2 channels single = False par_filename = 'None' @@ -255,7 +255,7 @@ def test_calspec3_config_miri(tmp_cwd, miri_full_coverage): pars_input['subchannel'] = [] pars_input['filter'] = [] pars_input['grating'] = [] - weighting = 'msm' + weighting = 'drizzle' output_type = 'band' single = False par_filename = 'None' @@ -321,14 +321,14 @@ def test_calspec3_config_miri(tmp_cwd, miri_full_coverage): def test_calspec3_config_miri_multi(tmp_cwd, miri_full_coverage): - """ Test CalSpec3 MIRI configuration default band cubes""" + """ Test CalSpec3 MIRI configuration default output type = multi""" pars_input = {} pars_input['channel'] = [] pars_input['subchannel'] = [] pars_input['filter'] = [] pars_input['grating'] = [] - weighting = 'msm' + weighting = 'drizzle' output_type = 'multi' single = False par_filename = 'None' @@ -356,6 +356,7 @@ def test_calspec3_config_miri_multi(tmp_cwd, miri_full_coverage): cubeinfo.instrument = this_instrument cubeinfo.determine_band_coverage(master_table) num_cubes, cube_pars = cubeinfo.number_cubes() + assert num_cubes == 1 assert cubeinfo.all_channel == ['1', '1', '1', '2', '2', '2', '3', '3', '3', '4', '4', '4'] @@ -374,6 +375,52 @@ def test_calspec3_config_miri_multi(tmp_cwd, miri_full_coverage): 'short', 'medium', 'long'] +def test_calspec3_config_miri_multi_ch1(tmp_cwd, miri_full_coverage): + """ Test CalSpec3 MIRI configuration output = band channel = 1""" + + pars_input = {} + pars_input['channel'] = [] + pars_input['subchannel'] = [] + pars_input['filter'] = [] + pars_input['grating'] = [] + weighting = 'drizzle' + output_type = 'multi' + channel = '1' + single = False + par_filename = 'None' + + pars = { + 'channel': pars_input['channel'], + 'subchannel': pars_input['subchannel'], + 'grating': pars_input['grating'], + 'filter': pars_input['filter'], + 'weighting': weighting, + 'single': single, + 'channel': channel, + 'output_type': output_type} + + cubeinfo = cube_build.CubeData( + miri_full_coverage, + par_filename, + **pars) + + master_table = file_table.FileTable() + this_instrument = master_table.set_file_table( + cubeinfo.input_models) + + assert this_instrument == 'MIRI' + + cubeinfo.instrument = this_instrument + cubeinfo.determine_band_coverage(master_table) + num_cubes, cube_pars = cubeinfo.number_cubes() + assert num_cubes == 1 + assert cubeinfo.all_channel == ['1', '1', '1'] + assert cubeinfo.all_subchannel == ['short', 'medium', 'long'] + + assert cube_pars['1']['par1'] == ['1', '1', '1'] + assert cube_pars['1']['par2'] == ['short', 'medium', 'long'] + + def test_calspec3_config_nirspec(tmp_cwd, nirspec_medium_coverage): """ Test CalSpec3 configuration for NIRSpec - default band cubes""" @@ -382,7 +429,7 @@ def test_calspec3_config_nirspec(tmp_cwd, nirspec_medium_coverage): pars_input['subchannel'] = [] pars_input['filter'] = [] pars_input['grating'] = [] - weighting = 'msm' + weighting = 'drizzle' output_type = 'band' single = False par_filename = 'None' @@ -465,3 +512,96 @@ def test_calspec3_config_nirspec_multi(tmp_cwd, nirspec_medium_coverage): assert cube_pars['1']['par1'] == ['g140m', 'g235m'] assert cube_pars['1']['par2'] == ['f100lp', 'f170lp'] + + +def test_calspec3_config_nirspec_grating(tmp_cwd, nirspec_medium_coverage): + """ Test CalSpec3 configuration for NIRSpec - select grating""" + + pars_input = {} + pars_input['channel'] = [] + pars_input['subchannel'] = [] + pars_input['filter'] = [] + pars_input['grating'] = [] + weighting = 'drizzle' + grating = 'g140m' + filter = 'f100lp' + output_type = 'multi' + single = False + par_filename = 'None' + + pars = { + 'channel': pars_input['channel'], + 'subchannel': pars_input['subchannel'], + 'grating': pars_input['grating'], + 'filter': pars_input['filter'], + 'weighting': weighting, + 'grating': grating, + 'filter': filter, + 'single': single, + 'output_type': output_type} + + cubeinfo = cube_build.CubeData( + nirspec_medium_coverage, + par_filename, + **pars) + + master_table = file_table.FileTable() + this_instrument = master_table.set_file_table( + cubeinfo.input_models) + + assert this_instrument == 'NIRSPEC' + + cubeinfo.instrument = this_instrument + cubeinfo.determine_band_coverage(master_table) + num_cubes, cube_pars = cubeinfo.number_cubes() + + assert num_cubes == 1 + assert cubeinfo.all_grating == ['g140m'] + assert cubeinfo.all_filter == ['f100lp'] + + assert cube_pars['1']['par1'] == ['g140m'] + assert cube_pars['1']['par2'] == ['f100lp'] + + + +def test_calspec3_config_nirspec_no_grating(tmp_cwd, nirspec_medium_coverage): + """ Test CalSpec3 configuration for NIRSpec - select grating incorrect""" + + pars_input = {} + pars_input['channel'] = [] + pars_input['subchannel'] = [] + pars_input['filter'] = [] + pars_input['grating'] = [] + weighting = 'drizzle' + grating = 'g140h' + filter = 'f100lp' + output_type = 'multi' + single = False + par_filename = 'None' + + pars = { + 'channel': pars_input['channel'], + 'subchannel': pars_input['subchannel'], + 'grating': pars_input['grating'], + 'filter': pars_input['filter'], + 'weighting': weighting, + 'grating': grating, + 'filter': filter, + 'single': single, + 'output_type': output_type} + + + cubeinfo = cube_build.CubeData( + nirspec_medium_coverage, + par_filename, + **pars) + + master_table = file_table.FileTable() + this_instrument = master_table.set_file_table( + cubeinfo.input_models) + + assert this_instrument == 'NIRSPEC' + cubeinfo.instrument = this_instrument + with pytest.raises(NoFiltersError): + cubeinfo.determine_band_coverage(master_table) + diff --git a/jwst/cube_build/tests/test_cube_build_step.py b/jwst/cube_build/tests/test_cube_build_step.py index f2b396c723..3d1ad674ea 100644 --- a/jwst/cube_build/tests/test_cube_build_step.py +++ b/jwst/cube_build/tests/test_cube_build_step.py @@ -12,8 +12,8 @@ from jwst import assign_wcs from jwst.cube_build import CubeBuildStep -from jwst.cube_build.file_table import ErrorNoAssignWCS -from jwst.cube_build.cube_build import ErrorNoChannels +from jwst.cube_build.file_table import NoAssignWCSError +from jwst.cube_build.cube_build import NoChannelsError from jwst.datamodels import ModelContainer @@ -96,7 +96,7 @@ def miri_image(): image = IFUImageModel((20, 20)) image.data = np.random.random((20, 20)) image.meta.instrument.name = 'MIRI' - image.meta.instrument.detector = 'MIRIFULONG' + image.meta.instrument.detector = 'MIRIFUSHORT' image.meta.exposure.type = 'MIR_MRS' image.meta.instrument.channel = '12' image.meta.instrument.band = 'SHORT' @@ -118,22 +118,22 @@ def test_call_cube_build(tmp_cwd, miri_cube_pars, miri_image, tmp_path, as_filen # the image needs to be a full image and this take too much time # in a unit test - # Test ErrorNoAssignWCS is raised - with pytest.raises(ErrorNoAssignWCS): + # Test NoAssignWCSError is raised + with pytest.raises(NoAssignWCSError): step = CubeBuildStep() step.override_cubepar = miri_cube_pars step.channel = '3' step.run(step_input) # Test some defaults to step are setup correctly and - # is user specifies channel is set up correctly + # if user specifies channel it is set up correctly step = CubeBuildStep() step.override_cubepar = miri_cube_pars step.channel = '1' try: step.run(step_input) - except ErrorNoAssignWCS: + except NoAssignWCSError: pass assert step.pars_input['channel'] == ['1'] @@ -146,7 +146,7 @@ def test_call_cube_build(tmp_cwd, miri_cube_pars, miri_image, tmp_path, as_filen # save file with modifications if as_filename: miri_image.save(step_input) - with pytest.raises(ErrorNoChannels): + with pytest.raises(NoChannelsError): step = CubeBuildStep() step.override_cubepar = miri_cube_pars step.channel = '3' @@ -223,3 +223,4 @@ def test_call_cube_build_nirspec_multi(tmp_cwd, nirspec_data, tmp_path, as_filen model = result[0] assert model.meta.cal_step.cube_build == 'COMPLETE' assert model.meta.filename == 'test_nirspec_s3d.fits' + diff --git a/jwst/cube_build/tests/test_miri_cubepars.py b/jwst/cube_build/tests/test_miri_cubepars.py index b51f8d919c..c485927db0 100644 --- a/jwst/cube_build/tests/test_miri_cubepars.py +++ b/jwst/cube_build/tests/test_miri_cubepars.py @@ -109,18 +109,18 @@ def test_miri_use_cubepars(tmp_cwd, miri_cube_pars): par1 = '1' par2 = 'medium' - ascale, bscale, wscale = instrument_info.GetScale(par1, par2) + ascale, bscale, wscale = instrument_info.get_scale(par1, par2) # check that the values are read in correctly assert math.isclose(ascale, 0.13, abs_tol=0.00001) assert math.isclose(bscale, 0.13, abs_tol=0.00001) assert math.isclose(wscale, 0.001, abs_tol=0.00001) - roiw = instrument_info.GetWaveRoi(par1, par2) - rois = instrument_info.GetSpatialRoi(par1, par2) - power = instrument_info.GetMSMPower(par1, par2) - wavemin = instrument_info.GetWaveMin(par1, par2) - wavemax = instrument_info.GetWaveMax(par1, par2) + roiw = instrument_info.get_wave_roi(par1, par2) + rois = instrument_info.get_spatial_roi(par1, par2) + power = instrument_info.get_msm_power(par1, par2) + wavemin = instrument_info.get_wave_min(par1, par2) + wavemax = instrument_info.get_wave_max(par1, par2) assert math.isclose(roiw, 0.001, abs_tol=0.00001) assert math.isclose(rois, 0.1, abs_tol=0.00001) @@ -208,17 +208,17 @@ def test_miri_cubepars_user_defaults(tmp_cwd, miri_cube_pars): # first check that it reads in correct values for this band # from the reference file - ascale, bscale, wscale = instrument_info.GetScale(par1, par2) + ascale, bscale, wscale = instrument_info.get_scale(par1, par2) assert math.isclose(ascale, 0.35, abs_tol=0.00001) assert math.isclose(bscale, 0.35, abs_tol=0.00001) assert math.isclose(wscale, 0.006, abs_tol=0.00001) - roiw = instrument_info.GetWaveRoi(par1, par2) - rois = instrument_info.GetSpatialRoi(par1, par2) - power = instrument_info.GetMSMPower(par1, par2) - wavemin = instrument_info.GetWaveMin(par1, par2) - wavemax = instrument_info.GetWaveMax(par1, par2) + roiw = instrument_info.get_wave_roi(par1, par2) + rois = instrument_info.get_spatial_roi(par1, par2) + power = instrument_info.get_msm_power(par1, par2) + wavemin = instrument_info.get_wave_min(par1, par2) + wavemax = instrument_info.get_wave_max(par1, par2) assert math.isclose(roiw, 0.006, abs_tol=0.00001) assert math.isclose(rois, 0.4, abs_tol=0.00001) @@ -371,17 +371,17 @@ def test_miri_cubepars_multiple_bands(tmp_cwd, miri_cube_pars): # first check that it reads in correct values for this band # from the reference file - ascale, bscale, wscale = instrument_info.GetScale(par1, par2) + ascale, bscale, wscale = instrument_info.get_scale(par1, par2) assert math.isclose(ascale, 0.2, abs_tol=0.00001) assert math.isclose(bscale, 0.2, abs_tol=0.00001) assert math.isclose(wscale, 0.003, abs_tol=0.00001) - roiw = instrument_info.GetWaveRoi(par1, par2) - rois = instrument_info.GetSpatialRoi(par1, par2) - power = instrument_info.GetMSMPower(par1, par2) - wavemin = instrument_info.GetWaveMin(par1, par2) - wavemax = instrument_info.GetWaveMax(par1, par2) + roiw = instrument_info.get_wave_roi(par1, par2) + rois = instrument_info.get_spatial_roi(par1, par2) + power = instrument_info.get_msm_power(par1, par2) + wavemin = instrument_info.get_wave_min(par1, par2) + wavemax = instrument_info.get_wave_max(par1, par2) assert math.isclose(roiw, 0.003, abs_tol=0.00001) assert math.isclose(rois, 0.2, abs_tol=0.00001) diff --git a/jwst/cube_build/tests/test_nirspec_cubepars.py b/jwst/cube_build/tests/test_nirspec_cubepars.py index 40a38b2933..c730786a84 100644 --- a/jwst/cube_build/tests/test_nirspec_cubepars.py +++ b/jwst/cube_build/tests/test_nirspec_cubepars.py @@ -137,16 +137,16 @@ def test_nirspec_cubepars(tmp_cwd, nirspec_cube_pars): par1 = 'prism' par2 = 'clear' - ascale, bscale, wscale = instrument_info.GetScale(par1, par2) + ascale, bscale, wscale = instrument_info.get_scale(par1, par2) assert math.isclose(ascale, 0.1, abs_tol=0.00001) assert math.isclose(bscale, 0.1, abs_tol=0.00001) assert math.isclose(wscale, 0.005, abs_tol=0.00001) - roiw = instrument_info.GetWaveRoi(par1, par2) - rois = instrument_info.GetSpatialRoi(par1, par2) - power = instrument_info.GetMSMPower(par1, par2) - wavemin = instrument_info.GetWaveMin(par1, par2) - wavemax = instrument_info.GetWaveMax(par1, par2) + roiw = instrument_info.get_wave_roi(par1, par2) + rois = instrument_info.get_spatial_roi(par1, par2) + power = instrument_info.get_msm_power(par1, par2) + wavemin = instrument_info.get_wave_min(par1, par2) + wavemax = instrument_info.get_wave_max(par1, par2) assert math.isclose(roiw, 0.011, abs_tol=0.00001) assert math.isclose(rois, 0.201, abs_tol=0.00001) diff --git a/jwst/cube_build/tests/test_wcs.py b/jwst/cube_build/tests/test_wcs.py index b319444245..7b5d166df1 100644 --- a/jwst/cube_build/tests/test_wcs.py +++ b/jwst/cube_build/tests/test_wcs.py @@ -245,10 +245,10 @@ def test_footprint_miri(): this_channel = '3' coord_system = 'skyalign' instrument_info = instrument_defaults.InstrumentInfo() - instrument_info.SetXSliceLimits(0, 101, this_channel) - x1, x2 = instrument_info.GetMIRISliceEndPts(this_channel) + instrument_info.set_xslice_limits(0, 101, this_channel) + x1, x2 = instrument_info.get_miri_slice_endpts(this_channel) - corners = cube_build_wcs_util.find_corners_MIRI(input_model, + corners = cube_build_wcs_util.find_corners_miri(input_model, this_channel, instrument_info, coord_system)