ENH: include exclusive values from IERS_B in IERS_Auto (⏰ wait for #16187)

astropy/utils/iers/iers.py

@@ -34,6 +34,7 @@
 from astropy import units as u
 from astropy import utils
 from astropy.table import MaskedColumn, QTable
+from astropy.table import join as table_join
 from astropy.time import Time, TimeDelta
 from astropy.utils.data import (
     clear_download_cache,
@@ -589,9 +590,25 @@
         # Since only 'P' and 'I' are possible and 'P' is guaranteed to come
         # after 'I', we can use searchsorted for 100 times speed up over
         # finding the first index where the flag equals 'P'.
+
+        def select_offset(array):
+            if not hasattr(array, "mask"):
+                return slice(None), 0
+
+            s = ~array.mask
+            # Get the position of first non-masked value.
+            # This is a bit of a hack since s isn't necessarily sorted, but we
+            # get the correct result by letting np.searchsorted *assuming* it is.
+            offset = np.searchsorted(s, True)
+
+            return s, offset
+
+        select1, offset1 = select_offset(table["UT1Flag_A"])
+        select2, offset2 = select_offset(table["PolPMFlag_A"])
+
         p_index = min(
-            np.searchsorted(table["UT1Flag_A"], "P"),
-            np.searchsorted(table["PolPMFlag_A"], "P"),
+            offset1 + np.searchsorted(table["UT1Flag_A"][select1], "P"),
+            offset2 + np.searchsorted(table["PolPMFlag_A"][select2], "P"),
         )
         table.meta["predictive_index"] = p_index
         table.meta["predictive_mjd"] = table["MJD"][p_index].value
@@ -932,30 +949,37 @@
         IERS-A has IERS-B values included, but for reasons unknown these
         do not match the latest IERS-B values (see comments in #4436).
         Here, we use the bundled astropy IERS-B table to overwrite the values
-        in the downloaded IERS-A table.
+        in the downloaded IERS-A table. Values from IERS-B that are not present
+        in IERS-A are inserted.
         """
         iers_b = IERS_B.open()
         # Substitute IERS-B values for existing B values in IERS-A table
         mjd_b = table["MJD"][np.isfinite(table["UT1_UTC_B"])]
         i0 = np.searchsorted(iers_b["MJD"], mjd_b[0], side="left")
         i1 = np.searchsorted(iers_b["MJD"], mjd_b[-1], side="right")
-        iers_b = iers_b[i0:i1]
-        n_iers_b = len(iers_b)
+        iers_b_overlap = iers_b[i0:i1]
+        n_iers_b = len(iers_b_overlap)
         # If there is overlap then replace IERS-A values from available IERS-B
         if n_iers_b > 0:
             # Sanity check that we are overwriting the correct values
-            if not u.allclose(table["MJD"][:n_iers_b], iers_b["MJD"]):
+            if not u.allclose(table["MJD"][:n_iers_b], iers_b_overlap["MJD"]):
                 raise ValueError(
                     "unexpected mismatch when copying IERS-B values into IERS-A table."
                 )
             # Finally do the overwrite
-            table["UT1_UTC_B"][:n_iers_b] = iers_b["UT1_UTC"]
-            table["PM_X_B"][:n_iers_b] = iers_b["PM_x"]
-            table["PM_Y_B"][:n_iers_b] = iers_b["PM_y"]
-            table["dX_2000A_B"][:n_iers_b] = iers_b["dX_2000A"]
-            table["dY_2000A_B"][:n_iers_b] = iers_b["dY_2000A"]
-
-        return table
+            table["UT1_UTC_B"][:n_iers_b] = iers_b_overlap["UT1_UTC"]
+            table["PM_X_B"][:n_iers_b] = iers_b_overlap["PM_x"]
+            table["PM_Y_B"][:n_iers_b] = iers_b_overlap["PM_y"]
+            table["dX_2000A_B"][:n_iers_b] = iers_b_overlap["dX_2000A"]
+            table["dY_2000A_B"][:n_iers_b] = iers_b_overlap["dY_2000A"]
+
+        return table_join(
+            table,
+            iers_b[:i0],
+            join_type="outer",
+            keys="MJD",
+            metadata_conflicts="silent",  # tmp filtering
+        )
 
 
 class earth_orientation_table(ScienceState):

astropy/utils/iers/tests/test_iers.py

@@ -13,6 +13,7 @@
     IERS_B_README,
     IERS_LEAP_SECOND_FILE,
 )
+from numpy.testing import assert_array_equal
 
 from astropy import units as u
 from astropy.config import set_temp_cache
@@ -319,7 +320,7 @@ def test_no_auto_download(self):
     def test_simple(self):
         with iers.conf.set_temp("iers_auto_url", self.iers_a_url_1):
             dat = iers.IERS_Auto.open()
-            assert dat["MJD"][0] == 57359.0 * u.d
+            assert dat["MJD"][0] == 37665.0 * u.d
             assert dat["MJD"][-1] == 57539.0 * u.d
 
             # Pretend we are accessing at a time 7 days after start of predictive data
@@ -386,22 +387,25 @@ def test_IERS_B_parameters_loading_into_IERS_Auto():
     B = iers.IERS_B.open()
 
     ok_A = A["MJD"] <= B["MJD"][-1]
+    assert np.sum(ok_A) == len(B)
     assert not np.all(ok_A), "IERS B covers all of IERS A: should not happen"
 
-    # We only overwrite IERS_B values in the IERS_A table that were already
-    # there in the first place.  Better take that into account.
-    ok_A &= np.isfinite(A["UT1_UTC_B"])
+    # see https://github.com/astropy/astropy/issues/13494
+    assert A["MJD"][0] == A["MJD"].min(), "IERS_Auto isn't properly sorted"
+    assert A["MJD"][-1] == A["MJD"].max(), "IERS_Auto isn't properly sorted"
+    assert A["MJD"][0] <= B["MJD"][0], "IERS_Auto doesn't go back as far as IERS B"
 
     i_B = np.searchsorted(B["MJD"], A["MJD"][ok_A])
+    assert_array_equal(i_B, np.arange(len(B)))
 
     assert np.all(np.diff(i_B) == 1), "Valid region not contiguous"
-    assert np.all(A["MJD"][ok_A] == B["MJD"][i_B])
+    assert np.all(A["MJD"][ok_A] == B["MJD"])
     # Check that values are copied correctly.  Since units are not
     # necessarily the same, we use allclose with very strict tolerance.
     for name in ("UT1_UTC", "PM_x", "PM_y", "dX_2000A", "dY_2000A"):
         assert_quantity_allclose(
             A[name][ok_A],
-            B[name][i_B],
+            B[name],
             rtol=1e-15,
             err_msg=(
                 f"Bug #9206 IERS B parameter {name} not copied over "