Return distances from GCP minimization

pyorbital/geoloc_avhrr.py

@@ -70,9 +70,10 @@
     """
     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))
-    logger.debug(f"GCP distances: median {np.median(distances)}, 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),
@@ -81,12 +82,52 @@
     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, "
+                 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))
-    logger.debug(f"Remaining GCP distances: median {np.median(distances)}, std {np.std(distances)}")
 
-    return time_diff, roll, pitch, yaw
+    minimized_median_distance = np.median(distances)
+    logger.debug(f"Remaining GCP distances: median {minimized_median_distance}, std {np.std(distances)}")
+
+    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(

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,
@@ -220,10 +224,34 @@ 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_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: