Fix clock drift rounding error

[carloshorn]

This PR closes #80.

[carloshorn]

Ohh... I did not create the branch from master.... so it also contains #81.

[mraspaud]

ok, we'll just have to be careful and merge them in the right order

[mraspaud]

why remove min_idx here?

[carloshorn]

because this leads to the index error. The max argument compares line numbers and is therefore in the line number domain. the - min_idx is the transformation to the complete_lons array index domain. In the case of the crashing orbits, we had exactly the situation, where the right side wins the max comparison.

[carloshorn]

Maybe better variable names could help to avoid this confusion. Something like min_line, max_line, line_range = max_line - min_line. Then we can use index or idx for the indices of complete_lons, complete_lats which are numbers in 0:line_range+1, where the interesting ones are self.scans["scan_line_number"] - min_line.

[carloshorn]

because this leads to the index error. The max argument compares line numbers and is therefore in the line number domain. the - min_idx is the transformation to the complete_lons array index domain. In the case of the crashing orbits, we had exactly the situation, where the right side wins the max comparison.

From reading again the code, I see, that the right hand side should win the max comparison whenever the clock drift is positive, because the array line_indices results from

pygac/pygac/pod_reader.py

Lines 438 to 444 in 5c59f2b

	offsets *= -2
	int_offsets = np.floor(offsets).astype(np.int)
	# filling out missing geolocations with computation from pyorbital.
	line_indices = (self.scans["scan_line_number"]
	+ int_offsets)

where offsets is a floating point number in seconds. The multiplication with -2 to convert the time offset into a line number offset is not well documented, but we know that the sign of the offset determines who wins the max comparison. I think the reason why it sill worked in most of the cases is the +2 in the following line:

pygac/pygac/pod_reader.py

Line 453 in 5c59f2b

idx_len = max_idx - min_idx + 2

the line I tested that crashed had an int_offset of -3, which in combination with the unexplained +2 appeared suspicious to me. Since the max index has already the +1, the range of valid indices should not require any +2. My first guess was that the +2 is required for the later interpolation, but slerp is linear and therefore only requires the next neighbours.

[carloshorn]

Now, I assume that it also requires this change to cover all possible indexes

min_idx = min(line_indices.min(), self.scans["scan_line_number"].min())

[carloshorn]

And without understanding the transformation, but arguing with symmetry, I would assume that

pygac/pygac/pod_reader.py

Line 440 in 5c59f2b

int_offsets = np.floor(offsets).astype(np.int)

should be changed to

int_offsets = np.where(offsets<0, np.floor(offsets), np.ceil(offsets)).astype(np.int) 

to cover the right index range, but here I am not yet 100% sure.

[carloshorn]

I need to further understand the underlying code...
It appears like the same offset adjustment is done again in Reader.compute_lonlat

pygac/pygac/reader.py

Lines 485 to 497 in 5c59f2b

	if clock_drift_adjust:
	from pygac.clock_offsets_converter import get_offsets
	try:
	offset_times, clock_error = get_offsets(self.spacecraft_name)
	except KeyError:
	LOG.info("No clock drift info available for %s",
	self.spacecraft_name)
	else:
	offset_times = np.array(offset_times, dtype='datetime64[ms]')
	offsets = np.interp(utcs.astype(np.uint64),
	offset_times.astype(np.uint64),
	clock_error)
	utcs = utcs - (offsets * 1000).astype('timedelta64[ms]')

where in this case there is no rounding error.
So the first get_offsets call is only used to get the index range. Then it is called again to do the actual shifting by the correct offset...

I will prepare a refactoring of the clock drift adjustment, so maybe #81 and #82 become obsolete.
Could you tell me the difference between reader.time and reader.utcs?
It appears to me like it is the same information, but stored in different data types.
Storing it twice has the disadvantage of synchronization. Maybe one of them should become a property which does the type conversion.