Refactor all calibration outside of the reader

pygac/calibration.py → pygac/noaa_calibration.py

@@ -40,6 +40,64 @@
 LOG = logging.getLogger(__name__)
 
 
+def calibrate(ds, custom_coeffs=None, coeffs_file=None):
+    channels = ds["channels"].data
+    times = ds.coords["times"].data
+    scan_line_numbers = ds["scan_line_index"].data
+
+    calibration_coeffs = Calibrator(
+            ds.attrs["spacecraft_name"],
+            custom_coeffs=custom_coeffs,
+            coeffs_file=coeffs_file
+        )
+
+    # `times` is in nanosecond resolution. However, convertion from datetime64[ns] to datetime objects
+    # does not work, and anyway only the date is needed.
+    start_time = times[0]
+    start_date = start_time.astype("datetime64[D]").item()
+    year = start_date.year
+
+    delta = (start_time.astype("datetime64[D]") - start_time.astype("datetime64[Y]")).astype(int)
+    jday = delta.astype(int) + 1
+
+    corr = ds.attrs['sun_earth_distance_correction_factor']
+
+    # how many reflective channels are there ?
+    tot_ref = channels.shape[2] - 3
+
+    channels[:, :, 0:tot_ref] = calibrate_solar(
+        channels[:, :, 0:tot_ref],
+        np.arange(tot_ref),
+        year, jday,
+        calibration_coeffs,
+        corr
+    )
+
+    prt = ds["prt_counts"].data
+    ict = ds["ict_counts"].data
+    space = ds["space_counts"].data
+
+    ir_channels_to_calibrate = [3, 4, 5]
+
+    for chan in ir_channels_to_calibrate:
+        channels[:, :, chan - 6] = calibrate_thermal(
+            channels[:, :, chan - 6],
+            prt,
+            ict[:, chan - 3],
+            space[:, chan - 3],
+            scan_line_numbers,
+            chan,
+            calibration_coeffs
+        )
+
+    new_ds = ds.copy()
+    new_ds["channels"].data = channels
+
+    new_ds.attrs["calib_coeffs_version"] = calibration_coeffs.version
+
+    return new_ds
+
+
 class CoeffStatus(Enum):
     """Indicates the status of calibration coefficients."""
     NOMINAL = 'nominal'
@@ -450,16 +508,18 @@ def calibrate_thermal(counts, prt, ict, space, line_numbers, channel, cal):
     if channel == 3:
         zeros = ict < ict_threshold
         nonzeros = np.logical_not(zeros)
-
-        ict[zeros] = np.interp((zeros).nonzero()[0],
-                               (nonzeros).nonzero()[0],
-                               ict[nonzeros])
+        try:
+            ict[zeros] = np.interp((zeros).nonzero()[0],
+                                (nonzeros).nonzero()[0],
+                                ict[nonzeros])
+        except ValueError: # 3b has no valid data
+            return counts
         zeros = space < space_threshold
         nonzeros = np.logical_not(zeros)
 
         space[zeros] = np.interp((zeros).nonzero()[0],
-                                 (nonzeros).nonzero()[0],
-                                 space[nonzeros])
+                                (nonzeros).nonzero()[0],
+                                space[nonzeros])
 
     # convolving and smoothing PRT, ICT and SPACE values
     if lines > 51:

pygac/pod_reader.py

@@ -38,20 +38,20 @@
 
 import datetime
 import logging
+
 try:
     from enum import IntFlag
 except ImportError:
     # python version < 3.6, use a simple object without nice representation
     IntFlag = object
 
 import numpy as np
-
-from pyorbital.geoloc_instrument_definitions import avhrr_gac
 from pyorbital.geoloc import compute_pixels, get_lonlatalt
+from pyorbital.geoloc_instrument_definitions import avhrr_gac
 
 from pygac.clock_offsets_converter import get_offsets
 from pygac.correct_tsm_issue import TSM_AFFECTED_INTERVALS_POD, get_tsm_idx
-from pygac.reader import Reader, ReaderError, NoTLEData
+from pygac.reader import NoTLEData, Reader, ReaderError
 from pygac.slerp import slerp
 from pygac.utils import file_opener
 
@@ -246,7 +246,7 @@ class PODReader(Reader):
     def correct_scan_line_numbers(self):
         """Correct the scan line numbers."""
         # Perform common corrections first.
