Merge 1da3a89 into bf83a7f

soar-telescope · Oct 9, 2018 · 4b423ac · 4b423ac
2 parents bf83a7f + 1da3a89
commit 4b423ac
Show file tree

Hide file tree

Showing 5 changed files with 106 additions and 28 deletions.
diff --git a/goodman_pipeline/core/__init__.py b/goodman_pipeline/core/__init__.py
@@ -34,6 +34,7 @@
                    normalize_master_flat,
                    ra_dec_to_deg,
                    read_fits,
+                   record_trace_information,
                    save_extracted,
                    search_comp_group,
                    setup_logging,

diff --git a/goodman_pipeline/core/core.py b/goodman_pipeline/core/core.py
@@ -1846,12 +1846,32 @@ def read_fits(full_path, technique='Unknown'):
         ccd.header.add_blank('-- GSP END --', after='GSP_WREJ')
 
     ccd.header.set('BUNIT', after='CCDSUM')
-    # ccd.header.set('', value='', comment='')
-    # ccd.header.set('', value='', comment='')
-    # ccd.header.set('', value='', comment='')
-    # ccd.header.set('', value='', comment='')
-    # ccd.header.set('', value='', comment='')
-    # ccd.header.set('', value='', comment='')
+
+    return ccd
+
+
+def record_trace_information(ccd, trace_info):
+    last_keyword = None
+    for info_key in trace_info:
+        info_value, info_comment = trace_info[info_key]
+        log.debug(
+            "Adding trace information: "
+            "{:s} = {:s} / {:s}".format(info_key,
+                                        str(info_value),
+                                        info_comment))
+
+        if last_keyword is None:
+            ccd.header.set(info_key,
+                           value=info_value,
+                           comment=info_comment)
+            last_keyword = info_key
+        else:
+            ccd.header.set(info_key,
+                           value=info_value,
+                           comment=info_comment,
+                           after=last_keyword)
+            last_keyword = info_key
+
     return ccd
 
 
