From 670686a76b100ee8f338156e3359f95498d92acb Mon Sep 17 00:00:00 2001 From: Sergio Pascual Date: Tue, 14 Apr 2020 17:54:29 +0200 Subject: [PATCH] Subtract SkyRss from image --- megaradrp/datamodel.py | 10 ++- megaradrp/processing/extractobj.py | 45 +------------ megaradrp/processing/sky.py | 88 ++++++++++++++++++++++++++ megaradrp/processing/tests/test_sky.py | 32 ++++++++++ megaradrp/recipes/scientific/base.py | 8 +-- megaradrp/tests/test_requirements.py | 8 +++ 6 files changed, 140 insertions(+), 51 deletions(-) create mode 100644 megaradrp/processing/sky.py create mode 100644 megaradrp/processing/tests/test_sky.py create mode 100644 megaradrp/tests/test_requirements.py diff --git a/megaradrp/datamodel.py b/megaradrp/datamodel.py index e729c5d4..3eea9c6a 100644 --- a/megaradrp/datamodel.py +++ b/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) diff --git a/megaradrp/processing/extractobj.py b/megaradrp/processing/extractobj.py index 74ad6950..aaa727a0 100644 --- a/megaradrp/processing/extractobj.py +++ b/megaradrp/processing/extractobj.py @@ -11,7 +11,6 @@ import math import uuid -import logging import numpy import astropy.wcs @@ -19,7 +18,6 @@ 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 diff --git a/megaradrp/processing/sky.py b/megaradrp/processing/sky.py new file mode 100644 index 00000000..274ff739 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/megaradrp/processing/tests/test_sky.py b/megaradrp/processing/tests/test_sky.py new file mode 100644 index 00000000..269db7ed --- /dev/null +++ b/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 diff --git a/megaradrp/recipes/scientific/base.py b/megaradrp/recipes/scientific/base.py index b7e3d10d..af351127 100644 --- a/megaradrp/recipes/scientific/base.py +++ b/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 ) \ No newline at end of file diff --git a/megaradrp/tests/test_requirements.py b/megaradrp/tests/test_requirements.py new file mode 100644 index 00000000..9b16f74e --- /dev/null +++ b/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) \ No newline at end of file