Skip to content

Fix the parallax correction code for the unit and computation of satellite altitude #3217

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 9 commits into
base: main
Choose a base branch
from
Open

Fix the parallax correction code for the unit and computation of satellite altitude #3217

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
3429bc6
add Earth radius to satellite altitude to compute the satellite carte…
ameraner Aug 25, 2025
da37b4b
fix docstring of get_satpos to indicate that the altitude is in meter…
ameraner Aug 25, 2025
cbab880
fix the automated _get_satpos_from_platform_name to return the height…
ameraner Aug 25, 2025
5a64be6
remove faulty meters convertion in _get_satpos_from_cth
ameraner Aug 25, 2025
7d5d3d2
fix test_get_surface_parallax_displacement value - displacement gets …
ameraner Aug 25, 2025
34b9ec0
fix test_correct_area_partlycloudy results
ameraner Aug 25, 2025
6e3aa52
fix unit of altitude in tests and fix failing tests
ameraner Aug 25, 2025
2d31272
fix test for tle-retrieved altitude
ameraner Aug 26, 2025
abccb51
Merge branch 'main' into fix_parallax
ameraner Aug 26, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions satpy/modifiers/parallax.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -155,7 +155,7 @@ def _get_parallax_shift_xyz(sat_lon, sat_lat, sat_alt, lon, lat, parallax_distan
Returns:
Parallax shift in cartesian coordinates in meter.
"""
sat_xyz = np.hstack(lonlat2xyz(sat_lon, sat_lat)) * sat_alt
sat_xyz = np.hstack(lonlat2xyz(sat_lon, sat_lat)) * (sat_alt + EARTH_RADIUS * 1e3)
cth_xyz = np.stack(lonlat2xyz(lon, lat), axis=-1) * EARTH_RADIUS*1e3 # km → m
delta_xyz = cth_xyz - sat_xyz
sat_distance = np.sqrt((delta_xyz*delta_xyz).sum(axis=-1))
Expand Down Expand Up @@ -542,12 +542,12 @@ def _get_corrector(self, base_area):


def _get_satpos_from_cth(cth_dataset):
"""Obtain satellite position from CTH dataset, height in meter.
"""Obtain satellite position from CTH dataset.

From a CTH dataset, obtain the satellite position lon, lat, altitude/m,
From a CTH dataset, obtain the satellite position lon, lat, altitude [m],
either directly from orbital parameters, or, when missing, from the
platform name using pyorbital and skyfield.
"""
(sat_lon, sat_lat, sat_alt_km) = get_satpos(
(sat_lon, sat_lat, sat_alt) = get_satpos(
cth_dataset, use_tle=True)
return (sat_lon, sat_lat, sat_alt_km * 1000)
return (sat_lon, sat_lat, sat_alt)
53 changes: 26 additions & 27 deletions satpy/tests/modifier_tests/test_parallax.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,11 +73,11 @@ def _get_fake_areas(center, sizes, resolution, code=4326): # noqa: D417
for size in sizes]


def _get_attrs(lat, lon, height=35_000):
def _get_attrs(lat, lon, height=35_000_000):
"""Get attributes for datasets in fake scene."""
return {
"orbital_parameters": {
"satellite_actual_altitude": height, # in km above surface
"satellite_actual_altitude": height, # in m above surface
"satellite_actual_longitude": lon,
"satellite_actual_latitude": lat},
"units": "m" # does not apply to orbital parameters, I think!
Expand Down Expand Up @@ -205,7 +205,7 @@ def test_get_surface_parallax_displacement(self):

val = get_surface_parallax_displacement(
0, 0, 36_000_000, 0, 10, 10_000)
np.testing.assert_allclose(val, 2141.2404451757875)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice, this tells us how far we were off :) 👍

np.testing.assert_allclose(val, 2075.6863144582703)


@pytest.mark.filterwarnings("ignore:Overlap checking not implemented:UserWarning")
Expand Down Expand Up @@ -241,7 +241,7 @@ def test_correct_area_clearsky(self, sat_pos, ar_pos, resolution, caplog):
{"CTH_clear": np.full((large, large), np.nan)},
daskify=False,
area=fake_area_large,
common_attrs=_get_attrs(sat_lat, sat_lon, 35_000))
common_attrs=_get_attrs(sat_lat, sat_lon, 35_000_000))

