Refactor clock drift #84

carloshorn · 2020-09-02T15:45:01Z

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

Enhancements

By adding more comment lines and proper variable names, the code should be more comprehensible.
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

It turned out that IndexError in PODReader._adjust_clock_drift #76 was a consequence of an error in computing the maximum line number of shifted line indices.

pygac/pygac/pod_reader.py

Line 452 in 5c59f2b

max(self.scans["scan_line_number"] - min_idx)) + 1

The - max_idx is wrong, see discussion in Fix clock drift rounding error #82.

This PR takes care of skipped scan lines which addresses the following TODO:

pygac/pygac/pod_reader.py

Line 418 in 5c59f2b

TODO: bad things might happen when scanlines are skipped.

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:

pygac/pygac/pod_reader.py

Lines 446 to 448 in 5c59f2b

    
           missed = sorted((set(line_indices) | 
        
                            set(line_indices + 1)) 
        
                           - set(self.scans["scan_line_number"]))

The current implementation computes wrong utc values for missed scan lines if the first scan line number is not equal to one, see

pygac/pygac/pod_reader.py

Lines 463 to 464 in 5c59f2b

missed_utcs = ((np.array(missed) - 1) * np.timedelta64(500, "ms")

+ self.utcs[0])

The correct time transformation is:

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

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:

pygac/pygac/pod_reader.py

Lines 466 to 468 in 5c59f2b

    
           mlons, mlats = self.compute_lonlat(width=self.lats.shape[1], 
        
                                              utcs=missed_utcs, 
        
                                              clock_drift_adjust=True)

otherwise, the values in complete_lons, complete_lats are not necessarily in the correct order.

It corrects Reduce rounding error in POD reader adjust clock drift #80, by making reader.times a property and taking into account milliseconds in the final assignment

pygac/pygac/pod_reader.py

Line 487 in 5c59f2b

self.utcs += offsets.astype('timedelta64[s]')

Furthermore, note that this line is also incorrect, because after

pygac/pygac/pod_reader.py

Line 438 in 5c59f2b

offsets *= -2

the variable offsets has the wrong units (1/scan_rate = 1/(0.5 seconds/line) = 2 lines/seconds)

Potential Side Effects:

This PR removes the method Reader.compute_lonlat, see

pygac/pygac/reader.py

Line 471 in 5c59f2b

def compute_lonlat(self, width, utcs=None, clock_drift_adjust=True):

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 · 2020-09-02T15:52:55Z

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.

mraspaud

Thanks for starting this. 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.
No need to delay tests any further now though, go ahead.

mraspaud · 2020-09-02T16:54:44Z

pygac/pod_reader.py

+        # TODO: Are we sure all satellites have this scan width in degrees ?
+        # TODO: Does this also work for LAC?


Good questions, they should be answered before we merge this.

Are we sure all satellites have this scan width in degrees ?

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.

Does this also work for LAC?

No, according to Table 8.3.1.3.1-1. (LAC/HRPT Data Characteristics.) the scan rate for LAC/HRPT is 360 scans per minute => 1/6 seconds/line instead of the 1/2 seconds/line which is currently assumed in PODReader._adjust_clock_drift. This can be fixed by making the scan rate a LAC/GAC reader attribute. This should be passed as frequency argument to pyorbita..geoloc_instrument_definitions.avhrr_gac.

Just noticed that the GAC/LAC readers already have the attribute reader.scan_freq... I will definitely use it!

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.

mraspaud · 2020-09-02T17:18:48Z

pygac/reader.py

