Refactor clock drift

[carloshorn]

This PR closes #76 and closes #80 and the associated PRs: closes #81 and closes #82.

Enhancements

1. By adding more comment lines and proper variable names, the code should be more comprehensible.
2. By making reader.times a property, there is no need to synchronize the array with reader.utcs anymore. This makes sense, because both arrays contain the same information but expressed in different data types.

Fixes

1. It turned out that #76 was a consequence of an error in computing the maximum line number of shifted line indices. https://github.com/pytroll/pygac/blob/5c59f2be9437ad7bfc0005038f45f8416d10828d/pygac/pod_reader.py#L452 The - max_idx is wrong, see discussion in #82.

2. This PR takes care of skipped scan lines which addresses the following TODO: https://github.com/pytroll/pygac/blob/5c59f2be9437ad7bfc0005038f45f8416d10828d/pygac/pod_reader.py#L418 by taking into account the complete line number range when sub-setting the missed lines, instead of only considering the shifted lines as in the current implementation: https://github.com/pytroll/pygac/blob/5c59f2be9437ad7bfc0005038f45f8416d10828d/pygac/pod_reader.py#L446-L448

3. The current implementation computes wrong utc values for missed scan lines if the first scan line number is not equal to one, see https://github.com/pytroll/pygac/blob/5c59f2be9437ad7bfc0005038f45f8416d10828d/pygac/pod_reader.py#L463-L464 The correct time transformation is:

   missed_utcs = ((missed - self.scans["scan_line_number"][0])*np.timedelta64(500, "ms")
              + self.utcs[0])

4. The computed longitudes and latitudes for the slerp interpolation should correspond to the utcs before clock adjustment. So the time difference between each line in the interpolation array is equal to the scan rate of 500 ms. Therefore, it should be clock_drift_adjust=False in the following call: https://github.com/pytroll/pygac/blob/5c59f2be9437ad7bfc0005038f45f8416d10828d/pygac/pod_reader.py#L466-L468 otherwise, the values in complete_lons, complete_lats are not necessarily in the correct order.

5. It corrects #80, by making reader.times a property and taking into account milliseconds in the final assignment https://github.com/pytroll/pygac/blob/5c59f2be9437ad7bfc0005038f45f8416d10828d/pygac/pod_reader.py#L487 Furthermore, note that this line is also incorrect, because after https://github.com/pytroll/pygac/blob/5c59f2be9437ad7bfc0005038f45f8416d10828d/pygac/pod_reader.py#L438 the variable offsets has the wrong units (1/scan_rate = 1/(0.5 seconds/line) = 2 lines/seconds)

Potential Side Effects:

1. This PR removes the method Reader.compute_lonlat, see https://github.com/pytroll/pygac/blob/5c59f2be9437ad7bfc0005038f45f8416d10828d/pygac/reader.py#L471 It is not considered as necessary for the public interface of Reader, because it is only called for clock drift adjustment, where the calling method reader._adjust_clock_drift is private already.

[carloshorn]

I noticed that there were not tests yet regarding the clock drift adjustment. I would postpone the creation of tests after the first code review.

[carloshorn]

I would have liked to see tests for the changed functionality or bugfixing before the refactoring took place, it would ensure that the refactoring is just a refactoring, ie that all the tests are still passing.

Good point for the future, I am not yet much into TDD, but I clearly see the benefit. How can we proceed on it? I could create a new branch where I implement the test and include the fixes and changes step by step. Or is the now implemented test meaningful enough to show the correctness of this PR?

[mraspaud]

@carloshorn looks go with the test you added! I don't think you need to rewind your work now.

As far as TDD goes, it's just as the name says :) Write a test first (that fails) and then proceed with the implementation so that the test passes. When the test passes, refactor so that the code becomes clean. See a demo of how to do this here: https://www.youtube.com/watch?v=58jGpV2Cg50&t=2628s

