{ "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." ] }, { "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", "execution_count": 5, "metadata": { "collapsed": false, "internals": { "frag_number": 2, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "import netCDF4\n", "import numpy as np" ] }, { "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", "execution_count": 6, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 4, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "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" ] } ], "source": [ "f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc')\n", "print(f) " ] }, { "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", "execution_count": 7, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 6, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "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" ] } ], "source": [ "print(f.variables.keys()) # get all variable names\n", "temp = f.variables['temperature'] # temperature variable\n", "print(temp) " ] }, { "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", "execution_count": 8, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 8 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "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" ] } ], "source": [ "for d in f.dimensions.items():\n", " print(d)" ] }, { "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", "execution_count": 9, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 10 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "(u'MT', u'Depth', u'Y', u'X')" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "temp.dimensions" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 11, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "(1, 10, 850, 712)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "temp.shape" ] }, { "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", "execution_count": 11, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 13, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "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" ] } ], "source": [ "mt = f.variables['MT']\n", "depth = f.variables['Depth']\n", "x,y = f.variables['X'], f.variables['Y']\n", "print(mt)\n", "print(x) " ] }, { "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", "execution_count": 12, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 15 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 41023.25]\n" ] } ], "source": [ "time = mt[:] # Reads the netCDF variable MT, array of one element\n", "print(time) " ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 16 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0. 100. 200. 400. 700. 1000. 2000. 3000. 4000. 5000.]\n" ] } ], "source": [ "dpth = depth[:] # examine depth array\n", "print(dpth) " ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 17, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "shape of temp variable: (1, 10, 850, 712)\n", "shape of temp slice: (6, 425, 356)\n" ] } ], "source": [ "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))" ] }, { "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", "execution_count": 15, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 19 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "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" ] } ], "source": [ "lat, lon = f.variables['Latitude'], f.variables['Longitude']\n", "print(lat)" ] }, { "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", "execution_count": 16, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 20, "slide_type": "subslide" }, "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "# 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)" ] }, { "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", "execution_count": 17, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 25, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 6.4631 degC\n", "32.6572 psu\n" ] } ], "source": [ "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))" ] }, { "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", "execution_count": 19, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 28, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "http://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p5deg/GFS_Global_0p5deg_20150711_0600.grib2/GC\n" ] } ], "source": [ "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 " ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 28, "slide_helper": "subslide_end", "slide_type": "subslide" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "float32 Temperature_surface(time2, lat, lon)\n", " long_name: Temperature @ Ground or water surface\n", " units: K\n", " abbreviation: TMP\n", " missing_value: nan\n", " grid_mapping: LatLon_Projection\n", " coordinates: reftime time2 lat lon \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 time2(time2)\n", " units: Hour since 2015-07-11T06: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" ] } ], "source": [ "# 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])" ] }, { "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", "execution_count": 21, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 31 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "shape=(361, 720), type=, missing_value=nan\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD7CAYAAAB37B+tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+MJkeZ378Ptmf5dbPLgM8/N7dWwCc23NhGwr6EO7Fc\nwNhSgi9ShNlVyAlOEQpZQAQns2sHZ+OLHQb5yEm3AkXHDzmEmeCYHwIFDtvEm3CKzE/bY1gc2zpW\nYn322vHCDhd0s177yR/dNVNvvVXV1d3V3VX9Ph9pNO/bb/+o7q769tNPPfUUMTMEQRCEfHnR0AUQ\nBEEQ2iFCLgiCkDki5IIgCJkjQi4IgpA5IuSCIAiZI0IuCIKQOWcPcVAikphHQRCEBjAz2RY6/wC8\nGMB3ADwI4EcADpXLDwE4DuCB8u9abZuDAB4D8AiAqx37Zd9xU/1T55/jn5Rdyj0LZc+13KFld2mn\n1yJn5r8hojcz86+I6GwAf0FE3wDAAD7OzB/X1yei3QCuB7AbwEUA7iWiS5n5hZAnjSAIglCfSh85\nM/+q/DgH4BwUIg4A0+Y9cB2AVWZ+jpmPAXgcwJURyikIgiA4qBRyInoRET0I4ASAu5n5u+VP7yei\nh4jo00S0o1x2IQqXi+I4Cst8LBwZugAtODJ0AVpwZOgCNOTI0AVowZGhC9CQI0MXoAVHmm4YYpG/\nwMyXA7gYwFVE9HcAfBLAJQAuB/AkgD/27aJp4VKDmY8MXYamSNn7J9dyA/mWPddyA+3KHhy1wsyn\niOg+ANcw86ZwE9GnAHyt/PoEgJ3aZheXy6YgokPa1yM53wBBEIQuIKI9APZUrlf2hLp28ioAZ5j5\nF0T0EgDfBPBRAD9k5qfKdT4E4A3MvK/s7FxB4Re/CMC9AF7NxkGIiNkWQiMIgpAg206ub2rYxsI8\nuZZ1jUs7qyzyCwDcQURnoXDDfIGZv05E/5mILkfhNvkpgPcCADMfJaI7ARwFcAbA+0wRFwRBAADa\nB547vN56PzFEtEqozXX1Yy4tLOPm0oW8fHKpN1HX8VrknR1ULHJBmBlo31Y/WQzhNtGF03UsU1xN\n4XaJtsnSwjJuueE24PYt/VrTBjgudqxrLu0UIR8ZekXmFWuIqCB0gl73gK36p0Ty9P75aEKuhLnN\nQ+KZ7ec6f7v9rBs2P9/w/O0Tvx07+3TlvheZ6WbcxDc8fzvmz9rYbIdt3TEi5DOCaVkM8ZonzB5K\nUE0xNeufWs80MlS9fWb7uZhfKoRyfXkOAHDuqWec+9O3A7C5rULtw0aIIMckhrU+GiE3X2PWLHlb\nun69SRnXK6IIutA1Pp/36f3zU8ts637vldsnvle15ZtxEwPAP6bbwgs6ALE0KVshV0KtLoQp5Eq4\nRKgmMS2k0688Vvywd1FcLkJtaB84pN64LHObkAPAR1ZuAjApxKboTblsLiSs+UauJMjMC7nCaXnf\nUC6/fXatcEGIgcvtUWf7qm1NUX5odXp1m5BPuWxu256FmMf2DoxOyGfZfdIUsxFhdQ3Mi3IdhU6x\nWemq47OOK8UUdHPbFHnDs6eieguaxpEngwh3GKZVNSXelvXnDq9PvfqK+2WchLpIYhwHcNejucPr\neGb7uTimLatq47wC2nZyqz7fxTfihudv773TMpTFD/d3rGws8tQZstO1SqxD2PSla4LeVYPvS0yE\nMFyiW3WffGJt1km9foW4U9pia49D8IZnTwGI14eXvWslZYZw+8QQ7ypEbMcD7QP36Urz1c9Tn9s2\nFSYYu4/rZtzEt9yw1YE6pD89ph6IkHdEn6O6dGILuYj2OKmqJzHvu2tA0JCoqLah/Om7zsxNDAhq\nS/Y+8uS4gVh/yncl4n1Y3nWPozfQtpEOQsesrgF7F50/V90/30jhuiGJLqoGEeWGPr5l+7s2wCvd\nH1Ms8oasP7+Nj519uujQqPlaaLNc9AbVl3h3QR1BDzlPeUA0x3V9m9Yxm5Db9ml70NdBRbNEGQlZ\nduYvLSwD6H7g0K4zc9j+rg3MHV7v5IHk0s7KiSUEO2rYsG8IsA2fiNt+19erX8r+Mcsfcm5CN/AK\nyKw3Vd+b7N8n8E1QogsUrpHQhFYuTu+fxy24lZZPLuEuvhF38Y1tdudkkZmULrgGQHXFaCzyIQRw\n/fltfPtZN+AW3Bp0jK6FzBZGOAuk+pDr0z9dlzZvgCFRKm3EXYUmzp+1semiWPwwQH/FzuOHsjkS\nvIMBReoNokt346g7O32VJKVGbitnDFeKHtplG3RhivsYBT/0PvcR+tjqfpaRJURrXCfCpGlHo2sc\ngY8mQu5b14USdKDIRrh8cgmn9883atP6fe9SzHedmcO5p57pzM8/qs7OOg1lCJGv45uMaaWHpvEc\nm4gD9fztXdWJKPdy72KxH/VfZ3Vtcx2FyxBwlUUXbfNcQwTdZWmb1zZU2H2W++n989iOjYnyNWWy\nfPPgC7t5luvpb/skK4u8S9dEVQVVy2wDZ0LwNZKYVrkqW9Ny9kFObwR1LNvOClEReVIH/dqrOuKy\nHl1vkFXr2HA9RNtY7rU61mmt2NfeRXxk5aZonZ6LHwZuvv3GYPdqW7J3rYRUmBgCoVc4s/KFZnTz\nlc01HL5LIUhZ1HOgSzEf+t70+VC11fNQC99Ev261QiBX1/AQLmtS/E2GTBeSrWslpsBVCaarUrle\n6WzLbbOgmBaQuY8if0R3jUkEvB0hVmTM8FFXveqCPutG3XDT0OtYq99j7yJ2fW6udn6WXWfqRaf1\nTfJCXgeXYKrlVbkhYuYsUcedO7zuFHDz+9gF13WOuZ97mxBLX0d0zPkt28Z2d1GONmxO9Yb653Lu\nqWfwPYSP9FTzHtTtwPQZgk07bZ3HSt210mTQiC3uVG8cthviO05Ig9L3rzdOU9irlgvD4xO9kIiN\ntuhuOGA89aOucMXulDaF1TdsP4b7ZNvJdT69f6tjVUXIqMgWoP4o1iwHBOk30mXZ6uvqVrVvail9\nPdoH74CDkMbks7yrGmGXjbSLGcubklJZXLgG0FSJhj4wJmR9E73u6P9P759PR8RX16b/GmC2vabF\nafMGra5pVwODFBsL8zR3eH0qxFElDIuap9xnkRPRiwH8TwDbULhh7mLmQ0S0AOALAH4DwDEA72Dm\nX5TbHATwHgDPA/gAM99t2a/TInf5Gn1pVpt2Mrlm9TY7Jc1j2ywml+D36e80yd1l0Qf6/YnZsEJn\njq+qO643vSa06awHYBdvFU1jiayJ5TqIaZnrkWdLC8vO6JWYHZoqnUeMBFqNLHJm/hsAb2bmywFc\nDuAaIroKwAEA9zDzpQC+VX4HEe0GcD2A3QCuAfAJIgqy+s0nrHnzXBMftH0q+ypzEyu7dWNx7KcJ\nqYt41VtWH8ePhRpKHjqk3OYH93WS20Q+5Pp1fn1bWOcmLis7pi+5al+LH44flTJ/1gYtMlPMLIgm\nlZ2dzPyr8uMcgHNQdC68HcCbyuV3ADiCQsyvA7DKzM8BOEZEjwO4EsD9vmOEirFrIEMTfBXcF6Wi\nW+uhQlnVeVV1vLHieuPpki4y7RVW3vTyKmvchms8gHmN6u67ahwDENAOTevbI+AxB1Y1cVf52FiY\nJ9fDdn15rsP4se6o7OwsLeofAvjbAA4z80Ei+jkzv6L8nQCcZOZXENGfArifmT9f/vYpAN9g5i8a\n+5x4PWib68HXuamvH2optbViXA0mNFwxdH8x6KvTtYuHktkpGSJIVYNgmuKqg1VuFeVCLKYxW596\nC43xhucrw8bCfJwxDJqoN528ou/cNHpmROViuYv7G9zThMadncz8QulauRjAVUT0OuN3hj8EKHpY\nTN1IljoZ1Lp4FW1jcbYRv7qv3V27OroKp1Pf1X02j6Ove3r//KD5rtWbnBJpvV6qcsXs5Ky65m0z\nC26ydxHYu5jVZN68AtpYmKdbcCtdtpdx2V5OWsR9BMeRM/MpIroPwNsAnCCi85n5KSK6AMDT5WpP\nANipbXZxuWwKIjqkfX0z9vJ9QeWwPJX1hqkqZtGgQ/YYH1fcuC+Spsl+ba4J/bPveC6hSM2VY7Oy\nbXVgShCNZP5dJ8qy1UHbOtsOb8Uj07754i1iZXKdybJPjoFQhHSQtiGHvPgxswymlFxPh4j2ANhT\nuV5F1MqrAJxh5l8Q0UsAfBPAR8sdP8vMy0R0AMAOZj5QdnauoPCLXwTgXgCvZuMgtteDqjzdQPVr\ncV3rwvfq28WgjLq4ImNCtnMxpBvHLEOsPg89EgEYfoYZXz1sMxDEPE+1Px9D3e82cd4h+26a7TEn\nNrM0avW56RD9CwDcQURnoXDDfIGZv05E9wO4k4j+EGX4IQAw81EiuhPAUQBnALzPFHEXthth3rDN\nSmdYW1Wjrlwi2KTHP4YPvQ5NHiipDjbyhZA2IVWL0ayLurC3ncihSOfQjBh1V2+TsTupQ98CXNke\nxyjmQNi5eYWcmR8G8HrL8pMA3uLY5jYArVKL+V6ZXB2dIYN6Qmlj0aaALv5mpE0qjLXRuYj1ljD1\nkDi8Fa9uw1YX1PcYuI5b5QoT/Oh6Vtwrf9tNbmTnUFaW2QmVyvDoJg3O1olZtS/z3Ju8wVQRsxNP\n79RM9eG67eQ60z7UztFRB3PfVR3W5u91r531DTlgm1TfnlKlbp1JTshduGJJ9cpkVixdnH0i4sqB\nYvs+BE3FvOp8becdep2aEjMm2IzpHzoixWRjYT5q/LMNZbn5xiK0iWWPSdeCPuaHRVX/X5JJs1QS\neFsoU9XNcg2oCInpNbcfgjox6G2OYRsY00dDiCVsesWOOVAsZ2ydoaF0Ff1SRdsR2iH7zgkz+mpq\nbMErtyP5iSVi53H2WdopEtM378sj49pfVw03ZFRhXcwKnpIl3id6ZIPKtgf4O+pjDTSKQZdCru8/\nB1x1Ws/bg1WyCnkyrpXYN7JuxfS5aOruw+bSGaKh1HGTdFG+0Nd65UtuepxZFXGgOPfQ8+96dG0K\n+xkb6uFcdd+SsMiJ1jjWnIQumligoVZLnY7Rqrj1utuF7qfPeOKQa21bp64gN0n2P3ZsGUN1unhg\ntx081JVVPnfYPvdAqkxY3nBcV4dFPqiQ225cH765VCJSTJqO/jRdSGbfgG9fXVnivjLq6+h+wCYz\nsIglN4kpBn260JrQdVscQ/2Y0MnUhBx73ce1+fFSqnxdYqvYIZ20IftRy6t88XWvdZP+B3V/2wi5\nYCdUzLtoX6HC3KVBoSfJG1udymry5Sr/bRIWtcr21rFLCKjuuFQ0DSF0hVw2cdOEjh7cfGCsVK4q\n1MTWflydnF0fuyuqLO1oycAyITmL3BSRkERRg6DnYu5BzIHq1+SmbimbZWbLh+KaTLaJO0ihLPKx\nWU4poI+QtoVrAsOmc2hybFtKalvdGWNUE+0DJx+1Akze2CYDGTon0kwoMbFZ0+afa11FiPWuRlJW\nNb6QdYR+sHUi6g/qLowhV7u1DQSrGx48d3h9KiV1iIjnjDkXsXO9VCzytlEZvWIKus0ij+R6qftK\n3LaTK3QAlbk/vZyhDXQMVlLqmHOHVrkp24q7+RDXJ85Qy0NT8vr2raiqQ2Ppf9m8ZjlY5FU3MmTA\njEs4msyObtsHr4CYF0kl0rcJNa+AXL/59u0qq0uUQyxiff0qzAbl26aOD95G7g0rR0JCQOsIeEgn\nvBIg3bI03xLNuqzXYf3NssnDpU6cfcq0yn7YF3Vvkm193YI0f1cXwZfjvOr4XYcxFf6vwoonLHpf\nk0JFtO51DWnoIdtWvUWMoWHliK8+NHVjmvsMsex1S91X50IHk+n7bHAK2cArIFq1u1iSda2YQhwy\nWMBVecxp34Dp4a+uY5vb+8oRul4IbV43Q7YLPX7d/fpe3wG3T1PEvRtc/uJQ11nbDnRfO7a1F1+9\nM+uR2YE7diEvclBdln4cucvPanba1PHjhdxcr5ivrk0k7woR8tg5Y2yEhiTGILTz0vSv28ro65wS\nMY9PlaFi3lvbd58VXyfaJGQWMNf+9bqlt7Oxi7eOL2olKSH3WbT6zYst5EB1ha8l0qtrhX88pFNU\nI1aUTmwxb/t2oJiVMLGUsA3Zd4Ws1nkLcz3cQ+qKWkfPuBly7FkQ7coHW2rZD5WQ+26OLtpmAw8R\n1Lo33vcaau7Le3xXmKIm5K7oEFW5QwfddBnFExIvbltXjz0POY5Y5N2h6pKtztSJSqqD603RjGcP\nmbBa37ZRYTKjUtdyiFox4RWQ6tE2YyltkR2+ikn7wCoG1YxFdVG7Iq+u+a1wi8CbZdZnvnGVp7GI\n14iDt4m4q1y2WPY6DW8skQUpYovmCnXdNcUWUaViwM11q9piEuHGQ1BzzMpgUSuuhm5OKmG6NHS/\nWIjFqCySucOWYxluGn0dMwSvsT/OEU9edzBESEhgbHyRQDEQS3xYTCvdFhYY0j/S1CVYZUzNkiXe\nljQtcldstgc1rZYeH+57TXGFRbnis6dQAt1gtKcZJ2uWRy9HcHksx5iiweAk2/F9I/TEuk4LVz03\nqeNG6wsR8XCSiCOfQHUSrjhGSzrEyBfC5qq4TUXSWq4AbC4IszzW/aprYjl3n5slZlSL6fc2O6Y3\n9ytJsJKi8EsX98fmXtHddXXfEhUhoYx1mNW6pLwRAGobXV6LnIh2EtF9RPRjIvoREX2gXH6IiI4T\n0QPl37XaNgeJ6DEieoSIrq53KuXoSctcnUDhbtEFZGNhnkJcDjYBnVhf+bZX1+xWgMX3zbxIdYfh\nVz04JsqpH8+w/qv6A2JgK2eVheTygwrDYrYRV3toWq+6qItijdfDG7VCROcDOJ+ZHySilwP4AYDf\nB/AOAL9k5o8b6+9G8Rx9A4CLANwL4FJmfsFYzzv5srUsDh+13itvDvyJ+npoxJOrY1dts4kh9iHD\nm13bh4RqBe2/Aa57oB8zxLUi/vF+8WWuDBXi0Le+tsxynZiwygFLsIR9QJDXtcLMTwF4qvz810T0\nExQCDcD6xLwOwCozPwfgGBE9DuBKAPdPFThS+GAxqGB+SuQ3FuaJsCXyXmxuC8sy2/B+575t7haP\na8jsWAQAHN6lLa/XYKoaaNPOSz01qlq2lVq0+M3n5hKLfVhsYYjmwJ++BpoJ0yhjcUrQS0OSyK6b\nwT5yItoF4AoUovxGAO8non8K4PsAPszMvwBwISZF+zi2hL8VPlF3/VYsn7f62rayiRmdluYT0OGX\n9lZsVwiixwUTGtrnI6Z15POZVllMIuLpEyLMrqil2GMYmvrmx8yUe9nWZ6gRNCCodKscAfDvmfkr\nRPTrAJ4pf/4jABcw8x8S0Z8CuJ+ZP19u9ykAX2fmLxn7Y+ChSt9yl36yKSH3UTfaI2DSidAKu7Sw\njOWTS/WOH0DdRhiSLyX0d996QneE5iuqGnDWScemhvjHHffK41qpDD8konMAfBHAf2HmrwAAMz/N\nJQA+hcJ9AgBPANipbX5xuczCJ4HVfw48fAg4cWS6wH1M4hB6DK0z1LptwEhOk9CImToibu6vTeii\nKxzNtT8VdugaraeW29YT0qDJ4KG2A4dMRMQNThwpNPLhQ8DrvuRcraqzkwDcAeBZZv6QtvwCZn6y\n/PwhAG9g5n1aZ+eV2OrsfDUbB9m0yG0YOUpcESxtmLDGldjWfXDYRmzqYYKeiBbXIIy21N2fz9ry\nve7KQI28qcrnYcOXT6VtndORejWdH2cy+KFB9kMi+h0A/wvAGrC58xsB7AVwebnspwDey8wnym1u\nBPAeAGcAfJCZv2nZr1vIHbgEXe8UCBV9Z0WusqxdfnS1zLa9R9Rto+dcVK3X9qFQZ9SeCPk4qCvo\nsUZ2uvYz63VqU8uMCLXT++c1bUktjW1NIcfexenEVWbPLsLEvDMh9+3D3Beqc3cPRZWlJW6RceBq\nB02jVprmcJl1AVcQrbHSB/cAwQbhh6mhz6ITYqF70YW3qjMzxOq2LXOIuqsjKYZlHgvfMWQiiJFh\nGCM2N0qMXDsSwljB3kV3u6swENMXciWISnDLyhYs2KGE+MirQhWrsDQUnRCRHlrEhfGgT+zAvEim\nb1YRo0NT6lQYzjQdFaQr5LqAd0FI1InGZqC+/jpqE/aWybRyYGlhGcCtQxdDiACvgKrmiK0TQ25G\nLW25TSbHc7SdRWtsTOmK543eRprZDwH3yMiU8ZXP5n5pcT5tsyM2pRBxYVSsrgXHmPsw66HP/SZ+\n8S3cMeMGHhdwuha5i5ZiXgxzremW2bsIojWe8su7OkGryhjhgTSU9a7i2m9ZGOTwQodUzaHpMxhE\nmJuhd3Bu4rLGPbqRjkXeIFd2E6yCHEJZPj3fuflbG3JxqwjjZajJjGfezRLBhZxe+GHTATo1aGqV\n2zIgKir3V5GUKychl4iVcaInRDNTKlTlFmpiHPmmbpwFmgVsNByi3zlaJEpfVjkQf8To5v7U+diE\nWxAywXxYe3OVN2y3syjeithRd8MLucKM+OhQ1BtdxLJs5rZEa6zPM6pPN8e8OPE3sb+K80s15nZj\nYX7KWhPGievNyxT1NoJsdVUKtRnWtRI5kqMOzry/LfdXB1t+85RRDVsGBM0Oroe2bmiIENejneak\n6loZmAmXSAt0y7xy3X1gWwePaYWnFOonwj2b+CZv7mPKQSGM4SzyveVxB/Qdx7bKzf36sI2iazoK\nrkt3hwi4YE4TZ45AbltHZmnav/Zak6JFnlIHYCSffKiLxZXv24VtsEWsvN7mPiRnuKDDK6A2MwS5\n3kAVUtfak59rJWIn6MTTcYCHSsxp3No2BH3ihzb7EcaHbpHbEmr53gbrxogr0Z/52PKaDCfkDWbV\nmdguoqDH8pPX6fBUVkjosGd9OH5MsdX3JSIuVNG1T3zsPvc2gxF9pGeR15l+LSKxxLwuSjyrJo1Q\n+CIE6gqxvt9tJ9dZwgoFG7wCqnKj2OpOk/q0sTBPekiiqpdSP/2kJ+R9snfRns2wAW0HGOkNpa9Z\nxZvuWxrV7BESYugT3aYhiq52MTMEZkGcbSHHpJg1tcpjjBKtFNXVtaDGYLPKzcZgWuJqO3GtCD6a\ndIK3DVHkFWweU+LV3aSXa6UPLEKtKslENrIaecrbsO3kOsealLbKUraFONYV8FkKFxPc1Hkrk7qy\nxeZMZ74BkS5jcpUyn7MzBo4Z7fXJhEOFPGauFr1BtJ3T0BTZ0MYW2tBExAUTVSdUO1LflxaWN9Me\nm8xq/WkcjVMx+fJwrhVdVG1JpuruoyWhF9iaN6UlMRPwu159Y01CMasNUHDjcn24RBxo1hE601Ro\n5LA+8hABd/1uy5jo25fxm813R/swmavc2Ca2gJt01ZkT0nkqDUuIhTzs+8cr5ES0k4juI6IfE9GP\niOgD5fIFIrqHiB4loruJaIe2zUEieoyIHiGiq2MUspaAOn1LbjdJVWdMF1Z4k3LExHxo1HmISEMV\nqlBW+kxGmgxA1VRvzwH4EDM/SEQvB/ADIroHwLsB3MPMHyOiJQAHABwgot0ArgewG8BFAO4lokuZ\n+QXnEQKmRZuYHDZ0wJBnXZ9gbnZEoHyArPQbVx6bjYV52nZ4ujN17AMvhDQo3C1F3aN94Fmvd7Z8\nNbY35roPQK+QM/NTAJ4qP/81Ef0EhUC/HcCbytXuAHAEhZhfB2CVmZ8DcIyIHgdwJYD7nQep4+f2\nuVkCp0uyVSTzojlnAdJmUKk8UE1U4it9GHTM42wNsLD3BRTXRawnoR9C3+oaT82YKLYEdy5BryPm\nwT5yItoF4AoA3wFwHjOfKH86AeC88vOFAI5rmx1HIfztadKxGTjAZ8pXbk4eMaLcD64GJO4SoWvq\nGibbTq6zGrQ3phwsrtQc+rK6FnmQkJdulS8C+CAz/1L/jYv4Rd/F7efC61OsuaZbq7EvZ8XpMLlW\nV2Jqy26oVxQRcaFvqhJtmalz9f9q5Cjtw+b/Pso8FCHuqCofOYjoHBQi/jlm/kq5+AQRnc/MTxHR\nBQCeLpc/AWCntvnF5bJpHj609fnX9wDn7aksbFt8FyT0adjXa14st4oabGTuj1dAtG9+kFnThdkl\nJHe+q536ltO+eZ47vJ6NUVJ1HZRr5ex3/necWf4L4OFt3v15BwQREaHwgT/LzB/Sln+sXLZMRAcA\n7GBm1dm5gsIvfhGAewG8mo2DTEwsgTAHf9NOANs+XEzs2zbyquFs4XXQ5/+M4ZPvegCPXhlzaUTC\n8Oj1Uk0d6JtWztV2zd9sRkvKKEPL1Xc3pXuOkZ1VrpU3AvgnAN5MRA+Uf9cA+CiAtxLRowB+r/wO\nZj4K4E4ARwF8A8D7TBE30U/A5Tdqkrc7yjRU+mTQPYg4MBnqGKOHX8RVSBF94FpVHa37Jp2bq6Xu\nG4iNwad6M29ELatZI2S70FSxmxj+8L57z3MYDp9DGYU8iDkoLaf6WOe8T79ye2JD9EvqukraWKk5\nDk6IVeY2qWclba3QBzHFN6f6GuO8B7PI9adKrSdSgNXu2842MMYplppFnnssa1NftnlvcrJ0hPzo\nUoBTr7vmudu0LlmLHJiebqyPnCO+ZWPATPDf5DxzsmoEAfDX89Trs6mDei72qjzwSQg5UG9ig6bi\nW9kBasaI9zztW0z0a9Sk88dW6VO3aIRxsbSwXHubqrf11DtCm07wUhlHPgRFTuNELvjexazCmRT6\nNVTxtbRvvrimK/5tRcSFITDrnS8NblPGmooiGYu8DjFiyq14Zg7KEV18Qy0R1+tn6q+lwjhpYpXP\nIkl0dtrwDQ4wCe38dGUcs5GzgJuYMxDZzs0cEm3DnHXIZqWrwR0Rii3MEH0bCrnWUZd2JulaCRFx\n20CAWCkyxyTiodTNtiYIuZKriPvI0rViI2RAkC7QUUZ+ZkLIDEFVYZ1mj3poFkVxyQhV9F1Hxjgu\nIkmL3EabGW1s27ks0DFa42bnsc39EdrBXCceXa0r7hYhRcY0KjlJi7yJSIdMLKx3+M2aG8E836rc\nzjEmarZtb8a3C8LQjKE+JtvZaU4LtZX9q9lUbiGiNEZrXMfs9AQmz9kUdnUdzU5OfVkVIdN7jcEi\nEtqRipDqOpGiHri0M0mLHCguorIKJ0RcF28jO6GJ2k5EvCAkUkf/c60fW3jHYBEJzaF9mJpTdmiK\nHOeJjGUccxLpAAARaUlEQVQJIFmLHHDEPttm6FFzdqrPLjzrzIKQA2GhiF0f14dY57OFemNrmkMp\nNrZ5M1PShuwscsAirqYQq+nc6k6/1uF0banTViibWimhQ49tgj+WuRoFO6lY467orhzqXtJC3giX\nSFuWq5uV0hO3D1RisiHOOyQpmi7m206uy3R0I0bd26Gs8Sn3rbZcJ3VjIn0h91nPum+8QYKrXGcU\niYFqQOYM5VUVNoao1tmHsuJFzIVQQi38qvVsKa9TJX0hB6bF3JWl0CXmeiepcsf49jej9FVRVedW\njBBHIX/0wIa+aBrinCrJDwhiXiSiNQ4SWyXYtomTAbfQZ5yuti2T4Yfzg76Z6BPwyiCi2SRWp6cr\nbUebFNgFaYp5Hha5z5oGqvOIh0SyCBPhh30eR7fGmmRsFPInZp3T02/k5B5pQ/IW+SZKjH3Ca4sx\nn2Fruw9oH1p1Rtq21S1zQWhKLNE2B8alSB4WuaKJ9RzokiFaE+FoQFfWu5rqqot9C4kS4e3Y5ToJ\n7pPJ9A29UsiJ6DNEdIKIHtaWHSKi40T0QPl3rfbbQSJ6jIgeIaKro5fY5V5pup4gCEmgJjiP2fFZ\naz9G35pu0adsjQNhFvlnAVxjLGMAH2fmK8q/bwAAEe0GcD2A3eU2nyCieFZ/E193JDEXf209iNaY\naK0IZ1T/HddQrdt3GYUE0dpxXUG3peQIdq849ENNkxhciIGoFFlm/jaAn1t+sp3cdQBWmfk5Zj4G\n4HEAV7Yp4EQDV3HjNQb9VFJD6EXMBaFblFVeB9eAHpegh8xFkINfXKeNtfx+InqIiD5NRDvKZRcC\nOK6tcxzARS2OUeATb9e6dbapEPOQDH7CFlbrus79EGYbo574kt/pYYamQKtcKepPj2Q5vX9+OpNq\nxq7YpkL+SQCXALgcwJMA/tiz7rBWbEQBETG3s+lGsblIQq592YjEvSK4aDOAxzlPb8bCbdJIyJn5\naS4B8ClsuU+eALBTW/XictkUZYep+tvTpBxebKM3Q0MXhTjINRUawLxIrrbqiguPPTI0lVGcRLRH\n10rXeo2EnIgu0L7+IwAqouWrAN5JRHNEdAmA1wD4rm0fzHxI+zvSpBxBBIqJGeq22UlHa5ybv6xP\n6lrR5voT/Q4i/EJJEzG3/eZaVkUqk5Ez8xFdK13rVQ4IIqJVAG8C8Coi+hmAfwtgDxFdjsJt8lMA\n7y0PepSI7gRwFMAZAO/jIRKeA/aZhBwVwxmvvCks4lIx2RRkVxphV/4bLWbf1rElQ/MFhaofRGts\n1iPbUH7faM6mYp7qkHyTpCeWKNaN4Dd1WHq+ASe6pSgW+TSVQq5wpVTwTNMn11kwcUWM6eJdJ69K\nqLinNigty4kloiERE91RdV1Dk52VpPA6KySIx81iRrX4IlnU8pDkXKmJuI98cq00pYE1Lvip/ZZk\nm9nJtk65vG3+FmGEBPafhGZPNMMUzW1zMyiSt8ibDBBwUgpKHZGQkMMtWo/A9I2+NRKeyWhPwUbs\n9ujaX27tPnkhb4XD8qsSCP13W3rVWWJi1iAluF1GlyjLPPMBGkJkSiNgaWG5schWhSrmlFvFJAsh\nb2SVRxKB3J7MMfE+8ERkhYxQ7djlN9fJza0CZCLkQqI0EfOanc6S30bQWT65BKBbAyvH/pn8hDxU\nPOrmW/GQ4xO6LZVD7dtc25rhiiLmgs7SwnL0fW4szJP6i77zHshPyOuixzk3dAfk+IRuTdV0eVXT\n70VGxFz4yMpNmyIeS8xzFm+dvIS8iXhYrL8qn/vm8ODVtZm0xoHAh5d+D0TMhS4x6pdysQgFecWR\ne0YEOlEWuQo9DOw4jRr2OAv4JsCWwVhCW1bXgJW4uxyDJa7IxiLfFFZbDhUbutUuoWzdEeIr18MJ\nG96LWY4eEgp0K3zqTXnGjYW8LPJQRLT7QX9DsvVDtOyb0JlVF5cwzR/tuxVAOhkKUyAbixwIdHd4\nREPcJQ2oM8uS7drLQ1VoCdEazz27C0CFVY7ZfXPLSsitiPtkOKqud8P7YcuDIQiA582srGuhFvqY\n/ONAhkJutaoDBGMmQwhbEu0Npqb/0pa8KGTCXGEGiDg+ZExkJ+RCz8R4y4lgmQszTM3645useSxx\n4yZZCnld61qs8bywWePA+F6HhQYEhrXOWifoOKNWNETE+yF0tqUQZq0RChWsrgGHd7XezZgNgSwt\n8lBExPuh6jq3vQ+n98+PuhEKflTEygSO4AZXR/nY68/ohJxXQOpv6LKMFf0ah15ntZ5rffGHC0EE\ndHLWmbtzLCQ/+bIwPMo1EuvhaHO1uPzip/fPy5vVjLPt5PrkRC+e0d1jrzcu7RQhFwZBF3PdgrLN\nuTj212LBjS7igGFdV+RemiUhr3StENFniOgEET2sLVsgonuI6FEiupuIdmi/HSSix4joESK6Ot4p\nCGPCbGRq+i1xsQhebHnxHcxStswQH/lnAVxjLDsA4B5mvhTAt8rvIKLdAK4HsLvc5hNENDo/vBAH\nU7Rtvs0u/ZvbTq6zafEJ6aFyj5/eP18rdfIYLXIXleGHzPxtItplLH47gDeVn+8AcASFmF8HYJWZ\nnwNwjIgeB3AlgPsjlVcYIfoADj0RUtcNUR2HELcPQBD6pmkc+XnMfKL8fALAeeXnCzEp2scBXNTw\nGMLI2ViYp+kc0/PR8067MDtYhUzRsmzO6sO49YAgZmYi8r2eyqurkDRzh9elQzVBtp1c59Ap3WY9\ns2lTIT9BROcz81NEdAGAp8vlTwDYqa13cblsCiI6pH09wsxHGpZFEBqxsTBP4iNPm+WTS51MtpwL\nRLQHwJ7K9ULCD0sf+deY+bfK7x8D8CwzLxPRAQA7mPlA2dm5gsIvfhGAewG8mo2DSPih4EIJa5dx\nwOoYSsjFGk8T10PWNjvQrFjkbcIPVwH8bwC/SUQ/I6J3A/gogLcS0aMAfq/8DmY+CuBOAEcBfAPA\n+0wRF4QhMcVBRDxd1L3xWeS8ApoVEfcRErWy1/HTWxzr3wbgtjaFEmYXFbnShTWui7gIeJ5Ix7Qd\nifEWkiN2Y9XjxWfZ3zoGZMCYndGnsRXyouvwseWTS2KNC6NDcq0Io0dcKnni6+yc1Xjxxp2dgiAI\nfSOusHqIa0UYPWKF54cvxr/wk0unp45Y5IIgJItrQJAM5JpELHJBEJJCF2mXa0XesiYRi1wQhKSo\nEmkR8WlEyAVBEDJHwg8FQUgOlw981q1xCT8UBCEbNhbmaWNhnpYWlmuFIM7S9G46YpELgpAsNst8\nlq1yl3ZK1IogCMmzsTBPN+OmUtRvHbYwCSIWuSAIySK+8knERy4IQnYoX/nQ5UgdscgFQRAyQSxy\nQRCEkSJCLgiCkDki5IIgCJkjQi4IgpA5IuSCIAiZI0IuCIKQOa1GdhLRMQDrAJ4H8BwzX0lECwC+\nAOA3ABwD8A5m/kXLcgqCIAgO2lrkDGAPM1/BzFeWyw4AuIeZLwXwrfK7IAiC0BExXCtmcPrbAdxR\nfr4DwO9HOIYgCILgoNXITiL6SwCnULhW/hMz/xkR/ZyZX1H+TgBOqu/adjKyUxCEqKwRTYjZ4gg1\npqvsh29k5ieJ6FwA9xDRI/qPzMxkXFxBEISYmAIOjFPEfbQScmZ+svz/DBF9GcCVAE4Q0fnM/BQR\nXQDgadu2RHRI+3qEmY+0KYuwhV6xVYVeI+JFZrL9Jgi5MnYRJ6I9APZUrtfUtUJELwVwFjP/kohe\nBuBuAP8OwFsAPMvMy0R0AMAOZj5gbJula8VWaRRdVp66x/Wt79uHa7s6xzAfFlXHFISmjF3Ebbi0\ns42QXwLgy+XXswF8npn/Qxl+eCeAvwVH+GFOQl5XFHWUqLWpXG2Onytjb4xCc6raw9jrTnQh76Iw\nqdCHeIZUuFkUcZOxN0yhPkO9GaeACDm2KkAMd0Rs9DINXZYcGHuDFbYw2+0sW+VZC7lPgEO3tSHi\nOR7G3HhnGbNzfpZFHEhQyB/q/ajCrDD2xjwrNDGuxnDvfed9GYAu4sgFITnqRN8IaTIrb8ixzlOE\nXJgZ2kYQCQXyoIxDzIeVCLkgCMH4xCcFKzrFh3Uf10WEXBg9qTXsronlW7aJYuiAr6FI7V73da1E\nyIXRklqj7oJYQuFKOCUpHeozxINOolaEUZKr6KRs7eZG33Wgj3snUSvCzJCjiIuAx2eW3iZEyAWh\nY2wD2kS4+yV2rvLU7p+4VgTBgiQ6E/SO3VQ6eV2uFRFyQRCETHAJeYw5OwVBEIQBESEXBEHIHBFy\nQRCEzBEhFwRByBwRckEQhMwRIRcEQcgcEXJBEITMESEXBEHIHBFyQRCEzOlEyInoGiJ6hIgeI6Kl\nLo4hCIIgFEQXciI6C8BhANcA2A1gLxG9NvZxhuB7QxegBVL2/sm13EC+Zc+13EC7sndhkV8J4HFm\nPsbMzwH4rwCu6+A4vfP9oQvQAil7/+RabiDfsudabqBd2bsQ8osA/Ez7frxcJgiCIHRAF0I+eKpH\nQRCEWSJ6Glsi+m0Ah5j5mvL7QQAvMPOyto6IvSAIQgN6yUdORGcD+D8A/j6AvwLwXQB7mfknUQ8k\nCIIgAOhgqjdmPkNE+wF8E8BZAD4tIi4IgtAdg8wQJAiCIMSj15GdqQ8UIqLPENEJInpYW7ZARPcQ\n0aNEdDcR7dB+O1ieyyNEdPUwpQaIaCcR3UdEPyaiHxHRBzIq+4uJ6DtE9GBZ9kO5lL0sy1lE9AAR\nfa38nku5jxHRWln275bLki87Ee0goruI6CdEdJSIrsqk3L9ZXmv1d4qIPhCt7Mzcyx8KN8vjAHYB\nOAfAgwBe29fxA8v4uwCuAPCwtuxjAP51+XkJwEfLz7vLczinPKfHAbxooHKfD+Dy8vPLUfRRvDaH\nspfleWn5/2wA9wO4KqOy/0sAnwfw1VzqS1menwJYMJYlX3YAdwB4j1ZftudQbuMcXgTgSQA7Y5W9\nz8L/XQB/rn0/AODA0BfVUs5dmBTyRwCcV34+H8Aj5eeDAJa09f4cwG8PXf6yLF8B8Jbcyg7gpQB+\ngGJQWfJlB3AxgHsBvBnA13KqL6WQv9JYlnTZS9H+S8vypMttKe/VAL4ds+x9ulZyHSh0HjOfKD+f\nAHBe+flCFOegSOJ8iGgXireK7yCTshPRi4joQRRlvJuZv4s8yv4fAfwrAC9oy3IoN1CM97iXiL5P\nRP+sXJZ62S8B8AwRfZaIfkhEf0ZEL0P65TZ5J4DV8nOUsvcp5Nn3qnLxaPSdx6DnSEQvB/BFAB9k\n5l/qv6VcdmZ+gZkvR2HhXkVErzN+T67sRPQPADzNzA8AmIrrBdIst8YbmfkKANcC+BdE9Lv6j4mW\n/WwArwfwCWZ+PYD/h+LNfqtQaZZ7EyKaA/APAfw387c2Ze9TyJ9A4RNS7MTkEydVThDR+QBARBcA\neLpcbp7PxeWyQSCic1CI+OeY+Svl4izKrmDmUwDuA/A2pF/2vwfg7UT0UxTW1e8R0eeQfrkBAMz8\nZPn/GQBfRuHOSr3sxwEcZ2aVX+ouFML+VOLl1rkWwA/K6w5EuuZ9Cvn3AbyGiHaVT6XrAXy1x+M3\n5asA/qD8/Aco/M9q+TuJaI6ILgHwGhSDn3qHiAjApwEcZeY/0X7KoeyvUj31RPQSAG8F8BMkXnZm\nvpGZdzLzJShelf8HM78r9XIDABG9lIh+rfz8MhQ+24eReNmZ+SkAPyOiS8tFbwHwYwBfQ8LlNtiL\nLbcKEOua9+zkvxZFRMXjAA4O3elgKd8qitGop1H4898NYAFFh9ajAO4GsENb/8byXB4B8LYBy/07\nKPy0DwJ4oPy7JpOy/xaAHwJ4CIWY/JtyefJl18rzJmxFrSRfbhS+5gfLvx+ptphJ2S9DkfH1IQBf\nQtEBmny5y7K8DMD/BfBr2rIoZZcBQYIgCJkjU70JgiBkjgi5IAhC5oiQC4IgZI4IuSAIQuaIkAuC\nIGSOCLkgCELmiJALgiBkjgi5IAhC5vx/oWJ9OHx0YTwAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "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)" ] }, { "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", "execution_count": 22, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 34 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "name of time dimension = time2\n", "units = Hour since 2015-07-11T06: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" ] } ], "source": [ "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[:]))" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 35, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['2015-07-11 06:00:00', '2015-07-11 09:00:00', '2015-07-11 12:00:00', '2015-07-11 15:00:00', '2015-07-11 18:00:00', '2015-07-11 21:00:00', '2015-07-12 00:00:00', '2015-07-12 03:00:00', '2015-07-12 06:00:00', '2015-07-12 09:00:00']\n" ] } ], "source": [ "dates = num2date(times[:], times.units)\n", "print([date.strftime('%Y-%m-%d %H:%M:%S') for date in dates[:10]]) # print only first ten..." ] }, { "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", "execution_count": 24, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 37 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2015-07-14 07:22:39.579246\n", "index = 24, date = 2015-07-14 06:00:00\n" ] } ], "source": [ "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]))" ] }, { "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", "execution_count": 25, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 39, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Boulder forecast valid at 2015-07-14 06:00:00 UTC = 296.8 K\n" ] } ], "source": [ "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))" ] }, { "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", "execution_count": 26, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 41 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-rw-r--r-- 1 jwhitaker staff 8985332 Jul 10 06:43 data/prmsl.2000.nc\r\n", "-rw-r--r-- 1 jwhitaker staff 8968789 Jul 10 06:43 data/prmsl.2001.nc\r\n", "-rw-r--r-- 1 jwhitaker staff 8972796 Jul 10 06:43 data/prmsl.2002.nc\r\n", "-rw-r--r-- 1 jwhitaker staff 8974435 Jul 10 06:43 data/prmsl.2003.nc\r\n", "-rw-r--r-- 1 jwhitaker staff 8997438 Jul 10 06:43 data/prmsl.2004.nc\r\n", "-rw-r--r-- 1 jwhitaker staff 8976678 Jul 10 06:43 data/prmsl.2005.nc\r\n", "-rw-r--r-- 1 jwhitaker staff 8969714 Jul 10 06:43 data/prmsl.2006.nc\r\n", "-rw-r--r-- 1 jwhitaker staff 8974360 Jul 10 06:43 data/prmsl.2007.nc\r\n", "-rw-r--r-- 1 jwhitaker staff 8994260 Jul 10 06:43 data/prmsl.2008.nc\r\n", "-rw-r--r-- 1 jwhitaker staff 8974678 Jul 10 06:43 data/prmsl.2009.nc\r\n", "-rw-r--r-- 1 jwhitaker staff 8970732 Jul 10 06:43 data/prmsl.2010.nc\r\n", "-rw-r--r-- 1 jwhitaker staff 8976285 Jul 10 06:43 data/prmsl.2011.nc\r\n" ] } ], "source": [ "!ls -l data/prmsl*nc" ] }, { "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", "execution_count": 27, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 43, "slide_helper": "subslide_end" }, "slide_helper": "slide_end", "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "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" ] } ], "source": [ "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))" ] }, { "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", "execution_count": 28, "metadata": { "collapsed": false, "internals": { "frag_helper": "fragment_end", "frag_number": 45 }, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "f.close()\n", "gfs.close()" ] }, { "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": { "celltoolbar": "Raw Cell Format", "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.9" } }, "nbformat": 4, "nbformat_minor": 0 }