From 5130cdbc0d7cdcfb545692b8cdaa4b2d10dcee71 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Thu, 28 Jul 2022 11:19:59 +0200 Subject: [PATCH] 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 | 3 + satpy/modifiers/parallax.py | 84 +++++++++++++-------- satpy/tests/modifier_tests/test_parallax.py | 44 ++++++----- 3 files changed, 81 insertions(+), 50 deletions(-) diff --git a/doc/source/modifiers.rst b/doc/source/modifiers.rst index 9682d4f851..23bb73d9d9 100644 --- a/doc/source/modifiers.rst +++ b/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 _`. + .. versionadded:: 0.37 diff --git a/satpy/modifiers/parallax.py b/satpy/modifiers/parallax.py index 07b7b61701..0464f1a5f6 100644 --- a/satpy/modifiers/parallax.py +++ b/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,12 +105,12 @@ 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) @@ -116,6 +118,23 @@ def forward_parallax(sat_lon, sat_lat, sat_alt, lon, lat, height): 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`` diff --git a/satpy/tests/modifier_tests/test_parallax.py b/satpy/tests/modifier_tests/test_parallax.py index fda0e1d718..70d51b49e9 100644 --- a/satpy/tests/modifier_tests/test_parallax.py +++ b/satpy/tests/modifier_tests/test_parallax.py @@ -80,25 +80,25 @@ 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() @@ -106,16 +106,16 @@ def test_forward_parallax_clearsky(self): @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,15 +137,15 @@ 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( @@ -153,9 +153,9 @@ def test_forward_parallax_cloudy_slant(self): 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: