Utility function for parallax displacement + docs

Add utility function for parallax displacement, as suggested by @ghiggi.
Rename the function `forward_parallax` to be called
`get_parallax_corrected_lonlats`.
Various smaller edits/changes in documentation and variable names.

doc/source/modifiers.rst

@@ -122,4 +122,7 @@ or incorrect results.  It does not yet perform any checks that the
 provided (cloud top) height covers the area of the dataset for which
 the parallax correction shall be applied.
 
+For more general background information and web routines related to the
+parallax effect, see also `this collection at the CIMSS website <https://cimss.ssec.wisc.edu/goes/webapps/parallax/>_`.
+
 .. versionadded:: 0.37

satpy/modifiers/parallax.py

@@ -46,6 +46,7 @@
 import xarray as xr
 from pyorbital.orbital import A as EARTH_RADIUS
 from pyorbital.orbital import get_observer_look
+from pyproj import Geod
 from pyresample.bucket import BucketResampler
 from pyresample.geometry import SwathDefinition
 
@@ -64,8 +65,8 @@ class IncompleteHeightWarning(UserWarning):
     """Raised when heights only partially overlap with area to be corrected."""
 
 
-def forward_parallax(sat_lon, sat_lat, sat_alt, lon, lat, height):
-    """Calculate forward parallax effect.
+def get_parallax_corrected_lonlats(sat_lon, sat_lat, sat_alt, lon, lat, height):
+    """Calculate parallax corrected lon/lats.
 
     Satellite geolocation generally assumes an unobstructed view of a smooth
     Earth surface.  In reality, this view may be obstructed by clouds or
@@ -84,12 +85,13 @@ def forward_parallax(sat_lon, sat_lat, sat_alt, lon, lat, height):
 
     .. note::
 
-        Be careful with units!  Heights may be either in m or km, and may
-        refer to either the Earth's surface on the Earth's centre.  Cloud top
-        height is usually reported in meters above the Earth's surface, rarely in
-        km.  Satellite altitude may be reported in either m or km, but orbital
-        parameters may be in relation the the Earths centre.  The Earth radius
-        from pyresample is reported in km.
+        Be careful with units!  This code expects ``sat_alt`` and
+        ``height`` to be in meter above the Earth's surface.  You may
+        have to convert your input correspondingly.  Cloud Top Height
+        is usually reported in meters above the Earth's surface, rarely
+        in km.  Satellite altitude may be reported in either m or km, but
+        orbital parameters are usually in relation the the Earth's centre.
+        The Earth radius from pyresample is reported in km.
 
     Args:
         sat_lon (number): Satellite longitude in geodetic coordinates [degrees]
@@ -103,19 +105,36 @@ def forward_parallax(sat_lon, sat_lat, sat_alt, lon, lat, height):
             will be based.  Typically this is the cloud top height. [m]
 
     Returns:
-        tuple[float, float]: New geolocation
-            New geolocation ``(lon, lat)`` for the longitude and
-            latitude that were to be corrected, in geodetic coordinates. [degrees]
+        tuple[float, float]: Corrected geolocation
+            Corrected geolocation ``(lon, lat)`` in geodetic coordinates for
+            the pixel(s) to be corrected, in geodetic coordinates. [degrees]
     """
     elevation = _get_satellite_elevation(sat_lon, sat_lat, sat_alt, lon, lat)
-    parallax_distance = _calculate_parallax_distance(height, elevation)
+    parallax_distance = _calculate_slant_cloud_distance(height, elevation)
     shifted_xyz = _get_parallax_shift_xyz(
             sat_lon, sat_lat, sat_alt, lon, lat, parallax_distance)
 
     return xyz2lonlat(
         shifted_xyz[..., 0], shifted_xyz[..., 1], shifted_xyz[..., 2])
 
 
+def get_surface_parallax_displacement(
+        sat_lon, sat_lat, sat_alt, lon, lat, height):
+    """Calculate surface parallax displacement.
+
+    Calculate the displacement due to parallax error.  Input parameters are
+    identical to :func:`get_parallax_corrected_lonlats`.
+
+    Returns:
+        number or array: parallax displacement in meter
+    """
+    (corr_lon, corr_lat) = get_parallax_corrected_lonlats(sat_lon, sat_lat, sat_alt, lon, lat, height)
+    # Get parallax displacement
+    geod = Geod(ellps="sphere")
+    _, _, parallax_dist = geod.inv(corr_lon, corr_lat, lon, lat)
+    return parallax_dist
+
+
 def _get_parallax_shift_xyz(sat_lon, sat_lat, sat_alt, lon, lat, parallax_distance):
     """Calculate the parallax shift in cartesian coordinates.
 
@@ -130,13 +149,11 @@ def _get_parallax_shift_xyz(sat_lon, sat_lat, sat_alt, lon, lat, parallax_distan
             in geodetic coordinates [degrees]
         lat (array or number): Latitudes of pixel/pixels to be corrected, in
             geodetic coordinates [degrees]
-        parallax_distance (number): Cloud to ground distance with parallax
+        parallax_distance (array or number): Cloud to ground distance with parallax
             effect [m].
 
     Returns:
-        Parallax shift in cartesian coordinates.  If the input parameters
-        ``sat_alt`` and ``parallax_distance`` are in meter, then the returned
-        coordinates will also be in meter.
+        Parallax shift in cartesian coordinates in meter.
     """
     sat_xyz = np.hstack(lonlat2xyz(sat_lon, sat_lat)) * sat_alt
     cth_xyz = np.stack(lonlat2xyz(lon, lat), axis=-1) * EARTH_RADIUS*1e3  # km → m
@@ -159,11 +176,11 @@ def _get_satellite_elevation(sat_lon, sat_lat, sat_alt, lon, lat):
     return elevation
 
 
-def _calculate_parallax_distance(height, elevation):
-    """Calculate cloud to ground distance with parallax effect.
+def _calculate_slant_cloud_distance(height, elevation):
+    """Calculate slant cloud to ground distance.
 
     From (cloud top) height and satellite elevation, calculate the
-    cloud-to-ground distance taking the parallax effect into consideration.
+    slant cloud-to-ground distance along the line of sight of the satellite.
     """
     if np.isscalar(elevation) and elevation == 0:
         raise NotImplementedError(
@@ -180,7 +197,7 @@ class ParallaxCorrection:
     """Parallax correction calculations.
 
     This class contains higher-level functionality to wrap the parallax
-    correction calculations in :func:`forward_parallax`.  The class is
+    correction calculations in :func:`get_parallax_corrected_lonlats`.  The class is
     initialised using a base area, which is the area for which a corrected
     geolocation will be calculated.  The resulting object is a callable.
     Calling the object with an array of (cloud top) heights returns a
@@ -197,10 +214,12 @@ class ParallaxCorrection:
     A note on the algorithm and the implementation.  Parallax correction
     is inherently an inverse problem.  The reported geolocation in
     satellite data files is the true location plus the parallax error.
-    Therefore, this class first calculates the parallax error (using
-    :func:`forward_parallax`), which gives a shifted longitude and
+    Therefore, this class first calculates the true geolocation (using
+    :func:`get_parallax_corrected_lonlats`), which gives a shifted longitude and
     shifted latitude on an irregular grid.  The difference between
     the original and the shifted grid is the parallax error or shift.
+    The magnitude of this error can be estimated with
+    :func:`get_surface_parallax_displacement`.
     With this difference, we need to invert the parallax correction to
     calculate the corrected geolocation.  Due to parallax correction,
     high clouds shift a lot, low clouds shift a little, and cloud-free
@@ -306,11 +325,11 @@ def corrected_area(self, cth_dataset,
 
         (base_lon, base_lat) = self.base_area.get_lonlats(chunks=lonlat_chunks)
         # calculate the shift/error due to the parallax effect
-        (shifted_lon, shifted_lat) = forward_parallax(
+        (corrected_lon, corrected_lat) = get_parallax_corrected_lonlats(
                 sat_lon, sat_lat, sat_alt_m,
                 base_lon, base_lat, cth_dataset.data)
 
-        shifted_area = self._get_swathdef_from_lon_lat(shifted_lon, shifted_lat)
+        shifted_area = self._get_swathdef_from_lon_lat(corrected_lon, corrected_lat)
 
         # But we are not actually moving pixels, rather we want a
         # coordinate transformation. With this transformation we approximately
@@ -375,13 +394,14 @@ def _check_overlap(self, cth_dataset):
     def _get_corrected_lon_lat(self, base_lon, base_lat, shifted_area):
         """Calculate the corrected lon/lat based from the shifted area.
 
-        After calculating the shifted area based on :func:`forward_parallax`,
+        After calculating the shifted area based on
+        :func:`get_parallax_corrected_lonlats`,
         we invert the parallax error and estimate where those pixels came from.
         For details on the algorithm, see the class docstring.
         """
-        (shifted_lon, shifted_lat) = shifted_area.get_lonlats(chunks=1024)
-        lon_diff = shifted_lon - base_lon
-        lat_diff = shifted_lat - base_lat
+        (corrected_lon, corrected_lat) = shifted_area.get_lonlats(chunks=1024)
+        lon_diff = corrected_lon - base_lon
+        lat_diff = corrected_lat - base_lat
         # We use the bucket resampler here, because parallax correction
         # inevitably means there will be 2 source pixels ending up in the same
         # destination pixel.  We want to choose the biggest shift (max abs in
@@ -398,15 +418,15 @@ def _get_corrected_lon_lat(self, base_lon, base_lat, shifted_area):
         #   so a cloud that was rectangular at the start may no longer be
         #   rectangular at the end
         bur = BucketResampler(self.base_area,
-                              da.array(shifted_lon), da.array(shifted_lat))
+                              da.array(corrected_lon), da.array(corrected_lat))
         inv_lat_diff = bur.get_abs_max(lat_diff)
         inv_lon_diff = bur.get_abs_max(lon_diff)
 
         inv_lon = base_lon - inv_lon_diff
         inv_lat = base_lat - inv_lat_diff
         if self.debug_mode:
-            self.diagnostics["shifted_lon"] = shifted_lon
-            self.diagnostics["shifted_lat"] = shifted_lat
+            self.diagnostics["corrected_lon"] = corrected_lon
+            self.diagnostics["corrected_lat"] = corrected_lat
             self.diagnostics["inv_lon"] = inv_lon
             self.diagnostics["inv_lat"] = inv_lat
             self.diagnostics["inv_lon_diff"] = inv_lon_diff
@@ -426,7 +446,7 @@ class ParallaxCorrectionModifier(ModifierBase):
 
     Apply parallax correction as a modifier.  Uses the
     :class:`ParallaxCorrection` class, which in turn uses the
-    :func:`forward_parallax` function.  See the documentation there for
+    :func:`get_parallax_corrected_lonlats` function.  See the documentation there for
     details on the behaviour.
 
     To use this, add to ``composites/visir.yaml`` within ``SATPY_CONFIG_PATH``

satpy/tests/modifier_tests/test_parallax.py

@@ -80,42 +80,42 @@ def _get_attrs(lat, lon, height=35_000):
 class TestForwardParallax:
     """Test the forward parallax function with various inputs."""
 
-    def test_forward_parallax_ssp(self):
+    def test_get_parallax_corrected_lonlats_ssp(self):
         """Test that at SSP, parallax correction does nothing."""
-        from ...modifiers.parallax import forward_parallax
+        from ...modifiers.parallax import get_parallax_corrected_lonlats
         sat_lat = sat_lon = lon = lat = 0.
         height = 5000.  # m
         sat_alt = 30_000_000.  # m
-        corr_lon, corr_lat = forward_parallax(
+        corr_lon, corr_lat = get_parallax_corrected_lonlats(
             sat_lon, sat_lat, sat_alt, lon, lat, height)
         assert corr_lon == corr_lat == 0
 
-    def test_forward_parallax_clearsky(self):
+    def test_get_parallax_corrected_lonlats_clearsky(self):
         """Test parallax correction for clearsky case (returns NaN)."""
-        from ...modifiers.parallax import forward_parallax
+        from ...modifiers.parallax import get_parallax_corrected_lonlats
         sat_lat = sat_lon = 0
         lat = np.linspace(-20, 20, 25).reshape(5, 5)
         lon = np.linspace(-20, 20, 25).reshape(5, 5).T
         height = np.full((5, 5), np.nan)  # no CTH --> clearsky
         sat_alt = 35_000_000.  # m above surface
-        (corr_lon, corr_lat) = forward_parallax(
+        (corr_lon, corr_lat) = get_parallax_corrected_lonlats(
             sat_lon, sat_lat, sat_alt, lon, lat, height)
         # clearsky becomes NaN
         assert np.isnan(corr_lon).all()
         assert np.isnan(corr_lat).all()
 
     @pytest.mark.parametrize("lat,lon", [(0, 0), (0, 40), (0, 179.9)])
     @pytest.mark.parametrize("resolution", [0.01, 0.5, 10])
-    def test_forward_parallax_cloudy_ssp(self, lat, lon, resolution):
+    def test_get_parallax_corrected_lonlats_cloudy_ssp(self, lat, lon, resolution):
         """Test parallax correction for fully cloudy scene at SSP."""
-        from ...modifiers.parallax import forward_parallax
+        from ...modifiers.parallax import get_parallax_corrected_lonlats
 
         N = 5
         lats = np.linspace(lat-N*resolution, lat+N*resolution, 25).reshape(N, N)
         lons = np.linspace(lon-N*resolution, lon+N*resolution, 25).reshape(N, N).T
         height = np.full((N, N), 10_000)  # constant high clouds at 10 km
         sat_alt = 35_000_000.  # satellite at 35 Mm
-        (corr_lon, corr_lat) = forward_parallax(
+        (corr_lon, corr_lat) = get_parallax_corrected_lonlats(
             lon, lat, sat_alt, lons, lats, height)
         # confirm movements behave as expected
         geod = Geod(ellps="sphere")
@@ -137,25 +137,25 @@ def test_forward_parallax_cloudy_ssp(self, lat, lon, resolution):
         assert (np.diff(np.diag(corr_delta)[:N//2+1]) < 0).all()
         assert (np.diff(np.diag(corr_delta)[N//2:]) > 0).all()
 
-    def test_forward_parallax_cloudy_slant(self):
+    def test_get_parallax_corrected_lonlats_cloudy_slant(self):
         """Test parallax correction for fully cloudy scene (not SSP)."""
-        from ...modifiers.parallax import forward_parallax
+        from ...modifiers.parallax import get_parallax_corrected_lonlats
         sat_lat = sat_lon = 0
         lat = np.linspace(-20, 20, 25).reshape(5, 5)
         lon = np.linspace(-20, 20, 25).reshape(5, 5).T
         height = np.full((5, 5), 10_000)  # constant high clouds at 10 km
         sat_alt = 35_000_000.  # satellite at 35 Mm
-        (corr_lon, corr_lat) = forward_parallax(
+        (corr_lon, corr_lat) = get_parallax_corrected_lonlats(
             sat_lon, sat_lat, sat_alt, lon, lat, height)
         # reference value from Simon Proud
         np.testing.assert_allclose(
             corr_lat[4, 4], 19.955, rtol=5e-4)
         np.testing.assert_allclose(
             corr_lon[4, 4], 19.960, rtol=5e-4)
 
-    def test_forward_parallax_mixed(self):
+    def test_get_parallax_corrected_lonlats_mixed(self):
         """Test parallax correction for mixed cloudy case."""
-        from ...modifiers.parallax import forward_parallax
+        from ...modifiers.parallax import get_parallax_corrected_lonlats
 
         sat_lon = sat_lat = 0
         sat_alt = 35_785_831.0  # m
@@ -167,7 +167,7 @@ def test_forward_parallax_mixed(self):
             [np.nan, 7000., 8000., 9000., np.nan],
             [np.nan, 7000., 7000., 7000., np.nan],
             [np.nan, 4000., 3000., np.nan, np.nan]])
-        (corrected_lon, corrected_lat) = forward_parallax(
+        (corrected_lon, corrected_lat) = get_parallax_corrected_lonlats(
             sat_lon, sat_lat, sat_alt, lon, lat, alt)
         assert corrected_lon.shape == lon.shape
         assert corrected_lat.shape == lat.shape
@@ -178,19 +178,27 @@ def test_forward_parallax_mixed(self):
         assert np.isfinite(corrected_lon[~np.isnan(alt)]).all()
         assert np.isfinite(corrected_lat[~np.isnan(alt)]).all()
 
-    def test_forward_parallax_horizon(self):
+    def test_get_parallax_corrected_lonlats_horizon(self):
         """Test that exception is raised if satellites exactly at the horizon.
 
         Test the rather unlikely case of a satellite elevation of exactly 0
         """
-        from ...modifiers.parallax import forward_parallax
+        from ...modifiers.parallax import get_parallax_corrected_lonlats
         sat_lat = sat_lon = lon = lat = 0.
         height = 5000.
         sat_alt = 30_000_000.
         with unittest.mock.patch("satpy.modifiers.parallax.get_observer_look") as smpg:
             smpg.return_value = (0, 0)
             with pytest.raises(NotImplementedError):
-                forward_parallax(sat_lon, sat_lat, sat_alt, lon, lat, height)
+                get_parallax_corrected_lonlats(sat_lon, sat_lat, sat_alt, lon, lat, height)
+
+    def test_get_surface_parallax_displacement(self):
+        """Test surface parallax displacement."""
+        from ...modifiers.parallax import get_surface_parallax_displacement
+
+        val = get_surface_parallax_displacement(
+                0, 0, 36_000_000, 0, 10, 10_000)
+        np.testing.assert_allclose(val, 2141.2404451757875)
 
 
 class TestParallaxCorrectionClass: