Repository: pytroll/pygac (language: python)

Issue:
Update documentation
Some ideas:

- Command line and library usage
- Explain config file contents
- Include API of reader classes (so that the satpy reader can refer to this for available `reader_kwargs`)
- [Installation instructions](https://pygac.readthedocs.io/en/develop/installation.html)
    - There's no need for using prefix with virtualenv or conda. Just `pip install pygac` for the latest release or `pip install git+...` for the latest development version
    - Development installation instructions (`git clone && pip install -e .`) would be nice, too


Links to POD/KLM user guides are broken
Couldn't find them on the web. I've asked NOAA NCDC where the pages are located now.


Discussion / hints:
Nice, we got a hacktoberfest tag.
I could also provide some instructions on the calibration coefficient creation and the json layout.
Does pytroll have a style guide for the documentation? Otherwise, I would just try to write it as similar as possible to the existing docs.
Yeah, but now it's over :scream:

Nevertheless, sounds great! ~~I think as long as it's not too much text, we can stick to the README. But at some point we should probably change to sphinx.~~ I'm not aware of a style guide for docs. But the existing pygac docs, plus maybe satpy & pyresample are probably a good reference.
Thanks. I think @abhaydd had a place, but it's not any official website.
http://webapp1.dlib.indiana.edu/virtual_disk_library/index.cgi/4284724/FID1385/podug/index.htm
http://webapp1.dlib.indiana.edu/virtual_disk_library/index.cgi/2790181/FID1497/klm/index.htm
I recommend to get offline PDF copies.  I can't attach one here due to the 10 MB file limit but we don't know if the University of Indiana will keep this one forever.

See #44 
Got a reply from NOAA: https://www1.ncdc.noaa.gov/pub/data/satellite/publications/podguides/
That's where I got the PDFs from.  I decided to download a local copy.  I found that the US government has a habit of shutting down whenever I need to consult those guides.

Fix patch:
diff --git a/README.md b/README.md
index 4f1aafa..4d7905a 100644
--- a/README.md
+++ b/README.md
@@ -5,6 +5,7 @@ pygac
 [![Coverage](https://codecov.io/gh/pytroll/pygac/branch/main/graph/badge.svg?token=DQMgf2LxuM)](https://codecov.io/gh/pytroll/pygac)
 
 
-A python package to read, calibrate and navigate NOAA and Metop AVHRR GAC and LAC data.
-
+Pygac is a Python package to read, calibrate and navigate data from the AVHRR 
+instrument onboard NOAA and MetOp satellites in GAC and LAC format.
 
+The documentation is available at https://pygac.readthedocs.io/.
diff --git a/doc/source/api.rst b/doc/source/api.rst
index 53d53f6..39c292a 100644
--- a/doc/source/api.rst
+++ b/doc/source/api.rst
@@ -1,607 +1,95 @@
-The :mod:`pygac` API
-====================
+The Pygac API
+=============
 
-*pygac* interface consists of a number of python modules. The following schematic shows a general structure and the processing flow of *pygac*.   
+.. inheritance-diagram:: pygac.gac_pod.GACPODReader pygac.gac_klm.GACKLMReader
+                         pygac.lac_pod.LACPODReader pygac.lac_klm.LACKLMReader
 
-It must be noted that *pygac* expects Level 1b file to contain normal GAC header and data records, the format of which are mentioned in the official NOAA POD and KLM Data User Guides. The user should not prepend any other header (e.g. when downloading GAC data from CLASS archive etc) to the L1b file. In the first pre-processing step, *pygac* determines whether the GAC data comes from the second (i.e. NOAA-14 and before) or the third generation (NOAA-15 and onwards) AVHRR instrument by "pygac-run".
-This is done by reading the first three bytes of the data set. If they contain the any of the following values, ["CMS", "NSS", "UKM", "DSS"], then the KLM reader from "gac_klm.py" file is invoked, otherwise the POD reader is invoked (gac_pod.py).
 
+Reader Base Classes
+-------------------
 
-GAC POD reader
---------------
-   
-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.
-
-In case of POD family of satellites, the first data record often doesn't start with first scan line number. In fact, the latter often has unambiguous value and may not be the continuation of the last record number (as expected due to overlap).
-
-Here is an example showing scanlines numbers for the first and the last 3 lines from one random NOAA-11 orbit:
-
-[12804 3 4 ..., 12806 12807 12808]
-You can see that the first two scanlines are missing and the orbit is actually starting with scanline number 3.
-
-*pygac* would rearrange the entire orbit in the increasing order of the scanline numbers.
-
-Sometimes the scanlines have erroneous time information. Here is an example of time stamp from few consecutive scanlines from a radom orbit from NOAA-11::
-
-       [datetime.datetime(2071, 4, 28, 15, 44, 46, 329000)],
-       [datetime.datetime(2067, 4, 27, 15, 44, 46, 336000)],
-       [datetime.datetime(2065, 4, 25, 15, 44, 46, 328000)],
-       [datetime.datetime(2059, 4, 22, 15, 44, 46, 331000)],
-       [datetime.datetime(2057, 4, 22, 15, 27, 17, 757000)],
-       [datetime.datetime(2055, 4, 20, 15, 27, 17, 753000)],
-       [datetime.datetime(2069, 7, 20, 14, 17, 19, 377000)],
-       [datetime.datetime(1981, 7, 25, 14, 17, 12, 379000)],
-       [datetime.datetime(1987, 7, 27, 14, 17, 12, 380000)],
-       [datetime.datetime(1991, 7, 29, 14, 17, 11, 359000)],
-       [datetime.datetime(1993, 7, 29, 14, 17, 11, 360000)],
+Common functionality shared by multiple readers.
 
-It is evident that year, month values are jumping and have non-sense values.
+Reader
+~~~~~~
 
-*pygac* would handle this in either of the following two ways. If the first scanline has valid time stamp, *pygac* would use this to calculate time stamp for remaining scanlines assuming the scanning rate of 0.5 sec/scanline and using actual scanline numbers (to account for any potential gaps in the data). If the first scanline doesn't have valid time stamp or missing altogether, then *pygac* computes the median value of valid time stamps in the data and extrapolates it based on scanline number and scanning rate to compute new time stamps for all scanlines.
-
-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
+.. automodule:: pygac.reader
    :members:
    :undoc-members:
 
-GAC KLM reader
---------------
+GAC format reader
+~~~~~~~~~~~~~~~~~
 
-.. automodule:: pygac.gac_klm
+.. automodule:: pygac.gac_reader
    :members:
    :undoc-members:
 
 
-GAC Reader
-----------
+LAC format reader
+~~~~~~~~~~~~~~~~~
 
-.. automodule:: pygac.gac_reader
+.. automodule:: pygac.lac_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
----------------------
-
-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* currently supports calibration of all GAC data from AVHRRs onboard
- NOAA-7 and onwards, including Metop satellites.
-
-At present, updated calibration coefficients provided by Andrew Heidinger (NOAA) under SCOPE-CM project are applied for all satellites. The solar channel calibration (Channels 1 and 2, and Channel 3a if available) takes into account inter-satellite differences and is derived using amalgamation of different calibration references including the most recent MODIS Collection 6 data, in-situ targets, and simultaneous nadir observations. The detailed methodology is presented in Heidinger et al. (2010). The resulting (inter)calibration coefficients are of highest quality and follow the Global Climate Observing System (GCOS) standards for deriving fundamental climate data records.
-
-The reflectances are normalized by a factor (as a function of day of a year) to account for changing Earth-Sun distance. However, it is left to the potential user to apply further normalization using cosine of solar zenith angle (depending on application in question).
-
-
-The thermal channel intercalibration is done from scratch, starting from obtaining Platinum Resistance Thermometer (PRT), space and Internal Calibration Target (ICT, blackbody) counts. ICT temperatures are obtained from PRT counts based on coefficients provided in the POD and KLM Data User Guides (Kidwell, 2000). For each thermal channel, a smoothing window of 51 successive PRT, ICT and space counts is used to obtain robust gain values and to dampen undue high frequency fluctuations in the count data (Trishchenko, 2002). This window size is easily configurable to suit user needs. Section 7.1.2.5 of KLM Data User Guide presents the summary of equations implemented in *pygac* to calibrate thermal channels (http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c7/sec7-1.htm), including non-linearity correction (Walton et al. 1998).
-
-The PRT readings are supposed to be present in a specific order, for example, the reset value followed by four readings from the PRTs. However, this may not always be the case for orbits that contain data gaps or due to any other unexplained reason. If not taken into account properly, this irregularity could result in the underestimation of brightness temperatures, when calibration information is smoothed over many scanlines. In *pygac*, this inconsistency is handled properly while calibration thermal channels.
+POD series reader
+~~~~~~~~~~~~~~~~~
 
-In some cases it was found that, apart from the reset values, even the readings from any one of the four PRTs could also have very low suspicious values. This could also seriously affect the computation of brightness temperatures. *pygac* detects such anomalies and corrects them using interpolation of nearby valid PRT readings.
-
-
-
-
-GAC I/O module
---------------
-
-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:
-
-``ECC_GAC_avhrr_satellitename_99999_yyyymmddThhmmsstZ_yyyymmddThhmmsstZ.h5``
-
-and
-
-``ECC_GAC_sunsatangles_satellitename_99999_yyyymmddThhmmsstZ_yyyymmddThhmmsstZ.h5``
+.. automodule:: pygac.pod_reader
+   :members:
+   :undoc-members:
+   :exclude-members: tsm_affected_intervals
 
-and
 
-``ECC_GAC_qualflags_satellitename_99999_yyyymmddThhmmsstZ_yyyymmddThhmmsstZ.h5``
+KLM series reader
+~~~~~~~~~~~~~~~~~
 
-where,
+.. automodule:: pygac.klm_reader
+   :members:
+   :undoc-members:
+   :exclude-members: tsm_affected_intervals
 
-``ECC``: ESA CCI Clouds (This prefix can be changed/specified by the user)
 
-``avhrr``: denoting that it contains reflectances and BTs
+Actual Reader Implementations
+-----------------------------
 
-``sunsatangles``: denoting that it contains angles
+Actual reader implementations building upon the base classes.
 
-``qualflags``: denoting that it contains quality flag information
+GAC POD reader
+~~~~~~~~~~~~~~
 
-``yyyymmddThhmmsstZ``: yy:year, mm:month, dd:day, hh:hour, mm:min, ss:sec, t:tenth of second (for the start and the end of the orbit).
+.. automodule:: pygac.gac_pod
+   :members:
+   :undoc-members:
 
-Letters ``T`` and ``Z`` are separators for time info.
 
-The value of ``99999`` is currently used instead of providing actual orbit number.
+GAC KLM reader
+~~~~~~~~~~~~~~
 
-Appendices A, B and C provide detailed format of these files, including variable names, scaling, etc.
+.. automodule:: pygac.gac_klm
+   :members:
+   :undoc-members:
 
 
-The start and end times in the header and in actual L1b data can be different for orbits. The mismatch can range from few milliseconds to days. It was decided to trust the time stamps in L1b data in *pygac*. After reorganizing based on scanline numbers (see issue highlighted above), the time stamps from the first and the last scanlines are taken as start and end times.
+LAC POD reader
+~~~~~~~~~~~~~~
 
+.. automodule:: pygac.lac_pod
+   :members:
+   :undoc-members:
 
-In some orbits, the latitude and longitude information contains corrupt values for a part/s of the orbit and these scanlines are not flagged in the corresponding scanline-by-scaline quality flags. Currently, *pygac* uses only a simple if_else construct to constrain valid range. Further improvement could be done using extra- and interpolation techniques.
 
-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.
+LAC KLM reader
+~~~~~~~~~~~~~~
 
-.. automodule:: pygac.gac_io
+.. automodule:: pygac.lac_klm
    :members:
    :undoc-members:
 
+Calibration
+-----------
 
-Supplement A: Structure of an output file containing reflectances and brightness temperatures
----------------------------------------------------------------------------------------------
-   
-Input: L1b file: ``NSS.GHRR.NN.D06279.S1800.E1955.B0711012.GC``
-
-Output::
-
-  Group how {
-     variables:
-       char channel_list(6, 9);
-         :_lastModified = "2014-10-29T13:55:05Z";
-
-     how:yaw_error = 0.0; // double
-     how:roll_error = 0.0; // double
-     how:pich_error = 0.0; // double
-     how:startepochs = 1160157611L; // long
-     how:endepochs = 1160164510L; // long
-     how:platform = "noaa18";
-     how:instrument = "avhrr";
-     how:orbit_number = 99999; // int
-     how:software = "pyGAC";
-     how:version = "1.0";
-   }
-
-   Group image1 {
-     variables:
-       short data(13686, 409);
-         :_lastModified = "2014-10-29T13:55:05Z";
-
-     Group how {
-       how:sun_earth_distance_correction_applied = "TRUE";
-       how:sun_earth_distance_correction_factor = 0.9982412208987179; // double
-     }
-
-     Group what {
-       what:product = "SATCH";
-       what:quantity = "REFL";
-       what:dataset_name = "Channel 1 reflectance";
-       what:units = "%";
-       what:gain = 0.01f; // float
-       what:offset = 0.0f; // float
-       what:missingdata = -32001; // int
-       what:nodata = -32001; // int
-       what:starttime = "180011";
-       what:endtime = "195510";
-       what:startdate = "20061006";
-       what:enddate = "20061006";
-     }
-
-     image1:channel = "1";
-     image1:description = "AVHRR ch1";
-   }
-
-   Group image2 {
-     variables:
-       short data(13686, 409);
-         :_lastModified = "2014-10-29T13:55:05Z";
-
-     Group how {
-       how:sun_earth_distance_correction_applied = "TRUE";
-       how:sun_earth_distance_correction_factor = 0.9982412208987179; // double
-     }
-
-     Group what {
-       what:product = "SATCH";
-       what:quantity = "REFL";
-       what:dataset_name = "Channel 2 reflectance";
-       what:units = "%";
-       what:gain = 0.01f; // float
-       what:offset = 0.0f; // float
-       what:missingdata = -32001; // int
-       what:nodata = -32001; // int
-       what:starttime = "180011";
-       what:endtime = "195510";
-       what:startdate = "20061006";
-       what:enddate = "20061006";
-     }
-
-     image2:channel = "2";
-     image2:description = "AVHRR ch2";
-   }
-
-   Group image3 {
-     variables:
-       short data(13686, 409);
-         :_lastModified = "2014-10-29T13:55:05Z";
-
-     Group how {
-     }
-
-     Group what {
-       what:product = "SATCH";
-       what:quantity = "TB";
-       what:dataset_name = "Channel 3b brightness temperature";
-       what:units = "K";
-       what:gain = 0.01f; // float
-       what:offset = 273.15f; // float
-       what:missingdata = -32001; // int
-       what:nodata = -32001; // int
-       what:starttime = "180011";
-       what:endtime = "195510";
-       what:startdate = "20061006";
-       what:enddate = "20061006";
-     }
-
-     image3:description = "AVHRR ch3b";
-     image3:channel = "3b";
-   }
-
-   Group image4 {
-     variables:
-       short data(13686, 409);
-         :_lastModified = "2014-10-29T13:55:05Z";
-
-     Group how {
-     }
-
-     Group what {
-       what:product = "SATCH";
-       what:quantity = "TB";
-       what:dataset_name = "Channel 4 brightness temperature";
-       what:units = "K";
-       what:gain = 0.01f; // float
-       what:offset = 273.15f; // float
-       what:missingdata = -32001; // int
-       what:nodata = -32001; // int
-       what:starttime = "180011";
-       what:endtime = "195510";
-       what:startdate = "20061006";
-       what:enddate = "20061006";
-     }
-
-     image4:channel = "4";
-     image4:description = "AVHRR ch4";
-   }
-
-   Group image5 {
-     variables:
-       short data(13686, 409);
-         :_lastModified = "2014-10-29T13:55:05Z";
-
-     Group how {
-     }
-
-     Group what {
-       what:product = "SATCH";
-       what:quantity = "TB";
-       what:dataset_name = "Channel 5 brightness temperature";
-       what:units = "K";
-       what:gain = 0.01f; // float
-       what:offset = 273.15f; // float
-       what:missingdata = -32001; // int
-       what:nodata = -32001; // int
-       what:starttime = "180011";
-       what:endtime = "195510";
-       what:startdate = "20061006";
-       what:enddate = "20061006";
-     }
-
-     image5:channel = "5";
-     image5:description = "AVHRR ch5";
-   }
-
-   Group image6 {
-     variables:
-       short data(13686, 409);
-         :_lastModified = "2014-10-29T13:55:05Z";
-
-     Group how {
-       how:sun_earth_distance_correction_applied = "TRUE";
-       how:sun_earth_distance_correction_factor = 0.9982412208987179; // double
-     }
-
-     Group what {
-       what:product = "SATCH";
-       what:quantity = "REFL";
-       what:dataset_name = "Channel 3a reflectance";
-       what:units = "%";
-       what:gain = 0.01f; // float
-       what:offset = 0.0f; // float
-       what:missingdata = -32001; // int
-       what:nodata = -32001; // int
-       what:starttime = "180011";
-       what:endtime = "195510";
-       what:startdate = "20061006";
-       what:enddate = "20061006";
-     }
-
-     image6:channel = "3a";
-     image6:description = "AVHRR ch3a";
-   }
-
-   Group what {
-     what:object = "SATP";
-     what:sets = 6; // int
-     what:version = "H5rad ?.?";
-     what:date = "20061006";
-     what:time = "180011";
-   }
-
-   Group where {
-
-     Group lat {
-       variables:
-         int data(13686, 409);
-           :_lastModified = "2014-10-29T13:55:05Z";
-
-       Group what {
-         what:dataset_name = "Latitude";
-         what:units = "Deg";
-         what:gain = 0.001f; // float
-         what:offset = 0.0f; // float
-         what:missingdata = -32001; // int
-         what:nodata = -32001; // int
-         what:starttime = "180011";
-         what:endtime = "195510";
-         what:startdate = "20061006";
-         what:enddate = "20061006";
-       }
-     }
-
-     Group lon {
-       variables:
-         int data(13686, 409);
-           :_lastModified = "2014-10-29T13:55:05Z";
-
-       Group what {
-         what:dataset_name = "Longitude";
-         what:units = "Deg";
-         what:gain = 0.001f; // float
-         what:offset = 0.0f; // float
-         what:missingdata = -32001; // int
-         what:nodata = -32001; // int
-         what:starttime = "180011";
-         what:endtime = "195510";
-         what:startdate = "20061006";
-         what:enddate = "20061006";
-       }
-     }
-
-     where:num_of_pixels = 409; // int
-     where:num_of_lines = 13686; // int
-     where:xscale = 0.0f; // float
-     where:yscale = 0.0f; // float
-   }
-  }
-
-
-
-Supplement B: Structure of an output file containing Sun and satellite positions
---------------------------------------------------------------------------------
-
-Input: L1b file: ``NSS.GHRR.NN.D06279.S1800.E1955.B0711012.GC``
-
-Output::
-
-  Group how {
-     how:yaw_error = 0.0; // double
-     how:roll_error = 0.0; // double
-     how:pich_error = 0.0; // double
-     how:startepochs = 1160157611L; // long
-     how:endepochs = 1160164510L; // long
-     how:platform = "noaa18";
-     how:instrument = "avhrr";
-     how:orbit_number = 99999; // int
-     how:software = "pyGAC";
-     how:version = "1.0";
-   }
-
-   Group image1 {
-     variables:
-       short data(13686, 409);
-         :_lastModified = "2014-10-29T13:55:06Z";
-
-     Group what {
-       what:product = "SUNZ";
-       what:quantity = "DEG";
-       what:dataset_name = "Solar zenith angle";
-       what:units = "Deg";
-       what:gain = 0.01f; // float
-       what:offset = 0.0f; // float
-       what:missingdata = -32001; // int
-       what:nodata = -32001; // int
-       what:starttime = "180011";
-       what:endtime = "195510";
-       what:startdate = "20061006";
-       what:enddate = "20061006";
-     }
-
-     image1:description = "Solar zenith angle";
-   }
-
-   Group image2 {
-     variables:
-       short data(13686, 409);
-         :_lastModified = "2014-10-29T13:55:06Z";
-
-     Group what {
-       what:product = "SATZ";
-       what:quantity = "DEG";
-       what:dataset_name = "Satellite zenith angle";
-       what:units = "Deg";
-       what:gain = 0.01f; // float
-       what:offset = 0.0f; // float
-       what:missingdata = -32001; // int
-       what:nodata = -32001; // int
-       what:starttime = "180011";
-       what:endtime = "195510";
-       what:startdate = "20061006";
-       what:enddate = "20061006";
-     }
-
-     image2:description = "Satellite zenith angle";
-   }
-
-   Group image3 {
-     variables:
-       short data(13686, 409);
-         :_lastModified = "2014-10-29T13:55:06Z";
-
-     Group what {
-       what:product = "SSAZD";
-       what:quantity = "DEG";
-       what:dataset_name = "Relative satellite-sun azimuth angle";
-       what:units = "Deg";
-       what:gain = 0.01f; // float
-       what:offset = 0.0f; // float
-       what:missingdata = -32001; // int
-       what:nodata = -32001; // int
-       what:starttime = "180011";
-       what:endtime = "195510";
-       what:startdate = "20061006";
-       what:enddate = "20061006";
-     }
-
-     image3:description = "Relative satellite-sun azimuth angle";
-   }
-
-   Group image4 {
-     variables:
-       short data(13686, 409);
-         :_lastModified = "2014-10-29T13:55:06Z";
-
-     Group what {
-       what:product = "SUNA";
-       what:quantity = "DEG";
-       what:dataset_name = "Solar azimuth angle";
-       what:units = "Deg";
-       what:gain = 0.01f; // float
-       what:offset = 180.0f; // float
-       what:missingdata = -32001; // int
-       what:nodata = -32001; // int
-       what:starttime = "180011";
-       what:endtime = "195510";
-       what:startdate = "20061006";
-       what:enddate = "20061006";
-     }
-
-     image4:description = "Solar azimuth angle";
-   }
-
-   Group image5 {
-     variables:
-       short data(13686, 409);
-         :_lastModified = "2014-10-29T13:55:06Z";
-
-     Group what {
-       what:product = "SATA";
-       what:quantity = "DEG";
-       what:dataset_name = "Satellite azimuth angle";
-       what:units = "Deg";
-       what:gain = 0.01f; // float
-       what:offset = 180.0f; // float
-       what:missingdata = -32001; // int
-       what:nodata = -32001; // int
-       what:starttime = "180011";
-       what:endtime = "195510";
-       what:startdate = "20061006";
-       what:enddate = "20061006";
-     }
-
-     image5:description = "Satellite azimuth angle";
-   }
-
-   Group what {
-     what:object = "SATP";
-     what:sets = 5; // int
-     what:version = "H5rad ?.?";
-     what:date = "20061006";
-     what:time = "180011";
-   }
-
-   Group where {
-
-     Group lat {
-       variables:
-         int data(13686, 409);
-           :_lastModified = "2014-10-29T13:55:06Z";
-
-       Group what {
-         what:dataset_name = "Latitude";
-         what:units = "Deg";
-         what:gain = 0.001f; // float
-         what:offset = 0.0f; // float
-         what:missingdata = -32001; // int
-         what:nodata = -32001; // int
-         what:starttime = "180011";
-         what:endtime = "195510";
-         what:startdate = "20061006";
-         what:enddate = "20061006";
-       }
-     }
-
-     Group lon {
-       variables:
-         int data(13686, 409);
-           :_lastModified = "2014-10-29T13:55:06Z";
-
-       Group what {
-         what:dataset_name = "Longitude";
-         what:units = "Deg";
-         what:gain = 0.001f; // float
-         what:offset = 0.0f; // float
-         what:missingdata = -32001; // int
-         what:nodata = -32001; // int
-         what:starttime = "180011";
-         what:endtime = "195510";
-         what:startdate = "20061006";
-         what:enddate = "20061006";
-       }
-     }
-
-     where:num_of_pixels = 409; // int
-     where:num_of_lines = 13686; // int
-     where:xscale = 0.0f; // float
-     where:yscale = 0.0f; // float
-   }
-  }
-
-
-Supplement C: Structure of an output file containing quality flags
-------------------------------------------------------------------
-
-
-The file that contains quality flags has following information.
-
-1) There is a variable called ``qual_flags/data``. It will have a dimension of (X,7), where X is the number of data records in the GAC orbit.
-
-The 7 columns contain the following information.
-
-=======  ==========================================================================================
-Col 1    Scan line number
-Col 2    Fatal error flag (scan line should not be used for analysis).
-Col 3    Insufficient data for calibration (scan line should not be used for analysis).
-Col 4    Insufficient data for navigation (scan line should not be used for analysis).
-Col 5-7  whether solar contamination of blackbody occurred in in channels 3, 4, and 5 respectively.
-=======  ==========================================================================================
-
-If the values for these flags are greater than zero, then the data should not be used. If everything is normal, then all values should be zero.
-
-2) There are also two important attributes that provide "last scan line number" and "total number of data records".
-
-By combining information from column 1 and these two attributes,  the user is able to figure out where the gap occurs and also exact time for each scan line.
+.. automodule:: pygac.calibration
+   :members:
+   :undoc-members:
diff --git a/doc/source/conf.py b/doc/source/conf.py
index 63e9bc5..388725e 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -25,7 +25,8 @@
 
 # Add any Sphinx extension module names here, as strings. They can be extensions
 # coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
-extensions = ['sphinx.ext.autodoc', 'sphinx.ext.doctest', 'sphinx.ext.todo']
+extensions = ['sphinx.ext.autodoc', 'sphinx.ext.doctest', 'sphinx.ext.todo',
+              'sphinx.ext.inheritance_diagram']
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
diff --git a/doc/source/index.rst b/doc/source/index.rst
index 020fabd..c682565 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -3,19 +3,27 @@
    You can adapt this file completely to your liking, but it should at least
    contain the root `toctree` directive.
 
-Welcome to pygac's documentation!
-=================================
+Pygac's documentation
+=====================
 
-Contents:
+Pygac is a Python package to read, calibrate and navigate data from the AVHRR
+instrument onboard NOAA and MetOp satellites in GAC and LAC format.
+
+
+Table of Contents
+-----------------
 
 .. toctree::
    :maxdepth: 2
 
+   introduction
    installation
    usage
+   methods
    api
-   
-   
+   legacy
+
+
 
 Indices and tables
 ==================
diff --git a/doc/source/installation.rst b/doc/source/installation.rst
index b9d7001..e5558fc 100644
--- a/doc/source/installation.rst
+++ b/doc/source/installation.rst
@@ -1,57 +1,57 @@
 Installation
 ------------
 
-pygac reads NOAA AVHRR Global Area Coverage (GAC) data, and performs state of
-the art calibration and navigation.
+You can install the latest stable release of the software via the python package index (pypi)
 
-You can install the latest stable release of the software via the python package index (pypi)::
+.. code-block:: bash
 
-  $> pip install --prefix=$PREFIX pygac
+  pip install pygac
 
-where ``$PREFIX`` denotes your desired installation prefix. If you want access to the full 
-source code and/or test the latest version under development you can download the pygac 
-source code from github_::
 
-  $> git clone git://github.com/pytroll/pygac.git
-  $> cd pygac
+TLE files
+~~~~~~~~~
+The pygac package requires Two-Line Element files stored per-satellite
+in files with names such as TLE_noaa19.txt. The contents should be the
+historical TLEs, i.e. a concatenation of just lines 1 and 2 without the
+satellite name. For example
 
-and then run::
+.. code-block::
 
-  $> pip install --prefix=$PREFIX .
+    1 23455U 94089A   01122.93455091  .00000622  00000-0  36103-3 0  7210
+    2 23455  99.1771 113.3063 0008405 277.6106  82.4106 14.12671703326608
+    1 23455U 94089A   01326.97611660  .00000739  00000-0  42245-3 0  9806
+    2 23455  99.1886 322.4670 0009980  66.2863 293.9354 14.12871991355419
+    etc
 
-or, if you want to hack the package::
+These can be downloaded from CelesTrak via the `special data request form`_.
 
