WaveLength Calibration has been modified to generate a JSON file to a…

…ddress more information than just the coeficient values of each aperture

megaradrp/core/recipe.py

@@ -19,6 +19,7 @@
 
 import logging
 
+import numpy as np
 from astropy.io import fits
 import numina.array.combine as combine
 from numina.flow import SerialFlow
@@ -191,3 +192,8 @@ def get_parameters(self, rinput):
         #     pass
 
         return parameters
+
+    def get_wlcalib(self, data):
+
+        wlcalib = [elem['aperture']['function']['coecifients'] for elem in data]
+        return np.array(wlcalib)

megaradrp/products.py

@@ -102,10 +102,6 @@ def _datatype_load(self, obj):
         return traces
 
 
-class WavelengthCalibration(ArrayType, DataProductTag):
-    pass
-
-
 class MasterWeights(DataProductType):
     def __init__(self, default=None):
         super(MasterWeights, self).__init__(ptype=dict, default=default)
@@ -147,4 +143,7 @@ def _datatype_load(self, obj):
                 data = json.load(fd)
         except IOError as e:
             raise e
-        return data
+        return data
+
+class WavelengthCalibration(JSONstorage):
+    pass

megaradrp/recipes/auxiliary/focusspec.py

@@ -110,7 +110,7 @@ def run(self, rinput):
         _logger.info('pair lines in images')
         line_fibers = self.filter_lines(ever)
 
