.bumpversion.cfg

 [bumpversion]
-current_version = 1.2.1
+current_version = 1.3.0
 commit = True
 tag = True
 

.stickler.yml

 linters:
   flake8:
     fixer: true
+    max-line-length: 120
 fixers:
   enable: true

CHANGELOG.md

+## Version 1.3.0 (2019/12/05)
+
+
+### Pull Requests Merged
+
+#### Features added
+
+* [PR 45](https://github.com/pytroll/pygac/pull/45) - Add LAC support ([5](https://github.com/pytroll/pygac/issues/5))
+* [PR 42](https://github.com/pytroll/pygac/pull/42) - Update documentation
+* [PR 41](https://github.com/pytroll/pygac/pull/41) - Add meta_data dictionary to reader class.
+
+In this release 3 pull requests were closed.
+
+
 ## Version 1.2.1 (2019/11/15)
 
 ### Pull Requests Merged

README.md

 pygac
 =====
 
-[![Build Status](https://travis-ci.org/adybbroe/pygac.png?branch=feature-clock)](https://travis-ci.org/adybbroe/pygac)
-[![Coverage Status](https://coveralls.io/repos/adybbroe/pygac/badge.png?branch=feature-clock)](https://coveralls.io/r/adybbroe/pygac?branch=develop)
-[![Code Health](https://landscape.io/github/pytroll/pygac/develop/landscape.png)](https://landscape.io/github/pytroll/pygac/develop)
+[![Build Status](https://travis-ci.org/adybbroe/pygac.png?branch=feature-clock)](https://travis-ci.org/pytroll/pygac)
+[![Coverage Status](https://coveralls.io/repos/adybbroe/pygac/badge.png?branch=feature-clock)](https://coveralls.io/r/pytroll/pygac?branch=develop)
 
 
-A python package to read, calibrate and navigate NOAA AVHRR GAC data.
+A python package to read, calibrate and navigate NOAA and Metop AVHRR GAC and LAC data.
+
 
-Abhay Devatshale, Martin Raspaud and Adam Dybbroe 
-April 2014, Norrkoping, Sweden

debian/changelog

+pygac (1.3.0-1) unstable; urgency=medium
+
+  * New upstream release.
+  * Update copyright file.
+  * Refresh all patches.
+
+ -- Antonio Valentino <antonio.valentino@tiscali.it>  Fri, 06 Dec 2019 06:58:30 +0000
+
 pygac (1.2.1-1) unstable; urgency=medium
 
   * New upstream release.

debian/copyright

@@ -10,11 +10,13 @@ Copyright: 2010-2019 Martin Raspaud <martin.raspaud@smhi.se>
                      Sara Hornquist <sara.hornquist@smhi.se>
                      Stephan Finkensieper <stephan.finkensieper@dwd.de>
                      Cornelia Schlundt <cornelia.schlundt@dwd.de>
+                     Sajid Pareeth <sajid.pareeth@fmach.it>
+                     Nina Hakansson <nina.hakansson@smhi.se>
                      Pytroll Developers
 License: GPL-3+
 
 Files: debian/*
-Copyright: 2018 Antonio Valentino <antonio.valentino@tiscali.it>
+Copyright: 2018-2019 Antonio Valentino <antonio.valentino@tiscali.it>
 License: GPL-3+
 
 License: GPL-3+

debian/patches/0002-Fix-config.patch

@@ -8,10 +8,10 @@ Subject: Fix config
  2 files changed, 3 insertions(+), 3 deletions(-)
 
 diff --git a/pygac/__init__.py b/pygac/__init__.py
-index 9225981..a30e971 100644
+index 0940aa4..cf52fa8 100644
 --- a/pygac/__init__.py
 +++ b/pygac/__init__.py
-@@ -30,7 +30,7 @@ try:
+@@ -31,7 +31,7 @@ try:
      CONFIG_FILE = os.environ['PYGAC_CONFIG_FILE']
  except KeyError:
      LOG.error('Environment variable PYGAC_CONFIG_FILE not set!')
@@ -21,10 +21,10 @@ index 9225981..a30e971 100644
  if not os.path.exists(CONFIG_FILE) or not os.path.isfile(CONFIG_FILE):
      LOG.warning(
 diff --git a/setup.py b/setup.py
-index cefd25e..5ffbda7 100644
+index 22bd50c..52788d8 100644
 --- a/setup.py
 +++ b/setup.py
-@@ -104,8 +104,8 @@ if __name__ == '__main__':
+@@ -107,8 +107,8 @@ if __name__ == '__main__':
                              'scipy>=0.8.0',
                              'python-geotiepoints>=1.1.8'],
            scripts=[os.path.join('bin', item) for item in os.listdir('bin')],

doc/source/api.rst

@@ -9,10 +9,6 @@ This is done by reading the first three bytes of the data set. If they contain t
 
 GAC POD reader
 --------------
-
-.. automodule:: pygac.gac_pod
-   :members:
-   :undoc-members:
    
 As the format of GAC POD header is changed twice (once in 1992 and again in 1994), there are currently three readers integrated in the *pygac* to read POD header.
 
@@ -45,6 +41,9 @@ It is evident that year, month values are jumping and have non-sense values.
 
 Whenever possible, *pygac* uses RPY corrections along with other orbital parameters to compute accurate satellite location (e.g. instead of assuming constant altitude). However, RPY corrections are not available for all NOAA satellites. In case of the majority of the POD family satellites, these corrections are set to zero.
 
+.. automodule:: pygac.gac_pod
+   :members:
+   :undoc-members:
 
 GAC KLM reader
 --------------
@@ -54,34 +53,29 @@ GAC KLM reader
    :undoc-members:
 
 
-Computation of geolocation
---------------
+GAC Reader
+----------
 
-.. automodule:: pygac.gac_klm
+.. automodule:: pygac.gac_reader
    :members:
    :undoc-members:
 
 
+Computation of geolocation
+--------------------------
+
 Each GAC row has total 409 pixels. But lat-lon values are provided for every eigth pixel starting from pixel 5 and ending at pixel 405. Using Numpy, Scipy and Pyresample packages, inter- and extrapolation is carried out to obtain geolocation for each pixel along the scanline.
 
 If the GAC data belongs to POD family, then clock drift errors are used to adjust existing Lat-Lon information. Here, *pygac* makes use of PyOrbital package, which is a part of PyTroll suite of Python interface developed to process meteorological satellite data (further information here: http://www.pytroll.org/ and https://github.com/mraspaud/pyorbital). *pygac* interpolates the clock offset and adjusts the nominal scan times to the actual scan times. Since the geolocation was computed using the nominal scan times, *pygac* interpolates the latitudes and longitudes to the actual scan times using spherical linear interpolation, aka slerp. However, in the case of a clock drift error greater than the scan rate of the dataset, the latitude and longitude for each pixel of the scan lines that cannot have an interpolated geolocation (typically at the start or end of the dataset) are recomputed. This is done using pyorbital, which in turn uses TLEs to compute the position of the satellite at each scan time and the instrument geometry compute the longitude and latitude of each pixel of the dataset. Since this operation can be quite costly, the interpolation is prefered whenever possible.
 
 Computation of angles
---------------
-
-.. automodule:: pygac.gac_reader
-   :members:
-   :undoc-members:
+---------------------
 
 The azimuth angles are calculated using get_alt_az and get_observer_look from pyorbital. In pyorbital documentation there is a link describing how calculations are done http://celestrak.com/columns/v02n02/. The azimuth described in the link is measured as clockwise from North instead of counter-clockwise from South. Counter clockwise from south would be the standard for a right-handed orthogonal coordinate system. Pygac was updated to use the same definition for angles as pyorbital (2019, September, version > 1.1.0). Previous versions used azimuth +/-180 degrees, which correspond to degrees clockwise from south. All angles are converted to degrees. All azimuth angles are converted to range ]-180, 180] (2019 October version > 1.1.0). Note that ]-180, 180] is an open interval.
 
 
 GAC calibration/inter-calibration
---------------
-
-.. automodule:: pygac.gac_klm
-   :members:
-   :undoc-members:
+---------------------------------
 
 *pygac* currently supports calibration of all GAC data from AVHRRs onboard
  NOAA-7 and onwards, including Metop satellites.
@@ -103,10 +97,6 @@ In some cases it was found that, apart from the reset values, even the readings
 GAC I/O module
 --------------
 
-.. automodule:: pygac.gac_pod
-   :members:
-   :undoc-members:
-
 The I/O module generates three HDF5 files, one containing reflectances, brightness temperatures, and lat/lon information. The other output file contains solar and satellite zenith and azimuth angles. And the third file contains quality flags.
 
 The output file name format is:
@@ -147,13 +137,13 @@ In some orbits, the latitude and longitude information contains corrupt values f
 
 For extremely warm and cold temperatures, the channel 3b is saturating producing irrelevant brightness temperatures. Such saturation is often not flagged in quality information. *pygac* currently uses a simple if_else construct to constrain valid range of Bts (170.0K<BT<350.0K). Such condition is also applied to split-window channels.
 
-
-Supplement A: Structure of an output file containing reflectances and brightness temperatures
---------------
-
-.. automodule:: pygac.gac_klm
+.. automodule:: pygac.gac_io
    :members:
    :undoc-members:
+
+
+Supplement A: Structure of an output file containing reflectances and brightness temperatures
+---------------------------------------------------------------------------------------------
    
 Input: L1b file: ``NSS.GHRR.NN.D06279.S1800.E1955.B0711012.GC``
 
@@ -402,11 +392,7 @@ Output::
 
 
 Supplement B: Structure of an output file containing Sun and satellite positions
---------------
-
-.. automodule:: pygac.gac_klm
-   :members:
-   :undoc-members:
+--------------------------------------------------------------------------------
 
 Input: L1b file: ``NSS.GHRR.NN.D06279.S1800.E1955.B0711012.GC``
 
@@ -597,11 +583,8 @@ Output::
 
 
 Supplement C: Structure of an output file containing quality flags
---------------
+------------------------------------------------------------------
 
-.. automodule:: pygac.gac_klm
-   :members:
-   :undoc-members:
 
 The file that contains quality flags has following information.
 

doc/source/usage.rst

 Usage
 -----
 
+Standalone
+~~~~~~~~~~
+
 Copy the template file ``etc/pygac.cfg.template`` to ``pygac.cfg`` and place
 it in a directory as you please. Set the environment variable PYGAC_CONFIG_FILE
 pointing to the file. E.g.::
@@ -26,3 +29,7 @@ portion of the GAC orbit that user wants to process.  The first scanline number
 starts at 0. If zeroes are specified at both locations, then the entire orbit
 will be processed.
 
+In Satpy
+~~~~~~~~
+It is also possible to use pygac as a library in Satpy, see this `example notebook
+<https://github.com/pytroll/pytroll-examples/blob/master/satpy/avhrr_l1b_gaclac.ipynb>`_.

pygac/__init__.py

@@ -22,6 +22,7 @@
 
 import logging
 import os
+import numpy as np
 
 from pygac.version import __version__  # noqa
 
@@ -52,3 +53,18 @@ def centered_modulus(array, divisor):
     arr = array % divisor
     arr[arr > divisor / 2] -= divisor
     return arr
+
+
+def calculate_sun_earth_distance_correction(jday):
+    """Calculate the sun earth distance correction.
+
+    In 2008 3-4 different equations of ESD were considered.
+    This one was chosen as it at the time gave reflectances most closely
+    matching the PATMOS-x data provided then by Andy Heidinger.
+
+    Formula might need to be reconsidered if jday is updated to a float.
+
+    """
+    # Earth-Sun distance correction factor
+    corr = 1.0 - 0.0334 * np.cos(2.0 * np.pi * (jday - 2) / 365.25)
+    return corr

pygac/gac_io.py

@@ -83,7 +83,11 @@ def save_gac(satellite_name,
              bt3, bt4, bt5,
              sun_zen, sat_zen, sun_azi, sat_azi, rel_azi,
              qual_flags, start_line, end_line,
-             gac_file, midnight_scanline, miss_lines):
+             gac_file, meta_data):
+
+    midnight_scanline = meta_data['midnight_scanline']
+    miss_lines = meta_data['missing_scanlines']
+    corr = meta_data['sun_earth_distance_correction_factor']
 
     last_scan_line_number = qual_flags[-1, 0]
 
@@ -199,10 +203,6 @@ def save_gac(satellite_name,
     starttime = start.strftime("%H%M%S%f")[:-5]
     enddate = end.strftime("%Y%m%d")
     endtime = end.strftime("%H%M%S%f")[:-5]
-    jday = int(start.strftime("%j"))
-
-    # Earth-Sun distance correction factor
-    corr = 1.0 - 0.0334 * np.cos(2.0 * np.pi * (jday - 2) / 365.25)
 
     # Apply scaling & offset
     bt3 -= 273.15

pygac/gac_pod.py

-#!/usr/bin/env python
-
-# -*- coding: utf-8 -*-
-
-# Copyright (c) 2014, 2015 Abhay Devasthale
+#!/usr/bin/python
+# Copyright (c) 2014-2016.
+#
 
 # Author(s):
 
 #   Abhay Devasthale <abhay.devasthale@smhi.se>
 #   Adam Dybbroe <adam.dybbroe@smhi.se>
+#   Sajid Pareeth <sajid.pareeth@fmach.it>
 #   Martin Raspaud <martin.raspaud@smhi.se>
 
 # This work was done in the framework of ESA-CCI-Clouds phase I
@@ -25,127 +24,24 @@
 
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-
-"""Read a gac file.
+"""Read a GAC POD file.
 
 Format specification can be found here:
 http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c2/sec2-0.htm
 http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c3/sec3-1.htm
-
 """
 
-from __future__ import print_function
-
-import numpy as np
-from pygac.gac_reader import GACReader, inherit_doc
-from pygac.correct_tsm_issue import TSM_AFFECTED_INTERVALS_POD, get_tsm_idx
-import datetime
-from pygac import gac_io
-
 import logging
-LOG = logging.getLogger(__name__)
-
-# common header
-header0 = np.dtype([("noaa_spacecraft_identification_code", ">u1"),
-                    ("data_type_code", ">u1"),
-                    ("start_time", ">u2", (3, )),
-                    ("number_of_scans", ">u2"),
-                    ("end_time", ">u2", (3, ))])
-
-# until 8 september 1992
-header1 = np.dtype([("noaa_spacecraft_identification_code", ">u1"),
-                    ("data_type_code", ">u1"),
-                    ("start_time", ">u2", (3, )),
-                    ("number_of_scans", ">u2"),
-                    ("end_time", ">u2", (3, )),
-                    ("processing_block_id", "S7"),
-                    ("ramp_auto_calibration", ">u1"),
-                    ("number_of_data_gaps", ">u2"),
-                    ("dacs_quality", ">u1", (6, )),
-                    ("calibration_parameter_id", ">i2"),
-                    ("dacs_status", ">u1"),
-                    ("spare1", ">i1", (5, )),
-                    ("data_set_name", "S44"),
-                    ("spare2", ">u2", (1568, ))])
 
-header2 = np.dtype([("noaa_spacecraft_identification_code", ">u1"),
-                    ("data_type_code", ">u1"),
-                    ("start_time", ">u2", (3, )),
-                    ("number_of_scans", ">u2"),
-                    ("end_time", ">u2", (3, )),
-                    ("processing_block_id", "S7"),
-                    ("ramp_auto_calibration", ">u1"),
-                    ("number_of_data_gaps", ">u2"),
-                    ("dacs_quality", ">u1", (6, )),
-                    ("calibration_parameter_id", ">i2"),
-                    ("dacs_status", ">u1"),
-                    ("spare1", ">i1", (5, )),
-                    ("data_set_name", "S42"),
-                    ("blankfill", "S2"),
-                    ("year_of_epoch", ">u2"),
-                    ("julian_day_of_epoch", ">u2"),
-                    ("millisecond_utc_epoch_time_of_day", ">u4"),
-                    # Keplerian orbital elements
-                    ("semi_major_axis", ">f8"),
-                    ("eccentricity", ">f8"),
-                    ("inclination", ">f8"),
-                    ("argument_of_perigee", ">f8"),
-                    ("right_ascension", ">f8"),
-                    ("mean_anomaly", ">f8"),
-                    # cartesian inertial true date of elements
-                    ("x_component_of_position_vector", ">f8"),
-                    ("y_component_of_position_vector", ">f8"),
-                    ("z_component_of_position_vector", ">f8"),
-                    ("x_dot_component_of_position_vector", ">f8"),
-                    ("y_dot_component_of_position_vector", ">f8"),
-                    ("z_dot_component_of_position_vector", ">f8"),
-                    ("spare2", ">u2", (1516, ))])
+import numpy as np
 
-header3 = np.dtype([("noaa_spacecraft_identification_code", ">u1"),
-                    ("data_type_code", ">u1"),
-                    ("start_time", ">u2", (3, )),
-                    ("number_of_scans", ">u2"),
-                    ("end_time", ">u2", (3, )),
-                    ("processing_block_id", "S7"),
-                    ("ramp_auto_calibration", ">u1"),
-                    ("number_of_data_gaps", ">u2"),
-                    ("dacs_quality", ">u1", (6, )),
-                    ("calibration_parameter_id", ">i2"),
-                    ("dacs_status", ">u1"),
-                    ("reserved_for_mounting_and_fixed_attitude_correction_indicator",
-                     ">i1"),
-                    ("nadir_earth_location_tolerance", ">i1"),
-                    ("spare1", ">i1"),
-                    ("start_of_data_set_year", ">u2"),
-                    ("data_set_name", "S44"),
-                    ("year_of_epoch", ">u2"),
-                    ("julian_day_of_epoch", ">u2"),
-                    ("millisecond_utc_epoch_time_of_day", ">u4"),
-                    # Keplerian orbital elements
-                    ("semi_major_axis", ">i4"),
-                    ("eccentricity", ">i4"),
-                    ("inclination", ">i4"),
-                    ("argument_of_perigee", ">i4"),
-                    ("right_ascension", ">i4"),
-                    ("mean_anomaly", ">i4"),
-                    # cartesian inertial true date of elements
-                    ("x_component_of_position_vector", ">i4"),
-                    ("y_component_of_position_vector", ">i4"),
-                    ("z_component_of_position_vector", ">i4"),
-                    ("x_dot_component_of_position_vector", ">i4"),
-                    ("y_dot_component_of_position_vector", ">i4"),
-                    ("z_dot_component_of_position_vector", ">i4"),
-                    # future use
-                    ("yaw_fixed_error_correction", ">i2"),
-                    ("roll_fixed_error_correction", ">i2"),
-                    ("pitch_fixed_error_correction", ">i2"),
-                    ("spare2", ">u2", (1537, ))])
+from pygac.pod_reader import PODReader, main_pod
+from pygac.gac_reader import GACReader
 
+LOG = logging.getLogger(__name__)
 
 scanline = np.dtype([("scan_line_number", ">i2"),
                      ("time_code", ">u2", (3, )),
-                     # ("time_code", ">u1", (6, )),
                      ("quality_indicators", ">u4"),
                      ("calibration_coefficients", ">i4", (10, )),
                      ("number_of_meaningful_zenith_angles_and_earth_location_appended",
@@ -159,280 +55,26 @@ scanline = np.dtype([("scan_line_number", ">i2"),
                      ("spare3", "u2", (11, ))])
 
 
-@inherit_doc
-class GACPODReader(GACReader):
-
-    spacecrafts_orbital = {25: 'tiros n',
-                           2: 'noaa 6',
-                           4: 'noaa 7',
-                           6: 'noaa 8',
-                           7: 'noaa 9',
-                           8: 'noaa 10',
-                           1: 'noaa 11',
-                           5: 'noaa 12',
-                           3: 'noaa 14',
-                           }
-    spacecraft_names = {25: 'tirosn',
-                        2: 'noaa6',
-                        4: 'noaa7',
-                        6: 'noaa8',
-                        7: 'noaa9',
-                        8: 'noaa10',
-                        1: 'noaa11',
-                        5: 'noaa12',
-                        3: 'noaa14',
-                        }
-
-    tsm_affected_intervals = TSM_AFFECTED_INTERVALS_POD
-
-    def correct_scan_line_numbers(self):
-        # Perform common corrections first.
-        super(GACPODReader, self).correct_scan_line_numbers()
-
-        # cleaning up the data
-        min_scanline_number = np.amin(np.absolute(self.scans["scan_line_number"][:]))
-        if self.scans["scan_line_number"][0] == self.scans["scan_line_number"][-1] + 1:
-            while self.scans["scan_line_number"][0] != min_scanline_number:
-                self.scans = np.roll(self.scans, -1)
-        else:
-            while self.scans["scan_line_number"][0] != min_scanline_number:
-                self.scans = self.scans[1:]
-
-        self.scans = self.scans[self.scans["scan_line_number"] != 0]
-
-    def read(self, filename):
-        super(GACPODReader, self).read(filename=filename)
-
-        # choose the right header depending on the date
-        with open(filename) as fd_:
-            head = np.fromfile(fd_, dtype=header0, count=1)[0]
-            year, jday, _ = self.decode_timestamps(head["start_time"])
-
-            start_date = (datetime.date(year, 1, 1) + datetime.timedelta(days=int(jday) - 1))
-
-            if start_date < datetime.date(1992, 9, 8):
-                header = header1
-            elif start_date <= datetime.date(1994, 11, 15):
-                header = header2
-            else:
-                header = header3
-
-        with open(filename) as fd_:
-            self.head = np.fromfile(fd_, dtype=header, count=1)[0]
-            self.scans = np.fromfile(fd_,
-                                     dtype=scanline,
-                                     count=self.head["number_of_scans"])
-        self.correct_scan_line_numbers()
-        self.spacecraft_id = self.head["noaa_spacecraft_identification_code"]
-        if self.spacecraft_id == 1 and start_date < datetime.date(1982, 1, 1):
-            self.spacecraft_id = 25
-        self.spacecraft_name = self.spacecraft_names[self.spacecraft_id]
-        LOG.info(
-            "Reading %s data", self.spacecrafts_orbital[self.spacecraft_id])
-
-        return self.head, self.scans
-
-    def get_header_timestamp(self):
-        year, jday, msec = self.decode_timestamps(self.head["start_time"])
-        try:
-            return self.to_datetime(self.to_datetime64(year=year, jday=jday,
-                                                       msec=msec))
-        except ValueError as err:
-            raise ValueError('Corrupt header timestamp: {0}'.format(err))
-
-    @staticmethod
-    def decode_timestamps(encoded):
-        """
-        Decode timestamps.
-        @return: year, day of year, milliseconds since 00:00
-        """
-        ndims = len(encoded.shape)
-        if ndims == 1:
-            # Single header timestamp
-            enc0 = encoded[0]
-            enc1 = encoded[1]
-            enc2 = encoded[2]
-        elif ndims == 2:
-            # Scanline timestamps
-            enc0 = encoded[:, 0]
-            enc1 = encoded[:, 1]
-            enc2 = encoded[:, 2]
-        else:
-            raise ValueError('Invalid timestamp dimension')
-
-        year = enc0 >> 9
-        year = np.where(year > 75, year + 1900, year + 2000)
-        jday = (enc0 & 0x1FF)
-        msec = ((np.uint32(enc1 & 2047) << 16) | (np.uint32(enc2)))
-
-        return year, jday, msec
+class GACPODReader(GACReader, PODReader):
+    """The GAC POD reader class.
 
-    def _get_times(self):
-        return self.decode_timestamps(self.scans["time_code"])
+    The `scan_points` attributes provides the position of the longitude and latitude points to
+    compute relative to the full swath width.
 
-    def _adjust_clock_drift(self):
-        """Adjust the geolocation to compensate for the clock error.
+    The offset attribute tells where in the file the scanline data starts.
+    """
 
-        TODO: bad things might happen when scanlines are skipped.
-        """
-        tic = datetime.datetime.now()
-        self.get_times()
-        from pygac.clock_offsets_converter import get_offsets
-        try:
-            offset_times, clock_error = get_offsets(self.spacecraft_name)
-        except KeyError:
-            LOG.info("No clock drift info available for %s",
-                     self.spacecraft_name)
-        else:
-            offset_times = np.array(offset_times, dtype='datetime64[ms]')
-            offsets = np.interp(self.utcs.astype(np.uint64),
-                                offset_times.astype(np.uint64),
-                                clock_error)
-            LOG.info("Adjusting for clock drift of %s to %s",
-                     str(min(offsets)),
-                     str(max(offsets)))
-            self.times = (self.utcs +
-                          offsets.astype('timedelta64[s]')).astype(datetime.datetime)
-            offsets *= -2
-
-            int_offsets = np.floor(offsets).astype(np.int)
-
-            # filling out missing geolocations with computation from pyorbital.
-            line_indices = (self.scans["scan_line_number"]
-                            + int_offsets)
-
-            missed = sorted((set(line_indices) |
-                             set(line_indices + 1))
-                            - set(self.scans["scan_line_number"]))
-
-            min_idx = min(line_indices)
-            max_idx = max(max(line_indices),
-                          max(self.scans["scan_line_number"] - min_idx)) + 1
-            idx_len = max_idx - min_idx + 2
-
-            complete_lons = np.full((idx_len, self.lats.shape[1]), np.nan,
-                                    dtype=np.float)
-            complete_lats = np.full((idx_len, self.lats.shape[1]), np.nan,
-                                    dtype=np.float)
-
-            complete_lons[self.scans["scan_line_number"] - min_idx] = self.lons
-            complete_lats[self.scans["scan_line_number"] - min_idx] = self.lats
-
-            missed_utcs = ((np.array(missed) - 1) * np.timedelta64(500, "ms")
-                           + self.utcs[0])
-
-            mlons, mlats = self.compute_lonlat(width=self.lats.shape[1],
-                                               utcs=missed_utcs,
-                                               clock_drift_adjust=True)
-
-            complete_lons[missed - min_idx] = mlons
-            complete_lats[missed - min_idx] = mlats
-
-            from pygac.slerp import slerp
-            off = offsets - np.floor(offsets)
-            res = slerp(complete_lons[line_indices - min_idx, :],
-                        complete_lats[line_indices - min_idx, :],
-                        complete_lons[line_indices - min_idx + 1, :],
-                        complete_lats[line_indices - min_idx + 1, :],
-                        off[:, np.newaxis, np.newaxis])
-
-            self.lons = res[:, :, 0]
-            self.lats = res[:, :, 1]
-            self.utcs += offsets.astype('timedelta64[s]')
-
-        toc = datetime.datetime.now()
-        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
-        return lons, lats
-
-    def get_telemetry(self):
-        number_of_scans = self.scans["telemetry"].shape[0]
-        decode_tele = np.zeros((int(number_of_scans), 105))
-        decode_tele[:, ::3] = (self.scans["telemetry"] >> 20) & 1023
-        decode_tele[:, 1::3] = (self.scans["telemetry"] >> 10) & 1023
-        decode_tele[:, 2::3] = self.scans["telemetry"] & 1023
-
-        prt_counts = np.mean(decode_tele[:, 17:20], axis=1)
-
-        # getting ICT counts
-
-        ict_counts = np.zeros((int(number_of_scans), 3))
-        ict_counts[:, 0] = np.mean(decode_tele[:, 22:50:3], axis=1)
-        ict_counts[:, 1] = np.mean(decode_tele[:, 23:51:3], axis=1)
-        ict_counts[:, 2] = np.mean(decode_tele[:, 24:52:3], axis=1)
-
-        # getting space counts
-
-        space_counts = np.zeros((int(number_of_scans), 3))
-        space_counts[:, 0] = np.mean(decode_tele[:, 54:100:5], axis=1)
-        space_counts[:, 1] = np.mean(decode_tele[:, 55:101:5], axis=1)
-        space_counts[:, 2] = np.mean(decode_tele[:, 56:102:5], axis=1)
-
-        return prt_counts, ict_counts, space_counts
-
-    def _get_corrupt_mask(self):
-        """Get mask for corrupt scanlines."""
-        mask = ((self.scans["quality_indicators"] >> 31) |
-                ((self.scans["quality_indicators"] << 4) >> 31) |
-                ((self.scans["quality_indicators"] << 5) >> 31))
-        return mask.astype(bool)
-
-    def get_qual_flags(self):
-        """Read quality flags."""
-        number_of_scans = self.scans["telemetry"].shape[0]
-        qual_flags = np.zeros((int(number_of_scans), 7))
-        qual_flags[:, 0] = self.scans["scan_line_number"]
-        qual_flags[:, 1] = (self.scans["quality_indicators"] >> 31)
-        qual_flags[:, 2] = ((self.scans["quality_indicators"] << 4) >> 31)
-        qual_flags[:, 3] = ((self.scans["quality_indicators"] << 5) >> 31)
-        qual_flags[:, 4] = ((self.scans["quality_indicators"] << 13) >> 31)
-        qual_flags[:, 5] = ((self.scans["quality_indicators"] << 14) >> 31)
-        qual_flags[:, 6] = ((self.scans["quality_indicators"] << 15) >> 31)
-
-        return qual_flags
-
-    def postproc(self, channels):
-        """No POD specific postprocessing to be done."""
-        pass
-
-    def get_tsm_pixels(self, channels):
-        """Determine pixels affected by the scan motor issue.
-
-        Uses channels 1, 2, 4 and 5. Neither 3a, nor 3b.
-        """
-        return get_tsm_idx(channels[:, :, 0], channels[:, :, 1],
-                           channels[:, :, 3], channels[:, :, 4])
+    def __init__(self, *args, **kwargs):
+        """Init the GAC POD reader."""
+        GACReader.__init__(self, *args, **kwargs)
+        self.scanline_type = scanline
+        self.offset = 3220
+        self.scan_points = np.arange(3.5, 2048, 5)
 
 
 def main(filename, start_line, end_line):
-    tic = datetime.datetime.now()
-    reader = GACPODReader()
-    reader.read(filename)
-    reader.get_lonlat()
-    channels = reader.get_calibrated_channels()
-    sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles()
-    qual_flags = reader.get_qual_flags()
-    if np.all(reader.mask):
-        print("ERROR: All data is masked out. Stop processing")
-        raise ValueError("All data is masked out.")
-
-    gac_io.save_gac(reader.spacecraft_name,
-                    reader.utcs,
-                    reader.lats, reader.lons,
-                    channels[:, :, 0], channels[:, :, 1],
-                    np.full_like(channels[:, :, 0], np.nan),
-                    channels[:, :, 2],
-                    channels[:, :, 3],
-                    channels[:, :, 4],
-                    sun_zen, sat_zen, sun_azi, sat_azi, rel_azi,
-                    qual_flags, start_line, end_line,
-                    reader.filename,
-                    reader.get_midnight_scanline(),
-                    reader.get_miss_lines())
-    LOG.info("pygac took: %s", str(datetime.datetime.now() - tic))
+    """Generate a l1c file."""
+    return main_pod(GACPODReader, filename, start_line, end_line)
 
 
 if __name__ == "__main__":

pygac/lac_klm.py

+#!/usr/bin/python
+# Copyright (c) 2014-2019
+#
+
+# Author(s):
+
+#   Abhay Devasthale <abhay.devasthale@smhi.se>
+#   Sajid Pareeth <sajid.pareeth@fmach.it>
+#   Martin Raspaud <martin.raspaud@smhi.se>
+#   Adam Dybbroe <adam.dybbroe@smhi.se>
+
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Reader for LAC KLM data."""
+
+import logging
+
+import numpy as np
+
+from pygac.klm_reader import KLMReader, main_klm
+from pygac.lac_reader import LACReader
+
+LOG = logging.getLogger(__name__)
+
+
+# video data object
+
+scanline = np.dtype([("scan_line_number", ">u2"),
+                     ("scan_line_year", ">u2"),
+                     ("scan_line_day_of_year", ">u2"),
+                     ("satellite_clock_drift_delta", ">i2"),
+                     ("scan_line_utc_time_of_day", ">u4"),
+                     ("scan_line_bit_field", ">u2"),
+                     ("zero_fill0", ">i2", (5, )),
+                     # QUALITY INDICATORS
+                     ("quality_indicator_bit_field", ">u4"),
+                     ("scan_line_quality_flags", [("reserved", ">u1"),
+                                                  ("time_problem_code", ">u1"),
+                                                  ("calibration_problem_code", ">u1"),
+                                                  ("earth_location_problem_code", ">u1")]),
+                     ("calibration_quality_flags", ">u2", (3, )),
+                     ("count_of_bit_errors_in_frame_sync", ">u2"),
+                     ("zero_fill1", ">i4", (2, )),
+                     # CALIBRATION COEFFICIENTS
+                     ("visible_operational_cal_ch_1_slope_1", ">i4"),
+                     ("visible_operational_cal_ch_1_intercept_1", ">i4"),
+                     ("visible_operational_cal_ch_1_slope_2", ">i4"),
+                     ("visible_operational_cal_ch_1_intercept_2", ">i4"),
+                     ("visible_operational_cal_ch_1_intersection", ">i4"),
+                     ("visible_test_cal_ch_1_slope_1", ">i4"),
+                     ("visible_test_cal_ch_1_intercept_1", ">i4"),
+                     ("visible_test_cal_ch_1_slope_2", ">i4"),
+                     ("visible_test_cal_ch_1_intercept_2", ">i4"),
+                     ("visible_test_cal_ch_1_intersection", ">i4"),
+                     ("visible_prelaunch_cal_ch_1_slope_1", ">i4"),
+                     ("visible_prelaunch_cal_ch_1_intercept_1", ">i4"),
+                     ("visible_prelaunch_cal_ch_1_slope_2", ">i4"),
+                     ("visible_prelaunch_cal_ch_1_intercept_2", ">i4"),
+                     ("visible_prelaunch_cal_ch_1_intersection", ">i4"),
+                     ("visible_operational_cal_ch_2_slope_1", ">i4"),
+                     ("visible_operational_cal_ch_2_intercept_1", ">i4"),
+                     ("visible_operational_cal_ch_2_slope_2", ">i4"),
+                     ("visible_operational_cal_ch_2_intercept_2", ">i4"),
+                     ("visible_operational_cal_ch_2_intersection", ">i4"),
+                     ("visible_test_cal_ch_2_slope_1", ">i4"),
+                     ("visible_test_cal_ch_2_intercept_1", ">i4"),
+                     ("visible_test_cal_ch_2_slope_2", ">i4"),
+                     ("visible_test_cal_ch_2_intercept_2", ">i4"),
+                     ("visible_test_cal_ch_2_intersection", ">i4"),
+                     ("visible_prelaunch_cal_ch_2_slope_1", ">i4"),
+                     ("visible_prelaunch_cal_ch_2_intercept_1", ">i4"),
+                     ("visible_prelaunch_cal_ch_2_slope_2", ">i4"),
+                     ("visible_prelaunch_cal_ch_2_intercept_2", ">i4"),
+                     ("visible_prelaunch_cal_ch_2_intersection", ">i4"),
+                     ("visible_operational_cal_ch_3a_slope_1", ">i4"),
+                     ("visible_operational_cal_ch_3a_intercept_1", ">i4"),
+                     ("visible_operational_cal_ch_3a_slope_2", ">i4"),
+                     ("visible_operational_cal_ch_3a_intercept_2", ">i4"),
+                     ("visible_operational_cal_ch_3a_intersection", ">i4"),
+                     ("visible_test_cal_ch_3a_slope_1", ">i4"),
+                     ("visible_test_cal_ch_3a_intercept_1", ">i4"),
+                     ("visible_test_cal_ch_3a_slope_2", ">i4"),
+                     ("visible_test_cal_ch_3a_intercept_2", ">i4"),
+                     ("visible_test_cal_ch_3a_intersection", ">i4"),
+                     ("visible_prelaunch_cal_ch_3a_slope_1", ">i4"),
+                     ("visible_prelaunch_cal_ch_3a_intercept_1", ">i4"),
+                     ("visible_prelaunch_cal_ch_3a_slope_2", ">i4"),
+                     ("visible_prelaunch_cal_ch_3a_intercept_2", ">i4"),
+                     ("visible_prelaunch_cal_ch_3a_intersection", ">i4"),
+                     ("ir_operational_cal_ch_3b_coefficient_1", ">i4"),
+                     ("ir_operational_cal_ch_3b_coefficient_2", ">i4"),
+                     ("ir_operational_cal_ch_3b_coefficient_3", ">i4"),
+                     ("ir_test_cal_ch_3b_coefficient_1", ">i4"),
+                     ("ir_test_cal_ch_3b_coefficient_2", ">i4"),
+                     ("ir_test_cal_ch_3b_coefficient_3", ">i4"),
+                     ("ir_operational_cal_ch_4_coefficient_1", ">i4"),
+                     ("ir_operational_cal_ch_4_coefficient_2", ">i4"),
+                     ("ir_operational_cal_ch_4_coefficient_3", ">i4"),
+                     ("ir_test_cal_ch_4_coefficient_1", ">i4"),
+                     ("ir_test_cal_ch_4_coefficient_2", ">i4"),
+                     ("ir_test_cal_ch_4_coefficient_3", ">i4"),
+                     ("ir_operational_cal_ch_5_coefficient_1", ">i4"),
+                     ("ir_operational_cal_ch_5_coefficient_2", ">i4"),
+                     ("ir_operational_cal_ch_5_coefficient_3", ">i4"),
+                     ("ir_test_cal_ch_5_coefficient_1", ">i4"),
+                     ("ir_test_cal_ch_5_coefficient_2", ">i4"),
+                     ("ir_test_cal_ch_5_coefficient_3", ">i4"),
+                     # NAVIGATION
+                     ("computed_yaw_steering", ">i2", (3,)),  # only in version 5
+                     ("total_applied_attitude_correction",
+                      ">i2", (3,)),  # only in version 5
+                     ("navigation_status_bit_field", ">u4"),
+                     ("time_associated_with_tip_euler_angles", ">u4"),
+                     ("tip_euler_angles", ">i2", (3, )),
+                     ("spacecraft_altitude_above_reference_ellipsoid", ">u2"),
+                     ("angular_relationships", ">i2", (153, )),
+                     ("zero_fill2", ">i2", (3, )),
+                     ("earth_location", ">i4", (102, )),
+                     ("zero_fill3", ">i4", (2, )),
+                     # HRPT MINOR FRAME TELEMETRY
+                     ("frame_sync", ">u2", (6, )),
+                     ("id", ">u2", (2, )),
+                     ("time_code", ">u2", (4, )),
+                     ('telemetry', [("ramp_calibration", '>u2', (5, )),
+                                    ("PRT", '>u2', (3, )),
+                                    ("ch3_patch_temp", '>u2'),
+                                    ("spare", '>u2'), ]),
+                     ("back_scan", ">u2", (30, )),
+                     ("space_data", ">u2", (50, )),
+                     ("sync_delta", ">u2"),
+                     ("zero_fill4", ">i2"),
+                     # EARTH OBSERVATIONS
+                     ("sensor_data", ">u4", (3414,)),
+                     ("zero_fill5", ">i4", (2,)),
+                     # DIGITAL B HOUSEKEEPING TELEMETRY
+                     ("digital_b_telemetry_update_flags", ">u2"),
+                     ("avhrr_digital_b_data", ">u2"),
+                     ("zero_fill6", ">i4", (3,)),
+                     # ANALOG HOUSEKEEPING DATA (TIP)
+                     ("analog_telemetry_update_flags", ">u4"),
+                     ("patch_temperature_range", ">u1"),
+                     ("patch_temperature_extended", ">u1"),
+                     ("patch_power", ">u1"),
+                     ("radiator_temperature", ">u1"),
+                     ("blackbody_temperature_1", ">u1"),
+                     ("blackbody_temperature_2", ">u1"),
+                     ("blackbody_temperature_3", ">u1"),
+                     ("blackbody_temperature_4", ">u1"),
+                     ("electronics_current", ">u1"),
+                     ("motor_current", ">u1"),
+                     ("earth_shield_position", ">u1"),
+                     ("electronics_temperature", ">u1"),
+                     ("cooler_housing_temperature", ">u1"),
+                     ("baseplate_temperature", ">u1"),
+                     ("motor_housing_temperature", ">u1"),
+                     ("a/d_converter_temperature", ">u1"),
+                     ("detector_#4_bias_voltage", ">u1"),
+                     ("detector_#5_bias_voltage", ">u1"),
+                     ("blackbody_temperature_channel3b", ">u1"),
+                     ("blackbody_temperature_channel4", ">u1"),
+                     ("blackbody_temperature_channel5", ">u1"),
+                     ("reference_voltage", ">u1"),
+                     ("zero_fill7", ">i2", (3,)),
+                     # CLOUDS FROM AVHRR (CLAVR)
+                     ("reserved_clavr_status_bit_field", ">u4"),
+                     ("reserved_clavr", ">u4"),
+                     ("reserved_clavr_ccm", ">u2", (256,)),
+                     # FILLER
+                     ("zero_fill8", ">i4", (94,))])
+
+
+class LACKLMReader(LACReader, KLMReader):
+    """The LAC KLM reader.
+
+    The offset attribute tells where in the file the scanline data starts.
+    """
+
+    def __init__(self, *args, **kwargs):
+        """Init the LAC KLM reader."""
+        LACReader.__init__(self, *args, **kwargs)
+        self.scanline_type = scanline
+        self.offset = 15872
+        # packed: 15872
+        # self.offset = 22528
+        # unpacked: 22528
+
+
+def main(filename, start_line, end_line):
+    """Generate a l1c file."""
+    return main_klm(LACKLMReader, filename, start_line, end_line)
+
+
+if __name__ == "__main__":
+    import sys
+    main(sys.argv[1], sys.argv[2], sys.argv[3])

pygac/lac_pod.py

+#!/usr/bin/python
+# Copyright (c) 2014-2019.
+#
+
+# Author(s):
+
+#   Abhay Devasthale <abhay.devasthale@smhi.se>
+#   Adam Dybbroe <adam.dybbroe@smhi.se>
+#   Sajid Pareeth <sajid.pareeth@fmach.it>
+#   Martin Raspaud <martin.raspaud@smhi.se>
+
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""The LAC POD reader routines."""
+
+import logging
+
+import numpy as np
+
+from pygac.pod_reader import PODReader, main_pod
+from pygac.lac_reader import LACReader
+
+LOG = logging.getLogger(__name__)
+
+scanline = np.dtype([("scan_line_number", ">i2"),
+                     ("time_code", ">u2", (3, )),
+                     ("quality_indicators", ">u4"),
+                     ("calibration_coefficients", ">i4", (10, )),
+                     ("number_of_meaningful_zenith_angles_and_earth_location_appended",
+                      ">u1"),
+                     ("solar_zenith_angles", "i1", (51, )),
+                     ("earth_location", ">i2", (102, )),
+                     ("telemetry", ">u4", (35, )),
+                     ("sensor_data", ">u4", (3414, )),
+                     ("add_on_zenith", ">u2", (10, )),
+                     ("clock_drift_delta", ">u2"),  # only in newest version
+                     ("spare3", "u2", (337, ))])
+
+
+class LACPODReader(LACReader, PODReader):
+    """The LAC POD reader.
+
+    The `scan_points` attributes provides the position of the longitude and latitude points to
+    compute relative to the full swath width.
+
+    The offset attribute tells where in the file the scanline data starts.
+    """
+
+    def __init__(self, *args, **kwargs):
+        """Init the LAC POD reader."""
+        LACReader.__init__(self, *args, **kwargs)
+        self.scanline_type = scanline
+        self.offset = 14800
+        self.scan_points = np.arange(2048)
+
+
+def main(filename, start_line, end_line):
+    """Generate a l1c file."""
+    return main_pod(LACPODReader, filename, start_line, end_line)
+
+
+if __name__ == "__main__":
+    import sys
+    main(sys.argv[1], sys.argv[2], sys.argv[3])

pygac/lac_reader.py

+#!/usr/bin/python
+# Copyright (c) 2014-2019.
+#
+
+# Author(s):
+
+#   Abhay Devasthale <abhay.devasthale@smhi.se>
+#   Adam Dybbroe <adam.dybbroe@smhi.se>
+#   Sajid Pareeth <sajid.pareeth@fmach.it>
+#   Martin Raspaud <martin.raspaud@smhi.se>
+
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""The LAC reader."""
+
+from pygac.reader import Reader
+import pygac.pygac_geotiepoints as gtp
+
+
+class LACReader(Reader):
+    """Reader for LAC data."""
+
+    # Scanning frequency (scanlines per millisecond)
+    scan_freq = 6.0 / 1000.0
+
+    def __init__(self, *args, **kwargs):
+        """Init the LAC reader."""
+        super(LACReader, self).__init__(*args, **kwargs)
+        self.scan_width = 2048
+        self.lonlat_interpolator = gtp.lac_lat_lon_interpolator

pygac/pod_reader.py

+#!/usr/bin/python
+# Copyright (c) 2014-2019 Pygac developers
+#
+
+# Author(s):
+
+#   Abhay Devasthale <abhay.devasthale@smhi.se>
+#   Adam Dybbroe <adam.dybbroe@smhi.se>
+#   Sajid Pareeth <sajid.pareeth@fmach.it>
+#   Martin Raspaud <martin.raspaud@smhi.se>
+
+# This work was done in the framework of ESA-CCI-Clouds phase I
+
+
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+"""POD file reading.
+
+Format specification can be found here:
+http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c2/sec2-0.htm
+http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c3/sec3-1.htm
+"""
+
+import datetime
+import logging
+
+import numpy as np
+
+from pygac.correct_tsm_issue import TSM_AFFECTED_INTERVALS_POD, get_tsm_idx
+from pygac.reader import Reader
+
+LOG = logging.getLogger(__name__)
+
+# common header
+header0 = np.dtype([("noaa_spacecraft_identification_code", ">u1"),
+                    ("data_type_code", ">u1"),
+                    ("start_time", ">u2", (3, )),
+                    ("number_of_scans", ">u2"),
+                    ("end_time", ">u2", (3, ))])
+
+
+# For L1B data before September 8, 1992
+# http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/k/app-k.htm
+header1 = np.dtype([("noaa_spacecraft_identification_code", ">u1"),
+                    ("data_type_code", ">u1"),
+                    ("start_time", ">u2", (3, )),
+                    ("number_of_scans", ">u2"),
+                    ("end_time", ">u2", (3, )),
+                    ("processing_block_id", "S7"),
+                    ("ramp_auto_calibration", ">u1"),
+                    ("number_of_data_gaps", ">u2"),
+                    ("dacs_quality", ">u1", (6, )),
+                    ("calibration_parameter_id", ">i2"),
+                    ("dacs_status", ">u1"),
+                    ("spare1", ">i1", (5, )),
+                    ("data_set_name", "S44")])
+
+# For L1B data between October 21, 1992 to November 15, 1994
+# http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/l/app-l.htm
+header2 = np.dtype([("noaa_spacecraft_identification_code", ">u1"),
+                    ("data_type_code", ">u1"),
+                    ("start_time", ">u2", (3, )),
+                    ("number_of_scans", ">u2"),
+                    ("end_time", ">u2", (3, )),
+                    ("processing_block_id", "S7"),
+                    ("ramp_auto_calibration", ">u1"),
+                    ("number_of_data_gaps", ">u2"),
+                    ("dacs_quality", ">u1", (6, )),
+                    ("calibration_parameter_id", ">i2"),
+                    ("dacs_status", ">u1"),
+                    ("spare1", ">i1", (5, )),
+                    ("data_set_name", "S42"),
+                    ("blankfill", "S2"),
+                    ("julian_year_of_epoch", ">u2"),
+                    ("julian_day_of_epoch", ">u2"),
+                    ("millisecond_utc_epoch_time_of_day", ">u4"),
+                    # Keplerian orbital elements
+                    ("semi_major_axis", ">f8"),
+                    ("eccentricity", ">f8"),
+                    ("inclination", ">f8"),
+                    ("argument_of_perigee", ">f8"),
+                    ("right_ascension", ">f8"),
+                    ("mean_anomaly", ">f8"),
+                    # cartesian inertial true date of elements
+                    ("x_component_of_position_vector", ">f8"),
+                    ("y_component_of_position_vector", ">f8"),
+                    ("z_component_of_position_vector", ">f8"),
+                    ("x_dot_component_of_position_vector", ">f8"),
+                    ("y_dot_component_of_position_vector", ">f8"),
+                    ("z_dot_component_of_position_vector", ">f8")])
+
+# For L1B data post November 15, 1994
+# http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c2/sec2-0.htm
+header3 = np.dtype([("noaa_spacecraft_identification_code", ">u1"),
+                    ("data_type_code", ">u1"),
+                    ("start_time", ">u2", (3, )),
+                    ("number_of_scans", ">u2"),
+                    ("end_time", ">u2", (3, )),
+                    ("processing_block_id", "S7"),
+                    ("ramp_auto_calibration", ">u1"),
+                    ("number_of_data_gaps", ">u2"),
+                    ("dacs_quality", ">u1", (6, )),
+                    ("calibration_parameter_id", ">i2"),
+                    ("dacs_status", ">u1"),
+                    ("reserved_for_mounting_and_fixed_attitude_correction_indicator",
+                     ">i1"),
+                    ("nadir_earth_location_tolerance", ">i1"),
+                    ("spare1", ">i1"),
+                    ("start_of_data_set_year", ">u2"),
+                    ("data_set_name", "S44"),
+                    ("year_of_epoch", ">u2"),
+                    ("julian_day_of_epoch", ">u2"),
+                    ("millisecond_utc_epoch_time_of_day", ">u4"),
+                    # Keplerian orbital elements
+                    ("semi_major_axis", ">i4"),
+                    ("eccentricity", ">i4"),
+                    ("inclination", ">i4"),
+                    ("argument_of_perigee", ">i4"),
+                    ("right_ascension", ">i4"),
+                    ("mean_anomaly", ">i4"),
+                    # cartesian inertial true date of elements
+                    ("x_component_of_position_vector", ">i4"),
+                    ("y_component_of_position_vector", ">i4"),
+                    ("z_component_of_position_vector", ">i4"),
+                    ("x_dot_component_of_position_vector", ">i4"),
+                    ("y_dot_component_of_position_vector", ">i4"),
+                    ("z_dot_component_of_position_vector", ">i4"),
+                    # future use
+                    ("yaw_fixed_error_correction", ">i2"),
+                    ("roll_fixed_error_correction", ">i2"),
+                    ("pitch_fixed_error_correction", ">i2")])
+
+# archive header
+tbm_header = np.dtype([('fill', 'S30'),
+                       ('data_set_name', 'S44'),
+                       ('select_flag', 'S1'),
+                       ('beginning_latitude', 'S3'),
+                       ('ending_latitude', 'S3'),
+                       ('beginning_longitude', 'S4'),
+                       ('ending_longitude', 'S4'),
+                       ('start_hour', 'S2'),
+                       ('start_minute', 'S2'),
+                       ('number_of_minutes', 'S3'),
+                       ('appended_data_flag', 'S1'),
+                       ('channel_select_flag', 'S1', (20, )),
+                       ('sensor_data_word_size', 'S2'),
+                       ('fill2', 'S3')])
+
+
+class PODReader(Reader):
+    """The POD reader."""
+
+    spacecrafts_orbital = {25: 'tiros n',
+                           2: 'noaa 6',
+                           4: 'noaa 7',
+                           6: 'noaa 8',
+                           7: 'noaa 9',
+                           8: 'noaa 10',
+                           1: 'noaa 11',
+                           5: 'noaa 12',
+                           3: 'noaa 14',
+                           }
+    spacecraft_names = {25: 'tirosn',
+                        2: 'noaa6',
+                        4: 'noaa7',
+                        6: 'noaa8',
+                        7: 'noaa9',
+                        8: 'noaa10',
+                        1: 'noaa11',
+                        5: 'noaa12',
+                        3: 'noaa14',
+                        }
+
+    tsm_affected_intervals = TSM_AFFECTED_INTERVALS_POD
+
+    def correct_scan_line_numbers(self):
+        """Correct the scan line numbers."""
+        # Perform common corrections first.
+        super(PODReader, self).correct_scan_line_numbers()
+
+        # cleaning up the data
+        min_scanline_number = np.amin(
+            np.absolute(self.scans["scan_line_number"][:]))
+        if self.scans["scan_line_number"][0] == self.scans["scan_line_number"][-1] + 1:
+            while self.scans["scan_line_number"][0] != min_scanline_number:
+                self.scans = np.roll(self.scans, -1)
+        else:
+            while self.scans["scan_line_number"][0] != min_scanline_number:
+                self.scans = self.scans[1:]
+
+        self.scans = self.scans[self.scans["scan_line_number"] != 0]
+
+    def read(self, filename):
+        """Read the data.
+
+        Args:
+            filename: string
+                The filename to read from.
+
+        Returns:
+            header: numpy record array
+                The header metadata
+            scans: numpy record array
+                The scanlines
+
+        """
+        super(PODReader, self).read(filename=filename)
+        # choose the right header depending on the date
+        with open(filename) as fd_:
+            # read archive header
+            self.tbm_head = np.fromfile(fd_, dtype=tbm_header, count=1)[0]
+            if ((not self.tbm_head['data_set_name'].startswith(self.creation_site + b'.')) and
+                    (self.tbm_head['data_set_name'] != b'\x00' * 42 + b'  ')):
+                fd_.seek(0)
+                self.tbm_head = None
+                tbm_offset = 0
+            else:
+                tbm_offset = tbm_header.itemsize
+
+            head = np.fromfile(fd_, dtype=header0, count=1)[0]
+            year, jday, _ = self.decode_timestamps(head["start_time"])
+
+            start_date = (datetime.date(year, 1, 1) +
+                          datetime.timedelta(days=int(jday) - 1))
+
+            if start_date < datetime.date(1992, 9, 8):
+                header = header1
+            elif start_date <= datetime.date(1994, 11, 15):
+                header = header2
+            else:
+                header = header3
+
+            fd_.seek(tbm_offset, 0)
+            self.head = np.fromfile(fd_, dtype=header, count=1)[0]
+            fd_.seek(self.offset + tbm_offset, 0)
+            self.scans = np.fromfile(fd_,
+                                     dtype=self.scanline_type,
+                                     count=self.head["number_of_scans"])
+
+        self.correct_scan_line_numbers()
+        self.spacecraft_id = self.head["noaa_spacecraft_identification_code"]
+        if self.spacecraft_id == 1 and start_date < datetime.date(1982, 1, 1):
+            self.spacecraft_id = 25
+        self.spacecraft_name = self.spacecraft_names[self.spacecraft_id]
+        LOG.info(
+            "Reading %s data", self.spacecrafts_orbital[self.spacecraft_id])
+
+        return self.head, self.scans
+
+    def get_header_timestamp(self):
+        """Get the timestamp from the header.
+
+        Returns:
+            A datetime object containing the timestamp from the header.
+
+        Raises:
+            A ValueError if the timestamp is corrupt.
+
+        """
+        year, jday, msec = self.decode_timestamps(self.head["start_time"])
+        try:
+            return self.to_datetime(self.to_datetime64(year=year, jday=jday,
+                                                       msec=msec))
+        except ValueError as err:
+            raise ValueError('Corrupt header timestamp: {0}'.format(err))
+
+    @staticmethod
+    def decode_timestamps(encoded):
+        """Decode timestamps.
+
+        Returns:
+            year
+            day of year
+            milliseconds since 00:00
+
+        """
+        ndims = len(encoded.shape)
+        if ndims == 1:
+            # Single header timestamp
+            enc0 = encoded[0]
+            enc1 = encoded[1]
+            enc2 = encoded[2]
+        elif ndims == 2:
+            # Scanline timestamps
+            enc0 = encoded[:, 0]
+            enc1 = encoded[:, 1]
+            enc2 = encoded[:, 2]
+        else:
+            raise ValueError('Invalid timestamp dimension')
+
+        year = enc0 >> 9
+        year = np.where(year > 75, year + 1900, year + 2000)
+        jday = (enc0 & 0x1FF)
+        msec = ((np.uint32(enc1 & 2047) << 16) | (np.uint32(enc2)))
+
+        return year, jday, msec
+
+    def _get_times(self):
+        return self.decode_timestamps(self.scans["time_code"])
+
+    def _adjust_clock_drift(self):
+        """Adjust the geolocation to compensate for the clock error.
+
+        TODO: bad things might happen when scanlines are skipped.
+        """
+        tic = datetime.datetime.now()
+        self.get_times()
+        from pygac.clock_offsets_converter import get_offsets
+        try:
+            offset_times, clock_error = get_offsets(self.spacecraft_name)
+        except KeyError:
+            LOG.info("No clock drift info available for %s",
+                     self.spacecraft_name)
+        else:
+            offset_times = np.array(offset_times, dtype='datetime64[ms]')
+            offsets = np.interp(self.utcs.astype(np.uint64),
+                                offset_times.astype(np.uint64),
+                                clock_error)
+            LOG.info("Adjusting for clock drift of %s to %s",
+                     str(min(offsets)),
+                     str(max(offsets)))
+            self.times = (self.utcs +
+                          offsets.astype('timedelta64[s]')).astype(datetime.datetime)
+            offsets *= -2
+
+            int_offsets = np.floor(offsets).astype(np.int)
+
+            # filling out missing geolocations with computation from pyorbital.
+            line_indices = (self.scans["scan_line_number"]
+                            + int_offsets)
+
+            missed = sorted((set(line_indices) |
+                             set(line_indices + 1))
+                            - set(self.scans["scan_line_number"]))
+
+            min_idx = min(line_indices)
+            max_idx = max(max(line_indices),
+                          max(self.scans["scan_line_number"] - min_idx)) + 1
+            idx_len = max_idx - min_idx + 2
+
+            complete_lons = np.full((idx_len, self.lats.shape[1]), np.nan,
+                                    dtype=np.float)
+            complete_lats = np.full((idx_len, self.lats.shape[1]), np.nan,
+                                    dtype=np.float)
+
+            complete_lons[self.scans["scan_line_number"] - min_idx] = self.lons
+            complete_lats[self.scans["scan_line_number"] - min_idx] = self.lats
+
+            missed_utcs = ((np.array(missed) - 1) * np.timedelta64(500, "ms")
+                           + self.utcs[0])
+            try:
+                mlons, mlats = self.compute_lonlat(width=self.lats.shape[1],
+                                                   utcs=missed_utcs,
+                                                   clock_drift_adjust=True)
+            except IndexError as err:
+                LOG.warning(
+                    'Cannot perform clock drift correction: %s', str(err))
+                return
+
+            complete_lons[missed - min_idx] = mlons
+            complete_lats[missed - min_idx] = mlats
+
+            from pygac.slerp import slerp
+            off = offsets - np.floor(offsets)
+            res = slerp(complete_lons[line_indices - min_idx, :],
+                        complete_lats[line_indices - min_idx, :],
+                        complete_lons[line_indices - min_idx + 1, :],
+                        complete_lats[line_indices - min_idx + 1, :],
+                        off[:, np.newaxis, np.newaxis])
+
+            self.lons = res[:, :, 0]
+            self.lats = res[:, :, 1]
+            self.utcs += offsets.astype('timedelta64[s]')
+
+        toc = datetime.datetime.now()
+        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
+        return lons, lats
+
+    def get_telemetry(self):
+        """Get the telemetry.
+
+        Returns:
+            prt_counts: np.array
+            ict_counts: np.array
+            space_counts: np.array
+
+        """
+        number_of_scans = self.scans["telemetry"].shape[0]
+        decode_tele = np.zeros((int(number_of_scans), 105))
+        decode_tele[:, ::3] = (self.scans["telemetry"] >> 20) & 1023
+        decode_tele[:, 1::3] = (self.scans["telemetry"] >> 10) & 1023
+        decode_tele[:, 2::3] = self.scans["telemetry"] & 1023
+
+        prt_counts = np.mean(decode_tele[:, 17:20], axis=1)
+
+        # getting ICT counts
+
+        ict_counts = np.zeros((int(number_of_scans), 3))
+        ict_counts[:, 0] = np.mean(decode_tele[:, 22:50:3], axis=1)
+        ict_counts[:, 1] = np.mean(decode_tele[:, 23:51:3], axis=1)
+        ict_counts[:, 2] = np.mean(decode_tele[:, 24:52:3], axis=1)
+
+        # getting space counts
+
+        space_counts = np.zeros((int(number_of_scans), 3))
+        space_counts[:, 0] = np.mean(decode_tele[:, 54:100:5], axis=1)
+        space_counts[:, 1] = np.mean(decode_tele[:, 55:101:5], axis=1)
+        space_counts[:, 2] = np.mean(decode_tele[:, 56:102:5], axis=1)
+
+        return prt_counts, ict_counts, space_counts
+
+    def _get_corrupt_mask(self):
+        """Get mask for corrupt scanlines."""
+        mask = ((self.scans["quality_indicators"] >> 31) |
+                ((self.scans["quality_indicators"] << 4) >> 31) |
+                ((self.scans["quality_indicators"] << 5) >> 31))
+        return mask.astype(bool)
+
+    def get_qual_flags(self):
+        """Read quality flags."""
+        number_of_scans = self.scans["telemetry"].shape[0]
+        qual_flags = np.zeros((int(number_of_scans), 7))
+        qual_flags[:, 0] = self.scans["scan_line_number"]
+        qual_flags[:, 1] = (self.scans["quality_indicators"] >> 31)
+        qual_flags[:, 2] = ((self.scans["quality_indicators"] << 4) >> 31)
+        qual_flags[:, 3] = ((self.scans["quality_indicators"] << 5) >> 31)
+        qual_flags[:, 4] = ((self.scans["quality_indicators"] << 13) >> 31)
+        qual_flags[:, 5] = ((self.scans["quality_indicators"] << 14) >> 31)
+        qual_flags[:, 6] = ((self.scans["quality_indicators"] << 15) >> 31)
+
+        return qual_flags
+
+    def postproc(self, channels):
+        """No POD specific postprocessing to be done."""
+        pass
+
+    def get_tsm_pixels(self, channels):
+        """Determine pixels affected by the scan motor issue.
+
+        Uses channels 1, 2, 4 and 5. Neither 3a, nor 3b.
+        """
+        return get_tsm_idx(channels[:, :, 0], channels[:, :, 1],
+                           channels[:, :, 3], channels[:, :, 4])
+
+
+def main_pod(reader_cls, filename, start_line, end_line):
+    """Generate a l1c file."""
+    from pygac import gac_io
+    tic = datetime.datetime.now()
+    reader = reader_cls()
+    reader.read(filename)
+    reader.get_lonlat()
+    channels = reader.get_calibrated_channels()
+    sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles()
+
+    qual_flags = reader.get_qual_flags()
+    if (np.all(reader.mask)):
+        print("ERROR: All data is masked out. Stop processing")
+        raise ValueError("All data is masked out.")
+    gac_io.save_gac(reader.spacecraft_name,
+                    reader.utcs,
+                    reader.lats, reader.lons,
+                    channels[:, :, 0], channels[:, :, 1],
+                    np.full_like(channels[:, :, 0], np.nan),
+                    channels[:, :, 2],
+                    channels[:, :, 3],
+                    channels[:, :, 4],
+                    sun_zen, sat_zen, sun_azi, sat_azi, rel_azi,
+                    qual_flags, start_line, end_line,
+                    reader.filename,
+                    reader.meta_data)
+    LOG.info("pygac took: %s", str(datetime.datetime.now() - tic))

pygac/pygac_geotiepoints.py

@@ -29,7 +29,7 @@ import geotiepoints as gtp
 import numpy as np
 
 
-def Gac_Lat_Lon_Interpolator(lons_subset, lats_subset):
+def gac_lat_lon_interpolator(lons_subset, lats_subset):
     """Interpolate lat-lon values in the AVHRR GAC data.
 
     Each GAC row has 409 pixels.
@@ -37,10 +37,20 @@ def Gac_Lat_Lon_Interpolator(lons_subset, lats_subset):
     ranging from 5 to 405. Interpolate to full resolution.
 
     """
-    # cols_subset = np.arange(0, 404, 8)
-    # cols_full = np.arange(405)
     cols_subset = np.arange(4, 405, 8)
     cols_full = np.arange(409)
+    return lat_lon_interpolator(lons_subset, lats_subset, cols_subset, 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)
+
+
+def lat_lon_interpolator(lons_subset, lats_subset, cols_subset, cols_full):
+    """Interpolate lat-lon values in the AVHRR data."""
     lines = lats_subset.shape[0]
     rows_subset = np.arange(lines)
     rows_full = np.arange(lines)