-  $> pip install -e --prefix=$PREFIX .
+.. _special data request form:
+    https://celestrak.com/NORAD/archives/request.php
 
-To uninstall the package run::
 
-  $> pip uninstall pygac
+Development
+~~~~~~~~~~~
 
+For development clone the repository from github and install pygac in editable mode.
 
-At runtime
-----------
+.. code-block:: bash
 
-Let $PREFIX be the prefix of your pygac installation. To make the python 
-interpreter aware of the pygac module, update the ``$PYTHONPATH`` environment
-variable::
+  git clone git://github.com/pytroll/pygac
+  cd pygac
+  pip install -e .[dev]
 
-  $> export PYTHONPATH=$PREFIX/lib/python2.7/site-packages:$PYTHONPATH
+It is recommended to activate `pre-commit`_ checks.
 
-Furthermore, update the ``$PATH`` environment variable to have the pygac 
-scripts available in your session::
+.. code-block:: bash
 
-  $> export PATH=$PREFIX/bin:$PATH
+  pre-commit install
 
+The test suite can be run using pytest.
 
-.. _github: http://github.com/pytroll/pygac
+.. code-block:: bash
+
+  pytest -vs pygac/tests
 
-TLE files
----------
-The pygac package requires Two-Line Element files stored per-satellite
-in files with names such as TLE_noaa19.txt. The name format and directory can be
-configured in the config file (see the Usage section). The contents should be the
-historical TLEs, i.e. a concatenation of just lines 1 and 2 without the satellite
-name. These can be downloaded from space-track.org. Once logged in use URL such as
-(replace CCCCC with the five-digit NORAD catalogue ID):
 