with caplog.at_level(logging.DEBUG):
new_area = corrector(sc["CTH_clear"])
Expand Down Expand Up @@ -274,7 +274,7 @@ def test_correct_area_ssp(self, lat, lon, resolution):
{"CTH_constant": np.full((large, large), 10000)},
daskify=False,
area=fake_area_large,
common_attrs=_get_attrs(lat, lon, 35_000))
common_attrs=_get_attrs(lat, lon, 35_000_000))
new_area = corrector(sc["CTH_constant"])
assert new_area.shape == fake_area_small.shape
old_lonlats = fake_area_small.get_lonlats()
Expand Down Expand Up @@ -323,29 +323,27 @@ def test_correct_area_partlycloudy(self, daskify):
])},
daskify=daskify,
area=fake_area_large,
common_attrs=_get_attrs(0, 0, 40_000))
common_attrs=_get_attrs(0, 0, 40_000_000))
new_area = corrector(sc["CTH"])
assert new_area.shape == fake_area_small.shape
(new_lons, new_lats) = new_area.get_lonlats()
assert fake_area_lons[3, 4] != new_lons[3, 4]

np.testing.assert_allclose(
new_lons,
np.array([
[np.nan, np.nan, 0.0, 0.1, 0.2],
[-0.20078652, -0.10044222, 0.0, 0.1, 0.2],
[-0.20068529, -0.10034264, 0.0, 0.1, 0.2],
[np.nan, np.nan, np.nan, np.nan, np.nan],
[-0.20048537, -0.10038778, 0., 0.10038778, 0.20058219]]),
np.array([[np.nan, np.nan, 0., 0.1, 0.2],
[-0.20077639, -0.10043652, 0., 0.1, 0.2],
[-0.20067643, -0.10033821, 0., 0.1, 0.2],
[np.nan, np.nan, np.nan, np.nan, np.nan],
[-0.20047906, -0.10038274, 0., 0.10038274, 0.20057462]]),
rtol=1e-5)
np.testing.assert_allclose(
new_lats,
np.array([
[np.nan, np.nan, 50.2, 50.2, 50.2],
[50.2110675, 50.22493181, 50.1, 50.1, 50.1],
[50.09680357, 50.09680346, 50.0, 50.0, 50.0],
[np.nan, np.nan, np.nan, np.nan, np.nan],
[49.86860622, 49.9097198, 49.90971976, 49.9097198, 49.88231496]]),
np.array([[np.nan, np.nan, 50.2, 50.2, 50.2],
[50.20963341, 50.22331809, 50.1, 50.1, 50.1],
[50.09555017, 50.09555005, 50., 50., 50.],
[np.nan, np.nan, np.nan, np.nan, np.nan],
[49.86771297, 49.90828968, 49.90828963, 49.90828968, 49.88124283]]),
rtol=1e-6)

@pytest.mark.parametrize(("res1", "res2"), [(0.08, 0.3), (0.3, 0.08)])
Expand Down Expand Up @@ -376,7 +374,7 @@ def test_correct_area_clearsky_different_resolutions(self, res1, res2):
{"CTH_clear": np.full(area1.shape, np.nan)},
daskify=False,
area=area1,
common_attrs=_get_attrs(0, 0, 35_000))
common_attrs=_get_attrs(0, 0, 35_000_000))

corrector = ParallaxCorrection(area2)
new_area = corrector(sc["CTH_clear"])
Expand All @@ -398,7 +396,7 @@ def test_correct_area_cloudy_no_overlap(self, ):
{"CTH_constant": np.full((9, 9), 10000)},
daskify=False,
area=fake_area_large,
common_attrs=_get_attrs(0, 0, 35_000))
common_attrs=_get_attrs(0, 0, 35_000_000))

