.bumpversion.cfg

 [bumpversion]
-current_version = 1.1.0
+current_version = 1.2.0
 commit = True
 tag = True
 

.stickler.yml

+linters:
+  flake8:
+    fixer: true
+fixers:
+  enable: true

.travis.yml

@@ -19,5 +19,5 @@ deploy:
     secure: FznD38SHxuj1RoFXUqLstF//O3+pypF84hCOHO8A3Poa+ygh7X4a+9aimCeiuY9d+5tbE0ZlyG7LMcCOFA3JCvH7lUChmQbbn3pcFvJdJzzYQgxYNfvCL9YtjRF/n648bVdpN265hm07rmOe7DHbysw4q8hMlxtUH87MYctRg90=
   on:
     tags: true
-    repo: adybbroe/pygac
+    repo: pytroll/pygac
     branch: master

CHANGELOG.md

+## Version 1.2.0 (2019/10/17)
+
+### Issues Closed
+
+* [Issue 7](https://github.com/pytroll/pygac/issues/7) - Project URL points to wrong domain
+
+In this release 1 issue was closed.
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 30](https://github.com/pytroll/pygac/pull/30) - Feature update angles computation
+
+#### Features added
+
+* [PR 34](https://github.com/pytroll/pygac/pull/34) - Use the geotiepoints library
+* [PR 31](https://github.com/pytroll/pygac/pull/31) - Refactor I/O
+
+In this release 3 pull requests were closed.
+
 ## Version 1.1.0 (2019/06/12)
 
 ### Issues Closed

RELEASING.md

+# Releasing Pygac
+
+1. checkout master
+2. pull from repo
+3. run the unittests
+4. run `loghub`.  Replace <github username> and <previous version> with proper
+   values.  To get the previous version run `git tag` and select the most
+   recent with highest version number.
+
+```
+loghub pytroll/pygac -u <github username> -st v<previous version> -plg bug "Bugs fixed" -plg enhancement "Features added" -plg documentation "Documentation changes" -plg backwards-incompatibility "Backwards incompatible changes"
+```
+
+This command will create a CHANGELOG.temp file which need to be added
+to the top of the CHANGLOG.md file.  The same content is also printed
+to terminal, so that can be copy-pasted, too.  Remember to update also
+the version number to the same given in step 5. Don't forget to commit
+CHANGELOG.md!
+
+5. Bump up the version with
+
+```
+bumpversion <level>
+```
+
+where level is one of `major`, `minor`, or `patch`.
+
+
+6. push changes to github `git push --follow-tags`
+7. Verify travis tests passed
+8. Deploy sdist (and wheel) to PyPI

debian/changelog

-pygac (1.1.0-4) UNRELEASED; urgency=medium
+pygac (1.2.0-1) unstable; urgency=medium
 
+  * New upstream release.
   * Bump Standards-Version to 4.4.1, no changes.
-  * Explicitly set Rules-Requires-Root: no.
+  * debian/control:
+    - explicitly set Rules-Requires-Root: no
+    - python3-geotiepoints is nor a dependency
+      (before it was only recommended)
+  * debian/patches:
+    - refresh all patches
+  * Remove obsolete fields Name from debian/upstream/metadata.
 
- -- Antonio Valentino <antonio.valentino@tiscali.it>  Mon, 30 Sep 2019 19:45:55 +0200
+ -- Antonio Valentino <antonio.valentino@tiscali.it>  Sun, 20 Oct 2019 06:50:06 +0000
 
 pygac (1.1.0-3) unstable; urgency=medium
 

debian/control

@@ -23,15 +23,15 @@ Homepage: https://github.com/pytroll/pygac
 
 Package: python3-pygac
 Architecture: any
-Depends: python3-h5py,
+Depends: python3-geotiepoints,
+         python3-h5py,
          python3-numpy,
          python3-pyorbital,
          python3-scipy,
          ${python3:Depends},
          ${shlibs:Depends},
          ${misc:Depends}
-Recommends: python3-geotiepoints,
-            ${python3:Recommends}
+Recommends: ${python3:Recommends}
 Suggests: ${python3:Suggests}
 Description: Python package to read and calibrate NOAA AVHRR GAC data
  pygac reads NOAA AVHRR Global Area Coverage (GAC) data, and performs

debian/patches/0001-Temporary-disable-broken-tests.patch

@@ -9,7 +9,7 @@ See also upstream issue https://github.com/pytroll/pygac/issues/16.
  2 files changed, 2 insertions(+)
 
 diff --git a/pygac/tests/test_calibrate_klm.py b/pygac/tests/test_calibrate_klm.py
-index 3ea36e8..7ed9015 100644
+index 8a489ce..1f2abd3 100644
 --- a/pygac/tests/test_calibrate_klm.py
 +++ b/pygac/tests/test_calibrate_klm.py
 @@ -33,6 +33,7 @@ from pygac.calibration import calibrate_solar, calibrate_thermal
@@ -21,7 +21,7 @@ index 3ea36e8..7ed9015 100644
  
          counts = np.array([[0, 0, 0, 0, 0,
 diff --git a/pygac/tests/test_calibrate_pod.py b/pygac/tests/test_calibrate_pod.py
-index 6027e35..fcf4cdb 100644
+index fe1bffc..aa35d11 100644
 --- a/pygac/tests/test_calibrate_pod.py
 +++ b/pygac/tests/test_calibrate_pod.py
 @@ -33,6 +33,7 @@ from pygac.calibration import calibrate_solar, calibrate_thermal

debian/patches/0002-Fix-config.patch

@@ -4,12 +4,11 @@ Subject: Fix config
 
 ---
  pygac/__init__.py | 2 +-
- pygac/gac_io.py   | 2 +-
  setup.py          | 4 ++--
- 3 files changed, 4 insertions(+), 4 deletions(-)
+ 2 files changed, 3 insertions(+), 3 deletions(-)
 
 diff --git a/pygac/__init__.py b/pygac/__init__.py
-index f7fee41..4425074 100644
+index 9225981..a30e971 100644
 --- a/pygac/__init__.py
 +++ b/pygac/__init__.py
 @@ -30,7 +30,7 @@ try:
@@ -20,27 +19,14 @@ index f7fee41..4425074 100644
 +    CONFIG_FILE = '/etc/pygac.cfg'
  
  if not os.path.exists(CONFIG_FILE) or not os.path.isfile(CONFIG_FILE):
-     LOG.warning(str(CONFIG_FILE) + " pointed to by the environment "
-diff --git a/pygac/gac_io.py b/pygac/gac_io.py
-index a2f0cf4..d613718 100644
---- a/pygac/gac_io.py
-+++ b/pygac/gac_io.py
-@@ -48,7 +48,7 @@ try:
-     CONFIG_FILE = os.environ['PYGAC_CONFIG_FILE']
- except KeyError:
-     LOG.exception('Environment variable PYGAC_CONFIG_FILE not set!')
--    raise
-+    CONFIG_FILE = '/etc/pygac.cfg'
- 
- if not os.path.exists(CONFIG_FILE) or not os.path.isfile(CONFIG_FILE):
-     raise IOError(str(CONFIG_FILE) + " pointed to by the environment " +
+     LOG.warning(
 diff --git a/setup.py b/setup.py
-index 803d379..2d08902 100644
+index cefd25e..5ffbda7 100644
 --- a/setup.py
 +++ b/setup.py
-@@ -105,8 +105,8 @@ if __name__ == '__main__':
-           extras_require={'geolocation interpolation': ['python-geotiepoints'],
-                           },
+@@ -104,8 +104,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')],
 -          data_files=[('etc', ['etc/pygac.cfg.template']),
 -                      ('gapfilled_tles', ['gapfilled_tles/TLE_noaa16.txt'])],

debian/patches/0003-fix-import-mock.patch

@@ -8,7 +8,7 @@ Subject: fix import mock
  2 files changed, 2 insertions(+), 2 deletions(-)
 
 diff --git a/pygac/tests/test_io.py b/pygac/tests/test_io.py
-index db278bb..69f96a3 100644
+index e3cea78..dcd462d 100644
 --- a/pygac/tests/test_io.py
 +++ b/pygac/tests/test_io.py
 @@ -19,7 +19,7 @@
@@ -18,10 +18,10 @@ index db278bb..69f96a3 100644
 -import mock
 +from unittest import mock
  import numpy as np
- from pygac.gac_io import update_start_end_line
- 
+ import numpy.testing
+ import pygac.gac_io as gac_io
 diff --git a/pygac/tests/test_reader.py b/pygac/tests/test_reader.py
-index 7574017..868385a 100644
+index 6fea775..be9447a 100644
 --- a/pygac/tests/test_reader.py
 +++ b/pygac/tests/test_reader.py
 @@ -20,7 +20,7 @@
@@ -31,5 +31,5 @@ index 7574017..868385a 100644
 -import mock
 +from unittest import mock
  import numpy as np
+ import numpy.testing
  from pygac.gac_reader import GACReader
- 

debian/upstream/metadata

----
 Bug-Database: https://github.com/pytroll/pygac/issues
 Bug-Submit: https://github.com/pytroll/pygac/issues/new
-Name: pygac
 Repository: https://github.com/pytroll/pygac.git
 Repository-Browse: https://github.com/pytroll/pygac

doc/source/api.rst

@@ -66,6 +66,15 @@ Each GAC row has total 409 pixels. But lat-lon values are provided for every eig
 
 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
 --------------

pygac/__init__.py

@@ -33,5 +33,22 @@ except KeyError:
     CONFIG_FILE = ''
 
 if not os.path.exists(CONFIG_FILE) or not os.path.isfile(CONFIG_FILE):
-    LOG.warning(str(CONFIG_FILE) + " pointed to by the environment "
+    LOG.warning(
+        str(CONFIG_FILE) + " pointed to by the environment "
         + "variable PYGAC_CONFIG_FILE is not a file or does not exist!")
+
+
+def get_absolute_azimuth_angle_diff(sat_azi, sun_azi):
+    """Calculates absolute azimuth difference angle. """
+    rel_azi = abs(sat_azi - sun_azi)
+    rel_azi = rel_azi % 360
+    # Not using np.where to avoid copying array
+    rel_azi[rel_azi > 180] = 360.0 - rel_azi[rel_azi > 180]
+    return rel_azi
+
+
+def centered_modulus(array, divisor):
+    """Transform array to half open range ]-divisor/2, divisor/2]."""
+    arr = array % divisor
+    arr[arr > divisor / 2] -= divisor
+    return arr

pygac/calibration.py

@@ -439,12 +439,17 @@ def calibrate_solar(counts, chan, year, jday, spacecraft, corr=1):
     sth = (cal.ah[chan] * (100.0 + cal.bh[chan] * t
                            + cal.ch[chan] * t * t)) / 100.0
     if cal.c_s is not None:
-        return np.where(counts <= cal.c_s[chan],
+        refl = np.where(counts <= cal.c_s[chan],
                         (counts - cal.c_dark[chan]) * stl * corr,
                         ((cal.c_s[chan] - cal.c_dark[chan]) * stl
                          + (counts - cal.c_s[chan]) * sth) * corr)
     else:
-        return (counts - cal.c_dark[chan]) * stl * corr
+        refl = (counts - cal.c_dark[chan]) * stl * corr
+
+    # Mask negative reflectances
+    refl[refl < 0] = np.nan
+
+    return refl
 
 
 def calibrate_thermal(counts, prt, ict, space, line_numbers, channel, spacecraft):
@@ -555,4 +560,7 @@ def calibrate_thermal(counts, prt, ict, space, line_numbers, channel, spacecraft
     if chan == 0:
         bt = np.where((counts - new_space) >= 0, 0.0, bt)
 
+    # Mask values outside valid range
+    bt = np.where(np.logical_or(bt < 170.0, bt > 350.0), np.nan, bt)
+
     return bt

pygac/gac_io.py

@@ -33,26 +33,27 @@ import time
 import h5py
 import numpy as np
 
-from .correct_tsm_issue import flag_pixels as flag_tsm_pixels
-
 try:
     import ConfigParser
 except ImportError:
     import configparser as ConfigParser
 
+from pygac import CONFIG_FILE
+from pygac.utils import slice_channel, strip_invalid_lat, check_user_scanlines
 
 LOG = logging.getLogger(__name__)
 
 
-try:
-    CONFIG_FILE = os.environ['PYGAC_CONFIG_FILE']
-except KeyError:
-    LOG.exception('Environment variable PYGAC_CONFIG_FILE not set!')
-    raise
+MISSING_DATA = -32001
+MISSING_DATA_LATLON = -999999
+
 
+def read_config():
+    """Read output dir etc from config file."""
     if not os.path.exists(CONFIG_FILE) or not os.path.isfile(CONFIG_FILE):
-    raise IOError(str(CONFIG_FILE) + " pointed to by the environment " +
-                  "variable PYGAC_CONFIG_FILE is not a file or does not exist!")
+        raise IOError('{} pointed to by the environment variable '
+                      'PYGAC_CONFIG_FILE is not a file or does not exist!'
+                      .format(str(CONFIG_FILE)))
 
     conf = ConfigParser.ConfigParser()
     try:
@@ -71,20 +72,8 @@ OUTPUT_FILE_PREFIX = options['output_file_prefix']
     SUNSATANGLES_DIR = os.environ.get('SM_SUNSATANGLES_DIR', OUTDIR)
     AVHRR_DIR = os.environ.get('SM_AVHRR_DIR', OUTDIR)
     QUAL_DIR = os.environ.get('SM_AVHRR_DIR', OUTDIR)
-MISSING_DATA = -32001
-MISSING_DATA_LATLON = -999999
-
-
-def update_start_end_line(start_line, end_line, temp_start_line, temp_end_line):
-    """Update user start/end lines after data has been sliced using temporary
-    start/end lines.
 
-    Returns:
-        Updated start_line, updated end_line
-    """
-    new_start_line = max(0, start_line - temp_start_line)
-    new_end_line = min(end_line, temp_end_line) - temp_start_line
-    return new_start_line, new_end_line
+    return OUTPUT_FILE_PREFIX, SUNSATANGLES_DIR, AVHRR_DIR, QUAL_DIR
 
 
 def save_gac(satellite_name,
@@ -93,100 +82,119 @@ def save_gac(satellite_name,
              ref1, ref2, ref3,
              bt3, bt4, bt5,
              sun_zen, sat_zen, sun_azi, sat_azi, rel_azi,
-             mask, qual_flags, start_line, end_line, tsmcorr,
-             gac_file, midnight_scanline, miss_lines, switch=None):
+             qual_flags, start_line, end_line,
+             gac_file, midnight_scanline, miss_lines):
 
-    along_track = lats.shape[0]
     last_scan_line_number = qual_flags[-1, 0]
 
-    # Determine scanline range requested by user
-    start_line = int(start_line)
-    end_line = int(end_line)
-    if end_line == 0:
-        # If the user specifies 0 as the last scanline, process all scanlines
-        end_line = along_track
-
-    bt3 = np.where(np.logical_or(bt3 < 170.0, bt3 > 350.0),
-                   MISSING_DATA, bt3 - 273.15)
-    bt4 = np.where(np.logical_or(bt4 < 170.0, bt4 > 350.0),
-                   MISSING_DATA, bt4 - 273.15)
-    bt5 = np.where(np.logical_or(bt5 < 170.0, bt5 > 350.0),
-                   MISSING_DATA, bt5 - 273.15)
-
-    lats = np.where(np.logical_or(lats < -90.00, lats > 90.00),
-                    MISSING_DATA_LATLON, lats)
-    lons = np.where(np.logical_or(lons < -180.00, lons > 180.00),
-                    MISSING_DATA_LATLON, lons)
-
-    sat_azi -= 180.0
-    rel_azi = abs(rel_azi)
-    rel_azi = 180.0 - rel_azi
-
-    for array in [bt3, bt4, bt5]:
-        array[array != MISSING_DATA] = 100 * array[array != MISSING_DATA]
-        array[mask] = MISSING_DATA
-    for array in [ref1, ref2, ref3,
-                  sun_zen, sat_zen, sun_azi, sat_azi, rel_azi]:
-        array *= 100
-        array[mask] = MISSING_DATA
-    for array in [lats, lons]:
-        array[array != MISSING_DATA_LATLON] = 1000.0 * \
-            array[array != MISSING_DATA_LATLON]
-        array[mask] = MISSING_DATA_LATLON
-
-    for ref in [ref1, ref2, ref3]:
-        ref[ref < 0] = MISSING_DATA
-
-    if switch is not None:
-        ref3[switch == 0] = MISSING_DATA
-        bt3[switch == 1] = MISSING_DATA
-        ref3[switch == 2] = MISSING_DATA
-        bt3[switch == 2] = MISSING_DATA
+    # Strip invalid coordinates
+    first_valid_lat, last_valid_lat = strip_invalid_lat(lats)
+    if first_valid_lat > start_line:
+        LOG.info('New start_line chosen (due to invalid lat/lon '
+                 'info) = ' + str(first_valid_lat))
+    if end_line > last_valid_lat:
+        LOG.info('New end_line chosen (due to invalid lat/lon '
+                 'info) = ' + str(last_valid_lat))
+
+    # Check user-defined scanlines
+    start_line, end_line = check_user_scanlines(
+        start_line=start_line,
+        end_line=end_line,
+        first_valid_lat=first_valid_lat,
+        last_valid_lat=last_valid_lat)
 
-    for array in [ref1, ref2, ref3, bt3, bt4, bt5]:
-        array[np.isnan(array)] = MISSING_DATA
+    # Slice data using new start/end lines
+    _, miss_lines, midnight_scanline = slice_channel(
+        np.zeros(lats.shape),
+        start_line=start_line,
+        end_line=end_line,
+        first_valid_lat=first_valid_lat,
+        last_valid_lat=last_valid_lat,
+        qual_flags=qual_flags,
+        miss_lines=miss_lines,
+        midnight_scanline=midnight_scanline)
 
-    # Choose new temporary start/end lines if lat/lon info is invalid
-    no_wrong_lat = np.where(lats != MISSING_DATA_LATLON)
-    temp_start_line = min(no_wrong_lat[0])
-    temp_end_line = max(no_wrong_lat[0])
-    if temp_start_line > start_line:
-        LOG.info('New start_line chosen (due to invalid lat/lon info) = ' + str(temp_start_line))
-    if end_line > temp_end_line:
-        LOG.info('New end_line chosen (due to invalid lat/lon info) = ' + str(temp_end_line))
-
-    # Slice data using temporary start/end lines
-    ref1 = ref1[temp_start_line:temp_end_line + 1, :].copy()
-    ref2 = ref2[temp_start_line:temp_end_line + 1, :].copy()
-    ref3 = ref3[temp_start_line:temp_end_line + 1, :].copy()
-    bt3 = bt3[temp_start_line:temp_end_line + 1, :].copy()
-    bt4 = bt4[temp_start_line:temp_end_line + 1, :].copy()
-    bt5 = bt5[temp_start_line:temp_end_line + 1, :].copy()
-    sun_zen = sun_zen[temp_start_line:temp_end_line + 1, :].copy()
-    sun_azi = sun_azi[temp_start_line:temp_end_line + 1, :].copy()
-    sat_zen = sat_zen[temp_start_line:temp_end_line + 1, :].copy()
-    sat_azi = sat_azi[temp_start_line:temp_end_line + 1, :].copy()
-    rel_azi = rel_azi[temp_start_line:temp_end_line + 1, :].copy()
-    lats = lats[temp_start_line:temp_end_line + 1, :].copy()
-    lons = lons[temp_start_line:temp_end_line + 1, :].copy()
-    miss_lines = np.sort(np.array(
-        qual_flags[0:temp_start_line, 0].tolist() +
-        miss_lines.tolist() +
-        qual_flags[temp_end_line+1:, 0].tolist()
-    ))
-    qual_flags = qual_flags[temp_start_line:temp_end_line+1, :].copy()
-    xutcs = xutcs[temp_start_line:temp_end_line+1].copy()
-
-    # Update user start/end lines to the new slice
-    start_line, end_line = update_start_end_line(
+    ref1, _, _ = slice_channel(ref1,
+                               start_line=start_line,
+                               end_line=end_line,
+                               first_valid_lat=first_valid_lat,
+                               last_valid_lat=last_valid_lat)
+    ref2, _, _ = slice_channel(ref2,
+                               start_line=start_line,
+                               end_line=end_line,
+                               first_valid_lat=first_valid_lat,
+                               last_valid_lat=last_valid_lat)
+    ref3, _, _ = slice_channel(ref3,
+                               start_line=start_line,
+                               end_line=end_line,
+                               first_valid_lat=first_valid_lat,
+                               last_valid_lat=last_valid_lat)
+    bt3, _, _ = slice_channel(bt3,
+                              start_line=start_line,
+                              end_line=end_line,
+                              first_valid_lat=first_valid_lat,
+                              last_valid_lat=last_valid_lat)
+    bt4, _, _ = slice_channel(bt4,
                               start_line=start_line,
                               end_line=end_line,
-        temp_start_line=temp_start_line,
-        temp_end_line=temp_end_line)
+                              first_valid_lat=first_valid_lat,
+                              last_valid_lat=last_valid_lat)
+    bt5, _, _ = slice_channel(bt5,
+                              start_line=start_line,
+                              end_line=end_line,
+                              first_valid_lat=first_valid_lat,
+                              last_valid_lat=last_valid_lat)
+    sun_zen, _, _ = slice_channel(sun_zen,
+                                  start_line=start_line,
+                                  end_line=end_line,
+                                  first_valid_lat=first_valid_lat,
+                                  last_valid_lat=last_valid_lat)
+    sun_azi, _, _ = slice_channel(sun_azi,
+                                  start_line=start_line,
+                                  end_line=end_line,
+                                  first_valid_lat=first_valid_lat,
+                                  last_valid_lat=last_valid_lat)
+    sat_zen, _, _ = slice_channel(sat_zen,
+                                  start_line=start_line,
+                                  end_line=end_line,
+                                  first_valid_lat=first_valid_lat,
+                                  last_valid_lat=last_valid_lat)
+    sat_azi, _, _ = slice_channel(sat_azi,
+                                  start_line=start_line,
+                                  end_line=end_line,
+                                  first_valid_lat=first_valid_lat,
+                                  last_valid_lat=last_valid_lat)
+    rel_azi, _, _ = slice_channel(rel_azi,
+                                  start_line=start_line,
+                                  end_line=end_line,
+                                  first_valid_lat=first_valid_lat,
+                                  last_valid_lat=last_valid_lat)
+    lons, _, _ = slice_channel(lons,
+                               start_line=start_line,
+                               end_line=end_line,
+                               first_valid_lat=first_valid_lat,
+                               last_valid_lat=last_valid_lat)
+    lats, _, _ = slice_channel(lats,
+                               start_line=start_line,
+                               end_line=end_line,
+                               first_valid_lat=first_valid_lat,
+                               last_valid_lat=last_valid_lat)
+    qual_flags, _, _ = slice_channel(qual_flags,
+                                     start_line=start_line,
+                                     end_line=end_line,
+                                     first_valid_lat=first_valid_lat,
+                                     last_valid_lat=last_valid_lat)
+    xutcs, _, _ = slice_channel(xutcs,
+                                start_line=start_line,
+                                end_line=end_line,
+                                first_valid_lat=first_valid_lat,
+                                last_valid_lat=last_valid_lat)
+
+    total_number_of_scan_lines = lats.shape[0]
 
     # Reading time from the body of the gac file
-    start = xutcs[start_line].astype(datetime.datetime)
-    end = xutcs[end_line].astype(datetime.datetime)
+    start = xutcs[0].astype(datetime.datetime)
+    end = xutcs[-1].astype(datetime.datetime)
     startdate = start.strftime("%Y%m%d")
     starttime = start.strftime("%H%M%S%f")[:-5]
     enddate = end.strftime("%Y%m%d")
@@ -196,53 +204,22 @@ def save_gac(satellite_name,
     # Earth-Sun distance correction factor
     corr = 1.0 - 0.0334 * np.cos(2.0 * np.pi * (jday - 2) / 365.25)
 
-    # Slice scanline range requested by user
-    ref1 = ref1[start_line:end_line+1, :].copy()
-    ref2 = ref2[start_line:end_line+1, :].copy()
-    ref3 = ref3[start_line:end_line+1, :].copy()
-    bt3 = bt3[start_line:end_line+1, :].copy()
-    bt4 = bt4[start_line:end_line+1, :].copy()
-    bt5 = bt5[start_line:end_line+1, :].copy()
-    sun_zen = sun_zen[start_line:end_line+1, :].copy()
-    sun_azi = sun_azi[start_line:end_line+1, :].copy()
-    sat_zen = sat_zen[start_line:end_line+1, :].copy()
-    sat_azi = sat_azi[start_line:end_line+1, :].copy()
-    rel_azi = rel_azi[start_line:end_line+1, :].copy()
-    lats = lats[start_line:end_line+1, :].copy()
-    lons = lons[start_line:end_line+1, :].copy()
-    qual_flags = qual_flags[start_line:end_line+1, :].copy()
-    xutcs = xutcs[start_line:end_line+1].copy()
-
-    if midnight_scanline is not None:
-        # Update midnight scanline to the final scanline range
-        midnight_scanline -= (temp_start_line + start_line)
-
-        # Set midnight scanline to None if it has been removed due to invalid
-        # lat/lon info (< 0) or lies outside the user defined scanline range
-        if midnight_scanline < 0 or (midnight_scanline > end_line or
-                                     midnight_scanline < start_line):
-            midnight_scanline = None
-
-    # Compute total number of scanlines
-    total_number_of_scan_lines = end_line - start_line + 1
-
-    # Correct for temporary scan motor issue.
-    #
-    # TODO: The thresholds in tsm.flag_pixels() were derived from the final
-    #       pygac output, that's why the correction is applied here. It would
-    #       certainly be more consistent to apply the correction in GACReader,
-    #       but that requires a new threshold analysis.
-    if tsmcorr:
-        LOG.info('Correcting for temporary scan motor issue')
-        tic = datetime.datetime.now()
-        (ref1, ref2, bt3, bt4, bt5, ref3) = flag_tsm_pixels(channel1=ref1,
-                                                            channel2=ref2,
-                                                            channel3b=bt3,
-                                                            channel4=bt4,
-                                                            channel5=bt5,
-                                                            channel3a=ref3,
-                                                            fillv=MISSING_DATA)
-        LOG.debug('TSM correction took: %s', str(datetime.datetime.now() - tic))
+    # Apply scaling & offset
+    bt3 -= 273.15
+    bt4 -= 273.15
+    bt5 -= 273.15
+    for array in [bt3, bt4, bt5, ref1, ref2, ref3, sun_zen, sat_zen, sun_azi,
+                  sat_azi, rel_azi]:
+        array *= 100.0
+    for array in [lats, lons]:
+        array *= 1000.0
+
+    # Replace NaN with fill values
+    for array in [ref1, ref2, ref3, bt3, bt4, bt5, sun_zen, sat_zen, sun_azi,
+                  sat_azi, rel_azi]:
+        array[np.isnan(array)] = MISSING_DATA
+    for array in [lats, lons]:
+        array[np.isnan(array)] = MISSING_DATA_LATLON
 
     avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime,
                 lats, lons, ref1, ref2, ref3, bt3, bt4, bt5,
@@ -260,6 +237,9 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime,
                 miss_lines):
     import os
 
+    # Read output dir etc from config file
+    OUTPUT_FILE_PREFIX, SUNSATANGLES_DIR, AVHRR_DIR, QUAL_DIR = read_config()
+
     # Calculate start and end time in sec1970
     t_obj = time.strptime(startdate + starttime[0:6], "%Y%m%d%H%M%S")
     starttime_sec1970 = calendar.timegm(t_obj)

pygac/gac_klm.py

@@ -32,9 +32,8 @@ http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c8/sec83142-1.htm
 from __future__ import print_function
 
 import numpy as np
-from .correct_tsm_issue import TSM_AFFECTED_INTERVALS_KLM
+from pygac.correct_tsm_issue import TSM_AFFECTED_INTERVALS_KLM, get_tsm_idx
 from pygac.gac_reader import GACReader, inherit_doc
-import pygac.geotiepoints as gtp
 import datetime
 from pygac import gac_io
 import logging
@@ -522,12 +521,10 @@ class GACKLMReader(GACReader):
 
         return prt_counts, ict_counts, space_counts
 
-    def get_lonlat(self):
-        arr_lat = self.scans["earth_location"][:, 0::2] / 1e4
-        arr_lon = self.scans["earth_location"][:, 1::2] / 1e4
-
-        self.lons, self.lats = gtp.Gac_Lat_Lon_Interpolator(arr_lon, arr_lat)
-        return self.lons, self.lats
+    def _get_lonlat(self):
+        lats = self.scans["earth_location"][:, 0::2] / 1e4
+        lons = self.scans["earth_location"][:, 1::2] / 1e4
+        return lons, lats
 
     def get_header_timestamp(self):
         year = self.head['start_of_data_set_year']
@@ -546,29 +543,62 @@ class GACKLMReader(GACReader):
         return year, jday, msec
 
     def get_ch3_switch(self):
-        """0 for 3b, 1 for 3a
+        """Channel 3 identification.
+
+        0: Channel 3b (Brightness temperature
+        1: Channel 3a (Reflectance)
+        2: Transition (No data)
         """
         return self.scans["scan_line_bit_field"][:] & 3
 
-    def get_corrupt_mask(self):
-
-        # corrupt scanlines
-
+    def _get_corrupt_mask(self):
+        """Get mask for corrupt scanlines."""
         mask = ((self.scans["quality_indicator_bit_field"] >> 31) |
                 ((self.scans["quality_indicator_bit_field"] << 3) >> 31) |
                 ((self.scans["quality_indicator_bit_field"] << 4) >> 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_indicator_bit_field"] >> 31)
-        qual_flags[:, 2] = ((self.scans["quality_indicator_bit_field"] << 3) >> 31)
-        qual_flags[:, 3] = ((self.scans["quality_indicator_bit_field"] << 4) >> 31)
-        qual_flags[:, 4] = ((self.scans["quality_indicator_bit_field"] << 24) >> 30)
-        qual_flags[:, 5] = ((self.scans["quality_indicator_bit_field"] << 26) >> 30)
-        qual_flags[:, 6] = ((self.scans["quality_indicator_bit_field"] << 28) >> 30)
+        qual_flags[:, 2] = (
+            (self.scans["quality_indicator_bit_field"] << 3) >> 31)
+        qual_flags[:, 3] = (
+            (self.scans["quality_indicator_bit_field"] << 4) >> 31)
+        qual_flags[:, 4] = (
+            (self.scans["quality_indicator_bit_field"] << 24) >> 30)
+        qual_flags[:, 5] = (
+            (self.scans["quality_indicator_bit_field"] << 26) >> 30)
+        qual_flags[:, 6] = (
+            (self.scans["quality_indicator_bit_field"] << 28) >> 30)
+
+        return qual_flags
+
+    def postproc(self, channels):
+        """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
 
-        return mask.astype(bool), qual_flags
+    def _adjust_clock_drift(self):
+        """Clock drift correction is only applied to POD satellites."""
+        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[:, :, 4], channels[:, :, 5])
 
 
 def main(filename, start_line, end_line):
@@ -578,9 +608,8 @@ def main(filename, start_line, end_line):
     reader.get_lonlat()
     channels = reader.get_calibrated_channels()
     sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles()
-
-    mask, qual_flags = reader.get_corrupt_mask()
-    if (np.all(mask)):
+    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.")
 
@@ -593,12 +622,10 @@ def main(filename, start_line, end_line):
                     channels[:, :, 4],
                     channels[:, :, 5],
                     sun_zen, sat_zen, sun_azi, sat_azi, rel_azi,
-                    mask, qual_flags, start_line, end_line,
-                    reader.is_tsm_affected(),
+                    qual_flags, start_line, end_line,
                     reader.filename,
                     reader.get_midnight_scanline(),
-                    reader.get_miss_lines(),
-                    reader.get_ch3_switch())
+                    reader.get_miss_lines())
     LOG.info("pygac took: %s", str(datetime.datetime.now() - tic))
 
 

pygac/gac_pod.py

@@ -39,8 +39,7 @@ from __future__ import print_function
 
 import numpy as np
 from pygac.gac_reader import GACReader, inherit_doc
-import pygac.geotiepoints as gtp
-from .correct_tsm_issue import TSM_AFFECTED_INTERVALS_POD
+from pygac.correct_tsm_issue import TSM_AFFECTED_INTERVALS_POD, get_tsm_idx
 import datetime
 from pygac import gac_io
 
@@ -186,9 +185,9 @@ class GACPODReader(GACReader):
 
     tsm_affected_intervals = TSM_AFFECTED_INTERVALS_POD
 
-    def correct_scan_line_numbers(self, plot=False):
+    def correct_scan_line_numbers(self):
         # Perform common corrections first.
-        super(GACPODReader, self).correct_scan_line_numbers(plot=plot)
+        super(GACPODReader, self).correct_scan_line_numbers()
 
         # cleaning up the data
         min_scanline_number = np.amin(np.absolute(self.scans["scan_line_number"][:]))
@@ -271,7 +270,7 @@ class GACPODReader(GACReader):
     def _get_times(self):
         return self.decode_timestamps(self.scans["time_code"])
 
-    def adjust_clock_drift(self):
+    def _adjust_clock_drift(self):
         """Adjust the geolocation to compensate for the clock error.
 
         TODO: bad things might happen when scanlines are skipped.
@@ -311,8 +310,10 @@ class GACPODReader(GACReader):
                           max(self.scans["scan_line_number"] - min_idx)) + 1
             idx_len = max_idx - min_idx + 2
 
-            complete_lons = np.zeros((idx_len, 409), dtype=np.float) * np.nan
-            complete_lats = np.zeros((idx_len, 409), dtype=np.float) * np.nan
+            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
@@ -320,7 +321,9 @@ class GACPODReader(GACReader):
             missed_utcs = ((np.array(missed) - 1) * np.timedelta64(500, "ms")
                            + self.utcs[0])
 
-            mlons, mlats = self.compute_lonlat(missed_utcs, True)
+            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
@@ -340,13 +343,10 @@ class GACPODReader(GACReader):
         toc = datetime.datetime.now()
         LOG.debug("clock drift adjustment took %s", str(toc - tic))
 
-    def get_lonlat(self):
-        # interpolating lat-on points using PYTROLL geotiepoints
-        arr_lat = self.scans["earth_location"][:, 0::2] / 128.0
-        arr_lon = self.scans["earth_location"][:, 1::2] / 128.0
-
-        self.lons, self.lats = gtp.Gac_Lat_Lon_Interpolator(arr_lon, arr_lat)
-        return self.lons, self.lats
+    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]
@@ -373,14 +373,15 @@ class GACPODReader(GACReader):
 
         return prt_counts, ict_counts, space_counts
 
-    def get_corrupt_mask(self):
-
-        # corrupt scanlines
-
+    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"]
@@ -391,7 +392,19 @@ class GACPODReader(GACReader):
         qual_flags[:, 5] = ((self.scans["quality_indicators"] << 14) >> 31)
         qual_flags[:, 6] = ((self.scans["quality_indicators"] << 15) >> 31)
 
-        return mask.astype(bool), qual_flags
+        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(filename, start_line, end_line):
@@ -399,12 +412,10 @@ def main(filename, start_line, end_line):
     reader = GACPODReader()
     reader.read(filename)
     reader.get_lonlat()
-    reader.adjust_clock_drift()
     channels = reader.get_calibrated_channels()
     sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles()
-
-    mask, qual_flags = reader.get_corrupt_mask()
-    if (np.all(mask)):
+    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.")
 
@@ -412,13 +423,12 @@ def main(filename, start_line, end_line):
                     reader.utcs,
                     reader.lats, reader.lons,
                     channels[:, :, 0], channels[:, :, 1],
-                    np.ones_like(channels[:, :, 0]) * -1,
+                    np.full_like(channels[:, :, 0], np.nan),
                     channels[:, :, 2],
                     channels[:, :, 3],
                     channels[:, :, 4],
                     sun_zen, sat_zen, sun_azi, sat_azi, rel_azi,
-                    mask, qual_flags, start_line, end_line,
-                    reader.is_tsm_affected(),
+                    qual_flags, start_line, end_line,
                     reader.filename,
                     reader.get_midnight_scanline(),
                     reader.get_miss_lines())

pygac/gac_reader.py

@@ -23,22 +23,25 @@
 """Generic reader for GAC data. Can't be used as is, has to be subclassed to add
 specific read functions.
 """
+from abc import ABCMeta, abstractmethod, abstractproperty
+import logging
 import numpy as np
-from pygac import CONFIG_FILE
+import os
+import six
+import types
+
+from pygac import (CONFIG_FILE, centered_modulus,
+                   get_absolute_azimuth_angle_diff)
 try:
     import ConfigParser
 except ImportError:
     import configparser as ConfigParser
 
-import os
-import logging
 from pyorbital.orbital import Orbital
 from pyorbital import astronomy
 import datetime
 from pygac.calibration import calibrate_solar, calibrate_thermal
-from abc import ABCMeta, abstractmethod, abstractproperty
-import six
-import types
+import pygac.geotiepoints as gtp
 
 LOG = logging.getLogger(__name__)
 
@@ -48,7 +51,23 @@ class GACReader(six.with_metaclass(ABCMeta)):
     scan_freq = 2.0/1000.0
     """Scanning frequency (scanlines per millisecond)"""
 
-    def __init__(self, tle_thresh=7):
+    def __init__(self, interpolate_coords=True, adjust_clock_drift=True,
+                 tle_dir=None, tle_name=None, tle_thresh=7):
+        """
+        Args:
+            interpolate_coords: Interpolate coordinates from every eighth pixel
+                to all pixels.
+            adjust_clock_drift: Adjust the geolocation to compensate for the
+                clock error (POD satellites only).
+            tle_dir: Directory holding TLE files
+            tle_name: Filename pattern of TLE files.
+            tle_thresh: Maximum number of days between observation and nearest
+                TLE
+        """
+        self.interpolate_coords = interpolate_coords
+        self.adjust_clock_drift = adjust_clock_drift
+        self.tle_dir = tle_dir
+        self.tle_name = tle_name
         self.tle_thresh = tle_thresh
         self.head = None
         self.scans = None
@@ -60,6 +79,7 @@ class GACReader(six.with_metaclass(ABCMeta)):
         self.times = None
         self.tle_lines = None
         self.filename = None
+        self._mask = None
 
     @abstractmethod
     def read(self, filename):
@@ -177,7 +197,14 @@ class GACReader(six.with_metaclass(ABCMeta)):
         """
         return (scan_line_number - 1) / self.scan_freq
 
-    def compute_lonlat(self, utcs=None, clock_drift_adjust=True):
+    def compute_lonlat(self, width, utcs=None, clock_drift_adjust=True):
+        """Compute lat/lon coordinates.
+
+        Args:
+            width: Number of coordinates per scanlines
+            utcs: Scanline timestamps
+            clock_drift_adjust: If True, adjust clock drift.
+        """
         tle1, tle2 = self.get_tle_lines()
 
         scan_points = np.arange(3.5, 2048, 5)
@@ -226,7 +253,7 @@ class GACReader(six.with_metaclass(ABCMeta)):
 
         lons, lats = pos_time[:2]
 
-        return lons.reshape(-1, 409), lats.reshape(-1, 409)
+        return lons.reshape(-1, width), lats.reshape(-1, width)
 
     def get_calibrated_channels(self):
         channels = self.get_counts()
@@ -256,8 +283,74 @@ class GACReader(six.with_metaclass(ABCMeta)):
                 self.scans["scan_line_number"],
                 chan,
                 self.spacecraft_name)
+
+        # Mask out corrupt values
+        channels[self.mask] = np.nan
+
+        # Apply KLM/POD specific postprocessing
+        self.postproc(channels)
+
+        # Mask pixels affected by scan motor issue
+        if self.is_tsm_affected():
+            LOG.info('Correcting for temporary scan motor issue')
+            self.mask_tsm_pixels(channels)
+
         return channels
 
+    def get_lonlat(self):
+        """Compute lat/lon coordinates.
+
+        TODO: Switch to faster interpolator?
+        """
+        if self.lons is None and self.lats is None:
+            self.lons, self.lats = self._get_lonlat()
+
+            # Interpolate from every eighth pixel to all pixels.
+            if self.interpolate_coords:
+                self.lons, self.lats = gtp.Gac_Lat_Lon_Interpolator(
+                    self.lons, self.lats)
+
+            # Adjust clock drift
+            if self.adjust_clock_drift:
+                self._adjust_clock_drift()
+
+            # Mask out corrupt scanlines
+            self.lons[self.mask] = np.nan
+            self.lats[self.mask] = np.nan
+
+            # Mask values outside the valid range
+            self.lats[np.fabs(self.lats) > 90.0] = np.nan
+            self.lons[np.fabs(self.lons) > 180.0] = np.nan
+
+        return self.lons, self.lats
+
+    @abstractmethod
+    def _get_lonlat(self):
+        """KLM/POD specific readout of lat/lon coordinates."""
+        raise NotImplementedError
+
+    @property
+    def mask(self):
+        """Mask for corrupt scanlines."""
+        if self._mask is None:
+            self._mask = self._get_corrupt_mask()
+        return self._mask
+
+    @abstractmethod
+    def _get_corrupt_mask(self):
+        """KLM/POD specific readout of corrupt scanline mask."""
+        raise NotImplementedError
+
+    @abstractmethod
+    def postproc(self, channels):
+        """Apply KLM/POD specific postprocessing."""
+        raise NotImplementedError
+
+    @abstractmethod
+    def _adjust_clock_drift(self):
+        """Adjust clock drift."""
+        raise NotImplementedError
+
     @staticmethod
     def tle2datetime64(times):
         """Convert TLE timestamps to numpy.datetime64
@@ -279,6 +372,11 @@ class GACReader(six.with_metaclass(ABCMeta)):
         return times64
 
     def get_tle_file(self):
+        """Find TLE file for the current satellite."""
+        tle_dir, tle_name = self.tle_dir, self.tle_name
+
+        # If user didn't specify TLE dir/name, try config file
+        if tle_dir is None or tle_name is None:
             conf = ConfigParser.ConfigParser()
             try:
                 conf.read(CONFIG_FILE)
@@ -287,13 +385,21 @@ class GACReader(six.with_metaclass(ABCMeta)):
                               str(CONFIG_FILE))
                 raise
 
-        values = {"satname": self.spacecraft_name, }
             options = {}
             for option, value in conf.items('tle', raw=True):
                 options[option] = value
-        tle_filename = os.path.join(options['tledir'],
-                                    options["tlename"] % values)
+
+            tle_dir = options['tledir']
+            tle_name = options['tlename']
+
+        values = {"satname": self.spacecraft_name, }
+        tle_filename = os.path.join(tle_dir, tle_name % values)
         LOG.info('TLE filename = ' + str(tle_filename))
+
+        return tle_filename
+
+    def read_tle_file(self, tle_filename):
+        """Read TLE file."""
         with open(tle_filename, 'r') as fp_:
             return fp_.readlines()
 
@@ -307,7 +413,7 @@ class GACReader(six.with_metaclass(ABCMeta)):
             return self.tle_lines
 
         self.get_times()
-        tle_data = self.get_tle_file()
+        tle_data = self.read_tle_file(self.get_tle_file())
         sdate = np.datetime64(self.times[0], '[ms]')
         dates = self.tle2datetime64(
             np.array([float(line[18:32]) for line in tle_data[::2]]))
@@ -347,7 +453,25 @@ class GACReader(six.with_metaclass(ABCMeta)):
         return tle1, tle2
 
     def get_angles(self):
+        """Get azimuth and zenith angles.
+
+        Azimuth angle definition is the same as in pyorbital, but with
+        different units (degrees not radians for sun azimuth angles)
+        and different ranges.
+
+        Returns:
+            sat_azi: satellite azimuth angle
+                degree clockwise from north in range ]-180, 180],
+            sat_zentih: satellite zenith angles in degrees in range [0,90],
+            sun_azi: sun azimuth angle
+                degree clockwise from north in range ]-180, 180],
+            sun_zentih: sun zenith angles in degrees in range [0,90],
+            rel_azi: absolute azimuth angle difference in degrees between sun
+                and sensor in range [0, 180]
+
+        """
         self.get_times()
+        self.get_lonlat()
         tle1, tle2 = self.get_tle_lines()
         orb = Orbital(self.spacecrafts_orbital[self.spacecraft_id],
                       line1=tle1, line2=tle2)
@@ -364,10 +488,15 @@ class GACReader(six.with_metaclass(ABCMeta)):
                                             self.lons, self.lats)
         del alt
         sun_azi = np.rad2deg(sun_azi)
-        sun_azi = np.where(sun_azi < 0, sun_azi + 180, sun_azi - 180)
+        rel_azi = get_absolute_azimuth_angle_diff(sun_azi, sat_azi)
 
-        rel_azi = abs(sat_azi - sun_azi)
-        rel_azi = np.where(rel_azi > 180.0, 360.0 - rel_azi, rel_azi)
+        # Scale angles range to half open interval ]-180, 180]
+        sat_azi = centered_modulus(sat_azi, 360.0)
+        sun_azi = centered_modulus(sun_azi, 360.0)
+
+        # Mask corrupt scanlines
+        for arr in (sat_azi, sat_zenith, sun_azi, sun_zenith, rel_azi):
+            arr[self.mask] = np.nan
 
         return sat_azi, sat_zenith, sun_azi, sun_zenith, rel_azi
 
@@ -426,7 +555,7 @@ class GACReader(six.with_metaclass(ABCMeta)):
 
         return year, jday, msec
 
-    def correct_scan_line_numbers(self, plot=False):
+    def correct_scan_line_numbers(self):
         """Remove scanlines with corrupted scanline numbers
 
         This includes:
@@ -439,17 +568,12 @@ class GACReader(six.with_metaclass(ABCMeta)):
             - NSS.GHRR.NJ.D96064.S0043.E0236.B0606162.WI
             - NSS.GHRR.NJ.D99286.S1818.E2001.B2466869.WI
 
-        Args:
-            plot (bool): If True, plot results.
+        Returns:
+            Intermediate and final results (for plotting purpose)
         """
         along_track = np.arange(1, len(self.scans["scan_line_number"])+1)
-
-        # Plot original scanline numbers
-        if plot:
-            import matplotlib.pyplot as plt
-            _, (ax0, ax1) = plt.subplots(nrows=2)
-            ax0.plot(along_track, self.scans["scan_line_number"], "b-",
-                     label="original")
+        results = {'along_track': along_track,
+                   'n_orig': self.scans['scan_line_number'].copy()}
 
         # Remove scanlines whose scanline number is outside the valid range
         within_range = np.logical_and(self.scans["scan_line_number"] < 15000,
@@ -495,30 +619,16 @@ class GACReader(six.with_metaclass(ABCMeta)):
         LOG.debug('Removed %s scanline(s) with corrupt scanline numbers',
                   str(len(along_track) - len(self.scans)))
 
-        # Plot corrected scanline numbers
-        if plot:
-            along_track = along_track[within_range]
-            along_track = along_track[diffs <= thresh]
-            ax0.plot(along_track, self.scans["scan_line_number"], "r--",
-                     label="corrected")
-            ax0.set_ylabel("Scanline Number")
-            ax0.set_xlabel("Along Track")
-            ax0.legend(loc="best")
-
-            ax1.plot(np.arange(len(nz_diffs)), nz_diffs)
-            ax1.axhline(thresh, color="r", label="thresh={0:.2f}"
-                        .format(thresh))
-            ax1.set_xlabel("Index")
-            ax1.set_ylabel("nonzero |n - n'|")
-            ax1.legend()
-
-            plt.tight_layout()
-            plt.savefig(self.filename + "_scanline_number_correction.png",
-                        bbox_inches='tight')
+        results.update({'n_corr': self.scans['scan_line_number'],
+                        'within_range': within_range,
+                        'diffs': diffs,
+                        'thresh': thresh,
+                        'nz_diffs': nz_diffs})
+        return results
 
     def correct_times_thresh(self, max_diff_from_t0_head=6*60*1000,
                              min_frac_near_t0_head=0.01,
-                             max_diff_from_ideal_t=10*1000, plot=False):
+                             max_diff_from_ideal_t=10*1000):
         """Correct corrupted timestamps using a threshold approach.
 
         The threshold approach is based on the scanline number and the header
@@ -555,20 +665,20 @@ class GACReader(six.with_metaclass(ABCMeta)):
                 If a scanline timestamp deviates more than max_diff_from_ideal_t
                 from the ideal timestamp, it is regarded as corrupt and will be
                 replaced with the ideal timestamp.
-            plot (bool): If True, plot results.
+        Returns:
+            Intermediate and final results (for plotting purpose)
         """
-        apply_corr = True
-        fail_reason = ""
+        results = {}
         dt64_msec = ">M8[ms]"
 
         # Check whether scanline number increases monotonically
         n = self.scans["scan_line_number"]
+        results.update({'t': self.utcs.copy(), 'n': n})
         if np.any(np.diff(n) < 0):
             LOG.error("Cannot perform timestamp correction. Scanline number "
                       "does not increase monotonically.")
-            apply_corr = False
-            fail_reason = "Scanline number jumps backwards"
-            # We could already return here, but continue for plotting purpose
+            results['fail_reason'] = "Scanline number jumps backwards"
+            return results
 
         # Convert time to milliseconds since 1970-01-01
         t = self.utcs.astype("i8")
@@ -598,60 +708,27 @@ class GACReader(six.with_metaclass(ABCMeta)):
         #    we do not have reliable information and cannot proceed.
         near_t0_head = np.where(
             np.fabs(offsets - t0_head) <= max_diff_from_t0_head)[0]
+        results.update({'offsets': offsets,
+                        't0_head': t0_head,
+                        'max_diff_from_t0_head': max_diff_from_t0_head})
         if near_t0_head.size / float(n.size) >= min_frac_near_t0_head:
             t0 = np.median(offsets[near_t0_head])
         else:
             LOG.error("Timestamp mismatch. Cannot perform correction.")
-            fail_reason = "Timestamp mismatch"
-            apply_corr = False
-            t0 = 0
+            results['fail_reason'] = "Timestamp mismatch"
+            return results
 
         # Add estimated offset to the ideal timestamps
         tn += t0
 
         # Replace timestamps deviating more than a certain threshold from the
         # ideal timestamp with the ideal timestamp.
-        if apply_corr:
         corrupt_lines = np.where(np.fabs(t - tn) > max_diff_from_ideal_t)
         self.utcs[corrupt_lines] = tn[corrupt_lines].astype(dt64_msec)
         LOG.debug("Corrected %s timestamp(s)", str(len(corrupt_lines[0])))
 
-        # Plot results
-        if plot:
-            import matplotlib.pyplot as plt
-            along_track = np.arange(n.size)
-            _, (ax0, ax1, ax2) = plt.subplots(nrows=3, sharex=True,
-                                              figsize=(8, 10))
-
-            ax0.plot(along_track, t, "b-", label="original")
-            if apply_corr:
-                ax0.plot(along_track, self.utcs, color="red", linestyle="--",
-                         label="corrected")
-            else:
-                ax0.set_title(fail_reason)
-            ax0.set_ylabel("Time [ms since 1970]")
-            ax0.set_ylim((self.utcs.min()-np.timedelta64(30, "m")).astype("i8"),
-                         (self.utcs.max()+np.timedelta64(30, "m")).astype("i8"))
-            ax0.legend(loc="best")
-            ax2.plot(along_track, n)
-            ax2.set_ylabel("Scanline number")
-            ax2.set_xlabel("Along Track")
-
-            ax1.plot(along_track, offsets)
-            ax1.fill_between(
-                along_track,
-                t0_head - np.ones(along_track.size)*max_diff_from_t0_head,
-                t0_head + np.ones(along_track.size)*max_diff_from_t0_head,
-                facecolor="g", alpha=0.33)
-            ax1.axhline(y=t0_head, color="g", linestyle="--",
-                        label="Header timestamp")
-            ax1.set_ylim(t0_head-5*max_diff_from_t0_head,
-                         t0_head+5*max_diff_from_t0_head)
-            ax1.set_ylabel("Offset t-tn [ms]")
-            ax1.legend(loc="best")
-
-            plt.savefig(self.filename+"_timestamp_correction.png",
-                        bbox_inches="tight", dpi=100)
+        results.update({'tn': tn, 'tcorr': self.utcs, 't0': t0})
+        return results
 
     @abstractproperty
     def tsm_affected_intervals(self):
@@ -717,6 +794,19 @@ class GACReader(six.with_metaclass(ABCMeta)):
         missing = sorted(ideal.difference(set(self.scans['scan_line_number'])))
         return np.array(missing)
 
+    def mask_tsm_pixels(self, channels):
+        """Mask pixels affected by the scan motor issue."""
+        idx = self.get_tsm_pixels(channels)
+        channels[idx] = np.nan
+
+    @abstractmethod
+    def get_tsm_pixels(self, channels):
+        """Determine pixels affected by the scan motor issue.
+
+        Channel selection is POD/KLM specific.
+        """
+        raise NotImplementedError
+
 
 def inherit_doc(cls):
     """Make a class method inherit its docstring from the parent class.