-https://www.space-track.org/basicspacedata/query/class/tle/EPOCH/1990-01-01--2029-01-01/NORAD_CAT_ID/CCCCC/orderby/EPOCH ASC/format/tle
+.. _pre-commit:
+    https://pre-commit.com/
diff --git a/doc/source/introduction.rst b/doc/source/introduction.rst
new file mode 100644
index 0000000..68c3129
--- /dev/null
+++ b/doc/source/introduction.rst
@@ -0,0 +1,47 @@
+Introduction
+============
+
+Supported Data Format
+---------------------
+
+Pygac reads AVHRR GAC (Global Area Coverage) and LAC (Local Area Coverage)
+level 1b data from NOAA, which is described in the `POD`_ (NOAA-14 and
+before) and `KLM`_ (NOAA-15 and following) user guides. The data can be
+obtained from `NOAA CLASS`_, where you can also find a comprehensive
+`introduction`_.
+
+.. note::
+    Pygac is currently not able to read files with the CLASS archive
+    header included.
+
+
+.. _NOAA CLASS:
+    https://www.class.noaa.gov/
+.. _POD:
+    https://www.ncei.noaa.gov/pub/data/satellite/publications/podguides/N-15%20thru%20N-19/
+.. _KLM:
+    https://www.ncei.noaa.gov/pub/data/satellite/publications/podguides/TIROS-N%20thru%20N-14/
+.. _introduction:
+    https://www.class.noaa.gov/release/data_available/avhrr/index.htm
+
+
+Supported Sensors
+-----------------
+Pygac currently supports AVHRR generations 1-3 onboard NOAA (TIROS-N, NOAA-6
+and onwards) and MetOp satellites.
+
+
+.. _here:
+    https://www.avl.class.noaa.gov/release/data_available/avhrr/index.htm
+
+
+Related Projects
+----------------
+
+- `pygac-fdr`_: Generate a fundamental data record of AVHRR GAC data using
+  Pygac.
+- `level1c4pps`_: Prepare AVHRR GAC data for `NWCSAF/PPS`_ using Pygac.
+
+.. _level1c4pps: https://github.com/foua-pps/level1c4pps
+.. _NWCSAF/PPS: https://www.nwcsaf.org/16
+.. _pygac-fdr: https://github.com/pytroll/pygac-fdr
\ No newline at end of file
diff --git a/doc/source/legacy.rst b/doc/source/legacy.rst
new file mode 100644
index 0000000..ce3cadf
--- /dev/null
+++ b/doc/source/legacy.rst
@@ -0,0 +1,535 @@
+Legacy Output
+=============
+
+GAC I/O module
+--------------
+
+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:
+
+``ECC_GAC_avhrr_satellitename_99999_yyyymmddThhmmsstZ_yyyymmddThhmmsstZ.h5``
+
+and
+
+``ECC_GAC_sunsatangles_satellitename_99999_yyyymmddThhmmsstZ_yyyymmddThhmmsstZ.h5``
+
+and
+
+``ECC_GAC_qualflags_satellitename_99999_yyyymmddThhmmsstZ_yyyymmddThhmmsstZ.h5``
+
+where,
+
+``ECC``: ESA CCI Clouds (This prefix can be changed/specified by the user)
+
+``avhrr``: denoting that it contains reflectances and BTs
+
+``sunsatangles``: denoting that it contains angles
+
+``qualflags``: denoting that it contains quality flag information
+
+``yyyymmddThhmmsstZ``: yy:year, mm:month, dd:day, hh:hour, mm:min, ss:sec, t:tenth of second (for the start and the end of the orbit).
+
+Letters ``T`` and ``Z`` are separators for time info.
+
+The value of ``99999`` is currently used instead of providing actual orbit
+number.
+
+Appendices A, B and C provide detailed format of these files, including
+variable names, scaling, etc.
+
+
+The start and end times in the header and in actual L1b data can be different
+for orbits. The mismatch can range from few milliseconds to days. It was
+decided to trust the time stamps in L1b data in Pygac. After reorganizing
+based on scanline numbers (see issue highlighted above), the time stamps
+from the first and the last scanlines are taken as start and end times.
+
+
+In some orbits, the latitude and longitude information contains corrupt
+values for a part/s of the orbit and these scanlines are not flagged in the
+corresponding scanline-by-scaline quality flags. Currently, Pygac uses only
+a simple if_else construct to constrain valid range. Further improvement could
+be done using extra- and interpolation techniques.
+
+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.
+
+
+.. 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``
+
+Output::
+
+  Group how {
+     variables:
+       char channel_list(6, 9);
+         :_lastModified = "2014-10-29T13:55:05Z";
+
+     how:yaw_error = 0.0; // double
+     how:roll_error = 0.0; // double
+     how:pich_error = 0.0; // double
+     how:startepochs = 1160157611L; // long
+     how:endepochs = 1160164510L; // long
+     how:platform = "noaa18";
+     how:instrument = "avhrr";
+     how:orbit_number = 99999; // int
+     how:software = "pyGAC";
+     how:version = "1.0";
+   }
+
+   Group image1 {
+     variables:
+       short data(13686, 409);
+         :_lastModified = "2014-10-29T13:55:05Z";
+
+     Group how {
+       how:sun_earth_distance_correction_applied = "TRUE";
+       how:sun_earth_distance_correction_factor = 0.9982412208987179; // double
+     }
+
+     Group what {
+       what:product = "SATCH";
+       what:quantity = "REFL";
+       what:dataset_name = "Channel 1 reflectance";
+       what:units = "%";
+       what:gain = 0.01f; // float
+       what:offset = 0.0f; // float
+       what:missingdata = -32001; // int
+       what:nodata = -32001; // int
+       what:starttime = "180011";
+       what:endtime = "195510";
+       what:startdate = "20061006";
+       what:enddate = "20061006";
+     }
+
+     image1:channel = "1";
+     image1:description = "AVHRR ch1";
+   }
+
+   Group image2 {
+     variables:
+       short data(13686, 409);
+         :_lastModified = "2014-10-29T13:55:05Z";
+
+     Group how {
+       how:sun_earth_distance_correction_applied = "TRUE";
+       how:sun_earth_distance_correction_factor = 0.9982412208987179; // double
+     }
+
+     Group what {
+       what:product = "SATCH";
+       what:quantity = "REFL";
+       what:dataset_name = "Channel 2 reflectance";
+       what:units = "%";
+       what:gain = 0.01f; // float
+       what:offset = 0.0f; // float
+       what:missingdata = -32001; // int
+       what:nodata = -32001; // int
+       what:starttime = "180011";
+       what:endtime = "195510";
+       what:startdate = "20061006";
+       what:enddate = "20061006";
+     }
+
+     image2:channel = "2";
+     image2:description = "AVHRR ch2";
+   }
+
+   Group image3 {
+     variables:
+       short data(13686, 409);
+         :_lastModified = "2014-10-29T13:55:05Z";
+
+     Group how {
+     }
+
+     Group what {
+       what:product = "SATCH";
+       what:quantity = "TB";
+       what:dataset_name = "Channel 3b brightness temperature";
+       what:units = "K";
+       what:gain = 0.01f; // float
+       what:offset = 273.15f; // float
+       what:missingdata = -32001; // int
+       what:nodata = -32001; // int
+       what:starttime = "180011";
+       what:endtime = "195510";
+       what:startdate = "20061006";
+       what:enddate = "20061006";
+     }
+
+     image3:description = "AVHRR ch3b";
+     image3:channel = "3b";
+   }
+
+   Group image4 {
+     variables:
+       short data(13686, 409);
+         :_lastModified = "2014-10-29T13:55:05Z";
+
+     Group how {
+     }
+
+     Group what {
+       what:product = "SATCH";
+       what:quantity = "TB";
+       what:dataset_name = "Channel 4 brightness temperature";
+       what:units = "K";
+       what:gain = 0.01f; // float
+       what:offset = 273.15f; // float
+       what:missingdata = -32001; // int
+       what:nodata = -32001; // int
+       what:starttime = "180011";
+       what:endtime = "195510";
+       what:startdate = "20061006";
+       what:enddate = "20061006";
+     }
+
+     image4:channel = "4";
+     image4:description = "AVHRR ch4";
+   }
+
+   Group image5 {
+     variables:
+       short data(13686, 409);
+         :_lastModified = "2014-10-29T13:55:05Z";
+
+     Group how {
+     }
+
+     Group what {
+       what:product = "SATCH";
+       what:quantity = "TB";
+       what:dataset_name = "Channel 5 brightness temperature";
+       what:units = "K";
+       what:gain = 0.01f; // float
+       what:offset = 273.15f; // float
+       what:missingdata = -32001; // int
+       what:nodata = -32001; // int
+       what:starttime = "180011";
+       what:endtime = "195510";
+       what:startdate = "20061006";
+       what:enddate = "20061006";
+     }
+
+     image5:channel = "5";
+     image5:description = "AVHRR ch5";
+   }
+
+   Group image6 {
+     variables:
+       short data(13686, 409);
+         :_lastModified = "2014-10-29T13:55:05Z";
+
+     Group how {
+       how:sun_earth_distance_correction_applied = "TRUE";
+       how:sun_earth_distance_correction_factor = 0.9982412208987179; // double
+     }
+
+     Group what {
+       what:product = "SATCH";
+       what:quantity = "REFL";
+       what:dataset_name = "Channel 3a reflectance";
+       what:units = "%";
+       what:gain = 0.01f; // float
+       what:offset = 0.0f; // float
+       what:missingdata = -32001; // int
+       what:nodata = -32001; // int
+       what:starttime = "180011";
+       what:endtime = "195510";
+       what:startdate = "20061006";
+       what:enddate = "20061006";
+     }
+
+     image6:channel = "3a";
+     image6:description = "AVHRR ch3a";
+   }
+
+   Group what {
+     what:object = "SATP";
+     what:sets = 6; // int
+     what:version = "H5rad ?.?";
+     what:date = "20061006";
+     what:time = "180011";
+   }
+
+   Group where {
+
+     Group lat {
+       variables:
+         int data(13686, 409);
+           :_lastModified = "2014-10-29T13:55:05Z";
+
+       Group what {
+         what:dataset_name = "Latitude";
+         what:units = "Deg";
+         what:gain = 0.001f; // float
+         what:offset = 0.0f; // float
+         what:missingdata = -32001; // int
+         what:nodata = -32001; // int
+         what:starttime = "180011";
+         what:endtime = "195510";
+         what:startdate = "20061006";
+         what:enddate = "20061006";
+       }
+     }
+
+     Group lon {
+       variables:
+         int data(13686, 409);
+           :_lastModified = "2014-10-29T13:55:05Z";
+
+       Group what {
+         what:dataset_name = "Longitude";
+         what:units = "Deg";
+         what:gain = 0.001f; // float
+         what:offset = 0.0f; // float
+         what:missingdata = -32001; // int
+         what:nodata = -32001; // int
+         what:starttime = "180011";
+         what:endtime = "195510";
+         what:startdate = "20061006";
+         what:enddate = "20061006";
+       }
+     }
+
+     where:num_of_pixels = 409; // int
+     where:num_of_lines = 13686; // int
+     where:xscale = 0.0f; // float
+     where:yscale = 0.0f; // float
+   }
+  }
+
+
+
+Supplement B: Structure of an output file containing Sun and satellite positions
+--------------------------------------------------------------------------------
+
+Input: L1b file: ``NSS.GHRR.NN.D06279.S1800.E1955.B0711012.GC``
+
+Output::
+
+  Group how {
+     how:yaw_error = 0.0; // double
+     how:roll_error = 0.0; // double
+     how:pich_error = 0.0; // double
+     how:startepochs = 1160157611L; // long
+     how:endepochs = 1160164510L; // long
+     how:platform = "noaa18";
+     how:instrument = "avhrr";
+     how:orbit_number = 99999; // int
+     how:software = "pyGAC";
+     how:version = "1.0";
+   }
+
+   Group image1 {
+     variables:
+       short data(13686, 409);
+         :_lastModified = "2014-10-29T13:55:06Z";
+
+     Group what {
+       what:product = "SUNZ";
+       what:quantity = "DEG";
+       what:dataset_name = "Solar zenith angle";
+       what:units = "Deg";
+       what:gain = 0.01f; // float
+       what:offset = 0.0f; // float
+       what:missingdata = -32001; // int
+       what:nodata = -32001; // int
+       what:starttime = "180011";
+       what:endtime = "195510";
+       what:startdate = "20061006";
+       what:enddate = "20061006";
+     }
+
+     image1:description = "Solar zenith angle";
+   }
+
+   Group image2 {
+     variables:
+       short data(13686, 409);
+         :_lastModified = "2014-10-29T13:55:06Z";
+
+     Group what {
+       what:product = "SATZ";
+       what:quantity = "DEG";
+       what:dataset_name = "Satellite zenith angle";
+       what:units = "Deg";
+       what:gain = 0.01f; // float
+       what:offset = 0.0f; // float
+       what:missingdata = -32001; // int
+       what:nodata = -32001; // int
+       what:starttime = "180011";
+       what:endtime = "195510";
+       what:startdate = "20061006";
+       what:enddate = "20061006";
+     }
+
+     image2:description = "Satellite zenith angle";
+   }
+
+   Group image3 {
+     variables:
+       short data(13686, 409);
+         :_lastModified = "2014-10-29T13:55:06Z";
+
+     Group what {
+       what:product = "SSAZD";
+       what:quantity = "DEG";
+       what:dataset_name = "Relative satellite-sun azimuth angle";
+       what:units = "Deg";
+       what:gain = 0.01f; // float
+       what:offset = 0.0f; // float
+       what:missingdata = -32001; // int
+       what:nodata = -32001; // int
+       what:starttime = "180011";
+       what:endtime = "195510";
+       what:startdate = "20061006";
+       what:enddate = "20061006";
+     }
+
+     image3:description = "Relative satellite-sun azimuth angle";
+   }
+
+   Group image4 {
+     variables:
+       short data(13686, 409);
+         :_lastModified = "2014-10-29T13:55:06Z";
+
+     Group what {
+       what:product = "SUNA";
+       what:quantity = "DEG";
+       what:dataset_name = "Solar azimuth angle";
+       what:units = "Deg";
+       what:gain = 0.01f; // float
+       what:offset = 180.0f; // float
+       what:missingdata = -32001; // int
+       what:nodata = -32001; // int
+       what:starttime = "180011";
+       what:endtime = "195510";
+       what:startdate = "20061006";
+       what:enddate = "20061006";
+     }
+
+     image4:description = "Solar azimuth angle";
+   }
+
+   Group image5 {
+     variables:
+       short data(13686, 409);
+         :_lastModified = "2014-10-29T13:55:06Z";
+
+     Group what {
+       what:product = "SATA";
+       what:quantity = "DEG";
+       what:dataset_name = "Satellite azimuth angle";
+       what:units = "Deg";
+       what:gain = 0.01f; // float
+       what:offset = 180.0f; // float
+       what:missingdata = -32001; // int
+       what:nodata = -32001; // int
+       what:starttime = "180011";
+       what:endtime = "195510";
+       what:startdate = "20061006";
+       what:enddate = "20061006";
+     }
+
+     image5:description = "Satellite azimuth angle";
+   }
+
+   Group what {
+     what:object = "SATP";
+     what:sets = 5; // int
+     what:version = "H5rad ?.?";
+     what:date = "20061006";
+     what:time = "180011";
+   }
+
+   Group where {
+
+     Group lat {
+       variables:
+         int data(13686, 409);
+           :_lastModified = "2014-10-29T13:55:06Z";
+
+       Group what {
+         what:dataset_name = "Latitude";
+         what:units = "Deg";
+         what:gain = 0.001f; // float
+         what:offset = 0.0f; // float
+         what:missingdata = -32001; // int
+         what:nodata = -32001; // int
+         what:starttime = "180011";
+         what:endtime = "195510";
+         what:startdate = "20061006";
+         what:enddate = "20061006";
+       }
+     }
+
+     Group lon {
+       variables:
+         int data(13686, 409);
+           :_lastModified = "2014-10-29T13:55:06Z";
+
+       Group what {
+         what:dataset_name = "Longitude";
+         what:units = "Deg";
+         what:gain = 0.001f; // float
+         what:offset = 0.0f; // float
+         what:missingdata = -32001; // int
+         what:nodata = -32001; // int
+         what:starttime = "180011";
+         what:endtime = "195510";
+         what:startdate = "20061006";
+         what:enddate = "20061006";
+       }
+     }
+
+     where:num_of_pixels = 409; // int
+     where:num_of_lines = 13686; // int
+     where:xscale = 0.0f; // float
+     where:yscale = 0.0f; // float
+   }
+  }
+
+
+Supplement C: Structure of an output file containing quality flags
+------------------------------------------------------------------
+
+
+The file that contains quality flags has following information.
+
+1) There is a variable called ``qual_flags/data``. It will have a dimension of
+(X,7), where X is the number of data records in the GAC orbit. The 7 columns
+contain the following information.
+
+=======  ==========================================================================================
+Col 1    Scan line number
+Col 2    Fatal error flag (scan line should not be used for analysis).
+Col 3    Insufficient data for calibration (scan line should not be used for analysis).
+Col 4    Insufficient data for navigation (scan line should not be used for analysis).
+Col 5-7  whether solar contamination of blackbody occurred in in channels 3, 4, and 5 respectively.
+=======  ==========================================================================================
+
+If the values for these flags are greater than zero, then the data should not
+be used. If everything is normal, then all values should be zero.
+
+2) There are also two important attributes that provide "last scan line number"
+and "total number of data records".
+
+By combining information from column 1 and these two attributes, the user is
+able to figure out where the gap occurs and also exact time for each scan line.
diff --git a/doc/source/methods.rst b/doc/source/methods.rst
new file mode 100644
index 0000000..e3ebb13
--- /dev/null
+++ b/doc/source/methods.rst
@@ -0,0 +1,151 @@
+Methods
+=======
+
+Calibration
+-----------
+
+At present, calibration coefficients provided by Andrew Heidinger
+(NOAA) under SCOPE-CM project are applied for all satellites.
+
+Pygac comes with one default calibration file, which can be found in
+`pygac/data/calibration.json`_. The current version is *PATMOS-x, v2017r1*,
+including provisional coefficients for MetOp-C. A record of all versions is
+kept in :class:`pygac.calibration.Calibrator.version_hashs`. Alternatively, it
+is possible to pass custom coefficients to the reader, see
+:class:`pygac.reader.Reader`. This could also be a previous default version.
+
+.. _pygac/data/calibration.json:
+    https://github.com/pytroll/pygac/blob/main/pygac/data/calibration.json
+
+The solar channel calibration (Channels 1 and 2, and Channel 3a if available)
+takes into account inter-satellite differences and is derived using
+amalgamation of different calibration references including the most recent
+MODIS Collection 6 data, in-situ targets, and simultaneous nadir
+observations. The detailed methodology is presented in Heidinger et al.
+(2010). The resulting (inter)calibration coefficients are of highest quality
+and follow the Global Climate Observing System (GCOS) standards for
+deriving fundamental climate data records.
+
+The reflectances are normalized by a factor (as a function of day of a year)
+to account for changing Earth-Sun distance. However, it is left to the
+user to apply further normalization using cosine of solar zenith
+angle for example (depending on application in question).
+
+The thermal channel intercalibration is done from scratch, starting from
+obtaining Platinum Resistance Thermometer (PRT), space and Internal
+Calibration Target (ICT, blackbody) counts. ICT temperatures are obtained
+from PRT counts based on coefficients provided in the POD and KLM Data User
+Guides (Kidwell, 2000). For each thermal channel, a smoothing window of 51
+successive PRT, ICT and space counts is used to obtain robust gain values and
+to dampen undue high frequency fluctuations in the count data (Trishchenko,
+2002). Section 7.1.2.5 of `KLM User Guide`_ presents the summary of equations
+implemented in Pygac to calibrate thermal channels, including non-linearity
+correction (Walton et al. 1998).
+
+.. _KLM User Guide:
+    https://www.ncei.noaa.gov/pub/data/satellite/publications/podguides/TIROS-N%20thru%20N-14/
+
+The PRT readings are supposed to be present in a specific order, for example,
+the reset value followed by four readings from the PRTs. However, this may
+not always be the case for orbits that contain data gaps or due to any other
+unexplained reason. If not taken into account properly, this irregularity
+could result in the underestimation of brightness temperatures, when
+calibration information is smoothed over many scanlines. In Pygac, this
+inconsistency is handled properly while calibrating thermal channels.
+
+In some cases it was found that, apart from the reset values, even the
+readings from any one of the four PRTs could also have very low suspicious
+values. This could also seriously affect the computation of brightness
+temperatures. Pygac detects such anomalies and corrects them using
+interpolation of nearby valid PRT readings.
+
+
+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. 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 preferred
+whenever possible.
+
+.. _PyOrbital:
+    https://pyorbital.readthedocs.io
+
+
+Computation of Angles
+---------------------
+
+The azimuth angles are calculated using `get_alt_az`_ and `get_observer_look`_
+from pyorbital. 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.
+
+
+.. _get_alt_az:
+    https://pyorbital.readthedocs.io/en/latest/#pyorbital.astronomy.get_alt_az
+.. _get_observer_look:
+    https://pyorbital.readthedocs.io/en/latest/#pyorbital.orbital.Orbital.get_observer_look
+
+
+Correction of Satellite Location
+--------------------------------
+
+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.
+
+
+Correction of Scanline Timestamps
+---------------------------------
+
+The geolocation in Pygac depends on accurate scanline timestamps. However,
+these may be corrupt, especially for older sensors. Assuming a constant
+scanning rate, Pygac attempts to fix them using extrapolation based on the scan
+line number and a reference time.
+
+Finding the right reference time is difficult due to the multitude of
+possible timestamp corruptions. But the combination of the following three
+options proved to be a robust reference in many situations:
+Timestamp of the first scanline, median time offset of all scanlines and header
+timestamp. See
+:meth:`pygac.reader.Reader.correct_times_median` and
+:meth:`pygac.reader.Reader.correct_times_thresh`
+for details.
+
+Finally, not only timestamps but also scanline numbers may be corrupt.
+Therefor lines with erroneous scanline numbers are removed before
+extrapolation, see :meth:`pygac.reader.Reader.correct_scan_line_numbers`.
+
+
+Scan-Motor-Issue
+----------------
+
+Between 2001 and 2004 GAC data from NOAA-14, NOAA-15, and NOAA-16 frequently
+contain a significant amount of noise towards an edge of the swath. As
+reported by `Schlundt et al (2017)`_, section 5.2, this is probably caused by a
+temporary scan-motor issue. Pygac tries to identify and mask affected pixels.
+
+.. _Schlundt et al (2017):
+    https://climate.esa.int/media/documents/Cloud_Technical-Report-AVHRR-GAC-FCDR-generation_v1.0.pdf
diff --git a/doc/source/usage.rst b/doc/source/usage.rst
index 4aa2532..48dc0f0 100644
--- a/doc/source/usage.rst
+++ b/doc/source/usage.rst
@@ -1,35 +1,106 @@
 Usage
 -----
 
-Standalone
+Through Satpy
+~~~~~~~~~~~~~
+
+The preferred way of using Pygac is through `Satpy`_. Results are returned as
+dask-friendly `xarray`_ DataArrays with proper dataset/coordinate names and
+additional metadata. It is also possible to select a user-defined range of
+scanlines. Furthermore, Satpy provides many options for resampling,
+visualizing and saving the data.
+
+.. code-block:: python
+
+    import satpy
+
+    # Channel set for KLM satellites. For POD satellites the channels are
+    # ['1', '2', '3', '4', '5'].
+    channels = ['1', '2', '3a', '3b', '4', '5']
+    ancillary = ['solar_zenith_angle',
+                 'sensor_zenith_angle',
+                 'solar_azimuth_angle',
+                 'sensor_azimuth_angle',
+                 'sun_sensor_azimuth_difference_angle',
+                 'qual_flags',
+                 'latitude',
+                 'longitude']
+
+    scene = satpy.Scene(filenames=['NSS.GHRR.NP.D15361.S0121.E0315.B3547172.SV'],
+                        reader='avhrr_l1b_gaclac',
+                        reader_kwargs={'tle_dir': '/path/to/tle/',
+                                       'tle_name': 'TLE_%(satname)s.txt'})
+    scene.load(channels + ancillary)
+
+
+For a list of Satpy reader keyword arguments see `satpy.readers.avhrr_l1b_gaclac`_
+and for further Pygac reader keyword arguments see :class:`pygac.reader.Reader`.
+Especially it is possible to choose a different version of calibration
+coefficients or even specify your own.
+
+.. _Satpy: https://satpy.readthedocs.io
+.. _xarray: https://xarray.pydata.org
+.. _satpy.readers.avhrr_l1b_gaclac:
+    https://satpy.readthedocs.io/en/stable/api/satpy.readers.avhrr_l1b_gaclac.html?highlight=avhrr_l1b_gaclac
+.. _example notebook:
+    https://github.com/pytroll/pytroll-examples/blob/main/satpy/avhrr_l1b_gaclac.ipynb
+
+
+Direct Usage
+~~~~~~~~~~~~
+
+Alternatively you can also use Pygac directly.
+
+.. code-block:: python
+
+    from pygac import get_reader_cls
+    
+    filename = 'NSS.GHRR.NP.D15361.S0121.E0315.B3547172.SV'
+    reader_cls = get_reader_cls(filename)
+    reader = reader_cls(tle_dir='/path/to/tle', tle_name='TLE_%(satname)s.txt')
+    reader.read(filename)
+
+    channels = reader.get_calibrated_channels()
+    lons, lats = reader.get_lonlat()
+    scanline_times = reader.get_times()
+    bad_quality_lines = reader.mask
+
+
+Legacy CLI
 ~~~~~~~~~~
 
