Skip to content

Instantly share code, notes, and snippets.

@jswhit

jswhit/reading_netCDF.ipynb

Created October 18, 2014 17:32
Show Gist options
  • Save jswhit/43dc116595ed66a92804 to your computer and use it in GitHub Desktop.
Save jswhit/43dc116595ed66a92804 to your computer and use it in GitHub Desktop.
Download ZIP
reading netcdf ipython notebook
Display the source blob
Display the rendered blob
Raw
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": {}
}
]
}
@majdo19
Copy link

majdo19 commented Mar 28, 2019

hi ,

i had question how you can plot all netcdf4 with the same variable prmsl i mean plot all years ? or choose what we want to plot !!?
thanks you

@iditbela2
Copy link

Thanks a lot for the code!
I think maybe in getclosest_ij you're missing a np.meshgrid first.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment