Extract calibration modules from the reader

.pre-commit-config.yaml

@@ -2,12 +2,12 @@ exclude: '^$'
 fail_fast: false
 
 repos:
-  - repo: https://github.com/astral-sh/ruff-pre-commit
+  #- repo: https://github.com/astral-sh/ruff-pre-commit
     # Ruff version.
-    rev: 'v0.3.2'
-    hooks:
-      - id: ruff
-        args: [--fix]
+    #rev: 'v0.3.2'
+    #hooks:
+      #- id: ruff
+        #args: [--fix]
   - repo: https://github.com/pre-commit/pre-commit-hooks
     rev: v4.5.0
     hooks:

continuous_integration/environment.yaml

@@ -11,6 +11,7 @@ dependencies:
   - packaging
   - pytest
   - pytest-cov
+  - xarray
   - pip
   - pip:
     - docutils

pygac/calibration.py → pygac/calibration/noaa.py

@@ -40,14 +40,80 @@
 LOG = logging.getLogger(__name__)
 
 
+def calibrate(ds, custom_coeffs=None, coeffs_file=None):
+    """Calibrate the dataset `ds`, possibly using custom_coeffs or a coeff file.
+
+    Args:
+        ds: xarray Dataset object containing the data to calibrate.
+        custom_calibration: dictionary with a subset of user defined satellite specific
+                            calibration coefficients
+        calibration_file: path to json file containing default calibrations
+    """
+    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'
     PROVISIONAL = 'provisional'
     EXPERIMENTAL = 'experimental'
 
 
-class Calibrator(object):
+class Calibrator:
     """Factory class to create namedtuples holding the calibration coefficients.
 
     Attributes:
@@ -183,7 +249,7 @@
     @classmethod
     def _get_coeffs_version(cls, coeff_file_content):
         """Determine coefficient version."""
         md5_hash = hashlib.md5(coeff_file_content)
         digest = md5_hash.hexdigest()
         version_dict = cls.version_hashs.get(
             digest,
@@ -450,16 +516,18 @@
     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:
@@ -491,7 +559,7 @@
     # TsBB = A + B*TBB    (7.1.2.4-3)
     # NBB = c1*nu_e^3/(exp(c2*nu_e/TsBB) - 1)    (7.1.2.4-3)
     # where the constants of the Planck function are defined as c1 = 2*h*c^2, c2 = h*c/k_B.
-    # constatns
+    # constants
     c1 = 1.1910427e-5  # mW/m^2/sr/cm^{-4}
     c2 = 1.4387752  # cm K
     # coefficients

pygac/gac_klm.py

@@ -127,7 +127,8 @@
                      ("spacecraft_altitude_above_reference_ellipsoid", ">u2"),
                      ("angular_relationships", ">i2", (153, )),
                      ("zero_fill3", ">i2", (3, )),
-                     ("earth_location", ">i4", (102, )),
+                     ("earth_location", [("lats", ">i2"),
+                                         ("lons", ">i2")], (51,)),
                      ("zero_fill4", ">i4", (2, )),
                      # HRPT MINOR FRAME TELEMETRY
                      ("frame_sync", ">u2", (6, )),

pygac/gac_pod.py

@@ -42,7 +42,8 @@
                      ("number_of_meaningful_zenith_angles_and_earth_location_appended",
                       ">u1"),
                      ("solar_zenith_angles", "i1", (51, )),
-                     ("earth_location", ">i2", (102, )),
+                     ("earth_location", [("lats", ">i2"),
+                                         ("lons", ">i2")], (51,)),
                      ("telemetry", ">u4", (35, )),
                      ("sensor_data", ">u4", (682, )),
                      ("add_on_zenith", ">u2", (10, )),

pygac/gac_reader.py

@@ -28,6 +28,8 @@
 
 import logging
 
+import pygac.pygac_geotiepoints as gtp
+from pygac.pygac_geotiepoints import GAC_LONLAT_SAMPLE_POINTS
 from pygac.reader import Reader, ReaderError
 import pygac.pygac_geotiepoints as gtp
 
@@ -42,10 +44,11 @@ class GACReader(Reader):
     scan_freq = 2.0 / 1000.0
     # Max scanlines
     max_scanlines = 15000
+    lonlat_sample_points = GAC_LONLAT_SAMPLE_POINTS
 
     def __init__(self, *args, **kwargs):
         """Init the GAC reader."""
-        super(GACReader, self).__init__(*args, **kwargs)
+        super().__init__(*args, **kwargs)
         self.scan_width = 409
         self.lonlat_interpolator = gtp.gac_lat_lon_interpolator
 

pygac/klm_reader.py

@@ -761,8 +761,8 @@ def get_telemetry(self):
 
     def _get_lonlat(self):
         """Get the longitudes and latitudes."""
-        lats = self.scans["earth_location"][:, 0::2] / 1e4
-        lons = self.scans["earth_location"][:, 1::2] / 1e4
+        lats = self.scans["earth_location"]["lats"] / 1e4
+        lons = self.scans["earth_location"]["lons"] / 1e4
         return lons, lats
 
     def get_header_timestamp(self):
@@ -806,16 +806,15 @@ def _get_ir_channels_to_calibrate(self):
             ir_channels_to_calibrate = [4, 5]
         return ir_channels_to_calibrate
 
-    def postproc(self, channels):
+    def postproc(self, ds):
         """Apply KLM specific postprocessing.
 
         Masking out 3a/3b/transition.
         """
         switch = self.get_ch3_switch()
-        channels[:, :, 2][switch == 0] = np.nan
-        channels[:, :, 3][switch == 1] = np.nan
-        channels[:, :, 2][switch == 2] = np.nan
-        channels[:, :, 3][switch == 2] = np.nan
+        ds["channels"].loc[dict(channel_name="3a", scan_line_index=((switch==0) | (switch==2)))] = np.nan
+        ds["channels"].loc[dict(channel_name="3b", scan_line_index=((switch==1) | (switch==2)))] = np.nan
+
 
     def _adjust_clock_drift(self):
         """Adjust the geolocation to compensate for the clock error.