+.. note::
+
+    Usage of the legacy command line program ``pygac-run`` is deprecated in
+    favour of the above options.
+
+There is also a legacy command line program ``pygac-run`` which saves the
+results to HDF5 and requires a configuration file.
+
 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.::
+it in a directory as you please. Set the environment variable ``PYGAC_CONFIG_FILE``
+pointing to the file. e.g.
+
+.. code-block:: bash
  
-  $> PYGAC_CONFIG_FILE=/home/user/pygac.cfg; export PYGAC_CONFIG_FILE
+  PYGAC_CONFIG_FILE=/home/user/pygac.cfg; export PYGAC_CONFIG_FILE
 
 Also adapt the configuration file to your needs. The ``tledir`` parameter should
 be set to where your Two Line Element (TLE) files are located.
 
-The main script is ``pygac-run``. It automatically checks for the type of file
-format and invokes either ``gac_pod.py`` (POD family, up to and including NOAA-14)
-or ``gac_klm.py`` (KLM family, NOAA-15 and onwards including Metop-A and -B).
+Then call ``pygac-run`` on a GAC/LAC file.
 
-You can test it directly on the testdata included in the package. The result
-will be three hdf5 files, one with the calibrated AVHRR data, the other with
-sun-satellite viewing geometry data and this third with scanline quality
-information::
+.. code-block:: bash
 
- $> pygac-run testdata/NSS.GHRR.NL.D02187.S1904.E2058.B0921517.GC 0 0
+  pygac-run testdata/NSS.GHRR.NL.D02187.S1904.E2058.B0921517.GC 0 0
  
 The last two digits are the start and end scanline numbers, thus specifying the