-        super(PODReader, self).correct_scan_line_numbers()
+        super().correct_scan_line_numbers()
 
         # cleaning up the data
         min_scanline_number = np.amin(

pygac/reader.py

@@ -24,26 +24,24 @@
 
 Can't be used as is, has to be subclassed to add specific read functions.
 """
-from abc import ABCMeta, abstractmethod
 import datetime
 import logging
-
-import numpy as np
 import os
 import re
-import six
 import types
 import warnings
-import pyorbital
+from abc import ABCMeta, abstractmethod
 
-from pygac.utils import (centered_modulus,
-                         calculate_sun_earth_distance_correction,
-                         get_absolute_azimuth_angle_diff)
-from pyorbital.orbital import Orbital
+import numpy as np
+import pyorbital
+import six
+from packaging.version import Version
 from pyorbital import astronomy
-from pygac.calibration import Calibrator, calibrate_solar, calibrate_thermal
+from pyorbital.orbital import Orbital
+
 from pygac import gac_io
-from packaging.version import Version
+from pygac.noaa_calibration import calibrate
+from pygac.utils import calculate_sun_earth_distance_correction, centered_modulus, get_absolute_azimuth_angle_diff
 
 LOG = logging.getLogger(__name__)
 
@@ -163,16 +161,6 @@ def filename(self, filepath):
             else:
                 self._filename = os.path.basename(filepath)
 
-    @property
-    def calibration(self):
-        """Get the property 'calibration'."""
-        calibration = Calibrator(
-            self.spacecraft_name,
-            custom_coeffs=self.custom_calibration,
-            coeffs_file=self.calibration_file
-        )
-        return calibration
-
     @abstractmethod
     def read(self, filename, fileobj=None):  # pragma: no cover
         """Read the GAC/LAC data.
@@ -509,55 +497,86 @@ def get_sun_earth_distance_correction(self):
         return calculate_sun_earth_distance_correction(jday)
 
     def update_meta_data(self):
-        """Add some metd data to the meta_data dicitonary."""
-        if 'sun_earth_distance_correction_factor' not in self.meta_data:
-            self.meta_data['sun_earth_distance_correction_factor'] = (
+        """Add some meta data to the meta_data dicitonary."""
+        meta_data = self.meta_data
+        self._update_meta_data_object(meta_data)
+        if 'gac_header' not in meta_data:
+            meta_data['gac_header'] = self.head.copy()
+
+    def _update_meta_data_object(self, meta_data):
+        if 'sun_earth_distance_correction_factor' not in meta_data:
+            meta_data['sun_earth_distance_correction_factor'] = (
                 self.get_sun_earth_distance_correction())
-        if 'midnight_scanline' not in self.meta_data:
-            self.meta_data['midnight_scanline'] = self.get_midnight_scanline()
-        if 'missing_scanlines' not in self.meta_data:
-            self.meta_data['missing_scanlines'] = self.get_miss_lines()
-        if 'gac_header' not in self.meta_data:
-            self.meta_data['gac_header'] = self.head.copy()
-        self.meta_data['calib_coeffs_version'] = self.calibration.version
+        if 'midnight_scanline' not in meta_data:
+            meta_data['midnight_scanline'] = self.get_midnight_scanline()
+        if 'missing_scanlines' not in meta_data:
+            meta_data['missing_scanlines'] = self.get_miss_lines()
+
+    def read_as_dataset(self, file_to_read):
+        self.read(file_to_read)
+        return self.create_dataset()
+
+    def create_dataset(self):
+        import xarray as xr
+
+        head = dict(zip(self.head.dtype.names, self.head.item()))
+        scans = self.scans
+
+        times = self.get_times()
+        line_numbers = scans["scan_line_number"]
+        counts = self.get_counts()
+
+        scan_size = counts.shape[1]
+
+        if scan_size == 409:
+            sub_columns = np.arange(4, 405, 8)
+        else:
+            sub_columns = np.arange(24, 2048, 40)
+
+        if counts.shape[-1] == 5:
+            channel_names = ["1", "2", "3", "4", "5"]
+            ir_channel_names = ["3", "4", "5"]
+        else:
+            channel_names = ["1", "2", "3a", "3b", "4", "5"]
+            ir_channel_names = ["3b", "4", "5"]
+        channels = xr.DataArray(counts, dims=["scan_line_index", "columns", "channel_name"],
+                                coords=dict(scan_line_index=line_numbers, channel_name=channel_names,
+                                            columns=np.arange(scan_size))
+                                )
+        channels = channels.assign_coords(times=("scan_line_index", times))
+        lons, lats = self._get_lonlat()
+        longitudes = xr.DataArray(lons, dims=["scan_line_index", "columns"],
+                                  coords=dict(scan_line_index=line_numbers,
+                                              columns=sub_columns))
+        latitudes = xr.DataArray(lats, dims=["scan_line_index", "columns"],
+                                 coords=dict(scan_line_index=line_numbers,
+                                             columns=sub_columns))
+        channels = channels.assign_coords(longitudes=(("scan_line_index", "columns"),
+                                                      longitudes.reindex_like(channels).data),
+                                          latitudes=(("scan_line_index", "columns"),
+                                                     latitudes.reindex_like(channels).data))
+
+        prt, ict, space = self.get_telemetry()
+        prt = xr.DataArray(prt, dims=["scan_line_index"])
+        ict = xr.DataArray(ict, dims=["scan_line_index", "ir_channel_name"],
+                           coords=dict(ir_channel_name=ir_channel_names))
+        space = xr.DataArray(space, dims=["scan_line_index", "ir_channel_name"])
+        ds = xr.Dataset(dict(channels=channels, prt_counts=prt, ict_counts=ict, space_counts=space), attrs=head)
+
+
+        ds.attrs["spacecraft_name"] = self.spacecraft_name
+        self._update_meta_data_object(ds.attrs)
+        return ds
 
     def get_calibrated_channels(self):
         """Calibrate and return the channels."""
-        channels = self.get_counts()
-        self.get_times()
-        times = self.times
-        self.update_meta_data()
-        year = times[0].year
-        delta = times[0].date() - datetime.date(year, 1, 1)
-        jday = delta.days + 1
-
-        corr = self.meta_data['sun_earth_distance_correction_factor']
-        calibration_coeffs = self.calibration
-
-        # how many reflective channels are there ?
-        tot_ref = channels.shape[2] - 3
-
-        channels[:, :, 0:tot_ref] = calibrate_solar(
-            channels[:, :, 0:tot_ref],
-            np.arange(tot_ref),
-            year, jday,
-            calibration_coeffs,
-            corr
-        )
-        prt, ict, space = self.get_telemetry()
+        ds = self.create_dataset()
 
-        ir_channels_to_calibrate = self._get_ir_channels_to_calibrate()
-
-        for chan in ir_channels_to_calibrate:
-            channels[:, :, chan - 6] = calibrate_thermal(
-                channels[:, :, chan - 6],
-                prt,
-                ict[:, chan - 3],
-                space[:, chan - 3],
-                self.scans["scan_line_number"],
-                chan,
-                calibration_coeffs
-            )
+        noaa_calibration_kwargs = dict(custom_coeffs=self.custom_calibration,
+                                       coeffs_file=self.calibration_file)
+        calibrated_ds = calibrate(ds, **noaa_calibration_kwargs)
+        channels = calibrated_ds["channels"].data
+        self.meta_data["calib_coeffs_version"] = calibrated_ds.attrs["calib_coeffs_version"]
 
         # Mask out corrupt values
         channels[self.mask] = np.nan

pygac/tests/test_calibrate_klm.py

@@ -25,13 +25,10 @@
 
 
 import unittest
-try:
-    import mock
-except ImportError:
-    from unittest import mock
+
 import numpy as np
 
-from pygac.calibration import Calibrator, calibrate_solar, calibrate_thermal
+from pygac.noaa_calibration import Calibrator, calibrate_solar, calibrate_thermal
 
 
 class TestGenericCalibration(unittest.TestCase):

pygac/tests/test_calibrate_pod.py

@@ -25,13 +25,10 @@
 
 
 import unittest
-try:
-    import mock
-except ImportError:
-    from unittest import mock
+
 import numpy as np
 
-from pygac.calibration import Calibrator, calibrate_solar, calibrate_thermal
+from pygac.noaa_calibration import Calibrator, calibrate_solar, calibrate_thermal
 
 
 class TestGenericCalibration(unittest.TestCase):

pygac/tests/test_calibration_coefficients.py

@@ -24,13 +24,14 @@
 
 import sys
 import unittest
+
 try:
     import mock
 except ImportError:
     from unittest import mock
 import numpy as np
 
-from pygac.calibration import Calibrator, calibrate_solar, CoeffStatus
+from pygac.noaa_calibration import Calibrator, CoeffStatus, calibrate_solar
 
 # dummy user json file including only noaa19 data with changed channel 1 coefficients
 user_json_file = b"""{
@@ -118,7 +119,7 @@
 
 class TestCalibrationCoefficientsHandling(unittest.TestCase):
 
-    @mock.patch('pygac.calibration.open', mock.mock_open(read_data=user_json_file))
+    @mock.patch('pygac.noaa_calibration.open', mock.mock_open(read_data=user_json_file))
     def test_user_coefficients_file(self):
         if sys.version_info.major < 3:
             cal = Calibrator('noaa19', coeffs_file="/path/to/unknow/defaults.json")