{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Global Average Annual Temperature Maps\n\nProduces maps of global temperature forecasts from the A1B and E1 scenarios.\n\nThe data used comes from the HadGEM2-AO model simulations for the A1B and E1\nscenarios, both of which were derived using the IMAGE Integrated Assessment\nModel (Johns et al. 2011; Lowe et al. 2009).\n\n## References\n\n    Johns T.C., et al. (2011) Climate change under aggressive mitigation: the\n    ENSEMBLES multi-model experiment. Climate Dynamics, Vol 37, No. 9-10,\n    doi:10.1007/s00382-011-1005-5.\n\n    Lowe J.A., C.D. Hewitt, D.P. Van Vuuren, T.C. Johns, E. Stehfest, J-F.\n    Royer, and P. van der Linden, 2009. New Study For Climate Modeling,\n    Analyses, and Scenarios. Eos Trans. AGU, Vol 90, No. 21,\n    doi:10.1029/2009EO210001.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os.path\n\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nimport iris\nimport iris.coords as coords\nimport iris.plot as iplt\n\n\ndef cop_metadata_callback(cube, field, filename):\n    \"\"\"\n    A function which adds an \"Experiment\" coordinate which comes from the\n    filename.\n    \"\"\"\n\n    # Extract the experiment name (such as A1B or E1) from the filename (in\n    # this case it is just the start of the file name, before the first \".\").\n    fname = os.path.basename(filename)  # filename without path.\n    experiment_label = fname.split(\".\")[0]\n\n    # Create a coordinate with the experiment label in it...\n    exp_coord = coords.AuxCoord(\n        experiment_label, long_name=\"Experiment\", units=\"no_unit\"\n    )\n\n    # ...and add it to the cube.\n    cube.add_aux_coord(exp_coord)\n\n\ndef main():\n    # Load E1 and A1B scenarios using the callback to update the metadata.\n    scenario_files = [\n        iris.sample_data_path(fname) for fname in [\"E1.2098.pp\", \"A1B.2098.pp\"]\n    ]\n    scenarios = iris.load(scenario_files, callback=cop_metadata_callback)\n\n    # Load the preindustrial reference data.\n    preindustrial = iris.load_cube(iris.sample_data_path(\"pre-industrial.pp\"))\n\n    # Define evenly spaced contour levels: -2.5, -1.5, ... 15.5, 16.5 with the\n    # specific colours.\n    levels = np.arange(20) - 2.5\n    red = (\n        np.array(\n            [\n                0,\n                0,\n                221,\n                239,\n                229,\n                217,\n                239,\n                234,\n                228,\n                222,\n                205,\n                196,\n                161,\n                137,\n                116,\n                89,\n                77,\n                60,\n                51,\n            ]\n        )\n        / 256.0\n    )\n    green = (\n        np.array(\n            [\n                16,\n                217,\n                242,\n                243,\n                235,\n                225,\n                190,\n                160,\n                128,\n                87,\n                72,\n                59,\n                33,\n                21,\n                29,\n                30,\n                30,\n                29,\n                26,\n            ]\n        )\n        / 256.0\n    )\n    blue = (\n        np.array(\n            [\n                255,\n                255,\n                243,\n                169,\n                99,\n                51,\n                63,\n                37,\n                39,\n                21,\n                27,\n                23,\n                22,\n                26,\n                29,\n                28,\n                27,\n                25,\n                22,\n            ]\n        )\n        / 256.0\n    )\n\n    # Put those colours into an array which can be passed to contourf as the\n    # specific colours for each level.\n    colors = np.stack([red, green, blue], axis=1)\n\n    # Make a wider than normal figure to house two maps side-by-side.\n    fig, ax_array = plt.subplots(1, 2, figsize=(12, 5))\n\n    # Loop over our scenarios to make a plot for each.\n    for ax, experiment, label in zip(\n        ax_array, [\"E1\", \"A1B\"], [\"E1\", \"A1B-Image\"]\n    ):\n        exp_cube = scenarios.extract_cube(\n            iris.Constraint(Experiment=experiment)\n        )\n        time_coord = exp_cube.coord(\"time\")\n\n        # Calculate the difference from the preindustial control run.\n        exp_anom_cube = exp_cube - preindustrial\n\n        # Plot this anomaly.\n        plt.sca(ax)\n        ax.set_title(f\"HadGEM2 {label} Scenario\", fontsize=10)\n        contour_result = iplt.contourf(\n            exp_anom_cube, levels, colors=colors, extend=\"both\"\n        )\n        plt.gca().coastlines()\n\n    # Now add a colourbar who's leftmost point is the same as the leftmost\n    # point of the left hand plot and rightmost point is the rightmost\n    # point of the right hand plot.\n\n    # Get the positions of the 2nd plot and the left position of the 1st plot.\n    left, bottom, width, height = ax_array[1].get_position().bounds\n    first_plot_left = ax_array[0].get_position().bounds[0]\n\n    # The width of the colorbar should now be simple.\n    width = left - first_plot_left + width\n\n    # Add axes to the figure, to place the colour bar.\n    colorbar_axes = fig.add_axes([first_plot_left, 0.18, width, 0.03])\n\n    # Add the colour bar.\n    cbar = plt.colorbar(\n        contour_result, colorbar_axes, orientation=\"horizontal\"\n    )\n\n    # Label the colour bar and add ticks.\n    cbar.set_label(preindustrial.units)\n    cbar.ax.tick_params(length=0)\n\n    # Get the time datetime from the coordinate.\n    time = time_coord.units.num2date(time_coord.points[0])\n    # Set a title for the entire figure, using the year from the datetime\n    # object. Also, set the y value for the title so that it is not tight to\n    # the top of the plot.\n    fig.suptitle(\n        f\"Annual Temperature Predictions for {time.year}\",\n        y=0.9,\n        fontsize=18,\n    )\n\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
}