Use centroid function in several places

megaradrp/datamodel.py

@@ -123,16 +123,20 @@ def gather_info_oresult(self, val):
         return [self.gather_info_dframe(f) for f in val.images]
 
     def fiber_scale_unit(self, img, unit=False):
-        funit = img['FIBERS'].header.get("FUNIT", "arcsec")
+        return fiber_scale_unit(img, unit=unit)
 
-        if funit == "arcsec":
-            scale = 1
-        else:
-            scale = self.PLATESCALE
-        if unit:
-            return scale, funit
-        else:
-            return scale
+
+def fiber_scale_unit(img, unit=False):
+    funit = img['FIBERS'].header.get("FUNIT", "arcsec")
+
+    if funit == "arcsec":
+        scale = 1
+    else:
+        scale = cons.GTC_FC_A_PLATESCALE.value
+    if unit:
+        return scale, funit
+    else:
+        return scale
 
 
 def get_fiberconf(img):
@@ -191,7 +195,6 @@ def read_fibers_extension(hdr, insmode='LCB'):
     return fp.FocalPlaneConf.from_header(hdr)
 
 
-
 def describe_hdulist_megara(hdulist):
     prim = hdulist[0].header
     instrument = prim.get("INSTRUME", "unknown")

megaradrp/ntypes.py

@@ -14,7 +14,7 @@
 from numina.core import DataFrameType
 from numina.types.product import DataProductMixin
 from numina.types.datatype import DataType
-from numina.types.array import ArrayType
+from numina.types.array import ArrayType, ArrayNType
 from numina.types.linescatalog import LinesCatalog
 from numina.exceptions import ValidationError
 
@@ -30,6 +30,12 @@ def validate_fiber_ext(header_f):
     _logger.debug('validate fiber extension')
 
 
+class Point2D(ArrayNType):
+    """A type of 2D cartesian coordinates"""
+    def __init__(self, *args, **kwds):
+        super(ArrayNType, self).__init__(2)
+
+
 class MegaraFrame(DataFrameType):
     """A processed frame"""
     DATATYPE = MegaraDataType.IMAGE_RAW

megaradrp/processing/centroid.py

@@ -0,0 +1,124 @@
+#
+# Copyright 2020 Universidad Complutense de Madrid
+#
+# This file is part of Megara DRP
+#
+# SPDX-License-Identifier: GPL-3.0+
+# License-Filename: LICENSE.txt
+#
+
+"""Compute centroids in MEGARA RSS images """
+
+import math
+import logging
+
+import numpy
+from scipy.spatial import KDTree
+from numina.constants import FWHM_G
+
+from megaradrp.instrument.focalplane import FocalPlaneConf
+
+
+_logger = logging.getLogger(__name__)
+
+
+def find_brightest_spaxel(final, extraction_region):
+    """Find the centroid around the brightest spaxel"""
+    fp_conf = FocalPlaneConf.from_img(final)
+    _logger.debug("LCB configuration is %s", fp_conf.conf_id)
+
+    rssdata = final[0].data
+    cut1, cut2 = extraction_region
+
+    flux_per_cell_all = rssdata[:, cut1:cut2].mean(axis=1)
+
+    max_cell = flux_per_cell_all.argmax() + 1
+    max_fiber_ = fp_conf.fibers[max_cell]
+
+    _logger.info("maximum flux in spaxel %d -- %s", max_cell, max_fiber_.name)
+    return max_fiber_
+
+
+def calc_centroid(final, extraction_region, point, nrings):
+    """compute the centroid around point"""
+    import megaradrp.instrument.constants as cons
+    import megaradrp.datamodel as dm
+
+    fp_conf = FocalPlaneConf.from_img(final)
+    _logger.debug("LCB configuration is %s", fp_conf.conf_id)
+
+    rssdata = final[0].data
+    cut1, cut2 = extraction_region
+
+    points = [point]
+
+    fibers = fp_conf.connected_fibers(valid_only=True)
+
+    grid_coords = []
+    for fiber in fibers:
+        grid_coords.append((fiber.x, fiber.y))
+    # setup kdtree for searching
+    kdtree = KDTree(grid_coords)
+
+    # Other possibility is
+    # query using radius instead
+    # radius = 1.2
+    # kdtree.query_ball_point(points, k=7, r=radius)
+    _logger.debug('adding %d nrings', nrings)
+    npoints = 1 + 3 * nrings * (nrings + 1)
+    _logger.debug('adding %d fibers', npoints)
+
+    dis_p, idx_p = kdtree.query(points, k=npoints)
+
+    _logger.info('Using %d nearest fibers', npoints)
+
+    scale, funit = dm.fiber_scale_unit(final, unit=True)
+    _logger.debug('unit is %s', funit)
+
+    platescale = cons.GTC_FC_A_PLATESCALE.value
+    positions = []
+    for diss, idxs, point in zip(dis_p, idx_p, points):
+        # For each point
+        value = [p * scale for p in point]
+        value_mm = [(v / platescale) for v in value]
+        _logger.info('For point %s arcsec', value)
+        _logger.info('For point %s mm', value_mm)
+        colids = []
+        coords = []
+        for dis, idx in zip(diss, idxs):
+            fiber = fibers[idx]
+            colids.append(fiber.fibid - 1)
+            coords.append((fiber.x, fiber.y))
+        _logger.debug('nearest fibers')
+        _logger.debug('%s', [col + 1 for col in colids])
+        coords = numpy.asarray(coords) * scale
+        # flux_per_cell = flux_per_cell_all[colids]
+        flux_per_cell = rssdata[colids, cut1:cut2].mean(axis=1)
+        flux_per_cell_total = flux_per_cell.sum()
+        flux_per_cell_norm = flux_per_cell / flux_per_cell_total
+        # centroid
+        scf = coords.T * flux_per_cell_norm
+        centroid = scf.sum(axis=1)
+        _logger.info('centroid: %s arcsec', list(centroid))
+        _logger.info('centroid: %s mm', list(centroid / platescale))
+        # central coords
+        c_coords = coords - centroid
+        scf0 = scf - centroid[:, numpy.newaxis] * flux_per_cell_norm
+        mc2 = numpy.dot(scf0, c_coords)
+        _logger.info('2nd order moments, x2=%f, y2=%f, xy=%f arcsec^2', mc2[0, 0], mc2[1, 1], mc2[0, 1])
+        if (mc2[0, 0] > 0) and (mc2[1, 1] > 0):
+            _logger.info('FWHM , x=%f, y=%f arcsec',
+                             FWHM_G * math.sqrt(mc2[0, 0]),
+                             FWHM_G * math.sqrt(mc2[1, 1])
+                             )
+        positions.append(centroid / platescale)
+    return positions[0]
+
+
+def calc_centroid_brightest(final, extraction_region, nrings):
+    """Compute the centroid around the brightest spaxel"""
+    max_fiber_ = find_brightest_spaxel(final, extraction_region)
+    point = (max_fiber_.x, max_fiber_.y)
+    # The brightest spaxel
+    position = calc_centroid(final, extraction_region, point, nrings)
+    return position

