jswhit/reading_netCDF.ipynb

## reading_netCDF.ipynb
{
 "metadata": {
  "celltoolbar": "Slideshow",
  "name": "",
  "signature": "sha256:4edf2a4e7e7449ca4ecbb1946583e5a8fb9c3634af1b7bbf57fce7bee30a6d99"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "# Reading netCDF data\n",
      "- requires [numpy](http://numpy.scipy.org) and netCDF/HDF5 C libraries.\n",
      "- Github site: https://github.com/Unidata/netcdf4-python\n",
      "- Online docs: http://unidata.github.io/netcdf4-python/\n",
      "- Based on Konrad Hinsen's old [Scientific.IO.NetCDF](http://dirac.cnrs-orleans.fr/plone/software/scientificpython/) API, with lots of added netcdf version 4 features.\n",
      "- Developed by Jeff Whitaker at NOAA, with many contributions from users."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Interactively exploring a netCDF File\n",
      "\n",
      "Let's explore a netCDF file from the *Atlantic Real-Time Ocean Forecast System*\n",
      "\n",
      "first, import netcdf4-python"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import netCDF4"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Create a netCDF4.Dataset object\n",
      "- **`f`** is an object, representing an open netCDF file.\n",
      "- printing the object gives you summary information, similar to *`ncdump -h`*.\n",
      "- the **`variables`** attribute of **`f`** is a dictionary containing netcdf variable objects."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')\n",
      "print f "
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<type 'netCDF4.Dataset'>\n",
        "root group (NETCDF4_CLASSIC data model, file format UNDEFINED):\n",
        "    Conventions: CF-1.0\n",
        "    title: HYCOM ATLb2.00\n",
        "    institution: National Centers for Environmental Prediction\n",
        "    source: HYCOM archive file\n",
        "    experiment: 90.9\n",
        "    history: archv2ncdf3z\n",
        "    dimensions(sizes): MT(1), Y(850), X(712), Depth(10)\n",
        "    variables(dimensions): float64 \u001b[4mMT\u001b[0m(MT), float64 \u001b[4mDate\u001b[0m(MT), float32 \u001b[4mDepth\u001b[0m(Depth), int32 \u001b[4mY\u001b[0m(Y), int32 \u001b[4mX\u001b[0m(X), float32 \u001b[4mLatitude\u001b[0m(Y,X), float32 \u001b[4mLongitude\u001b[0m(Y,X), float32 \u001b[4mu\u001b[0m(MT,Depth,Y,X), float32 \u001b[4mv\u001b[0m(MT,Depth,Y,X), float32 \u001b[4mtemperature\u001b[0m(MT,Depth,Y,X), float32 \u001b[4msalinity\u001b[0m(MT,Depth,Y,X)\n",
        "    groups: \n",
        "\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Access a netCDF variable\n",
      "- variable objects stored by name in **`variables`** dict.\n",
      "- print the variable yields summary info (including all the attributes).\n",
      "- no actual data read yet (just have a reference to the variable object with metadata)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "temp = f.variables['temperature']  # temperature variable\n",
      "print temp "
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<type 'netCDF4.Variable'>\n",
        "float32 temperature(MT, Depth, Y, X)\n",
        "    coordinates: Longitude Latitude Date\n",
        "    standard_name: sea_water_potential_temperature\n",
        "    units: degC\n",
        "    _FillValue: 1.26765e+30\n",
        "    valid_range: [ -5.07860279  11.14989948]\n",
        "    long_name:   temp [90.9H]\n",
        "unlimited dimensions: MT\n",
        "current shape = (1, 10, 850, 712)\n",
        "filling on\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## List the Dimensions\n",
      "\n",
      "- All variables in a netCDF file have an associated shape, specified by a list of dimensions.\n",
      "- Let's list all the dimensions in this netCDF file.\n",
      "- Note that the **`MT`** dimension is special (*`unlimited`*), which means it can be appended to."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for dname, d in f.dimensions.items():\n",
      "    print(d)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<type 'netCDF4.Dimension'> (unlimited): name = 'MT', size = 1\n",
        "\n",
        "<type 'netCDF4.Dimension'>: name = 'Y', size = 850\n",
        "\n",
        "<type 'netCDF4.Dimension'>: name = 'X', size = 712\n",
        "\n",
        "<type 'netCDF4.Dimension'>: name = 'Depth', size = 10\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Each variable has a **`dimensions`** and a **`shape`** attribute."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "temp.dimensions"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "(u'MT', u'Depth', u'Y', u'X')"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "temp.shape"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "(1, 10, 850, 712)"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "### Each dimension typically has a variable associated with it (called a *coordinate* variable).\n",
      "- *Coordinate variables* are 1D variables that have the same name as dimensions.\n",
      "- Coordinate variables and *auxiliary coordinate variables* (named by the *coordinates* attribute) locate values in time and space."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mt = f.variables['MT']\n",
      "depth = f.variables['Depth']\n",
      "x,y = f.variables['X'], f.variables['Y']\n",
      "print mt \n",
      "print x                 "
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<type 'netCDF4.Variable'>\n",
        "float64 MT(MT)\n",
        "    long_name: time\n",
        "    units: days since 1900-12-31 00:00:00\n",
        "    calendar: standard\n",
        "    axis: T\n",
        "unlimited dimensions: MT\n",
        "current shape = (1,)\n",
        "filling on, default _FillValue of 9.96920996839e+36 used\n",
        "\n",
        "<type 'netCDF4.Variable'>\n",
        "int32 X(X)\n",
        "    point_spacing: even\n",
        "    axis: X\n",
        "unlimited dimensions: \n",
        "current shape = (712,)\n",
        "filling on, default _FillValue of -2147483647 used\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Accessing data from a netCDF variable object\n",
      "\n",
      "- netCDF variables objects behave much like numpy arrays.\n",
      "- slicing a netCDF variable object returns a numpy array with the data.\n",
      "- Boolean array and integer sequence indexing behaves differently for netCDF variables than for numpy arrays. Only 1-d boolean arrays and integer sequences are allowed, and these indices work independently along each dimension (similar to the way vector subscripts work in fortran)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "time = mt[:]  # Reads the netCDF variable MT, array of one element\n",
      "print time "
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 41023.25]\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dpth = depth[:] # examine depth array\n",
      "print dpth "
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[    0.   100.   200.   400.   700.  1000.  2000.  3000.  4000.  5000.]\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xx,yy = x[:],y[:]\n",
      "print 'shape of temp variable:',temp.shape\n",
      "tempslice = temp[0, dpth > 400, xx > xx.max()/2, yy < yy.max()/2]\n",
      "print 'shape of temp slice:',tempslice.shape"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "shape of temp variable: (1, 10, 850, 712)\n",
        "shape of temp slice:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " (6, 356, 424)\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## What is the sea surface temperature and salinity at 50N, 140W?\n",
      "### Finding the latitude and longitude indices of 50N, 140W\n",
      "\n",
      "- The `X` and `Y` dimensions don't look like longitudes and latitudes\n",
      "- Use the auxilary coordinate variables named in the `coordinates` variable attribute, `Latitude` and `Longitude`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "lat, lon = f.variables['Latitude'], f.variables['Longitude']\n",
      "print lat"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<type 'netCDF4.Variable'>\n",
        "float32 Latitude(Y, X)\n",
        "    standard_name: latitude\n",
        "    units: degrees_north\n",
        "unlimited dimensions: \n",
        "current shape = (850, 712)\n",
        "filling on, default _FillValue of 9.96920996839e+36 used\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Aha!  So we need to find array indices `iy` and `ix` such that `Latitude[iy, ix]` is close to 50.0 and `Longitude[iy, ix]` is close to -140.0 ..."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "latvals = lat[:]; lonvals = lon[:] # extract lat/lon values (in degrees) to numpy arrays\n",
      "import numpy as np\n",
      "# a function to find the index of the point closest (in squared distance) to give lat/lon value.\n",
      "def getclosest_ij(lats,lons,latpt,lonpt):\n",
      "    dist_sq = (lats-latpt)**2 + (lons-lonpt)**2  # find squared distance of every point on grid\n",
      "    minindex_flattened = dist_sq.argmin() # 1D index of minimum dist_sq element\n",
      "    # Get 2D index for latvals and lonvals arrays from 1D index\n",
      "    return np.unravel_index(minindex_flattened, latvals.shape)\n",
      "iy_min, ix_min = getclosest_ij(latvals, lonvals, 50., -140)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "### Now we have all the information we need to find our answer.\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "```\n",
      "|----------+--------|\n",
      "| Variable |  Index |\n",
      "|----------+--------|\n",
      "| MT       |      0 |\n",
      "| Depth    |      0 |\n",
      "| Y        | iy_min |\n",
      "| X        | ix_min |\n",
      "|----------+--------|\n",
      "```"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "### What is the sea surface temperature and salinity at the specified point?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sal = f.variables['salinity']\n",
      "# Read values out of the netCDF file for temperature and salinity\n",
      "print temp[0,0,iy_min,ix_min], temp.units\n",
      "print sal[0,0,iy_min,ix_min], sal.units"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "6.46312 degC\n",
        "32.6572 psu\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Remote data access via openDAP\n",
      "\n",
      "- Remote data can be accessed seamlessly with the netcdf4-python API\n",
      "- Access happens via the DAP protocol and DAP servers, such as TDS.\n",
      "- many formats supported, like GRIB, are supported \"under the hood\"."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "The following example showcases some nice netCDF features:\n",
      "\n",
      "1. We are seamlessly accessing **remote** data, from a TDS server.\n",
      "2. We are seamlessly accessing **GRIB2** data, as if it were netCDF data.\n",
      "3. We are generating **metadata** on-the-fly.\n",
      "4. We are using [pyUDL](https://github.com/Unidata/pyudl) to get URL from **TDS data catalog**."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# import pyUDL\n",
      "from pyudl.tds import get_latest_dods_url\n",
      "# specify URL for TDS data catalog\n",
      "gfs_catalog_url = \"http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/Global_0p5deg/catalog.xml\" \n",
      "# get latest forecast time available from catalog\n",
      "latest = get_latest_dods_url(gfs_catalog_url)\n",
      "# open the remote dataset\n",
      "gfs = netCDF4.Dataset(latest)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Look at metadata for a specific variable\n",
      "sfctmp = gfs.variables['Temperature_surface']\n",
      "print sfctmp\n",
      "for dname in sfctmp.dimensions:\n",
      "    print gfs.variables[dname]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<type 'netCDF4.Variable'>\n",
        "float32 Temperature_surface(time, lat, lon)\n",
        "    long_name: Temperature @ Ground or water surface\n",
        "    units: K\n",
        "    missing_value: nan\n",
        "    abbreviation: TMP\n",
        "    coordinates: time \n",
        "    Grib_Variable_Id: VAR_0-0-0_L1\n",
        "    Grib2_Parameter: [0 0 0]\n",
        "    Grib2_Parameter_Discipline: Meteorological products\n",
        "    Grib2_Parameter_Category: Temperature\n",
        "    Grib2_Parameter_Name: Temperature\n",
        "    Grib2_Level_Type: Ground or water surface\n",
        "    Grib2_Generating_Process_Type: Forecast\n",
        "unlimited dimensions: \n",
        "current shape = (32, 361, 720)\n",
        "filling off\n",
        "\n",
        "<type 'netCDF4.Variable'>\n",
        "float64 time(time)\n",
        "    units: Hour since 2014-10-18T12:00:00Z\n",
        "    standard_name: time\n",
        "    long_name: GRIB forecast or observation time\n",
        "    _CoordinateAxisType: Time\n",
        "unlimited dimensions: \n",
        "current shape = (32,)\n",
        "filling off\n",
        "\n",
        "<type 'netCDF4.Variable'>\n",
        "float32 lat(lat)\n",
        "    units: degrees_north\n",
        "    _CoordinateAxisType: Lat\n",
        "unlimited dimensions: \n",
        "current shape = (361,)\n",
        "filling off\n",
        "\n",
        "<type 'netCDF4.Variable'>\n",
        "float32 lon(lon)\n",
        "    units: degrees_east\n",
        "    _CoordinateAxisType: Lon\n",
        "unlimited dimensions: \n",
        "current shape = (720,)\n",
        "filling off\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "##Missing values\n",
      "- when `data == var.missing_value` somewhere, a masked array is returned.\n",
      "- illustrate with soil moisture data (only defined over land)\n",
      "- white areas on plot are masked values over water."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "soilmvar = gfs.variables['Volumetric_Soil_Moisture_Content_depth_below_surface_layer']\n",
      "soilm = soilmvar[0,0,::-1,:] # flip the data in latitude so North Hemisphere is up on the plot\n",
      "print soilm.shape, type(soilm), soilmvar.missing_value\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline\n",
      "cs = plt.contourf(soilm)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(361, 720) <class 'numpy.ma.core.MaskedArray'> nan\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAVMAAAD7CAYAAADXc3dDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnW2MJVeZ3/9/23Tztm1vGxi/zTKWwIkHMdhs7JDAhiYh\nxpYSzIcVZkYbkYCSlWAAGbyZGa/WNIt4aWKT1e6IVbIY5HXoDl4DDhYL2CZuxH4wL1nbbTw4YC0j\nMY7ddtwwA0E747GffLhVd84991TVqapTdU9VPT+p1ffWrZdzb53zr+ec5znPoYhAURRFqcdpsy6A\noihKH1AxVRRFCYCKqaIoSgBUTBVFUQKgYqooihIAFVNFUZQAnDGLi5LUeCxFUTqHiDDrs1zLlOTz\nSX6X5AMkD5H8RLJ9meQRkvcnf1cZxxwg+ROSj5C8IqdQ0f99+MMfnnkZtJxaTi1nHOUsItcyFZG/\nJ/kmEfk1yTMA/A3JNwAQAJ8WkU9b4rsTwDUAdgI4H8A9JC8SkecKS6IoitJhCsdMReTXycs5AKcD\n+Hny3mXuXg1gTUSeEZHDAB4FcHmAciqKokRNoZiSPI3kAwA2AdwrIg8nH72P5IMkbyZ5VrLtPABH\njMOPYGShdpKlpaVZF8ELLWdYtJxhGUo56TMWAAAkzwTwTQD7ARwC8FTy0UcBnCsi7yb5ZwDuE5Ev\nJMd8FsBfi8iXrXOJ73UVRVFigCQkxwHl7c0XkaMkvwbgH4nIunGBzwK4M3n7GIDtxmEXJNumWF5e\nHr9eWlrqzNNLUZRhsL6+jvX1de/9cy1Tki8BcFJEfkHyBRhZph8B8LCIPJHscy2Ay0RkT+KAWsVo\nnPR8APcAeIVthqplqihKFbhn8r2stnjtmpbpuQBuIXkaRuOrt4rIt0j+JclLMPLq/xTA7wOAiBwi\neRtGwwAnAbxHVVNRuoctWj7MHTwGADixd2G8zUfsuGd6P9/rp8em+7cprlNlmYXWqWWqKHFRRTx9\nSMUt6/ymEJbl6K3zWNh3ArixHS0JNmaqzB6z0s3yCax0lyLhSq1LYNLC9CEVRlMg57fS87nPle43\nd/AYnjrzpaWud/iME6MXN2XqG3aJ4Niz8+P3N55+HVa29o3fH18s9x3zUMu0Q8xvHSvdhVKUlFPC\ndooTexcy65FLeF1deWCyLprHmeLsIj3PH63+IX6XH8/dNxSXPX20kogWWabRiukGR2XeJTJ+nb63\n97G395m0QZzYuzCuqCGfrkp/cYkpcErQ8rrcD675tbX0GvsWV1oTR1/qakTnxdRFKrC3y/X4Y3ws\ndPE6Afe4n/oqrP3HFEXzfpt1IqseZNWb9Bj7YW2LoilI9rm+f/aZFb5NewxWTMdcR2zcNLkpFdNd\nHwJuuHG4gqoMkywxLcLlNc87t0scTUFyWbrHP37mVHuNgRA91yIxjT+fqeWpS3+UXSKtefFiZ37r\nmPNP6TeheiF2ndm3uFIopOn107/xua4/GqRMIdn1oXau0wlvft5TZUhWaVmBnN8addNMhvR79Qn7\n3qcCNr91LFNU02PSz02r1P4sj7JW3WVPH42uy3/s2XksnH680WvE382PnA2yUedXSAtz3+JKI2Ka\nNui8hq3UwyWmPt39rHuSVa9c3fQQ9fsGtOetd7Hj5FxtMe3+mGlkmFEGTVJHRPctrkzE0pmo2MVP\n3tima6ZPlYdYVv3at7iC6569cWJbKIvOjPccx4i2RBtjpiqmnrQdhtXEmKcKadzYYUmmp7zqvTO7\n8zGE1aVleOrMl7YqqG2IaSfGTGfOddlhWnUxB/2zrMnQ17KxG1SZ8TSlHfLuSdZn9v12zWgyZyDZ\nQwcu6tYJn2uEYsfJORw+4wRul+uxq4XrqWXqwbFn53H4jBOlxl2yHAZZn4fCtDpCncunAfl8HxXn\nbEyr0aToXuaNnWaJqWu2kq+Yhhh3byOwf8fJOdx4+nUAwjldtZsfAHPw3Le7kOUwSP+bXvaVrX2t\nPrGrcHxxYWK8rswc75S8qYuhiCF7UFVcoli2TuT1MkzBTn+fdIqy2e3Pu6bppa86LdNuB9c9e2PQ\nLv+Ok3N46dGnxu9DPcR7L6bOwOEmLKC0q+8R25pVJpeINkV6naxrZCWxyLOEyia+yDqveZ65g2Hz\nDYQQ0zwxsR8qTVP1AWu3AXu2UpboFlnGtrf/drl+XMeqCisQdvbUZU9Pxrq2JaY6ZupLgAkCVYTU\njhN1YZ8r75giQbSFrgp2Mgz7PKaATl+vXsVvSkRTTFGa35r+vMqQSN4xIXsspgVqRwCkr11KYZZh\n/vqj+P5N9YUvvX567tvl+mBd/u+ffeZY5NscWuqUZRqiUoX+cfOsUHvbDfhDAOUt0jxxrGLd+liY\ntgDWsUp9SUU2b+44UGw5lqXJ4ZUsL3oepjOySEyzHliuOftZ57L3TbEfTE04p5qwTFNCR930ppvf\ndKiQ/aQObQ3UtTBcgmoKqUvsTAsxawwzizL718mBWUQdZ1oox1ksuJxNrsQlVepxlgHgItTDLLSQ\npt37pqzRXnTz61T4PBErOzukCmZFLxrHzCPvmCoCVvaYMkKalfOyDCEiElLnSl7Xv6qjJyXvIRYa\nu4zzW8cgqwtTww1NPSDsrnleGYsELUQZ7eiaZieLFhO9mLZpOTR1rVOVr9k40ixsJ09b16l63ZBi\nNHfw2ITYuBr5aMbRpJXn0yWv+hDz+X6++4Wqs6FD28rEKqfxoL7sODkHINzMrFBELaYhKkpRsHMb\nYt3kNVze16Iuf5PYlqrPdYvGSpvCXFIjHSv0+X1c45S+5c87xowDnW/5t0ip00UuM5Rl7vvSo0/h\n+IeKU/eNsz+VFNHQY+xZRCumIyshf586Y5xdGiuriq+YNYUpFD5WqpmQGHAPbWQ9PLKmR/okAzHL\nVkUg7eN8cUU5AO4ogTbYt7iCG4BCL3jZYTMfbrjxevzuTW5vfghH0vGPT47LHluZjEUNQbQOKHNm\nRhZ1x7t8sAWpbevJhyxLtKyAuo7J++5ZY6RZv5FvbKuPmJZZZaBMl9M1a6iNez7rumVGmwD5ghrS\nq2+eK2++fh1Bnd86NuXkSgP7SyXXruPNJ/l8AN8GMA9gDsD/EJEDJBcBfBHAywEcBvB2EflFcswB\nAO8C8CyA94vIXY7zOsU0bwEvmyxPZh2KxCAWqlqZTVipLoEtO+5Y9vfOOiZEly11WgGzG35o65p1\nczL45nrwPUfe9NIQ1mndjG+1vPki8vck3yQivyZ5BoC/IfkGAG8FcLeIfIrkPgD7AewnuRPANQB2\nAjgfwD0kLxKR5woL6hBGV6NPvbMhhDQvZKgoNGgWgpvGH1YRxSbK6zqn3bXPIv2NQ3SP6zD9AJ9t\n/gDXZIaUputcUde+ickFPseFihdtOtubdzef5AsxslL/LYAvAXijiGySPAfAuoj8w8QqfU5EVpJj\nvgFgWUTus841YZlm/ZiuAfo0zKXs3PCswX7Xvj5liYGySTFcx9iUEeomrM4q1LFIs+qRbZUWZfWy\nhyaqhsHl/e6yOruJKyGsUJ/z293xmFYdrr0GFMnTSD4AYBPAvSLyMIBtIrKZ7LIJYFvy+jwAR4zD\nj2BkoZbGfEKnfyayCqdlk4ptGg+XN6XRFA77/Ob7Np03bVqbVS1D32vn7WNOQrDXEfJpoK71h0JR\n9jfxmfJbtQzpPYo1cUsof0V6H+159V2ijGV6JoBvAjgA4Msi8pvGZ1siskjyzwDcJyJfSLZ/FsBf\ni8iXrXNNWKZ2lz3PUWFnvDFxJXdIj8ujaizkrMdYm47hLDtDKo/UsrMtPN/UhGZ3MKSAlo0tNcmr\ni1UJHYfqIvSU29APNNMRFtOaZcFmQInIUZJfA/DbADZJniMiT5A8F8CTyW6PAdhuHHZBsm2K5eXl\n8et7/8MS3oLXAph0AhWJRVYAtklR176q1ZllyYYQ1qa6x2XPGdJh9cf4GFZwLFdITVyfNWGFjsRw\ndF67ldj31HR6Hl9cAAxrMUS+2jJCWmZ/k6q/oW/9DvHAi0VA19fXsb6+7r1/kTf/JQBOisgvSL4A\nI8v0IwDeAuBpEVkhuR/AWSKSOqBWAVyOxAEF4BW2697lzTe9qOnYkEv8qlqfQH0L1GeftqzUKmFP\noa9XxaPukzHe5zxNYNYl1/0sW6ayCU5CzbryuVbZ71Lko3Bdu2/JwOuGRr0awC0Yja2eBuBWEflP\nSWjUbQB+C9OhUddjFBp1EsAHROSbjvMWxpm68iva40bpwmNlrQCf7muocdLQ4po2ON/yFY1blnGQ\n+Di8qgpO3fM0Qd7CdmUJOfe/DHmhhWWoIqZFeRFixzbwOpc1yuxG2TfQte53HeqOCYayEnyoG1ua\nRxVRtc9bVfyGIqaAX53NizopS2gLtayg9kFMU44vLtT35rdN1hS7UGR1V+0/H0KNj7qiCfKuWbZ8\n+xZXMj3Oru1F3ulQQhoraSOahRA0Vf+z6lianyKEcWLXzdhCCctQ5feIzjJ15RXNcjQVBdb7UuWm\n51lzRTGsZeJDs8ZifeNFm6aukJqLq9VZ/qJLhByWKiKUgyrPKi26RlfvZ1nLNDoxzVpfx7eLUYWm\nhCdr7nrZ69tdP1vIq5436/y+hOjGlVnCo08UCWrZsfEs6kzVHS9lkjOxwfecKbF3+20BNd+fOPvM\nbiSHtm/YtJe+uUZWV4yyCDX25WJiPamD+7wSlGSJZp2ZYVmhaEMRxSZperKI78Pepk67CD0OHRLf\nJYiyiMIybdLqDEWo8dHQ5zWnLbpmiZlL+eZZoFWccXnfp0qyjKEJsM/aTnUoctgWTV11HVumF5NV\nP2K+z520TNMnVBeEFKieGzR0qFUW+xZXpixUuxGk/8sIeNkY2vH1S1gfXff6hiJU1z7FtALLnrNO\n8H2b069DUtcJNzPLFLunr+sSnrbEaFa4xkJ9k2RUdYJlWQyuxlzkHHMdF7v1ERNZFmNIUS1LUU8k\nrz7llbtrD8yp3LYFlunMQqPsRA52SEWZEKAuY1e8stmG0v1DOKHyzlN0HNC9xhIDroeO+YCy20ke\nIUL6mppJ2OW64ZtQZ2aW6VyN7DCuJ2TXLVczJtQkFVUfS9O24kOFjvle05yNppapP1m5fG3HYVt1\nvMq4ekqeaNZZ1iQGOhe0n4fZgG0hbdyKXdto9PR5FTXrM1cgdpXzVMUVoN1UWrw+I6vTv5v928Zs\nLKTt0Wdqd1uLWIbGp8zRhEbVwewKxVzpijixdwErB09171Mr1f5eWWPKru8e6vcoso5O7F0o5XRS\nppmVyNgJ18s4HE3yLM++CqhJp7r5ZZ0qecTk2MrqjvsMY8xqqKPvU0pngWsxv6aw67+Zl9XV0/OJ\nV07pW10YJ12K1QFVhqoJRezxG1k9laHfh3T/kKTnNM9bx6qs6jjyJcs5aNK3xhMDbQTsm9fgHneY\noqtu1jVkuoZv/Z5ZN7/MTBvfbrzrc1sMfU13HxGtEiebl7wamK6YeRZp00La9LmV2dHkrCczfWaX\nvfgmxxcXppKH28xMTH1uju9ULtd4Xp2baB6bJZYhKskfrY6WZ3DNXioiFpFzZZhXS7UdQvkIzLhg\n7hn9H783ruMzHGW22dH2ftQFH6MpagfUaAbHpKD6xMiVWZI2b8yx6dlZpojmhR7NUjjzrq+iGQf2\nhIkqM/SA6fo+v3UMcwezj7XHUVNv/lDrxcy7+fbNACa74i5BzYrB8+2KZFm8dlIVuxs/WfGMCmOH\nTe3elXntuYPHJmYu2WVOP8tLSDILcS0zTqqxpvWo00spOtY/8N89Tuo6v7keVtfveR0n4MwdUGl8\nndltzlpMLd2eFcx8Yq9/jGNePKQpoIXd+VRIMwTUnM2VljtdMCwNfcpKxnxi7wJOnH148v3ehcZj\nXrPK4rO8cpNLMA+BMo236Yeq6z7adVnv8ylmOmbqG2BuJmwokxILCJsvc6q8LovUQ+jmt45NptAz\nuvsfLVqZcW0j1/Jtir44EmKnyWGdPKvU1aXvYmxoaMrcj5lZprY1mmIuGVHkCMr6onlLMZif2Sug\nmueduJ5LIM1tprhZQpdak67rfHTPxzI/y8QlpGsblaxVM0pCiYMYHIsn9hYbLFntty/ktdssoshn\nmpI31mZap1k5MMss+5ySFUNnjuE6HVFZYuo4l891M61cXyu0psWaF+bl22iG7HwIQTpHv8oElKqx\n2FXo8z3OHTNd69iyJb4N0jWumeUsyiOvEjrF1CWi1ripd6KIImvSU6TrWjMusTSTb/S58cRE1eiR\nsjPl6tLn+mAbZGXENLebT3I7yXtJPkzyhyTfn2xfJnmE5P3J31XGMQdI/oTkIySvKPtlQt0o00R3\nmexFJrwpMBNis3vXqT+PMlRi9y7MPb3De/c66QpnvS68coq696LNJDd9p/TwG4rHTJ8BcK2IvArA\n6wC8l+TFAATAp0Xk0uTv6wBAcieAawDsBHAlgM+QrDUuG2oJWtc2e3tezJ2XWCWWZtaNmLqmLciJ\nSLumb2Zdv2hZ5iJ8rCEVynjxDZfTLr4fdQyTUt18kncAOAjg9QB+JSI3WZ8fAPCciKwk778BYFlE\n7rP2c2baN6n6lDa7p0Vzik18KmPmPnZ4lDF+aceV+l4zK3Y2ZLIXm7yptxpX2h5FyypnPaiL6lMI\nhnKP7XsgqwHzmZLcAeBSAKkwvo/kgyRvJnlWsu08AEeMw44AON/3GiEwc0OaSUXMtXBcc9orVzZX\nnKlj7NRlQRZZnS6qHNMkarW2gDG2XvWBqd38cri0owgvMSX5YgC3A/iAiPwKwJ8DuBDAJQAeB3BT\nzuHte7g8sMdTgXqWax4rW/vGlqkpqiGsTPN8psiWFdsyPQFX6NlQLJZYUbEMS5pFy/wrojBon+Tz\nAHwJwH8TkTsAQESeND7/LIA7k7ePAdhuHH5Bsm2ah5ZPvX7ZErBtqbi0NSk1QO8TalQw+8nH8vXt\nmmVNl/VZUK8qWUlMVDibZdSlzK5/s8jX0KcMUN5srgNPrnvvnjtmSpIAbgHwtIhca2w/V0QeT15f\nC+AyEdmTOKBWAVyOUff+HgCvsOOgfMZMgbCzbiqFndiCmiWeBaKaYo6fNhkXWDXpxeAaSwcY11vH\nw71OUhMfQmVh6xK5OlEnNAojR9PvAXiTFQa1QnKD5IMA3gjgWgAQkUMAbgNwCMDXAbwnM6B0opDu\neEtyY2Ric/rzsfnt+CwoZWYWFcxEOrF3ITPBSdb+VSwQu7sfy/iqUp6xiDke1HbPp6n7PBQhrcvM\ngvaBB6eD3gF3IPzaBkR2GcdvTO1nfu68Zoh0elaZcinIHOWDz1huGXyEWRtOfIRKBen7ADcZWn2o\nY5nGl8807c5Y3vH0S45vriVmpoUqsmvSYp1BYpAxNaZ5NmlRzjpPqhIfdmrLoVH3oTVbMbUt0ixr\nr6Qg1er6F5XJMzOU85zJdwhtcValTI5SZXakeRNyc0VUJJa6GDWe7X12YpolSlnOneQ11ypcy+6e\nZwmzNWTAPQXlKbrmLC3iAlzRAZqoJF5CdLddyXtOnXchSfC80PgKE52ihOE02zHTLFxjqaEoENIU\nkV3uBCemIFcpn8Mj66LM7C0ffBNgqJjGj5mMo2x33FeUp8W23xRmhgMAvKZjY6YpAYV0PIZa0sqc\nSEtX1UKtSBtOJxXS7pF297mn3LjmUESxCqEs8ZkvW+IksFBljqHmdMPTYzIrYdkuvEdgf9sM0cnQ\nR3S8szqZQlpBg+IT04bGGUV2Vc5I75ynG6CLn0cIoSu7yKASP2n9S/83FZUxBEt2SkhTfahozMUn\npg0yjkVNf6yKP5qIkdPURyAL9snKAhQCn26+Cmk3kdVyD0wlhwC94fjEtCEv+FRX3+PHs48hNya6\n/+M/2TXxV4U2utza6PpL1r2tc8+HYJ2GJD4HVFUvuSdTAf0FVIlZnZitVTC43WSQtJkT1ZV0xQyV\n0bCobjJyRi2Mp5MGSSs5AJqYhh6fmM5gTfiyjK3THCs0a+XTlLYqelaSaUAtj74gq8D81qxLocTX\nzW8QHxEMRZGV5xLaLPGtYzHa616l57PPqVZpt8lLilPnwd3H5N9NJUcalJhO4Os8CkyZih0id6i5\nwkCZrOFKt3Dd1yZy2irZDFdMW8AMY3GNabmy/JvWY13Mc6iI9p+yFqi5YoKPBdpHKzUkcU4nbZAp\n51Cd8KgSmBXRZ0E0l/iVqcymkKpzaRgU1Q/XygkpQ8ukX62rnz+ddNiWaQTOLpc3P2SlVmtC8aFs\nnStr1Q6B+Lz5TeLKi1qBJh1YeeU6vrigFVfJpKh+aN1plsFZpua4UhVRnJWQ+mCOwZpr2ft08bWh\n9QP7fruWGA95HV1g8RTDEVMrMbPrs6ltLXj70/I0MV6lAjlc2hDUXlDUxktoQP+7+caP4RyfTPOW\numZeea4v5UNeF6yKkGqXXyliJKgfm6gnWfXG3qfPTDif8vwmJY2p/oppwYJ281sZlcY6rolufdY0\nTx/KVnT15CtaZ9qhm918x4J7zu0Grhi8qSe0eay1ImpomkqdVnW2izae4aD3GtP6EWBILzfOlOR2\nAH8J4GUABMB/FZE/JbkI4IsAXg7gMIC3i8gvkmMOAHgXgGcBvF9E7nKc1y/OdPeuUSKH1Cyv8YVt\nkcm0DBsUUVf3qk58n283XxuPkkVWHep7nSmdXX9tA3XjTJ8BcK2IvArA6wC8l+TFAPYDuFtELgLw\nreQ9SO4EcA2AnQCuBPAZkuWtX+NJMR7PDCikKa41wtuYtx+KtMJrxnylKkP0xje1YGCu0InIEyLy\nQPL6VwB+BOB8AG8FcEuy2y0A3pa8vhrAmog8IyKHATwK4HLnyV1mdotz5W2BdVmHoX90V6UN4cWX\n1f5bEkpccE9zotQ0ZptLp3kX4qFN3g4okjsAXArguwC2ichm8tEmgG3J6/MA3GccdgQj8fUvYGBB\nrTJf+dSxANAdkXJ5alVkFR/KRIfMbx1L2sZk6r8u1bVR2kKzredP9/bBqwtO8sUAvgTgAyLyy4lC\njQZd8yb4B5/87/M08X7iJLQ1LS5khSs7F1tRfMirN0UzrMy/rlquVZ24hZYpyedhJKS3isgdyeZN\nkueIyBMkzwXwZLL9MQDbjcMvSLZN89DyqdcvWwK2LXkV2P6idXI1Npnl3ocmxC61MFRIlbKkdSZU\n/Rm1r3jrYaE1vrkOPLnufb5cy5QkAdwM4JCI/Inx0VcBvDN5/U4Adxjb30FyjuSFAF4J4HvOk796\nefw391ev9S5wKnyhBbBMsuaQNGEN120ImrxCsalaJ7pqnQIAti1h7tsfxNy3PzjSqgKKLNPXA/g9\nABsk70+2HQDwSQC3kXw3ktAoABCRQyRvA3AIwEkA75GCHH+pdei7fo2Z/zOU2GWdq6mUZK5UaFUE\nUGdBKU0QulfDPfHm081qQ1V6vLPLZ7p7dN1YFv1yrZPURlfZTEgSCzpMoNiEeGjHWqd8v9uJs8+M\nO59pLDGSqXXa1DBCk1TpgmlXXmmbWOtbKJGfmWVqXtfnR85aYTMUrgD+rlAlScWQElso9QkthDFm\n9i/6jtFbpiazXBp5Is9pZDfZpm6G81gtBGU4zB2ML3TKNCrsfK0+BkcUYpoWNnYRiwV7facyaFC/\nEhMxCmrVNhGFmJrMunF3UdCrdO0VZdbE4nwORXRiqviRZZ1mCaYKqVKVputObNZpVfqbHLoAp8Op\nY1apHSPnctIVNYQ0BjDGEC1l9uhD2J/oLNM2bl6Xwp7KMHewfHxoF4c1lG7T5ciZPKIT0zbwSb83\nNHQVU6Ut7PbX5XR+Jr0XU5dATCyo13EhLRLAvM/tZaGLhNJ3P6UfNClwrt5h10W102JqzlhSRhQG\nHmc8SOrErKq4KlWwc3GkFmtXRTWKGVA2duP0EcwyYRZNrlU/K2xHlEsos8aqysyGMh1VLhFVB1Z/\n4J5w4Uumc7SKARRDWyXZnRlQKU02yL5asnm/mTkpIlSlzAvBUku1H4SqK0X5LszMcVni3QVLNdrQ\nqBCiZ6brC3XOPmA3klSIuQeF4WFmAuE8NPNU9xlZpuHPm7VESFH+jZhT+QGRWqZFlF2SxDW2GvNN\nqUpd8Srzm/h4/12Cq1brsCgyYFxLsOcxXvY9QqIVU1fDNi1M2+pURrRtfZcVVLVWO8TaRuX65OMc\nTvcplRFu965ou/zRiqkP5lhL3ud5+/SNdFzUtbhZrJVQiRORaisF+wiwvY/PKhsmMdblTolp3k3K\nW8OpSHT7iG0BNp0LVsPU+omslgtB9LFGM1nb8DpHrEQtpmlXv2iMNO0m2DdhSOLpIh3XbHNIJG14\nddIEKt0nr806xXLNfyw01nYdZZzpxL4F5nwqoraH0PWD92nmU4yY98qOYY0xs7riB7kB7B51+YuM\nmqZJrz+LsfdOxpnWwTc4OMYxl1io+tuk47WmaKaVPlZrQikmHTvNMlBcwwBB7vfaxpTFGvMQQLRx\npiZF4pjlaMq8oekNWq02wN53QluQ6sHvL43HcO+ebqOx1qdCy5Tk50huknzI2LZM8gjJ+5O/q4zP\nDpD8CclHSF4RsrBlnoyFN9dxkxRFySAJkwrliPK5Xtfw6eZ/HsCV1jYB8GkRuTT5+zoAkNwJ4BoA\nO5NjPkMyyFCCb7hF1g2f2OYhpDoMMA25MQqx4sb4z7lPxIHVSnja7HrHPO5eKHQi8h0AP3d85BqI\nvRrAmog8IyKHATwK4PKqhUsFrcpa9qljyvRk+47jqJBO4y2QavH3EjPm1DZYGhsPN+rS3MFjUQsp\nUM8B9T6SD5K8meRZybbzABwx9jkC4Pwa1wBQ8mZZsWpVb7Q6TDIwul+ZQd0qqP3E0fUOGl9snr+D\ndaiqmP45gAsBXALgcQA35exbO/Zq4malHr70h7ffp9sy8BFJV8zq0DC78mOrtEQFV+u+f7genkGN\njg4KqEklb76IPJm+JvlZAHcmbx8DsN3Y9YJk2xTLy8vj10tLS1haWpr4fKoxugQy70m2NoqNKzXv\nd3xdFVLf/cwGZsYjKj1lLfse980IWV9fx/r6uvf+XkH7JHcAuFNEXp28P1dEHk9eXwvgMhHZkzig\nVjEaJz3H+FoJAAANRElEQVQfwD0AXmFH6PsE7Y/FtIxXb/euyZu9lm9RmWMw3IOJY6ssTtcXJsTU\n/C3T39ciFVSXmMY+zqWUJ++hWSpUKkeYXcy6LhUF7RdapiTXALwRwEtI/gzAhwEskbwEoy78TwH8\nPgCIyCGStwE4BOAkgPd4T3XKwhbGov2KtiU4b8zArSqnRZpWePO3se5JlWEApbvkPTxLWaY9qy+F\nYioiux2bP5ez/8cBfLxOoSrhacH6PN3U+VQR+x70rLEok4jEmw5vFsQ7nbROwwwQ8NunsZ9apF17\nl9Ov6Dil/1Rtax0Myi8iWjF1ht3Y3U0XZbv2AfbtNTUrvQbwK056+LCNem6+rHp69YHCsdUy4nhi\n70LhWkh9obbYmU4ps4EY22Jfu0eZESUcUF2oP9FapqUxG7TjBuWN7Qxx3GcqfrSOpeDh/BvibzwI\ndu8qvSabeWyfiF5Mu/BE6hpjL2xdEVUUNONfMMW5KxoQvZhm0qAIdOXmVSH4GKbtmMq7tlqn/cMx\nfdt3lqG5nJDrddeIesx0jGtspaEwnNGNVE9+aUrG+Sr9oaz45eUfTkVZVtE5v0X3LFNtoPWY8e+n\n1mnPKBG072N12is1dInuiWlF8m6Qa5mN3pLXFTfHUeuKbs51VFD7i8+y632lP2Ka03h9nnRdfiKW\nwWstdJ+gfNcxymDZt7gyfl11PaiuGzL9EVPFi1aD6FVgB0XQ3KYdpDti2vHEsbFQaJmGFEC9T73G\nfDCvbO2bYUnioBNiOhaAovE+g7RrMYSuewjGzoE6AqgPvMFRZuXgvtMJMQUcFlXaWDOcJUPubnjR\nRBdcBVSpSNfHSwHP5NDBL+qRHNp5nIcX2HwS9uEGNYX9W5q/W+6DKFB8r/YYug/3GJbp2YdHG61F\n8HzpQlutnRxaGQZll3dxovlMh40j5tSnTnVBSH3oTDcf8LNmUo9iX25QVOStw+U5bDCkMbShMyQh\nBXpqmWoXsh7ptL6p33HVWK4CmE6/55lSTe9PD9FeSLcsUx+0oVbHHCvNnTFmR1eUSHaijsGe4rjv\nQ7vXnbdM0wTSKqL1OPX7NdsA9D4NBzNxiavL36cuPtAxb74SJ1OzqhwxvxPZgJReML/lv6xzH5ZO\nL/Lm966br8wIM+7XwMxTqfSHKslqUvHtK4ViSvJzJDdJPmRsWyR5N8kfk7yL5FnGZwdI/oTkIySv\naKrgSjyMx1DVCTEYyj4chzBv38cy/TyAK61t+wHcLSIXAfhW8h4kdwK4BsDO5JjPkFTrd8CYjaiJ\nLn7frZ0YMX/zvgtkGQqFTkS+A+Dn1ua3ArgleX0LgLclr68GsCYiz4jIYQCPArg8TFGVmCnKF9vk\neBn3aI7UWVFkoZoJofs+Xl7Vm79NRDaT15sAtiWvzwNwn7HfEQDnV7yG0jHabixqlc6G44sL3r/9\nOF6550IKBAiNEhEhmeeaV7e90ih98BR3ifmtY+Nk0B/d87HcfftujZpUFdNNkueIyBMkzwXwZLL9\nMQDbjf0uSLZNsby8PH69tLSEpaWlikVRhkpqIamQKk2wvr6O9fV17/294kxJ7gBwp4i8Onn/KQBP\ni8gKyf0AzhKR/YkDahWjcdLzAdwD4BV2UKnGmQ4DMw4xpIWSZitSMZ0Ndhc/ywnVN6u0dtYokmsA\n3gjgJSR/BuAGAJ8EcBvJdwM4DODtACAih0jeBuAQgJMA3qOqqTSJCmn7mGOmQxFSH3QGlNIYaYML\nKXhNnFMpzxDFVPOZKjNj3NACNCz13HeHPgqpDyqmSmM01ajUKlViRMVUiR61SuMiq4s/VIs0Rad6\nKp1CrdJ40OQ1k6iYKtFzfHFh/KfEy9Cn9KqYKoqiBEDFVFEUb4ZufeahYqooSiXUATWJiqmiKN6Y\ngqkOqElUTBVFqYTvarZDQcVUUZTSmEKqFuoIFVNFUUphL9+sS5eM0EQniqKUJvXqy+rk6z6jSz0r\nitIoqYU69LApnZuvKEpp+m6FVkG7+YqiKB5oN19RFKUFVEwVRVECoGKqKIoSABVTRVGUAKiYKoqi\nBEDFVFEUJQC14kxJHgZwDMCzAJ4RkctJLgL4IoCXAzgM4O0i8oua5VQURYmaupapAFgSkUtF5PJk\n234Ad4vIRQC+lbxXFEXpNSG6+XYQ61sB3JK8vgXA2wJcQ1EUJWpqzYAi+XcAjmLUzf8vIvIXJH8u\nIr+ZfE4AW+l74zidAaUoA2CDp2ytXR1v80UzoOrOzX+9iDxO8qUA7ib5iPmhiAjJbv+CiqJUwhTS\nIVBLTEXk8eT/UyS/AuByAJskzxGRJ0ieC+BJ17HLy8vj10tLS1haWqpTlEGRVtL0SW9X2qYtAPv6\nimLTdp1sgvX1dayvr3vvX7mbT/KFAE4XkV+SfBGAuwB8BMCbATwtIisk9wM4S0T2W8d2opuf9WRt\nomIUXWuDxC6RUk/7LLG1PzfL4LqGz3W72FiU8LjqSV/qRlE3v46YXgjgK8nbMwB8QUQ+kYRG3Qbg\nt5ARGhWrmNbtllStNH3uDvWlISnFtGl8zILGxLQOMYppU4JWVJH6LKQmfWlQSjYqpj0X0yxv4ixE\nbNbXj5W+NLahYY6dD2EYqGlvfuNUcXZk3dh0THBWqIC6UYdWvxnKfZ2ZZfpg61dVusRQGmAXKWsU\ndPFeur7ja4BuW6bKMOmzV7jL9LF3Feo7qZgqnWHWwzRdoUlHUNnwvJgJ/T1UTBVlIPRFBKvS9PdX\nMVU6QV8t0rpecHXeFdPWQ0TFVImaPolElUbte0ysQyBDip5RMVWiJUZxcBFL9zmWcpjMYo7+rH4H\nFVMlSrogpDGKV+z0eVhCxVRRSqACGoaQFmss90TFVIkSTTozLMpmQ4vxPqmYKr0kxsamhCHWe6tL\nPSuKogRAxVRRFCUAKqaKoigBUDFVFEUJgIqpoihKAFRMFUVRAqBiqiiKEgAVU0VRlAComCqKogSg\nETEleSXJR0j+hOS+Jq6hKIoSE8HFlOTpAA4CuBLATgC7SV4c+jpt8P1ZF8ATLWdYtJxhGUo5m7BM\nLwfwqIgcFpFnAPx3AFc3cJ3G+cGsC+CJljMsWs6wDKWcTYjp+QB+Zrw/kmxTFEXpLU2Iaf+yviqK\nohRACZzxmuTrACyLyJXJ+wMAnhORFWMfFVxFUTqHiGTm/2tCTM8A8L8B/AsA/wfA9wDsFpEfBb2Q\noihKRARPDi0iJ0nuBfBNAKcDuFmFVFGUvhPcMlUURRkirc6AiimYn+TnSG6SfMjYtkjybpI/JnkX\nybOMzw4k5X6E5BUtlnM7yXtJPkzyhyTfH2NZST6f5HdJPkDyEMlPxFhO49qnk7yf5J2xlpPkYZIb\nSTm/F3E5zyJ5O8kfJff+H8dWTpL/IPkd07+jJN8ftJwi0sofRl3+RwHsAPA8AA8AuLit6zvK8zsA\nLgXwkLHtUwD+Y/J6H4BPJq93JuV9XlL+RwGc1lI5zwFwSfL6xRiNR18caVlfmPw/A8B9AN4QYzmT\n638QwBcAfDXie/9TAIvWthjLeQuAdxn3/swYy2mU9zQAjwPYHrKcbX6BfwLgG8b7/QD2t/kjOsq0\nA5Ni+giAbcnrcwA8krw+AGCfsd83ALxuRmW+A8CbYy4rgBdiNKHkVTGWE8AFAO4B8CYAd8Z67xMx\nPdvaFlU5E+H8O8f2qMpple0KAN8JXc42u/ldCObfJiKbyetNANuS1+dhVN6UmZSd5A6MrOnvIsKy\nkjyN5ANJee4VkYdjLCeA/wzgDwA8Z2yLsZwC4B6SPyD575NtsZXzQgBPkfw8yb8l+RckXxRhOU3e\nAWAteR2snG2Kaac8XTJ6HOWVudXvQ/LFAL4E4AMi8suJgkRSVhF5TkQuwcjy+2ck32R9PvNykvxX\nAJ4UkfsBOGMGYyhnwutF5FIAVwF4L8nfmShEHOU8A8BrAXxGRF4L4P9h1Os8VYg4ygkAIDkH4F8D\n+KupQtQsZ5ti+hhGYxQp2zGp/DGwSfIcACB5LoAnk+122S9ItrUCyedhJKS3isgdMZcVAETkKICv\nAfjtCMv5TwG8leRPMbJO/jnJWyMsJ0Tk8eT/UwC+glHei9jKeQTAERFJ84TcjpG4PhFZOVOuAvC/\nkt8UCPh7timmPwDwSpI7kqfDNQC+2uL1ffgqgHcmr9+J0fhkuv0dJOdIXgjglRhNRmgckgRwM4BD\nIvInsZaV5EtSTyjJFwD4lwDuj62cInK9iGwXkQsx6u79TxH5N7GVk+QLSf5G8vpFGI3zPRRbOUXk\nCQA/I3lRsunNAB4GcGdM5TTYjVNd/LQ8YcrZ8sDvVRh5ox8FcKDNazvKsobRDK0TGI3l/jsAixg5\nJn4M4C4AZxn7X5+U+xEAb2mxnG/AaGzvAYzE6X6M0htGVVYArwbwt0k5NwD8QbI9qnJaZX4jTnnz\noyonRmORDyR/P0zbS2zlTK77Gowcjg8C+DJGTqkYy/kiAP8XwG8Y24KVU4P2FUVRAqDLliiKogRA\nxVRRFCUAKqaKoigBUDFVFEUJgIqpoihKAFRMFUVRAqBiqiiKEgAVU0VRlAD8f9xUIQkX3i0tAAAA\nAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x11221e890>"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "##Packed integer data\n",
      "There is a similar feature for variables with `scale_factor` and `add_offset` attributes.\n",
      "\n",
      "- short integer data will automatically be returned as float data, with the scale and offset applied.  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Dealing with dates and times\n",
      "- time variables usually measure relative to a fixed date using a certain calendar, with units specified like ***`hours since YY:MM:DD hh-mm-ss`***.\n",
      "- **`num2date`** and **`date2num`** convenience functions provided to convert between these numeric time coordinates and handy python datetime instances.  \n",
      "- **`date2index`** finds the time index corresponding to a datetime instance."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from netCDF4 import num2date, date2num, date2index\n",
      "timedim = sfctmp.dimensions[0]\n",
      "print timedim\n",
      "times = gfs.variables[timedim]\n",
      "print times[:]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "time\n",
        "[  0.   3.   6.   9.  12.  15.  18.  21.  24.  27.  30.  33.  36.  39.  42.\n",
        "  45.  48.  51.  54.  57.  60.  63.  66.  69.  72.  75.  78.  81.  84.  87.\n",
        "  90.  93.]\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dates = num2date(times[:], times.units)\n",
      "print [date.strftime('%Y%m%d%H') for date in dates]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "['2014101812', '2014101815', '2014101818', '2014101821', '2014101900', '2014101903', '2014101906', '2014101909', '2014101912', '2014101915', '2014101918', '2014101921', '2014102000', '2014102003', '2014102006', '2014102009', '2014102012', '2014102015', '2014102018', '2014102021', '2014102100', '2014102103', '2014102106', '2014102109', '2014102112', '2014102115', '2014102118', '2014102121', '2014102200', '2014102203', '2014102206', '2014102209']\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "###Get index associated with a specified date, extract forecast data for that date."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from datetime import datetime, timedelta\n",
      "date = datetime.now() + timedelta(days=3)\n",
      "print date\n",
      "ntime = date2index(date,times,select='nearest')\n",
      "print ntime, dates[ntime]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2014-10-21 09:55:04.939928\n",
        "23 2014-10-21 09:00:00\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "###Get temp forecast for Boulder (near 40N, -105W)\n",
      "- use function **`getcloses_ij`** we created before..."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "lats, lons = gfs.variables['lat'][:], gfs.variables['lon'][:]\n",
      "# lats, lons are 1-d. Make them 2-d using numpy.meshgrid.\n",
      "lons, lats = np.meshgrid(lons,lats)\n",
      "j, i = getclosest_ij(lats,lons,40,-105)\n",
      "fcst_temp = sfctmp[ntime,j,i]\n",
      "print 'Boulder forecast valid at %s UTC = %5.1f %s' % (dates[ntime],fcst_temp,sfctmp.units)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Boulder forecast valid at 2014-10-21 09:00:00 UTC = 290.8 K\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "##Exercise\n",
      "\n",
      "Get probability of precip (3 hour accumulation > 0.25mm) at Boulder for tomorrow from [SREF](http://www.emc.ncep.noaa.gov/mmb/SREF/SREF.html).\n",
      "\n",
      "(catalog URL = 'http://thredds.ucar.edu/thredds/catalog/grib/NCEP/SREF/CONUS_40km/ensprod/catalog.html')"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "##Simple multi-file aggregation\n",
      "\n",
      "What if you have a bunch of netcdf files, each with data for a different year, and you want to access all the data as if it were in one file?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!ls -l data/prmsl*nc"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-rw-r--r--  1 jsw  staff  8985332 Oct 17 19:53 data/prmsl.2000.nc\r\n",
        "-rw-r--r--  1 jsw  staff  8968789 Oct 17 19:53 data/prmsl.2001.nc\r\n",
        "-rw-r--r--  1 jsw  staff  8972796 Oct 17 19:53 data/prmsl.2002.nc\r\n",
        "-rw-r--r--  1 jsw  staff  8974435 Oct 17 19:53 data/prmsl.2003.nc\r\n",
        "-rw-r--r--  1 jsw  staff  8997438 Oct 17 19:53 data/prmsl.2004.nc\r\n",
        "-rw-r--r--  1 jsw  staff  8976678 Oct 17 19:53 data/prmsl.2005.nc\r\n",
        "-rw-r--r--  1 jsw  staff  8969714 Oct 17 19:53 data/prmsl.2006.nc\r\n",
        "-rw-r--r--  1 jsw  staff  8974360 Oct 17 19:53 data/prmsl.2007.nc\r\n",
        "-rw-r--r--  1 jsw  staff  8994260 Oct 17 19:53 data/prmsl.2008.nc\r\n",
        "-rw-r--r--  1 jsw  staff  8974678 Oct 17 19:53 data/prmsl.2009.nc\r\n",
        "-rw-r--r--  1 jsw  staff  8970732 Oct 17 19:53 data/prmsl.2010.nc\r\n",
        "-rw-r--r--  1 jsw  staff  8976285 Oct 17 19:53 data/prmsl.2011.nc\r\n"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "**`MFDataset`** uses file globbing to patch together all the files into one big Dataset.\n",
      "You can also pass it a list of specific files.\n",
      "\n",
      "Limitations:\n",
      "\n",
      "- It can only  aggregate the data along the leftmost dimension of each variable.\n",
      "- only works with `NETCDF3`, or `NETCDF4_CLASSIC` formatted files.\n",
      "- kind of slow."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mf = netCDF4.MFDataset('data/prmsl*nc')\n",
      "times = mf.variables['time']\n",
      "dates = num2date(times[:],times.units)\n",
      "print 'starting date = ',dates[0]\n",
      "print 'ending date = ',dates[-1]\n",
      "prmsl = mf.variables['prmsl']\n",
      "print times.shape, prmsl.dimensions, prmsl.shape"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "starting date =  2000-01-01 00:00:00\n",
        "ending date =  2011-12-31 00:00:00\n",
        "(4383,) (u'time', u'lat', u'lon') (4383, 91, 180)\n"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Closing your netCDF file\n",
      "\n",
      "It's good to close netCDF files, but not actually necessary when Dataset is open for read access only.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f.close()\n",
      "gfs.close()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##That's it!\n",
      "\n",
      "Now you're ready to start exploring your data interactively.\n",
      "\n",
      "To be continued with **Writing netCDF data** ...."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}