{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "722c095d-cc78-4683-b73f-1c57aa78940e",
      "metadata": {},
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "ea5821e8-9579-4537-aa51-49e736c2aa2f",
      "metadata": {},
      "outputs": [],
      "source": [
        "from marine_qc import (\n",
        "    Climatology,\n",
        "    do_bayesian_buddy_check,\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e932294d-d925-4a69-9b89-5dea456a9a13",
      "metadata": {},
      "outputs": [],
      "source": [
        "from data import get_climatology_data, get_grouped_data"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0f6009ea-ba4e-4098-82df-657958aa806a",
      "metadata": {},
      "source": [
        "# How to use quality control checks with grouped reports"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "74f5740d-a913-4179-965b-93911e8338eb",
      "metadata": {},
      "source": [
        "The ``marine_qc`` toolbox provides some quality control (QC) checks that work a group of marine reports. Creating some test dataset we can show how to use them on \"real\" data. There is one check that is shown here:\n",
        "\n",
        "* a buddy check that tests that each observation is reasonably close to the average of its near neighbours in time and space\n",
        "\n",
        "For more information of all available QC checks on a group of individual reports see [the overview of QC functions for grouped reports](https://marine-qc.readthedocs.io/en/stable/overview_grp.html)."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "dbb94aa5-e9d5-4c70-9e40-44eed8bbd20b",
      "metadata": {},
      "source": [
        "The test dataset we provide is a ``pandas.DataFrame`` including marine reports representing the track of five ships. The QC functions return a QC flag for each individual report. The flags are:\n",
        "\n",
        "* `0`: the check has passed\n",
        "* `1`: the check has failed\n",
        "* `2`: the check is untestable"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "3244744c-acd5-4c1e-bd9c-cd6769721dbe",
      "metadata": {},
      "outputs": [],
      "source": [
        "data = get_grouped_data()\n",
        "data"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c047af51-ef10-40f9-860c-87b998e9760a",
      "metadata": {},
      "source": [
        "The ship data include six different parameters (`lat`, `lon`, `date`, and `sst`)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "52bffd46-9270-4568-9d5b-0cfe29432b46",
      "metadata": {},
      "outputs": [],
      "source": [
        "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5), constrained_layout=True)\n",
        "\n",
        "colors = plt.cm.tab10.colors\n",
        "\n",
        "i = 0\n",
        "for color, (platform, group) in zip(colors, data.groupby(\"platform\"), strict=False):\n",
        "    ax1.plot(\n",
        "        group[\"lon\"],\n",
        "        group[\"lat\"],\n",
        "        \"-o\",\n",
        "        ms=5,\n",
        "        label=platform,\n",
        "        color=color,\n",
        "    )\n",
        "\n",
        "    ax2.plot(\n",
        "        group[\"date\"],\n",
        "        group[\"sst\"],\n",
        "        \"-o\",\n",
        "        ms=5,\n",
        "        label=platform,\n",
        "        color=color,\n",
        "    )\n",
        "\n",
        "    outlier = data.loc[(data.platform == platform) & (data.date == pd.Timestamp(f\"2026-07-01 {13 + i}:00\"))]\n",
        "\n",
        "    ax1.scatter(\n",
        "        outlier.lon,\n",
        "        outlier.lat,\n",
        "        s=150,\n",
        "        facecolors=\"none\",\n",
        "        edgecolors=\"red\",\n",
        "        linewidth=2,\n",
        "        zorder=10,\n",
        "    )\n",
        "\n",
        "    ax2.scatter(\n",
        "        outlier.date,\n",
        "        outlier.sst,\n",
        "        s=150,\n",
        "        facecolors=\"none\",\n",
        "        edgecolors=\"red\",\n",
        "        linewidth=2,\n",
        "        zorder=10,\n",
        "    )\n",
        "\n",
        "    i += 1\n",
        "\n",
        "ax1.set_xlabel(\"Longitude (°)\")\n",
        "ax1.set_ylabel(\"Latitude (°)\")\n",
        "ax1.set_title(\"Buddy observations\")\n",
        "ax1.legend(loc=\"best\")\n",
        "\n",
        "ax2.set_xlabel(\"Time\")\n",
        "ax2.set_ylabel(\"SST (°C)\")\n",
        "ax2.set_title(\"Sea surface temperature\")\n",
        "ax2.tick_params(axis=\"x\", rotation=30)\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b6f59ff4-d2f8-4689-b7df-7427b595669b",
      "metadata": {},
      "source": [
        "We see five different ships. Each of them includes one value that is much higher than the others. These values should be detected and flagged as `1` (failed)."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6501f9c6-7ac7-4fc8-9196-8a0c4e445499",
      "metadata": {},
      "source": [
        "Firtly, we import some external data. This is a `xarray.Dataset` containing [CF](https://cfconventions.org/)-compliant data:\n",
        "\n",
        "* a land-sea mask where `1` denotes a land point and `0` denotes a sea point (`land_sea_mask`)\n",
        "* a simple, latitudinally and longitudinally dependent sea-surface temperature field (`sst`)\n",
        "* a sea-surface temperature standard deviation of grid cell minus neighbourhood mean (`sst_stde1`)\n",
        "* a sea-surface temperature standard deviation of point observation minus grid cell mean (`sst_stde2`)\n",
        "* a sea-surface temperature standard deviation of neighbourhood mean uncertainty (`sst_stde3`)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "05c15942-3a83-41e2-b4d0-220947661a16",
      "metadata": {},
      "outputs": [],
      "source": [
        "climatology_data = get_climatology_data()\n",
        "climatology_data"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4a4c9b93-7589-4114-bbf7-da7768524ea7",
      "metadata": {},
      "source": [
        "For the buddy check, we need some external data:\n",
        "\n",
        "* `climatology`: The climatological average(s) used to calculate anomalies.\n",
        "* `stdev1`: Field of standard deviations representing standard deviation of difference between target gridcell and complete neighbour average (grid area to neighbourhood difference)\n",
        "* `stdev2`: Field of standard deviations representing standard deviation of difference between a single observation and the target gridcell average (point to grid area difference)\n",
        "* `stdev3`: Field of standard deviations representing standard deviation of difference between random neighbour gridcell and full neighbour average (uncertainty in neighbour average)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "d563eabe-d8db-41f2-88e6-ac036047bc39",
      "metadata": {},
      "outputs": [],
      "source": [
        "fig, axes = plt.subplots(\n",
        "    2,\n",
        "    2,\n",
        "    figsize=(14, 8),\n",
        "    constrained_layout=True,\n",
        ")\n",
        "\n",
        "plots = [\n",
        "    (\"sst\", \"sea-surface temperature\", \"RdYlBu_r\"),\n",
        "    (\"sst_stdev1\", \"stdev1\", \"viridis\"),\n",
        "    (\"sst_stdev2\", \"stdev2\", \"viridis\"),\n",
        "    (\"sst_stdev3\", \"stdev3\", \"viridis\"),\n",
        "]\n",
        "\n",
        "for ax, (var, title, cmap) in zip(axes.flat, plots, strict=False):\n",
        "    climatology_data[var].isel(time=0).plot(\n",
        "        ax=ax,\n",
        "        cmap=cmap,\n",
        "        add_colorbar=True,\n",
        "    )\n",
        "    ax.set_title(title)\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b0df8f00-c509-4762-9c11-0b0ef397df07",
      "metadata": {},
      "source": [
        "Unfortunately, the standard deviations have to be converted manually to `marine_qc.Climatology` classes."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "0500d952-6ae4-427d-8490-71fad195d1ca",
      "metadata": {},
      "outputs": [],
      "source": [
        "sst_stdev1 = Climatology(climatology_data.sst_stdev1)\n",
        "sst_stdev2 = Climatology(climatology_data.sst_stdev2)\n",
        "sst_stdev3 = Climatology(climatology_data.sst_stdev3)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a4bc4780-0c3d-44c7-8e08-94aa6fe25301",
      "metadata": {},
      "source": [
        "Besiedes the external climatologies some additional parameters needs to be set:\n",
        "\n",
        "* `prior_probability_of_gross_error`: prior probability of gross error, which is the background rate of gross errors (`0.05`)\n",
        "* `quantization_interval`: smallest possible increment in the input values (`0.1`)\n",
        "* `one_sigma_measurement_uncertainty`: estimated one sigma measurement uncertainty (`1.0`)\n",
        "* `limits`: list with three members which specify the search range for the buddy check (`[2, 2, 4]`)\n",
        "* `noise_scaling`: tuning parameter used to multiply `stdev2` (`3.0`)\n",
        "* `maximum_anomaly`: largest absolute anomaly, assumes that the maximum and minimum anomalies have the same magnitude (`8.0`)\n",
        "* `fail_probability`: probability of gross error that corresponds to a failed test. Anything with a probability of gross error greater than fail_probability will be considered failing (`0.3`)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "55d0c7fe-9f31-4740-8a84-7b8eedbca766",
      "metadata": {},
      "outputs": [],
      "source": [
        "qc_buddy = do_bayesian_buddy_check(\n",
        "    value=data.sst,\n",
        "    lat=data.lat,\n",
        "    lon=data.lon,\n",
        "    date=data.date,\n",
        "    climatology=climatology_data.sst,\n",
        "    stdev1=sst_stdev1,\n",
        "    stdev2=sst_stdev2,\n",
        "    stdev3=sst_stdev3,\n",
        "    prior_probability_of_gross_error=0.05,\n",
        "    quantization_interval=0.1,\n",
        "    one_sigma_measurement_uncertainty=1.0,\n",
        "    limits=[2, 2, 4],\n",
        "    noise_scaling=3.0,\n",
        "    maximum_anomaly=8.0,\n",
        "    fail_probability=0.3,\n",
        ")\n",
        "pd.DataFrame({**data, \"qc_flag\": qc_buddy})"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "51380051-d3da-4cf4-9d76-d2c9a1714273",
      "metadata": {},
      "source": [
        "We get the expected results."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b0302488-db38-401f-a944-4e5f7f969b3d",
      "metadata": {},
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.13.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}