megaradrp/recipes/auxiliary/acquisitionlcb.py

@@ -20,8 +20,9 @@
 from numina.constants import FWHM_G
 
 from megaradrp.instrument.focalplane import FocalPlaneConf
-from megaradrp.recipes.scientific.base import ImageRecipe
 from megaradrp.ntypes import ProcessedRSS, ProcessedImage
+from megaradrp.processing.centroid import calc_centroid_brightest
+from megaradrp.recipes.scientific.base import ImageRecipe
 
 
 class AcquireLCBRecipe(ImageRecipe):
@@ -108,79 +109,7 @@ def run(self, rinput):
         fp_conf = FocalPlaneConf.from_img(final)
         self.logger.debug("LCB configuration is %s", fp_conf.conf_id)
 
-        rssdata = final[0].data
-
-        scale, funit = self.datamodel.fiber_scale_unit(final, unit=True)
-        self.logger.debug('unit is %s', funit)
-        platescale = self.datamodel.PLATESCALE
-
-        cut1, cut2 = rinput.extraction_region
-
-        # points = [(0, 0)] # Center of fiber 313
-        points = list(rinput.points)
-
-        flux_per_cell_all = rssdata[:, cut1:cut2].mean(axis=1)
-
-        max_cell = flux_per_cell_all.argmax() + 1
-        max_fiber_ = fp_conf.fibers[max_cell]
-
-        self.logger.info("maximum flux in spaxel %d -- %s", max_cell, max_fiber_.name)
-        # Extend points with the brightest spaxel
-        points.append((max_fiber_.x, max_fiber_.y))
-
-        fibers = fp_conf.connected_fibers(valid_only=True)
-
-        grid_coords = []
-        for fiber in fibers:
-            grid_coords.append((fiber.x, fiber.y))
-        # setup kdtree for searching
-        kdtree = KDTree(grid_coords)
-
-        # Other posibility is
-        # query using radius instead
-        # radius = 1.2
-        # kdtree.query_ball_point(points, k=7, r=radius)
-        self.logger.debug('adding %d nrings', rinput.nrings)
-        npoints = 1 + 3 * rinput.nrings * (rinput.nrings + 1)
-        self.logger.debug('adding %d fibers', npoints)
-
-        dis_p, idx_p = kdtree.query(points, k=npoints)
-
-        self.logger.info('Using %d nearest fibers', npoints)
-        for diss, idxs, point in zip(dis_p, idx_p, points):
-            # For each point
-            value = [p * scale for p in point]
-            value_mm = [(v / platescale) for v in value]
-            self.logger.info('For point %s arcsec', value)
-            self.logger.info('For point %s mm', value_mm)
-            colids = []
-            coords = []
-            for dis, idx in zip(diss, idxs):
-                fiber = fibers[idx]
-                colids.append(fiber.fibid - 1)
-                coords.append((fiber.x, fiber.y))
-            self.logger.debug('nearest fibers')
-            self.logger.debug('%s', [col + 1 for col in colids])
-            coords = np.asarray(coords) * scale
-            # flux_per_cell = flux_per_cell_all[colids]
-            flux_per_cell = rssdata[colids, cut1:cut2].mean(axis=1)
-            flux_per_cell_total = flux_per_cell.sum()
-            flux_per_cell_norm = flux_per_cell / flux_per_cell_total
-            # centroid
-            scf = coords.T * flux_per_cell_norm
-            centroid = scf.sum(axis=1)
-            self.logger.info('centroid: %s arcsec', list(centroid))
-            self.logger.info('centroid: %s mm', list(centroid / platescale))
-            # central coords
-            c_coords = coords - centroid
-            scf0 = scf - centroid[:, np.newaxis] * flux_per_cell_norm
-            mc2 = np.dot(scf0, c_coords)
-            self.logger.info('2nd order moments, x2=%f, y2=%f, xy=%f arcsec^2', mc2[0,0], mc2[1,1], mc2[0,1])
-            if (mc2[0,0] > 0) and (mc2[1,1] > 0):
-                self.logger.info('FWHM , x=%f, y=%f arcsec',
-                                 FWHM_G * math.sqrt(mc2[0,0]),
-                                 FWHM_G * math.sqrt(mc2[1,1])
-                                 )
+        centroid = calc_centroid_brightest(final, rinput.extraction_region, rinput.nrings)
 
         if False:
             self.compute_dar(final)