@@ -2089,6 +2109,37 @@ def trace(ccd, model, trace_model, model_fitter, sampling_step, nsigmas=2, plots
 
     rms_error = np.sqrt(np.sum(np.array(sampling_differences))/len(sampling_differences))
 
+    log.debug("RMS Error of unclipped trace differences {:.3f}".format(rms_error))
+
+    clipped_values = sigma_clip(sampling_differences,
+                                sigma=2,
+                                iters=3,
+                                cenfunc=np.ma.median)
+    if np.ma.is_masked(clipped_values):
+        _sampling_axis = list(sampling_axis)
+        _sample_values = list(sample_values)
+
+        sampling_axis = []
+        sample_values = []
+        for i in range(len(clipped_values)):
+            if clipped_values[i] is not np.ma.masked:
+                sampling_axis.append(_sampling_axis[i])
+                sample_values.append(_sample_values[i])
+
+        log.debug("Re-fitting the trace for a better trace.")
+
+        fitted_trace = model_fitter(trace_model, sampling_axis, sample_values)
+
+        sampling_differences = [
+            (fitted_trace(sampling_axis[i]) - sample_values[i]) ** 2 for i in
+            range(len(sampling_axis))]
+
+        rms_error = np.sqrt(
+            np.sum(np.array(sampling_differences)) / len(sampling_differences))
+
+        log.debug(
+            "RMS Error after sigma-clipping trace differences {:.3f}".format(rms_error))
+
     trace_info = collections.OrderedDict()
 
     trace_info['GSP_TMOD'] = [fitted_trace.__class__.__name__,

diff --git a/goodman_pipeline/core/tests/test_core.py b/goodman_pipeline/core/tests/test_core.py
@@ -7,6 +7,7 @@
 from astropy.modeling import (models,
                               fitting)
 import astropy.units as u
+import collections
 import numpy as np
 import os
 import pandas
@@ -47,6 +48,7 @@
                     normalize_master_flat,
                     ra_dec_to_deg,
                     read_fits,
+                    record_trace_information,
                     save_extracted,
                     search_comp_group,
                     setup_logging,
@@ -679,6 +681,42 @@ def test_check_comp_group__lamp_does_not_exist(self):
         self.assertTrue(comp_group.equals(new_group))
 
 
+class RecordTraceInformationTest(TestCase):
+
+    def setUp(self):
+        self.ccd = CCDData(data=np.ones((800, 2000)),
+                           meta=fits.Header(),
+                           unit='adu')
+
+        self.all_keywords = ['GSP_TMOD',
+                             'GSP_TORD',
+                             'GSP_TC00',
+                             'GSP_TC01',
+                             'GSP_TC02',
+                             'GSP_TERR']
+
+        self.trace_info = collections.OrderedDict()
+
+        self.trace_info['GSP_TMOD'] = ['Polinomial1D',
+                                       'Model name used to fit trace']
+
+        self.trace_info['GSP_TORD'] = [2, 'Degree of the model used to fit '
+                                          'target trace']
+
+        self.trace_info['GSP_TC00'] = [500, 'Parameter c0']
+        self.trace_info['GSP_TC01'] = [1, 'Parameter c1']
+        self.trace_info['GSP_TC02'] = [2, 'Parameter c2']
+        self.trace_info['GSP_TERR'] = [0.5, 'RMS error of target trace']
+
+    def test_record_trace_information(self):
+        ccd = record_trace_information(ccd=self.ccd, trace_info=self.trace_info)
+        new_keys = [key for key in ccd.header.keys()]
+
+        self.assertTrue(all([key in new_keys for key in self.all_keywords]))
+        self.assertEqual(ccd.header['GSP_TMOD'], 'Polinomial1D')
+        self.assertEqual(ccd.header['GSP_TORD'], 2)
+
+
 class SpectroscopicModeTest(TestCase):
 
     def setUp(self):

diff --git a/goodman_pipeline/spectroscopy/redspec.py b/goodman_pipeline/spectroscopy/redspec.py
@@ -21,6 +21,7 @@
                     identify_targets,
                     trace_targets,
                     extraction,
+                    record_trace_information,
                     save_extracted)
 
 from ..core import (NoMatchFound,
@@ -389,27 +390,8 @@ def _run(self, data_container, extraction_type):
                         else:
                             target_number = 0
                         try:
-                            last_keyword = None
-                            for info_key in trace_info:
-                                info_value, info_comment = trace_info[info_key]
-                                self.log.debug(
-                                    "Adding trace information: "
-                                    "{:s} = {:s} / {:s}".format(info_key,
-                                                                str(info_value),
-                                                                info_comment))
-
-                                if last_keyword is None:
-                                    ccd.header.set(info_key,
-                                                  value=info_value,
-                                                  comment=info_comment)
-                                    last_keyword = info_key
-                                else:
-                                    ccd.header.set(info_key,
-                                                   value=info_value,
-                                                   comment=info_comment,
-                                                   after=last_keyword)
-                                    last_keyword = info_key
-
+                            ccd = record_trace_information(ccd=ccd,
+                                                           trace_info=trace_info)
                             # target extraction
                             extracted = extraction(
                                 ccd=ccd,
@@ -429,6 +411,8 @@ def _run(self, data_container, extraction_type):
                                                          value=saved_ccd.header['GSP_FNAM'],
                                                          comment='Science target file the lamp was extracted for.')
 
+                                    comp_lamp = record_trace_information(ccd=comp_lamp,
+                                                                         trace_info=trace_info)
                                     extracted_lamp = extraction(
                                         ccd=comp_lamp,
                                         target_trace=single_trace,
@@ -438,7 +422,9 @@ def _run(self, data_container, extraction_type):
                                         ccd=extracted_lamp,
                                         destination=self.args.destination,
                                         target_number=target_number)
+
                                     all_lamps.append(extracted_lamp)
+
                             extracted_target_and_lamps.append([extracted,
                                                                all_lamps])
 

diff --git a/goodman_pipeline/spectroscopy/wavelength.py b/goodman_pipeline/spectroscopy/wavelength.py
@@ -755,7 +755,9 @@ def _evaluate_solution(self, clipped_differences):
         self.rms_error = np.sqrt(
             np.sum(square_differences) / len(square_differences))
 
-        self.log.info('RMS Error : {:.3f}'.format(self.rms_error))
+        self.log.info('Wavelength solution RMS Error : {:.3f}'.format(
+            self.rms_error))
+
         return self.rms_error, self.n_points, self.n_rejections
 
     def _get_lines_in_lamp(self, ccddata_lamp=None):