{ "metadata": { "name": "", "signature": "sha256:dd72246b5115e614f58e416dbaf20bdc9b9fd21cd509e3765fca519ab90d3930" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": { "internals": { "slide_helper": "subslide_end", "slide_type": "subslide" }, "slide_helper": "slide_end", "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.\n", "\n", "**Important Note**: To run this notebook, you will need the data files from the github repository (the data is not included in the source tarball release). Please go to https://github.com/Unidata/netcdf4-python and follow the instructions for cloning the repository." ] }, { "cell_type": "markdown", "metadata": { "internals": { "slide_type": "subslide" }, "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 and numpy" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from __future__ import print_function # make sure print behaves the same in 2.7 and 3.x\n", "import netCDF4\n", "import numpy as np" ], "language": "python", "metadata": { "internals": { "frag_number": 2, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 2, "slide_type": "subslide" }, "slideshow": { "slide_type": "slide" } }, "source": [ "## Create a netCDF4.Dataset object\n", "- **`f`** is a `Dataset` object, representing an open netCDF file.\n", "- printing the object gives you summary information, similar to *`ncdump -h`*." ] }, { "cell_type": "code", "collapsed": false, "input": [ "f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc')\n", "print(f) " ], "language": "python", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 4, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "\n", "root group (NETCDF4_CLASSIC data model, file format HDF5):\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": { "internals": { "frag_helper": "fragment_end", "frag_number": 4, "slide_type": "subslide" }, "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": [ "print(f.variables.keys()) # get all variable names\n", "temp = f.variables['temperature'] # temperature variable\n", "print(temp) " ], "language": "python", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 6, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[u'MT', u'Date', u'Depth', u'Y', u'X', u'Latitude', u'Longitude', u'u', u'v', u'temperature', u'salinity']\n", "\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": { "internals": { "frag_helper": "fragment_end", "frag_number": 6, "slide_type": "subslide" }, "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 d in f.dimensions.items():\n", " print(d)" ], "language": "python", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 8 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "(u'MT', (unlimited): name = 'MT', size = 1\n", ")\n", "(u'Y', : name = 'Y', size = 850\n", ")\n", "(u'X', : name = 'X', size = 712\n", ")\n", "(u'Depth', : name = 'Depth', size = 10\n", ")\n" ] } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 9 }, "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": { "internals": { "frag_helper": "fragment_end", "frag_number": 10 }, "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": { "internals": { "frag_helper": "fragment_end", "frag_number": 11, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 6, "text": [ "(1, 10, 850, 712)" ] } ], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 11, "slide_type": "subslide" }, "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": { "internals": { "frag_helper": "fragment_end", "frag_number": 13, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "\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", "\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": { "internals": { "frag_helper": "fragment_end", "frag_number": 13, "slide_type": "subslide" }, "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": { "internals": { "frag_helper": "fragment_end", "frag_number": 15 }, "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": { "internals": { "frag_helper": "fragment_end", "frag_number": 16 }, "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: %s' % repr(temp.shape))\n", "tempslice = temp[0, dpth > 400, yy > yy.max()/2, xx > xx.max()/2]\n", "print('shape of temp slice: %s' % repr(tempslice.shape))" ], "language": "python", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 17, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "shape of temp variable: (1, 10, 850, 712)\n", "shape of temp slice: (6, 425, 356)" ] }, { "output_type": "stream", "stream": "stdout", "text": [ "\n" ] } ], "prompt_number": 10 }, { "cell_type": "markdown", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 17, "slide_type": "subslide" }, "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": { "internals": { "frag_helper": "fragment_end", "frag_number": 19 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "\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": { "internals": { "frag_helper": "fragment_end", "frag_number": 20, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "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": [ "# extract lat/lon values (in degrees) to numpy arrays\n", "latvals = lat[:]; lonvals = lon[:] \n", "# a function to find the index of the point closest pt\n", "# (in squared distance) to give lat/lon value.\n", "def getclosest_ij(lats,lons,latpt,lonpt):\n", " # find squared distance of every point on grid\n", " dist_sq = (lats-latpt)**2 + (lons-lonpt)**2 \n", " # 1D index of minimum dist_sq element\n", " minindex_flattened = dist_sq.argmin() \n", " # Get 2D index for latvals and lonvals arrays from 1D index\n", " return np.unravel_index(minindex_flattened, lats.shape)\n", "iy_min, ix_min = getclosest_ij(latvals, lonvals, 50., -140)" ], "language": "python", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 20, "slide_type": "subslide" }, "slideshow": { "slide_type": "slide" } }, "outputs": [], "prompt_number": 12 }, { "cell_type": "markdown", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 22 }, "slideshow": { "slide_type": "fragment" } }, "source": [ "### Now we have all the information we need to find our answer.\n" ] }, { "cell_type": "markdown", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 23 }, "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": { "internals": { "frag_helper": "fragment_end", "frag_number": 24 }, "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('%7.4f %s' % (temp[0,0,iy_min,ix_min], temp.units))\n", "print('%7.4f %s' % (sal[0,0,iy_min,ix_min], sal.units))" ], "language": "python", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 25, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 6.4631 degC\n", "32.6572 psu\n" ] } ], "prompt_number": 13 }, { "cell_type": "markdown", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 25, "slide_type": "subslide" }, "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": { "internals": { "frag_helper": "fragment_end", "frag_number": 27 }, "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." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import datetime\n", "date = datetime.datetime.now()\n", "# build URL for latest synoptic analysis time\n", "URL = 'http://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p5deg/GFS_Global_0p5deg_%04i%02i%02i_%02i%02i.grib2/GC' %\\\n", "(date.year,date.month,date.day,6*(date.hour//6),0)\n", "# keep moving back 6 hours until a valid URL found\n", "validURL = False; ncount = 0\n", "while (not validURL and ncount < 10):\n", " print(URL)\n", " try:\n", " gfs = netCDF4.Dataset(URL)\n", " validURL = True\n", " except RuntimeError:\n", " date -= datetime.timedelta(hours=6)\n", " ncount += 1 " ], "language": "python", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 28, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "http://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p5deg/GFS_Global_0p5deg_20150310_1200.grib2/GC\n" ] } ], "prompt_number": 14 }, { "cell_type": "code", "collapsed": false, "input": [ "# Look at metadata for a specific variable\n", "# gfs.variables.keys() will show all available variables.\n", "sfctmp = gfs.variables['Temperature_surface']\n", "# get info about sfctmp\n", "print(sfctmp)\n", "# print coord vars associated with this variable\n", "for dname in sfctmp.dimensions: \n", " print(gfs.variables[dname])" ], "language": "python", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 28, "slide_helper": "subslide_end", "slide_type": "subslide" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "slide" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "\n", "float32 Temperature_surface(time1, lat, lon)\n", " long_name: Temperature @ Ground or water surface\n", " units: K\n", " missing_value: nan\n", " abbreviation: TMP\n", " coordinates: time1 \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 = (93, 361, 720)\n", "filling off\n", "\n", "\n", "float64 time1(time1)\n", " units: Hour since 2015-03-10T12:00:00Z\n", " standard_name: time\n", " long_name: GRIB forecast or observation time\n", " calendar: proleptic_gregorian\n", " _CoordinateAxisType: Time\n", "unlimited dimensions: \n", "current shape = (93,)\n", "filling off\n", "\n", "\n", "float32 lat(lat)\n", " units: degrees_north\n", " _CoordinateAxisType: Lat\n", "unlimited dimensions: \n", "current shape = (361,)\n", "filling off\n", "\n", "\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": { "internals": { "frag_helper": "fragment_end", "frag_number": 28, "slide_type": "subslide" }, "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", "# flip the data in latitude so North Hemisphere is up on the plot\n", "soilm = soilmvar[0,0,::-1,:] \n", "print('shape=%s, type=%s, missing_value=%s' % \\\n", " (soilm.shape, type(soilm), soilmvar.missing_value))\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "cs = plt.contourf(soilm)" ], "language": "python", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 31 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "shape=(361, 720), type=, missing_value=nan\n" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD7CAYAAAB37B+tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+sZVd13z8L228SIGPzgNjGdjpWwSlueAFUm7Qk4iVN\njJFaTNsImFFbVNIoamJAFLczHhd7MIEyyKFRMyJqA0ROxJtCSKC4isF2wlMTVZgQsAfHuGCFkTIU\njx0PeEhRZ2y8+sc9+71999v7nH1+3XvOvesjPb17z899z9n7e9ZZe+21RVUxDMMwxssz5l0AwzAM\nox0m5IZhGCPHhNwwDGPkmJAbhmGMHBNywzCMkWNCbhiGMXLOncdJRcRiHg3DMBqgqhIuK7XIReQH\nROReEblPRB4QkUPF8kMickJEvlz8vcbb50YR+bqIPCQi15QUZnR/t9xyy9zLYGUfz99Yyz3mso+1\n3LllT1Fqkavq/xORn1bV74nIucCfisidgAIfUNUPBMJ/JfAG4ErgEuAeEblCVZ8uO49hGIbRnEof\nuap+r/i4ApzHRMQBdpj3wHXAUVV9UlWPAw8DV3dQTsMwDCNBpZCLyDNE5D7gJHCXqn6hWPUWEblf\nRD4sIhcUy14AnPB2P8HEMl8I1tfX512ExljZZ89Yyw3jLftYyw3tyi5lfpepDUXOBz4JvAV4rPgD\neDdwsar+goj8BvB5Vf1osc+HgD9U1T8IjqW55zUMwzAmiAga6ezMjlpR1SdE5HPAtar6a96BPwTc\nUXz9JnCZt9ulxbIdHDp0aOvz+vr6qJ+khmEYfbC5ucnm5mbldqUWuYg8D3hKVb8jIj8IfBZ4H/Al\nVX2k2ObtwFWquq/o7Nxg4he/BLgHeGFofptFbhjGmNh16vTW5zOru5PL+qapRX4xcLuInMPEn/4x\nVf1DEfkdEXkpk47PbwC/BKCqD4rIx4EHgaeAXzbFNgwjhi+EbehCRKuEOtzWP+f+1cPcXHw+fGr/\nzETdJ9tH3ulJzSI3jKWjK+EO8YUzdY5QXEPhzi3b/tXD3HrDe+G2bf06JtsG8lrPupayyE3IDcPo\nhFAMy8SzLVVWcw6Pnf/85Lrbzrlh6vsN379t6/Pxc89WHntNldPf3wXA7nPOTK0LLfo6mJAvEa5y\nz+MVz1hO6lrCqeWPnf98du/fFsrTh1d4/hOPbX1PHu+950fPf/rwSrLMOYLcJV1Y6wsj5OFrjP/d\nX77MxF4bTdSNvimzjs9ev7P+rRzZuf2fPXdakKva8s3cBMDPy3tzijg3utKk1uGH88IJdepCXPX4\nE4BZnz7htfBfQ+06GWWUWcypurPr1GnOXr8b3cjzNfsCvn/1MDAtxGFbl33T+5157/kcKwKgf77y\nbMvBaCzypOV9Q7H8tuW2wg2jLbLPE8saD3wntAC6Ub5tKPShBQ5xIQ+td1/Mh0zX3oHRu1ZCIV92\n90kTog3CLHSjoI4g1z1uzI0C9VwpOQ+BoXHV40902sZG61pxmHDnsf2au/0dtn2Uoa9SiDcyE/jF\nRPZNi7QvjitH8u971WCYsB7GeOz853Pc+17VxsMwwU/oQW74/m0z77TMZe0dMxwoNBaLfOjMs9M1\nNwQr1uHkmKWYm69+OJSFDIai71PlhvGte5iuXznulLbE2uM8mJVrZTQW+ZCZh9unKkIg9SqbWj5L\nTMRnjxPWUJjL7kWVe2XLMMhww7htV46cZs9TK1Mhhn0QnmNe/vRjIjPRA7PIWzLLUV0+dazwHPE2\ncV1MyizjKtdH2/N1fewmuHYyL3/62jvoNBBj9J2dg+MGmXrK9yXiXQt2F8SGRNuDYJiEQp6izIWS\n2iZ0vVRtm6LPDvh5CLk/vsWiVgbO6e/v4vi5Zxs9cWN+SV8Q+8pJ0Sfu1bmOFVb1O/uwGJeJlHjq\nRnpdaLGH+6WOX3afch8m7tz7Vw/z8/LebkZCFr78WLx6H+x5anskah/GTUrIK2cIMuKEuRhyKRPx\n2Hp/uzEQNlr/e26DdszDn7/r1OlRPkhj6MZOga367ov3ypHTW39lx3fHkH3bf2VlSuHO7ec1uZmb\ntkZvtuFW3sPhU/v5hB7kE3qw9fFirKlOpROYJQtjkc9FAG8QTh9e2ZEUJ8WiCEQOoSB0dawYXVrt\nXbqKUve7ydtL1/idnzHhDe+Zfw9i5Q6Pkdqm7G3A8c6Nm7YyDG65KN4Buw4+0egtLfa228eAIvcG\n0ae7caFdK2UCOSRLNlbOsbpSZkmVkEO+KJaF1HVFm/vphKpuOV3cdt2RmbGyuuOkrntTIU9tG8P9\nDj9DYVOXRSisfYq5c630pTsLJeRDSkgfo+ztwES7GVVinrL6Q0syzNsR0qZOdHlv/XC9soRTbeqU\n2zfnQZl6QIRhjWV++XCfcHnZ/v79anKP3NiFMFti10L+CT3Y6+QSCyHkfYpgTrrNXacmnSaHT+3v\n/Nwm8DvJERifmDj7gtjk3LlWcZllW3WOsg7GcH0fpETHH/STuh51o2LaWu5NBd2/P64ztQvW3gE3\n33aQW3lPJ8erYvRCPiuhi2UKDCuBoytBNyHfSa6Iz6JDtI2bYoyE/Rs5I4Jz3V9VQp6y8Mvuc879\n8e9N21DEeaYLGa2Qdz2rSK4/PSXeVRw+tX9r+64t97EQs0TrxLlXWaVtOlKbWug+qdwibt3YBb3u\nm1AbcizyLga0+ffksfOfXzs/y56ntieoyA1u6AMTco+yhhijjpD7+KLuvlex7A+BMdDnFGZ9kTOv\nJfQv4itH8vPslGVN3NFxWUPIoZ5VvqbaKD+QO2dMP5r60Rc610qfDauJiLcRYRPwtPXehTUdw/V7\nuP9lb2a5vtA++lLakNsmyq57W3QDdp1Ki2+YiKvL1BLhm9JVjz9RKuah+6SJ6DrtuPWGiT/edaxO\nBhHCzdCZb33QA4LqCLIbyJGzj79d7j65jTLcbkiNeUyk3Ct9+cSdSMdmV6pqxG4b/88/bo4x0PSt\nbxb4A4Nyr3/OtqmO0pwHRxcDt/oaGORwIn3s1+LRMV12kJYKuYj8gIjcKyL3icgDInKoWL4qIneL\nyNdE5C4RucDb50YR+bqIPCQi19QtUJ0b5BpN0xvqv/qEDS72PSXKZWLdVwMdcsPPJdVg+xTsppxZ\n3V274YX1p4xZP/D7qD+h2J9Z3b1j9KfDrfPRjenlOQ/ROsxjTMmtvGdigTPxs6+pwm3aeZRLpY9c\nRJ6pqt8TkXOBPwXeBvwz4K9V9f0ish94jqoeEJErmSS1vAq4BLgHuEJVnw6OucNHXkeMYxMntCGs\n1KFvO1wX2ye2Tdnx27IsvvS6vzPXpRG7f20aV9Uw8lSZ2tzHoblvQtoO3Imta3sO/1ixMMSusxV2\nTWMfuap+r/i4ApwHKPBa4FXF8tuBTeAAcB1wVFWfBI6LyMPA1cDny87RVMSbcvb63bxzI93wykQ6\n9zW5TgOremiUbbfogt5ElKs6mP31XVhGOblAwvuUuudlZe9LuGNhtl0dswmxfERdW+ap33r68ArD\nGQueT45F/gzgS8DfBo6o6o0i8m1VfU6xXoBTqvocEfkN4POq+tFi3YeAO1X194NjTlnkdWe4CYU8\nN6eFa3BVFndXr505Qpxzvtxtlg2/cd/MTVuinGMdd/2anZvYKRTyW3nPVkREWD+7vKept8I+wiWb\nXtuqcnR9z9z5/HDEPU/l506aB42zH6rq06r6UuBS4BUi8mPBemVipScPUbewMXx/as5osjDyYNep\n01lukS59h6HvvU7DrNr28Kn9U9vE/PyxsgydVDnDZWGjvpX3JDPlucx3qX37ouweunW+ePnl7ErE\nq+rEUER8HjhLf/c5Z7ayIg5ZxMuoFUcuIu8Evgf8IrCuqo+IyMXA51T174jIAQBVfV+x/WeAW1T1\n3uA4esstt2x9X19f59VrL4+eMxz8UTdeNGSeYtaltd+2HFXMy/8aPvRir/2xOuAL+KyGS6fwLWvY\nvpbO+m3j13Wk7k9O/cq5r00t9T5z1eSk0Vg0Njc32dzc3Pr+rne9q/6AIBF5HvCUqn5HRH4Q+Czw\nPmAdeFxVDxfifUHQ2Xk1252dLwx7Nqs6O50fPEzoU5XZLafi5foqc6z3ugxFyKG6Mc/K914nNjtF\naIUPTcgdsTj1OsRGG+fexxh93tsuHlZlxy6bNHpRiD2oGo3sFJGXMOnMPIeJG+ZjqvqrIrIKfBz4\nEeA48HpV/U6xz0HgzcBTwNtU9bOR49Ya2Vk1U0mY2ayKsEOpju96Vv71vklZc2GKgb4t8zYNMOWX\nnreQh1S9TbQ5Xgr/vtXtgym737ntrC8/eR/nHCruWvgBHqMaou8nzgmFPLejs4zcjsOmozqHKuZV\nwtx0XRO6aHRjEPK+X/9v5qbGLpYqYsdt4nJpmqWwLosk5Kk3jtEIeWi9VGVLa3rjY9anvy61rEtm\nLfpVscxl23RN20YXCljXIYVdMCsfbllOD2hez/pI11yFCfk2sTe50Qh5mG+hqjHUufEpge5bsGOk\nzhlzbXR5znkPIumqsTkhD3/PojXmKppm6cyhr3pSdo/aRtEs0v2PWeWjEXLYmZM4tg7qJ0+q8o3P\n0kKuM1JxFuedVfrVLoXcEUa4LCtdTFKcom9Rn1WW06ETWuHhdTn73POHn/1wu9Dpm7tyJP94YRa9\nNj38dehjCH5V9EHbss8qFWsf0QZdRL4sAn1HYvWRaqJPmoZ6zos2c/oOJvuhX9g6lvbZ63dv/aXW\nxxhqhyR011HVpKHNclIBR2ogTxXu95mI16NpXpfYsq6EckyCO0tyDaxBCLns2ykgsRlmckWmzrZD\nxIlwrgUVWu9+Q617rD4oO2fTNwEn3ibi2zgxnGWn9s3cxJnV3Y3rVZ9vgkM21rpmrj7yMCKlyhJP\niX3ujOM+fUWhtKWOS6Wrc6amZoPZTPzrrPG6ojy2V+dZEBsBGpIaqJTbJurEnVfR98CzRagfU5Ng\nJ3zkcxPylcefyNo2FWVRZsHnzg/ZV6xtF5QN6Gh6zqq5NHNmfW8q7OHv8Tsomwq5EadKzMtCbeuO\nsSiLviqjT3+7E28/kdqisHBTvZUJSq7YtB34k0vsYVQnXWmZtZRrObvtQuH2ybluTS31ebp2lhV3\nr8P6UxaCmmuZx/ZvKsZN9quytPuM4Bkig/CRl+EqVp8CkNMxWLeyhdkI/f9N3CV1zu93AMf6C5r0\nIXQVXZL6vmiW0zxxWf10Y3Kvz6zu3mEYhIT1tGy73D6bMpqO3XC/LbeDvGlH+hAp608YvGsF6vvR\nYom2+qTtQ6aqc6rKYod6bo8yq7qJ+yRnn5Svtk1GQKOacMBQmYiXta9cV0jqQZCKF2/jjkk9/FOp\njMfI6OLI6yT5qcJ1es5z3semvyUV5uUfN3z9ffe+92T7uSHfhVLlL2+La9wm4rOhKlFamag2HasQ\nxkGH39uMgUj5wP1lY7fGcwdMDdIir8q+V2Zxu+1DX7BuTI8KbTLnZzhZRaoCullfYr8nXBbbx+8A\nDKnKCVOHssmPm+yb2i92ncZqIY2NWP3220ZstHNOfcqp01UJ2MreOGPERH+Z6tGuU6eTFvlcfeQ5\nvtpYZXCVMCYc/vbu+E60m8716fyMTYn54MNlLsvj2evTI7ncA8pdM3eM8DqGn7vwkYf7Vy0vO/4y\nNb4h0sQNFuK/cab85mWCnGpPda1z5wOf1cjkeVKmQXN1rfgVyn86h2Ls0tmm9i07djoyI39wkf8A\niL0Wtsn1kWP9+Nv4sfIxgfa3CUMy+8Y/X+r+xER8EcPEFoWqCBe3jb8s1qnt94eAGwNQDOparTcp\nDJA9P+siUTbF5WB85CkmgrC7cVKn0Arv2h8bq+RNrYNs18TG9v+Y2ykl9E3IfSCk4vjDRn7ravw8\nJubDIhVhVGZ4pEQctn3j4XGcEO8P6kVoJPn7+X5j/2Gw6KwcOc3Zo/F1c/WRhxZkyvKFaad/rh/P\n36+K8MkeeyuA6VFWMcpGmbZ5pS0biFFWjrbUFfEyqubaBHO7dE0XswnVHTORIhR6N2F2atvwuMtQ\nN6ru1yBHdobiHVI2RVYf002lEvTHfHruNaeO+yJ0JaU6b2KVO2d0qzvHLPPMlIl43YfoMjTUWRPO\nhds0TNSnzEApE/KwbudGl5T51BeNpkI+t87OHLEJo0TKhL2qk8Tt7/5S/qaYrz517BxLtEzoY8eN\nTdzrd2hWdVSWXtejxyrLG7Jy5HTyd3Yh4jBp0Cbi/dMmvQK07ySPiXXVgJ2+UmAMnVgAQ9m1H+zE\nEmVTuuXGVub41ctyTYSVaGo268gEF3VC8srK4MhqNE6c967lbZuxXdnvCbdJUWdar2WxtuZBKvwQ\n4v0aYeivwxkTKfpITLdMlrjDn3AZIm/YR0eSa2XyQ3bevCpRjrleUtuHlRb6m7ey7HW2Kn/KFKE1\n7QQ5IcxtXCxVr+CTKKLpZX5kQh0m18Cs8XmT+3ZZp061GewD9oCH/L6uwQl5U3KGeTcRthy/Yu5x\nczozp5++CVdIhVUdi7M/e/3uPKvd3z6Ce1Py+zfaxPCWRbIY7akb7VU1iCdFqp3ELPsy/DeCZasX\nk7f8Zg+vUiEXkcuA3wF+GFDgv6rqfxaRQ8C/Bh4rNj2oqncW+9wIvBn4PvBWVb2rToHOrO7eDq+L\nrQu+1xWR2NRvfiU8s7obVrdnaH/nxk3RyuePFG1j/YbHrjzW3rW0371wncQaVY67pArf3RW6vtzv\nsLwpw8fVsXCEtBNRP+VDV1TlcXHr24zJGD2Zrs8YpT5yEbkIuEhV7xORZwN/DrwOeD3wXVX9QLD9\nlUxk+CrgEuAe4ApVfTrYrtRHHiMlEKkO0KZWYll4X27kTI6vvMpK2TpG4AOvyieeOl9V+XKJRRk1\nmYfT/OOzZUfypUjkSVVu8jI/eZfCv8x1omzQD5D0kZdGrajqI6p6X/H5b4CvMhFogB0HA64Djqrq\nk6p6HHgYuDp27DCKJPaXQyrxUpPKcPb63VPRKmGlDcuX8rHXtYZjLpetaJG9a6w8vicZPZIaTZkT\nClkWkZJC9u2sbC61qL9Nijr31ugO3wLfv3o46QapSncbY57J6RYN3UinEikL1c72kYvIHuBlwOeB\nVwJvEZF/CXwReIeqfgd4QbHecYJt4W9FmTCn1pXt4wb25OQN8UmlEijLCRN+r9MBVGaB52wf0kWM\neVXOmtR6E/D54btTDh8pj0CJkbLGTcT7oW5eqKzww8Ktsgn8qqp+SkR+mG3/+LuBi1X1F0TkN4DP\nq+pHi/0+BPyhqv5BcDxl7/R5Y4LVNMlVDpWvMAG5YYS5mdzKqDP4pwltxTx2X/wRr1UuMMcyv0LP\ng6o6H6ZCrkrP0JeIW70oMTQTrpVKIReR84D/Adypqr8eWb8HuENVXyIiBwBU9X3Fus8At6jqvcE+\nyo/dsr3gh9fhwvUd5+5TyF22wa6oG0dbJfhhigD/PCmrvMsImzJy70vYr2EdofMlV8gdbcYRNMXq\nx4Ste3VyEx7d3F7xwLvq+8hFRIAPAw/6Ii4iF3ub/RPgK8XnTwNvFJEVEbkceBHwhejBX3IIHvin\nk/8REZ/6MR0zxFd8v9fe99HndmD666veApr4xn1y70sX/RZGd5T5X6F+fTP6YWqw4e+9fKKR7i+1\nT0XUyk8C/xM4xiT8EOAgsBd4abHsG8AvqerJYp+DTMIPnwLepqqfjRxX4f7tUBs/XjoRfpOqgLHJ\nIqoIR09V0SQdbK5lnJoYI3W+No2rzGLP+V3hRB3GuMl5INeNiKqLf9xlr1MiGeGHTV0rfVAq5JAt\n5jsq4tFjqFbHYTa19H23hv+/bHvIT6YVE9pUeGHqGFA9eCNGnYRfZlkvBrlC3mfIoRkF22Sl/Rj8\nEP2UqHtUWt9719JPtTo5STIIR0/GKnuu9Z6yzsNlqTeD0KXSJMVt2GDLGqn5upeHsrTMTY/XxeC0\nRSYt4mltnOtUbxw9Fi9cRpa+Umsitn9uUqkS3IVtmwWu6vg+MX93leDGpt7yy5wKi2waW26Mk7KR\nuiFdGCpGHjv0JaWTHsOxyB1NLOaqh0HZMcNhsYlhsv7EEo6+xLyukLp9yibCNYwY4aTkbYlGKW3A\nVg6R4nNfgQxjJXk9MlNPz9dHnksH1nRjinPHZgjKyv2dKrt7YFTkV6jTOemGWYeC/u595dkFm1jg\n5lpZHOqIapULrk69cOdddh95tojvXRugj7zMj12XGvupFn50V4ZUR2uwXGS7I7UzSzyj3KlGkxLf\n0D+e8peb+8QICRPBQXqC7xB7sDcjFrCRZLA+8hxyrPGK3Nw+viA3KUdVLO4Owosf+rv8B0eGjz72\nF9vOMLIp6p0/oUtZB3zXLLWbxbX5jPZfxvxcK26Ifk7hQ4Euc0dkHG/KKvfPUbZv4GIJiVbGo4Hl\nX/YGMGPaNFCzvhYTX8hzXIhtQ1KbjAFZJLZ+fy0B//FhzdlZ6wnkto35ncN1meywyssENiM0svOK\n2PIJXUVfkTfGeAmjWPquH8so3o6u30KGF7WSQ0rgcnt4m7pXIvs6y151bWfF3Fib+rzj5pVZ41Vv\nCHPCTeZhVvliEpsvt8sOzpBlFvMuGUfUSg84MU52fIZWf+bx6jAU32COmyWc2s2EfDkI62isrlhd\nqMcOt24t4q6VcVrkHRL1lzfAt8yztp2Tbzykrogby40l1eqQDJdtLktrkUNglde0wHOOW0ZOStEd\n071BcrBSH/nVwawtI123yvLPG3GmjMZGmjO0zs4FpZaLJXLz3BD5qo4mFwa5ZSn7x6rZURpLNxtO\n32YsL6m3sdzO0Nj0gEZAyzf0pRbyqEtlVh2MifNEE+Vk3OToA6SBmJt4GyG7TpXPF1uW37+ugDvR\nX0rhbyHmSy3kDtW1TsL96ljjqmvJ7afEPBJDn+WvzqwUvnCbiBsxZjk4aNGZavO+5rTUHxPygqbh\niK3PmyPK3g0uK+fWscIcLkEFiTXIpbWCjEpSMeX+6OJY3WlSn0KXoauXu06dtvpZwnIL+d5IbHdD\nen0Q1H3lCv3lOa6ZGqkHTPSXj6oBQitHTrPr1PSfbzC0jXpa6MFrHbhzlz780Ed1rdGw2VlZ8znn\n0Q2Qo8UX77Vt5fE9wM685LEBIIYRw68ns3qQW93MY7kt8oKpStl2qH9LUhM75J4nmkMmEea00FaO\n0SuhwKZSPphffTYst5C3eKUp66xsgp+gKGwQTaySZPk8UXcPjTrWleWQNhw73HFHj3FmNS7ovstl\nmd1yXQw+jLHcQu6RW7m6FnAoF8W6gpkqX+wB0SRxlgm4EbLVQenVO39axFTnutEdy+0jDzoBZV/5\nkP15Rba0wR9cVDbQyHzlRle4kca+mIdMxN1ce11RapGLyGUi8jkR+QsReUBE3losXxWRu0XkayJy\nl4hc4O1zo4h8XUQeEpFruijkvAW0Dys8pC9/tX/cuhO6pjDBN6rYcrv4cdJkTpNo1KbKtfIk8HZV\n/bvATwC/IiIvBg4Ad6vqFcAfFd8RkSuBNwBXAtcCHxSR1u6bvvxK0XP5ye5nIOA+fXQMxVKSAjs7\nQQ2jB1x6ZzfozgS8H0pFVlUfUdX7is9/A3wVuAR4LXB7sdntwOuKz9cBR1X1SVU9DjwMXN1DuesT\nTgdXImIpARc51ttDJRan26Xl6yyklM/SMGbFypHT+WMWZmjEzYIdulLHmCrZNttaFpE9wMuAe4EL\nVfVkseokcGHx+QXACW+3E0yEfxhkzu0ZVp5ZDoBJhR92RWpgh7lLjL6p+3Yr+9gatLdIOVgai3mJ\nKzSrs1NEng38PvA2Vf2uyHYWRVXVSVraJP3kyc1NAVnnieeNgmwTW94U3YBdp7r3H4YiHaa9NRE3\nZsnZ63ez60h6cpIqsXZJuvx2Mro63HIim5BKIReR85iI+O+q6qeKxSdF5CJVfURELgYeLZZ/E7jM\n2/3SYlmE3/Q+/z3gqnold1Tk6q5FxsWcVadnVxVT9jFJtBWU24n56BqAMWqmIloa1r1UnDqMJ/Gb\nboCQ0BJf5H/mFPzxxyqPVxW1IsCHgQdV9de9VZ8G3lR8fhPwKW/5G0VkRUQuB14EfCF+9H+z/bf3\nF/JFOLXdAnbcdfYqmbg2XYj4Ir3yGrMhlhDLfQ9ztNR1M46pHma1vwvXYe9vsq2Xcap85K8E/jnw\n0yLy5eLvWuB9wM+JyNeAnym+o6oPAh8HHgTuBH5Zq6Ygis1WX7VtnX0iNPVDzyKCpU7yqtzj9Yl1\nnBpN8Ov52eceB2gs4GFI45jEPGRLm2pGlc1vqre9sz+vo3IGnjm5VULG8LpoQ/aNriiboKIqj0sY\njTXkNhNS68FzVAY41dusZuMJGFMsa1dldXkumu47ZivHGAc54uss1iqLvWldnweduDjnZZH75w17\noWOjv6qGmucORU9S8VCZ9+jStvjXuE7F8RtE3X0Noy4xAQ51IEbOLEZDt9Jj0TgQ6OEgLfIC/wKX\nJbBPCbz/ubHPdsE6S8NMc02m6xqTVWMsLjlWeM6k5UOvz6EOuj83EXqZETWYpFlnVqfDkXyfV5V/\nrJcOt8Q0aWNjkqa2/mjOsNKbNW4MmTCnUKq+7zqVjl8fAqEO5jII10qMslessJMjXBbbJ4tYTHok\nBnss+NfwzOpudp06nSXIsWs/5MpvLAZtLebct84x12WRAbtWmpIj0lmWaEkmwLGKOExXWCfiVRNJ\npBrT0F9LDWOZGY1FXiXaZYJdKfiRGefHLOAhOR2WOULtHgxlYZFDf3U1hkmq/u1fPQzA4VP7Oz3f\nWOvoqCzy2E11Qu13evifuwwpXCQRD7EBPMaYOHxqv4l4BoMU8ipinRk5PvKkb3jBIlb6wK/8rhe9\najswl4xRzazrSJsxFUNldEJe1/KuG5K4iNZ4jrjmWin+JLpVOF/8ojUaYzFYJEEfjZDnRKdUiXzl\n6MSRhxrWYV6jNes8CAxjFixCfRxMHLmPC5WLUTahaytGHGaYQ9k19Yk9MMNOTn9ZGS6/erivYwy5\nZIz+ya2bfROWYUz1crAWeXgRe+mk8yzwRRbxFOE1daPH3GiyFF1X8EWwiIzmDDWPz5jq5GCFHLqz\nupPHWbIixqpeAAARF0lEQVROzroCHHZqlnVy5h6jjDE1HKMbZF/D3Eg94pdlLHVy0EKeijKJhR+W\nYSF327S1pptaT7mCnnLBjKVBGc0Ychsd6huDz6CFHNKZwEKaVAS3z7LlEHGC2kTU216rJmI+Jl+l\nUY+yJHmzJhVQsXLk9NyCA3IZvJD7hBc4vPBNZhWBcTxxu8bvwPQzJc7C+m0izCbmC8xAosWG/FZQ\nxSCjVqpo8gQf802aFbO6RuYmMXxU16ZyAbWlq+OMicFb5HVe5esK/FBe6eaJ813P0uJNXfcwyZex\nPHTZFnNEPGfcidtuDAxeyGO0mkAiOI6xTRNRbyK4frL8WBnaHNsYJ7Pup0q5aVPfh84oXSs+vp/c\nhHl8pMTcRNzom5iY+6GQ4f8hB0WM0iKH7p6Y7mYuY4dnF/Tlkpm1u8cYBvO2hJtMiTgEKoVcRD4i\nIidF5CveskMickJEvlz8vcZbd6OIfF1EHhKRa/oqeJeYNW8Y88dZvLMW85wQ5yFb45Bnkf82cG2w\nTIEPqOrLir87AUTkSuANwJXFPh8UkZlY/X3ffHvVr4cfylgV1mgDfoyQ3M7ILqgy4oYu4pAh5Kr6\nJ8C3I6t2zFIBXAccVdUnVfU48DBwdZsCDsnlYWJTj3BC3Hm/NhvDJxTNWdeZcFL3MYg4tPORv0VE\n7heRD4vIBcWyFwAnvG1OAJe0OEeU1M3t66abgNfDf/hO3ZOSgR8m8sasidW5sRocTYX8N4HLgZcC\n3wJ+rWTb1pOC1rmwY7wJY8e5TmRf8zeoqkmhDaPrtt1qnt+B0Sj8UFUfdZ9F5EPAHcXXbwKXeZte\nWizbwaFDh7Y+r6+vs76+XnrOJhc2to91avZD1gTXhpGBbqSzIvY1ajN1LubsWtnc3GRzc7NyO6ma\nzR5ARPYAd6jqS4rvF6vqt4rPbweuUtV9RWfnBhO/+CXAPcALNTiJiISL0uduYeH5VIUVxaZDCyuN\nhcPtxF2nLAox93O/h/d3LD5Jo398MZ+XATa0Ni8iqOqO/smc8MOjwP8CflRE/kpE3gwcFpFjInI/\n8Crg7QCq+iDwceBB4E7gl7MVO0FWaODRYzssvtDX5Y4Re+qmbpZZ72l8N0r2dSryv4scQyRuoVt/\nhOFITXAyVj92n2RZ5J2ftIZF7hr2jhsXe1XPnCjCCU/Z07bNBMXLQEzEa1nmiXu1cuS0XWdjB357\nDC30PgbxDHU0Z2OLfCh0dZOGNhvJWHFvSslr6d6Saj5w7d4YMXITXHU54ntoIl7G4IXc3ZwdfrKY\nGGR0qLkbbVZfc2TfdKOpbEApQS85vmGU0Xem07EZFIMX8uRTsUUUxJietEOiVnjh3rVWc6K2CWU0\nFpvYDD6pdVXLc88xdAYv5D5TT8mUSGQIfKqjbWt9RECW1YJ3gtqoEzL0hVsIotEhqQE9dfeJMbb2\nPgohdxb0jqdkC4uvjLE9jfuitUW8d83E2+iM3Pl7y8R6Udv2KITckbxBPQm6MaFW5a8S7kxht/TC\nxhRFP8ss0syOzRqHMQl5TADcsnBdzc61kLF1dPRBmC+l9jUpuwc1742JuQHUMtjq1td5THnYJeMR\n8hixG+uWdWilj/Xmzo0e3Ckm5oYjV6RzthuzePuMW8jL8K31QFj8IeIxdGM67HEZCSN7sl5l6wh4\nxbax625ibsRGeqdyKi2qPzzGuIW8zN3iPgdhcLmhh7qx/Wf05G6KucQMowNyrfFFYTRCrhqJgChz\nn8TimE0o+qVO7Hhq28Q9WibryuiGZXqbHo2QV+ILQyAQJgIDJGaNu/sWEfNlapRGnFTiO2NkQh61\nykNKLMIq37gxY/z+i/AhbG9PRoHI5CEf9lnl9GEtSz/X4LMf7tzXa/gZrhY/sZP5u+sxda1zqSvA\n4cjPkqyIsFh+TSMPJ+R1KXsTH2s9Gn32Q0fUqk74W/0baSJen61rXScu392LqoYXWt7B8Zct6sCo\nSWaCvGWxyEcn5MD0a3iFYCzDTZwZdQda5VpRYaRRap2xnJSNF6nJosSNh4zOtbJ1jBoxxWaNN6cq\nwVhlg0rlI/eXV7jJzK2y3Mg+St1uuSzCpCUp10qjyZfHhIn4bEhe5421nQ+DHCvb2+YsRQO2e7m8\n7F1rHbUydhEvY5yulUxMxHumsKSrrnNltFBk4NbUObB7aUxoKuKLXn8WziJf9Bs2BKaEeSPvdVd1\nYpm7/1lYVkvDcfTY9puZsYPRCrkJ9mzwhbdtHH6t/RukVTAWmNQo4LJxI0tUbxbatWJ0g+pap4Op\ndhwrbIxmiRsFlW9vFtUEZAi5iHxERE6KyFe8ZasicreIfE1E7hKRC7x1N4rI10XkIRG5pq+CG+Nm\nS8z9ePIOIhOMJSQh5suULTPHIv9t4Npg2QHgblW9Avij4jsiciXwBuDKYp8PiohZ/UacWBjiDLEJ\nnkdAVWK8Esy14qGqfwJ8O1j8WuD24vPtwOuKz9cBR1X1SVU9DjwMXN1NUY2FxIl5KOpFNMwsGqMJ\n+oixNziguY/8QlU9WXw+CVxYfH4BcMLb7gRwScNzGAvOVs53XZv+vzGbBGdTKQDM17oQLOs8Aq2j\nVlRVRaRsmObsh44aRiYuqZplxhweu06dZuXI5PPZ5x6PjvZ1y5dNuEOaCvlJEblIVR8RkYuBR4vl\n3wQu87a7tFi2g0OHDm19Xl9fZ319vWFRDKMZZ1Z3m0tlLPid4kVK26nlC8rm5iabm5uV22XlWhGR\nPcAdqvqS4vv7gcdV9bCIHAAuUNUDRWfnBhO/+CXAPcALw8QqXeRaMRaPXaemh2D3ZWW58zghX3Zr\nbqiE9QHYmb++YFnuYSrXSqWQi8hR4FXA85j4w28G/jvwceBHgOPA61X1O8X2B4E3A08Bb1PVz0aO\naUJuRPEbbx+N01ngi5BAaRmIvjF5YarLIuCOxkLeU2FMyI0oruH2JeJbPvElE4CxUuX6Wrb7uDAT\nSxhGXfzwQstPPz5sgpFqTMiNQdFX6JjNFrUAWIhoEnOtGAuP87vDYuekXjSinZ0Fy/owNteKYRij\nwX/47sAs8x2YRW4YxiDxxTy0ys0in8YscsMwBklZx7QN5JpmtBNLGIaxuPjD82Msq0WewixywzAG\njblVqjEhNwxj0FgceTUm5IZhjAazxuNY1IphGIMlDEOsGgew6EnQLNeKYRijIxZPvsyDulJCblEr\nhmEMFj+lcekgoSXHfOSGYQwef7i+xZDvxCxywzAGi/N3y77ldafkYD5ywzCMkWBD9A3DMBYUE3LD\nMIyRY0JuGIYxckzIDcMwRo4JuWEYxsgxITcMwxg5reLIReQ4cBr4PvCkql4tIqvAx4C/BRwHXq+q\n32lZTsMwDCNBW4tcgXVVfZmqXl0sOwDcrapXAH9UfDcMwzB6ogvXShic/lrg9uLz7cDrOjiHYRiG\nkaDVyE4R+UvgCSaulf+iqr8lIt9W1ecU6wU45b57+9nITsMwOuWYTNuUawuoMX1lP3ylqn5LRJ4P\n3C0iD/krVVVFZPGupmEYgyEUcFhMES+jlZCr6reK/4+JyCeBq4GTInKRqj4iIhcDj8b2PXTo0Nbn\n9fV11tfX2xTF8PArtqvQx0RYU42uM4yxsugivrm5yebmZuV2jV0rIvJM4BxV/a6IPAu4C3gX8LPA\n46p6WEQOABeo6oFg31G6VmKVxtFn5al73rLty46R2q/OOcKHRdU5DaMpiy7iMTqfIUhELgc+WXw9\nF/ioqv7HIvzw48CPkAg/HJOQ1xVFHydqbSpXm/OPlUVvjEZzqtrDotcdm+qtBrMQz5wKt4wiHrLo\nDdOoz7zejIeACTnbFaALd0TX+GWad1nGwKI3WGObsN0us1U+aiEvE+DcfWOYeC4Oi9x4l5mwc36Z\nRRwGKOT3z/ysxrKw6I15WWhiXC3CvS/73T8OvcSRG8bgqBN9YwyTZXlD7up3mpAbS0PbCCJjgj0o\nu6HLh5UJuWEY2ZSJzxCs6CE+rGdxXUzIjYVnaA27b7ryLcdEMXfA17wY2r2e1bUyITcWlqE16j7o\nSihSCacspUN95vGgs6gVYyEZq+gM2dodG7OuA7O4dxa1YiwNYxRxE/DuWaa3CRNyw+iZ2IA2E+7Z\n0nWu8qHdP3OtGEYES3Rm+B27Q+nkTblWTMgNwzBGQkrIu5iz0zAMw5gjJuSGYRgjx4TcMAxj5JiQ\nG4ZhjBwTcsMwjJFjQm4YhjFyTMgNwzBGjgm5YRjGyDEhNwzDGDm9CLmIXCsiD4nI10Vkfx/nMAzD\nMCZ0LuQicg5wBLgWuBLYKyIv7vo88+DP5l2AFljZZ89Yyw3jLftYyw3tyt6HRX418LCqHlfVJ4H/\nBlzXw3lmzhfnXYAWWNlnz1jLDeMt+1jLDe3K3oeQXwL8lff9RLHMMAzD6IE+hHyxM7gbhmEMjM7T\n2IrITwCHVPXa4vuNwNOqetjbxsTeMAyjATPJRy4i5wL/G/iHwP8BvgDsVdWvdnoiwzAMA+hhqjdV\nfUpErgc+C5wDfNhE3DAMoz/mMkOQYRiG0R0zHdk59IFCIvIRETkpIl/xlq2KyN0i8jURuUtELvDW\n3Vj8lodE5Jr5lBpE5DIR+ZyI/IWIPCAibx1R2X9ARO4VkfuKsh8aS9mLspwjIl8WkTuK72Mp93ER\nOVaU/QvFssGXXUQuEJFPiMhXReRBEXnFSMr9o8W1dn9PiMhbOyu7qs7kj4mb5WFgD3AecB/w4lmd\nP7OMPwW8DPiKt+z9wL8vPu8H3ld8vrL4DecVv+lh4BlzKvdFwEuLz89m0kfx4jGUvSjPM4v/5wKf\nB14xorL/W+CjwKfHUl+K8nwDWA2WDb7swO3Am736cv4Yyh38hmcA3wIu66rssyz83wc+430/AByY\n90WNlHMP00L+EHBh8fki4KHi843Afm+7zwA/Me/yF2X5FPCzYys78Ezgz5kMKht82YFLgXuAnwbu\nGFN9KYT8ucGyQZe9EO2/jCwfdLkj5b0G+JMuyz5L18pYBwpdqKoni88ngQuLzy9g8hscg/g9IrKH\nyVvFvYyk7CLyDBG5j0kZ71LVLzCOsv8n4N8BT3vLxlBumIz3uEdEvigiv1gsG3rZLwceE5HfFpEv\nichvicizGH65Q94IHC0+d1L2WQr56HtVdfJoLPsdc/2NIvJs4PeBt6nqd/11Qy67qj6tqi9lYuG+\nQkR+LFg/uLKLyD8CHlXVLwM74nphmOX2eKWqvgx4DfArIvJT/sqBlv1c4OXAB1X15cD/ZfJmv12o\nYZZ7CxFZAf4x8HvhujZln6WQf5OJT8hxGdNPnKFyUkQuAhCRi4FHi+Xh77m0WDYXROQ8JiL+u6r6\nqWLxKMruUNUngM8Br2b4Zf8HwGtF5BtMrKufEZHfZfjlBkBVv1X8fwz4JBN31tDLfgI4oaouv9Qn\nmAj7IwMvt89rgD8vrjt0dM1nKeRfBF4kInuKp9IbgE/P8PxN+TTwpuLzm5j4n93yN4rIiohcDryI\nyeCnmSMiAnwYeFBVf91bNYayP8/11IvIDwI/B3yVgZddVQ+q6mWqejmTV+U/VtV/MfRyA4jIM0Xk\nh4rPz2Lis/0KAy+7qj4C/JWIXFEs+lngL4A7GHC5A/ay7VaBrq75jJ38r2ESUfEwcOO8Ox0i5TvK\nZDTqWSb+/H8FrDLp0PoacBdwgbf9weK3PAS8eo7l/kkmftr7gC8Xf9eOpOwvAb4E3M9ETP5DsXzw\nZffK8yq2o1YGX24mvub7ir8HXFscSdl/nEnG1/uBP2DSATr4chdleRbw18APecs6KbsNCDIMwxg5\nNtWbYRjGyDEhNwzDGDkm5IZhGCPHhNwwDGPkmJAbhmGMHBNywzCMkWNCbhiGMXJMyA3DMEbO/wfA\n2s/4AWDW9AAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 16 }, { "cell_type": "markdown", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 32, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "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": { "internals": { "frag_helper": "fragment_end", "frag_number": 32, "slide_type": "subslide" }, "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] # time dim name\n", "print('name of time dimension = %s' % timedim)\n", "times = gfs.variables[timedim] # time coord var\n", "print('units = %s, values = %s' % (times.units, times[:]))" ], "language": "python", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 34 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "name of time dimension = time1\n", "units = Hour since 2015-03-10T12:00:00Z, values = [ 0. 3. 6. 9. 12. 15. 18. 21. 24. 27. 30. 33.\n", " 36. 39. 42. 45. 48. 51. 54. 57. 60. 63. 66. 69.\n", " 72. 75. 78. 81. 84. 87. 90. 93. 96. 99. 102. 105.\n", " 108. 111. 114. 117. 120. 123. 126. 129. 132. 135. 138. 141.\n", " 144. 147. 150. 153. 156. 159. 162. 165. 168. 171. 174. 177.\n", " 180. 183. 186. 189. 192. 195. 198. 201. 204. 207. 210. 213.\n", " 216. 219. 222. 225. 228. 231. 234. 237. 240. 252. 264. 276.\n", " 288. 300. 312. 324. 336. 348. 360. 372. 384.]\n" ] } ], "prompt_number": 17 }, { "cell_type": "code", "collapsed": false, "input": [ "dates = num2date(times[:], times.units)\n", "print([date.strftime('%Y-%m-%d %H:%M:%S') for date in dates[:10]]) # print only first ten..." ], "language": "python", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 35, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "['2015-03-10 12:00:00', '2015-03-10 15:00:00', '2015-03-10 18:00:00', '2015-03-10 21:00:00', '2015-03-11 00:00:00', '2015-03-11 03:00:00', '2015-03-11 06:00:00', '2015-03-11 09:00:00', '2015-03-11 12:00:00', '2015-03-11 15:00:00']\n" ] } ], "prompt_number": 18 }, { "cell_type": "markdown", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 35, "slide_type": "subslide" }, "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('index = %s, date = %s' % (ntime, dates[ntime]))" ], "language": "python", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 37 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "2015-03-13 14:47:13.043237\n", "index = 25, date = 2015-03-13 15:00:00\n" ] } ], "prompt_number": 19 }, { "cell_type": "markdown", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 38 }, "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' % \\\n", " (dates[ntime],fcst_temp,sfctmp.units))" ], "language": "python", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 39, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Boulder forecast valid at 2015-03-13 15:00:00 UTC = 292.7 K\n" ] } ], "prompt_number": 20 }, { "cell_type": "markdown", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 39, "slide_type": "subslide" }, "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": { "internals": { "frag_helper": "fragment_end", "frag_number": 41 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "-rw-r--r-- 1 jwhitaker climate 8985332 Mar 10 09:57 data/prmsl.2000.nc\r\n", "-rw-r--r-- 1 jwhitaker climate 8968789 Mar 10 09:57 data/prmsl.2001.nc\r\n", "-rw-r--r-- 1 jwhitaker climate 8972796 Mar 10 09:57 data/prmsl.2002.nc\r\n", "-rw-r--r-- 1 jwhitaker climate 8974435 Mar 10 09:57 data/prmsl.2003.nc\r\n", "-rw-r--r-- 1 jwhitaker climate 8997438 Mar 10 09:57 data/prmsl.2004.nc\r\n", "-rw-r--r-- 1 jwhitaker climate 8976678 Mar 10 09:57 data/prmsl.2005.nc\r\n", "-rw-r--r-- 1 jwhitaker climate 8969714 Mar 10 09:57 data/prmsl.2006.nc\r\n", "-rw-r--r-- 1 jwhitaker climate 8974360 Mar 10 09:57 data/prmsl.2007.nc\r\n", "-rw-r--r-- 1 jwhitaker climate 8994260 Mar 10 09:57 data/prmsl.2008.nc\r\n", "-rw-r--r-- 1 jwhitaker climate 8974678 Mar 10 09:57 data/prmsl.2009.nc\r\n", "-rw-r--r-- 1 jwhitaker climate 8970732 Mar 10 09:57 data/prmsl.2010.nc\r\n", "-rw-r--r-- 1 jwhitaker climate 8976285 Mar 10 09:57 data/prmsl.2011.nc\r\n" ] } ], "prompt_number": 21 }, { "cell_type": "markdown", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 42 }, "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 = %s' % dates[0])\n", "print('ending date = %s'% dates[-1])\n", "prmsl = mf.variables['prmsl']\n", "print('times shape = %s' % times.shape)\n", "print('prmsl dimensions = %s, prmsl shape = %s' %\\\n", " (prmsl.dimensions, prmsl.shape))" ], "language": "python", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 43, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "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", "times shape = 4383\n", "prmsl dimensions = (u'time', u'lat', u'lon'), prmsl shape = (4383, 91, 180)\n" ] } ], "prompt_number": 22 }, { "cell_type": "markdown", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 43, "slide_type": "subslide" }, "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": { "internals": { "frag_helper": "fragment_end", "frag_number": 45 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "prompt_number": 23 }, { "cell_type": "markdown", "metadata": { "internals": { "frag_helper": "fragment_end", "frag_number": 45, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "-" } }, "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": {} } ] }