-        return (((year - 1970).astype('datetime64[Y]')
-                 + (jday - 1).astype('timedelta64[D]')).astype('datetime64[ms]')
+        return (year.astype(str).astype('datetime64[Y]')


Could you explain the logic here?

I usually try to avoid hard coded numbers as much as possible. In this case, I did not like that PyGAC hard codes the numpy.datetime64 epoch (1970). I know that there is a very low chance that numpy will ever change the epoch, but still I think PyGAC should not make any assumptions on it.
The shortest way I found to become independent of the epoch was the string conversion, because something like epoch = int(str(np.datetime64(0, 'Y'))) looks ugly to me and numpy also does not have any attribute to access its epoch directly (or at least I could not easily find it).

Makes sense, thanks for explaining!

carloshorn · 2020-09-03T13:01:20Z

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 · 2020-09-03T15:34:40Z

@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 · 2020-09-04T09:17:18Z

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:

pygac/pygac/reader.py

Lines 523 to 524 in 5c59f2b

    
           sgeom = avhrr_gac(utcs.astype(datetime.datetime), 
        
                             self.scan_points, 55.385)

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 · 2020-09-04T19:08:22Z

@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 · 2020-09-07T07:25:47Z

@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 · 2020-09-07T10:29:39Z

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 · 2020-09-07T10:42:18Z

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

ninahakansson · 2020-09-07T10:56:37Z

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 · 2020-09-07T10:58:00Z

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 · 2020-09-07T11:02:54Z

@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 · 2020-09-07T11:05:16Z

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

carloshorn · 2020-09-07T11:08:50Z

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

mraspaud

LGTM, thanks for digging into this.

…tions

mraspaud · 2020-09-07T12:22:36Z

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

ghost · 2020-09-07T12:23:59Z

Congratulations 🎉. DeepCode analyzed your code in 2.218 seconds and we found no issues. Enjoy a moment of no bugs ☀️.

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

carloshorn · 2020-09-07T12:24:38Z

@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 · 2020-09-07T12:29:41Z

Retry DeepCode

sfinkens · 2020-09-08T10:14:16Z

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

mraspaud · 2020-09-11T08:02:36Z

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

sfinkens · 2020-09-11T09:27:49Z

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 · 2020-09-11T12:59:47Z

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

mraspaud · 2020-09-11T13:03:36Z

Thanks for the hard work everyone!

horn and others added 6 commits September 1, 2020 13:30

convert Reader.times into a property

6b2f556

refactored adjust clock drift logic

a76df6e

fix undefined variable

b777678

refactor adjust clock drift code layout

917e417

Fixing style errors.

d801dcb

fix variable names

5216187

carloshorn requested a review from mraspaud September 2, 2020 15:53

carloshorn self-assigned this Sep 2, 2020

carloshorn added bug enhancement labels Sep 2, 2020

mraspaud reviewed Sep 2, 2020

View reviewed changes

fix floor call in POD reader and add test

5f7d2fc

LAC/HRPT compatibility.

5c34bee

Understood why clock drift adjustment is not needed for KLM

856b23e

mraspaud approved these changes Sep 7, 2020

View reviewed changes

removed unexplained try except clause on calculating missing geo loca…

3f04306

…tions

sfinkens approved these changes Sep 11, 2020

View reviewed changes

mraspaud merged commit aa841ef into pytroll:master Sep 11, 2020

carloshorn deleted the refactor_clock_drift branch September 11, 2020 13:15

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Refactor clock drift #84

Refactor clock drift #84

carloshorn commented Sep 2, 2020 •

edited by mraspaud

Loading

carloshorn commented Sep 2, 2020

mraspaud left a comment

mraspaud Sep 2, 2020

carloshorn Sep 3, 2020

carloshorn Sep 3, 2020

mraspaud Sep 3, 2020

mraspaud Sep 2, 2020

carloshorn Sep 3, 2020 •

edited

Loading

mraspaud Sep 3, 2020

carloshorn commented Sep 3, 2020

mraspaud commented Sep 3, 2020

carloshorn commented Sep 4, 2020

ninahakansson commented Sep 4, 2020

carloshorn commented Sep 7, 2020

mraspaud commented Sep 7, 2020

mraspaud commented Sep 7, 2020

ninahakansson commented Sep 7, 2020

carloshorn commented Sep 7, 2020

carloshorn commented Sep 7, 2020

mraspaud commented Sep 7, 2020 •

edited

Loading

carloshorn commented Sep 7, 2020

mraspaud left a comment

mraspaud commented Sep 7, 2020

ghost commented Sep 7, 2020 •

edited by ghost

Loading

carloshorn commented Sep 7, 2020

mraspaud commented Sep 7, 2020

sfinkens commented Sep 8, 2020

mraspaud commented Sep 11, 2020

sfinkens commented Sep 11, 2020 •

edited

Loading

sfinkens commented Sep 11, 2020

mraspaud commented Sep 11, 2020

	missed = sorted((set(line_indices) \|
	set(line_indices + 1))
	- set(self.scans["scan_line_number"]))

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

	mlons, mlats = self.compute_lonlat(width=self.lats.shape[1],
	utcs=missed_utcs,
	clock_drift_adjust=True)

		# TODO: Are we sure all satellites have this scan width in degrees ?
		# TODO: Does this also work for LAC?

Refactor clock drift #84

Refactor clock drift #84

Conversation

carloshorn commented Sep 2, 2020 • edited by mraspaud Loading

Enhancements

Fixes

Potential Side Effects:

carloshorn commented Sep 2, 2020

mraspaud left a comment

Choose a reason for hiding this comment

mraspaud Sep 2, 2020

Choose a reason for hiding this comment

carloshorn Sep 3, 2020

Choose a reason for hiding this comment

carloshorn Sep 3, 2020

Choose a reason for hiding this comment

mraspaud Sep 3, 2020

Choose a reason for hiding this comment

mraspaud Sep 2, 2020

Choose a reason for hiding this comment

carloshorn Sep 3, 2020 • edited Loading

Choose a reason for hiding this comment

mraspaud Sep 3, 2020

Choose a reason for hiding this comment

carloshorn commented Sep 3, 2020

mraspaud commented Sep 3, 2020

carloshorn commented Sep 4, 2020

ninahakansson commented Sep 4, 2020

carloshorn commented Sep 7, 2020

mraspaud commented Sep 7, 2020

mraspaud commented Sep 7, 2020

ninahakansson commented Sep 7, 2020

carloshorn commented Sep 7, 2020

carloshorn commented Sep 7, 2020

mraspaud commented Sep 7, 2020 • edited Loading

carloshorn commented Sep 7, 2020

mraspaud left a comment

Choose a reason for hiding this comment

mraspaud commented Sep 7, 2020

ghost commented Sep 7, 2020 • edited by ghost Loading

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

carloshorn commented Sep 7, 2020

mraspaud commented Sep 7, 2020

sfinkens commented Sep 8, 2020

mraspaud commented Sep 11, 2020

sfinkens commented Sep 11, 2020 • edited Loading

sfinkens commented Sep 11, 2020

mraspaud commented Sep 11, 2020

carloshorn commented Sep 2, 2020 •

edited by mraspaud

Loading

carloshorn Sep 3, 2020 •

edited

Loading

mraspaud commented Sep 7, 2020 •

edited

Loading

ghost commented Sep 7, 2020 •

edited by ghost

Loading

sfinkens commented Sep 11, 2020 •

edited

Loading