diff --git a/jwst/cube_build/cube_build_wcs_util.py b/jwst/cube_build/cube_build_wcs_util.py index 987e53998f6..6bebd30a9f4 100644 --- a/jwst/cube_build/cube_build_wcs_util.py +++ b/jwst/cube_build/cube_build_wcs_util.py @@ -113,20 +113,11 @@ def find_corners_MIRI(input, this_channel, instrument_info, coord_system): if coord_system != 'internal_cal': # 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 + 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 # ***************************************************************************** diff --git a/jwst/cube_build/ifu_cube.py b/jwst/cube_build/ifu_cube.py index bf0a510ade1..1c95517f9c4 100644 --- a/jwst/cube_build/ifu_cube.py +++ b/jwst/cube_build/ifu_cube.py @@ -546,11 +546,6 @@ def build_ifucube(self): self.spaxel_iflux = np.zeros(total_num, dtype=np.float64) self.spaxel_dq = np.zeros(total_num, dtype=np.uint32) - for i in range(total_num): - self.spaxel_flux[i] = 0.0 - self.spaxel_iflux[i] = 0.0 - self.spaxel_weight[i] = 0.0 - self.spaxel_var[i] = 0.0 # ______________________________________________________________________________ subtract_background = True @@ -578,7 +573,6 @@ 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) @@ -586,18 +580,16 @@ def build_ifucube(self): coord1, coord2, wave, flux, err, slice_no, rois_pixel, roiw_pixel, weight_pixel,\ softrad_pixel, scalerad_pixel = pixelresult - # 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 + # by default flag the dq plane based on the FOV of the detector projected to sky flag_dq_plane = 1 - if wave.size == 0: - log.warning(f'No valid data found on file {input_model.meta.filename}') - else: - flag_dq_plane = 0 if self.skip_dqflagging: flag_dq_plane = 0 - t1 = time.time() - log.debug("Time to transform pixels to output frame = %.1f s" % (t1 - t0,)) + # check that there is valid data returned + # If all the data is flagged as DO_NOT_USE - not common- then log warning + if wave.size == 0: + log.warning(f'No valid data found on file {input_model.meta.filename}') + flag_dq_plane = 0 t0 = time.time() roiw_ave = np.mean(roiw_pixel) @@ -1391,10 +1383,10 @@ def map_detector_to_outputframe(self, this_par1, if self.linear_wavelength: min_wave_tolerance = self.zcoord[0] - self.roiw - max_wave_tolerance = (self.zcoord[-1] + self.roiw) + max_wave_tolerance = self.zcoord[-1] + self.roiw else: min_wave_tolerance = self.zcoord[0] - np.max(self.roiw_table) - max_wave_tolerance = (self.zcoord[-1] + np.max(self.roiw_table)) + max_wave_tolerance = self.zcoord[-1] + np.max(self.roiw_table) valid_min = np.where(wave_all >= min_wave_tolerance) not_mapped_low = wave_all.size - len(valid_min[0]) @@ -1414,8 +1406,8 @@ def map_detector_to_outputframe(self, this_par1, # all_flags = (dqflags.pixel['DO_NOT_USE'] + # dqflags.pixel['NON_SCIENCE']) - valid3 = np.logical_and((wave_all >= min_wave_tolerance), - (wave_all <= max_wave_tolerance)) + valid3 = np.logical_and(wave_all >= min_wave_tolerance, + wave_all <= max_wave_tolerance) # find the location of good data # good_data = np.where((np.bitwise_and(dq_all, all_flags) == 0) & @@ -1427,13 +1419,16 @@ def map_detector_to_outputframe(self, this_par1, # good data holds the location of pixels we want to map to cube # define variables as numpy arrays (numba needs this defined) - flux = np.zeros(flux_all[good_data].shape, dtype=np.float64) - err = np.zeros(flux_all[good_data].shape, dtype=np.float64) - coord1 = np.zeros(flux_all[good_data].shape, dtype=np.float64) - coord2 = np.zeros(flux_all[good_data].shape, dtype=np.float64) - wave = np.zeros(flux_all[good_data].shape, dtype=np.float64) - slice_no = np.zeros(flux_all[good_data].shape) - flux[:] = flux_all[good_data] + 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] diff --git a/jwst/cube_build/src/cube_match_internal.c b/jwst/cube_build/src/cube_match_internal.c index cd2563180cc..f9a7f2145f9 100644 --- a/jwst/cube_build/src/cube_match_internal.c +++ b/jwst/cube_build/src/cube_match_internal.c @@ -1,10 +1,13 @@ /* -The detector pixels are represented by a 'point could' 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 e modified shephard weighting method to determine how to weight each point clold member -in the spaxel. +Thie 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 area of the detector pixel overlapped onto an output +spaxel is used to weight the pixel flux for determining the spaxel flux. + Main function for Python: cube_wrapper_internal @@ -91,52 +94,8 @@ extern double sh_find_overlap(const double xcenter, const double ycenter, extern double find_area_quad(double MinX, double MinY, double Xcorner[], double Ycorner[]); -//_______________________________________________________________________ -// Allocate the memory to for the spaxel arrays in the cube -//_______________________________________________________________________ - -int mem_alloc(int nelem, double **fluxv, double **weightv, double **varv, double **ifluxv) { - - double *f, *w, *v, *i; - const char *msg = "Couldn't allocate memory for output arrays."; - - // flux: - f = (double*)malloc(nelem* sizeof(double)); - if (f) { - *fluxv = f; - } else { - PyErr_SetString(PyExc_MemoryError, msg); - return 1; - } - - // weight: - w = (double*)malloc(nelem* sizeof(double)); - if (w) { - *weightv = w; - } else { - PyErr_SetString(PyExc_MemoryError, msg); - return 1; - } - - // variance: - v = (double*)malloc(nelem* sizeof(double)); - if (v) { - *varv = v; - } else { - PyErr_SetString(PyExc_MemoryError, msg); - return 1; - } - // iflux - i = (double*)malloc(nelem* sizeof(double)); - if (i) { - *ifluxv = i; - } else { - PyErr_SetString(PyExc_MemoryError, msg); - return 1; - } +extern alloc_flux_arrays(int nelem, double **fluxv, double **weightv, double **varv, double **ifluxv); - return 0; -} //_______________________________________________________________________ // based on a 2-D drizzling type of algorithm find the contribution of @@ -151,20 +110,11 @@ int match_detector_cube(int instrument, int naxis1, int naxis2, int nz, int npt, double *pixel_flux, double *pixel_err, double **spaxel_flux, double **spaxel_weight, double **spaxel_var,double **spaxel_iflux){ - double *fluxv = NULL, *weightv=NULL, *varv=NULL ; // vectors for spaxel - double *ifluxv = NULL; // vector for spaxel + double *fluxv, *weightv, *varv, *ifluxv ; // vectors for spaxel // allocate memory to hold output - if (mem_alloc(ncube, &fluxv, &weightv, &varv, &ifluxv)) return 1; + if (alloc_flux_arrays(ncube, &fluxv, &weightv, &varv, &ifluxv)) return 1; - double set_zero=0.0; - // Set all data to zero - for (int i = 0; i < ncube; i++){ - varv[i] = set_zero; - fluxv[i] = set_zero; - ifluxv[i] = set_zero; - weightv[i] = set_zero; - } // loop over each valid point on detector and find match to IFU cube based // on along slice coordinate and wavelength diff --git a/jwst/cube_build/src/cube_match_sky.c b/jwst/cube_build/src/cube_match_sky.c index 2e931e4fa0a..27d1fcef58e 100644 --- a/jwst/cube_build/src/cube_match_sky.c +++ b/jwst/cube_build/src/cube_match_sky.c @@ -1,9 +1,9 @@ /* -The detector pixels are represented by a 'point could' on the sky. The IFU cube is +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) +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 @@ -103,69 +103,24 @@ extern double sh_find_overlap(const double xcenter, const double ycenter, const double xlength, const double ylength, double xPixelCorner[],double yPixelCorner[]); -// Allocate the memory to for the spaxel arrays in the cube +extern alloc_flux_arrays(int nelem, double **fluxv, double **weightv, double **varv, double **ifluxv); -int mem_alloc(int nelem, double **fluxv, double **weightv, double **varv, double **ifluxv) { - - double *f, *w, *v, *i; - const char *msg = "Couldn't allocate memory for output arrays."; - - // flux: - f = (double*)malloc(nelem* sizeof(double)); - if (f) { - *fluxv = f; - } else { - PyErr_SetString(PyExc_MemoryError, msg); - return 1; - } - - // weight: - w = (double*)malloc(nelem* sizeof(double)); - if (w) { - *weightv = w; - } else { - PyErr_SetString(PyExc_MemoryError, msg); - return 1; - } - - // variance: - v = (double*)malloc(nelem* sizeof(double)); - if (v) { - *varv = v; - } else { - PyErr_SetString(PyExc_MemoryError, msg); - return 1; - } - // iflux - i = (double*)malloc(nelem* sizeof(double)); - if (i) { - *ifluxv = i; - } else { - PyErr_SetString(PyExc_MemoryError, msg); - return 1; - } - - return 0; -} //________________________________________________________________________________ // allocate the memory for the spaxel DQ array -int mem_alloc_dq(long nelem, int **idqv) { - int *i; +int mem_alloc_dq(int nelem, int **idqv) { + const char *msg = "Couldn't allocate memory for output arrays."; - i = (int*)malloc(nelem * sizeof(int)); - if (i) { - *idqv = i; - } else { - PyErr_SetString(PyExc_MemoryError, msg); - return 1; + if (!(*idqv = (int*)calloc(nelem, sizeof(int)))) { + PyErr_SetString(PyExc_MemoryError, msg); + return 1; } + return 0; } - //________________________________________________________________________________ // Routine for MIRI DQ plane assignment @@ -438,10 +393,8 @@ int overlap_fov_with_spaxels(int overlap_partial, int overlap_full, } // end loop over ix return status ; - } - //________________________________________________________________________________ // Routine to setting NIRSpec dq plane for each wavelength plane @@ -511,9 +464,7 @@ int slice_wave_plane_nirspec(int w, int slicevalue, 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, @@ -573,15 +524,13 @@ int overlap_slice_with_spaxels(int overlap_partial, } // Swap start and end points if necessary and store swap state - //bool swapped; - //swapped = false; + if (x1 > x2){ x1 = x2; x2 = x1; y1 = y2; y2 = y1; - //swapped = true; } // Recalculate differences @@ -604,7 +553,7 @@ int overlap_slice_with_spaxels(int overlap_partial, yuse = x; xuse = y; } - //coord = (y, x) if is_steep else (x, y); + int index = (yuse * naxis1) + xuse; wave_slice_dq[index] = overlap_partial; error -= abs(dy); @@ -612,7 +561,6 @@ int overlap_slice_with_spaxels(int overlap_partial, y += ystep; error += dx; } - } return status; } @@ -621,17 +569,11 @@ int overlap_slice_with_spaxels(int overlap_partial, // Set the spaxel dq = 0. This is used when not determining the FOV on the sky for // setting the DQ plane. This is case for internalCal type cubes -int dq_set_zero(int ncube, int **spaxel_dq){ - int *idqv = NULL; // int vector for spaxel +int set_dqplane_to_zero(int ncube, int **spaxel_dq){ + int *idqv; // int vector for spaxel if (mem_alloc_dq(ncube, &idqv)) return 1; - - // Set all data to zero - for (long i = 0; i < ncube; i++){ - idqv[i] = 0; - } *spaxel_dq = idqv; - return 0; } @@ -647,66 +589,60 @@ int dq_miri(int start_region, int end_region, int overlap_partial, int overlap_f long ncube, int npt, int **spaxel_dq) { - int *idqv = NULL; // int vector for spaxel - - if (mem_alloc_dq(ncube, &idqv)) return 1; - - // Set all data to zero - for (long i = 0; i < ncube; i++){ - idqv[i] = 0; - } + 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]; + 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 + // 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 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 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]; + 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]; + 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]; - } + 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 + } + } // end loop over wavelength - *spaxel_dq = idqv; + *spaxel_dq = idqv; - return 0; + return 0; } //________________________________________________________________________________ @@ -745,14 +681,9 @@ int dq_nirspec(int overlap_partial, */ - int *idqv = NULL; // int vector for spaxel - + int *idqv ; // int vector for spaxel if (mem_alloc_dq(ncube, &idqv)) return 1; - - // Set all data to zero - for (long i = 0; i < ncube; i++){ - idqv[i] = 0; - } + // for each of the 30 slices - find the projection of this slice // onto each of the IFU wavelength planes. @@ -803,6 +734,8 @@ int dq_nirspec(int overlap_partial, // 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, @@ -814,23 +747,11 @@ int match_point_emsm(double *xc, double *yc, double *zc, double **spaxel_iflux) { - /* -return values: spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux -*/ - double *fluxv = NULL, *weightv=NULL, *varv=NULL ; // vectors for spaxel - double *ifluxv = NULL; // vector for spaxel + double *fluxv, *weightv, *varv, *ifluxv; // vector for spaxel - // allocate memory to hold output - if (mem_alloc(ncube, &fluxv, &weightv, &varv, &ifluxv)) return 1; - - double set_zero=0.0; - // Set all data to zero - for (int i = 0; i < ncube; i++){ - varv[i] = set_zero; - fluxv[i] = set_zero; - ifluxv[i] = set_zero; - weightv[i] = set_zero; - } + // 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 @@ -969,6 +890,7 @@ return values: spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux // 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, @@ -981,34 +903,22 @@ int match_point_msm(double *xc, double *yc, double *zc, double **spaxel_flux, double **spaxel_weight, double **spaxel_var, double **spaxel_iflux) { - /* -return values: spaxel_flux, spaxel_weight, spaxel_var, spaxel_iflux -*/ - double *fluxv = NULL, *weightv=NULL, *varv=NULL ; // vectors for spaxel - double *ifluxv = NULL; // vector for spaxel - - // allocate memory to hold output - if (mem_alloc(ncube, &fluxv, &weightv, &varv, &ifluxv)) return 1; - - double set_zero=0.0; - // Set all data to zero - for (int i = 0; i < ncube; i++){ - varv[i] = set_zero; - fluxv[i] = set_zero; - ifluxv[i] = set_zero; - weightv[i] = set_zero; - } - - // 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) { + 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){ @@ -1289,9 +1199,8 @@ static PyObject *cube_wrapper(PyObject *module, PyObject *args) { &spaxel_dq); } } else{ // set dq plane to 0 - - status1 = dq_set_zero(ncube, &spaxel_dq); - + status1 = set_dqplane_to_zero(ncube, &spaxel_dq); + } //______________________________________________________________________ @@ -1329,7 +1238,8 @@ static PyObject *cube_wrapper(PyObject *module, PyObject *args) { nxx, nyy, nwave, ncube, npt, cdelt1, cdelt2, &spaxel_flux, &spaxel_weight, &spaxel_var, &spaxel_iflux); } - + + if (status || status1) { goto fail; diff --git a/setup.cfg b/setup.cfg index 1dae02845de..fc49ca894ed 100644 --- a/setup.cfg +++ b/setup.cfg @@ -35,7 +35,6 @@ install_requires = gwcs>=0.16.1 jsonschema>=3.0.2 numpy>=1.17 - numba>=0.50.0 photutils>=1.1.0 psutil>=5.7.2 poppy>=0.9.0