FIX: Remove obvious solar radio burst F10.7 data

There are some anomalous F10.7 values over 400 which were happening
during a solar radio burst and don't allow the model to run when values
are over 750 (MSIS returns nans for everything after that point).
For now, lets flag values over 400 and fill it with F10.7a values
on that same date. We can change this later if there is a more definitive
way to justify these values, but this at least allows users to run over
the entire dataset without MSIS crashing.

pymsis/utils.py

@@ -115,7 +115,13 @@ def _load_f107_ap_data() -> dict[str, npt.NDArray]:
     ap[ap < 0] = np.nan
     daily_ap[daily_ap < 0] = np.nan
     # There are also some non-physical f10.7 values
-    arr["f107"][arr["f107"] <= 0] = np.nan
+    # F10.7 > 400 is unrealistic indicating a solar radio burst
+    solar_burst = 400
+    bad_f107 = (arr["f107"] <= 0) | (arr["f107"] > solar_burst)
+    # Set the data to the 81-day average as a temporary fix
+    arr["f107"][bad_f107] = arr["f107a"][bad_f107]
+    # flag it as interpolated or predicted values so we warn the user
+    arr["f107-type"][bad_f107] = b"INT"
 
     ap_data = np.ones((len(ap), 7)) * np.nan
     # daily Ap
@@ -148,6 +154,9 @@ def _load_f107_ap_data() -> dict[str, npt.NDArray]:
     interpolated = b"INT"
     predicted = b"PRD"
     warn_data = (arr["f107-type"] == interpolated) | (arr["f107-type"] == predicted)
+    # Because we use the F107 from the day before, we need to shift the warning
+    # to the following day when we would actually use the value
+    warn_data[1:] = warn_data[:-1]
 
     # Set the global module-level data variable
     data = {

tests/f107_ap_test_data.txt

@@ -362,7 +362,7 @@ DATE,BSRN,ND,KP1,KP2,KP3,KP4,KP5,KP6,KP7,KP8,KP_SUM,AP1,AP2,AP3,AP4,AP5,AP6,AP7,
 2000-12-26,2285,16,17,13,10,10,3,10,20,27,110,6,5,4,4,2,4,7,12,6,0.2,1,166,188.8,182.6,OBS,173.6,172.9,168.4,169.3
 2000-12-27,2285,17,20,13,10,10,23,33,23,17,150,7,5,4,4,9,18,9,6,8,0.4,2,162,187.6,181.4,OBS,173.8,173.3,168.5,169.6
 2000-12-28,2285,18,23,10,13,10,30,13,10,3,113,9,4,5,4,15,5,4,2,6,0.3,1,162,185.4,179.3,OBS,173.9,173.8,168.6,170.0
-2000-12-29,2285,19,10,33,17,10,3,10,13,20,117,4,18,6,4,2,4,5,7,6,0.3,1,151,181.5,175.5,OBS,173.7,174.3,168.4,170.4
+2000-12-29,2285,19,10,33,17,10,3,10,13,20,117,4,18,6,4,2,4,5,7,6,0.3,1,151,999,175.5,OBS,173.7,174.3,168.4,170.4
 2000-12-30,2285,20,13,7,7,10,10,7,0,10,63,5,3,3,4,4,3,0,4,3,0.1,0,152,182.1,176.1,INT,173.5,174.8,168.2,170.9
 2000-12-31,2285,21,13,0,3,10,7,0,0,3,37,5,0,2,4,3,0,0,2,2,0.0,0,119,169.5,163.9,PRD,173.3,175.0,168.1,171.0
 2040-12-01,2825,21,,,,,,,,,,,,,,,,,,,,,0,69.7,67.73,PRM,68.5,68.3,66.8,67.7

tests/test_utils.py

@@ -79,8 +79,17 @@ def test_loading_data(monkeypatch, tmp_path):
             [186.3, 186.3],
             [[7, 4, 5, 9, 4, 5.25, 5.75], [7, 4, 5, 9, 4, 5.25, 5.75]],
         ),
+        # Bad F10.7 value artificially input on 2000-12-19,
+        # should be replaced with F10.7a value on that same day
+        (
+            np.datetime64("2000-12-30T12:00"),
+            [173.7],
+            [173.5],
+            [3, 4, 4, 3, 3, 6.375, 5.375],
+        ),
     ],
 )
+@pytest.mark.filterwarnings("ignore:There is data that was either interpolated")
 def test_get_f107_ap(dates, expected_f107, expected_f107a, expected_ap):
     f107, f107a, ap = utils.get_f107_ap(dates)
     assert_array_equal(f107, expected_f107)
@@ -111,7 +120,8 @@ def test_get_f107_ap_out_of_range(dates):
 @pytest.mark.parametrize(
     "dates",
     [
-        (np.datetime64("2000-12-30T00:00")),
+        (np.datetime64("2000-12-30T00:00")),  # Bad data inserted
+        (np.datetime64("2000-12-31T00:00")),  # Interpolated data
         ([np.datetime64("2000-12-01T00:00"), np.datetime64("2000-12-31T00:00")]),
     ],
 )