Subtract SkyRss from image

megaradrp/datamodel.py

@@ -313,8 +313,8 @@ def get_fiberconf(img):
         return get_fiberconf_default(main_insmode)
 
 
-def get_fiberconf_default(insmode):
-    """Obtain default FiberConf object"""
+def create_default_fiber_header(insmode):
+    """Obtain default FIBER header"""
     if insmode == 'LCB':
         slit_file = 'lcb_default_header.txt'
     elif insmode == 'MOS':
@@ -326,6 +326,12 @@ def get_fiberconf_default(insmode):
     data = pkgutil.get_data('megaradrp.instrument.configs', slit_file)
     default_hdr = StringIO(data.decode('utf8'))
     hdr_fiber = fits.header.Header.fromfile(default_hdr)
+    return hdr_fiber
+
+
+def get_fiberconf_default(insmode):
+    """Obtain default FiberConf object"""
+    hdr_fiber = create_default_fiber_header(insmode)
     return read_fibers_extension(hdr_fiber)
 
 

megaradrp/processing/extractobj.py

@@ -11,15 +11,13 @@
 
 import math
 import uuid
-import logging
 
 import numpy
 import astropy.wcs
 import astropy.io.fits as fits
 import astropy.units as u
 from scipy.spatial import KDTree
 from scipy.ndimage.filters import gaussian_filter
-from numina.frame.utils import copy_img
 from numina.array.wavecalib.crosscorrelation import periodic_corr1d
 
 import megaradrp.datamodel as dm