pygac/lac_klm.py

@@ -126,7 +126,8 @@
                      ("spacecraft_altitude_above_reference_ellipsoid", ">u2"),
                      ("angular_relationships", ">i2", (153, )),
                      ("zero_fill2", ">i2", (3, )),
-                     ("earth_location", ">i4", (102, )),
+                     ("earth_location", [("lats", ">i2"),
+                                         ("lons", ">i2")], (51,)),
                      ("zero_fill3", ">i4", (2, )),
                      # HRPT MINOR FRAME TELEMETRY
                      ("frame_sync", ">u2", (6, )),

pygac/lac_pod.py

@@ -37,10 +37,10 @@
                      ("time_code", ">u2", (3, )),
                      ("quality_indicators", ">u4"),
                      ("calibration_coefficients", ">i4", (10, )),
-                     ("number_of_meaningful_zenith_angles_and_earth_location_appended",
-                      ">u1"),
+                     ("number_of_meaningful_zenith_angles_and_earth_location_appended", ">u1"),
                      ("solar_zenith_angles", "i1", (51, )),
-                     ("earth_location", ">i2", (102, )),
+                     ("earth_location", [("lats", ">i2"),
+                                         ("lons", ">i2")], (51,)),
                      ("telemetry", ">u4", (35, )),
                      ("sensor_data", ">u4", (3414, )),
                      ("add_on_zenith", ">u2", (10, )),

pygac/lac_reader.py

@@ -26,6 +26,8 @@
 
 import logging
 
+import pygac.pygac_geotiepoints as gtp
+from pygac.pygac_geotiepoints import LAC_LONLAT_SAMPLE_POINTS
 from pygac.reader import Reader, ReaderError
 import pygac.pygac_geotiepoints as gtp
 
@@ -40,6 +42,7 @@ class LACReader(Reader):
     scan_freq = 6.0 / 1000.0
     # Max scanlines
     max_scanlines = 65535
+    lonlat_sample_points = LAC_LONLAT_SAMPLE_POINTS
 
     def __init__(self, *args, **kwargs):
         """Init the LAC reader."""

pygac/pod_reader.py

@@ -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(
@@ -341,8 +341,8 @@ def read_header(cls, filename, fileobj=None, header_date="auto"):
 
     @classmethod
     def _validate_tbm_header(cls, potential_tbm_header):
-        data_set_name = potential_tbm_header['data_set_name']
-        allowed_empty = (42*b'\x00' + b'  ')
+        data_set_name = potential_tbm_header["data_set_name"]
+        allowed_empty = (42*b"\x00" + b"  ")
         if data_set_name == allowed_empty:
             return potential_tbm_header.copy()
 
@@ -445,10 +445,9 @@ def _get_times(self):
     def _compute_missing_lonlat(self, missed_utcs):
         """Compute lon lat values using pyorbital."""
         tic = datetime.datetime.now()
-
         scan_rate = datetime.timedelta(milliseconds=1/self.scan_freq).total_seconds()
         sgeom = avhrr_gac(missed_utcs.astype(datetime.datetime),
-                          self.scan_points, frequency=scan_rate)
+                          self.lonlat_sample_points, frequency=scan_rate)
         t0 = missed_utcs[0].astype(datetime.datetime)
         s_times = sgeom.times(t0)
         tle1, tle2 = self.get_tle_lines()
@@ -546,8 +545,8 @@ def _adjust_clock_drift(self):
         LOG.debug("clock drift adjustment took %s", str(toc - tic))
 
     def _get_lonlat(self):
-        lats = self.scans["earth_location"][:, 0::2] / 128.0
-        lons = self.scans["earth_location"][:, 1::2] / 128.0
+        lats = self.scans["earth_location"]["lats"] / 128.0
+        lons = self.scans["earth_location"]["lons"] / 128.0
         return lons, lats
 
     def get_telemetry(self):

pygac/pygac_geotiepoints.py

@@ -28,6 +28,9 @@
 
 import numpy as np
 
+GAC_LONLAT_SAMPLE_POINTS = np.arange(4, 405, 8)
+LAC_LONLAT_SAMPLE_POINTS = np.arange(24, 2048, 40)
+
 
 def gac_lat_lon_interpolator(lons_subset, lats_subset):
     """Interpolate lat-lon values in the AVHRR GAC data.
@@ -37,16 +40,15 @@ def gac_lat_lon_interpolator(lons_subset, lats_subset):
     ranging from 5 to 405. Interpolate to full resolution.
 
     """
-    cols_subset = np.arange(4, 405, 8)
     cols_full = np.arange(409)
-    return lat_lon_interpolator(lons_subset, lats_subset, cols_subset, cols_full)
+    return lat_lon_interpolator(lons_subset, lats_subset, GAC_LONLAT_SAMPLE_POINTS, cols_full)
+
 
 
 def lac_lat_lon_interpolator(lons_subset, lats_subset):
     """Interpolate lat-lon values in the AVHRR LAC data."""
-    cols_subset = np.arange(24, 2048, 40)
     cols_full = np.arange(2048)
-    return lat_lon_interpolator(lons_subset, lats_subset, cols_subset, cols_full)
+    return lat_lon_interpolator(lons_subset, lats_subset, LAC_LONLAT_SAMPLE_POINTS, cols_full)
 
 
 def lat_lon_interpolator(lons_subset, lats_subset, cols_subset, cols_full):