Merge pull request #1876 from pypeit/hotfix_coadd3d

Fixes the WCS used for weighted combination

pypeit/core/datacube.py

@@ -238,8 +238,8 @@ def extract_point_source(wave, flxcube, ivarcube, bpmcube, wcscube, exptime,
 
     # Convert from counts/s/Ang/arcsec**2 to counts. The sensitivity function expects counts as input
     numxx, numyy, numwave = flxcube.shape
-    arcsecSQ = (wcscube.wcs.cdelt[0] * wcscube.wcs.cunit[0].to(units.arcsec)) * \
-               (wcscube.wcs.cdelt[1] * wcscube.wcs.cunit[1].to(units.arcsec))
+    arcsecSQ = abs((wcscube.wcs.cdelt[0] * wcscube.wcs.cunit[0].to(units.arcsec)) * \
+               (wcscube.wcs.cdelt[1] * wcscube.wcs.cunit[1].to(units.arcsec)))
     if fluxed:
         # The datacube is flux calibrated, in units of 10^-17 erg/s/cm**2/Ang/arcsec**2
         # Scale the flux and ivar cubes to be in units of erg/s/cm**2/Ang
@@ -1320,8 +1320,12 @@ def compute_weights(raImg, decImg, waveImg, sciImg, ivarImg, slitidImg,
     spec_bins = np.arange(1 + numwav) - 0.5
     bins = (xbins, ybins, spec_bins)
 
+    # Grab cos(dec) for convenience. Use the average of the min and max dec.
+    cosdec = np.cos(0.5 * (_dec_min + _dec_max) * np.pi / 180.0)
+    # Number of spaxels in the RA direction
+    numra = int((_ra_max - _ra_min) * cosdec / dspat)
+
     # Generate a 2D WCS to register all frames
-    numra = xbins.size - 1
     coord_min = [_ra_min, _dec_min, _wave_min]
     coord_dlt = [-dspat, dspat, dwv]
     whitelightWCS = generate_WCS(coord_min, coord_dlt, numra)
@@ -1357,7 +1361,7 @@ def compute_weights(raImg, decImg, waveImg, sciImg, ivarImg, slitidImg,
                                        sn_smooth_npix=sn_smooth_npix, weight_method=weight_method)
 
     # Because we pass back a weights array, we need to interpolate to assign each detector pixel a weight
-    all_wghts = [np.ones(_sciImg[0].shape) for _ in range(numframes)]
+    all_wghts = numframes*[np.ones(_sciImg[0].shape)]
     for ff in range(numframes):
         ww = (slitidImg[ff] > 0)
         all_wghts[ff][ww] = interp1d(wave_spec, weights[ff], kind='cubic',