{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Deriving the Coriolis Frequency Over the Globe\n\nThis code computes the Coriolis frequency and stores it in a cube with\nassociated metadata. It then plots the Coriolis frequency on an orthographic\nprojection.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import cartopy.crs as ccrs\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nimport iris\nfrom iris.coord_systems import GeogCS\nimport iris.plot as iplt\n\n\ndef main():\n    # Start with arrays for latitudes and longitudes, with a given number of\n    # coordinates in the arrays.\n    coordinate_points = 200\n    longitudes = np.linspace(-180.0, 180.0, coordinate_points)\n    latitudes = np.linspace(-90.0, 90.0, coordinate_points)\n    lon2d, lat2d = np.meshgrid(longitudes, latitudes)\n\n    # Omega is the Earth's rotation rate, expressed in radians per second\n    omega = 7.29e-5\n\n    # The data for our cube is the Coriolis frequency,\n    # `f = 2 * omega * sin(phi)`, which is computed for each grid point over\n    # the globe from the 2-dimensional latitude array.\n    data = 2.0 * omega * np.sin(np.deg2rad(lat2d))\n\n    # We now need to define a coordinate system for the plot.\n    # Here we'll use GeogCS; 6371229 is the radius of the Earth in metres.\n    cs = GeogCS(6371229)\n\n    # The Iris coords module turns the latitude list into a coordinate array.\n    # Coords then applies an appropriate standard name and unit to it.\n    lat_coord = iris.coords.DimCoord(\n        latitudes, standard_name=\"latitude\", units=\"degrees\", coord_system=cs\n    )\n\n    # The above process is repeated for the longitude coordinates.\n    lon_coord = iris.coords.DimCoord(\n        longitudes, standard_name=\"longitude\", units=\"degrees\", coord_system=cs\n    )\n\n    # Now we add bounds to our latitude and longitude coordinates.\n    # We want simple, contiguous bounds for our regularly-spaced coordinate\n    # points so we use the guess_bounds() method of the coordinate. For more\n    # complex coordinates, we could derive and set the bounds manually.\n    lat_coord.guess_bounds()\n    lon_coord.guess_bounds()\n\n    # Now we input our data array into the cube.\n    new_cube = iris.cube.Cube(\n        data,\n        standard_name=\"coriolis_parameter\",\n        units=\"s-1\",\n        dim_coords_and_dims=[(lat_coord, 0), (lon_coord, 1)],\n    )\n\n    # Now let's plot our cube, along with coastlines, a title and an\n    # appropriately-labelled colour bar:\n    ax = plt.axes(projection=ccrs.Orthographic())\n    ax.coastlines(resolution=\"10m\")\n    mesh = iplt.pcolormesh(new_cube, cmap=\"seismic\")\n    tick_levels = [-0.00012, -0.00006, 0.0, 0.00006, 0.00012]\n    plt.colorbar(\n        mesh,\n        orientation=\"horizontal\",\n        label=\"s-1\",\n        ticks=tick_levels,\n        format=\"%.1e\",\n    )\n    plt.title(\"Coriolis frequency\")\n    plt.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.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}