The hardest part is to know what to test. My rule of thumb is that every line should be tested, but we test only public functionality, not implementation details. We can talk more on slack about it if you want (in #random for example).

[carloshorn]

According to Table J-1 (Satellite Scanning Instrument Parameters for TIROS-N through NOAA-14) in the KLM guide, the max scan angle for AVHRR/2 is given by 55.37 deg. I don't know where the 55.385 deg in the code come from, but the 55.37 deg is also the default in pyorbita..geoloc_instrument_definitions.avhrr_gac.

I will have to check, but I think this was calculated somehow. It's somewhere in pyorbital maybe.

I am still trying to understand the given scan angle of 55.385 deg which is passed to pyorbital in the current implementation: https://github.com/pytroll/pygac/blob/5c59f2be9437ad7bfc0005038f45f8416d10828d/pygac/reader.py#L523-L524

Unfortunately, I could not find any clear hint in pyorbital. My first guess, that it is the scan midpoint plus half a scan step, did not lead to this number. Then, the following line got my attention: https://github.com/pytroll/pyorbital/blob/42768fb839f22fe05254bb98aa8a774199462a2a/pyorbital/geoloc.py#L284 But still, I could not reproduce the 55.385 degrees.

Can you think of anything else to try, or maybe some verbal ideas/explanations why the scan angle needs an adjustment?

[ninahakansson]

@carloshorn Nice work with the PR! About the 55.385 this might be inherited from the PPS-ahamap package. I found a note from 2002. I can ask KG if he remembers :).

   /*    Replaced with calculations taken from 6S Users' Guide Version 2 /KG */

   for(y=0;y<self->podlac->podHeader.numberOfScans;y++) {
      for(x=0;x<LAC_MXPIX;x++) {
     scan_angle=55.385*abs(x-1024)/1024;
     satzenang[y*LAC_MXPIX+x]=asin((1+sat_alt/earth_radius)*sin(scan_angle*DEG_TO_RAD))*RAD_TO_DEG;
      }
   }

[carloshorn]

@ninahakansson, please ask KG. I hope he remembers, otherwise, I would vote for using the default value 55.37 from the KLM Guide, because I would prefer to have a reference/justification/explanation for any assumption/choice/parameter that appears in the code. If someone remembers or comes up with an explanation later on, we can still open an issue, change back to 55.385 and add a comment line with the explanation. What do you think?

[mraspaud]

From what I can find, pyorbital is also defaulting to 55.37 for avhrr: https://github.com/pytroll/pyorbital/blob/master/pyorbital/geoloc_instrument_definitions.py#L52-L102

[mraspaud]

Also notice that the scan angle is 55.25 for noaa 16: https://www.ospo.noaa.gov/data/ppp/notice_files/mar0101.html

[ninahakansson]

KG had nothing to add and @abhaydd showed me that in one part of the KLM-guide that value is 55.4 and that the value varies between sources, versions of sources and actually between instruments. I have no opinion about updating it or not. As it is a small change PPS will be OK either way.

[carloshorn]

Also notice that the scan angle is 55.25 for noaa 16: https://www.ospo.noaa.gov/data/ppp/notice_files/mar0101.html

Currently, there is just a clock drift adjustment for POD satellites, but it is truly a thing to remember when adding clock drifts for KLM satellites in the future.

[carloshorn]

@ninahakansson, thanks for checking. In this case, I would stick on one source, which could be Table J-1 (Satellite Scanning Instrument Parameters for TIROS-N through NOAA-14) in the KLM guide and use it to justify the choice in the code. In the current PR, this is done indirectly by passing the responsibility to the default value of pyorbital.geoloc_instrument_definitions.avhrr_gac, which refers to exactly this table in its doc-string.

[mraspaud]

@carloshorn in principle clock drift adjustment shouldn't be needed as the clock is updated daily on the klm series.

[carloshorn]

@mraspaud, cool, then we are lucky and do not need to implement an extra case for NOAA-16.

[mraspaud]

@sfinkens I think you need to run your tests on this one

[ghost]

Congratulations :tada:. DeepCode analyzed your code in 2.218 seconds and we found no issues. Enjoy a moment of no bugs :sunny:.

👉 View analysis in DeepCode’s Dashboard | Configure the bot

[carloshorn]

@mraspaud, sorry I have not noticed that you were about to approving this. I went through the code a last time for my-self.

[mraspaud]

Retry DeepCode

[sfinkens]

Awesome, now I'm able to understand what's happening :smile: Thanks for digging into this! I'll run the tests ASAP.

[mraspaud]

@sfinkens did you run the tests on this one too?

[sfinkens]

I'm on it. There are differences, as expected, but I think not critical. Currently hunting a memory leak occuring while plotting the differences...

[sfinkens]

Regression tests show minor differences for the angles, which is expected because the scanline timestamps were corrected.

[mraspaud]

Thanks for the hard work everyone!