From c865cfb91bef609199bb6ae60658b8f61ca4e909 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Mon, 2 Jun 2025 15:49:42 +0200 Subject: [PATCH 1/4] Add a threshold for the minimization to fail --- pyorbital/geoloc_avhrr.py | 11 ++++++++--- pyorbital/tests/test_geoloc.py | 21 +++++++++++++++++++++ 2 files changed, 29 insertions(+), 3 deletions(-) diff --git a/pyorbital/geoloc_avhrr.py b/pyorbital/geoloc_avhrr.py index f4f8d040..a9e36d7f 100644 --- a/pyorbital/geoloc_avhrr.py +++ b/pyorbital/geoloc_avhrr.py @@ -62,7 +62,7 @@ def compute_avhrr_gcps_lonlatalt(gcps, max_scan_angle, rpy, start_time, tle) -> return get_lonlatalt(pixels_pos, s_times) -def estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, start_time, tle, max_scan_angle): +def estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, start_time, tle, max_scan_angle, max_accepted_distance=5000): """Estimate time offset and attitude deviations from gcps. Provided reference longitudes and latitudes for the gcps, this function minimises the attitude and time offset @@ -72,7 +72,8 @@ def estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, start_time, distances = compute_gcp_distances_to_reference_lonlats((0, 0, 0, 0), gcps, start_time, tle, max_scan_angle, (ref_lons, ref_lats)) - logger.debug(f"GCP distances: median {np.median(distances)}, std {np.std(distances)}") + original_median_distance = np.median(distances) + logger.debug(f"GCP distances: median {original_median_distance}, std {np.std(distances)}") # we need to work in seconds*1e3 to avoid the nanosecond precision issue res = minimize(compute_gcp_accumulated_squared_distances_to_reference_lonlats, x0=(0, 0, 0, 0), @@ -84,7 +85,11 @@ def estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, start_time, logger.debug(f"Estimated time difference to {time_diff} seconds, attitude to {roll}, {pitch}, {yaw} degrees") distances = compute_gcp_distances_to_reference_lonlats(res.x, gcps, start_time, tle, max_scan_angle, (ref_lons, ref_lats)) - logger.debug(f"Remaining GCP distances: median {np.median(distances)}, std {np.std(distances)}") + + minimized_median_distance = np.median(distances) + logger.debug(f"Remaining GCP distances: median {minimized_median_distance}, std {np.std(distances)}") + if minimized_median_distance > max_accepted_distance: + raise RuntimeError("Minimization did not provide significant improvement.") return time_diff, roll, pitch, yaw diff --git a/pyorbital/tests/test_geoloc.py b/pyorbital/tests/test_geoloc.py index c4d59e49..bad25b84 100644 --- a/pyorbital/tests/test_geoloc.py +++ b/pyorbital/tests/test_geoloc.py @@ -226,6 +226,27 @@ def test_minimize_geoloc_error(): assert yaw == pytest.approx(ref_yaw, abs=1e-2) +def test_minimize_geoloc_error_fails(): + """Test minimizing the distance to a set of gcps.""" + # Couple of example Two Line Elements + tle1 = "1 33591U 09005A 12345.45213434 .00000391 00000-0 24004-3 0 6113" + tle2 = "2 33591 098.8821 283.2036 0013384 242.4835 117.4960 14.11432063197875" + tle = (tle1, tle2) + + # Choosing a specific time, this should be relatively close to the issue date of the TLE + t = dt.datetime(2012, 12, 12, 4, 16, 1, 575000) + + ref_time_displacement = 60 + ref_time = t + dt.timedelta(seconds=ref_time_displacement) + rpy = (0, 0, 0) + max_scan_angle = 55.37 + # gcps are line/col + gcps = np.array([[2, 500], [1500, 700], [20, 1000], [500, 1100], [100, 2000]]) + ref_lons, ref_lats, _ = compute_avhrr_gcps_lonlatalt(gcps, max_scan_angle, rpy, ref_time, tle) + with pytest.raises(RuntimeError): + _ = estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, t, + tle, max_scan_angle, 5000) + class TestGeolocDefs: """Test the instrument definitions.""" From a6cf07f12e7217c68c2d245240218f440392924c Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Mon, 2 Jun 2025 16:04:36 +0200 Subject: [PATCH 2/4] Return gcp distances before and after minimization --- pyorbital/geoloc_avhrr.py | 14 ++++++-------- pyorbital/tests/test_geoloc.py | 26 +++----------------------- 2 files changed, 9 insertions(+), 31 deletions(-) diff --git a/pyorbital/geoloc_avhrr.py b/pyorbital/geoloc_avhrr.py index a9e36d7f..cdbb6051 100644 --- a/pyorbital/geoloc_avhrr.py +++ b/pyorbital/geoloc_avhrr.py @@ -62,7 +62,7 @@ def compute_avhrr_gcps_lonlatalt(gcps, max_scan_angle, rpy, start_time, tle) -> return get_lonlatalt(pixels_pos, s_times) -def estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, start_time, tle, max_scan_angle, max_accepted_distance=5000): +def estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, start_time, tle, max_scan_angle): """Estimate time offset and attitude deviations from gcps. Provided reference longitudes and latitudes for the gcps, this function minimises the attitude and time offset @@ -70,10 +70,10 @@ def estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, start_time, """ from scipy.optimize import minimize - distances = compute_gcp_distances_to_reference_lonlats((0, 0, 0, 0), gcps, start_time, tle, max_scan_angle, - (ref_lons, ref_lats)) - original_median_distance = np.median(distances) - logger.debug(f"GCP distances: median {original_median_distance}, std {np.std(distances)}") + original_distances = compute_gcp_distances_to_reference_lonlats((0, 0, 0, 0), gcps, start_time, tle, max_scan_angle, + (ref_lons, ref_lats)) + original_median_distance = np.median(original_distances) + logger.debug(f"GCP distances: median {original_median_distance}, std {np.std(original_distances)}") # we need to work in seconds*1e3 to avoid the nanosecond precision issue res = minimize(compute_gcp_accumulated_squared_distances_to_reference_lonlats, x0=(0, 0, 0, 0), @@ -88,10 +88,8 @@ def estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, start_time, minimized_median_distance = np.median(distances) logger.debug(f"Remaining GCP distances: median {minimized_median_distance}, std {np.std(distances)}") - if minimized_median_distance > max_accepted_distance: - raise RuntimeError("Minimization did not provide significant improvement.") - return time_diff, roll, pitch, yaw + return time_diff, (roll, pitch, yaw), (original_distances, distances) def compute_gcp_accumulated_squared_distances_to_reference_lonlats( diff --git a/pyorbital/tests/test_geoloc.py b/pyorbital/tests/test_geoloc.py index bad25b84..55930cfe 100644 --- a/pyorbital/tests/test_geoloc.py +++ b/pyorbital/tests/test_geoloc.py @@ -220,33 +220,13 @@ def test_minimize_geoloc_error(): # gcps are line/col gcps = np.array([[2, 500], [1500, 700], [20, 1000], [500, 1100], [100, 2000]]) ref_lons, ref_lats, _ = compute_avhrr_gcps_lonlatalt(gcps, max_scan_angle, rpy, ref_time, tle) - time_diff, roll, pitch, yaw = estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, t, - tle, max_scan_angle) + time_diff, (roll, pitch, yaw), (do, dm) = estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, t, + tle, max_scan_angle) assert time_diff == pytest.approx(ref_time_displacement, abs=1e-2) assert yaw == pytest.approx(ref_yaw, abs=1e-2) + assert min(do) > max(dm) -def test_minimize_geoloc_error_fails(): - """Test minimizing the distance to a set of gcps.""" - # Couple of example Two Line Elements - tle1 = "1 33591U 09005A 12345.45213434 .00000391 00000-0 24004-3 0 6113" - tle2 = "2 33591 098.8821 283.2036 0013384 242.4835 117.4960 14.11432063197875" - tle = (tle1, tle2) - - # Choosing a specific time, this should be relatively close to the issue date of the TLE - t = dt.datetime(2012, 12, 12, 4, 16, 1, 575000) - - ref_time_displacement = 60 - ref_time = t + dt.timedelta(seconds=ref_time_displacement) - rpy = (0, 0, 0) - max_scan_angle = 55.37 - # gcps are line/col - gcps = np.array([[2, 500], [1500, 700], [20, 1000], [500, 1100], [100, 2000]]) - ref_lons, ref_lats, _ = compute_avhrr_gcps_lonlatalt(gcps, max_scan_angle, rpy, ref_time, tle) - with pytest.raises(RuntimeError): - _ = estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, t, - tle, max_scan_angle, 5000) - class TestGeolocDefs: """Test the instrument definitions.""" From fe59c0b969b0e715b24a836961024d8626757559 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Thu, 5 Jun 2025 15:11:35 +0200 Subject: [PATCH 3/4] Estimate time offset from gcps --- pyorbital/geoloc_avhrr.py | 38 +++++++++++++++++++++++++++++++++- pyorbital/tests/test_geoloc.py | 29 +++++++++++++++++++++++++- 2 files changed, 65 insertions(+), 2 deletions(-) diff --git a/pyorbital/geoloc_avhrr.py b/pyorbital/geoloc_avhrr.py index cdbb6051..f54b47f3 100644 --- a/pyorbital/geoloc_avhrr.py +++ b/pyorbital/geoloc_avhrr.py @@ -82,7 +82,7 @@ def estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, start_time, if not res.success: raise RuntimeError("Time and attitude estimation did not converge") time_diff, roll, pitch, yaw = res.x * [1e3, 1, 1, 1] - logger.debug(f"Estimated time difference to {time_diff} seconds, attitude to {roll}, {pitch}, {yaw} degrees") + logger.debug(f"Estimated time difference to {time_diff} seconds, attitude to {np.rad2deg(roll)}, {np.rad2deg(pitch)}, {np.rad2deg(yaw)} degrees") distances = compute_gcp_distances_to_reference_lonlats(res.x, gcps, start_time, tle, max_scan_angle, (ref_lons, ref_lats)) @@ -92,6 +92,42 @@ def estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, start_time, return time_diff, (roll, pitch, yaw), (original_distances, distances) +def estimate_time_offset(gcps, ref_lons, ref_lats, start_time, tle, max_scan_angle): + """Estimate time offset from gcps. + + Provided reference longitudes and latitudes for the gcps, this function minimises the time offset + needed to match the gcp coordinates to the reference coordinates. + """ + from scipy.optimize import minimize + + original_distances = compute_gcp_distances_to_reference_lonlats((0, 0, 0, 0), gcps, start_time, tle, max_scan_angle, + (ref_lons, ref_lats)) + original_median_distance = np.median(original_distances) + logger.debug(f"GCP distances: median {original_median_distance}, std {np.std(original_distances)}") + + def gcp_distance_for_time(time): + dist = compute_gcp_accumulated_squared_distances_to_reference_lonlats((time[0], 0, 0, 0), gcps, start_time, tle, max_scan_angle, (ref_lons, ref_lats)) + return dist + + # we need to work in seconds*1e3 to avoid the nanosecond precision issue + res = minimize(gcp_distance_for_time, + x0=(0,), + bounds=((-0.03, 0.03),), + options=dict(ftol=1e-1), + ) + if not res.success: + raise RuntimeError("Time offset estimation did not converge") + time_diff, = res.x * [1e3,] + logger.debug(f"Estimated time difference to {time_diff} seconds") + distances = compute_gcp_distances_to_reference_lonlats((res.x[0], 0, 0, 0), gcps, start_time, tle, max_scan_angle, + (ref_lons, ref_lats)) + + minimized_median_distance = np.median(distances) + logger.debug(f"Remaining GCP distances: median {minimized_median_distance}, std {np.std(distances)}") + + return time_diff, (original_distances, distances) + + def compute_gcp_accumulated_squared_distances_to_reference_lonlats( variables, gcps, start_time, tle, max_scan_angle, refs): """Compute the summed squared distance fot gcps to reference lonlats. diff --git a/pyorbital/tests/test_geoloc.py b/pyorbital/tests/test_geoloc.py index 55930cfe..63ceb9ae 100644 --- a/pyorbital/tests/test_geoloc.py +++ b/pyorbital/tests/test_geoloc.py @@ -29,7 +29,11 @@ import pytest from pyorbital.geoloc import ScanGeometry, geodetic_lat, qrotate, subpoint -from pyorbital.geoloc_avhrr import compute_avhrr_gcps_lonlatalt, estimate_time_and_attitude_deviations +from pyorbital.geoloc_avhrr import ( + compute_avhrr_gcps_lonlatalt, + estimate_time_and_attitude_deviations, + estimate_time_offset, +) from pyorbital.geoloc_instrument_definitions import ( amsua, ascat, @@ -227,6 +231,29 @@ def test_minimize_geoloc_error(): assert min(do) > max(dm) +def test_minimize_time_error(): + """Test minimizing the distance to a set of gcps using only time offset.""" + # Couple of example Two Line Elements + tle1 = "1 33591U 09005A 12345.45213434 .00000391 00000-0 24004-3 0 6113" + tle2 = "2 33591 098.8821 283.2036 0013384 242.4835 117.4960 14.11432063197875" + tle = (tle1, tle2) + + # Choosing a specific time, this should be relatively close to the issue date of the TLE + t = dt.datetime(2012, 12, 12, 4, 16, 1, 575000) + + ref_time_displacement = 20 + ref_time = t + dt.timedelta(seconds=ref_time_displacement) + rpy = (0, 0, 0) + max_scan_angle = 55.37 + # gcps are line/col + gcps = np.array([[2, 500], [1500, 700], [20, 1000], [500, 1100], [100, 2000]]) + ref_lons, ref_lats, _ = compute_avhrr_gcps_lonlatalt(gcps, max_scan_angle, rpy, ref_time, tle) + time_diff, (do, dm) = estimate_time_offset(gcps, ref_lons, ref_lats, t, + tle, max_scan_angle) + assert time_diff == pytest.approx(ref_time_displacement, abs=1e-2) + assert min(do) > max(dm) + + class TestGeolocDefs: """Test the instrument definitions.""" From 0a7ede2d838460b94c8049eab9471bc2208e5080 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Fri, 6 Jun 2025 10:46:39 +0200 Subject: [PATCH 4/4] Fix linting --- pyorbital/geoloc_avhrr.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pyorbital/geoloc_avhrr.py b/pyorbital/geoloc_avhrr.py index f54b47f3..43bac2ce 100644 --- a/pyorbital/geoloc_avhrr.py +++ b/pyorbital/geoloc_avhrr.py @@ -82,7 +82,8 @@ def estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, start_time, if not res.success: raise RuntimeError("Time and attitude estimation did not converge") time_diff, roll, pitch, yaw = res.x * [1e3, 1, 1, 1] - logger.debug(f"Estimated time difference to {time_diff} seconds, attitude to {np.rad2deg(roll)}, {np.rad2deg(pitch)}, {np.rad2deg(yaw)} degrees") + logger.debug(f"Estimated time difference to {time_diff} seconds, " + f"attitude to {np.rad2deg(roll)}, {np.rad2deg(pitch)}, {np.rad2deg(yaw)} degrees") distances = compute_gcp_distances_to_reference_lonlats(res.x, gcps, start_time, tle, max_scan_angle, (ref_lons, ref_lats)) @@ -106,7 +107,8 @@ def estimate_time_offset(gcps, ref_lons, ref_lats, start_time, tle, max_scan_ang logger.debug(f"GCP distances: median {original_median_distance}, std {np.std(original_distances)}") def gcp_distance_for_time(time): - dist = compute_gcp_accumulated_squared_distances_to_reference_lonlats((time[0], 0, 0, 0), gcps, start_time, tle, max_scan_angle, (ref_lons, ref_lats)) + dist = compute_gcp_accumulated_squared_distances_to_reference_lonlats((time[0], 0, 0, 0), gcps, start_time, tle, + max_scan_angle, (ref_lons, ref_lats)) return dist # we need to work in seconds*1e3 to avoid the nanosecond precision issue