{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Calculating a Custom Statistic\n\nThis example shows how to define and use a custom\n:class:`iris.analysis.Aggregator`, that provides a new statistical operator for\nuse with cube aggregation functions such as :meth:`~iris.cube.Cube.collapsed`,\n:meth:`~iris.cube.Cube.aggregated_by` or\n:meth:`~iris.cube.Cube.rolling_window`.\n\nIn this case, we have a 240-year sequence of yearly average surface temperature\nover North America, and we want to calculate in how many years these exceed a\ncertain temperature over a spell of 5 years or more.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport iris\nfrom iris.analysis import Aggregator\nimport iris.plot as iplt\nimport iris.quickplot as qplt\nfrom iris.util import rolling_window\n\n\n# Define a function to perform the custom statistical operation.\n# Note: in order to meet the requirements of iris.analysis.Aggregator, it must\n# do the calculation over an arbitrary (given) data axis.\ndef count_spells(data, threshold, axis, spell_length):\n    \"\"\"\n    Function to calculate the number of points in a sequence where the value\n    has exceeded a threshold value for at least a certain number of timepoints.\n\n    Generalised to operate on multiple time sequences arranged on a specific\n    axis of a multidimensional array.\n\n    Args:\n\n    * data (array):\n        raw data to be compared with value threshold.\n\n    * threshold (float):\n        threshold point for 'significant' datapoints.\n\n    * axis (int):\n        number of the array dimension mapping the time sequences.\n        (Can also be negative, e.g. '-1' means last dimension)\n\n    * spell_length (int):\n        number of consecutive times at which value > threshold to \"count\".\n\n    \"\"\"\n    if axis < 0:\n        # just cope with negative axis numbers\n        axis += data.ndim\n    # Threshold the data to find the 'significant' points.\n    data_hits = data > threshold\n    # Make an array with data values \"windowed\" along the time axis.\n    hit_windows = rolling_window(data_hits, window=spell_length, axis=axis)\n    # Find the windows \"full of True-s\" (along the added 'window axis').\n    full_windows = np.all(hit_windows, axis=axis + 1)\n    # Count points fulfilling the condition (along the time axis).\n    spell_point_counts = np.sum(full_windows, axis=axis, dtype=int)\n    return spell_point_counts\n\n\ndef main():\n    # Load the whole time-sequence as a single cube.\n    file_path = iris.sample_data_path(\"E1_north_america.nc\")\n    cube = iris.load_cube(file_path)\n\n    # Make an aggregator from the user function.\n    SPELL_COUNT = Aggregator(\n        \"spell_count\", count_spells, units_func=lambda units: 1\n    )\n\n    # Define the parameters of the test.\n    threshold_temperature = 280.0\n    spell_years = 5\n\n    # Calculate the statistic.\n    warm_periods = cube.collapsed(\n        \"time\",\n        SPELL_COUNT,\n        threshold=threshold_temperature,\n        spell_length=spell_years,\n    )\n    warm_periods.rename(\"Number of 5-year warm spells in 240 years\")\n\n    # Plot the results.\n    qplt.contourf(warm_periods, cmap=\"RdYlBu_r\")\n    plt.gca().coastlines()\n    iplt.show()\n\n\nif __name__ == \"__main__\":\n    main()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.10.8"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}