An example of fixing nearest neighbour interpolation in iris

nearest_interpolation.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:ce0f3b988dd7e3d3fcbd212fc16cd56431d210fc7353bbcd8fe3a98ddf5568da"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Question about nearest neighbour interpolation being different with different iris functions.\n",
      "\n",
      "See https://groups.google.com/forum/#!topic/scitools-iris/tF-LgxoIBoI."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from iris.cube import Cube\n",
      "from iris.coords import DimCoord\n",
      "from iris.analysis import Nearest\n",
      "from iris.analysis.interpolate import nearest_neighbour_data_value"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "latitude = DimCoord(np.arange(-10, 11, 5), standard_name='latitude', units='degrees')\n",
      "longitude = DimCoord(np.arange(-5, 6, 5), standard_name='longitude', units='degrees')\n",
      "data = np.reshape(np.arange(15)+1.0,(5,3))\n",
      "cube = Cube(data, dim_coords_and_dims=[(latitude, 0), (longitude, 1)])\n",
      "\n",
      "print cube\n",
      "print longitude\n",
      "print latitude\n",
      "print data"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "unknown / (unknown)                 (latitude: 5; longitude: 3)\n",
        "     Dimension coordinates:\n",
        "          latitude                           x             -\n",
        "          longitude                          -             x\n",
        "DimCoord(array([-5,  0,  5]), standard_name='longitude', units=Unit('degrees'))\n",
        "DimCoord(array([-10,  -5,   0,   5,  10]), standard_name='latitude', units=Unit('degrees'))\n",
        "[[  1.   2.   3.]\n",
        " [  4.   5.   6.]\n",
        " [  7.   8.   9.]\n",
        " [ 10.  11.  12.]\n",
        " [ 13.  14.  15.]]\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sample_points = [[['latitude', 2.51], ['longitude', 2.5]],\n",
      "                 [['latitude', -2.51], ['longitude', 2.5]],\n",
      "                 [['latitude', 2.5], ['longitude', -2.51]],\n",
      "                 [['latitude', -2.51], ['longitude', -2.51]]]\n",
      "\n",
      "for point in sample_points:\n",
      "    print 'New interface: {}'.format(cube.interpolate(point,scheme=Nearest()).data)\n",
      "    print 'Original: {}'.format(nearest_neighbour_data_value(cube, point))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "New interface: 8.0\n",
        "Original: 11.0\n",
        "New interface: 8.0\n",
        "Original: 5.0\n",
        "New interface: 8.0\n",
        "Original: 7.0\n",
        "New interface: 8.0\n",
        "Original: 4.0\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If we set the coordinates to be floating point values though: "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "latitude = DimCoord(np.arange(-10., 11, 5), standard_name='latitude', units='degrees')\n",
      "longitude = DimCoord(np.arange(-5., 6, 5), standard_name='longitude', units='degrees')\n",
      "data = np.reshape(np.arange(15)+1.0,(5,3))\n",
      "cube = Cube(data, dim_coords_and_dims=[(latitude, 0), (longitude, 1)])\n",
      "\n",
      "print cube\n",
      "print longitude\n",
      "print latitude\n",
      "print data"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "unknown / (unknown)                 (latitude: 5; longitude: 3)\n",
        "     Dimension coordinates:\n",
        "          latitude                           x             -\n",
        "          longitude                          -             x\n",
        "DimCoord(array([-5.,  0.,  5.]), standard_name='longitude', units=Unit('degrees'))\n",
        "DimCoord(array([-10.,  -5.,   0.,   5.,  10.]), standard_name='latitude', units=Unit('degrees'))\n",
        "[[  1.   2.   3.]\n",
        " [  4.   5.   6.]\n",
        " [  7.   8.   9.]\n",
        " [ 10.  11.  12.]\n",
        " [ 13.  14.  15.]]\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for point in sample_points:\n",
      "    print 'New interface: {}'.format(cube.interpolate(point,scheme=Nearest()).data)\n",
      "    print 'Original: {}'.format(nearest_neighbour_data_value(cube, point))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "New interface: 11.0\n",
        "Original: 11.0\n",
        "New interface: 5.0\n",
        "Original: 5.0\n",
        "New interface: 7.0\n",
        "Original: 7.0\n",
        "New interface: 4.0\n",
        "Original: 4.0\n"
       ]
      }
     ],
     "prompt_number": 5
    }
   ],
   "metadata": {}
  }
 ]
}