@@ -443,6 +441,7 @@ def generate_sensitivity(final, spectrum, star_interp, extinc_interp,
         # FIXME: add history
         sens = fits.PrimaryHDU(s_response, header=final[0].header)
         # delete second axis keywords
+        # FIXME: delete axis with wcslib
         for key in ['CRPIX2', 'CRVAL2', 'CDELT2', 'CTYPE2']:
             if key in sens.header:
                 del sens.header[key]
@@ -455,45 +454,3 @@ def generate_sensitivity(final, spectrum, star_interp, extinc_interp,
         return sens
 
 
-def subtract_sky(img, ignored_sky_bundles=None, logger=None):
-    # Sky subtraction
-
-    if logger is None:
-        logger = logging.getLogger(__name__)
-
-    logger.info('obtain fiber information')
-    sky_img = copy_img(img)
-    final_img = copy_img(img)
-    fiberconf = dm.get_fiberconf(sky_img)
-    # Sky fibers
-    skyfibs = fiberconf.sky_fibers(valid_only=True,
-                                   ignored_bundles=ignored_sky_bundles)
-    logger.debug('sky fibers are: %s', skyfibs)
-    # Create empty sky_data
-    target_data = img[0].data
-
-    target_map = img['WLMAP'].data
-    sky_data = numpy.zeros_like(img[0].data)
-    sky_map = numpy.zeros_like(img['WLMAP'].data)
-    sky_img[0].data = sky_data
-
-    for fibid in skyfibs:
-        rowid = fibid - 1
-        sky_data[rowid] = target_data[rowid]
-        sky_map[rowid] = target_map[rowid]
-    # Sum
-    coldata = sky_data.sum(axis=0)
-    colsum = sky_map.sum(axis=0)
-
-    # Divide only where map is > 0
-    mask = colsum > 0
-    avg_sky = numpy.zeros_like(coldata)
-    avg_sky[mask] = coldata[mask] / colsum[mask]
-
-    # This should be done only on valid fibers
-    logger.info('ignoring invalid fibers: %s', fiberconf.invalid_fibers())
-    for fibid in fiberconf.valid_fibers():
-        rowid = fibid - 1
-        final_img[0].data[rowid, mask] = img[0].data[rowid, mask] - avg_sky[mask]
-
-    return final_img, img, sky_img

megaradrp/processing/sky.py

@@ -0,0 +1,88 @@
+#
+# Copyright 2019-2020 Universidad Complutense de Madrid
+#
+# This file is part of Megara DRP
+#
+# SPDX-License-Identifier: GPL-3.0+
+# License-Filename: LICENSE.txt
+#
+
+import logging
+
+import numpy
+
+import megaradrp.datamodel as dm
+from numina.frame.utils import copy_img
+
+
+def subtract_sky(img, ignored_sky_bundles=None, logger=None):
+    # Sky subtraction
+
+    if logger is None:
+        logger = logging.getLogger(__name__)
+
+    logger.info('obtain fiber information')
+    sky_img = copy_img(img)
+    final_img = copy_img(img)
+    fiberconf = dm.get_fiberconf(sky_img)
+    # Sky fibers
+    skyfibs = fiberconf.sky_fibers(valid_only=True,
+                                   ignored_bundles=ignored_sky_bundles)
+    logger.debug('sky fibers are: %s', skyfibs)
+    # Create empty sky_data
+    target_data = img[0].data
+
+    target_map = img['WLMAP'].data
+    sky_data = numpy.zeros_like(img[0].data)
+    sky_map = numpy.zeros_like(img['WLMAP'].data)
+    sky_img[0].data = sky_data
+
+    for fibid in skyfibs:
+        rowid = fibid - 1
+        sky_data[rowid] = target_data[rowid]
+        sky_map[rowid] = target_map[rowid]
+    # Sum
+    coldata = sky_data.sum(axis=0)
+    colsum = sky_map.sum(axis=0)
+
+    # Divide only where map is > 0
+    mask = colsum > 0
+    avg_sky = numpy.zeros_like(coldata)
+    avg_sky[mask] = coldata[mask] / colsum[mask]
+
+    # This should be done only on valid fibers
+    logger.info('ignoring invalid fibers: %s', fiberconf.invalid_fibers())
+    for fibid in fiberconf.valid_fibers():
+        rowid = fibid - 1
+        final_img[0].data[rowid, mask] = img[0].data[rowid, mask] - avg_sky[mask]
+    # Update headers
+    #
+    return final_img, img, sky_img
+
+
+def subtract_sky_rss(img, sky_img, ignored_sky_bundles=None, logger=None):
+    """Subtract a sky image from an image"""
+    # Sky subtraction
+
+    if logger is None:
+        logger = logging.getLogger(__name__)
+
+    #logger.info('obtain fiber information')
+    final_img = copy_img(img)
+    # fiberconf_sky = dm.get_fiberconf(sky_img)
+    # fiberconf_target = dm.get_fiberconf(img)
+
+    logger.debug('using WLMAP extension to compute valid regions')
+
+    v_map = img['WLMAP'].data > 0
+    sky_map = numpy.zeros_like(img['WLMAP'].data)
+    sky_data = sky_img[0].data
+    sky_map[:] = v_map[:]
+
+    # This should be done only on valid fibers
+    #logger.info('ignoring invalid fibers: %s', fiberconf_target.invalid_fibers())
+    final_img[0].data[v_map] = img[0].data[v_map] - sky_data[v_map]
+    final_img[0].data[~v_map] = 0.0
+    # Update headers
+    #
+    return final_img, img, sky_img

megaradrp/processing/tests/test_sky.py

@@ -0,0 +1,32 @@
+import numpy
+
+from megaradrp.datamodel import create_default_fiber_header
+from ..sky import subtract_sky_rss
+
+
+def create_rss(value, wlmap):
+    import astropy.io.fits as fits
+    data1 = value + numpy.zeros((623, 4300), dtype='float32')
+    hdu = fits.PrimaryHDU(data1)
+    hdrf = create_default_fiber_header('LCB')
+    fibers = fits.ImageHDU(header=hdrf, name='FIBERS')
+    rss_map = fits.ImageHDU(wlmap, name='WLMAP')
+    return fits.HDUList([hdu, fibers, rss_map])
+
+
+def test_subtract_sky_rss():
+
+    wlmap = numpy.zeros((623, 4300), dtype='float32')
+    wlmap[:,350:4105] = 1.0
+    wlmap[622,:] = 0
+    img1 = create_rss(1000, wlmap)
+    img2 = create_rss(400, wlmap)
+
+    final_img, img, sky_img = subtract_sky_rss(img1, img2)
+    assert img is img1
+    # In final image, regions outside WLMAP must be at zero
+    assert final_img[0].data[:, 100:200].min() == 0
+    assert final_img[0].data[:, 100:200].max() == 0
+
+    assert final_img[0].data[622, :].max() == 0
+    assert final_img[0].data[622, :].min() == 0

megaradrp/recipes/scientific/base.py

@@ -25,7 +25,7 @@
 from megaradrp.processing.fiberflat import FlipLR, FiberFlatCorrector
 from megaradrp.processing.twilight import TwilightCorrector
 from megaradrp.processing.extractobj import compute_centroid, compute_dar
-from megaradrp.processing.extractobj import subtract_sky
+from megaradrp.processing.sky import subtract_sky, subtract_sky_rss
 
 
 class ImageRecipe(MegaraBaseRecipe):
@@ -100,7 +100,7 @@ def run_reduction_1d(self, img, tracemap, wlcalib, fiberflat, twflat=None, offse
 
         flow2 = SerialFlow(correctors)
 
-        reduced_rss =  flow2(img)
+        reduced_rss = flow2(img)
         return reduced_rss
 
     def compute_dar(self, img):
@@ -125,8 +125,6 @@ def run_sky_subtraction(self, img, sky_rss=None, ignored_sky_bundles=None):
                                 )
         else:
             self.logger.info('use sky RSS image')
-            return subtract_sky(img,
-                                self.datamodel,
-                                ignored_sky_bundles=ignored_sky_bundles,
+            return subtract_sky_rss(img, sky_img=sky_rss,
                                 logger=self.logger
                                 )

megaradrp/tests/test_requirements.py

@@ -0,0 +1,8 @@
+
+from numina.core import Requirement
+from ..requirements import DiffuseLightRequirement
+
+
+def test_requires_df_l():
+    req = DiffuseLightRequirement()
+    assert isinstance(req, Requirement)