megaradrp/recipes/calibration/lcbstdstar.py

@@ -14,7 +14,6 @@
 import astropy.io.fits as fits
 import astropy.units as u
 import astropy.wcs
-
 from scipy.interpolate import interp1d
 
 from numina.array.numsplines import AdaptiveLSQUnivariateSpline
@@ -25,12 +24,14 @@
 from numina.types.array import ArrayType
 
 from megaradrp.instrument.focalplane import FocalPlaneConf
-from megaradrp.processing.extractobj import extract_star, generate_sensitivity
-from megaradrp.processing.extractobj import mix_values, compute_broadening
-from megaradrp.recipes.scientific.base import ImageRecipe
+from megaradrp.ntypes import Point2D
 from megaradrp.ntypes import ProcessedRSS, ProcessedFrame, ProcessedSpectrum
 from megaradrp.ntypes import ReferenceSpectrumTable, ReferenceExtinctionTable
 from megaradrp.ntypes import MasterSensitivity
+from megaradrp.processing.extractobj import extract_star, generate_sensitivity
+from megaradrp.processing.extractobj import mix_values, compute_broadening
+from megaradrp.processing.centroid import calc_centroid_brightest
+from megaradrp.recipes.scientific.base import ImageRecipe
 
 
 class LCBStandardRecipe(ImageRecipe):
@@ -71,7 +72,7 @@ class LCBStandardRecipe(ImageRecipe):
     the central spaxel containing the star and returned as `star_spectrum`.
 
     """
-    position = Requirement(list, "Position of the reference object", default=(0, 0))
+    position = Requirement(Point2D, "Position of the reference object", optional=True)
     nrings = Parameter(3, "Number of rings to extract the star",
                        validator=range_validator(minval=1))
     reference_spectrum = Requirement(ReferenceSpectrumTable, "Spectrum of reference star")
@@ -140,12 +141,24 @@ def run(self, rinput):
         # 1 + 6  + 12  + 18 for third ring
         # 1 + 6 * Sum_i=0^n =  1 + 3 * n * (n +1)
         # Using three rings around central point
+
+        # If position is None, find the brightest spaxel
+        # and use the centroid
+        if rinput.position is None:
+            self.logger.info('finding centroid of brightest spaxel')
+            extraction_region = [1000, 3000]
+            nrings = rinput.nrings
+            position = calc_centroid_brightest(final, extraction_region, nrings)
+        else:
+            position = rinput.position
+        self.logger.info('central position is %s', position)
+
         self.logger.debug('adding %d nrings', rinput.nrings)
         npoints = 1 + 3 * rinput.nrings * (rinput.nrings +1)
         self.logger.debug('adding %d fibers', npoints)
 
         fp_conf = FocalPlaneConf.from_img(final)
-        spectra_pack = extract_star(final, rinput.position, npoints,
+        spectra_pack = extract_star(final, position, npoints,
                                     fp_conf, logger=self.logger)
 
         spectrum, colids, wl_cover1, wl_cover2 = spectra_pack