-portion of the GAC orbit that user wants to process.  The first scanline number
+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/main/satpy/avhrr_l1b_gaclac.ipynb>`_.
+The result will be three hdf5 files, one with the calibrated AVHRR data,
+the other with sun-satellite viewing geometry data and this third with
+scanline quality information.
+
+
diff --git a/pygac/gac_klm.py b/pygac/gac_klm.py
index 293ac1f..081da8d 100644
--- a/pygac/gac_klm.py
+++ b/pygac/gac_klm.py
@@ -22,13 +22,7 @@
 # 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.
-
-Reads L1b GAC data from KLM series of satellites (NOAA-15 and later) and does most of the computations.
-Format specification can be found here:
-http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c8/sec83142-1.htm
-
-"""
+"""Reader for GAC KLM data."""
 
 from __future__ import print_function
 
diff --git a/pygac/gac_pod.py b/pygac/gac_pod.py
index 9cabb56..7700229 100644
--- a/pygac/gac_pod.py
+++ b/pygac/gac_pod.py
@@ -24,12 +24,7 @@
 
 # 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 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
-"""
+"""Reader for GAC POD data."""
 
 import logging
 
diff --git a/pygac/klm_reader.py b/pygac/klm_reader.py
index b927251..a0c5297 100644
--- a/pygac/klm_reader.py
+++ b/pygac/klm_reader.py
@@ -26,9 +26,12 @@
 
 """Read KLM data.
 
-Reads L1b GAC/LAC data from KLM series of satellites (NOAA-15 and later) and does most of the computations.
-Format specification can be found here:
-http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c8/sec83142-1.htm
+Reads L1b GAC/LAC data from KLM series of satellites (NOAA-15 and later).
+Format specification can be found in section 8 of the `KLM user guide`_.
+
+.. _KLM user guide:
+    https://www.ncei.noaa.gov/pub/data/satellite/publications/podguides/N-15%20thru%20N-19/
+
 """
 
 import datetime
@@ -51,21 +54,26 @@
 class KLM_QualityIndicator(IntFlag):
     """Quality Indicators.
 
-    Source:
-        KLM guide
-        Table 8.3.1.3.3.1-1. Format of packed LAC/HRPT Data Sets (Version 2, pre-April 28, 2005).
-        Table 8.3.1.3.3.2-1. Format of LAC/HRPT Data Record for NOAA-N (Version 5, post-November 14,
-                             2006, all spacecraft).
-        Table 8.3.1.4.3.1-1. Format of packed GAC Data Record for NOAA KLM (Version 2, pre-April 28, 2005).
-        Table 8.3.1.4.3.2-1. Format of GAC Data Record for NOAA-N (Version 4, post-January 25, 2006,
-                             all spacecraft).
-
-    Note:
-        Table 8.3.1.3.3.1-1. and Table 8.3.1.4.3.1-1. define bit: 21 as
-        "frame sync word not valid"
-        Table 8.3.1.3.3.2-1. and Table 8.3.1.4.3.2-1. define bit: 21 as
-        "flywheeling detected during this frame"
+    Source: KLM guide
+
+    - Table 8.3.1.3.3.1-1. Format of packed LAC/HRPT Data Sets (Version 2,
+      pre-April 28, 2005).
+    - Table 8.3.1.3.3.2-1. Format of LAC/HRPT Data Record for NOAA-N
+      (Version 5, post-November 14, 2006, all spacecraft).
+    - Table 8.3.1.4.3.1-1. Format of packed GAC Data Record for NOAA KLM
+      (Version 2, pre-April 28, 2005).
+    - Table 8.3.1.4.3.2-1. Format of GAC Data Record for NOAA-N (Version 4,
+      post-January 25, 2006, all spacecraft).
+
+    Notes:
+
+    - Table 8.3.1.3.3.1-1. and Table 8.3.1.4.3.1-1. define bit: 21 as
+      "frame sync word not valid"
+    - Table 8.3.1.3.3.2-1. and Table 8.3.1.4.3.2-1. define bit: 21 as
+      "flywheeling detected during this frame"
+
     """
+
     FATAL_FLAG = 2**31  # Data should not be used for product generation
     TIME_ERROR = 2**30  # Time sequence error detected within this scan
     DATA_GAP = 2**29  # Data gap precedes this scan
diff --git a/pygac/lac_pod.py b/pygac/lac_pod.py
index 760e09b..49f6103 100644
--- a/pygac/lac_pod.py
+++ b/pygac/lac_pod.py
@@ -22,7 +22,7 @@
 # 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."""
+"""Reader for LAC POD data."""
 
 import logging
 
diff --git a/pygac/pod_reader.py b/pygac/pod_reader.py
index 5b28471..ac58ec0 100644
--- a/pygac/pod_reader.py
+++ b/pygac/pod_reader.py
@@ -29,9 +29,11 @@
 
 """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
+Reads L1b GAC/LAC data from POD series of satellites (NOAA-14 and earlier).
+Format specification can be found in chapters 2 & 3 of the `POD user guide`_.
+
+.. _POD user guide:
+    https://www.ncei.noaa.gov/pub/data/satellite/publications/podguides/TIROS-N%20thru%20N-14/
 """
 
 import datetime
diff --git a/pygac/reader.py b/pygac/reader.py
index 90812c1..5b81265 100644
--- a/pygac/reader.py
+++ b/pygac/reader.py
@@ -85,7 +85,7 @@ class NoTLEData(IndexError):
 
 
 class Reader(six.with_metaclass(ABCMeta)):
-    """Reader for Gac and Lac, POD and KLM data."""
+    """Reader for GAC and LAC format, POD and KLM platforms."""
 
     # data set header format, see _validate_header for more details
     data_set_pattern = re.compile(
@@ -166,7 +166,7 @@ def calibration(self):
         return calibration
 
     @abstractmethod
-    def read(self, filename, fileobj=None):
+    def read(self, filename, fileobj=None):  # pragma: no cover
         """Read the GAC/LAC data.
 
         Args:
@@ -177,7 +177,7 @@ def read(self, filename, fileobj=None):
 
     @classmethod
     @abstractmethod
-    def read_header(cls, filename, fileobj=None):
+    def read_header(cls, filename, fileobj=None):  # pragma: no cover
         """Read the file header.
 
         Args:
@@ -354,7 +354,7 @@ def save(self, start_line, end_line, output_file_prefix="PyGAC", output_dir="./"
         )
 
     @abstractmethod
-    def get_header_timestamp(self):
+    def get_header_timestamp(self):  # pragma: no cover
         """Read start timestamp from the header.
 
         Returns:
@@ -400,7 +400,7 @@ def get_counts(self):
             return channels
 
     @abstractmethod
-    def _get_times(self):
+    def _get_times(self):  # pragma: no cover
         """Specify how to read scanline timestamps from GAC data.
 
         Returns:
@@ -500,8 +500,8 @@ def update_meta_data(self):
     def get_calibrated_channels(self):
         """Calibrate and return the channels."""
         channels = self.get_counts()
-        times = self.times
         self.get_times()
+        times = self.times
         self.update_meta_data()
         year = times[0].year
         delta = times[0].date() - datetime.date(year, 1, 1)
@@ -545,6 +545,11 @@ def get_calibrated_channels(self):
 
         return channels
 
+    @abstractmethod
+    def get_telemetry(self):  # pragma: no cover
+        """KLM/POD specific readout of telemetry."""
+        raise NotImplementedError
+
     def get_lonlat(self):
         """Compute lat/lon coordinates.
 
@@ -574,7 +579,7 @@ def get_lonlat(self):
         return self.lons, self.lats
 
     @abstractmethod
-    def _get_lonlat(self):
+    def _get_lonlat(self):  # pragma: no cover
         """KLM/POD specific readout of lat/lon coordinates."""
         raise NotImplementedError
 
@@ -585,6 +590,15 @@ def mask(self):
             self._mask = self._get_corrupt_mask()
         return self._mask
 
+    @abstractproperty
+    def QFlag(self):  # pragma: no cover
+        """KLM/POD specific quality indicators."""
+        raise NotImplementedError
+
+    @abstractproperty
+    def _quality_indicators_key(self):  # pragma: no cover
+        raise NotImplementedError
+
     def _get_corrupt_mask(self, flags=None):
         """Readout of corrupt scanline mask.
 
@@ -615,12 +629,12 @@ def get_qual_flags(self):
         return qual_flags
 
     @abstractmethod
-    def postproc(self, channels):
+    def postproc(self, channels):  # pragma: no cover
         """Apply KLM/POD specific postprocessing."""
         raise NotImplementedError
 
     @abstractmethod
-    def _adjust_clock_drift(self):
+    def _adjust_clock_drift(self):  # pragma: no cover
         """Adjust clock drift."""
         raise NotImplementedError
 
@@ -761,14 +775,19 @@ def get_angles(self):
         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],
+
+            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]
+            and sensor in range [0, 180]
 
         """
         self.get_times()
@@ -851,7 +870,7 @@ def correct_times_median(self, year, jday, msec):
                 msec0 = np.median(msec - msec_lineno)
                 msec = msec0 + msec_lineno
 
-        return year, jday, msec
+        return year.astype(int), jday.astype(int), msec
 
     def correct_scan_line_numbers(self):
         """Remove scanlines with corrupted scanline numbers.
@@ -859,7 +878,7 @@ def correct_scan_line_numbers(self):
         This includes:
             - Scanline numbers outside the valide range
             - Scanline numbers deviating more than a certain threshold from the
-            ideal case (1,2,3,...N)
+              ideal case (1,2,3,...N)
 
         Example files having corrupt scanline numbers:
             - NSS.GHRR.NJ.D96144.S2000.E2148.B0720102.GC
@@ -1031,7 +1050,7 @@ def correct_times_thresh(self, max_diff_from_t0_head=6*60*1000,
         return results
 
     @abstractproperty
-    def tsm_affected_intervals(self):
+    def tsm_affected_intervals(self):  # pragma: no cover
         """Specify time intervals being affected by the scan motor problem.
 
         Returns:
@@ -1104,7 +1123,7 @@ def mask_tsm_pixels(self, channels):
         channels[idx] = np.nan
 
     @abstractmethod
-    def get_tsm_pixels(self, channels):
+    def get_tsm_pixels(self, channels):  # pragma: no cover
         """Determine pixels affected by the scan motor issue.
 
         Channel selection is POD/KLM specific.
diff --git a/setup.py b/setup.py
index f0f00d6..5174ed8 100644
--- a/setup.py
+++ b/setup.py
@@ -43,6 +43,10 @@
                     'scipy>=0.8.0',
                     'python-geotiepoints>=1.1.8',
                     'bottleneck>=1.0.0']
+    extras_require = {
+        'dev': ['pytest', 'pre-commit', 'flake8']
+    }
+
     if sys.version_info < (3, 7):
         # To parse ISO timestamps in calibration.py
         requirements.append('python-dateutil>=2.8.0')
@@ -69,6 +73,7 @@
           # Project should use reStructuredText, so ensure that the docutils get
           # installed or upgraded on the target machine
           install_requires=requirements,
+          extras_require=extras_require,
           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'])],


Test patch:
diff --git a/pygac/tests/test_reader.py b/pygac/tests/test_reader.py
index b13ed91..29388f5 100644
--- a/pygac/tests/test_reader.py
+++ b/pygac/tests/test_reader.py
@@ -31,6 +31,8 @@
 import numpy as np
 import numpy.testing
 from pygac.gac_reader import GACReader, ReaderError
+from pygac.pod_reader import POD_QualityIndicator
+from pygac.gac_pod import scanline
 from pygac.reader import NoTLEData
 
 
@@ -42,6 +44,58 @@ def __fspath__(self):
         return self.path
 
 
+class FakeGACReader(GACReader):
+    QFlag = POD_QualityIndicator
+    _quality_indicators_key = "quality_indicators"
+    tsm_affected_intervals = {None: []}
+    along_track = 3
+    across_track = 4
+
+    def __init__(self):
+        super(FakeGACReader, self).__init__()
+        self.scan_width = self.across_track
+        scans = np.zeros(self.along_track, dtype=scanline)
+        scans["scan_line_number"] = np.arange(self.along_track)
+        scans["sensor_data"] = 128
+        self.scans = scans
+        self.head = {'foo': 'bar'}
+        self.spacecraft_name = 'noaa6'
+
+    def _get_times(self):
+        year = np.full(self.along_track, 1970, dtype=int)
+        jday = np.full(self.along_track, 1, dtype=int)
+        msec = 1000 * np.arange(1, self.along_track+1, dtype=int)
+        return year, jday, msec
+
+    def get_header_timestamp(self):
+        return datetime.datetime(1970, 1, 1)
+
+    def get_telemetry(self):
+        prt = 51 * np.ones(self.along_track)  # prt threshold is 50
+        ict = 101 * np.ones((self.along_track, 3))  # ict threshold is 100
+        space = 101 * np.ones((self.along_track, 3))  # space threshold is 100
+        return prt, ict, space
+
+    def _adjust_clock_drift(self):
+        pass
+
+    def _get_lonlat(self):
+        pass
+
+    def postproc(self, channels):
+        pass
+
+    def read(self, filename, fileobj=None):
+        pass
+
+    @classmethod
+    def read_header(cls, filename, fileobj=None):
+        pass
+
+    def get_tsm_pixels(self, channels):
+        pass
+
+
 class TestGacReader(unittest.TestCase):
     """Test the common GAC Reader."""
 
@@ -137,6 +191,21 @@ def test__get_calibrated_channels_uniform_shape(self, get_channels):
         with self.assertRaises(AssertionError):
             self.reader._get_calibrated_channels_uniform_shape()
 
+    def test_get_calibrated_channels(self):
+        reader = FakeGACReader()
+        res = reader.get_calibrated_channels()
+        expected = np.full(
+            (
+                FakeGACReader.along_track,
+                FakeGACReader.across_track,
+                5  # number of channels
+            ),
+            np.nan
+        )
+        expected[:, 1, 0] = 8.84714652
+        expected[:, 2, 1] = 10.23511303
+        np.testing.assert_allclose(res, expected)
+
     def test_to_datetime64(self):
         """Test conversion from (year, jday, msec) to datetime64."""
         t0 = GACReader.to_datetime64(year=np.array(1970), jday=np.array(1),