.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples/plot_robust_kalman.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_plot_robust_kalman.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_plot_robust_kalman.py:


Robust fitting and temporal smoothing
-------------------------------------

A large disturbance at a single station -- a spike, a baseline jump,
local interference -- is inconsistent with any smooth ionospheric
current system, but an ordinary least-squares fit will spread it across
the entire interpolated map. This example compares three defenses:

1. the ordinary snapshot fit (contaminated),
2. robust IRLS weighting (``fit(robust='bisquare')``), which rejects
   the outlier spatially at each time step, and
3. :class:`pysecs.KalmanSECS` with innovation gating, which couples the
   time steps and rejects observations inconsistent with the recent
   past.

.. GENERATED FROM PYTHON SOURCE LINES 17-79



.. image-sg:: /examples/images/sphx_glr_plot_robust_kalman_001.png
   :alt: Validation point response to a single-station spike at t=30
   :srcset: /examples/images/sphx_glr_plot_robust_kalman_001.png
   :class: sphx-glr-single-img





.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np

    from pysecs import SECS, KalmanSECS, T_df


    R_E = 6378e3
    rng = np.random.default_rng(11)

    # A small grid of divergence-free SECs and a ground station network
    lat_g, lon_g = np.meshgrid(np.linspace(-3, 3, 4), np.linspace(-3, 3, 4), indexing="ij")
    sec_locs = np.column_stack([lat_g.ravel(), lon_g.ravel(), np.full(16, R_E + 110e3)])
    nobs = 12
    obs_locs = np.column_stack(
        [rng.uniform(-4, 4, nobs), rng.uniform(-4, 4, nobs), np.full(nobs, R_E)]
    )

    # Static true current system observed with noise for one hour at 1 min
    true_amps = rng.normal(0, 1e4, 16)
    B_true = np.tensordot(true_amps, T_df(obs_locs, sec_locs), (0, 2))
    B_scale = np.max(np.abs(B_true))

    ntimes = 60
    times = np.arange(ntimes) * 60.0
    noise = 0.05 * B_scale
    B_obs = B_true[np.newaxis] + rng.normal(0, noise, (ntimes, nobs, 3))
    obs_std = np.full_like(B_obs, noise)

    # A 50x spike at station 0 halfway through the series
    B_obs[30, 0, :] += 50 * B_scale

    # Predict at a validation location far from the spiked station
    pred_loc = np.array([[1.5, 1.5, R_E]])
    B_val_true = np.tensordot(true_amps, T_df(pred_loc, sec_locs), (0, 2))[0]

    # 1. Ordinary snapshot fit
    snap = SECS(sec_df_loc=sec_locs)
    snap.fit(obs_locs, B_obs, obs_std=obs_std, epsilon=1e-6)
    pred_snap = snap.predict_B(pred_loc)

    # 2. Robust IRLS fit
    robust = SECS(sec_df_loc=sec_locs)
    robust.fit(obs_locs, B_obs, obs_std=obs_std, epsilon=1e-6, robust="bisquare")
    pred_robust = robust.predict_B(pred_loc)

    # 3. Kalman smoother with innovation gating
    kalman = KalmanSECS(SECS(sec_df_loc=sec_locs), tau=3600.0, gate_sigma=4.0)
    kalman.fit(obs_locs, B_obs, times, obs_std=obs_std)
    pred_kalman = kalman.predict_B(pred_loc)

    fig, ax = plt.subplots()
    minutes = times / 60
    ax.plot(minutes, pred_snap[:, 0] * 1e9, label="snapshot fit", alpha=0.8)
    ax.plot(minutes, pred_robust[:, 0] * 1e9, label="robust (bisquare)", alpha=0.8)
    ax.plot(minutes, pred_kalman[:, 0] * 1e9, label="Kalman + gating", alpha=0.8)
    ax.axhline(B_val_true[0] * 1e9, color="k", ls=":", label="truth")
    ax.set_xlabel("time (minutes)")
    ax.set_ylabel(r"predicted B$_x$ (nT)")
    ax.set_title("Validation point response to a single-station spike at t=30")
    ax.legend()
    plt.show()


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.244 seconds)


.. _sphx_glr_download_examples_plot_robust_kalman.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_robust_kalman.ipynb <plot_robust_kalman.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_robust_kalman.py <plot_robust_kalman.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_robust_kalman.zip <plot_robust_kalman.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_