corrector = ParallaxCorrection(fake_area_small)
with pytest.raises(MissingHeightError):
Expand All @@ -418,7 +416,7 @@ def test_correct_area_cloudy_partly_shifted(self, ):
{"CTH_constant": np.full((9, 9), 10000)},
daskify=False,
area=fake_area_large,
common_attrs=_get_attrs(0, 0, 35_000))
common_attrs=_get_attrs(0, 0, 35_000_000))

corrector = ParallaxCorrection(fake_area_small)

Expand All @@ -436,7 +434,7 @@ def test_correct_area_cloudy_same_area(self, ):
{"CTH_constant": np.full((9, 9), 10000)},
daskify=False,
area=area,
common_attrs=_get_attrs(0, 0, 35_000))
common_attrs=_get_attrs(0, 0, 35_000_000))

corrector = ParallaxCorrection(area)
corrector(sc["CTH_constant"])
Expand Down Expand Up @@ -487,11 +485,11 @@ def test_parallax_modifier_interface(self):
fake_bt = xr.DataArray(
np.linspace(220, 230, 25).reshape(5, 5),
dims=("y", "x"),
attrs={"area": area_small, **_get_attrs(0, 0, 35_000)})
attrs={"area": area_small, **_get_attrs(0, 0, 35_000_000)})
cth_clear = xr.DataArray(
np.full((9, 9), np.nan),
dims=("y", "x"),
attrs={"area": area_large, **_get_attrs(0, 0, 35_000)})
attrs={"area": area_large, **_get_attrs(0, 0, 35_000_000)})
modif = ParallaxCorrectionModifier(
name="parallax_corrected_dataset",
prerequisites=[fake_bt, cth_clear],
Expand Down Expand Up @@ -661,7 +659,7 @@ def test_modifier_interface_cloud_moves_to_observer(self, cth, use_dask, test_ar
cloud_location = {
"foroyar": {
7500: (197, 202, 152, 172),
15000: (239, 244, 165, 184)},
15000: (238, 243, 164, 184)},
"ouagadougou": {
7500: (159, 164, 140, 160),
15000: (163, 168, 141, 161)}}
Expand Down Expand Up @@ -691,8 +689,9 @@ def test_modifier_interface_cloud_moves_to_observer(self, cth, use_dask, test_ar
# cloud may shrink.
assert ((res.attrs["area"].get_lonlats()[1][dest_mask]).mean() <
fake_bt.attrs["area"].get_lonlats()[1][cma].mean())
# verify that all pixels at the new cloud location are indeed cloudy
assert (res.data[dest_mask] < 250).all()
# verify that most pixels at the new cloud location are indeed cloudy (calculate mean since some warm
# surface pixels nearby a non-square cloud could be included in the square dest_mask
assert res.data[dest_mask].mean() < 210


_test_yaml_code = """
Expand Down
2 changes: 1 addition & 1 deletion satpy/tests/test_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -224,7 +224,7 @@ def test_get_satpos_from_satname(self, caplog):
assert "Orbital parameters missing from metadata" in caplog.text
np.testing.assert_allclose(
(lon, lat, alt),
(119.39533705010592, -1.1491628298731498, 35803.19986408156),
(119.39533705010592, -1.1491628298731498, 35803199.86408156),
rtol=1e-4,
)

Expand Down
4 changes: 2 additions & 2 deletions satpy/utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -346,7 +346,7 @@ def get_satpos(
the dataset metadata contain the satellite position directly.

Returns:
Geodetic longitude, latitude, altitude [km]
Geodetic longitude, latitude, altitude [m]

"""
if preference is not None and preference not in ("nadir", "actual", "nominal", "projection"):
Expand Down Expand Up @@ -429,7 +429,7 @@ def _get_satpos_from_platform_name(cth_dataset):
gc = es.at(ts.from_datetime(
cth_dataset.attrs["start_time"].replace(tzinfo=datetime.timezone.utc)))
(lat, lon) = wgs84.latlon_of(gc)
height = wgs84.height_of(gc).to("km")
height = wgs84.height_of(gc).to("m")
return (lon.degrees, lat.degrees, height.value)


Expand Down
Loading