A new vorticity module.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A new vorticity calculation for TCRM\n",
    "\n",
    "Based on Johnny Lin's 'gemath' package (http://www.johnny-lin.com/py_pkgs/gemath/doc/index.html). The existing vorticity module was "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "import sys, cStringIO\n",
    "\n",
    "from Utilities import metutils\n",
    "\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def can_use_sphere(longitude, latitude):\n",
    "    \"\"\"Test if can use sphere package.\n",
    "   \n",
    "    Calling Sequence:\n",
    "        Result = can_use_sphere(longitude, latitude)\n",
    "\n",
    "    Note that Result is a 3-element tuple, not a scalar.  See the\n",
    "    description below for details.\n",
    "\n",
    "    Test if the longitude-latitude domain and Python package configur-\n",
    "    ation allows use of the NCAR SPHEREPACK 3.0 package.  Specifically \n",
    "    the function tests:\n",
    "\n",
    "    (1) Can you import sphere?\n",
    "    (2) Is the longitude vector evenly spaced?\n",
    "    (3) Does the longitude vector entirely span the globe, but doesn't\n",
    "        repeat?\n",
    "    (4) Is the latitude vector gaussian?  If not, is it evenly spaced \n",
    "        and includes the poles?\n",
    "\n",
    "    If the answer to (1) is no, then there's no use to go through the \n",
    "    other tests and the function return tuple element [0] is 0.  If (2) \n",
    "    or (3) is no, return tuple element [0] is 0.  If (4) is both not \n",
    "    gaussian and not evenly spaced with the poles, return tuple \n",
    "    element [0] is 0.  Otherwise, return tuple element [0] is 1.  Note \n",
    "    that (4) uses the latitude checker in the sphere package.\n",
    "\n",
    "\n",
    "    Method Arguments:\n",
    "    * longitude:  Vector of longitudes of the domain [deg].  Can be \n",
    "      in any order, and a Numeric array or regular list/tuple.\n",
    "\n",
    "    * latitude:  Vector of latitudes of the domain [deg].  Can be in\n",
    "      any order, and a Numeric array or regular list/tuple.\n",
    "\n",
    "\n",
    "    Output Result:\n",
    "    * A 3-element tuple:\n",
    "      [0]:  1 if the longitude-latitude domain and Python package \n",
    "            configuration passes tests to allow use of the NCAR \n",
    "            SPHEREPACK 3.0 package.  0 if does not pass those tests.\n",
    "      [1]:  String containing messages written to stdout by calls to \n",
    "            module sphere used for checking latitude (e.g. if it is\n",
    "            gaussian, etc.).  If there are no messages or there were \n",
    "            no calls needed to sphere for checking latitude, value is \n",
    "            empty string.\n",
    "      [2]:  Same as [1] except contains stderr messages.\n",
    "\n",
    "\n",
    "    Examples:\n",
    "    >>> from can_use_sphere import can_use_sphere\n",
    "    >>> import numpy as np\n",
    "    >>> lon = np.arange(36)*10\n",
    "    >>> lat = [-90, -60, -30, 0, 30, 60, 90]\n",
    "    >>> can_use_sphere(lon, lat)[0]\n",
    "    1\n",
    "    >>> lat = [-90, -60, -30, 0, 30, 60, 87]\n",
    "    >>> can_use_sphere(lon, lat)[0]\n",
    "    0\n",
    "    >>> can_use_sphere(lon, lat)[1].splitlines()[1]\n",
    "    'CANNOT PROCESS THE DATA - Latitude values are incorrect'\n",
    "    \"\"\"\n",
    "\n",
    "\n",
    "    #- Test if the sphere package exists:\n",
    "\n",
    "    try:  import sphere\n",
    "    except ImportError:  return (0,'','')\n",
    "\n",
    "\n",
    "    #- Convert input to Numeric and sort:\n",
    "   \n",
    "    lon = np.sort(np.array(longitude))\n",
    "    lat = np.sort(np.array(latitude))\n",
    "\n",
    "\n",
    "    #- Check if either vector is less than 2 elements.  If so,\n",
    "    #  return false:\n",
    "\n",
    "    if (len(lon) < 2) or (len(lat) < 2):  return (0,'','')\n",
    "\n",
    "\n",
    "    #- Is the longitude vector evenly spaced?  If not, return false:\n",
    "\n",
    "    diff_lon = lon[1:] - lon[0:-1]\n",
    "    if not np.allclose(diff_lon, diff_lon[0]):  return (0,'','')\n",
    "\n",
    "\n",
    "    #- Does the longitude vector exactly spans the globe but without\n",
    "    #  repeating?  If not, return false:\n",
    "\n",
    "    if not np.allclose( (lon[-1]+diff_lon[0]-lon[0]), 360.0 ):\n",
    "       return (0,'','')\n",
    "\n",
    "\n",
    "    #- Check latitude (e.g. whether it is gaussian, whether it includes\n",
    "    #  the pole points if they are evenly spaced) and any other bad \n",
    "    #  conditions for sphere using the sphere package checker.  If \n",
    "    #  sphere can't execute correctly, return 0.  Also return standard \n",
    "    #  error and standard out:\n",
    "\n",
    "    temp_out = sys.stdout\n",
    "    temp_err = sys.stderr\n",
    "    sys.stdout = cStringIO.StringIO()\n",
    "    sys.stderr = cStringIO.StringIO()\n",
    "\n",
    "    try:\n",
    "        try:\n",
    "            sph_obj = sphere.Sphere(lon, lat)\n",
    "            sph_obj_exception = 0\n",
    "        except:\n",
    "            sph_obj_exception = 1\n",
    "    finally:\n",
    "        stdout = sys.stdout.getvalue()\n",
    "        stderr = sys.stderr.getvalue()\n",
    "        sys.stdout.close()\n",
    "        sys.stderr.close()\n",
    "        sys.stdout = temp_out\n",
    "        sys.stderr = temp_err\n",
    "        if sph_obj_exception == 1:  return (0, stdout, stderr)\n",
    "\n",
    "\n",
    "    #- Return true if haven't returned yet:\n",
    "\n",
    "    return (1, stdout, stderr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lat = [-90, -60, -30, 0, 30, 60, 90]\n",
    "lon = np.arange(35)*10\n",
    "print can_use_sphere(lon, lat)[0]\n",
    "lon = [-180, 0, 180]\n",
    "print can_use_sphere(lon, lat)[0]\n",
    "lon = [-180, 180]\n",
    "print can_use_sphere(lon, lat)[0]\n",
    "lon = [-180, 170]\n",
    "print can_use_sphere(lon, lat)[0]\n",
    "lat = [45]\n",
    "print can_use_sphere(lon, lat)[0]\n",
    "lon = np.arange(36)*10\n",
    "lat = [-87.86380005,-85.0965271 ,-82.31291199,-79.52560425 \\\n",
    "              ,-76.73690033,-73.9475174 , -71.15775299,-68.3677597 \\\n",
    "              ,-65.5776062 ,-62.78735352,-59.99702072, -57.20663071 \\\n",
    "              ,-54.41619873,-51.62573242,-48.83523941,-46.04472733 \\\n",
    "              ,-43.25419617,-40.46364975,-37.67308807,-34.88252258 \\\n",
    "              ,-32.09194565, -29.30135918,-26.51076889,-23.72017479 \\\n",
    "              ,-20.92957497,-18.13897133, -15.34836483,-12.55775642 \\\n",
    "              , -9.76714516, -6.97653341, -4.18592072, -1.39530694 \\\n",
    "              ,  1.39530694,  4.18592072,  6.97653341,  9.76714516 \\\n",
    "              , 12.55775642, 15.34836483, 18.13897133, 20.92957497 \\\n",
    "              , 23.72017479, 26.51076889, 29.30135918, 32.09194565 \\\n",
    "              , 34.88252258, 37.67308807, 40.46364975, 43.25419617 \\\n",
    "              , 46.04472733, 48.83523941, 51.62573242, 54.41619873 \\\n",
    "              , 57.20663071, 59.99702072, 62.78735352, 65.5776062 \\\n",
    "              , 68.3677597 , 71.15775299, 73.9475174 , 76.73690033 \\\n",
    "              , 79.52560425, 82.31291199, 85.0965271 , 87.86380005 ]\n",
    "print can_use_sphere(lon, lat)[0]\n",
    "print can_use_sphere(lon, lat)[1]\n",
    "lat[-1] = 88.0\n",
    "print can_use_sphere(lon, lat)[0]\n",
    "print can_use_sphere(lon, lat)[1].splitlines()\n",
    "print can_use_sphere(lon, lat)[2]\n",
    "lat = [-90, -60, -30, 0, 30, 60, 90]\n",
    "print can_use_sphere(lon, lat)[0]\n",
    "print can_use_sphere(lon, lat)[1].splitlines()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def has_close(data, value, rtol=1.e-5, atol=1.e-8):\n",
    "    \"\"\"Test if data has any values \"equal\" to value.\n",
    "\n",
    "    Returns 1 if any of the elements of argument data has a value \n",
    "    \"equal\" to argument value; returns 0 otherwise.  If data or value\n",
    "    is floating point, \"equal\" means where abs(data-value) <= atol + \n",
    "    rtol * abs(value).  This is essentially the same algorithm used \n",
    "    in the Numeric function allclose.  If data and value are integer,\n",
    "    \"equal\" means strict equality.\n",
    "\n",
    "    Positional Input Arguments:\n",
    "    * data:   Data.  Scalar or Numeric array, Python list/tuple of\n",
    "              any size and shape.  Floating or integer type.\n",
    "    * value:  Test value.  Scalar or 1 element Numeric array, list,\n",
    "              or tuple.  Floating or integer type.\n",
    "\n",
    "    Keyword Input Arguments:\n",
    "    * rtol:   \"Relative\" tolerance.  Default is 1.e-5.  Used in the\n",
    "              comparison between data and value only if the two are \n",
    "              floating point.  Floating or integer type.\n",
    "    * atol:   \"Absolute\" tolerance.  Default is 1.e-8.  Used in the\n",
    "              comparison between data and value only if the two are \n",
    "              floating point.  Floating or integer type.\n",
    "\n",
    "    Examples:\n",
    "    >>> from has_close import has_close\n",
    "    >>> data = [20., -32., -1., 2., 5., 29.]\n",
    "    >>> has_close(data, -1.)\n",
    "    1\n",
    "    >>> has_close(data, 10.)\n",
    "    0\n",
    "    \"\"\"\n",
    "\n",
    "\n",
    "    #- Make sure data is Numeric type and value is a scalar within the\n",
    "    #  function:\n",
    "\n",
    "    dataN  = np.array(data)\n",
    "    valueS = np.array(value)\n",
    "\n",
    "\n",
    "    #- Safe compare if floating.  Strict compare if integer.  Any other\n",
    "    #  type returns an error:\n",
    "\n",
    "    if (dataN.dtype == float) or (valueS.dtype == type(1.)):\n",
    "        closemask = np.less_equal(np.abs(dataN-valueS), atol+rtol*np.abs(valueS))\n",
    "    elif (dataN.dtype == np.int32) and (valueS.dtype == type(1)):\n",
    "        closemask = np.where(dataN == valueS, 1, 0)\n",
    "    else:\n",
    "        print(dataN.dtype)\n",
    "        print(type(valueS))\n",
    "        raise ValueError, \"has_close:  Inputs must be float or integer\"\n",
    "\n",
    "\n",
    "    #- Return true if any elements of data has value:\n",
    "\n",
    "    if ma.maximum(closemask) == 1:\n",
    "        return 1\n",
    "    else:\n",
    "        return 0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = np.arange(20., dtype=np.float64)\n",
    "result = has_close(data, -1.)\n",
    "print result\n",
    "\n",
    "data = np.array([3, 2, 5, 10, -2, -3, 63])\n",
    "print has_close(data, -2)\n",
    "\n",
    "data = np.array([3, 2, 5, 10, -2, -3, 63])\n",
    "print has_close(data, -22)\n",
    "\n",
    "data = np.array([3, 2, 5, 10, -2, -3, 63])\n",
    "print has_close(data, -1.9999999)\n",
    "\n",
    "data = np.array([3, 2, 5, 10, -2, -3, 63])\n",
    "print has_close(data, -2.01)\n",
    "\n",
    "data = np.arange(20.)\n",
    "print has_close(data, 10)\n",
    "\n",
    "data = np.arange(20.)\n",
    "print has_close(data, [10])\n",
    "\n",
    "data = np.arange(20.)\n",
    "print has_close(data, np.array(15.))\n",
    "\n",
    "data = np.arange(20.)\n",
    "print has_close(data, np.array(50.))\n",
    "\n",
    "data = 4.\n",
    "print has_close(data, 2.)\n",
    "\n",
    "data = 4.\n",
    "print has_close(data, 4.)\n",
    "\n",
    "data = 4.\n",
    "print has_close(data, 3.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def deriv(*positional_inputs, **keyword_inputs):\n",
    "    \"\"\"Calculate the derivative along a single dimension.\n",
    "\n",
    "\n",
    "    Calling Sequence:\n",
    "        Result = deriv([x_in,] y_in, missing=1e+20, algorithm='default')\n",
    "\n",
    "\n",
    "    Positional Input Arguments:\n",
    "    * x_in:  Abscissa values of y_in to take with respect to.  If \n",
    "      not set, the derivative of y_in is take with respect to unit \n",
    "      abscissa intervals.  Numeric array of same shape and size as \n",
    "      y_in.  Must be monotonic and with no duplicate values.\n",
    "      Optional.  First positional argument out of two, if present.\n",
    "\n",
    "    * y_in:  Ordinate values, to take the derivative with respect \n",
    "      to.  Numeric array vector of rank 1.  Required.  Second posi-\n",
    "      tional argument, if x_in is present; only positional argument \n",
    "      if x_in is absent.\n",
    "\n",
    "\n",
    "    Keyword Input Arguments:\n",
    "    * missing:  If y_in and/or x_in has missing values, this is the \n",
    "      missing value value.  Scalar.  Default is 1e+20.\n",
    "\n",
    "    * algorithm:  Name of the algorithm to use.  String scalar.\n",
    "      Default is 'default'.  Possible values include:\n",
    "      + 'default':  Default method (currently set to 'order1').\n",
    "      + 'order1':  First-order finite-differencing (backward and\n",
    "        forward differencing used at the endpoints, and centered\n",
    "        differencing used everywhere else).  If abscissa intervals\n",
    "        are irregular, differencing will be correspondingly asym-\n",
    "        metric.\n",
    "\n",
    "\n",
    "    Output Result:\n",
    "    * Derivative of y_in with respect to x_in (or unit interval \n",
    "      abscissa, if x_in is not given).  Numeric array of same shape \n",
    "      and size as y_in.  If there are missing values, those elements \n",
    "      in the output are set to the value in |missing|.  For instance, \n",
    "      if y_in is only one element, a one-element vector is returned \n",
    "      as the derivative with the value of |missing|.  If there are \n",
    "      missing values in the output due to math errors and |missing| \n",
    "      is set to None, output will fill those missing values with the \n",
    "      MA default value of 1e+20.  \n",
    "\n",
    "\n",
    "    References:\n",
    "    * Press, W. H., et al. (1992):  Numerical Recipes in Fortran \n",
    "      77:  The Art of Scientific Computing.  New York, NY:  Cambridge\n",
    "      University Press, pp. 180-184.\n",
    "\n",
    "    * Wang, Y. (1999):  \"Numerical Differentiation,\" Introduction to \n",
    "      MHD Numerical Simulation in Space, ESS265: Instrumentation, \n",
    "      Data Processing and Data Analysis in Space Physics (UCLA).\n",
    "      URL:  http://www-ssc.igpp.ucla.edu/personnel/russell/ESS265/\n",
    "      Ch10/ylwang/node21.html.\n",
    "\n",
    "\n",
    "    Example with one argument, no missing values, using the default\n",
    "    method:\n",
    "    >>> from deriv import deriv\n",
    "    >>> import Numeric as N\n",
    "    >>> y = N.sin(N.arange(8))\n",
    "    >>> dydx = deriv(y)\n",
    "    >>> ['%.7g' % dydx[i] for i in range(4)]\n",
    "    ['0.841471', '0.4546487', '-0.3501755', '-0.83305']\n",
    "    >>> true = N.cos(N.arange(8))  #- Compare with exact solution\n",
    "    >>> ['%.7g' % true[i] for i in range(4)]  \n",
    "    ['1', '0.5403023', '-0.4161468', '-0.9899925']\n",
    "\n",
    "    Example with two arguments with missing values, using first-\n",
    "    order differencing:\n",
    "    >>> x = N.arange(8)/(2.*N.pi)\n",
    "    >>> y = N.sin(x)\n",
    "    >>> y[3] = 1e20            #- Set an element to missing value\n",
    "    >>> dydx = deriv(x, y, missing=1e20, algorithm='order1')\n",
    "    >>> ['%.7g' % dydx[i] for i in range(5)]\n",
    "    ['0.9957836', '0.9831985', '1e+20', '0.8844179', '1e+20']\n",
    "    >>> true = N.cos(x)       #- Compare with exact solution\n",
    "    >>> ['%.7g' % true[i] for i in range(5)]  \n",
    "    ['1', '0.9873616', '0.9497657', '0.8881628', '0.8041098']\n",
    "    \"\"\"\n",
    "\n",
    "    #- Establish y_in and x_in from *positional_inputs:\n",
    "\n",
    "    if len(positional_inputs) == 1:\n",
    "        y_in = positional_inputs[0]\n",
    "        x_in = np.arange(len(y_in), dtype=y_in.dtype)\n",
    "    elif len(positional_inputs) == 2:\n",
    "        x_in = positional_inputs[0]\n",
    "        y_in = positional_inputs[1]\n",
    "    else:\n",
    "        raise ValueError, \"deriv:  Bad inputs\"\n",
    "\n",
    "\n",
    "    #- Establish missing and algorithm from *keyword_inputs:\n",
    "\n",
    "    if keyword_inputs.has_key('missing') == 1:\n",
    "        missing = keyword_inputs['missing']\n",
    "    else:\n",
    "        missing = 1e+20\n",
    "\n",
    "    if keyword_inputs.has_key('algorithm') == 1:\n",
    "        algorithm = keyword_inputs['algorithm']\n",
    "    else:\n",
    "        algorithm = 'default'\n",
    "\n",
    "\n",
    "    #- Check positional and keyword inputs for possible errors:\n",
    "\n",
    "    if (len(y_in.shape) != 1) or (len(x_in.shape) != 1):\n",
    "        raise ValueError, \"deriv:  Inputs not a vector\"\n",
    "    if type(algorithm) != type(''):\n",
    "        raise ValueError, \"deriv:  algorithm not str\"\n",
    "\n",
    "\n",
    "    #- Set algorithm_to_use variable, based on the algorithm keyword.\n",
    "    #  The algorithm_to_use tells which algorithm below to actually\n",
    "    #  use (so here is where we set what algorithm to use for default):\n",
    "\n",
    "    if algorithm == 'default':\n",
    "        algorithm_to_use = 'order1'\n",
    "    else:\n",
    "        algorithm_to_use = algorithm\n",
    "\n",
    "\n",
    "    #- Change input to MA:  just set to input value unless there are\n",
    "    #  missing values, in which case add mask:\n",
    "\n",
    "    if missing == None:\n",
    "        x = ma.masked_array(x_in)\n",
    "        y = ma.masked_array(y_in)\n",
    "    else:\n",
    "        x = ma.masked_values(x_in, missing, copy=0)\n",
    "        y = ma.masked_values(y_in, missing, copy=0)\n",
    "\n",
    "\n",
    "    #- Calculate and return derivative:\n",
    "\n",
    "    #  * Create working arrays that are consistent with a 3-point\n",
    "    #    stencil in the interior and 2-point stencil on the ends:\n",
    "    #    *im1 means the point before index i, *ip1 means the point \n",
    "    #    after index i, and the i index array is just plain x or \n",
    "    #    y; the endpadded arrays replicate the ends of x and y.\n",
    "    #    I use an MA array filled approach instead of concatentation\n",
    "    #    because the MA concatenation routine doesn't work right\n",
    "    #    when the endpoint element is a missing value:\n",
    "\n",
    "    x_endpadded = ma.zeros(x.size+2, dtype=x.dtype)\n",
    "    x_endpadded[0]    = x[0]\n",
    "    x_endpadded[1:-1] = x \n",
    "    x_endpadded[-1]   = x[-1]\n",
    "\n",
    "    y_endpadded = ma.zeros(y.size+2, dtype=y.dtype)\n",
    "    y_endpadded[0]    = y[0]\n",
    "    y_endpadded[1:-1] = y\n",
    "    y_endpadded[-1]   = y[-1]\n",
    "\n",
    "    y_im1 = y_endpadded[:-2]\n",
    "    y_ip1 = y_endpadded[2:]\n",
    "    x_im1 = x_endpadded[:-2]\n",
    "    x_ip1 = x_endpadded[2:]\n",
    "\n",
    "\n",
    "    #  * Option 1:  First-order differencing (interior points use\n",
    "    #    centered differencing, and end points use forward or back-\n",
    "    #    ward differencing, as applicable):\n",
    "\n",
    "    if algorithm_to_use == 'order1':\n",
    "        dydx = (y_ip1 - y_im1) / (x_ip1 - x_im1) \n",
    "\n",
    "\n",
    "    #  * Option 2:  Bad algorithm specified:\n",
    "\n",
    "    else:\n",
    "        raise ValueError, \"deriv:  bad algorithm\"\n",
    "\n",
    "\n",
    "    #- Return derivative as Numeric array:\n",
    "\n",
    "    return ma.filled( dydx, missing )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y = np.array([3.,4.5])\n",
    "dydx = deriv(y)\n",
    "print dydx\n",
    "\n",
    "y = np.array([3.])\n",
    "dydx = deriv(y)\n",
    "print dydx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def curl_2d( x, y, Fx, Fy, missing=1e+20 \\\n",
    "           , algorithm='default', R_sphere=6.37122e+6):\n",
    "    \"\"\"Curl of a vector F on a 2-D \"rectangular\" grid.\n",
    "\n",
    "    The 2-D grid F is defined on is rectangular, meaning that while\n",
    "    the grid spacings do not have to be even, each grid box is\n",
    "    rectangular and grid boundaries are rectilinear and parallel to\n",
    "    each other.  If the 2-D grid is on a sphere, we assume that lines\n",
    "    of constant longitude and latitude form grids that are rectangular\n",
    "    (though in reality this is not strictly true).\n",
    "\n",
    "    Since F does not have a z-component, the curl of F is positive\n",
    "    pointing as a normal out of the plane defined by the basis for \n",
    "    x and y.\n",
    "\n",
    "\n",
    "    Positional Input Arguments:\n",
    "    * x:  x-coordinate of each F.  1-D Numeric array with number of \n",
    "      elements equal to the number of columns in array F.  Typically \n",
    "      x is in km, or on a sphere the longitude in deg.  Floating or\n",
    "      integer type.\n",
    "\n",
    "    * y:  y-coordinate of each F.  1-D Numeric array with number of \n",
    "      elements equal to the number of rows in array F.  Typically \n",
    "      y is in km, or on a sphere the latitude in deg.  Floating or\n",
    "      integer type.\n",
    "\n",
    "    * Fx:  x-component of vector quantity F.  2-D Numeric array of \n",
    "      shape (len(y), len(x)).  Floating or integer type.\n",
    "\n",
    "    * Fy:  y-component of vector quantity F.  2-D Numeric array of \n",
    "      same shape and size as Fx.  Floating or integer type.\n",
    "\n",
    "\n",
    "    Keyword Input Arguments:\n",
    "    * missing:  If Fx and/or Fy has missing values, this is the \n",
    "      missing value value.  Scalar.  Default is 1e+20.\n",
    "\n",
    "    * algorithm:  Name of the algorithm to use.  String scalar.\n",
    "      Default is 'default'.  Possible values include:\n",
    "\n",
    "      + 'default':  Default method.  Algorithm 'order1_cartesian' \n",
    "        is used.\n",
    "\n",
    "      + 'default_spherical':  If 'spherepack' can be used, chooses\n",
    "        that method.  Otherwise 'order1_spherical' is used.  Either\n",
    "        way this option assumes the domain is on a sphere of radius\n",
    "        equal to keyword R_sphere.\n",
    "\n",
    "      + 'order1_cartesian':  First-order finite-differencing (back-\n",
    "        ward and forward differencing used at the endpoints, and \n",
    "        centered differencing used everywhere else).  If abscissa \n",
    "        intervals are irregular, differencing will be correspon-\n",
    "        dingly asymmetric.  Grid is assumed Cartesian.  The curl\n",
    "        returned is in units of F divided by the units of x and y\n",
    "        (x and y must be in the same units).\n",
    "\n",
    "      + 'order1_spherical':  First-order finite-differencing (back-\n",
    "        ward and forward differencing used at the endpoints, and \n",
    "        centered differencing used everywhere else).  If abscissa \n",
    "        intervals are irregular, differencing will be correspon-\n",
    "        dingly asymmetric.  Grid is assumed spherical and x and y\n",
    "        are longitude and latitude (in deg).  The curl returned is \n",
    "        in units of F divided by units of R_sphere.  Note because\n",
    "        of singularities in the algorithm near the poles, the\n",
    "        curl for all points north of 88N or south of 88S latitude \n",
    "        are set to missing.\n",
    "\n",
    "      + 'spherepack':  The NCAR package SPHEREPACK 3.0 is used for \n",
    "        calculating the curl.  The spatial domain is the global \n",
    "        sphere; x and y are longitude and latitude (in deg), \n",
    "        respectively; x must be evenly spaced and y must be gauss-\n",
    "        ian or evenly spaced.  The domain must cover the entire \n",
    "        sphere and x and y must be longitude and latitude, respec-\n",
    "        tively, in deg.  There can be no missing values.  The \n",
    "        algorithm also assumes that x and y are monotonically in-\n",
    "        creasing from index 0.  The curl returned is in units of F \n",
    "        divided by m.  Thus, if F is velocity in m/s, the curl of \n",
    "        F is in units 1/s (and is the vorticity).  Note that in \n",
    "        this function's implementation, the final curl output using \n",
    "        algorithm 'spherepack' is adjusted for R_sphere not equal \n",
    "        to the default mean earth radius, if appropriate, so this \n",
    "        algorithm can be used on other planets.\n",
    "\n",
    "        Note that SPHEREPACK operates on single precision floating\n",
    "        point (Float32) values.  This function silently converts\n",
    "        input to Float32 for the computations, as needed (the input \n",
    "        parameters are not altered, however).\n",
    "\n",
    "    * R_sphere:  If the grid is the surface of a sphere, this is \n",
    "      the radius of the sphere.  This keyword is only used when \n",
    "      the algorithm is 'order1_spherical' or 'spherepack' (either \n",
    "      chosen explicitly by the algorithm keyword or automatically \n",
    "      when 'default_spherical' is used).  Default value is the \n",
    "      \"mean\" radius of the Earth, in meters.  \"Mean\" is in quotes\n",
    "      because this value changes depending on what you consider the\n",
    "      mean; the default in this function is the same Earth radius\n",
    "      used in module sphere.  Note that if the level is a distance \n",
    "      z above the \"mean\" radius of the earth, R_sphere should be \n",
    "      set to the \"mean\" radius plus z.  \n",
    "      \n",
    "      Value can be a scalar or a Numeric or MA array of the same \n",
    "      size and shape as Fx.  If keyword is an array, the value of\n",
    "      R_sphere for each element corresponds to the same location\n",
    "      as in the Fx and Fy arrays; there should also be no missing \n",
    "      values in R_sphere (function doesn't check for this, however).\n",
    "\n",
    "\n",
    "    Output:\n",
    "    * Curl of F.  Units depend on algorithm used to calculate curl\n",
    "      (see above discussion of keyword algorithm).  Numeric array of \n",
    "      same shape as Fx and Fy inputs.  \n",
    "      \n",
    "      If there are missing values, those elements in the output that \n",
    "      used missing values in the calculation of the curl are set to \n",
    "      the value in keyword missing.  If there are missing values in \n",
    "      the output due to math errors and keyword missing is set to \n",
    "      None, output will fill those missing values with the MA default \n",
    "      value of 1e+20.\n",
    "\n",
    "\n",
    "    References:\n",
    "    * Glickman, T. S. (Ed.) (2000):  \"Curl,\" Glossary of Meteorology,\n",
    "      Boston, MA:  American Meteorological Society, ISBN 1-878220-\n",
    "      34-9.\n",
    "    * Haltiner, G. J, and R. T. Williams (1980):  Numerical Prediction\n",
    "      and Dynamic Meteorology, New York:  John Wiley & Sons, pp. 8-9.\n",
    "    * Holton, J. R. (1992):  An Introduction to Dynamic Meteorology,\n",
    "      San Diego, CA:  Academic Press, p. 482.\n",
    "    * Kauffman, B. G., and W. G. Large (2002):  The CCSM Coupler,\n",
    "      Version 5.0, User's Guide, Source Code Reference, and Scientific\n",
    "      Description.  Boulder, CO:  NCAR.  Value for earth radius is \n",
    "      given at URL:\n",
    "      http://www.ccsm.ucar.edu/models/ccsm2.0/cpl5/users_guide/node10.html\n",
    "    * Overview of the CDAT Interface to the NCAR SPHEREPACK 3.0.  \n",
    "      SPHEREPACK is by J. C. Adams and P. N. Swarztrauber.  More\n",
    "      information on SPHEREPACK can be found at:  \n",
    "         http://www.scd.ucar.edu/css/software/spherepack/\n",
    "      while information on the CDAT implementation used here is at:\n",
    "         http://esg.llnl.gov/cdat/getting_started/cdat.htm\n",
    "\n",
    "\n",
    "    Example of a synthetic dataset of winds.  The curl is the vorti-\n",
    "    city.  Images of the wind field and the vorticity are online at\n",
    "    http://www.johnny-lin.com/py_pkgs/gemath/doc/test_curl_2d_add.html.\n",
    "    Velocity at y greater than 60 and less than -60 are set to 0, to\n",
    "    prevent spectral algorithms from giving errors at the poles.  Tiny \n",
    "    values less than 1e-10 are also set to 0:\n",
    "\n",
    "    (a) Import statements and create data:\n",
    "\n",
    "    >>> import Numeric as N\n",
    "    >>> from curl_2d import curl_2d\n",
    "    >>> nx = 181\n",
    "    >>> ny = 91\n",
    "    >>> x = N.arange(nx) * 2.0 - 180.0\n",
    "    >>> y = N.arange(ny) * 2.0 - 90.0\n",
    "    >>> y_2d = N.reshape( N.repeat(y, nx), (ny, nx) )\n",
    "    >>> x_2d = N.reshape( N.repeat(N.reshape(x,(1,nx)), ny), (ny, nx) )\n",
    "    >>> u = N.sin(x_2d*N.pi/180.)*N.sin((y_2d-90.)*N.pi/30.)\n",
    "    >>> v = N.sin((x_2d-30.)*N.pi/180.)*N.sin((y_2d-90.)*N.pi/30.)\n",
    "    >>> N.putmask( u, N.where(N.absolute(u) < 1e-10, 1, 0), 0. )\n",
    "    >>> N.putmask( v, N.where(N.absolute(v) < 1e-10, 1, 0), 0. )\n",
    "    >>> N.putmask( u, N.where(N.absolute(y_2d) > 60., 1, 0), 0. )\n",
    "    >>> N.putmask( v, N.where(N.absolute(y_2d) > 60., 1, 0), 0. )\n",
    "\n",
    "    (b) Use order 1 Cartesian algorithm:  If x and y are in meters,\n",
    "        and u and v are in m/s, the curl is in 1/s:\n",
    "\n",
    "    >>> curl = curl_2d(x, y, u, v, algorithm='order1_cartesian')\n",
    "    >>> ['%.7g' % curl[45,i] for i in range(88,92)]\n",
    "    ['-0.007251593', '-0.003628007', '0', '0.003628007']\n",
    "    >>> ['%.7g' % curl[0,i]  for i in range(8,12)]\n",
    "    ['0', '0', '0', '0']\n",
    "\n",
    "    (c) Use spectral method from SPHEREPACK (need to first remove the\n",
    "        repeating dateline point from the domain).  Here we take x and\n",
    "        y in deg, while u and v are in m/s.  If so the curl is in 1/s.\n",
    "        These results as much less than for example (b) above since\n",
    "        the domain is much larger:\n",
    "\n",
    "    >>> x = x[0:-1]\n",
    "    >>> u = u[:,0:-1]\n",
    "    >>> v = v[:,0:-1]\n",
    "    >>> curl = curl_2d(x, y, u, v, algorithm='spherepack')\n",
    "    >>> ['%.7g' % curl[45,i] for i in range(88,92)]\n",
    "    ['-6.436402e-08', '-3.220152e-08', '1.33852e-13', '3.220165e-08']\n",
    "    >>> ['%.7g' % curl[0,i]  for i in range(8,12)]\n",
    "    ['-2.299736e-20', '-4.806512e-21', '-9.283145e-22', '-1.368445e-20']\n",
    "\n",
    "    (d) Use order 1 spherical algorithm:  x and y are in deg and u and\n",
    "        v are in m/s.  The curl is in 1/s.  Note that these results are\n",
    "        nearly identical to those from SPHEREPACK, except near the poles\n",
    "        the values are set to missing:\n",
    "\n",
    "    >>> curl = curl_2d(x, y, u, v, algorithm='order1_spherical')\n",
    "    >>> ['%.7g' % curl[45,i] for i in range(88,92)]\n",
    "    ['-6.517317e-08', '-3.260645e-08', '0', '3.260645e-08']\n",
    "    >>> ['%.7g' % curl[0,i]  for i in range(8,12)]\n",
    "    ['1e+20', '1e+20', '1e+20', '1e+20']\n",
    "\n",
    "    (e) Use \"default\" algorithm, which for this domain chooses the \n",
    "        order 1 cartesian algorithm:\n",
    "\n",
    "    >>> curl = curl_2d(x, y, u, v, algorithm='default')\n",
    "    >>> ['%.7g' % curl[45,i] for i in range(88,92)]\n",
    "    ['-0.007251593', '-0.003628007', '0', '0.003628007']\n",
    "\n",
    "    (f) Use \"default_spherical\" algorithm with missing data (equal to\n",
    "        1e+20), which because of the missing data will choose the order\n",
    "        1 spherical algorithm.  We also do the calculation on Mars, and\n",
    "        illustrate passing in the radius as a constant as well as an\n",
    "        array of the same size and shape as u and v:\n",
    "\n",
    "    >>> u[30:46,90:114] = 1e+20\n",
    "    >>> curl = curl_2d( x, y, u, v, algorithm='default_spherical' \\\n",
    "                      , R_sphere=3390e+3)\n",
    "    >>> ['%.7g' % curl[45,i] for i in range(88,92)]\n",
    "    ['-1.224875e-07', '-6.128107e-08', '1e+20', '1e+20']\n",
    "\n",
    "    >>> R_mars = N.zeros(u.shape, typecode=N.Float) + 3390e+3\n",
    "    >>> curl = curl_2d( x, y, u, v, algorithm='default_spherical' \\\n",
    "                      , R_sphere=R_mars )\n",
    "    >>> ['%.7g' % curl[45,i] for i in range(88,92)]\n",
    "    ['-1.224875e-07', '-6.128107e-08', '1e+20', '1e+20']\n",
    "    \"\"\"\n",
    "    #import MA\n",
    "    #import Numeric as N\n",
    "    #from has_close import has_close\n",
    "    #from can_use_sphere import can_use_sphere\n",
    "\n",
    "\n",
    "\n",
    "#------- Nested Function:  Curl By Order 1 Cartesian Finite Diff. ------\n",
    "\n",
    "    def _order1_cartesian_curl():\n",
    "        \"\"\"\n",
    "        Calculate curl using Cartesian first-order differencing.\n",
    "\n",
    "        An ma array is returned.\n",
    "\n",
    "        Algorithm:  First-order differencing (interior points use \n",
    "        centered differencing, and end points use forward or backward \n",
    "        differencing, as applicable) in Cartesian coordinates (see \n",
    "        Glickman [2000], p. 194).\n",
    "        \"\"\"\n",
    "        dFy_dx_N = np.zeros( (len(y), len(x)), dtype=float )\n",
    "        dFx_dy_N = np.zeros( (len(y), len(x)), dtype=float )\n",
    "\n",
    "        for iy in range(len(y)):\n",
    "            dFy_dx_N[iy,:] = deriv( x, np.ravel(Fy[iy,:]) \\\n",
    "                                  , missing=missing, algorithm='order1')\n",
    "\n",
    "        for ix in range(len(x)):\n",
    "            dFx_dy_N[:,ix] = deriv( y, np.ravel(Fx[:,ix]) \\\n",
    "                                  , missing=missing, algorithm='order1')\n",
    "\n",
    "        dFy_dx = ma.masked_values(dFy_dx_N, missing, copy=0)\n",
    "        dFx_dy = ma.masked_values(dFx_dy_N, missing, copy=0)\n",
    "\n",
    "        return dFy_dx - dFx_dy\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#------- Nested Function:  Curl By Order 1 Spherical Finite Diff. ------\n",
    "\n",
    "    def _order1_spherical_curl():\n",
    "        \"\"\"\n",
    "        Calculate curl using spherical first-order differencing.\n",
    "\n",
    "        An ma array is returned.\n",
    "\n",
    "        Algorithm:  First-order differencing (interior points use \n",
    "        centered differencing, and end points use forward or backward \n",
    "        differencing, as applicable) in spherical coordinates (see \n",
    "        Holton 1992).\n",
    "\n",
    "        Note because of singularities in the algorithm near the poles, \n",
    "        the return value for all points north of 88N or south of 88S \n",
    "        latitude are set to mask true.\n",
    "\n",
    "        Key to some variables:\n",
    "        * xrad is x in radians (x is assumed in deg)\n",
    "        * yrad is y in radians (y is assumed in deg)\n",
    "        * ddFy is d(Fy) / d(xrad) as a masked array.  The version with\n",
    "          a \"_N\" suffix means the Numeric version.\n",
    "        * ddFx is d(Fx * cos(yrad)) / d(yrad) as a masked array.  The \n",
    "          version with a \"_N\" suffix means the Numeric version.\n",
    "        \"\"\"\n",
    "\n",
    "        xrad = x * np.pi/180.0\n",
    "        yrad = y * np.pi/180.0\n",
    "\n",
    "\n",
    "        #- Derivative preliminaries:\n",
    "\n",
    "        ddFy_N = np.zeros( (len(yrad), len(xrad)), dtype=float )\n",
    "        ddFx_N = np.zeros( (len(yrad), len(xrad)), dtype=float )\n",
    "\n",
    "        for iy in range(len(yrad)):\n",
    "            ddFy_N[iy,:] = deriv( xrad, np.ravel(Fy[iy,:]) \\\n",
    "                                , missing=missing, algorithm='order1')\n",
    "\n",
    "        for ix in range(len(xrad)):\n",
    "            tmp_MA = ma.masked_values(np.ravel(Fx[:,ix]), missing, copy=0)\n",
    "            tmp_N  = ma.filled(tmp_MA * np.cos(yrad), missing)\n",
    "            ddFx_N[:,ix] = deriv( yrad, tmp_N \\\n",
    "                                , missing=missing, algorithm='order1')\n",
    "\n",
    "        ddFy = ma.masked_values(ddFy_N, missing, copy=0)\n",
    "        ddFx = ma.masked_values(ddFx_N, missing, copy=0)\n",
    "\n",
    "\n",
    "        #- Calculate the curl:\n",
    "\n",
    "        tmpr = ( ddFy - ddFx ) / \\\n",
    "               ( R_sphere * np.reshape( np.repeat(np.cos(yrad),len(xrad)) \\\n",
    "                                     , (len(yrad),len(xrad)) ) )\n",
    "\n",
    "\n",
    "        #- Make points \"near\" poles be missing and return:\n",
    "\n",
    "        np_mask_lat =  88.0\n",
    "        sp_mask_lat = -88.0\n",
    "        yarr = np.reshape( np.repeat(y,len(x)), (len(y),len(x)) )\n",
    "\n",
    "        np_mask = ma.make_mask( np.where(yarr > np_mask_lat, 1, 0) )\n",
    "        sp_mask = ma.make_mask( np.where(yarr < sp_mask_lat, 1, 0) )\n",
    "\n",
    "        tmpr_mask = tmpr.mask\n",
    "        tmpr_mask = ma.mask_or(tmpr_mask, np_mask)\n",
    "        tmpr_mask = ma.mask_or(tmpr_mask, sp_mask)\n",
    "\n",
    "        return ma.masked_array(tmpr, mask=tmpr_mask)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#--------- Nested Function:  Curl By SPHEREPACK Spectral Method --------\n",
    "\n",
    "    def _spherepack_curl():\n",
    "        \"\"\"\n",
    "        Calculate curl using SPHEREPACK spectral methods.\n",
    "        \n",
    "        If the value of R_sphere is not the same as the earth radius \n",
    "        used in package sphere, adjust curl to reflect R_sphere.  A \n",
    "        Numeric array is returned.\n",
    "\n",
    "        SPHEREPACK outputs a number of diagnostic messages to\n",
    "        stdout (and perhaps to stderr).  The one that's irritating\n",
    "        is the one saying data will be converted to Float32.  Thus,\n",
    "        we silently ensure that calculations are used only on Float32 \n",
    "        data, to surpress the warnings.  I did this because redirect-\n",
    "        ing stdout/stderr turned out to be too difficult; just\n",
    "        altering the sys attributes didn't work.\n",
    "        \"\"\"\n",
    "        import sphere\n",
    "\n",
    "        sph_obj = sphere.Sphere( x.astype(ma.Float32) \\\n",
    "                               , y.astype(ma.Float32) )\n",
    "        if np.allclose(np.array(R_sphere), sphere.radius):\n",
    "            curl = sph_obj.vrt( Fx.astype(ma.Float32) \\\n",
    "                              , Fy.astype(ma.Float32) )\n",
    "        else:\n",
    "            curl = sph_obj.vrt( Fx.astype(ma.Float32) \\\n",
    "                              , Fy.astype(ma.Float32) ) \\\n",
    "                 * sphere.radius / R_sphere\n",
    "\n",
    "        return curl\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#-------------- Overall Function:  Settings, Body, Return --------------\n",
    "\n",
    "    #- Choose algorithm to compute curl:\n",
    "    \n",
    "    if algorithm == 'default':\n",
    "        _calculate_curl = _order1_cartesian_curl\n",
    "\n",
    "    elif algorithm == 'default_spherical':\n",
    "        if (can_use_sphere(x,y)[0] == 1) and \\\n",
    "           (not has_close(Fx, missing)) and \\\n",
    "           (not has_close(Fy, missing)):\n",
    "            _calculate_curl = _spherepack_curl\n",
    "        else:\n",
    "            _calculate_curl = _order1_spherical_curl\n",
    "\n",
    "    elif algorithm == 'order1_cartesian':\n",
    "        _calculate_curl = _order1_cartesian_curl\n",
    "\n",
    "    elif algorithm == 'order1_spherical':\n",
    "        _calculate_curl = _order1_spherical_curl\n",
    "\n",
    "    elif algorithm == 'spherepack':\n",
    "        if has_close(Fx, missing) or has_close(Fy, missing):\n",
    "            raise ValueError, \"curl_2d:  has missing values\"\n",
    "        else:\n",
    "            _calculate_curl = _spherepack_curl\n",
    "\n",
    "    else:\n",
    "        raise ValueError, \"curl_2d:  bad algorithm\"\n",
    "\n",
    "\n",
    "    #- Calculate curl and return from function:\n",
    "\n",
    "    return ma.filled( _calculate_curl(), missing )\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nx = 181\n",
    "ny = 91\n",
    "\n",
    "x = np.arange(nx) * 2.0 - 180.\n",
    "y = np.arange(ny) * 2.0 - 90.\n",
    "y_2d = np.reshape(np.repeat(y, nx), (ny, nx))\n",
    "x_2d = np.reshape(np.repeat(np.reshape(x, (1, nx)), ny), (ny, nx))\n",
    "\n",
    "xx, yy = np.meshgrid(x, y)\n",
    "u = np.sin(xx*np.pi/180.)*np.sin((yy-90.)*np.pi/30.)\n",
    "v = np.sin((xx-30.)*np.pi/180.)*np.sin((yy-90.)*np.pi/30.)\n",
    "np.putmask( u, np.where(np.absolute(u) < 1e-10, 1, 0), 0. )\n",
    "np.putmask( v, np.where(np.absolute(v) < 1e-10, 1, 0), 0. )\n",
    "np.putmask( u, np.where(np.absolute(yy) > 60., 1, 0), 0. )\n",
    "np.putmask( v, np.where(np.absolute(yy) > 60., 1, 0), 0. )\n",
    "#x = x[0:-1]\n",
    "#u = u[:,0:-1]\n",
    "#v = v[:,0:-1]\n",
    "\n",
    "R_earth = ma.zeros(u.shape, dtype=float) + 6371220.0\n",
    "R_earth[45,89] = 6371220.0 + 2000.0\n",
    "curl = curl_2d(x, y, u, v, algorithm='order1_spherical', R_sphere=R_earth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "['%.7g' % curl[45,i] for i in range(88,92)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "['%.7g' % curl[0,i]  for i in range(8,12)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.contourf(x[1:-1], y[1:-1], curl[1:-1,1:-1])\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "xx, yy = np.meshgrid(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(x)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(xx)\n",
    "print(yy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(x_2d)\n",
    "print(y_2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(u.shape)\n",
    "print(xx.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def mg_rel(u, v, lon, lat):\n",
    "    dx = np.zeros((len(lat), len(lon)-2))\n",
    "    dy = np.zeros((len(lat)-2, len(lon)))\n",
    "    du = np.zeros((len(lat)-2, len(lon)))\n",
    "    dv = np.zeros((len(lat), len(lon)-2))\n",
    "    zeta = np.zeros((len(lat)-2, len(lon)-2))\n",
    "\n",
    "    for i in xrange(1, len(lon)-1):\n",
    "        for j in xrange(0, len(lat)):\n",
    "            dx[j, i-1] = metutils.convert((lon[i+1] - lon[i-1]) * \\\n",
    "                                          np.cos(np.pi*lat[j]/180.),\n",
    "                                          \"deg\", \"m\")\n",
    "            dv[j, i-1] = v[j, i+1] - v[j, i-1]\n",
    "\n",
    "    for i in xrange(0, len(lon)):\n",
    "        for j in xrange(1, len(lat)-1):\n",
    "            dy[j-1, i] = metutils.convert((lat[j+1] - lat[j-1]), \"deg\", \"m\")\n",
    "            du[j-1, i] = u[j+1, i] - u[j-1, i]\n",
    "\n",
    "    for i in xrange(len(lon) - 2):\n",
    "        for j in xrange(len(lat) - 2):\n",
    "            zeta[j, i] = dv[j, i]/dx[j, i] - du[j, i]/dy[j, i]\n",
    "    return zeta\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "zeta = mg_rel(u,v,x,y)\n",
    "plt.contourf(x[1:-1], y[1:-1], zeta)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "['%.7g' % zeta[45,i]  for i in range(88,92)]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}