Compute spectral coverage in detail (fixes #212)

megaradrp/processing/extractobj.py

@@ -24,6 +24,38 @@
 from megaradrp.processing.fluxcalib import update_flux_limits
 
 
+def coverage_det(arr):
+    """Compute coverage"""
+
+    ma = numpy.max(arr)
+    mn = numpy.min(arr)
+
+    maxcover = len(arr) + 1
+
+    if ma == 0:
+        return slice(0, 0)
+
+    if mn == 1:
+        return slice(0, maxcover)
+
+    # We have 0s and 1s
+    diffa = numpy.diff(arr)
+
+    # For the left, find first +1
+    val_f, = numpy.where(diffa==1)
+    pix_f = 0
+    if len(val_f) > 0:
+        pix_f = val_f[0] + 1
+
+    # For the rigth, find first -1
+    val_l, = numpy.where(diffa==-1)
+    pix_l = maxcover
+    if len(val_l) > 0:
+        pix_l = val_l[0] + 1
+
+    return slice(pix_f, pix_l)
+
+
 def extract_star(rssimage, position, npoints, fiberconf, logger=None):
     """Extract a star given its center and the number of fibers to extract"""
 
@@ -77,12 +109,14 @@ def extract_star(rssimage, position, npoints, fiberconf, logger=None):
         valid_region = max_value_region
 
         # Interval with maximum coverage
-        nz_max, = numpy.nonzero(numpy.diff(max_value_region))
+        #nz_max, = numpy.nonzero(numpy.diff(max_value_region))
         # Interval with at least 1 fiber
-        nz_some, = numpy.nonzero(numpy.diff(some_value_region))
+        #nz_some, = numpy.nonzero(numpy.diff(some_value_region))
+        nz_max_slice = coverage_det(max_value_region)
+        nz_some_slice = coverage_det(some_value_region)
 
         # Collapse the flux in the optimal region
-        perf = flux_fiber[:, nz_max[0] + 1: nz_max[1] + 1].sum(axis=1)
+        perf = flux_fiber[:, nz_max_slice].sum(axis=1)
         # Contribution of each fiber to the total flux, 1D
         perf_norm = perf / perf.sum()
         contributions = numpy.zeros(shape=(rssdata.shape[0],))
@@ -113,7 +147,7 @@ def extract_star(rssimage, position, npoints, fiberconf, logger=None):
         # ax2.plot(coverage_total)
         #
         # plt.show()
-        pack = flux_total_c, colids, nz_max, nz_some
+        pack = flux_total_c, colids, nz_max_slice, nz_some_slice
         # pack = flux_total_c
         totals.append(pack)
 
@@ -267,7 +301,7 @@ def compute_dar(img, datamodel, logger=None, debug_plot=False):
 
 
 def generate_sensitivity(final, spectrum, star_interp, extinc_interp,
-                         cover1, cover2, sigma=20.0):
+                         wl_coverage1, wl_coverage2, sigma=20.0):
 
         wcsl = astropy.wcs.WCS(final[0].header)
 
@@ -292,32 +326,31 @@ def generate_sensitivity(final, spectrum, star_interp, extinc_interp,
         r0 = response_0 / r0max
         r1 = response_1 / r1max
 
-        pixm1, pixm2 = cover1
-        pixr1, pixr2 = cover2
+        pixm1, pixm2 = wl_coverage1
+        pixr1, pixr2 = wl_coverage2
 
         pixlims = {}
-        pixlims['PIXLIMR1'] = pixr1
+        pixlims['PIXLIMR1'] = pixr1 + 1  # Convert to 1-ref
         pixlims['PIXLIMR2'] = pixr2
-        pixlims['PIXLIMM1'] = pixm1
+        pixlims['PIXLIMM1'] = pixm1 + 1  # Convert to 1-ref
         pixlims['PIXLIMM2'] = pixm2
 
         max_valid = numpy.zeros_like(valid)
-        max_valid[pixm1:pixm2 + 1] = True
+        max_valid[wl_coverage1] = True
 
         partial_valid = numpy.zeros_like(valid)
-        partial_valid[pixr1:pixr2 + 1] = True
+        partial_valid[wl_coverage2] = True
 
-        valid = numpy.ones_like(response_0)
-        valid[pixm2:] = 0
-        valid[:pixm1+1] = 0
+        valid = numpy.zeros_like(response_0)
+        valid[wl_coverage1] = 1
 
         pixf1, pixf2 = int(math.floor(pixm1 +  2* sigma)), int(math.ceil(pixm2 - 2 * sigma))
 
-        pixlims['PIXLIMF1'] = pixf1
+        pixlims['PIXLIMF1'] = pixf1 + 1
         pixlims['PIXLIMF2'] = pixf2
 
         flux_valid = numpy.zeros_like(valid, dtype='bool')
-        flux_valid[pixf1:pixf2 + 1] = True
+        flux_valid[pixf1:pixf2] = True
 
         r0_ens = gaussian_filter(r0, sigma=sigma)
 

megaradrp/processing/tests/test_extractobj.py

@@ -0,0 +1,32 @@
+
+import megaradrp.processing.extractobj as eobj
+
+
+def test_cover_1():
+
+    arr = [1, 1, 1, 1, 1, 1]
+    res = eobj.coverage_det(arr)
+    assert res == slice(0, 7)
+
+    arr = [0, 0, 0, 0, 0, 0]
+    res = eobj.coverage_det(arr)
+    assert res == slice(0, 0)
+
+
+def test_cover_2():
+
+    arr = [0, 1, 1, 1, 0, 0]
+    res = eobj.coverage_det(arr)
+    assert res == slice(1, 4, None)
+
+    arr = [0, 1, 1, 1, 1, 0]
+    res = eobj.coverage_det(arr)
+    assert res == slice(1, 5)
+
+    arr = [0, 1, 1, 1, 1, 1]
+    res = eobj.coverage_det(arr)
+    assert res == slice(1, 7)
+
+    arr = [1, 1, 1, 1, 1, 0]
+    res = eobj.coverage_det(arr)
+    assert res == slice(0, 5)

megaradrp/recipes/calibration/lcbstdstar.py

@@ -104,7 +104,7 @@ def run(self, rinput):
         spectra_pack = extract_star(final, rinput.position, npoints,
                                     fiberconf, logger=self.logger)
 
-        spectrum, colids, cover1, cover2 = spectra_pack
+        spectrum, colids, wl_cover1, wl_cover2 = spectra_pack
         star_spectrum = fits.PrimaryHDU(spectrum, header=final[0].header)
 
         rad_vel = rinput.reference_spectrum_velocity * u.km / u.s
@@ -117,7 +117,7 @@ def run(self, rinput):
 
         fiber_ids = [colid + 1 for colid in colids]
         sigma = rinput.sigma_resolution
-        sens = generate_sensitivity(final, spectrum, star_interp, extinc_interp, cover1, cover2, sigma)
+        sens = generate_sensitivity(final, spectrum, star_interp, extinc_interp, wl_cover1, wl_cover2, sigma)
         self.logger.info('end LCBStandardRecipe reduction')
 
         return self.create_result(