diff --git a/CHANGES.rst b/CHANGES.rst index 9d430d4438..75917d1397 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -17,6 +17,8 @@ cube_build ---------- - Fix bug when creating cubes using output_type=channel. [#6138] +- Move computationally intensive routines to c extensions and + removed miri psf weight function. [#6093] datamodels ---------- diff --git a/jwst/cube_build/blot_cube.py b/jwst/cube_build/blot_cube.py new file mode 100644 index 0000000000..fe19f855b0 --- /dev/null +++ b/jwst/cube_build/blot_cube.py @@ -0,0 +1,50 @@ +""" Map the detector pixels to the cube coordinate system. +This is where the weight functions are used. +""" +import numpy as np +import logging +log = logging.getLogger(__name__) +log.setLevel(logging.DEBUG) + + +def blot_overlap(ipt, xstart, + xcenter, ycenter, + x_cube, y_cube, + flux_cube, + blot_xsize, + blot_flux, blot_weight): + + """ Blot the median sky image back to the detector + + ipt is the median element. + xcenter, ycenter are the detector pixel arrays + xstart is only valid for MIRI and is the start of the x detector value for channel + 0 for channels on left and ~512 for channels on right + x_cube, y_cube: median cube IFU mapped backwards to detector + flux_cube: median flux + blot_flux & blot_weight: blotted values of flux_cube to detector + """ + + roi_det = 1.0 # Just large enough that we don't get holes + xdistance = np.absolute(x_cube[ipt] - xcenter) + ydistance = np.absolute(y_cube[ipt] - ycenter) + + index_x = np.where(xdistance <= roi_det) + index_y = np.where(ydistance <= roi_det) + + if len(index_x[0]) > 0 and len(index_y[0]) > 0: + + d1pix = x_cube[ipt] - xcenter[index_x] + d2pix = y_cube[ipt] - ycenter[index_y] + + dxy = [(dx * dx + dy * dy) for dy in d2pix for dx in d1pix] + dxy = np.array(dxy) + dxy = np.sqrt(dxy) + weight_distance = np.exp(-dxy) + weighted_flux = weight_distance * flux_cube[ipt] + + index2d = [iy * blot_xsize + ix for iy in index_y[0] for ix in (index_x[0] + xstart)] + index2d = np.array(index2d) + + blot_flux[index2d] = blot_flux[index2d] + weighted_flux + blot_weight[index2d] = blot_weight[index2d] + weight_distance diff --git a/jwst/cube_build/blot_cube_build.py b/jwst/cube_build/blot_cube_build.py index b0fcff326f..281a861e4b 100644 --- a/jwst/cube_build/blot_cube_build.py +++ b/jwst/cube_build/blot_cube_build.py @@ -8,6 +8,8 @@ from gwcs import wcstools from . import instrument_defaults +from . import blot_cube + log = logging.getLogger(__name__) log.setLevel(logging.DEBUG) @@ -199,16 +201,15 @@ def blot_images_miri(self): # pixel. # the Regular grid is on the x,y detector - blot_flux = self.blot_overlap_quick(model, xcenter, ycenter, - xstart, x_cube, y_cube, - flux_cube) + blot_flux = self.blot_overlap_miri(model, xcenter, ycenter, + xstart, x_cube, y_cube, + flux_cube) blot.data = blot_flux blot_models.append(blot) return blot_models - # ________________________________________________________________________________ - def blot_overlap_quick(self, model, xcenter, ycenter, xstart, - x_cube, y_cube, flux_cube): + def blot_overlap_miri(self, model, xcenter, ycenter, xstart, + x_cube, y_cube, flux_cube): # blot_overlap_quick finds to overlap between the blotted sky values # (x_cube, y_cube) and the detector pixels. @@ -223,31 +224,16 @@ def blot_overlap_quick(self, model, xcenter, ycenter, xstart, blot_weight = np.ndarray.flatten(blot_weight) # loop over valid points in cube (not empty edge pixels) - roi_det = 1.0 # Just large enough that we don't get holes - ivalid = np.nonzero(np.absolute(flux_cube)) + ivalid = np.nonzero(np.absolute(flux_cube)) t0 = time.time() for ipt in ivalid[0]: - # search xcenter and ycenter seperately. These arrays are smallsh. - # xcenter size = naxis1 on detector (for MIRI only 1/2 array) - # ycenter size = naxis2 on detector - xdistance = np.absolute(x_cube[ipt] - xcenter) - ydistance = np.absolute(y_cube[ipt] - ycenter) - - index_x = np.where(xdistance <= roi_det) - index_y = np.where(ydistance <= roi_det) - - if len(index_x[0]) > 0 and len(index_y[0]) > 0: - d1pix = np.array(x_cube[ipt] - xcenter[index_x]) - d2pix = np.array(y_cube[ipt] - ycenter[index_y]) - - dxy = [(dx * dx + dy * dy) for dy in d2pix for dx in d1pix] - dxy = np.sqrt(dxy) - weight_distance = np.exp(-dxy) - weighted_flux = weight_distance * flux_cube[ipt] - index2d = [iy * blot_xsize + ix for iy in index_y[0] for ix in (index_x[0] + xstart)] - blot_flux[index2d] = blot_flux[index2d] + weighted_flux - blot_weight[index2d] = blot_weight[index2d] + weight_distance + blot_cube.blot_overlap(ipt, xstart, + xcenter, ycenter, + x_cube, y_cube, + flux_cube, + blot_xsize, + blot_flux, blot_weight) t1 = time.time() log.debug(f"Time to blot median image to input model = {t1-t0:.1f}") @@ -257,7 +243,6 @@ def blot_overlap_quick(self, model, xcenter, ycenter, xstart, blot_flux = blot_flux.reshape((blot_ysize, blot_xsize)) return blot_flux - # ************************************************************************ def blot_images_nirspec(self): """ Core blotting routine for NIRSPEC @@ -299,7 +284,7 @@ def blot_images_nirspec(self): nslices = 30 log.info('Looping over 30 slices on NIRSPEC detector, this takes a little while') t0 = time.time() - roi_det = 1.0 # Just large enough that we don't get holes + for ii in range(nslices): ts0 = time.time() # for each slice pull out the blotted values that actually fall on the slice region @@ -338,30 +323,17 @@ def blot_images_nirspec(self): y_slice = y_slice[fuse] flux_slice = flux_slice[fuse] + xstart = 0 nn = flux_slice.size for ipt in range(nn): - # search xcenter and ycenter seperately. These arrays are smallish. - # xcenter size = naxis1 on detector - # ycenter size = naxis2 on detector - xdistance = np.absolute(x_slice[ipt] - xcenter) - ydistance = np.absolute(y_slice[ipt] - ycenter) - - index_x = np.where(xdistance <= roi_det) - index_y = np.where(ydistance <= roi_det) - - if len(index_x[0]) > 0 and len(index_y[0]) > 0: - d1pix = np.array(x_slice[ipt] - xcenter[index_x]) - d2pix = np.array(y_slice[ipt] - ycenter[index_y]) - - dxy = [(dx * dx + dy * dy) for dy in d2pix for dx in d1pix] - dxy = np.sqrt(dxy) - weight_distance = np.exp(-dxy) - weighted_flux = weight_distance * flux_slice[ipt] - index2d = [iy * blot_xsize + ix for iy in index_y[0] for ix in (index_x[0])] - blot_flux[index2d] = blot_flux[index2d] + weighted_flux - blot_weight[index2d] = blot_weight[index2d] + weight_distance + blot_cube.blot_overlap(ipt, xstart, + xcenter, ycenter, + x_slice, y_slice, + flux_slice, + blot_xsize, + blot_flux, blot_weight) ts1 = time.time() - log.debug(f"Time to map 1 slice = {ts1-ts0:.1f}") + log.debug(f"Time to blot 1 slice on NIRspec = {ts1-ts0:.1f}") # done mapping median cube to this input model t1 = time.time() log.debug(f"Time to blot median image to input model = {t1-t0:.1f}") diff --git a/jwst/cube_build/cube_build.py b/jwst/cube_build/cube_build.py index 5d75a4e555..4878c1bbc4 100644 --- a/jwst/cube_build/cube_build.py +++ b/jwst/cube_build/cube_build.py @@ -18,7 +18,6 @@ def __init__(self, input_models, input_filenames, par_filename, - resol_filename, **pars): """ Initialize the high level of information for the ifu cube @@ -36,15 +35,12 @@ def __init__(self, list of fits filenames par_filename: str cube parameter reference filename - resol_filename: str - miri resolution reference filename pars : dictionary holding top level cube parameters """ self.input_models = input_models self.input_filenames = input_filenames self.par_filename = par_filename - self.resol_filename = resol_filename self.single = pars.get('single') self.channel = pars.get('channel') self.subchannel = pars.get('subchannel') @@ -73,7 +69,6 @@ def setup(self): Read in necessary reference data: * cube parameter reference file - * if miripsf weighting parameter is set then read in resolution file This routine fills in the instrument_info dictionary, which holds the default spatial and spectral size of the output cube, as well as, @@ -117,15 +112,6 @@ def setup(self): self.all_filter, instrument_info) # ------------------------------------------------------------------------------- -# Read the miri resolution reference file - if self.weighting == 'miripsf': - log.info('Reading default MIRI cube resolution file %s', - self.resol_filename) - cube_build_io_util.read_resolution_file(self.resol_filename, - self.all_channel, - self.all_subchannel, - instrument_info) -# _______________________________________________________________________________ self.instrument_info = instrument_info # _______________________________________________________________________________ # Set up values to return and acess for other parts of cube_build diff --git a/jwst/cube_build/cube_build_step.py b/jwst/cube_build/cube_build_step.py index 30ec223985..fd8551b07d 100755 --- a/jwst/cube_build/cube_build_step.py +++ b/jwst/cube_build/cube_build_step.py @@ -41,7 +41,7 @@ class CubeBuildStep (Step): scale1 = float(default=0.0) # cube sample size to use for axis 1, arc seconds scale2 = float(default=0.0) # cube sample size to use for axis 2, arc seconds scalew = float(default=0.0) # cube sample size to use for axis 3, microns - weighting = option('emsm','msm','miripsf',default = 'emsm') # Type of weighting function + weighting = option('emsm','msm',default = 'emsm') # Type of weighting function coord_system = option('skyalign','world','internal_cal','ifualign',default='skyalign') # Output Coordinate system. rois = float(default=0.0) # region of interest spatial size, arc seconds roiw = float(default=0.0) # region of interest wavelength size, microns @@ -57,7 +57,7 @@ class CubeBuildStep (Step): output_use_model = boolean(default=true) # Use filenames in the output models """ - reference_file_types = ['cubepar', 'resol'] + reference_file_types = ['cubepar'] # ________________________________________________________________________________ def process(self, input): @@ -140,7 +140,7 @@ def process(self, input): # if interpolation is point cloud then weighting can be # 1. MSM: modified shepard method # 2. EMSM - # 3. miripsf - weighting for MIRI based on PSF and LSF + if self.coord_system == 'skyalign': self.interpolation = 'pointcloud' @@ -221,16 +221,6 @@ def process(self, input): self.log.warning('No default cube parameters reference file found') return # ________________________________________________________________________________ -# If miripsf weight is set then set up reference file - resol_filename = None - if self.weighting == 'miripsf': - resol_filename = self.get_reference_file(self.input_models[0], 'resol') - self.log.info(f'MIRI resol reference file {resol_filename}') - if resol_filename == 'N/A': - self.log.warning('No spectral resolution reference file found') - self.log.warning('Run again and turn off miripsf') - return -# ________________________________________________________________________________ # shove the input parameters in to pars to pull out in general cube_build.py pars = { @@ -269,7 +259,6 @@ def process(self, input): self.input_models, self.input_filenames, par_filename, - resol_filename, **pars) # ________________________________________________________________________________ # cubeinfo.setup: diff --git a/jwst/cube_build/cube_build_wcs_util.py b/jwst/cube_build/cube_build_wcs_util.py index ebff509fbf..6bebd30a9f 100644 --- a/jwst/cube_build/cube_build_wcs_util.py +++ b/jwst/cube_build/cube_build_wcs_util.py @@ -111,27 +111,19 @@ def find_corners_MIRI(input, this_channel, instrument_info, coord_system): lambda_min = np.nanmin(lam) lambda_max = np.nanmax(lam) - # before returning, ra should be between 0 to 360 - if a_min < 0: - a_min = a_min + 360 - if a_max >= 360.0: - a_max = a_max - 360.0 - - if a1 < 0: - a1 = a1 + 360 - if a1 > 360.0: - a1 = a1 - 360.0 - - if a2 < 0: - a2 = a2 + 360 - if a2 > 360.0: - a2 = a2 - 360.0 + if coord_system != 'internal_cal': + # before returning, ra should be between 0 to 360 + a_min %= 360 + a_max %= 360 + + a1 %= 360 + a2 %= 360 return a_min, b1, a_max, b2, a1, b_min, a2, b_max, lambda_min, lambda_max # ***************************************************************************** -def find_corners_NIRSPEC(input, this_channel, instrument_info, coord_system): +def find_corners_NIRSPEC(input, instrument_info, coord_system): """Find the sky footprint of a slice of a NIRSpec exposure For each slice find: @@ -241,5 +233,5 @@ def find_corners_NIRSPEC(input, this_channel, instrument_info, coord_system): lambda_min = min(lambda_slice) lambda_max = max(lambda_slice) + return a_min, b1, a_max, b2, a1, b_min, a2, b_max, lambda_min, lambda_max -# ______________________________________________________________________________ diff --git a/jwst/cube_build/cube_cloud.py b/jwst/cube_build/cube_cloud.py deleted file mode 100644 index e8a0f133d7..0000000000 --- a/jwst/cube_build/cube_cloud.py +++ /dev/null @@ -1,382 +0,0 @@ -""" Map the detector pixels to the cube coordinate system. -This is where the weight functions are used. -""" -import numpy as np -import logging - -log = logging.getLogger(__name__) -log.setLevel(logging.DEBUG) - - -def match_det2cube_msm(naxis1, naxis2, naxis3, - cdelt1, cdelt2, - zcdelt3, - xcenters, ycenters, zcoord, - spaxel_flux, - spaxel_weight, - spaxel_iflux, - spaxel_var, - flux, - err, - coord1, coord2, wave, - weighting_type, - rois_pixel, roiw_pixel, weight_pixel, - softrad_pixel, scalerad_pixel, - cube_debug, debug_file): - """ Map the detector pixels to the cube spaxels using the MSM parameters - - - Match the Point Cloud members to the spaxel centers that fall in the ROI. - For each spaxel the coord1,coord1 and wave point cloud members are weighed - according to modified shepard method of inverse weighting based on the - distance between the point cloud member and the spaxel center. - - Note that this routine does NOT build the cube by looping over spaxels and - looking for pixels that contribute to those spaxels. The runtime is significantly - better to instead loop over pixels, and look for spaxels that they contribute to. - This way we can just keep a running sum of the weighted fluxes in a 1-d representation - of the output cube, which can be normalized by the weights at the end. - - Parameters - ---------- - naxis1 : int - size of the ifucube in 1st axis - naxis2 : int - size of the ifucube in 2nd axis - naxis3 : int - size of the ifucube in 3rd axis - cdelt1 : float - ifucube spaxel size in axis 1 dimension - cdelt2 : float - ifucube spaxel size in axis 2 dimension - cdelt3_normal :float - ifu spectral size at wavelength - rois_pixel : float - region of influence size in spatial dimension - roiw_pixel : float - region of influence size in spectral dimension - weight_power : float - msm weighting parameter - xcenter : numpy.ndarray - spaxel center locations 1st dimensions. - ycenter : numpy.ndarray - spaxel center locations 2nd dimensions. - zcoord : numpy.ndarray - spaxel center locations in 3rd dimensions - spaxel_flux : numpy.ndarray - contains the weighted summed detector fluxes that fall - within the roi - spaxel_weight : numpy.ndarray - contains the summed weights assocated with the detector fluxes - spaxel_iflux : numpy.ndarray - number of detector pixels falling with roi of spaxel center - spaxel_var: numpy.ndarray - contains the weighted summed variance within the roi - flux : numpy.ndarray - array of detector fluxes associated with each position in - coorr1, coord2, wave - err: numpy.ndarray - array of detector errors associated with each position in - coorr1, coord2, wave - coord1 : numpy.ndarray - contains the spatial coordinate for 1st dimension for the mapped - detector pixel - coord2 : numpy.ndarray - contains the spatial coordinate for 2nd dimension for the mapped - detector pixel - wave : numpy.ndarray - contains the spectral coordinate for the mapped detector pixel - - Returns - ------- - spaxel_flux, spaxel_weight, spaxel_ifux, and spaxel_var updated with the information - from the detector pixels that fall within the roi of the spaxel center. - """ - nplane = naxis1 * naxis2 - - # now loop over the pixel values for this region and find the spaxels that fall - # within the region of interest. - nn = coord1.size - # print('looping over n points mapping to cloud',nn) - # ________________________________________________________________________________ - for ipt in range(0, nn - 1): - # xcenters, ycenters is a flattened 1-D array of the 2 X 2 xy plane - # cube coordinates. - # find the spaxels that fall withing ROI of point cloud defined by - # coord1,coord2,wave - lower_limit = softrad_pixel[ipt] - xdistance = (xcenters - coord1[ipt]) - ydistance = (ycenters - coord2[ipt]) - radius = np.sqrt(xdistance * xdistance + ydistance * ydistance) - - indexr = np.where(radius <= rois_pixel[ipt]) - indexz = np.where(abs(zcoord - wave[ipt]) <= roiw_pixel[ipt]) - - # Find the cube spectral planes that this input point will contribute to - if len(indexz[0]) > 0: - d1 = np.array(coord1[ipt] - xcenters[indexr]) / cdelt1 - d2 = np.array(coord2[ipt] - ycenters[indexr]) / cdelt2 - d3 = np.array(wave[ipt] - zcoord[indexz]) / zcdelt3[indexz] - - dxy = (d1 * d1) + (d2 * d2) - - # shape of dxy is #indexr or number of overlaps in spatial plane - # shape of d3 is #indexz or number of overlaps in spectral plane - # shape of dxy_matrix & d3_matrix (#indexr, #indexz) - # rows = number of overlaps in spatial plane - # cols = number of overlaps in spectral plane - dxy_matrix = np.tile(dxy[np.newaxis].T, [1, d3.shape[0]]) - d3_matrix = np.tile(d3 * d3, [dxy_matrix.shape[0], 1]) - - # wdistance is now the spatial distance squared plus the spectral distance squared - wdistance = dxy_matrix + d3_matrix - if weighting_type == 'msm': - weight_distance = np.power(np.sqrt(wdistance), weight_pixel[ipt]) - weight_distance[weight_distance < lower_limit] = lower_limit - weight_distance = 1.0 / weight_distance - elif weighting_type == 'emsm': - weight_distance = np.exp(-wdistance / (scalerad_pixel[ipt] / cdelt1)) - - weight_distance = weight_distance.flatten('F') - weighted_flux = weight_distance * flux[ipt] - weighted_var = (weight_distance * err[ipt]) * (weight_distance * err[ipt]) - - # Identify all of the cube spaxels (ordered in a 1d vector) that this input point contributes to - icube_index = [iz * nplane + ir for iz in indexz[0] for ir in indexr[0]] - - if cube_debug in icube_index: - log.info('cube_debug %i %d %d', ipt, flux[ipt], weight_distance[icube_index.index(cube_debug)]) - - # Add the weighted flux and variance to running 1d cubes, along with the weights - # (for later normalization), and point count (for information) - spaxel_flux[icube_index] = spaxel_flux[icube_index] + weighted_flux - spaxel_weight[icube_index] = spaxel_weight[icube_index] + weight_distance - spaxel_iflux[icube_index] = spaxel_iflux[icube_index] + 1 - spaxel_var[icube_index] = spaxel_var[icube_index] + weighted_var - -# _______________________________________________________________________ - - -def match_det2cube_miripsf(alpha_resol, beta_resol, wave_resol, - naxis1, naxis2, naxis3, - xcenters, ycenters, zcoord, - spaxel_flux, - spaxel_weight, - spaxel_iflux, - spaxel_var, - spaxel_alpha, spaxel_beta, spaxel_wave, - flux, - err, - coord1, coord2, wave, alpha_det, beta_det, - weighting_type, - rois_pixel, roiw_pixel, - weight_pixel, - softrad_pixel, - scalerad_pixel): - """ Map the detector pixels to the cube spaxels using miri PSF weighting - - Map coordinates coord1,coord2, and wave of the point cloud to which - spaxels they overlap with in the ifucube. For each spaxel the coord1,coord2 - and wave point cloud members are weighting according to the miri psf and lsf. - The weighting function is based on the distance the point cloud member - and spaxel center in the alph-beta coordinate system. The alpha and beta - value of each point cloud member and spaxel center is passed to this - routine. - - Parameters - ---------- - alpha_resol : numpy.ndarray - alpha resolution table - beta_resol : numpy.ndarray - beta resolution table - wave_resol : numpy.ndarray - wavelength resolution table - naxis1 : int - size of the ifucube in 1st axis - naxis2 : int - size of the ifucube in 2nd axis - naxis3 : int - size of the ifucube in 3rd axis - xcenter : numpy.ndarray - spaxel center locations 1st dimensions. - ycenter : numpy.ndarray - spaxel center locations 2nd dimensions. - zcoord : numpy.ndarray - spaxel center locations in 3rd dimensions - spaxel_flux : numpy.ndarray - contains the weighted summed detector fluxes that fall withi the roi - spaxel_weight : numpy.ndarray - contains the summed weights assocated with the detector fluxes - spaxel_iflux : numpy.ndarray - number of detector pixels falling with roi of spaxel center - spaxel_var: numpy.ndarray - contains the weighted summed variance within the roi - flux : numpy.ndarray - array of detector fluxes associated with each position in - coorr1, coord2, wave - err: numpy.ndarray - array of detector errors associated with each position in - coorr1, coord2, wave - spaxel_alpha : numpy.ndarray - alpha value of spaxel centers - spaxel_beta : numpy.ndarray - beta value of spaxel centers - spaxel_wave : numpy.ndarray - alpha value of spaxel centers - coord1 : numpy.ndarray - contains the spatial coordinate for 1st dimension for the mapped - detector pixel - coord2 : numpy.ndarray - contains the spatial coordinate for 2nd dimension for the mapped - detector pixel - wave : numpy.ndarray - contains the spectral coordinate for the mapped detector pixel - alpha_det : numpy.ndarray - alpha coordinate of mapped detector pixels - beta_det : numpy.ndarray - beta coordinate of mapped detector pixels - rois_pixel : float - region of influence size in spatial dimension - roiw_pixel : float - region of influence size in spectral dimension - weight_power : float - msm weighting parameter - softrad_pxiel :float - weighting paramter - - Returns - ------- - spaxel_flux, spaxel_weight, spaxel_ifux, spaxel_var updated with the information - from the detector pixels that fall within the roi if the spaxel center. - - """ - - nplane = naxis1 * naxis2 - # now loop over the pixel values for this region and find the spaxels - # that fall within the region of interest. - nn = coord1.size - # print('looping over n points mapping to cloud',nn) - # _______________________________________________________________________ - for ipt in range(0, nn - 1): - lower_limit = softrad_pixel[ipt] - - # ________________________________________________________ - # if weight is miripsf -distances determined in alpha-beta - # coordinate system - - weights = FindNormalizationWeights(wave[ipt], - wave_resol, - alpha_resol, - beta_resol) - # ___________________________________________________________________ - # xcenters, ycenters is a flattened 1-D array of the 2 X 2 xy plane - # cube coordinates. - # find the spaxels that fall withing ROI of point cloud defined by - # coord1,coord2,wave - xdistance = (xcenters - coord1[ipt]) - ydistance = (ycenters - coord2[ipt]) - radius = np.sqrt(xdistance * xdistance + ydistance * ydistance) - - indexr = np.where(radius <= rois_pixel[ipt]) - indexz = np.where(abs(zcoord - wave[ipt]) <= roiw_pixel[ipt]) - - # _______________________________________________________________ - # TODO if this method is used replace two for loops with quicker - # list comprehension - # loop over the points in the ROI - for iz, zz in enumerate(indexz[0]): - istart = zz * nplane - for ir, rr in enumerate(indexr[0]): - cube_index = istart + rr - - alpha_distance = alpha_det[ipt] - spaxel_alpha[cube_index] - beta_distance = beta_det[ipt] - spaxel_beta[cube_index] - wave_distance = abs(wave[ipt] - spaxel_wave[cube_index]) - - xn = alpha_distance / weights[0] - yn = beta_distance / weights[1] - wn = wave_distance / weights[2] - - # only included the spatial dimensions - wdistance = (xn * xn + yn * yn + wn * wn) -# ________________________________________________________________________________ - # MSM weighting based on 1/r**power - if weighting_type == 'msm': - weight_distance = np.power(np.sqrt(wdistance), weight_pixel[ipt]) - if weight_distance < lower_limit: - weight_distance = lower_limit - weight_distance = 1.0 / weight_distance - elif weighting_type == 'emsm': - weight_distance = scalerad_pixel[ipt] * np.exp(1.0 / wdistance) - - weighted_flux = weight_distance * flux[ipt] - weighted_var = (weight_distance * err[ipt]) * (weight_distance * err[ipt]) - - spaxel_flux[cube_index] = spaxel_flux[cube_index] + weighted_flux - spaxel_weight[cube_index] = spaxel_weight[cube_index] + weight_distance - spaxel_iflux[cube_index] = spaxel_iflux[cube_index] + 1 - spaxel_var[cube_index] = spaxel_var[cube_index] + weighted_var -# _______________________________________________________________________ - - -def FindNormalizationWeights(wavelength, - wave_resol, - alpha_resol, - beta_resol): - """ Routine used in MIRI PSF weighting to normalize data - - - Normalize weighting to each point cloud that contributes to the - ROI of a spaxel. The normalization of weighting is determined from - width of PSF as wellas wavelength resolution - - Parameters - ---------- - wave_resol : numpy.ndarray - wavelength resolution array - alpha_resol : numpy.ndarray - alpha psf resolution array - beta_resol : numpy.ndarray - beta psf resolution array - - Returns - ------- - normalized weighting for 3 dimension - """ - alpha_weight = 1.0 - beta_weight = 1.0 - lambda_weight = 1.0 - - # alpha psf weighting - alpha_wave_cutoff = alpha_resol[0] - alpha_a_short = alpha_resol[1] - alpha_b_short = alpha_resol[2] - alpha_a_long = alpha_resol[3] - alpha_b_long = alpha_resol[4] - if wavelength < alpha_wave_cutoff: - alpha_weight = alpha_a_short + alpha_b_short * wavelength - else: - alpha_weight = alpha_a_long + alpha_b_long * wavelength - - # beta psf weighting - beta_wave_cutoff = beta_resol[0] - beta_a_short = beta_resol[1] - beta_b_short = beta_resol[2] - beta_a_long = beta_resol[3] - beta_b_long = beta_resol[4] - if wavelength < beta_wave_cutoff: - beta_weight = beta_a_short + beta_b_short * wavelength - else: - beta_weight = beta_a_long + beta_b_long * wavelength - - # wavelength weighting - wavecenter = wave_resol[0] - a_ave = wave_resol[1] - b_ave = wave_resol[2] - c_ave = wave_resol[3] - - wave_diff = wavelength - wavecenter - resolution = a_ave + b_ave * wave_diff + c_ave * wave_diff * wave_diff - lambda_weight = wavelength / resolution - weight = [alpha_weight, beta_weight, lambda_weight] - return weight diff --git a/jwst/cube_build/cube_cloud_quick.py b/jwst/cube_build/cube_cloud_quick.py deleted file mode 100644 index 4ce379abb7..0000000000 --- a/jwst/cube_build/cube_cloud_quick.py +++ /dev/null @@ -1,140 +0,0 @@ -# Routines used in Spectral Cube Building -import numpy as np -import math -import logging - -log = logging.getLogger(__name__) -log.setLevel(logging.DEBUG) -# ________________________________________________________________________________ - - -def match_det2cube_msm(naxis1, naxis2, naxis3, - cdelt1, cdelt2, cdelt3, - rois, roiw, - msm_weight_power, - xcoord, ycoord, zcoord, - spaxel_flux, - spaxel_weight, - spaxel_iflux, - spaxel_var, - flux, - err, - coord1, coord2, wave): - """ - Short Summary - ------------- - Match the Point Cloud members to the spaxel centers that fall in the ROI. - For each spaxel the coord1,coord1 and wave point cloud members are weighting - according to modified shepard method of inverse weighting based on the distance between - the point cloud member and the spaxel center. - - Parameters - ---------- - naxis1,naxis2,naxis3: size of the ifucube - cdelt1,cdelt2,cdelt3: ifucube spaxel size in the 3 dimensions - rois, roiw: region of influence size in spatial and spectral dimension - weight_power: msm weighting parameter - xcoord,ycoord: spaxel center locations in 1st and 2nd dimensions. These values 1 D - zcoord: spaxel center locations in 3rd dimensions - spaxel_flux: contains the weighted summed detector fluxes that fall withi the roi - spaxel_weight: contains the summed weights assocated with the detector fluxes - spaxel_iflux: number of detector pixels falling with roi of spaxel center - spaxel_var: contains the weighted summed detector variances that fall withi the roi - flux: array of detector fluxes associated with each position in coorr1, coord2, wave - err: array of detector errors associated with each position in coorr1, coord2, wave - coord1, coord2, wave - Returns - ------- - spaxel_flux, spaxel_weight, spaxel_ifux, spaxel_var updated with the information from the - detector pixels that fall within the roi if the spaxel center. - """ - - nplane = naxis1 * naxis2 - lower_limit = 0.01 - - # now loop over the pixel values for this region and find the spaxels that fall - # withing the region of interest. - nn = coord1.size - - # transform to ifu cube grid system - c1 = (coord1 - xcoord[0]) / cdelt1 - c2 = (coord2 - ycoord[0]) / cdelt2 - c3 = (wave - zcoord[0]) / cdelt3 - - spatial_box = math.ceil(rois / cdelt1) - spectral_radius = math.ceil(roiw / cdelt3) - - c1_min = (c1 - spatial_box).astype(int) - c1_max = (c1 + spatial_box).astype(int) - c1_min[c1_min < 0] = 0 - c1_max[c1_max >= naxis1] = naxis1 - 1 - - c2_min = (c2 - spatial_box).astype(int) - c2_min[c2_min < 0] = 0 - c2_max = (c2 + spatial_box).astype(int) - c2_max[c2_max >= naxis2] = naxis2 - 1 - - c3_min = (c3 - spectral_radius).astype(int) - c3_min[c3_min < 0] = 0 - c3_max = (c3 + spectral_radius).astype(int) - c3_max[c3_max >= naxis3] = naxis3 - 1 - - # print('looping over n points mapping to cloud',nn) - for ipt in range(0, nn - 1): - # xcoord is 1 d array of center of spaxel in x direction - # ycoord is 1 d array of center of spaxel in y direction - # find the spaxels that fall withing ROI of point cloud defined by - # coord1,coord2,wave - ixrange = np.arange(c1_min[ipt], c1_max[ipt] + 1) - iyrange = np.arange(c2_min[ipt], c2_max[ipt] + 1) - izrange = np.arange(c3_min[ipt], c3_max[ipt] + 1) - - xdistance = (xcoord[ixrange] - coord1[ipt]) - ydistance = (ycoord[iyrange] - coord2[ipt]) - zdistance = (zcoord[izrange] - wave[ipt]) - - # form 2D array of xy plane indices - ixindex = np.tile(ixrange[np.newaxis].T, [1, iyrange.shape[0]]) - iyindex = np.tile(iyrange, [ixindex.shape[0], 1]) - xyindex = (iyindex * naxis1) + ixindex - - # for 2D array of x distance y distance values - xvalues = np.tile(xdistance[np.newaxis].T, [1, ydistance.shape[0]]) - yvalues = np.tile(ydistance, [xvalues.shape[0], 1]) - radius = np.sqrt(xvalues * xvalues + yvalues * yvalues) - - # pull out values that fall in the roi - indexr = np.where(radius <= rois) - indexz = np.where(abs(zdistance) <= roiw) - - # grab the index of values that fall in the roi - xyindex_good = xyindex[indexr] - zindex_good = izrange[indexz] - - d1 = np.array(xvalues[indexr]) / cdelt1 - d2 = np.array(yvalues[indexr]) / cdelt2 - d3 = np.array(zdistance[indexz]) / cdelt3 - dxy = (d1 * d1) + (d2 * d2) - - # shape of dxy is #indexr or number of overlaps in spatial plane - # shape of d3 is #indexz or number of overlaps in spectral plane - # shape of dxy_matrix & d3_matrix (#indexr, #indexz) - # rows = number of overlaps in spatial plane - # cols = number of overlaps in spectral plane - dxy_matrix = np.tile(dxy[np.newaxis].T, [1, d3.shape[0]]) - d3_matrix = np.tile(d3 * d3, [dxy_matrix.shape[0], 1]) - - wdistance = dxy_matrix + d3_matrix - weight_distance = np.power(np.sqrt(wdistance), msm_weight_power) - weight_distance[weight_distance < lower_limit] = lower_limit - weight_distance = 1.0 / weight_distance - weight_distance = weight_distance.flatten('F') - weighted_flux = weight_distance * flux[ipt] - weighted_var = (weight_distance * err[ipt]) * (weight_distance * err[ipt]) - - icube_index = [iz * nplane + ixy for iz in zindex_good for ixy in xyindex_good] - spaxel_flux[icube_index] = spaxel_flux[icube_index] + weighted_flux - spaxel_weight[icube_index] = spaxel_weight[icube_index] + weight_distance - spaxel_iflux[icube_index] = spaxel_iflux[icube_index] + 1 - spaxel_var[icube_index] = spaxel_var[icube_index] + weighted_var -# _______________________________________________________________________ diff --git a/jwst/cube_build/cube_internal_cal.py b/jwst/cube_build/cube_internal_cal.py new file mode 100644 index 0000000000..c8ad539af1 --- /dev/null +++ b/jwst/cube_build/cube_internal_cal.py @@ -0,0 +1,143 @@ +""" Routines for creating single band, single exposure IFU Cubes with +the interoplation method = area , coord_system = internal_cal +""" +import numpy as np +from ..datamodels import dqflags +from jwst.transforms.models import _toindex +from .cube_match_internal import cube_wrapper_internal # c extension + + +def match_det2cube(instrument, + x, y, sliceno, + input_model, + transform, + spaxel_flux, + spaxel_weight, + spaxel_iflux, + spaxel_var, + 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. + 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 + sliceno : int + slice number + input_model : datamodel + input slope 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. + + Returns + ------- + spaxel filled in with needed information on overlapping detector pixels + """ + + x = _toindex(x) + y = _toindex(y) + + pixel_dq = input_model.dq[y, x] + + 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 holds the location of pixels we want to map to cube + x = x[good_data] + y = y[good_data] + + coord1, coord2, lam = transform(x, y) + valid = ~np.isnan(coord2) + x = x[valid] + y = y[valid] + + yy_bot = y + yy_top = y + 1 + xx_left = x + xx_right = x + 1 + # along slice dimension is second coordinate returned from transform + 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) + b4, a4, lam4 = transform(xx_left, yy_top) + # check if units are in microns or meters, if meters convert to microns + # only need to check one of the wavelenghts + lmax = np.nanmax(lam1.flatten()) + if lmax < 0.0001: + lam1 = lam1 * 1.0e6 + lam2 = lam2 * 1.0e6 + lam3 = lam3 * 1.0e6 + lam4 = lam4 * 1.0e6 + + 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) + a4, b4, lam4 = transform(xx_left, yy_top) + + # corners are returned Nanned if outside range of slice + # fixed the nanned corners when a1,lam1 is valid but + # adding 1 to x,y pushes data outside BB valid region + + # on the edge out of bounds + index_good2 = ~np.isnan(a2) + index_good3 = ~np.isnan(a3) + index_good4 = ~np.isnan(a4) + + # we need the cases of only all corners valid numbers + good = np.where(index_good2 & index_good3 & index_good4) + + a1 = a1[good] + a2 = a2[good] + a3 = a3[good] + a4 = a4[good] + lam1 = lam1[good] + lam2 = lam2[good] + lam3 = lam3[good] + lam4 = lam4[good] + x = x[good] + y = y[good] + + # center of first pixel, x,y = 1 for Adrian's equations + # but we want the pixel corners, x,y values passed into this + # routine to start at 0 + pixel_flux = input_model.data[y, x] + pixel_err = input_model.err[y, x] + + # 1-1 mapping in across slice direction (x for NIRSPEC, y for MIRI) + if instrument == 'NIRSPEC': + ss = sliceno + instrument_no = 1 + 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) + + spaxel_flux_slice, spaxel_weight_slice, spaxel_var_slice, spaxel_iflux_slice = result + spaxel_flux = spaxel_flux + np.asarray(spaxel_flux_slice, np.float64) + spaxel_weight = spaxel_weight + np.asarray(spaxel_weight_slice, np.float64) + spaxel_var = spaxel_var + np.asarray(spaxel_var_slice, np.float64) + spaxel_iflux = spaxel_iflux + np.asarray(spaxel_flux_slice,np.float64) + result = None diff --git a/jwst/cube_build/cube_overlap.py b/jwst/cube_build/cube_overlap.py deleted file mode 100644 index 825a08dfc5..0000000000 --- a/jwst/cube_build/cube_overlap.py +++ /dev/null @@ -1,596 +0,0 @@ -""" Routines for creating single band, single exposure IFU Cubes with -the interoplation method = area -""" -import numpy as np -import math -from ..datamodels import dqflags -from jwst.transforms.models import _toindex - - -def find_area_poly(nvertices, xpixel, ypixel): - """ Find the area of the polygon - - Parameters - ---------- - nvertices : int - number of Vertices of polygon - xpixel : numpy.ndarray - x coordinate of vertices - ypixel : numpy.ndarray - y coordinate of vertices - - Returns - ------- - area of polygon - - """ - areaPoly = 0.0 - xmin = min(xpixel) - ymin = min(ypixel) - - for i in range(0, nvertices - 1): - area = (xpixel[i] - xmin) * (ypixel[i + 1] - ymin) - (xpixel[i + 1] - xmin) * (ypixel[i] - ymin) - areaPoly = areaPoly + area - - areaPoly = abs(0.5 * areaPoly) - return areaPoly -# _____________________________________________________________________________ - - -def find_area_quad(MinX, MinY, Xcorner, Ycorner): - """ Find the area of an quadrilateral - - Parameters - ---------- - MinX : float - 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) - - Returns - ------- - Area - """ - PX = [] - PY = [] - - PX.append(Xcorner[0] - MinX) - PX.append(Xcorner[1] - MinX) - PX.append(Xcorner[2] - MinX) - PX.append(Xcorner[3] - MinX) - PX.append(PX[0]) - - PY.append(Ycorner[0] - MinY) - PY.append(Ycorner[1] - MinY) - PY.append(Ycorner[2] - MinY) - PY.append(Ycorner[3] - MinY) - PY.append(PY[0]) - - Area = 0.5 * ((PX[0] * PY[1] - PX[1] * PY[0]) + - (PX[1] * PY[2] - PX[2] * PY[1]) + - (PX[2] * PY[3] - PX[3] * PY[2]) + - (PX[3] * PY[4] - PX[4] * PY[3])) - - return abs(Area) -# _______________________________________________________________________ - - -def calcCondition(edge, x1, y1, x2, y2, left, right, top, bottom): - """ Determine if a point is inside a polygon - - Parameters - ---------- - edge : float - edge of spaxel - x1 : float - x min coordinate of pixel - y1 : float - y min coordinate of pixel - x2 : float - x max coordinate of pixel - y2 : float - y max coordinate of pixel - left : float - left side of spaxel - right : float - right side of spaxel - top : float - top of spaxel - bottom : float - bottom of spaxel - - Returns - ------- - where the detector pixel is in relation to a side of the spaxel - """ - - stat1 = insideWindow(edge, x1, y1, left, right, top, bottom) - stat2 = insideWindow(edge, x2, y2, left, right, top, bottom) - - if not stat1 and stat2: - return 1 - if stat1 and stat2: - return 2 - if stat1 and not stat2: - return 3 - if not stat1 and not stat2: - return 4 - return 0 # never executed -# _______________________________________________________________________ - - -def insideWindow(edge, x, y, left, right, top, bottom): - """Function used in determined overlap of detector pixel and spaxel - - Given the pixel edge and cener and the left,right,top bottom of - sides of spaxel return if detector edge is inside spaxel - - Parameters - ---------- - edge : float - edge of pixel - x : float - x center of pixel - y : float - y center of pixel - left : float - left side of spaxel - right : float - right side of spaxel - top : float - top of spaxel - bottom : float - bottom of spaxel - - Returns - ------- - returns true or false values if detector point is on "correct" - side of spaxel - """ - - CP_LEFT = 0 - CP_RIGHT = 1 - CP_BOTTOM = 2 - CP_TOP = 3 - - if edge == CP_LEFT: - return (x > left) - elif edge == CP_RIGHT: - return (x < right) - elif edge == CP_BOTTOM: - return (y > bottom) - elif edge == CP_TOP: - return (y < top) - else: - return 0 -# _______________________________________________________________________ - - -def solve_intersection(edge, x1, y1, x2, y2, - left, right, top, bottom): - """ Finds the intersection of a polygon and rectangular pixel - - Parameters - ---------- - edge : float - one of the 4 edges of spaxel - x1 : float - x min of pixel - y1 : float - y min of pixel - x2 : float - x max of pixel - y2 : float - y max of pixel - left : float - left side of spaxel - right : float - right side of spaxel - top : float - top of spaxel - bottom : float - bottom of spaxel - - Returns - ------- - returns x,y inside spaxel (detector detector region inside spaxel) - """ - x = 0 - y = 0 - CP_LEFT = 0 - CP_RIGHT = 1 - CP_BOTTOM = 2 - CP_TOP = 3 - m = 0 - if x2 != x1: - m = (y2 - y1) / (x2 - x1) - if edge == CP_LEFT: - x = left - y = y1 + m * (x - x1) - elif edge == CP_RIGHT: - x = right - y = y1 + m * (x - x1) - elif edge == CP_BOTTOM: - y = bottom - if(x1 != x2): - x = x1 + (1.0 / m) * (y - y1) - else: - x = x1 - elif edge == CP_TOP: - y = top - if(x1 != x2): - x = x1 + (1.0 / m) * (y - y1) - else: - x = x1 - return x, y -# _______________________________________________________________________ - - -def addpoint(x, y, xnew, ynew, nvertices2): - """ adds a point to vertices of the detector pixel region inside the spaxel - - Parameters - ---------- - x : float - x value to add - y : float - y value to add - xnew : numpy.ndarray - new x vertices - ynew : numpy.ndarray - new y vertices - nvertices2 : int - number of vertices - Returns - ------- - adds a vertice to the polygon describing the region of the detector pixel - inside the spaxel - """ - xnew[nvertices2] = x - ynew[nvertices2] = y - - nvertices2 = nvertices2 + 1 - - return nvertices2 -# ________________________________________________________________________________ - - -def sh_find_overlap(xcenter, ycenter, xlength, ylength, xp_corner, yp_corner): - """ Find overlap between pixel and spaxel - - Using the Sutherland_hedgeman Polygon Clipping Algorithm to solve the - overlap region first clip the x-y detector plane by the cube's xy rectangle - then find the overlap area - - Parameters - --------- - xcenter : float - center grid point in x dimension for cube (along slice- alpha) - ycenter : float - center grid point in y dimension for cube (lambda) - xlength : float - width of spaxel in x dimesion (along slice- alpha) - ylength : float - width of spaxel in y dimesion (lambda) - xp_corner : float - alpha pixel corner values - yp_Corner : float - lambda pixel corner values - - Returns - ------- - AreaOverlap - """ - - area_clipped = 0.0 - top = ycenter + 0.5 * ylength - bottom = ycenter - 0.5 * ylength - - left = xcenter - 0.5 * xlength - right = xcenter + 0.5 * xlength - - nvertices = 4 # input detector pixel vertices - max_vertices = 9 - # initialize xPixel, yPixel to the detector pixel corners. - # xPixel,yPixel will become the clipped polygon vertices - # inside the cube pixel - # xnew,ynew xpixel and ypixel of size max_vertices - - xPixel = [] - yPixel = [] - - xnew = [] - ynew = [] - - for j in range(0, 9): - xnew.append(0.0) - ynew.append(0.0) - xPixel.append(0.0) - yPixel.append(0.0) - - # Xpixel, YPixel closed (5 corners) - for i in range(0, 4): - xPixel[i] = xp_corner[i] - yPixel[i] = yp_corner[i] - xPixel[4] = xp_corner[0] - yPixel[4] = yp_corner[0] - - for i in range(0, 4): # 0:left, 1: right, 2: bottom, 3: top - nvertices2 = 0 - for j in range(0, nvertices): - x1 = xPixel[j] - y1 = yPixel[j] - x2 = xPixel[j + 1] - y2 = yPixel[j + 1] - condition = calcCondition(i, x1, y1, x2, y2, - left, right, top, bottom) - x = 0 - y = 0 - - if condition == 1: - x, y = solve_intersection(i, x1, y1, x2, y2, - left, right, top, bottom) - nvertices2 = addpoint(x, y, xnew, ynew, nvertices2) - nvertices2 = addpoint(x2, y2, xnew, ynew, nvertices2) - - elif condition == 2: - nvertices2 = addpoint(x2, y2, xnew, ynew, nvertices2) - elif condition == 3: - x, y = solve_intersection(i, x1, y1, x2, y2, - left, right, top, bottom) - nvertices2 = addpoint(x, y, xnew, ynew, nvertices2) - - # condition == 4: points outside - # Done looping over J corners - nvertices2 = addpoint(xnew[0], ynew[0], xnew, ynew, nvertices2) - - if nvertices2 > max_vertices: - raise Error2DPolygon(" Failure in finding the clipped polygon, nvertices2 > 9 ") - - nvertices = nvertices2 - 1 - - for k in range(0, nvertices2): - xPixel[k] = xnew[k] - yPixel[k] = ynew[k] - - # done loop over top,bottom,left,right - nvertices = nvertices + 1 - - if nvertices > 0: - area_clipped = find_area_poly(nvertices, xPixel, yPixel) - - return area_clipped -# _____________________________________________________________________________ - - -def match_det2cube(instrument, - x, y, sliceno, - input_model, - transform, - spaxel_flux, - spaxel_weight, - spaxel_iflux, - spaxel_var, - 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. - 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 - sliceno : int - slice number - input_model : datamodel - input slope 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. - - Returns - ------- - spaxel filled in with needed information on overlapping detector pixels - """ - - x = _toindex(x) - y = _toindex(y) - - pixel_dq = input_model.dq[y, x] - - 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 holds the location of pixels we want to map to cube - x = x[good_data] - y = y[good_data] - - coord1, coord2, lam = transform(x, y) - valid = ~np.isnan(coord2) - x = x[valid] - y = y[valid] - - yy_bot = y - yy_top = y + 1 - xx_left = x - xx_right = x + 1 - # along slice dimension is second coordinate returned from transform - 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) - b4, a4, lam4 = transform(xx_left, yy_top) - # check if units are in microns or meters, if meters convert to microns - # only need to check one of the wavelenghts - lmax = np.nanmax(lam1.flatten()) - if lmax < 0.0001: - lam1 = lam1 * 1.0e6 - lam2 = lam2 * 1.0e6 - lam3 = lam3 * 1.0e6 - lam4 = lam4 * 1.0e6 - - 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) - a4, b4, lam4 = transform(xx_left, yy_top) - - # corners are returned Nanned if outside range of slice - # fixed the nanned corners when a1,lam1 is valid but - # adding 1 to x,y pushes data outside BB valid region - - # index_bad2 = np.isnan(a2) - # index_bad3 = np.isnan(a3) - # index_bad4 = np.isnan(a4) - - # on the edge out of bounds - index_good2 = ~np.isnan(a2) - index_good3 = ~np.isnan(a3) - index_good4 = ~np.isnan(a4) - - # we need the cases of only all corners valid numbers - good = np.where(index_good2 & index_good3 & index_good4) - - a1 = a1[good] - a2 = a2[good] - a3 = a3[good] - a4 = a4[good] - lam1 = lam1[good] - lam2 = lam2[good] - lam3 = lam3[good] - lam4 = lam4[good] - x = x[good] - y = y[good] - # approximate what out of bound value should be - # w12 = np.mean(lam1[index_good2] - lam2[index_good2]) - # w13 = np.mean(lam1[index_good3] - lam3[index_good3]) - # w14 = np.mean(lam1[index_good4] - lam4[index_good4]) - - # a12 = np.mean(a1[index_good2] - a2[index_good2]) - # a13 = np.mean(a1[index_good3] - a3[index_good3]) - # a14 = np.mean(a1[index_good4] - a4[index_good4]) - - # a2[index_bad2] = a1[index_bad2] - a12 - # a3[index_bad3] = a1[index_bad3] - a13 - # a4[index_bad4] = a1[index_bad4] - a14 - - # lam2[index_bad2] = lam1[index_bad2] - w12 - # lam3[index_bad3] = lam1[index_bad3] - w13 - # lam4[index_bad4] = lam1[index_bad4] - w14 - - # center of first pixel, x,y = 1 for Adrian's equations - # but we want the pixel corners, x,y values passed into this - # routine to start at 0 - pixel_flux = input_model.data[y, x] - pixel_err = input_model.err[y, x] - - nzc = len(zcoord) - nac = len(acoord) - # 1-1 mapping in across slice direction (x for NIRSPEC, y for MIRI) - if instrument == 'NIRSPEC': - xx = sliceno - - if instrument == 'MIRI': - yy = sliceno - - # Loop over all pixels in slice - nn = len(x) - for ipixel in range(0, nn - 1): - # detector pixel -> 4 corners - # In along slice,wave space - # along slice: a - - along_corner = [] - wave_corner = [] - - along_corner.append(a1[ipixel]) - along_corner.append(a2[ipixel]) - along_corner.append(a3[ipixel]) - along_corner.append(a4[ipixel]) - - wave_corner.append(lam1[ipixel]) - wave_corner.append(lam2[ipixel]) - wave_corner.append(lam3[ipixel]) - wave_corner.append(lam4[ipixel]) -# ________________________________________________________________________________ -# Now it does not matter the WCS method used - along_min = min(along_corner) - along_max = max(along_corner) - wave_min = min(wave_corner) - wave_max = max(wave_corner) - - Area = find_area_quad(along_min, wave_min, along_corner, wave_corner) - - # estimate where the pixel overlaps in the cube - # find the min and max values in the cube xcoord,ycoord and zcoord - # along_min = -2.0 - # along_max = -1.86 - MinA = (along_min - crval_along) / cdelt_along - MaxA = (along_max - crval_along) / cdelt_along - ia1 = max(0, int(MinA)) - # ia2 = int(math.ceil(MaxA)) - ia2 = int(MaxA) - if ia2 >= nac: - ia2 = nac - 1 - - MinW = (wave_min - crval3) / cdelt3 - MaxW = (wave_max - crval3) / cdelt3 - iz1 = int(math.trunc(MinW)) - iz2 = int(math.ceil(MaxW)) - if iz2 >= nzc: - iz2 = nzc - 1 - - # loop over possible overlapping cube pixels - # noverlap = 0 - nplane = naxis1 * naxis2 - - for zz in range(iz1, iz2 + 1): - zcenter = zcoord[zz] - istart = zz * nplane - - for aa in range(ia1, ia2 + 1): - if instrument == 'NIRSPEC': - cube_index = istart + aa * naxis1 + xx # xx = slice # - - if instrument == 'MIRI': - cube_index = istart + yy * naxis1 + aa # xx = slice # - - acenter = acoord[aa] - area_overlap = sh_find_overlap(acenter, zcenter, - cdelt_along, cdelt3, - along_corner, wave_corner) - if area_overlap > 0.0: - AreaRatio = area_overlap / Area - - spaxel_flux[cube_index] = spaxel_flux[cube_index] + \ - (AreaRatio * pixel_flux[ipixel]) - spaxel_weight[cube_index] = spaxel_weight[cube_index] + \ - AreaRatio - spaxel_iflux[cube_index] = spaxel_iflux[cube_index] + 1 - spaxel_var[cube_index] = spaxel_var[cube_index] + \ - (AreaRatio * pixel_err[ipixel]) * (AreaRatio * pixel_err[ipixel]) -# ________________________________________________________________________________ - - -class Error2DPolygon(Exception): - """ Exception raise when no overlap between the detector pixel and - output plane is found - """ - pass diff --git a/jwst/cube_build/ifu_cube.py b/jwst/cube_build/ifu_cube.py index 5ec0dcdf35..c544ce1bc5 100644 --- a/jwst/cube_build/ifu_cube.py +++ b/jwst/cube_build/ifu_cube.py @@ -8,7 +8,6 @@ import logging import math from ..model_blender import blendmeta - from .. import datamodels from ..assign_wcs import pointing from jwst.transforms.models import _toindex @@ -19,10 +18,10 @@ from ..assign_wcs.util import wrap_ra from ..datamodels import dqflags from . import cube_build_wcs_util -from . import cube_overlap -from . import cube_cloud +from . import cube_internal_cal from . import coord from ..mrs_imatch.mrs_imatch_step import apply_background_2d +from .cube_match_sky import cube_wrapper # c extension log = logging.getLogger(__name__) log.setLevel(logging.DEBUG) @@ -136,9 +135,7 @@ def check_ifucube(self): IncorrectInput Interpolation = area was selected for when input data is more than one file or model - AreaInterpolation - If Interpolation = area then no user selected value can be set for - beta dimension of the output cube + """ num1 = len(self.list_par1) num_files = 0 @@ -148,8 +145,8 @@ def check_ifucube(self): n = len(self.master_table.FileMap[self.instrument][this_a][this_b]) num_files = num_files + n self.num_files = num_files -# do some basic checks on the cubes + # 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" + @@ -158,8 +155,8 @@ def check_ifucube(self): raise IncorrectInput("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 """ @@ -206,7 +203,7 @@ def define_cubename(self): if self.output_type == 'single': newname = self.output_name_base + ch_name + '-' + b_name + \ '_single_s3d.fits' -# ________________________________________________________________________________ + # ________________________________________________________________________________ elif self.instrument == 'NIRSPEC': # Check to see if the output base name already has a grating/prism @@ -227,11 +224,11 @@ def define_cubename(self): newname = self.output_name_base + fg_name + '_single_s3d.fits' if self.coord_system == 'internal_cal': newname = self.output_name_base + fg_name + '_internal_s3d.fits' -# ______________________________________________________________________________ + # ______________________________________________________________________________ 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 @@ -280,9 +277,9 @@ def set_geometry(self, corner_a, corner_b, lambda_min, lambda_max): 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 + # ________________________________________________________________________________ + # 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 systemtric cube centered on xi,eta = 0 xilimit = max(np.abs(xi_min), np.abs(xi_max)) etalimit = max(np.abs(eta_min), np.abs(eta_max)) @@ -306,42 +303,20 @@ def set_geometry(self, corner_a, corner_b, lambda_min, lambda_max): self.a_max = xi_max self.b_min = eta_min self.b_max = eta_max -# center of spaxels + # center of spaxels self.xcoord = np.zeros(self.naxis1) xstart = xi_min + self.cdelt1 / 2.0 - for i in range(self.naxis1): - self.xcoord[i] = xstart - xstart = xstart + 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 - for i in range(self.naxis2): - self.ycoord[i] = ystart - ystart = ystart + self.cdelt2 - - ygrid = np.zeros(self.naxis2 * self.naxis1) - xgrid = np.zeros(self.naxis2 * self.naxis1) - - k = 0 - ystart = self.ycoord[0] - for i in range(self.naxis2): - xstart = self.xcoord[0] - for j in range(self.naxis1): - xgrid[k] = xstart - ygrid[k] = ystart - xstart = xstart + self.cdelt1 - k = k + 1 - ystart = ystart + self.cdelt2 - -# ycube,xcube = np.mgrid[0:self.naxis2, -# 0:self.naxis1] -# xcube = xcube.flatten() -# ycube = ycube.flatten() + self.ycoord = np.arange(start=ystart, stop=ystart + self.naxis2 * self.cdelt2, step=self.cdelt2) - self.xcenters = xgrid - self.ycenters = ygrid - -# _______________________________________________________________________ + 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: @@ -360,9 +335,7 @@ def set_geometry(self, corner_a, corner_b, lambda_min, lambda_max): self.crval3 = self.lambda_min + self.cdelt3 / 2.0 self.crpix3 = 1.0 zstart = self.lambda_min + self.cdelt3 / 2.0 - for i in range(self.naxis3): - self.zcoord[i] = zstart - zstart = zstart + self.cdelt3 + self.zcoord = np.arange(start=zstart, stop=self.lambda_max, step=self.cdelt3) else: self.naxis3 = len(self.wavelength_table) @@ -376,7 +349,7 @@ 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 @@ -545,15 +518,12 @@ def build_ifucube(self): associated with the band 2. map_detector_to_output_frame: Maps the detector data to the cube output coordinate system 3. For each mapped detector pixel the ifu cube spaxel located in the region of - interest. There are three different routines to do this step each of them use - a slighly different weighting function in how to combine the detector fluxs that - fall within a region of influence from the spaxel center - a. cube_cloud:match_det2_cube_msm: This routine uses the modified + interest. There are two different routines to do this step, both of which use a c extension + to combine the detector fluxes that fall within a region of influence from the spaxel center + a. src/cube_match_sky: This routine uses the modified shepard method to determing the weighting function, which weights the detector fluxes based on the distance between the detector center and spaxel center. - b. cube_cloud:match_det2_cube_miripsf - the weighting function based on the width of the - psf and lsf. - c. cube_overlap.match_det2cube is only for single exposure, single band cubes and + b. src/cube_match_internalis only for single exposure, single band cubes and the ifucube in created in the detector plane. The weighting function is based on the overlap of between the detector pixel and spaxel. This method is simplified to determine the overlap in the along slice-wavelength plane. @@ -569,65 +539,33 @@ def build_ifucube(self): self.output_name = self.define_cubename() total_num = self.naxis1 * self.naxis2 * self.naxis3 - self.spaxel_flux = np.zeros(total_num) - self.spaxel_weight = np.zeros(total_num) - self.spaxel_iflux = np.zeros(total_num) - self.spaxel_var = np.zeros(total_num) - self.spaxel_dq = np.zeros((self.naxis3, self.naxis2 * self.naxis1), dtype=np.uint32) - - spaxel_ra = None - spaxel_dec = None - spaxel_wave = None -# ______________________________________________________________________________ -# Only performed if weighting = MIRIPSF, first convert xi,eta cube to -# v2,v3,wave. This information if past to cube_cloud and for each -# input_model the v2,v3, wave is converted to alpha,beta in detector plane -# ra,dec, wave is independent of input_model -# v2,v3, alpha,beta depends on the input_model - if self.weighting == 'miripsf': - spaxel_ra = np.zeros(total_num) - spaxel_dec = np.zeros(total_num) - spaxel_wave = np.zeros(total_num) - - nxy = self.xcenters.size - nz = self.zcoord.size - for iz in range(nz): - istart = iz * nxy - for ixy in range(nxy): - ii = istart + ixy - spaxel_ra[ii], spaxel_dec[ii] = coord.std2radec(self.crval1, - self.crval2, - self.xcenters[ixy], - self.ycenters[ixy]) - spaxel_wave[ii] = self.zcoord[iz] -# ______________________________________________________________________________ + + self.spaxel_flux = np.zeros(total_num,dtype=np.float64) + self.spaxel_weight = np.zeros(total_num, dtype=np.float64) + 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) + + # ______________________________________________________________________________ subtract_background = True - # now need to loop over every file that covers this - # channel/subchannel (MIRI) - # or Grating/filter(NIRSPEC) + # loop over every file that covers this channel/subchannel (MIRI) or + # Grating/filter(NIRSPEC) # and map the detector pixels to the cube spaxel - cube_debug = None - if self.xdebug is not None: - nplane = self.naxis1 * self.naxis2 - xydebug = self.ydebug * self.naxis1 + self.xdebug - cube_debug = (self.zdebug * nplane) + xydebug - log.info('Cube index debug %d', cube_debug) - log.info('%i %i %i %', xydebug, self.zdebug, self.ydebug, self.xdebug) - number_bands = len(self.list_par1) - for i in range(number_bands): - this_par1 = self.list_par1[i] - this_par2 = self.list_par2[i] + + for ib in range(number_bands): + this_par1 = self.list_par1[ib] + this_par2 = self.list_par2[ib] 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 + # ________________________________________________________________________________ + # 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] # set up input_model to be first file used to copy in basic header info # to ifucube meta data - if i == 0 and k == 0: + if ib == 0 and k == 0: input_model_ref = input_model log.debug(f"Working on Band defined by: {this_par1} {this_par2}") @@ -635,179 +573,141 @@ def build_ifucube(self): # POINTCLOUD used for skyalign and IFUalign # -------------------------------------------------------------------------------- if self.interpolation == 'pointcloud': - t0 = time.time() pixelresult = self.map_detector_to_outputframe(this_par1, subtract_background, input_model) coord1, coord2, wave, flux, err, slice_no, rois_pixel, roiw_pixel, weight_pixel,\ - softrad_pixel, scalerad_pixel, alpha_det, beta_det = pixelresult + softrad_pixel, scalerad_pixel = pixelresult + + # by default flag the dq plane based on the FOV of the detector projected to sky + flag_dq_plane = 1 + if self.skip_dqflagging: + flag_dq_plane = 0 + # check that there is valid data returned - # If all the data is flagged as DO_NOT_USE - not common then log warning and skip data + # If all the data is flagged as DO_NOT_USE - not common- then log warning if wave.size == 0: - no_data = True - else: - no_data = False + log.warning(f'No valid data found on file {input_model.meta.filename}') + flag_dq_plane = 0 - t1 = time.time() - log.debug("Time to transform pixels to output frame = %.1f s" % (t1 - t0,)) + t0 = time.time() + roiw_ave = np.mean(roiw_pixel) + # ______________________________________________________________________ + # C extension setup + # ______________________________________________________________________ + start_region = 0 + end_region = 0 - # If setting the DQ plane of the IFU - if self.skip_dqflagging or no_data: - log.info("Skipping setting DQ flagging") - else: - t0 = time.time() - roiw_ave = np.mean(roiw_pixel) - self.map_fov_to_dqplane(this_par1, coord1, coord2, wave, roiw_ave, slice_no) - t1 = time.time() - log.debug("Time to set initial dq values = %.1f s" % (t1 - t0,)) + if self.instrument == 'MIRI': + instrument = 0 + start_region = self.instrument_info.GetStartSlice(this_par1) + end_region = self.instrument_info.GetEndSlice(this_par1) + + else: # NIRSPEC + instrument = 1 + + result = None + weight_type = 0 # default to emsm + if self.weighting == 'msm': + weight_type = 1 + + 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(result[0], np.float64) + self.spaxel_weight = self.spaxel_weight + np.asarray(result[1], np.float64) + self.spaxel_var = self.spaxel_var + np.asarray(result[2], np.float64) + self.spaxel_iflux = self.spaxel_iflux + np.asarray(result[3],np.float64) + spaxel_dq.astype(np.uint) + self.spaxel_dq = np.bitwise_or(self.spaxel_dq, spaxel_dq) + result = None - if no_data: - log.warning(f'No valid data found on file {input_model.meta.filename}') - if self.weighting == 'msm' or self.weighting == 'emsm': - t0 = time.time() - cube_cloud.match_det2cube_msm(self.naxis1, self.naxis2, self.naxis3, - self.cdelt1, self.cdelt2, - self.cdelt3_normal, - self.xcenters, self.ycenters, self.zcoord, - self.spaxel_flux, - self.spaxel_weight, - self.spaxel_iflux, - self.spaxel_var, - flux, - err, - coord1, coord2, wave, - self.weighting, - rois_pixel, roiw_pixel, - weight_pixel, - softrad_pixel, - scalerad_pixel, - cube_debug, - self.debug_file) - - t1 = time.time() - log.debug("Time to match pixels to cube spaxels = %.1f s" % (t1 - t0,)) - # ________________________________________________________________________________ - elif self.weighting == 'miripsf': - wave_resol = self.instrument_info.Get_RP_ave_Wave(this_par1, - this_par2) - - alpha_resol = self.instrument_info.Get_psf_alpha_parameters() - beta_resol = self.instrument_info.Get_psf_beta_parameters() - wcsobj = input_model.meta.wcs - worldtov23 = wcsobj.get_transform(wcsobj.output_frame.name, "v2v3") - v2ab_transform = wcsobj.get_transform('v2v3', 'alpha_beta') - - spaxel_v2, spaxel_v3, zl = worldtov23(spaxel_ra, - spaxel_dec, - spaxel_wave) - - spaxel_alpha, spaxel_beta, spaxel_wave = v2ab_transform(spaxel_v2, - spaxel_v3, - zl) - cube_cloud.match_det2cube_miripsf(alpha_resol, - beta_resol, - wave_resol, - self.naxis1, self.naxis2, self.naxis3, - self.xcenters, self.ycenters, self.zcoord, - self.spaxel_flux, - self.spaxel_weight, - self.spaxel_iflux, - self.spaxel_var, - spaxel_alpha, spaxel_beta, spaxel_wave, - flux, - err, - coord1, coord2, wave, - alpha_det, beta_det, - 'emsm', - rois_pixel, roiw_pixel, - weight_pixel, - softrad_pixel, - scalerad_pixel) + # -------------------------------------------------------------------------------- + # # 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) + # MIRI # -------------------------------------------------------------------------------- - 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) - regions = list(range(start_region, end_region + 1)) - t0 = time.time() - - for i in regions: - log.info('Working on Slice # %d', i) - y, x = (det2ab_transform.label_mapper.mapper == i).nonzero() + 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) + regions = list(range(start_region, end_region + 1)) + t0 = time.time() + + for i in regions: + 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 - cube_overlap.match_det2cube(self.instrument, - x, y, slice, - input_model, - det2ab_transform, - self.spaxel_flux, - self.spaxel_weight, - self.spaxel_iflux, - self.spaxel_var, - self.xcoord, self.zcoord, - self.crval1, self.crval3, - self.cdelt1, self.cdelt3, - self.naxis1, self.naxis2) + index = np.where(y < 1023) + y = y[index] + x = x[index] + slice = i - start_region + cube_internal_cal.match_det2cube(self.instrument, + x, y, slice, + input_model, + det2ab_transform, + self.spaxel_flux, + self.spaxel_weight, + self.spaxel_iflux, + self.spaxel_var, + self.xcoord, self.zcoord, + self.crval1, self.crval3, + self.cdelt1, self.cdelt3, + self.naxis1, self.naxis2) t1 = time.time() - log.debug("Time to Map All slices on Detector to Cube = %.1f s" % (t1 - t0,)) - # -------------------------------------------------------------------------------- - # 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] - - for i in range(nslices): - # print('slice and slice map',i ,slicemap[i]) - slice_wcs = nirspec.nrs_wcs_set_input(input_model, i) - x, y = wcstools.grid_from_bounding_box(slice_wcs.bounding_box, step=(1, 1), center=True) - detector2slicer = slice_wcs.get_transform('detector', 'slicer') - - cube_overlap.match_det2cube(self.instrument, - x, y, slicemap[i], - input_model, - detector2slicer, - self.spaxel_flux, - self.spaxel_weight, - self.spaxel_iflux, - self.spaxel_var, - self.ycoord, self.zcoord, - self.crval2, self.crval3, - self.cdelt2, self.cdelt3, - self.naxis1, self.naxis2) - -# _______________________________________________________________________ -# Mapped all data to cube or Point Cloud -# now determine Cube Spaxel flux + # -------------------------------------------------------------------------------- + # 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] + + for i in range(nslices): + + slice_wcs = nirspec.nrs_wcs_set_input(input_model, i) + x, y = wcstools.grid_from_bounding_box(slice_wcs.bounding_box, step=(1, 1), center=True) + detector2slicer = slice_wcs.get_transform('detector', 'slicer') + + cube_internal_cal.match_det2cube(self.instrument, + x, y, slicemap[i], + input_model, + detector2slicer, + self.spaxel_flux, + self.spaxel_weight, + self.spaxel_iflux, + self.spaxel_var, + self.ycoord, self.zcoord, + self.crval2, self.crval3, + self.cdelt2, self.cdelt3, + self.naxis1, self.naxis2) + # _______________________________________________________________________ + # done looping over files self.find_spaxel_flux() self.set_final_dq_flags() - # result consist of ifu_cube and status + # shove Flux and iflux in the final IFU cube result = self.setup_final_ifucube_model(input_model_ref) -# _______________________________________________________________________ -# shove Flux and iflux in the final IFU cube - return result -# ******************************************************************************** + # ******************************************************************************** def build_ifucube_single(self): """ Build a set of single mode IFU cubes used for outlier detection and background matching @@ -820,26 +720,29 @@ def build_ifucube_single(self): # loop over input models single_ifucube_container = datamodels.ModelContainer() + weight_type = 0 # default to emsm + 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 j = 0 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 + # ________________________________________________________________________________ + # 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] - cube_debug = None + log.debug("Working on next Single IFU Cube = %i" % (j + 1)) t0 = time.time() # for each new data model create a new spaxel total_num = self.naxis1 * self.naxis2 * self.naxis3 - self.spaxel_flux = np.zeros(total_num) - self.spaxel_weight = np.zeros(total_num) + self.spaxel_flux = np.zeros(total_num,dtype=np.float64) + self.spaxel_weight = np.zeros(total_num, dtype=np.float64) self.spaxel_iflux = np.zeros(total_num) - self.spaxel_dq = np.zeros((self.naxis3, self.naxis2 * self.naxis1), dtype=np.uint32) - self.spaxel_var = np.zeros(total_num) + self.spaxel_dq = np.zeros(total_num, dtype=np.uint32) + self.spaxel_var = np.zeros(total_num,dtype=np.float64) subtract_background = False @@ -848,35 +751,41 @@ def build_ifucube_single(self): input_model) coord1, coord2, wave, flux, err, slice_no, rois_pixel, roiw_pixel, weight_pixel, \ - softrad_pixel, scalerad_pixel, alpha_det, beta_det = pixelresult - - cube_cloud.match_det2cube_msm(self.naxis1, - self.naxis2, - self.naxis3, - self.cdelt1, self.cdelt2, - self.cdelt3_normal, - self.xcenters, - self.ycenters, - self.zcoord, - self.spaxel_flux, - self.spaxel_weight, - self.spaxel_iflux, - self.spaxel_var, - flux, - err, - coord1, coord2, wave, - self.weighting, - rois_pixel, roiw_pixel, - weight_pixel, - softrad_pixel, - scalerad_pixel, - cube_debug, - self.debug_file) -# _______________________________________________________________________ -# shove Flux and iflux in the final ifucube + softrad_pixel, scalerad_pixel = pixelresult + + # 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': + instrument = 0 + else: # NIRSPEC + instrument = 1 + + result = None + 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(result[0], np.float64) + self.spaxel_weight = self.spaxel_weight + np.asarray(result[1], np.float64) + self.spaxel_var = self.spaxel_var + np.asarray(result[2], np.float64) + self.spaxel_iflux = self.spaxel_iflux + np.asarray(result[3],np.float64) + result = None + # ______________________________________________________________________ + # shove Flux and iflux in the final ifucube self.find_spaxel_flux() - # determine Cube Spaxel flux + # determine Cube Spaxel flux status = 0 result = self.setup_final_ifucube_model(input_model) ifucube_model, status = result @@ -888,8 +797,7 @@ 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 @@ -1068,15 +976,26 @@ def determine_cube_parameters(self): self.wavemax > table_wavelength[imax]): imax = imax + 1 - self.roiw_table = table_wroi[imin:imax + 1] - self.rois_table = table_sroi[imin:imax + 1] + num_table = imax - imin + 1 + self.scalerad_table = np.zeros(num_table) + self.scalerad_table[:] = table_scalerad[imin:imax + 1] + + self.softrad_table = np.zeros(num_table) + self.softrad_table[:] = table_softrad[imin:imax + 1] + + self.roiw_table = np.zeros(num_table) + self.roiw_table[:] = table_wroi[imin:imax + 1] + + self.rois_table = np.zeros(num_table) + self.rois_table[:] = table_sroi[imin:imax + 1] if self.num_files < 4: - self.rois_table = [i * 1.5 for i in self.rois_table] + self.rois_table = self.rois_table * 1.5 - self.softrad_table = table_softrad[imin:imax + 1] - self.weight_power_table = table_power[imin:imax + 1] - self.scalerad_table = table_scalerad[imin:imax + 1] - self.wavelength_table = table_wavelength[imin:imax + 1] + self.weight_power_table = np.zeros(num_table) + self.weight_power_table[:] = table_power[imin:imax + 1] + + self.wavelength_table = np.zeros(num_table) + self.wavelength_table[:] = table_wavelength[imin:imax + 1] # check if using default values from the table (not user set) if self.rois == 0.0: @@ -1278,8 +1197,10 @@ def setup_ifucube_wcs(self): # Find the footprint of the image spectral_found = hasattr(input_model.meta.wcsinfo,'spectral_region') spatial_found = hasattr(input_model.meta.wcsinfo,'s_region') - - if(spectral_found & spatial_found): + world = False + if self.coord_system == 'skyalign': + world = True + if(spectral_found & spatial_found & world): [lmin,lmax] = input_model.meta.wcsinfo.spectral_region spatial_box = input_model.meta.wcsinfo.s_region s = spatial_box.split(' ') @@ -1292,12 +1213,11 @@ def setup_ifucube_wcs(self): ca4 = float(s[9]) cb4 = float(s[10]) else: - log.info('Using old method of 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, - this_a, self.instrument_info, self.coord_system) ca1, cb1, ca2, cb2, ca3, cb3, ca4, cb4, lmin, lmax = ch_corners @@ -1323,8 +1243,8 @@ def setup_ifucube_wcs(self): lambda_min.append(lmin) lambda_max.append(lmax) - # ________________________________________________________________________________ - # done looping over files determine final size of cube + # ________________________________________________________________________________ + # done looping over files determine final size of cube corner_a = np.array(corner_a) corner_a = wrap_ra(corner_a) @@ -1367,16 +1287,16 @@ def setup_ifucube_wcs(self): 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 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}') -# ________________________________________________________________________________ - # Based on Scaling and Min and Max values determine naxis1, naxis2, naxis3 - # set cube CRVALs, CRPIXs + # ________________________________________________________________________________ + # 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': self.set_geometry(corner_a, corner_b, final_lambda_min, final_lambda_max) @@ -1386,7 +1306,6 @@ def setup_ifucube_wcs(self): self.print_cube_geometry() # ************************************************************************** - def map_detector_to_outputframe(self, this_par1, subtract_background, input_model): @@ -1433,30 +1352,20 @@ def map_detector_to_outputframe(self, this_par1, weighting parameter assocation with coord1,coord2 softrad_det : numpy.ndarray weighting parameter assocation with coord1,coord2 - alpha_det : numpy.ndarray - alpha coordinate of pixels - beta_det : numpy.ndarray - beta coordinate of pixel """ # intitalize alpha_det and beta_det to None. These are filled in # if the instrument is MIRI and the weighting is miripsf - alpha_det = None - beta_det = None - coord1 = None - coord2 = None - flux = None - err = None - wave = None - slice_no = None # Slice number + wave_all = None + slice_no_all = None # Slice number if self.instrument == 'MIRI': sky_result = self.map_miri_pixel_to_sky(input_model, this_par1, subtract_background) - (x, y, ra, dec, wave, slice_no, alpha, beta) = sky_result + (x, y, ra, dec, wave_all, slice_no_all) = sky_result elif self.instrument == 'NIRSPEC': sky_result = self.map_nirspec_pixel_to_sky(input_model) - (x, y, ra, dec, wave, slice_no) = sky_result + (x, y, ra, dec, wave_all, slice_no_all) = sky_result # ______________________________________________________________________________ # The following is for both MIRI and NIRSPEC @@ -1471,20 +1380,18 @@ def map_detector_to_outputframe(self, this_par1, # of interest would have them fall within one of the cube spectral planes # Note that the cube lambda refer to the spaxel midpoints, so we must account for both # the spaxel width and the ROI size + if self.linear_wavelength: - min_wave_tolerance = self.crval3 - np.absolute(self.zcoord[1] - self.zcoord[0]) / 2 - self.roiw - max_wave_tolerance = (self.zcoord[-1] + np.absolute(self.zcoord[-1] - self.zcoord[-2]) / 2 - + self.roiw) + min_wave_tolerance = self.zcoord[0] - self.roiw + max_wave_tolerance = self.zcoord[-1] + self.roiw else: - min_wave_tolerance = self.crval3 - np.absolute(self.zcoord[1] - self.zcoord[0]) / 2 - np.max(self.roiw_table) - max_wave_tolerance = (self.zcoord[-1] + np.absolute(self.zcoord[-1] - self.zcoord[-2]) / 2 - + np.max(self.roiw_table)) + min_wave_tolerance = self.zcoord[0] - np.max(self.roiw_table) + max_wave_tolerance = self.zcoord[-1] + np.max(self.roiw_table) - valid_min = np.where(wave >= min_wave_tolerance) - not_mapped_low = wave.size - len(valid_min[0]) - - valid_max = np.where(wave <= max_wave_tolerance) - not_mapped_high = wave.size - len(valid_max[0]) + valid_min = np.where(wave_all >= min_wave_tolerance) + not_mapped_low = wave_all.size - len(valid_min[0]) + 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}') @@ -1496,21 +1403,37 @@ def map_detector_to_outputframe(self, this_par1, # ______________________________________________________________________________ # using the DQFlags from the input_image find pixels that should be excluded # from the cube mapping - all_flags = (dqflags.pixel['DO_NOT_USE'] + - dqflags.pixel['NON_SCIENCE']) + # all_flags = (dqflags.pixel['DO_NOT_USE'] + + # dqflags.pixel['NON_SCIENCE']) - valid3 = np.logical_and((wave >= min_wave_tolerance), - (wave <= max_wave_tolerance)) + valid3 = np.logical_and(wave_all >= min_wave_tolerance, + wave_all <= max_wave_tolerance) # find the location of good data - good_data = np.where((np.bitwise_and(dq_all, all_flags) == 0) & - (valid2) & (valid3)) + # good_data = np.where((np.bitwise_and(dq_all, all_flags) == 0) & + # (valid2) & (valid3)) + + 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) # good data holds the location of pixels we want to map to cube - flux = flux_all[good_data] - err = err_all[good_data] - wave = wave[good_data] - slice_no = slice_no[good_data] + # define variables as numpy arrays (numba needs this defined) + flux_all_good = flux_all[good_data] + good_shape = flux_all_good.shape + flux = np.zeros(good_shape, dtype=np.float64) + err = np.zeros(good_shape, dtype=np.float64) + coord1 = np.zeros(good_shape, dtype=np.float64) + coord2 = np.zeros(good_shape, dtype=np.float64) + wave = np.zeros(good_shape, dtype=np.float64) + slice_no = np.zeros(good_shape) + + flux[:] = flux_all_good + err[:] = err_all[good_data] + wave[:] = wave_all[good_data] + slice_no[:] = slice_no_all[good_data] + + log.info(f'After removing pixels based on criteria 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 @@ -1528,14 +1451,20 @@ 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): - ifound = (np.abs(self.wavelength_table - w)).argmin() - rois_det[iw] = self.rois_table[ifound] - roiw_det[iw] = self.roiw_table[ifound] - softrad_det[iw] = self.softrad_table[ifound] - weight_det[iw] = self.weight_power_table[ifound] - scalerad_det[iw] = self.scalerad_table[ifound] + 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) ra_use = ra[good_data] dec_use = dec[good_data] coord1, coord2 = coord.radec2std(self.crval1, @@ -1543,12 +1472,8 @@ def map_detector_to_outputframe(self, this_par1, ra_use, dec_use, self.rot_angle) - if self.weighting == 'miripsf': - alpha_det = alpha[good_data] - beta_det = beta[good_data] - return coord1, coord2, wave, flux, err, slice_no, rois_det, roiw_det, weight_det, \ - softrad_det, scalerad_det, alpha_det, beta_det + softrad_det, scalerad_det # ______________________________________________________________________ def map_miri_pixel_to_sky(self, input_model, this_par1, subtract_background): @@ -1571,11 +1496,8 @@ def map_miri_pixel_to_sky(self, input_model, this_par1, subtract_background): Returns ------- x, y, ra, dec, lambda, slice_no of valid slice pixels - if weighting = mirispf returns alpha, beta """ - alpha = None - beta = None wave = None slice_no = None # Slice number @@ -1627,16 +1549,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.weighting == 'miripsf': - alpha, beta, lam = det2ab_transform(x, y) - - valid1 = ~np.isnan(alpha) - alpha = alpha[valid1] - beta = beta[valid1] - wave = wave[valid1] - x = x[valid1] - y = y[valid1] - sky_result = (x, y, ra, dec, wave, slice_no, alpha, beta) + sky_result = (x, y, ra, dec, wave, slice_no) return sky_result # ______________________________________________________________________ @@ -1645,13 +1558,12 @@ def map_nirspec_pixel_to_sky(self, 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. Parameter ---------- input: datamodel - input data model + input data model Returns ------- @@ -1674,9 +1586,12 @@ def map_nirspec_pixel_to_sky(self, input_model): nslices = 30 log.info("Mapping each NIRSpec slice to sky for input file") - t80 = time.time() + # t80 = time.time() for ii in range(nslices): + # t10 = time.time() slice_wcs = nirspec.nrs_wcs_set_input(input_model, ii) + # t11 = time.time() + # log.debug(f'Time to run nirspec nrs_wcs_set_input {t11-t10} ') x, y = wcstools.grid_from_bounding_box(slice_wcs.bounding_box) ra, dec, lam = slice_wcs(x, y) @@ -1700,10 +1615,9 @@ def map_nirspec_pixel_to_sky(self, input_model): lam_det[yind, xind] = lam flag_det[yind, xind] = 1 slice_det[yind, xind] = ii + 1 - # after looping over slices - pull out valid values - t81 = time.time() - log.debug(f'Time to map slices {t81-t80} ') + # t81 = time.time() + # log.debug(f'Time to map NIRSpec slices to sky {t81-t80} ') valid_data = np.where(flag_det == 1) y, x = valid_data @@ -1715,361 +1629,37 @@ def map_nirspec_pixel_to_sky(self, input_model): sky_result = (x, y, ra, dec, wave, slice_no) return sky_result - # ______________________________________________________________________ - def map_fov_to_dqplane(self, this_par1, coord1, coord2, wave, roiw_ave, slice_no): - """ Set an initial DQ flag for the IFU cube based on FOV of input data - - Map the FOV of each exposure channel (MIRI) or slice (NIRSPEC) to the DQ plane - and set an initial DQ flagging. The process is different for MIRI and NIRSpec. - In the MIRI case all the slices map roughly to the same FOV across the - wavelength range covered by the IFU. This is not the case for NIRSpec the - 30 different slices map to different FOV on the range of wavelengths. - - Paramteter - --------- - this_par1: Channel (MIRI) or Grating (NIRSpec) - 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 - slice_no: integer slice value of input data (used in MIRI case to find - the points of the edge slices.) - """ - - # MIRI mapping: - # The FOV is roughly the same for all the wavelength ranges. - # The offset in the slices makes the calculation of the four corners - # of the FOV more complicated. So we only use the two slices at - # the edges of the FOV to define the 4 corners. - # Note we can not use the wcs.footprint because this footprint only - # consists of 4 values ra min, ra max, dec min, dec max and we - # need 4 corners made of 8 different values. - - if self.instrument == 'MIRI': - - # find the wavelength boundaries of the band - use two extreme slices - wavemin = np.amin(wave) - wavemax = np.amax(wave) - - # self.zcoord holds the center of the wavelength bin - iwavemin = np.absolute(np.array(wavemin - self.zcoord) / self.cdelt3_normal) - iwavemax = np.absolute(np.array(wavemax - self.zcoord) / self.cdelt3_normal) - - imin = np.where(iwavemin == np.amin(iwavemin))[0] - imax = np.where(iwavemax == np.amin(iwavemax))[0] - - # for each wavelength plane - find the 2 extreme slices to set the FOV - for w in range(imin[0], imax[0]): - wave_distance = np.absolute(self.zcoord[w] - wave) - - # pull out the two extreme slices in the FOV. - # use these points to set the FOV - start_region = self.instrument_info.GetStartSlice(this_par1) - end_region = self.instrument_info.GetEndSlice(this_par1) - - istart = start_region - coord1_start = None - index_use = np.where((wave_distance < roiw_ave) & (slice_no == istart)) - if len(index_use[0]) > 0: - coord2_start = coord2[index_use] - coord1_start = coord1[index_use] - - iend = end_region - coord1_end = None - index_use = np.where((wave_distance < roiw_ave) & (slice_no == iend)) - if len(index_use[0]) > 0: - coord2_end = coord2[index_use] - coord1_end = coord1[index_use] - - # if there is valid data on this wavelength plane (not in a gap) - if coord1_start is not None and coord1_end is not None: - coord1_total = np.concatenate((coord1_start, coord1_end), axis=0) - coord2_total = np.concatenate((coord2_start, coord2_end), axis=0) - - # from an array of x and y values (contained in coord1_total and coord2_total) - # determine the footprint - footprint_all = self.four_corners(coord1_total, coord2_total) - - isline, footprint = footprint_all - (xi1, eta1, xi2, eta2, xi3, eta3, xi4, eta4) = footprint - - # find the overlap of FOV footprint and with IFU Cube - xi_corner = np.array([xi1, xi2, xi3, xi4]) - eta_corner = np.array([eta1, eta2, eta3, eta4]) - self.overlap_fov_with_spaxels(xi_corner, eta_corner, w, w) - - # NIRSpec Mapping: - # The FOV of each NIRSpec slice varies across the wavelength range. - # Each slice is mapped to each IFU wavelength plane - # The FOV of the slice is really just a line, so instead of using - # the routines the finds the overlap between a polygon and regular grid- - # which is used for MIRI - an algorithm that determines the spaxels that - # the slice line intersects is used instead. - - elif self.instrument == 'NIRSPEC': - # for each of the 30 slices - find the projection of this slice - # onto each of the IFU wavelength planes. - for islice in range(30): - index_slice = np.where(slice_no == islice + 1) - - # find the smaller set of wavelengths to search over for this slice - wavemin = np.amin(wave[index_slice]) - wavemax = np.amax(wave[index_slice]) - iwavemin = np.absolute(np.array(wavemin - self.zcoord) / (self.cdelt3_normal)) - iwavemax = np.absolute(np.array(wavemax - self.zcoord) / (self.cdelt3_normal)) - imin = np.where(iwavemin == np.amin(iwavemin))[0] - imax = np.where(iwavemax == np.amin(iwavemax))[0] - - # loop over valid wavelengths for slice and find the projection of the - # slice on the wavelength plane - - for w in range(imin[0], imax[0]): - wave_distance = np.absolute(self.zcoord[w] - wave) - index_use = np.where((wave_distance < roiw_ave) & (slice_no == islice + 1)) - # using only coord values for wavelength slice - if len(index_use[0]) > 0: - coord2_use = coord2[index_use] - coord1_use = coord1[index_use] - footprint_all = self.four_corners(coord1_use, coord2_use) - isline, footprint = footprint_all - # isline is true if slice values fall on line (which the will for NIRSpec) - # TODO for setting the dqfalg - should we be looking at footprint + sroi ? - (xi1, eta1, xi2, eta2, xi3, eta3, xi4, eta4) = footprint - # find the overlap with IFU Cube - xi_corner = np.array([xi1, xi2, xi3, xi4]) - eta_corner = np.array([eta1, eta2, eta3, eta4]) - - if isline: - self.overlap_slice_with_spaxels(xi_corner, eta_corner, w) - else: - self.overlap_fov_with_spaxels(xi_corner, eta_corner, w, w) -# ******************************************************************************** - - def overlap_slice_with_spaxels(self, xi_corner, eta_corner, w): - """ Set the initial dq plane of 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 - which IFU spaxels the slice falls on to set an initial dq flag. - - Parameters - --------- - xi_corner: holds the x starting and ending points of the slice - eta_corner: holds the y starting and ending points of the slice - wavelength: the wavelength bin of the IFU cube working with - - Sets - ---- - self.spaxel_dq : numpy.ndarray containing intermediate dq flag - - """ - - points = self.findpoints_on_slice(xi_corner, eta_corner) - num = len(points) - - for i in range(num): - xpt, ypt = points[i] - index = (ypt * self.naxis1) + xpt - self.spaxel_dq[w, index] = self.overlap_partial - -# ******************************************************************************** - - def findpoints_on_slice(self, xi_corner, eta_corner): - """ Bresenham's Line Algorithm to find points a line intersects with grid. - - Given the endpoints of a line find the spaxels this line intersects. - - Parameters - ----------- - xi_corner: holds the started in ending x values - eta_corner: holds the started in ending y values - - - Returns - ------- - Points: a tuple of x,y spaxel values that this line intersects - - """ - - # set up line - convert to integer values - x1 = int((xi_corner[0] - self.xcoord[0]) / self.cdelt1) - y1 = int((eta_corner[0] - self.ycoord[0]) / self.cdelt2) - x2 = int((xi_corner[1] - self.xcoord[0]) / self.cdelt1) - y2 = int((eta_corner[1] - self.ycoord[0]) / self.cdelt2) - - dx = x2 - x1 - dy = y2 - y1 - - # how steep is it - is_steep = abs(dy) > abs(dx) - - # Rotate line - if is_steep: - x1, y1 = y1, x1 - x2, y2 = y2, x2 - - # Swap start and end points if necessary and store swap state - swapped = False - if x1 > x2: - x1, x2 = x2, x1 - y1, y2 = y2, y1 - swapped = True - - # Recalculate differences - dx = x2 - x1 - dy = y2 - y1 - - # calculate error - error = int(dx / 2.0) - ystep = -1 - if y1 < y2: - ystep = 1 - - # iterate over grid to generate points between the start and end of line - y = y1 - points = [] - for x in range(x1, x2 + 1): - coord = (y, x) if is_steep else (x, y) - points.append(coord) - error -= abs(dy) - if error < 0: - y += ystep - error += dx - - # If coords were swapped then reverse - if swapped: - points.reverse() - return points -# ******************************************************************************** - - def overlap_fov_with_spaxels(self, xi_corner, eta_corner, wmin, wmax): - """find the amount of overlap of FOV with each spaxel - - Given the corners of the FOV find the spaxels that - 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. self.overlap_partial = overlap partial - b self.overlap_full = overlap_full - bit_wise combination of these values is allowed to account for - dithered FOVs. - - Parameter - ---------- - xi_corner: xi coordinates of the 4 corners of the FOV on the wavelenghth plane - eta_corner: eta coordinates of the 4 corners of the FOV on the wavelength plane - wmin: minimum wavelength bin in the IFU cube that this data covers - wmax: maximum wavelength bin in the IFU cube that this data covers - - Sets - ------- - self.spaxel_dq : numpy.ndarray containing intermediate dq flag + # ******************************************************************************** + 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 """ - - # loop over spaxels in the wavelength plane and set slice_dq - # lets roughly find the spaxels that might be overlapped - - ximin = np.amin(xi_corner) - ximax = np.amax(xi_corner) - etamin = np.amin(eta_corner) - etamax = np.amax(eta_corner) - index = np.where(((self.xcenters - self.cdelt1 / 2) > ximin) & - ((self.xcenters + self.cdelt1 / 2) < ximax) & - ((self.ycenters - self.cdelt2 / 2) > etamin) & - ((self.ycenters + self.cdelt2 / 2) < etamax)) - - wave_slice_dq = np.zeros(self.naxis2 * self.naxis1, dtype=np.int32) - area_box = self.cdelt1 * self.cdelt2 - for ixy in range(len(index[0])): - area_overlap = cube_overlap.sh_find_overlap(self.xcenters[index][ixy], - self.ycenters[index][ixy], - self.cdelt1, self.cdelt2, - xi_corner, eta_corner) - - overlap_coverage = area_overlap / area_box - - if overlap_coverage > self.tolerance_dq_overlap: - if overlap_coverage > 0.95: - wave_slice_dq[ixy] = self.overlap_full - else: - wave_slice_dq[ixy] = self.overlap_partial - - # set for a range of wavelengths - if wmin != wmax: - self.spaxel_dq[wmin:wmax, :] = np.bitwise_or(self.spaxel_dq[wmin:wmax, :], - wave_slice_dq) - - # set for a single wavelength - else: - self.spaxel_dq[wmin, :] = np.bitwise_or(self.spaxel_dq[wmin, :], - wave_slice_dq) -# ******************************************************************************* - - def four_corners(self, coord1, coord2): - """ helper function to compute the four corners of the FOV - - From an array of x and y values find the 4 corners enclosing these points - This routine defines the four corners as - corner 1: location of min coord2 - corner 2: location of min coord1 - corner 3: location of max coord2 - corner 4: location of max coord1 - - Parameter - ---------- - coord1: array of 4 x corners - coord2: array of 4 y corners - - Returns - ------- - Footprint of 4 corners - Is the data contained in coor1 and coord2 represented by a line if yes: - isline = True if not isline = False - - """ - - isline = False - index = np.where(coord2 == np.amin(coord2)) - xi_corner1 = coord1[index[0]] - eta_corner1 = coord2[index[0]] - - index = np.where(coord1 == np.amax(coord1)) - xi_corner2 = coord1[index[0]] - eta_corner2 = coord2[index[0]] - - index = np.where(coord2 == np.amax(coord2)) - xi_corner3 = coord1[index[0]] - eta_corner3 = coord2[index[0]] - - index = np.where(coord1 == np.amin(coord1)) - xi_corner4 = coord1[index[0]] - eta_corner4 = coord2[index[0]] - footprint = (xi_corner1[0], eta_corner1[0], - xi_corner2[0], eta_corner2[0], - xi_corner3[0], eta_corner3[0], - xi_corner4[0], eta_corner4[0]) - - distance_min_points = math.sqrt((xi_corner1 - xi_corner4)**2 + - (eta_corner1 - eta_corner4)**2) - - distance_max_points = math.sqrt((xi_corner2 - xi_corner3)**2 + - (eta_corner2 - eta_corner3)**2) - dist_tolerance = 0.0001 # tolerance used if points fall on a line - if ((distance_min_points < dist_tolerance) and - (distance_max_points < dist_tolerance)): - isline = True - footprint_all = (isline, footprint) - - return footprint_all -# ******************************************************************************** - + ifound = (np.abs(wavelength_table - w)).argmin() + rois_det[iw] = rois_table[ifound] + roiw_det[iw] = roiw_table[ifound] + softrad_det[iw] = softrad_table[ifound] + 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 """ -# 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. + # 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': good = self.spaxel_iflux > 0 self.spaxel_flux[good] = self.spaxel_flux[good] / self.spaxel_weight[good] @@ -2077,13 +1667,15 @@ def find_spaxel_flux(self): elif self.interpolation == 'pointcloud': # 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]) -# ******************************************************************************** + # ******************************************************************************** def set_final_dq_flags(self): + """ Set up the final dq flags, Good data(0) , NON_SCIENCE or DO_NOT_USE """ @@ -2185,28 +1777,28 @@ def set_final_dq_flags(self): 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 return IFUCube model """ - status = 0 # loop over the wavelength planes to confirm each plane has some data # for initial or final planes that do not have any data - eliminated them # from the IFUcube # Rearrange values from 1d vectors into 3d cubes - temp_flux = self.spaxel_flux.reshape((self.naxis3, - self.naxis2, self.naxis1)) - temp_wmap = self.spaxel_iflux.reshape((self.naxis3, - self.naxis2, self.naxis1)) - temp_dq = self.spaxel_dq.reshape((self.naxis3, + + flux = self.spaxel_flux.reshape((self.naxis3, + self.naxis2, self.naxis1)) + wmap = self.spaxel_iflux.reshape((self.naxis3, self.naxis2, self.naxis1)) - temp_var = self.spaxel_var.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)) # clean up empty wavelength planes except for single case if self.output_type != 'single': @@ -2214,7 +1806,7 @@ def setup_final_ifucube_model(self, model_ref): k = 0 found = 0 while (k < self.naxis3 and found == 0): - flux_at_wave = temp_flux[k, :, :] + flux_at_wave = flux[k, :, :] sum = np.nansum(flux_at_wave) if sum == 0.0: remove_start = remove_start + 1 @@ -2227,7 +1819,7 @@ def setup_final_ifucube_model(self, model_ref): found = 0 k = self.naxis3 - 1 while (k > 0 and found == 0): - flux_at_wave = temp_flux[k, :, :] + flux_at_wave = flux[k, :, :] sum = np.nansum(flux_at_wave) if sum == 0.0: remove_final = remove_final + 1 @@ -2247,10 +1839,10 @@ def setup_final_ifucube_model(self, model_ref): log.info('Number of wavelength planes removed with no data: %i', remove_total) - temp_flux = temp_flux[remove_start:self.naxis3 - remove_final, :, :] - temp_wmap = temp_wmap[remove_start:self.naxis3 - remove_final, :, :] - temp_dq = temp_dq[remove_start:self.naxis3 - remove_final, :, :] - temp_var = temp_var[remove_start:self.naxis3 - remove_final, :, :] + flux = flux[remove_start:self.naxis3 - remove_final, :, :] + wmap = wmap[remove_start:self.naxis3 - remove_final, :, :] + dq = dq[remove_start:self.naxis3 - remove_final, :, :] + var = var[remove_start:self.naxis3 - remove_final, :, :] if self.linear_wavelength: self.crval3 = self.zcoord[remove_start] @@ -2260,19 +1852,12 @@ def setup_final_ifucube_model(self, model_ref): self.naxis3 = self.naxis3 - (remove_start + remove_final) # end removing empty wavelength planes - naxis1 = self.naxis1 - naxis2 = self.naxis2 - naxis3 = self.naxis3 - - data = np.zeros((naxis3, naxis2, naxis1)) - idata = np.zeros((naxis3, naxis2, naxis1)) - dq_cube = np.zeros((naxis3, naxis2, naxis1)) - err_cube = np.zeros((naxis3, naxis2, naxis1)) + var = np.sqrt(var) if self.linear_wavelength: - ifucube_model = datamodels.IFUCubeModel(data=data, dq=dq_cube, - err=err_cube, - weightmap=idata) + 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) @@ -2281,19 +1866,14 @@ def setup_final_ifucube_model(self, model_ref): dtype=[('wavelength', ' 1: saved_model_type = ifucube_model.meta.model_type @@ -2402,125 +1982,8 @@ def setup_final_ifucube_model(self, model_ref): if self.instrument == 'MIRI': ifucube_model.meta.wcsinfo.cunit1 = 'arcsec' ifucube_model.meta.wcsinfo.cunit2 = 'arcsec' -# we only need to check list_par1[0] and list_par2[0] because these types -# of cubes are made from 1 exposures (setup_cube checks this at the start -# of cube_build). - if self.list_par1[0] == '1' and self.list_par2[0] == 'short': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL1A' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE1A' - if self.list_par1[0] == '2' and self.list_par2[0] == 'short': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL2A' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE2A' - if self.list_par1[0] == '3' and self.list_par2[0] == 'short': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL3A' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE3A' - if self.list_par1[0] == '4' and self.list_par2[0] == 'short': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL4A' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE4A' - - if self.list_par1[0] == '1' and self.list_par2[0] == 'medium': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL1B' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE1B' - if self.list_par1[0] == '2' and self.list_par2[0] == 'medium': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL2B' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE2B' - if self.list_par1[0] == '3' and self.list_par2[0] == 'medium': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL3B' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE3B' - if self.list_par1[0] == '4' and self.list_par2[0] == 'medium': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL4B' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE4B' - - if self.list_par1[0] == '1' and self.list_par2[0] == 'long': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL1C' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE1C' - if self.list_par1[0] == '2' and self.list_par2[0] == 'long': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL2C' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE2C' - if self.list_par1[0] == '3' and self.list_par2[0] == 'long': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL3C' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE3C' - if self.list_par1[0] == '4' and self.list_par2[0] == 'long': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL4C' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE4C' - - if self.list_par1[0] == '1' and self.list_par2[0] == 'short-medium': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL1A' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE1A' - if self.list_par1[0] == '2' and self.list_par2[0] == 'short-medium': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL2B' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE2B' - if self.list_par1[0] == '3' and self.list_par2[0] == 'short-medium': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL3B' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE3B' - if self.list_par1[0] == '4' and self.list_par2[0] == 'short-medium': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL4A' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE4A' - - if self.list_par1[0] == '1' and self.list_par2[0] == 'short-long': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL1A' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE1A' - if self.list_par1[0] == '2' and self.list_par2[0] == 'short-long': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL2C' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE2C' - if self.list_par1[0] == '3' and self.list_par2[0] == 'short-long': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL3C' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE3C' - if self.list_par1[0] == '4' and self.list_par2[0] == 'short-long': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL4A' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE4A' - - if self.list_par1[0] == '1' and self.list_par2[0] == 'medium-short': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL1B' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE1B' - if self.list_par1[0] == '2' and self.list_par2[0] == 'medium-short': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL2A' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE2A' - if self.list_par1[0] == '3' and self.list_par2[0] == 'medium-short': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL3A' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE3A' - if self.list_par1[0] == '4' and self.list_par2[0] == 'medium-short': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL4B' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE4B' - - if self.list_par1[0] == '1' and self.list_par2[0] == 'medium-long': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL1B' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE1B' - if self.list_par1[0] == '2' and self.list_par2[0] == 'medium-long': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL2C' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE2C' - if self.list_par1[0] == '3' and self.list_par2[0] == 'medium-long': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL3C' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE3C' - if self.list_par1[0] == '4' and self.list_par2[0] == 'medium-long': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL4B' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE4B' - - if self.list_par1[0] == '1' and self.list_par2[0] == 'long-short': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL1C' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE1C' - if self.list_par1[0] == '2' and self.list_par2[0] == 'long-short': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL2A' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE2A' - if self.list_par1[0] == '3' and self.list_par2[0] == 'long-short': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL3A' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE3A' - if self.list_par1[0] == '4' and self.list_par2[0] == 'long-short': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL4C' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE4C' - - if self.list_par1[0] == '1' and self.list_par2[0] == 'long-medium': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL1C' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE1C' - if self.list_par1[0] == '2' and self.list_par2[0] == 'long-medium': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL2B' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE2B' - if self.list_par1[0] == '3' and self.list_par2[0] == 'long-medium': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL3B' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE3B' - if self.list_par1[0] == '4' and self.list_par2[0] == 'long-medium': - ifucube_model.meta.wcsinfo.ctype1 = 'MRSAL4C' - ifucube_model.meta.wcsinfo.ctype2 = 'MRSBE4C' + ifucube_model.meta.wcsinfo.ctype1 = 'MRSALPHA' + ifucube_model.meta.wcsinfo.ctype2 = 'MRSBETA' if self.instrument == 'NIRSPEC': ifucube_model.meta.wcsinfo.cunit1 = 'meter' @@ -2528,12 +1991,12 @@ def setup_final_ifucube_model(self, model_ref): ifucube_model.meta.wcsinfo.ctype1 = 'NRSLICERX' ifucube_model.meta.wcsinfo.ctype2 = 'NRSSLICERY' -# set WCS information + # set WCS information wcsobj = pointing.create_fitswcs(ifucube_model) ifucube_model.meta.wcs = wcsobj - ifucube_model.meta.wcs.bounding_box = ((0, naxis1 - 1), - (0, naxis2 - 1), - (0, naxis3 - 1)) + ifucube_model.meta.wcs.bounding_box = ((0, self.naxis1 - 1), + (0, self.naxis2 - 1), + (0, self.naxis3 - 1)) ifucube_model.meta.cal_step.cube_build = 'COMPLETE' # problem with cube_build - contains only 0 data @@ -2542,8 +2005,8 @@ def setup_final_ifucube_model(self, model_ref): result = (ifucube_model, status) return result -# ******************************************************************************** + # ******************************************************************************** def blend_output_metadata(self, IFUCube): """Create new output metadata based on blending all input metadata.""" # Run fitsblender on output product @@ -2563,11 +2026,3 @@ class IncorrectParameter(Exception): """ Raises an exception if cube building parameter is nan """ pass - - -class AreaInterpolation(Exception): - """ Raises an exception if input parameter, Interpolation, is set to area - and the input parameter of the spatial scale of second dimension of cube - is also set. - """ - pass diff --git a/jwst/cube_build/src/cube_match_internal.c b/jwst/cube_build/src/cube_match_internal.c new file mode 100644 index 0000000000..887a1c9de3 --- /dev/null +++ b/jwst/cube_build/src/cube_match_internal.c @@ -0,0 +1,433 @@ +/* +This module forms the IFU cube in the "detector 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 +breaks the 3-D representation of the pixel's coordinates in the output plane to 2D. Standard +drizzling type routines are used to calculate the percentage of each pixel's overlap with +the output spaxels. The percentage of the pixel area overlapped onto an output +spaxel is used to weight the pixel flux for determining the spaxel flux. + + +Main function for Python: cube_wrapper_internal + +Python signature: + 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) +provide more details + +The output of this function is a tuple of 5 arrays:(spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux, spaxel_dq) +example output + +Parameters +---------- +instrument: int + 0 = MIRI, 1 = NIRSPEC. Used for set the dq plane +naxis1 : int + axis 1 of IFU cube +nasix2 : int + axis 2 of IFU cube +crval_along : double + along slice reference value in IFU cube +cdelt_along : double + along slice sampling in IFU cube +crval3 : double + wavelength reference value in IFU cube +cdetl3 : double + wavelength sampling in IFU cube +a1, a2, a3, a4: double array + 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 + 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 +zcoord : double array + Array holds the wavelength coordinates in the IFU cube +ss : double array + Array holds the slice number +pixel_flux : double array + Array that holds detector pixel flux +pixel_error : double array + Array that holds detector pixel flux error + +Returns +------- + +spaxel_flux : array +spaxel_weight : array +spaxel_iflux : array +spaxel_var : array +spaxel_dq : array + +*/ + +#include +#include +#include +#include +#include +#include +#include + +#define PY_ARRAY_UNIQUE_SYMBOL _jwst_cube_match_sky_numpy_api //WHAT IS THIS AND WHERE IS IT USED??? +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION + + +// routines used from cube_utils.c +extern double sh_find_overlap(const double xcenter, const double ycenter, + const double xlength, const double ylength, + double xPixelCorner[],double yPixelCorner[]); + + +extern double find_area_quad(double MinX, double MinY, double Xcorner[], double Ycorner[]); + + +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, + double *a1, double *a2,double *a3, double*a4, + double *lam1, double *lam2, double*lam3, double *lam4, + double *acoord, double *zcoord, int ss, + double *pixel_flux, double *pixel_err, + double **spaxel_flux, double **spaxel_weight, double **spaxel_var,double **spaxel_iflux){ + + double *fluxv, *weightv, *varv, *ifluxv ; // vectors for spaxel + + // allocate memory to hold output + if (alloc_flux_arrays(ncube, &fluxv, &weightv, &varv, &ifluxv)) return 1; + + + // loop over each valid point on detector and find match to IFU cube based + // on along slice coordinate and wavelength + for (int ipixel= 0; ipixel< npt; ipixel++){ + double along_corner[4]; + double wave_corner[4]; + + along_corner[0] = a1[ipixel]; + along_corner[1] = a2[ipixel]; + along_corner[2] = a3[ipixel]; + along_corner[3] = a4[ipixel]; + + wave_corner[0] = lam1[ipixel]; + wave_corner[1] = lam2[ipixel]; + wave_corner[2] = lam3[ipixel]; + wave_corner[3] = lam4[ipixel]; + + double along_min = 10000; + double wave_min = 10000; + double along_max = -10000; + double wave_max = -10000; + for (int j = 0; j< 4; j++){ + if(along_corner[j] < along_min) { along_min = along_corner[j];} + if(along_corner[j] > along_max) { along_max = along_corner[j];} + if(wave_corner[j] < wave_min) { wave_min = wave_corner[j];} + if(wave_corner[j] > wave_max) { wave_max = wave_corner[j];} + } + + double Area = find_area_quad(along_min, wave_min, along_corner, wave_corner); + + // estimate where the pixel overlaps in the cube + // find the min and max values in the cube appropriate xcoord,ycoord or zcoord + + int ia1 = (along_min - crval_along) / cdelt_along; + int ia2 = (along_max - crval_along) / cdelt_along; + if (ia1 < 0){ ia1 = 0;} + if (ia2 >= na){ ia2 = na -1;} + + double MinW = (wave_min - crval3) / cdelt3; + double MaxW = (wave_max - crval3) / cdelt3; + int iz1= (int) MinW; + int iz2 = round(MaxW); + + if(iz1 < 0){ iz1 = 0;} + if (iz2 >= nz){iz2 = nz - 1;} + + int nplane = naxis1 * naxis2; + // loop over possible overlapping cube pixels + for(int zz =iz1; zz < iz2+1; zz++){ + double zcenter = zcoord[zz]; + int istart = zz * nplane; + for (int aa= ia1; aa< ia2 + 1; aa++){ + int cube_index = 0; + if(instrument == 1) { // NIRSPec + cube_index = istart + aa * naxis1 + ss; // ss = slice # + } else { + cube_index = istart + ss * naxis1 + aa; // yy = slice # + } + double acenter = acoord[aa]; + double area_overlap = sh_find_overlap(acenter, zcenter, + cdelt_along, cdelt3, + along_corner, wave_corner); + + + if (area_overlap > 0.0) { + double AreaRatio = area_overlap / Area; + fluxv[cube_index] = fluxv[cube_index] + (AreaRatio * pixel_flux[ipixel]); + weightv[cube_index] = weightv[cube_index] + AreaRatio; + ifluxv[cube_index] = ifluxv[cube_index] + 1; + double err = (AreaRatio * pixel_err[ipixel]) * (AreaRatio * pixel_err[ipixel]); + varv[cube_index] = varv[cube_index] + err; + } + } + } + } + + *spaxel_flux = fluxv; + *spaxel_weight = weightv; + *spaxel_var = varv; + *spaxel_iflux = ifluxv; + + return 0; +} + + + +// C extension SETUP + +PyArrayObject * ensure_array(PyObject *obj, int *is_copy) { + if (PyArray_CheckExact(obj) && + PyArray_IS_C_CONTIGUOUS((PyArrayObject *) obj) && + PyArray_TYPE((PyArrayObject *) obj) == NPY_DOUBLE) { + *is_copy = 0; + return (PyArrayObject *) obj; + } else { + *is_copy = 1; + return (PyArrayObject *) PyArray_FromAny( + obj, PyArray_DescrFromType(NPY_DOUBLE), 0, 0, + NPY_ARRAY_CARRAY | NPY_ARRAY_FORCECAST, NULL + ); + } +} + + +static PyObject *cube_wrapper_internal(PyObject *module, PyObject *args) { + + PyObject *result = NULL, *a1o, *a2o, *a3o, *a4o, *lam1o, *lam2o, *lam3o, *lam4o, + *fluxo, *erro, *acoordo, *zcoordo; + + double crval_along, cdelt_along, crval3, cdelt3; + int nz, na, npt, naxis1, naxis2, ncube; + int ss; + int instrument; + + double *spaxel_flux=NULL, *spaxel_weight=NULL, *spaxel_var=NULL; + double *spaxel_iflux=NULL; + + int free_a1=0, free_a2=0, free_a3=0, free_a4 =0, free_lam1=0, free_lam2 =0, free_lam3=0, free_lam4=0; + int free_acoord=0, free_zcoord=0, status=0; + int free_flux=0, free_err=0; + + PyArrayObject *a1, *a2, *a3, *a4, *lam1, *lam2, *lam3, *lam4, *flux, *err, *acoord, *zcoord; + + PyArrayObject *spaxel_flux_arr=NULL, *spaxel_weight_arr=NULL, *spaxel_var_arr=NULL; + PyArrayObject *spaxel_iflux_arr=NULL; + npy_intp npy_ncube = 0; + + + if (!PyArg_ParseTuple(args, "iiiddddOOOOOOOOOOiOO:cube_wrapper_internal", + &instrument, &naxis1, &naxis2, &crval_along, &cdelt_along, + &crval3, &cdelt3, + &a1o, &a2o, &a3o, &a4o,&lam1o, &lam2o, &lam3o, &lam4o, + &acoordo, &zcoordo, &ss, &fluxo, &erro)){ + + return NULL; + } + + // check that input parameters are valid: + if ((cdelt3 < 0) || (cdelt_along < 0)) { + PyErr_SetString(PyExc_ValueError, + "'cdelt3' and 'cdelt_along' must be a strictly positive number."); + return NULL; + } + + // ensure we are working with numpy arrays and avoid creating new ones + // if possible: + if ((!(a1 = ensure_array(a1o, &free_a1))) || + (!(a2 = ensure_array(a2o, &free_a2))) || + (!(a3 = ensure_array(a3o, &free_a3))) || + (!(a4 = ensure_array(a4o, &free_a4))) || + (!(lam1 = ensure_array(lam1o, &free_lam1))) || + (!(lam2 = ensure_array(lam2o, &free_lam2))) || + (!(lam3 = ensure_array(lam3o, &free_lam3))) || + (!(lam4 = ensure_array(lam4o, &free_lam4))) || + (!(acoord = ensure_array(acoordo, &free_acoord))) || + (!(zcoord = ensure_array(zcoordo, &free_zcoord))) || + (!(flux = ensure_array(fluxo, &free_flux))) || + (!(err = ensure_array(erro, &free_err))) ) + + { + goto cleanup; + + } + + npt = (int) PyArray_Size((PyObject *) flux); + nz = (int) PyArray_Size((PyObject *) zcoord); + na = (int) PyArray_Size((PyObject *) acoord); + int n1 = (int) PyArray_Size((PyObject *) a1); + int n2 = (int) PyArray_Size((PyObject *) lam1); + + if (n1 != npt || n2 != npt ) { + PyErr_SetString(PyExc_ValueError, + "Input coordinate arrays of unequal size."); + goto cleanup; + + } + + ncube = naxis1 * naxis2 * nz; + + if (ncube ==0) { + // 0-length input arrays. Nothing to clip. Return 0-length arrays + spaxel_flux_arr = (PyArrayObject*) PyArray_EMPTY(1, &npy_ncube, NPY_DOUBLE, 0); + if (!spaxel_flux_arr) goto fail; + + spaxel_weight_arr = (PyArrayObject*) PyArray_EMPTY(1, &npy_ncube, NPY_DOUBLE, 0); + if (!spaxel_weight_arr) goto fail; + + spaxel_var_arr = (PyArrayObject*) PyArray_EMPTY(1, &npy_ncube, NPY_DOUBLE, 0); + if (!spaxel_var_arr) goto fail; + + spaxel_iflux_arr = (PyArrayObject*) PyArray_EMPTY(1, &npy_ncube, NPY_DOUBLE, 0); + if (!spaxel_iflux_arr) goto fail; + + result = Py_BuildValue("(OOOO)", spaxel_flux_arr, spaxel_weight_arr, spaxel_var_arr, + spaxel_iflux_arr); + + goto cleanup; + } + + //______________________________________________________________________ + // Match the point cloud elements to the spaxels they fail within the roi + //______________________________________________________________________ + status = match_detector_cube(instrument, naxis1, naxis2, nz, npt, ncube, na, + crval_along, cdelt_along, crval3, cdelt3, + (double *) PyArray_DATA(a1), + (double *) PyArray_DATA(a2), + (double *) PyArray_DATA(a3), + (double *) PyArray_DATA(a4), + (double *) PyArray_DATA(lam1), + (double *) PyArray_DATA(lam2), + (double *) PyArray_DATA(lam3), + (double *) PyArray_DATA(lam4), + (double *) PyArray_DATA(acoord), + (double *) PyArray_DATA(zcoord), + ss, + (double *) PyArray_DATA(flux), + (double *) PyArray_DATA(err), + &spaxel_flux, &spaxel_weight, &spaxel_var, &spaxel_iflux); + + if (status ) { + goto fail; + + } else { + // create return tuple: + npy_ncube = (npy_intp) ncube; + + spaxel_flux_arr = (PyArrayObject*) PyArray_SimpleNewFromData(1, &npy_ncube, NPY_DOUBLE, spaxel_flux); + if (!spaxel_flux_arr) goto fail; + spaxel_flux = NULL; + + spaxel_weight_arr = (PyArrayObject*) PyArray_SimpleNewFromData(1, &npy_ncube, NPY_DOUBLE, spaxel_weight); + if (!spaxel_weight_arr) goto fail; + spaxel_weight = NULL; + + spaxel_var_arr = (PyArrayObject*) PyArray_SimpleNewFromData(1, &npy_ncube, NPY_DOUBLE, spaxel_var); + if (!spaxel_var_arr) goto fail; + spaxel_var = NULL; + + spaxel_iflux_arr = (PyArrayObject*) PyArray_SimpleNewFromData(1, &npy_ncube, NPY_DOUBLE, spaxel_iflux); + if (!spaxel_iflux_arr) goto fail; + spaxel_iflux = NULL; + + PyArray_ENABLEFLAGS(spaxel_flux_arr, NPY_ARRAY_OWNDATA); + PyArray_ENABLEFLAGS(spaxel_weight_arr, NPY_ARRAY_OWNDATA); + PyArray_ENABLEFLAGS(spaxel_var_arr, NPY_ARRAY_OWNDATA); + PyArray_ENABLEFLAGS(spaxel_iflux_arr, NPY_ARRAY_OWNDATA); + + result = Py_BuildValue("(OOOO)", spaxel_flux_arr, spaxel_weight_arr, spaxel_var_arr, + spaxel_iflux_arr); + + goto cleanup; + } + + fail: + Py_XDECREF(spaxel_flux_arr); + Py_XDECREF(spaxel_weight_arr); + Py_XDECREF(spaxel_var_arr); + Py_XDECREF(spaxel_iflux_arr); + + free(spaxel_flux); + free(spaxel_weight); + free(spaxel_var); + free(spaxel_iflux); + + + if (!PyErr_Occurred()) { + PyErr_SetString(PyExc_MemoryError, + "Unable to allocate memory for output arrays."); + } + + cleanup: + if (free_a1)Py_XDECREF(a1); + if (free_a2)Py_XDECREF(a2); + if (free_a3)Py_XDECREF(a3); + if (free_a4)Py_XDECREF(a4); + if (free_lam1)Py_XDECREF(lam1); + if (free_lam2)Py_XDECREF(lam2); + if (free_lam3)Py_XDECREF(lam3); + if (free_lam4)Py_XDECREF(lam4); + if (free_acoord) Py_XDECREF(acoord); + if (free_zcoord) Py_XDECREF(zcoord); + if (free_flux) Py_XDECREF(flux); + if (free_err) Py_XDECREF(err); + + return result; +} + +static PyMethodDef cube_methods[] = +{ + { + "cube_wrapper_internal", + cube_wrapper_internal, + METH_VARARGS, + "cube_wrapper_internal(put in doc string)" + }, + {NULL, NULL, 0, NULL} /* sentinel */ +}; + +static struct PyModuleDef moduledef = { + PyModuleDef_HEAD_INIT, + "cube_match_internal", /* m_name */ + "find point cloud matches for each spaxel center", /* m_doc */ + -1, /* m_size */ + cube_methods, /* m_methods */ + NULL, /* m_reload */ + NULL, /* m_traverse */ + NULL, /* m_clear */ + NULL, /* m_free */ +}; + +PyMODINIT_FUNC PyInit_cube_match_internal(void) +{ + PyObject* m; + import_array(); + m = PyModule_Create(&moduledef); + return m; +} diff --git a/jwst/cube_build/src/cube_match_sky.c b/jwst/cube_build/src/cube_match_sky.c new file mode 100644 index 0000000000..07910ae52f --- /dev/null +++ b/jwst/cube_build/src/cube_match_sky.c @@ -0,0 +1,1351 @@ +/* +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 zcdelt3. +This module uses the modified shephard 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 + +Python signature: result = cube_wrapper(instrument, flag_dq_plane, weight_type, 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) +provide more details + +The output of this function is a tuple of 5 arrays:(spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux, spaxel_dq) +example output + +Parameters +---------- +instrument : int + 0 = MIRI, 1 = NIRSPEC. Used for set the dq plane +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 +weight_type : int + 0: use emsm weighting + 1: use msm weighting +start_region : int + starting slice number for detector region used in dq flagging +end_region: int + 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 +overlap_full : int + 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 +ycoord : double array + 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 + size: point cloud elements. Flux of each point cloud memeber +err : double array + size: point cloud elements. err of each point cloud memeber +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. Wavelegnth of each point cloud memeber +rois_pixel : double array + size: point cloud elements. Roi in spatial dimension to use for point cloud memeber +roiw_pixel : double array + size: point cloud elements. Roi in wavelength dimension to use point cloud memeber +scalerad_pixel : double array + size: point cloud elements. 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 +roiw_ave : double + Average roiw for all the wavelength planes. Used in dq flagging +cdelt1 : double + Naxis 1 scale for cube +cdelt2 : double + Naxis 2 scale for cube + + +Returns +------- +spaxel_flux : numpy.ndarray + IFU spaxel cflux +spaxel_weight : numpy.ndarray + IFU spaxel weight +spaxel_iflux : numpy.ndarray + IFU spaxel weight map (number of overlaps) +spaxel_var : numpy.ndarray + IFU spaxel error +spaxel_dq : numpy.ndarray + IFU spaxel dq +*/ + +#include +#include +#include +#include +#include +#include +#include + +#define PY_ARRAY_UNIQUE_SYMBOL _jwst_cube_match_sky_numpy_api //WHAT IS THIS AND WHERE IS IT USED??? +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION + +// routines used from cube_utils.c + +extern double sh_find_overlap(const double xcenter, const double ycenter, + const double xlength, const double ylength, + double xPixelCorner[],double yPixelCorner[]); + +extern int alloc_flux_arrays(int nelem, double **fluxv, double **weightv, double **varv, double **ifluxv); + + +//________________________________________________________________________________ +// allocate the memory for the spaxel DQ array + +int mem_alloc_dq(int nelem, int **idqv) { + + const char *msg = "Couldn't allocate memory for output arrays."; + + if (!(*idqv = (int*)calloc(nelem, sizeof(int)))) { + PyErr_SetString(PyExc_MemoryError, msg); + return 1; + } + + return 0; +} + +//________________________________________________________________________________ +// Routine for MIRI DQ plane assignment + +int corner_wave_plane_miri(int w, int start_region, int end_region, + double roiw_ave, + double *zc, + double *coord1, double *coord2, double *wave, + double *sliceno, + long ncube, int 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 + corners of the FOV. + + w : wavelength plane + start_region : starting slice # for channel (slice # in nirspec) + end_region : ending slice # for chanel (slice # in nirspec) + roiw_ave : average roiw for all wavelengths + zc : array of wavelenths + 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 + */ + + int status = 0; + int ic1_start_min = -1; + int ic1_start_max = -1; + int ic2_start_min = -1; + int ic2_start_max = -1; + + float c1_start_min = 10000; + float c2_start_min = 10000; + float c1_start_max = -10000; + float c2_start_max = -10000; + + int ic1_end_min = -1; + int ic1_end_max = -1; + int ic2_end_min = -1; + int ic2_end_max = -1; + + float c1_end_min = 10000; + float c2_end_min = 10000; + float c1_end_max = -10000; + float c2_end_max = -10000; + + // Loop over every point cloud member and pull out values + // 1. Fail withing roiw_ave of wavelength plane + // and + // 2. Are for either of the 2 extreme slices + + for (int ipt =0; ipt< npt ; ipt++){ + int slice = (int)sliceno[ipt]; + double wave_distance = fabs(zc[w] - wave[ipt]); + + float c11 = -1; // coord1 for start region + float c21 = -1; // coord2 for start region + float c12 = -1; // coord1 for end region + float c22 = -1; // corrd2 for end region + + if(slice == start_region || slice==end_region){ + // Find all the coordinates on wave slice with slice = start region + // These points will define corner 1 (min c2) and corner 2 (max c2) + if (wave_distance < roiw_ave && slice == start_region){ + c11 = coord1[ipt]; + c21 = coord2[ipt]; + } + + // Find all the coordinates on wave slice with slice = start region + // These points will define corner 4 (min c2) and corner 3 (max c2) + if (wave_distance < roiw_ave && slice == end_region){ + c12 = coord1[ipt]; + c22 = coord2[ipt]; + } + + // for the start region find min and max c1,c2 + if( c11 !=-1){ + // check c1 points for min + if (c11 < c1_start_min) { + c1_start_min = c11; + ic1_start_min = ipt; + } + // check c1 points for max + if (c11 > c1_start_max) { + c1_start_max= c11; + ic1_start_max = ipt; + } + // check c2 points for min + if (c21 < c2_start_min) { + c2_start_min = c21; + ic2_start_min = ipt; + } + // check c2 points for max + if (c21 > c2_start_max) { + c2_start_max= c21; + ic2_start_max = ipt; + } + } + // for the end region find min and max c1,c2 + if( c12 !=-1){ + if (c12 < c1_end_min) { + c1_end_min = c12; + ic1_end_min = ipt; + } + if (c12 > c1_end_max) { + c1_end_max = c12; + ic1_end_max = ipt; + } + if (c22 < c2_end_min) { + c2_end_min = c22; + ic2_end_min = ipt; + } + if (c22 > c2_end_max) { + c2_end_max = c22; + ic2_end_max = ipt; + } + } + } + } // end looping over point cloud + + // Make sure the 2 extreme slices are found on the FOV. Not finding both can occur for edge wavelength planes + // or empty wavelength planes between channels + + if( ic1_start_min ==-1 || ic1_start_max ==-1 || ic1_end_min == -1 || ic1_end_max == -1){ + status = 1; + return status; + } else { + // Find the length in the start and end regions for c1 and c2. This will help define which coordinates to use to set corners. + // Because we do not know the orientation on the sky pick the longest length to set how to pick corners. + float length_c1_start = c1_start_max - c1_start_min; + float length_c2_start = c2_start_max - c2_start_min; + //float length_c1_end = c1_end_max - c1_end_min; + // float length_c2_end = c2_end_max - c2_end_min; + + int c1_use = 1; // use the c1 coords to set corners + if(length_c1_start < length_c2_start){ + c1_use = 0; // use the c2 coords to set the corners + } + + if(c1_use ==0) { + corner1[0] = coord1[ic2_start_min]; + corner1[1] = coord2[ic2_start_min]; + + corner2[0] = coord1[ic2_start_max]; + corner2[1] = coord2[ic2_start_max]; + + corner3[0] = coord1[ic2_end_max]; + corner3[1] = coord2[ic2_end_max]; + + corner4[0] = coord1[ic2_end_min]; + corner4[1] = coord2[ic2_end_min]; + } else{ + corner1[0] = coord1[ic1_start_min]; + corner1[1] = coord2[ic1_start_min]; + + corner2[0] = coord1[ic1_start_max]; + corner2[1] = coord2[ic1_start_max]; + + corner3[0] = coord1[ic1_end_max]; + corner3[1] = coord2[ic1_end_max]; + + corner4[0] = coord1[ic1_end_min]; + corner4[1] = coord2[ic1_end_min]; + } + + status = 0; + return status; + } +} + +//________________________________________________________________________________ +// 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, + int naxis1, int naxis2, + double xcenters[], double ycenters[], + double xi_corner[], double eta_corner[], + int wave_slice_dq[]) { + + /* find the amount of overlap of FOV each spaxel + + Given the corners of the FOV find the spaxels that + 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 + dithered FOVs. + + Parameters + ---------- + overlap_partial : int + overlap_full : int + cdelt1 : double + IFU cube naxis 1 spatial scale + cdelt2 : double + IFU cube naxis 1 spatial scale + naxis1 : int + IFU cube naxis 1 size + naxis2 : int + IFU cube naxis 1 size + xcenters : double array + 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 wavelenghth plane + eta_corner: eta coordinates of the 4 corners of the FOV on the wavelength plane + + + Sets + ------- + wave_slice_dq: array containing intermediate dq flag + + */ + + // loop over spaxels in the wavelength plane and set slice_dq + // roughly find the spaxels that might be overlapped + + int status = 0; + double ximin = 1000.0; + double etamin = 1000.0; + double ximax = -10000.0; + double etamax = -1000.0; + for (int i=0; i<4 ; i++){ + if (xi_corner[i] < ximin){ ximin = xi_corner[i];} + if (xi_corner[i] > ximax){ ximax = xi_corner[i];} + if (eta_corner[i] < etamin){ etamin = eta_corner[i];} + if (eta_corner[i] < etamax){ etamax = eta_corner[i];} + } + + double area_box = cdelt1 * cdelt2; + float tolerance_dq_overlap = 0.05; //spaxel has to have 5% overlap to flag in FOV + // loop over cube xcenters and cube ycenters + for (int ix = 0; ix < naxis1; ix++){ + double x1 = (xcenters[ix] - cdelt1)/2; + double x2 = (xcenters[ix] + cdelt1)/2; + if(x1 > ximin && x2 < ximax){ + for (int iy = 0; iy< naxis2; iy++){ + + double y1 = (ycenters[ix] - cdelt2)/2; + double y2 = (ycenters[ix] + cdelt2)/2; + if(y1 > etamin && y2 < etamax){ + int ixy = iy* naxis1 + ix; + double area_overlap = sh_find_overlap(xcenters[ix],ycenters[iy], + cdelt1, cdelt2, + xi_corner, eta_corner); + + double overlap_coverage = area_overlap / area_box; + + if (overlap_coverage > tolerance_dq_overlap){ + if (overlap_coverage > 0.95) { + wave_slice_dq[ixy] = overlap_full; + }else{ + wave_slice_dq[ixy] = overlap_partial; + } + } + }// end y1, y2 test + }// end loop over iy + } // end x1, x2 test + } // end loop over ix + + return status ; +} + +//________________________________________________________________________________ +// Routine to setting NIRSpec dq plane for each wavelength plane + +int slice_wave_plane_nirspec(int w, int slicevalue, + double roiw_ave, + double *zc, + double *coord1, double *coord2, double *wave, + double *sliceno, + long ncube, int npt, + double *c1_min, double *c1_max, double *c2_min, double *c2_max) { + /* + + NIRSpec dq plane is set by mapping each slice to IFU wavelength plane + This routine maps each slice to sky and finds the min and max coordinates on the sky + of the slice. + + w : wavelength plane + slicevalue : slice # 1 to 30 + roiw_ave : average roiw for all wavelengths + zc : array of wavelenths + coord1 : point cloud xi values + coord2 : point cloud eta values + wave : point cloud wavelength values + sliceno : point cloud slice no - starts at 1 + ncube : number of cube values + npt: number of point cloud elements + + return: + c1_min, c1_max, c2_min, c2_max + */ + + double dvalue = 10000; + *c1_min = dvalue; + *c2_min = dvalue; + *c1_max = -dvalue; + *c2_max = -dvalue; + + for (int ipt =0; ipt< npt ; ipt++){ + int slice = (int)sliceno[ipt]; + double wave_distance = fabs(zc[w] - wave[ipt]); + + // Find all the coordinates on wave slice with slice = start region + + if (wave_distance < roiw_ave && slice == slicevalue){ + + double c1 = coord1[ipt]; + double c2 = coord2[ipt]; + // find min, max of xi eta + if (c1 < *c1_min) {*c1_min = c1;} + if (c1 > *c1_max) {*c1_max = c1;} + if (c2 < *c2_min) {*c2_min = c2;} + if (c2 > *c2_max) {*c2_max = c2;} + } + + } // end looping over point cloud + + int status = 0; + if(*c1_min == dvalue || *c2_min == dvalue || *c1_max ==-dvalue || *c2_max == -dvalue){ + // Problem finding limits of slice for wavelength plane + // This is likely caused the no valid data on wavelength plane + // The two ends of wavelengths can have DQ detector data set to DO_NOT_USE - setting up no + // valid data on thee planes in the IFU Cube. + + status = 1; + } + + return status; +} + +//________________________________________________________________________________ +// 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, + double xcenters[], double ycenters[], + double xi_min, double eta_min, + double xi_max, double eta_max, + int wave_slice_dq[]) { + + /* + Set the 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 + which IFU spaxels the slice falls on to set an initial dq flag. + + Parameters + --------- + overlap_partial: intermediate dq flag + + 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 + ------- + Points: a tuple of x,y spaxel values that this line intersects + + */ + + int status = 0; + + //set up line - convert to integer values + int x1 = (int)((xi_min - xcenters[0]) / cdelt1); + int y1 = (int)((eta_min - ycenters[0]) / cdelt2); + int x2 = (int)((xi_max - xcenters[0]) / cdelt1); + int y2 = (int)((eta_max - ycenters[0]) / cdelt2); + + int dx = x2 - x1; + int dy = y2 - y1; + + bool is_steep; + is_steep = abs(dy) > abs(dx); + + // if is_steep switch x and y + if (is_steep){ + x1 = y1; + y1 = x1; + x2 = y2; + y2 = x2; + } + + // Swap start and end points if necessary and store swap state + + if (x1 > x2){ + x1 = x2; + x2 = x1; + + y1 = y2; + y2 = y1; + } + + // Recalculate differences + dx = x2 - x1; + dy = y2 - y1; + + //calculate error + int error = (int)(dx / 2.0); + int ystep = -1; + if (y1 < y2){ + ystep = 1; + } + + // iterate over grid to generate points between the start and end of line + int y = y1; + for (int x = x1; x< (x2 + 1); x++){ + int yuse = y; + int xuse = x ; + if (is_steep){ + yuse = x; + xuse = y; + } + + int index = (yuse * naxis1) + xuse; + wave_slice_dq[index] = overlap_partial; + error -= abs(dy); + if (error < 0){ + y += ystep; + error += dx; + } + } + return status; +} + +//________________________________________________________________________________ +// 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){ + + int *idqv; // int vector for spaxel + if (mem_alloc_dq(ncube, &idqv)) return 1; + *spaxel_dq = idqv; + return 0; +} + +//________________________________________________________________________________ +// 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, + double cdelt1, double cdelt2, double roiw_ave, + double *xc, double *yc, double *zc, + double *coord1, double *coord2, double *wave, + double *sliceno, + long ncube, int npt, + int **spaxel_dq) { + + int *idqv ; // int vector for spaxel + if (mem_alloc_dq(ncube, &idqv)) return 1; + + double corner1[2]; + double corner2[2]; + double corner3[2]; + double corner4[2]; + + // for each wavelength plane find the 2 extreme slices to set FOV. Use these two extreme slices to set up the + // corner of the FOV for each wavelength + + int nxy = nx * ny; + // Loop over the wavelength planes and set DQ plane + for (int w = 0; w < nz; w++) { + + int wave_slice_dq[nxy]; + for( int i = 0; i < nxy; i ++){ + wave_slice_dq[i] = 0; + } + + int status = corner_wave_plane_miri( w, start_region, end_region, roiw_ave, zc, + coord1, coord2, wave, sliceno, ncube, npt, + corner1, corner2, corner3, corner4); + if( status == 0){ // found min and max slice on wavelengh plane + double xi_corner[4]; + double eta_corner[4]; + xi_corner[0] = corner1[0]; + xi_corner[1] = corner2[0]; + xi_corner[2] = corner3[0]; + xi_corner[3] = corner4[0]; + + eta_corner[0] = corner1[1]; + eta_corner[1] = corner2[1]; + eta_corner[2] = corner3[1]; + eta_corner[3] = corner4[1]; + + status = overlap_fov_with_spaxels(overlap_partial, overlap_full, + cdelt1,cdelt2, + nx, ny, + xc, yc, + xi_corner, eta_corner, + wave_slice_dq); + long istart = nxy*w; + long iend = istart + nxy; + for( long in = istart; in < iend; in ++){ + long ii = in - istart; + idqv[in] = wave_slice_dq[ii]; + } + } + } // end loop over wavelength + + *spaxel_dq = idqv; + + 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, + double *xc, double *yc, double *zc, + double *coord1, double *coord2, double *wave, + double *sliceno, + long ncube, int npt, + int **spaxel_dq) { + + /* + Set an initial 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 + range of wavelengths. The FOV of the slice is really just a line, so instead of using + the routines the finds the overlap between a polygon and regular grid- + which is used for MIRI - an algorithm that determines the spaxels that + the slice line intersects is used instead. + + Paramteter + --------- + overlap_partial - intermediate DQ flag + 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 + slice_no: integer slice value of input data (used in MIRI case to find + the points of the edge slices.) + + */ + + + int *idqv ; // int vector for spaxel + if (mem_alloc_dq(ncube, &idqv)) return 1; + + + // for each of the 30 slices - find the projection of this slice + // onto each of the IFU wavelength planes. + + for (int w = 0; w < nz; w++) { + int n_slice_found = 0; + for (int islice =1; islice< 31 ; islice++){ + double c1_min; + double c2_min; + double c1_max; + double c2_max; + int status = 0; + status = slice_wave_plane_nirspec( w, islice, roiw_ave, zc, + coord1, coord2, wave, sliceno, ncube, npt, + &c1_min, &c2_min, &c1_max, &c2_max); + + if( status ==0){ + + n_slice_found++; + int nxy = nx * ny; + int wave_slice_dq[nxy]; + for (int j =0; j< nxy; j++){ + wave_slice_dq[j] = 0; + } + status = overlap_slice_with_spaxels(overlap_partial, + cdelt1,cdelt2, + nx, ny, + xc, yc, + c1_min, c2_min, c1_max, c2_max, + wave_slice_dq); + + long istart = nxy*w; + long iend = istart + nxy; + for( long in = istart; in < iend; in ++){ + long ii = in - istart; + idqv[in] = wave_slice_dq[ii]; + } + } // end loop over status + + } // end loop over slices + + } // end of wavelength + *spaxel_dq = idqv; + + return 0; +} + +// 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 + +// return values: spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux + +int match_point_emsm(double *xc, double *yc, double *zc, + double *coord1, double *coord2, double *wave, + double *flux, double *err, + double *rois_pixel, double *roiw_pixel, double *scalerad_pixel, + double *zcdelt3, + int nx, int ny, int nwave, int ncube, int npt, + double cdelt1, double cdelt2, + double **spaxel_flux, double **spaxel_weight, double **spaxel_var, + double **spaxel_iflux) { + + + double *fluxv, *weightv, *varv, *ifluxv; // vector for spaxel + + // allocate memory to hold output + if (alloc_flux_arrays(ncube, &fluxv, &weightv, &varv, &ifluxv)) return 1; + + + // loop over each point cloud member and find which roi spaxels it is found + + for (int k = 0; k < npt; k++) { + // Search wave and find match + int iwstart = -1; + int iwend = -1; + int ii = 0; + int done_search_w = 0; + + while (ii < nwave && done_search_w == 0) { + float wdiff = fabs(zc[ii] - wave[k]); + if(wdiff <= roiw_pixel[k]){ + if (iwstart == -1){ + iwstart = ii; + } + } else{ + if(iwstart != -1 && iwend ==-1){ + iwend = ii; + done_search_w = 1; + } + } + ii = ii + 1; + } + // catch the case of iwstart near nwave and becomes = nwave before iwend can be set. + if(iwstart !=-1 && iwend == -1){ + iwend = nwave; + done_search_w = 1; + } + + + // Search xcenters and find match + int ixstart = -1; + int ixend = -1; + ii = 0; + int done_search_x = 0; + while (ii < nx && done_search_x == 0) { + double xdiff = fabs(xc[ii] - coord1[k]); + if(xdiff <= rois_pixel[k]){ + if (ixstart == -1){ + ixstart = ii; + } + } else{ + if(ixstart != -1 && ixend ==-1){ + ixend = ii; + done_search_x = 1; + } + } + ii = ii + 1; + } + // catch the case of ixstart near nx and becomes = nx before ixend can be set. + if(ixstart !=-1 && ixend == -1){ + ixend = nx; + done_search_x = 1; + } + + // Search ycenters and find match + int iystart = -1; + int iyend = -1; + ii = 0; + int done_search_y = 0; + while (ii < ny && done_search_y == 0) { + double ydiff = fabs(yc[ii] - coord2[k]); + if(ydiff <= rois_pixel[k]){ + if (iystart == -1){ + iystart = ii; + } + } else{ + if(iystart != -1 && iyend ==-1){ + iyend = ii; + done_search_y = 1; + } + } + ii = ii + 1; + } + // catch the case of iystart near ny and becomes = ny before iyend can be set. + if(iystart !=-1 && iyend == -1){ + iyend = ny; + done_search_y = 1; + } + + // set up the values for fluxv, weightv, ifluxv, varv + int nxy = nx * ny; + if(done_search_x == 1 && done_search_y ==1 && done_search_w ==1){ + // The search above for x,y was a crude search - now narrow the search using the distance between + // the spaxel center and point cloud + for (int ix = ixstart; ix< ixend; ix ++){ + for ( int iy = iystart; iy < iyend; iy ++){ + double ydist = fabs(yc[iy] - coord2[k]); + double xdist = fabs(xc[ix] - coord1[k]); + double radius = sqrt( xdist*xdist + ydist*ydist); + + if (radius <= rois_pixel[k]){ + // Find the index for this in spatial plane + int index_xy = iy* nx + ix; + for (int iw = iwstart; iw< iwend; iw++){ + int index_cube = iw*nxy + index_xy; + + double d1 = xdist/cdelt1; + double d2 = ydist/cdelt2; + double dxy = (d1 * d1) + (d2 * d2); + double d3 = (wave[k] - zc[iw])/ zcdelt3[iw]; + //double wdist = fabs( wave[k] - zc[iw]); + double d32 = d3 * d3; + double w = d32 + dxy; + double wn = -w/(scalerad_pixel[k]/cdelt1); + double ww = exp(wn); + + double weighted_flux = flux[k]* ww; + double weighted_var = (err[k]* ww) * (err[k]*ww); + fluxv[index_cube] = fluxv[index_cube] + weighted_flux; + weightv[index_cube] = weightv[index_cube] + ww; + varv[index_cube] = varv[index_cube] + weighted_var; + ifluxv[index_cube] = ifluxv[index_cube] +1.0; + + } + } + } // end loop over iy + } // end loop over ix + + } // end done_search_x, done_search_y, done_search_w + } // end loop over point cloud + + + // assign output values: + + *spaxel_flux = fluxv; + *spaxel_weight = weightv; + *spaxel_var = varv; + *spaxel_iflux = ifluxv; + + return 0; +} + + + +// 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 + +int match_point_msm(double *xc, double *yc, double *zc, + double *coord1, double *coord2, double *wave, + double *flux, double *err, + double *rois_pixel, double *roiw_pixel, + double *weight_pixel, double *softrad_pixel, + double *zcdelt3, + int nx, int ny, int nwave, int ncube, int npt, + double cdelt1, double cdelt2, + double **spaxel_flux, double **spaxel_weight, double **spaxel_var, + double **spaxel_iflux) { + + + double *fluxv, *weightv, *varv, *ifluxv; // vector for spaxel + + // allocate memory to hold output + if (alloc_flux_arrays(ncube, &fluxv, &weightv, &varv, &ifluxv)) return 1; + + // loop over each point cloud member and find which roi spaxels it is found + + for (int k = 0; k < npt; k++) { + // Search wave and find match + int iwstart = -1; + int iwend = -1; + int ii = 0; + int done_search_w = 0; + + while (ii < nwave && done_search_w == 0) { + float wdiff = fabs(zc[ii] - wave[k]); + if(wdiff <= roiw_pixel[k]){ + if (iwstart == -1){ + iwstart = ii; + } + } else{ + if(iwstart != -1 && iwend ==-1){ + iwend = ii; + done_search_w = 1; + } + } + ii = ii + 1; + } + // catch the case of iwstart near nwave and becomes = nwave before iwend can be set. + if(iwstart !=-1 && iwend == -1){ + iwend = nwave; + done_search_w = 1; + } + + + // Search xcenters and find match + int ixstart = -1; + int ixend = -1; + ii = 0; + int done_search_x = 0; + while (ii < nx && done_search_x == 0) { + double xdiff = fabs(xc[ii] - coord1[k]); + if(xdiff <= rois_pixel[k]){ + if (ixstart == -1){ + ixstart = ii; + } + } else{ + if(ixstart != -1 && ixend ==-1){ + ixend = ii; + done_search_x = 1; + } + } + ii = ii + 1; + } + // catch the case of ixstart near nx and becomes = nx before ixend can be set. + if(ixstart !=-1 && ixend == -1){ + ixend = nx; + done_search_x = 1; + } + + // Search ycenters and find match + int iystart = -1; + int iyend = -1; + ii = 0; + int done_search_y = 0; + while (ii < ny && done_search_y == 0) { + double ydiff = fabs(yc[ii] - coord2[k]); + if(ydiff <= rois_pixel[k]){ + if (iystart == -1){ + iystart = ii; + } + } else{ + if(iystart != -1 && iyend ==-1){ + iyend = ii; + done_search_y = 1; + } + } + ii = ii + 1; + } + // catch the case of iystart near ny and becomes = ny before iyend can be set. + if(iystart !=-1 && iyend == -1){ + iyend = ny; + done_search_y = 1; + } + + // set up the values for fluxv, weightv, ifluxv, varv + int nxy = nx * ny; + if(done_search_x == 1 && done_search_y ==1 && done_search_w ==1){ + // The search above for x,y was a crude search - now narrow the search using the distance between + // the spaxel center and point cloud + for (int ix = ixstart; ix< ixend; ix ++){ + for ( int iy = iystart; iy < iyend; iy ++){ + double ydist = fabs(yc[iy] - coord2[k]); + double xdist = fabs(xc[ix] - coord1[k]); + double radius = sqrt( xdist*xdist + ydist*ydist); + + if (radius <= rois_pixel[k]){ + // Find the index for this in spatial plane + int index_xy = iy* nx + ix; + for (int iw = iwstart; iw< iwend; iw++){ + int index_cube = iw*nxy + index_xy; + + double d1 = xdist/cdelt1; + double d2 = ydist/cdelt2; + double dxy = (d1 * d1) + (d2 * d2); + double d3 = (wave[k] - zc[iw])/ zcdelt3[iw]; + + double d32 = d3 * d3; + double w = d32 + dxy; + double wn = pow(sqrt(w), weight_pixel[k]); + if( wn < softrad_pixel[k]){ + wn = softrad_pixel[k]; + } + + double ww = 1.0/wn; + double weighted_flux = flux[k]* ww; + double weighted_var = (err[k]* ww) * (err[k]*ww); + fluxv[index_cube] = fluxv[index_cube] + weighted_flux; + weightv[index_cube] = weightv[index_cube] + ww; + varv[index_cube] = varv[index_cube] + weighted_var; + ifluxv[index_cube] = ifluxv[index_cube] +1.0; + + } + } + } // end loop over iy + } // end loop over ix + + } // end done_search_x, done_search_y, done_search_w + } // end loop over point cloud + + + // assign output values: + + *spaxel_flux = fluxv; + *spaxel_weight = weightv; + *spaxel_var = varv; + *spaxel_iflux = ifluxv; + + return 0; +} + + +PyArrayObject * ensure_array(PyObject *obj, int *is_copy) { + if (PyArray_CheckExact(obj) && + PyArray_IS_C_CONTIGUOUS((PyArrayObject *) obj) && + PyArray_TYPE((PyArrayObject *) obj) == NPY_DOUBLE) { + *is_copy = 0; + return (PyArrayObject *) obj; + } else { + *is_copy = 1; + return (PyArrayObject *) PyArray_FromAny( + obj, PyArray_DescrFromType(NPY_DOUBLE), 0, 0, + NPY_ARRAY_CARRAY | NPY_ARRAY_FORCECAST, NULL + ); + } +} + + +static PyObject *cube_wrapper(PyObject *module, PyObject *args) { + PyObject *result = NULL, *xco, *yco, *zco, *fluxo, *erro, *coord1o, *coord2o, *waveo, *slicenoo; + PyObject *rois_pixelo, *roiw_pixelo, *scalerad_pixelo, *zcdelt3o, *softrad_pixelo, *weight_pixelo; + + double cdelt1, cdelt2, roiw_ave; + int nwave, npt, nxx, nyy, ncube; + + int instrument, flag_dq_plane,start_region, end_region, overlap_partial, overlap_full, weight_type; + double *spaxel_flux=NULL, *spaxel_weight=NULL, *spaxel_var=NULL; + double *spaxel_iflux=NULL; + int *spaxel_dq=NULL; + + int free_xc=0, free_yc=0, free_zc=0, free_coord1=0, free_coord2 =0 , free_wave=0, status=0; + int free_rois_pixel=0, free_roiw_pixel=0, free_scalerad_pixel=0, free_flux=0, free_err=0, free_zcdelt3=0; + int free_sliceno=0, free_softrad_pixel, free_weight_pixel=0; + + PyArrayObject *xc, *yc, *zc, *flux, *err, *coord1, *coord2, *wave, *rois_pixel, *roiw_pixel, *scalerad_pixel; + PyArrayObject *zcdelt3, *sliceno, *softrad_pixel, *weight_pixel; + PyArrayObject *spaxel_flux_arr=NULL, *spaxel_weight_arr=NULL, *spaxel_var_arr=NULL; + PyArrayObject *spaxel_iflux_arr=NULL, *spaxel_dq_arr=NULL; + npy_intp npy_ncube = 0; + + int ny,nz; + + if (!PyArg_ParseTuple(args, "iiiiiiiOOOOOOOOOOOOOOOddd:cube_wrapper", + &instrument, &flag_dq_plane, &weight_type, &start_region, &end_region, &overlap_partial, &overlap_full, + &xco, &yco, &zco, &coord1o, &coord2o, &waveo, &fluxo, &erro, &slicenoo, + &rois_pixelo, &roiw_pixelo, &scalerad_pixelo, &weight_pixelo, &softrad_pixelo, &zcdelt3o, &roiw_ave, + &cdelt1, &cdelt2)) { + return NULL; + } + + + // check that input parameters are valid: + + if ((cdelt1 < 0) || (cdelt2 < 0)) { + PyErr_SetString(PyExc_ValueError, + "'cdelt1' and 'cdelt2' must be a strictly positive number."); + return NULL; + } + + // ensure we are working with numpy arrays and avoid creating new ones + // if possible: + if ((!(xc = ensure_array(xco, &free_xc))) || + (!(yc = ensure_array(yco, &free_yc))) || + (!(zc = ensure_array(zco, &free_zc))) || + (!(coord1 = ensure_array(coord1o, &free_coord1))) || + (!(coord2 = ensure_array(coord2o, &free_coord2))) || + (!(wave = ensure_array(waveo, &free_wave))) || + (!(flux = ensure_array(fluxo, &free_flux))) || + (!(err = ensure_array(erro, &free_err))) || + (!(sliceno = ensure_array(slicenoo, &free_sliceno))) || + (!(rois_pixel = ensure_array(rois_pixelo, &free_rois_pixel))) || + (!(roiw_pixel = ensure_array(roiw_pixelo, &free_roiw_pixel))) || + (!(zcdelt3 = ensure_array(zcdelt3o, &free_zcdelt3))) || + (!(scalerad_pixel = ensure_array(scalerad_pixelo, &free_scalerad_pixel))) || + (!(softrad_pixel = ensure_array(softrad_pixelo, &free_softrad_pixel))) || + (!(weight_pixel = ensure_array(weight_pixelo, &free_weight_pixel))) ) + { + goto cleanup; + + } + + npt = (int) PyArray_Size((PyObject *) coord1); + ny = (int) PyArray_Size((PyObject *) coord2); + nz = (int) PyArray_Size((PyObject *) wave); + if (ny != npt || nz != npt ) { + PyErr_SetString(PyExc_ValueError, + "Input coordinate arrays of unequal size."); + goto cleanup; + + } + + nxx = (int) PyArray_Size((PyObject *) xc); + nyy = (int) PyArray_Size((PyObject *) yc); + nwave = (int) PyArray_Size((PyObject *) zc); + + ncube = nxx * nyy * nwave; + + if (ncube ==0) { + // 0-length input arrays. Nothing to clip. Return 0-length arrays + spaxel_flux_arr = (PyArrayObject*) PyArray_EMPTY(1, &npy_ncube, NPY_DOUBLE, 0); + if (!spaxel_flux_arr) goto fail; + + spaxel_weight_arr = (PyArrayObject*) PyArray_EMPTY(1, &npy_ncube, NPY_DOUBLE, 0); + if (!spaxel_weight_arr) goto fail; + + spaxel_var_arr = (PyArrayObject*) PyArray_EMPTY(1, &npy_ncube, NPY_DOUBLE, 0); + if (!spaxel_var_arr) goto fail; + + spaxel_iflux_arr = (PyArrayObject*) PyArray_EMPTY(1, &npy_ncube, NPY_DOUBLE, 0); + if (!spaxel_iflux_arr) goto fail; + + spaxel_dq_arr = (PyArrayObject*) PyArray_EMPTY(1, &npy_ncube, NPY_INT, 0); + if (!spaxel_dq_arr) goto fail; + + result = Py_BuildValue("(OOOOO)", spaxel_flux_arr, spaxel_weight_arr, spaxel_var_arr, + spaxel_iflux_arr, spaxel_dq_arr); + + goto cleanup; + + } + + //______________________________________________________________________ + // if flag_dq_plane = 1, Set up the dq plane + //______________________________________________________________________ + int status1 = 0; + if(flag_dq_plane){ + if (instrument == 0){ + status1 = dq_miri(start_region, end_region,overlap_partial, overlap_full, + nxx,nyy,nwave, + cdelt1, cdelt2, roiw_ave, + (double *) PyArray_DATA(xc), + (double *) PyArray_DATA(yc), + (double *) PyArray_DATA(zc), + (double *) PyArray_DATA(coord1), + (double *) PyArray_DATA(coord2), + (double *) PyArray_DATA(wave), + (double *) PyArray_DATA(sliceno), + ncube, npt, + &spaxel_dq); + + } else{ + status1 = dq_nirspec(overlap_partial, + nxx,nyy,nwave, + cdelt1, cdelt2, roiw_ave, + (double *) PyArray_DATA(xc), + (double *) PyArray_DATA(yc), + (double *) PyArray_DATA(zc), + (double *) PyArray_DATA(coord1), + (double *) PyArray_DATA(coord2), + (double *) PyArray_DATA(wave), + (double *) PyArray_DATA(sliceno), + ncube, npt, + &spaxel_dq); + } + } else{ // set dq plane to 0 + status1 = set_dqplane_to_zero(ncube, &spaxel_dq); + + } + + //______________________________________________________________________ + // Match the point cloud elements to the spaxels they fail within the roi + //______________________________________________________________________ + if(weight_type ==0){ + status = match_point_emsm((double *) PyArray_DATA(xc), + (double *) PyArray_DATA(yc), + (double *) PyArray_DATA(zc), + (double *) PyArray_DATA(coord1), + (double *) PyArray_DATA(coord2), + (double *) PyArray_DATA(wave), + (double *) PyArray_DATA(flux), + (double *) PyArray_DATA(err), + (double *) PyArray_DATA(rois_pixel), + (double *) PyArray_DATA(roiw_pixel), + (double *) PyArray_DATA(scalerad_pixel), + (double *) PyArray_DATA(zcdelt3), + nxx, nyy, nwave, ncube, npt, cdelt1, cdelt2, + &spaxel_flux, &spaxel_weight, &spaxel_var, &spaxel_iflux); + } else{ + status = match_point_msm((double *) PyArray_DATA(xc), + (double *) PyArray_DATA(yc), + (double *) PyArray_DATA(zc), + (double *) PyArray_DATA(coord1), + (double *) PyArray_DATA(coord2), + (double *) PyArray_DATA(wave), + (double *) PyArray_DATA(flux), + (double *) PyArray_DATA(err), + (double *) PyArray_DATA(rois_pixel), + (double *) PyArray_DATA(roiw_pixel), + (double *) PyArray_DATA(weight_pixel), + (double *) PyArray_DATA(softrad_pixel), + (double *) PyArray_DATA(zcdelt3), + nxx, nyy, nwave, ncube, npt, cdelt1, cdelt2, + &spaxel_flux, &spaxel_weight, &spaxel_var, &spaxel_iflux); + } + + + if (status || status1) { + goto fail; + + } else { + // create return tuple: + npy_ncube = (npy_intp) ncube; + + spaxel_flux_arr = (PyArrayObject*) PyArray_SimpleNewFromData(1, &npy_ncube, NPY_DOUBLE, spaxel_flux); + if (!spaxel_flux_arr) goto fail; + spaxel_flux = NULL; + + spaxel_weight_arr = (PyArrayObject*) PyArray_SimpleNewFromData(1, &npy_ncube, NPY_DOUBLE, spaxel_weight); + if (!spaxel_weight_arr) goto fail; + spaxel_weight = NULL; + + spaxel_var_arr = (PyArrayObject*) PyArray_SimpleNewFromData(1, &npy_ncube, NPY_DOUBLE, spaxel_var); + if (!spaxel_var_arr) goto fail; + spaxel_var = NULL; + + spaxel_iflux_arr = (PyArrayObject*) PyArray_SimpleNewFromData(1, &npy_ncube, NPY_DOUBLE, spaxel_iflux); + if (!spaxel_iflux_arr) goto fail; + spaxel_iflux = NULL; + + spaxel_dq_arr = (PyArrayObject*) PyArray_SimpleNewFromData(1, &npy_ncube, NPY_INT, spaxel_dq); + if (!spaxel_dq_arr) goto fail; + spaxel_dq = NULL; + + PyArray_ENABLEFLAGS(spaxel_flux_arr, NPY_ARRAY_OWNDATA); + PyArray_ENABLEFLAGS(spaxel_weight_arr, NPY_ARRAY_OWNDATA); + PyArray_ENABLEFLAGS(spaxel_var_arr, NPY_ARRAY_OWNDATA); + PyArray_ENABLEFLAGS(spaxel_iflux_arr, NPY_ARRAY_OWNDATA); + PyArray_ENABLEFLAGS(spaxel_dq_arr, NPY_ARRAY_OWNDATA); + result = Py_BuildValue("(OOOOO)", spaxel_flux_arr, spaxel_weight_arr, spaxel_var_arr, + spaxel_iflux_arr, spaxel_dq_arr); + + goto cleanup; + + } + + fail: + Py_XDECREF(spaxel_flux_arr); + Py_XDECREF(spaxel_weight_arr); + Py_XDECREF(spaxel_var_arr); + Py_XDECREF(spaxel_iflux_arr); + Py_XDECREF(spaxel_dq_arr); + + free(spaxel_flux); + free(spaxel_weight); + free(spaxel_var); + free(spaxel_iflux); + free(spaxel_dq); + + if (!PyErr_Occurred()) { + PyErr_SetString(PyExc_MemoryError, + "Unable to allocate memory for output arrays."); + } + + cleanup: + if (free_xc)Py_XDECREF(xc); + if (free_yc) Py_XDECREF(yc); + if (free_zc) Py_XDECREF(zc); + if (free_coord1) Py_XDECREF(coord1); + if (free_coord2) Py_XDECREF(coord2); + if (free_wave) Py_XDECREF(wave); + if (free_flux) Py_XDECREF(flux); + if (free_err) Py_XDECREF(err); + if (free_rois_pixel) Py_XDECREF(rois_pixel); + if (free_roiw_pixel) Py_XDECREF(roiw_pixel); + if (free_scalerad_pixel) Py_XDECREF(scalerad_pixel); + if (free_softrad_pixel) Py_XDECREF(softrad_pixel); + if (free_weight_pixel) Py_XDECREF(weight_pixel); + if (free_zcdelt3) Py_XDECREF(zcdelt3); + if (free_sliceno) Py_XDECREF(sliceno); + + return result; +} + + +static PyMethodDef cube_methods[] = +{ + { + "cube_wrapper", + cube_wrapper, + METH_VARARGS, + "cube_wrapper(put in doc string)" + }, + {NULL, NULL, 0, NULL} /* sentinel */ +}; + + +static struct PyModuleDef moduledef = { + PyModuleDef_HEAD_INIT, + "cube_match_sky", /* m_name */ + "find point cloud matches for each spaxel center", /* m_doc */ + -1, /* m_size */ + cube_methods, /* m_methods */ + NULL, /* m_reload */ + NULL, /* m_traverse */ + NULL, /* m_clear */ + NULL, /* m_free */ +}; + +PyMODINIT_FUNC PyInit_cube_match_sky(void) +{ + PyObject* m; + import_array(); + m = PyModule_Create(&moduledef); + return m; +} diff --git a/jwst/cube_build/src/cube_utils.c b/jwst/cube_build/src/cube_utils.c new file mode 100644 index 0000000000..584643c8b2 --- /dev/null +++ b/jwst/cube_build/src/cube_utils.c @@ -0,0 +1,287 @@ +// This library contains c functions to support building IFU cubes + +#include +#include +#include +#include +#include + +#define CP_LEFT 0 +#define CP_RIGHT 1 +#define CP_BOTTOM 2 +#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) { + + const char *msg = "Couldn't allocate memory for output arrays."; + + // flux: + if (!(*fluxv = (double*)calloc(nelem, sizeof(double)))) { + PyErr_SetString(PyExc_MemoryError, msg); + goto failed_mem_alloc; + } + + //weight + if (!(*weightv = (double*)calloc(nelem, sizeof(double)))) { + PyErr_SetString(PyExc_MemoryError, msg); + goto failed_mem_alloc; + } + + //variance + if (!(*varv = (double*)calloc(nelem, sizeof(double)))) { + PyErr_SetString(PyExc_MemoryError, msg); + goto failed_mem_alloc; + } + + //iflux + if (!(*ifluxv = (double*)calloc(nelem, sizeof(double)))) { + PyErr_SetString(PyExc_MemoryError, msg); + goto failed_mem_alloc; + } + + return 0; + + failed_mem_alloc: + free(*fluxv); + free(*weightv); + free(*varv); + free(*ifluxv); + +} + +// support function for sh_find_overlap +void addpoint (double x, double y, double xnew[], double ynew[], int *nVertices2){ + 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){ + switch(edge) + { + case CP_LEFT: + return (x > left); + case CP_RIGHT: + return (x < right); + case CP_BOTTOM: + return (y > bottom); + case CP_TOP: + return (y < top); + } + 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) { + 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; + if(stat1 && stat2) return 2; + if(stat1 && !stat2) return 3; + if(!stat1 && !stat2) return 4; + return 0; //never executed + +} + +// support function for sh_find_overlap +void solveIntersection(int edge ,double x1,double y1,double x2,double y2, + double *x,double *y, + double left, double right, double top, double bottom){ + float m = 0; + if(x2 != x1) m = ((double)(y2-y1)/(double)(x2-x1)); + switch(edge) + { + case CP_LEFT: + *x = left; + *y = y1 + m * (*x - x1); + break; + case CP_RIGHT: + *x = right; + *y = y1 + m * (*x - x1); + break; + case CP_BOTTOM: + *y = bottom; + if(x1 != x2) + *x = x1 + (double)(1/m) * (*y - y1); + else + *x = x1; + break; + case CP_TOP: + *y = top; + if(x1 != x2) + *x = x1 + (double)(1/m) * (*y - y1); + else + *x = x1; + break; + } +} + + +double find_area_quad(double MinX, double MinY, double Xcorner[], double Ycorner[]){ + /* Find the area of an quadrilateral + + Parameters + ---------- + MinX : float + 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) + + Returns + ------- + Area + */ + + double PX[5]; + double PY[5]; + + PX[0] = Xcorner[0] - MinX; + PX[1] = Xcorner[1] - MinX; + PX[2] = Xcorner[2] - MinX; + PX[3] = Xcorner[3] - MinX; + PX[4]= PX[0]; + + PY[0] = Ycorner[0] - MinY; + PY[1] = Ycorner[1] - MinY; + PY[2] = Ycorner[2] - MinY; + PY[3] = Ycorner[3] - MinY; + PY[4] = PY[0]; + + double Area = 0.5 * ((PX[0] * PY[1] - PX[1] * PY[0]) + + (PX[1] * PY[2] - PX[2] * PY[1]) + + (PX[2] * PY[3] - PX[3] * PY[2]) + + (PX[3] * PY[4] - PX[4] * PY[3])); + + return fabs(Area); +} + + +// find the area of a closed polygon +double find_area_poly(int nVertices,double xPixel[],double yPixel[]){ + + double areaPoly = 0.0; + double xmin = xPixel[0]; + double ymin = yPixel[0]; + + for (int i = 1; i < nVertices; i++){ + if(xPixel[i] < xmin) xmin = xPixel[i]; + if(yPixel[i] < ymin) ymin = yPixel[i]; + } + + for (int i = 0; i < nVertices-1; i++){ + double area = ( xPixel[i]- xmin)*(yPixel[i+1]-ymin) - (xPixel[i+1]-xmin)*(yPixel[i]-ymin); + areaPoly = areaPoly + area; + } + areaPoly = 0.5* areaPoly; + return fabs(areaPoly); +} + + +//________________________________________________________________________________ +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 + + double areaClipped = 0.0; + double top = ycenter + 0.5*ylength; + double bottom = ycenter - 0.5*ylength; + double left = xcenter - 0.5*xlength; + double right = xcenter + 0.5*xlength; + + int nVertices = 4; // input detector pixel vertices + + int MaxVertices = 9; + double xPixel[9]= {0.0}; + double yPixel[9] = {0.0}; + double xnew[9]= {0.0}; + double ynew[9]= {0.0}; + + // initialize xPixel, yPixel to the detector pixel corners. + // xPixel,yPixel is become the clipped polygon vertices inside the cube pixel + for (int i = 0; i < 4; i++) { + xPixel[i] = xPixelCorner[i]; + yPixel[i] = yPixelCorner[i]; + } + xPixel[4] = xPixelCorner[0]; + yPixel[4] = yPixelCorner[0]; + + for (int i = 0 ; i < 4; i++) { // 0:left, 1: right, 2: bottom, 3: top + int nVertices2 = 0; + for (int j = 0; j< nVertices; j++){ + double x1 = xPixel[j]; + double y1 = yPixel[j]; + double x2 = xPixel[j+1]; + double y2 = yPixel[j+1]; + + int condition = calcCondition(i,x1,y1,x2,y2, + left,right,top,bottom); + + double x = 0; + double y = 0; + + switch(condition) + { + case 1: + solveIntersection(i,x1,y1,x2,y2, + &x, &y, + left,right,top,bottom); + + + addpoint (x, y, xnew, ynew, &nVertices2); + addpoint (x2, y2, xnew, ynew, &nVertices2); + break; + case 2: + addpoint (x2, y2, xnew, ynew, &nVertices2); + break; + case 3: + solveIntersection(i,x1,y1,x2,y2, + &x, &y, + left,right,top,bottom); + addpoint (x, y, xnew, ynew, &nVertices2); + break; + case 4: + break; + } + }// loop over j corners + + + addpoint (xnew[0], ynew[0], xnew, ynew, &nVertices2); // closed polygon + + //if(nVertices2 > MaxVertices ) { + // printf( " sh_find_overlap:: failure in finding the clipped polygon, nVertices2 > 9 \n "); + // exit(EXIT_FAILURE); + //} + + nVertices = nVertices2-1; + + for (int k = 0; k< nVertices2; k++){ + xPixel[k] = xnew[k]; + yPixel[k] = ynew[k]; + } + + //update + + } // loop over top,bottom,left,right + + nVertices++; + if(nVertices > 0) { + areaClipped = find_area_poly(nVertices,xPixel,yPixel); + } + + return areaClipped; +} diff --git a/jwst/cube_build/tests/test_configuration.py b/jwst/cube_build/tests/test_configuration.py index 8ddb1cf0eb..2bc9375041 100644 --- a/jwst/cube_build/tests/test_configuration.py +++ b/jwst/cube_build/tests/test_configuration.py @@ -213,7 +213,6 @@ def test_calspec2_config(_jail, miri_ifushort_short): output_type = 'multi' # calspec 2 setup. Only 1 cube create from 2 chanels single = False par_filename = 'None' - resol_filename = 'None' input_file = 'test.fits' input_models = [] @@ -234,7 +233,6 @@ def test_calspec2_config(_jail, miri_ifushort_short): input_models, input_filenames, par_filename, - resol_filename, **pars) master_table = file_table.FileTable() @@ -266,7 +264,6 @@ def test_calspec3_config_miri(_jail, miri_full_coverage): output_type = 'band' single = False par_filename = 'None' - resol_filename = 'None' input_file = 'test.fits' num_files = len(miri_full_coverage) @@ -287,7 +284,6 @@ def test_calspec3_config_miri(_jail, miri_full_coverage): miri_full_coverage, input_filenames, par_filename, - resol_filename, **pars) master_table = file_table.FileTable() @@ -348,7 +344,6 @@ def test_calspec3_config_miri_multi(_jail, miri_full_coverage): output_type = 'multi' single = False par_filename = 'None' - resol_filename = 'None' input_file = 'test.fits' num_files = len(miri_full_coverage) @@ -369,7 +364,6 @@ def test_calspec3_config_miri_multi(_jail, miri_full_coverage): miri_full_coverage, input_filenames, par_filename, - resol_filename, **pars) master_table = file_table.FileTable() @@ -411,7 +405,6 @@ def test_calspec3_config_nirspec(_jail, nirspec_medium_coverage): output_type = 'band' single = False par_filename = 'None' - resol_filename = 'None' input_file = 'test.fits' num_files = len(nirspec_medium_coverage) @@ -432,7 +425,6 @@ def test_calspec3_config_nirspec(_jail, nirspec_medium_coverage): nirspec_medium_coverage, input_filenames, par_filename, - resol_filename, **pars) master_table = file_table.FileTable() @@ -468,7 +460,6 @@ def test_calspec3_config_nirspec_multi(_jail, nirspec_medium_coverage): output_type = 'multi' single = False par_filename = 'None' - resol_filename = 'None' input_file = 'test.fits' num_files = len(nirspec_medium_coverage) @@ -489,7 +480,6 @@ def test_calspec3_config_nirspec_multi(_jail, nirspec_medium_coverage): nirspec_medium_coverage, input_filenames, par_filename, - resol_filename, **pars) master_table = file_table.FileTable() diff --git a/setup.py b/setup.py index bc5eeaa7fe..adfef8044b 100644 --- a/setup.py +++ b/setup.py @@ -34,6 +34,7 @@ # Include C extensions "jwst.lib.src": ["*.c"], + "jwst.cube_build.src": ["*.c"], # Include the transforms schemas "jwst.transforms": ["schemas/stsci.edu/jwst_pipeline/*.yaml"], @@ -58,6 +59,18 @@ ['jwst/lib/src/winclip.c'], include_dirs=include_dirs, define_macros=define_macros + ), + Extension( + 'jwst.cube_build.cube_match_internal', + ['jwst/cube_build/src/cube_match_internal.c','jwst/cube_build/src/cube_utils.c'], + include_dirs=include_dirs, + define_macros=define_macros + ), + Extension( + 'jwst.cube_build.cube_match_sky', + ['jwst/cube_build/src/cube_match_sky.c','jwst/cube_build/src/cube_utils.c'], + include_dirs=include_dirs, + define_macros=define_macros ) ], )