{
 "metadata": {
  "name": "area_weights"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Shapefile area averaging\n",
      "\n",
      "### Introduction\n",
      "\n",
      "This example demonstrates how to find the area weighting of a cube's data points that fall within a given geometry from a shapefile. \n",
      "\n",
      "Area weighting allows us to account for data cells that only fall partially within a given geometry, by determining the percentage of the cell's area that falls within the geometry.\n",
      "\n",
      "We can use that percentage to define its overall contribution when we perform a mathematical operation on the values within the geometry (for example by using a function from `iris.analysis` that may be used in conjunction with `cube.collapsed`).\n",
      "\n",
      "Let's look at this pictorially:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "![Illustration of a shapefile geometry with cube cells overlapping as indicated by the grid](files/img/area_weighting.png)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the image the black grid is the grid of data cells from the cube, the blue line is the extent of the shapefile geometry and the green filled area is the area of the grid that sits within the shapefile. Clearly there are numerous cells toward the centre of the grid that entirely sit within the geometry, so the values of these cells will contribute 100% to the mathematical operation to be performed. However the top centre cell only overlaps the geometry by about 10%, so will only contribute 10% of its value, and so on."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "An example of when this might be valuable is determining the total rainfall within a UK county over a given time period. This would be done by collapsing the cube over the time dimension and finding the sum of all the values that lie within the geometry of the county we are interested in. As in the image above, not all cells will be fully contained within the geometry and so in order to determine the amount of rainfall within the county's geometry, we need to find the area-weighted amount of total rainfall for that cell that fell within the geometry."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Worked Example\n",
      "\n",
      "In this example we will determine the area-weighted average air temperature of the UK in the year 2098 for one of our climate future scenarios (A1B)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's start by importing the libraries we will need for this example, and checking the version of Iris being used."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import cartopy.io.shapereader as shpreader\n",
      "import iris\n",
      "import iris.analysis.geometry as iag"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print('Iris version: {}'.format(iris.__version__))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Iris version: 1.5.1\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we will import some data; namely the A1B scenario global air temperatures file."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a1b_cube = iris.load_cube(iris.sample_data_path('A1B.2098.pp'))\n",
      "lats = a1b_cube.coord('latitude')\n",
      "lons = a1b_cube.coord('longitude')\n",
      "if not lats.has_bounds():\n",
      "    lats.guess_bounds()\n",
      "if not lons.has_bounds():\n",
      "    lons.guess_bounds()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Calculating the geometry weights will require known bounds to our grid cells so if either of our cube's lateral coordinates are not bounded, we use the `guess_bounds` method to set appropriate bounds on each lateral coordinate.\n",
      "\n",
      "Now we will import a shapefile to use, specifically the Natural Earth global country boundaries shapefile."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "filename = shpreader.natural_earth(resolution='110m',\n",
      "                                   category='cultural',\n",
      "                                   name='admin_0_countries')\n",
      "borders = shpreader.Reader(filename)\n",
      "uk_geom, = [g for g in borders.records() if g.attributes['name_long']=='United Kingdom']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here the shapefile is specified and read using the `cartopy.io.shapereader` module. The UK's geometry is determined by iterating over all the records in the shapefile, choosing the record where the attribute `name_long` matches 'United Kingdom'.\n",
      "\n",
      "This Natural Earth shapefile was chosen here for simplicity's sake, but any other shapefile can be accessed using the methods demonstrated here. It is worth noting that whilst `shpreader.natural_earth()` will download Natural Earth shapefiles as required, all it returns is a filepath to a shapefile.\n",
      "\n",
      "A very similar process should be able to be followed with any standard shapefile, as all shapefiles should be loaded with both `geometries()` and `records()` methods. The former returns a generator object of each geometry within the shapefile while the latter returns a generator object of each record within the shapefile.\n",
      "\n",
      "A record is made up of a geometry (one of the geometries from the geometries method) and an attributes dictionary. The keys within the attributes dictionary are dependant on the shapefile loaded, but are simple to locate:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print uk_geom.attributes.keys()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "['labelrank', 'gu_a3', 'long_len', 'scalerank', 'wikipedia', 'woe_id', 'brk_group', 'pop_est', 'sovereignt', 'pop_year', 'homepart', 'region_wb', 'adm0_a3', 'continent', 'region_un', 'adm0_a3_us', 'brk_name', 'mapcolor8', 'name_sort', 'wb_a3', 'formal_en', 'sov_a3', 'tiny', 'geounit', 'adm0_a3_is', 'iso_n3', 'brk_a3', 'abbrev_len', 'subunit', 'type', 'adm0_a3_un', 'gdp_year', 'economy', 'mapcolor13', 'fips_10', 'name_alt', 'gdp_md_est', 'brk_diff', 'wb_a2', 'income_grp', 'lastcensus', 'iso_a3', 'iso_a2', 'featurecla', 'postal', 'geou_dif', 'name_len', 'adm0_a3_wb', 'name', 'su_a3', 'level', 'admin', 'mapcolor7', 'note_adm0', 'adm0_dif', 'abbrev', 'subregion', 'formal_fr', 'su_dif', 'mapcolor9', 'un_a3', 'note_brk', 'name_long']\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "These operations can be performed on any shapefile record.\n",
      "\n",
      "Note, however, that the keys within an attributes dictionary are dependent on the shapefile, and that code that expects a certain key to appear within the attributes may fail if a different shapefile is processed."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "With the UK's geometry selected from the loaded shapefile, it is ready to be passed to the area weights calculating function:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "weights = iag.geometry_area_weights(a1b_cube, uk_geom.geometry)\n",
      "cube_collapsed = a1b_cube.collapsed(['latitude', 'longitude'],\n",
      "                                    iris.analysis.MEAN,\n",
      "                                    weights=weights)\n",
      "print('Average UK air temp in 2098: {:.3f}K'.format(cube_collapsed.data[0]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Average UK air temp in 2098: 288.180K\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "These steps calculate the area-weighted values for the cube based on the UK's geometries. These are passed to the collapse call, which collapses to find the area-weighted mean over both lateral coordinates. Finally, the resultant cube's single data point (average UK air temperature in 2098 according to the A1B scenario) is printed."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Something to be aware of\n",
      "\n",
      "A note of caution to end on. The example shown here works in Euclidean space; curvature of the earth is not accounted for in the weighted area calculations performed."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}