Add parallax correction via new ParallaxCorrectionModifier

[gerritholl]

Add a parallax correction to Satpy.

This is a team effort based on work by Simon Proud (@simonrp84), Xin Zhang (@zxdawn), Ulrich Hamann (@uhamann), Jussi Leinonen (@jleinonen), myself, and others.

We have made good progress during the Pytroll Contributor Week 29 November - 3 December 2021. More work was done in December 2021 and January 2022.

• Closes Add parallax correction #1845
• Tests added
• Fully documented
• Add your name to AUTHORS.md if not there already
• Add more tests where
  • source and target resolutions are very different
  • areas are only partly overlapping
  • areas are not overlapping at all
  • areas are the same
• get reliable reference values for unit tests where this value is not trivial
• Investigate if other resamplers resolve the problems at poles and antimeridian
• Make resampler and search radius configurable
• Add easy to use modifier interface
• Get orbital parameters from somewhere else than the metadata of the CTH dataset. This information is not included with nwcsaf-geo and arguably should come from the dataset to be corrected, not from the CTH (which theoretically doesn't even need to come from the same satellite). See Add parallax correction via new ParallaxCorrectionModifier #1904 (comment).
• Make sure heights are in metre
• Add test to cover the case where clouds are moved out of the area of interest completely.
• Improve speed (although there are no formal performance tests yet, small tests show the correction is very slow).
• Fix bug when loading corrected and uncorrected datasets simultaneously.
• Check why brightness is different for parallax corrected image and fix this.
• Make internal resampling of CTH more flexible (rather than hardcoding method and parameters)
• Make resampling in modifier more flexible, actually using the provided parameters.

Deferred:

• Not critical: implement checks for warning or error if areas do not or only partially overlap. This relies on intersection calculation never completes, infinite loop with 100% CPU pyresample#329 to be fixed first. Without those checks, applying parallax correction may yield unexpected results.
• Add diagnostic tool to get the displacement (distance & direction)
• Expose and document interface to intermediate swathdefinition as public API
• Add test for when parallax correction shifts a cloud across the antimeridian or pole (or document that this is not supported)
• Allow explicitly passing satellite location for datasets where this information is not included (either when no information provided at all, or when estimate based on TLE and skyfield is undesirable)
• Allow obtaining satellite location or viewing angles from the metadata not only of CTH but also of other datasets. This would require those to be passed to the modifier as optional prerequisites.
• Add performance tests (in benchmark directory?) with real scenes (could be a different PR)
  • Check dask-awareness
  • Check/control for what area parallax is calculated when scene is resampled, but modifier interface is used
  • Confirm that dask reuses calculations, such as when multiple channels are corrected
  • Check if further optimisations are needed for multiple channels, such as avoiding recalculating the dask graphs

[gerritholl]

Test to calculate parallax correction for New England, 2021-12-29 20:01:17Z.

Uncorrected RGB:

Corrected RGB:

Uncorrected C13:

Corrected C13:

Difference C13:

Code for test:

[gerritholl]

My understanding so far.

Consider:

• For location A", the satellite reports a cloud that is actually at location A. The location A" is actually clear-sky, but we have no way of knowing that.
• For location A, the satellite reports a cloud that is actually at location A'.
• For location A', the satellite observes clear-sky, but location A' is actually cloudy.

The forward parallax correction can calculate that for position A", a cloud seen at height h is in reality position A, and that at location A, a cloud at height h is in reality at location A'. We can create a SwathDefinition that assigns a new lat/lon for each and use that. Or we can — or must? — do another step and reassign the pixel values. Meaning:

• For location A, we know that this is the actual location of the cloud reported at location A" — so at this location, we assign the pixel value that was reported for location A".
• Likewise, for location A', we know that A' corresponds to the true location of the pixel reported at location A, so here we assign the pixel value that was reported for location A.

This reassignment of pixels, if I understand correctly, is done by this PR by calculating an inverse transformation using a convolution (I don't understand the details yet here), resampling to this, and then reassigning the original base area:

corrector = ParallaxCorrection(base_area)
plax_corr_area = corrector(cloud_top_height)  # SwathDefinition of same size as base_area
local_scene = global_scene.resample(plax_corr_area)
local_scene[some_dataset].attrs["area"] = base_area

This somewhat abuses the concept of resampling and I don't understand it fully yet, but it appears to work for a high resolution area, as the images in the previous comment show...

Thoughts welcome!

[gerritholl]

I selected a clear sky area in the range 71.9-73.2°W, 43.6–45.0°N, and plotted the difference in channel 13 before and after applying parallax correction in the way I understand it to be correct (see previous comment). This exhibits unphysical artefacts with differences of up to 4–5 K:

[jleinonen]

I pushed an alternative implementation as parallax2.py in order to avoid merge conflicts. If we decide to use it then we should eventually replace parallax.py with it.

Changes:

• Eliminate the dependence on projection.py (this too can be deleted if we adopt this version)
• Replace the RBF interpolation with NumpyBilinearResampler from Pyresample
• Do the interpolation on the correction to (lon,lat) instead of (lon,lat) itself

Things remaining to do/consider:

• The radius_of _influence parameter to the resampler is hardcoded - can/should we detect this automatically?
• Testing different scenarios: source resolution much higher/lower than the target resolution, source and target areas not fully overlapping
• Testing on real scenes and evaluating accuracy and performance
• Should we replace NumpyBilinearResampler with XArrayBilinearResampler?
• What about nearest neighbor interpolation?

[gerritholl]

With the latest update, the Δ looks much better:

Note the change in the scale. The high values near the bottom may be due to a cloud that is just outside the field of view.

[gerritholl]

We also need to test what happens at the antimeridian. And at the poles? Although most critical for geostationary applications in (northern) Europe, people may want to use it for polar orbiters too.

[mraspaud]

I just have three minor comments (on top of the one left from before), then I think we can merge. Great job!

[mraspaud]

LGTM

[ghiggi]

Hey @gerritholl. Here are the other comments. I forgot to submit the review :P

[djhoese]

LGTM

[djhoese]

It. is. done. 🎉

[mraspaud]

Thank you @gerritholl for driving this to its completion, and thanks to everyone involved for taking the time to make this PR better!