-        focus_wavelength = self.generateJSON(ever, rinput.wlcalib,
+        focus_wavelength = self.generateJSON(ever, self.get_wlcalib(rinput.wlcalib),
                                              rinput.obresult.images)
 
         _logger.info('fit FWHM of lines')

megaradrp/recipes/calibration/arc.py

@@ -26,23 +26,19 @@
 import numpy
 from astropy.io import fits
 from scipy.interpolate import interp1d
-
 from numina.core import Requirement, Product, Parameter
 from numina.core import DataFrameType
-from numina.core.products import ArrayType
 from numina.core.requirements import ObservationResultRequirement
 from numina.core.products import LinesCatalog
 from numina.array.wavecal.arccalibration import arccalibration_direct
 from numina.array.wavecal.arccalibration import fit_solution
 from numina.array.wavecal.arccalibration import gen_triplets_master
 from numina.array.wavecal.statsummary import sigmaG
 # FIXME: still using findpeaks1D functions
-from numina.array.peaks.findpeaks1D import findPeaks_spectrum
-from numina.array.peaks.findpeaks1D import refinePeaks_spectrum
 from numina.array.peaks.peakdet import find_peaks_indexes, refine_peaks
 
 from megaradrp.core.recipe import MegaraBaseRecipe
-from megaradrp.products import TraceMap, WavelengthCalibration
+from megaradrp.products import TraceMap, WavelengthCalibration, JSONstorage
 from megaradrp.requirements import MasterBiasRequirement, MasterBPMRequirement
 from megaradrp.requirements import MasterDarkRequirement
 from megaradrp.core.processing import apextract_tracemap
@@ -65,6 +61,7 @@ class ArcCalibrationRecipe(MegaraBaseRecipe):
     arc_image = Product(DataFrameType)
     arc_rss = Product(DataFrameType)
     master_wlcalib = Product(WavelengthCalibration)
+    data_wlcalib = Product(JSONstorage)
 
     def __init__(self):
         super(ArcCalibrationRecipe, self).__init__("0.1.0")
@@ -83,12 +80,12 @@ def run(self, rinput):
         _logger.info('extracted %i fibers', rssdata.shape[0])
 
         # Skip any other inputs for the moment
-        coeff_table = self.calibrate_wl(rssdata, rinput.lines_catalog,
+        coeff_table, data_wlcalib = self.calibrate_wl(rssdata, rinput.lines_catalog,
                                         rinput.polynomial_degree)
 
         # WL calibration goes here
         return self.create_result(arc_image=reduced, arc_rss=rss,
-                                  master_wlcalib=coeff_table)
+                                  master_wlcalib=coeff_table, data_wlcalib = data_wlcalib)
 
     def pintar_todas(self, diferencia_final, figure):
         import matplotlib.pyplot as plt
@@ -106,6 +103,9 @@ def calibrate_wl(self, rss, lines_catalog, poldeg, times_sigma=50.0):
         #
         # read master table (TBM) and generate auxiliary parameters (valid for
         # all the slits) for the wavelength calibration
+        lista_solution = []
+        lista_xpeaks_refined = []
+
         wv_master = lines_catalog[:, 0]
         ntriplets_master, ratios_master_sorted, triplets_master_sorted_list = \
             gen_triplets_master(wv_master)
@@ -119,11 +119,11 @@ def calibrate_wl(self, rss, lines_catalog, poldeg, times_sigma=50.0):
             # find peaks (initial search providing integer numbers)
             threshold = numpy.median(row) + times_sigma * sigmaG(row)
 
-            ipeaks_int = findPeaks_spectrum(row, nwinwidth, threshold)
-            ipeaks_float = refinePeaks_spectrum(row, ipeaks_int, nwinwidth, method=1)
+            # ipeaks_int = findPeaks_spectrum(row, nwinwidth, threshold)
+            # ipeaks_float = refinePeaks_spectrum(row, ipeaks_int, nwinwidth, method=1)
 
-            # ipeaks_int = find_peaks_indexes(row, nwinwidth, threshold)
-            # ipeaks_float = refine_peaks(row, ipeaks_int, nwinwidth)[0]
+            ipeaks_int = find_peaks_indexes(row, nwinwidth, threshold)
+            ipeaks_float = refine_peaks(row, ipeaks_int, nwinwidth)[0]
 
             # self.pintar_todas(ipeaks_float,'refinado_nico')
             # self.pintar_todas(ipeaks_float2,'refinado_nuevo')
@@ -142,7 +142,7 @@ def calibrate_wl(self, rss, lines_catalog, poldeg, times_sigma=50.0):
             finterp_channel = interp1d(range(xchannel.size), xchannel,
                                        kind='linear')
             xpeaks_refined = finterp_channel(ipeaks_float)
-
+            lista_xpeaks_refined.append(xpeaks_refined)
             try:
                 solution = arccalibration_direct(wv_master,
                                                  ntriplets_master,
@@ -174,10 +174,67 @@ def calibrate_wl(self, rss, lines_catalog, poldeg, times_sigma=50.0):
                              cdelt1_approx)
 
                 _logger.info('fitted coefficients %s', numpy_array_with_coeff)
-
+                lista_solution.append(solution)
                 coeff_table[idx] = numpy_array_with_coeff
 
             except TypeError as error:
                 _logger.error("%s", error)
 
-        return coeff_table
+        data_wlcalib = self.generateJSON(coeff_table, lista_solution,
+                          lista_xpeaks_refined, poldeg, lines_catalog)
+
+        _logger.info('Fin')
+
+        return coeff_table, data_wlcalib
+
+    def generateJSON(self, coeff_table, lista_solution,
+                     lista_xpeaks_refined, poldeg, lines_catalog):
+        '''
+            Final format of the features field is:{
+                  "features": [[<x_position>,
+                                calculated_wavelength,
+                                <original_lambda>,
+                                <original_flux>
+                                ]]
+        '''
+
+        from numpy.polynomial.polynomial import polyval
+        _logger.info('start JSON generation')
+
+        result = []
+        counter = 0
+        for ind, xpeaks in enumerate(lista_xpeaks_refined):
+            res = polyval(xpeaks, coeff_table[ind])
+            _logger.info('indice: %s', ind)
+            record = {}
+
+            features = []
+            feature = {'xpos':None,
+                       'wavelength':None,
+                       'reference':None,
+                       'flux':None,
+                       'category':None
+              }
+            for aux, elem in enumerate(xpeaks):
+                feature['xpos'] = xpeaks[aux]
+                feature['wavelength'] = res[aux]
+                feature['reference'] = lines_catalog[aux][0]
+                feature['flux'] = lines_catalog[aux][1]
+                feature['category'] = lista_solution[ind][aux]['type']
+                features.append(feature)
+
+            function = {
+                'method':'least squares',
+                'order':poldeg,
+                'coecifients': coeff_table[ind].tolist()
+            }
+
+            record['aperture'] = {'id': ind + 1,
+                                  'features': features,
+                                  'function': function}
+            result.append(record)
+            counter += 1
+
+        _logger.info('end JSON generation')
+
+        return result

megaradrp/recipes/calibration/flat.py

@@ -73,7 +73,7 @@ def run(self, rinput):
         rss = numpy.array(self.load_files_from_directory(rinput.master_weights, img=reduced[0].data))
 
         _logger.info('Starting: resample_rss_flux')
-        final, wcsdata = self.resample_rss_flux(rss.T, rinput.wlcalib)
+        final, wcsdata = self.resample_rss_flux(rss.T, self.get_wlcalib(rinput.wlcalib))
 
         _logger.info('Compute mean and resampling again')
 

megaradrp/recipes/calibration/twilight.py

@@ -63,7 +63,7 @@ def run(self, rinput):
 
         return self.run_args(rinput.obresult,
                              rinput.master_weights,
-                             rinput.wlcalib,
+                             self.get_wlcalib(rinput.wlcalib),
                              parameters
                             )
 
@@ -145,7 +145,7 @@ def run(self, rinput):
         self.logger.info('extraction completed')
 
         _logger.info('resampling spectra')
-        final, wcsdata = resample_rss_flux(rsshdu.data, rinput.wlcalib)
+        final, wcsdata = resample_rss_flux(rsshdu.data, self.get_wlcalib(rinput.wlcalib))
         # This value was ~0.4% and now is 4e-6 %
         # (abs(final.sum()-hdu_t.data.sum())/hdu_t.data.sum()*100)