golive4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# golive learning notebook\n",
    "\n",
    "## Recommendation\n",
    "\n",
    "> ***locate*** and ***work through*** the \n",
    "[Python Data Science Handbook](https://jakevdp.github.io/PythonDataScienceHandbook/)\n",
    "by Jake VanderPlas, available free online. \n",
    "\n",
    "\n",
    "## golive notebooks 1, 2, 3\n",
    "\n",
    "\n",
    "The first three **golive** notebooks cover access to and analysis of a large geospatial dataset using the\n",
    "Python stack: Python 3, numpy, matplotlib, pandas, xarray, dask, and ultimately the pangeo \n",
    "framework; with digressions to widgets and ipyleaflet maps and so on. \n",
    "\n",
    "\n",
    "Data and a scientific problem:\n",
    "\n",
    "\n",
    "> * The *golive* dataset is land ice velocity fields from 2013 through present\n",
    ">   * For our purposes: mesoscale valley glacier systems\n",
    "> * What sort of speed variability patterns are discernible in the **golive** data?\n",
    "\n",
    "\n",
    "This notebook is focused on mechanics of Python tools with\n",
    "the *golive* data as the subject matter.\n",
    "\n",
    "Some prefatory notes: \n",
    "\n",
    "- `xarray` is abbreviated `xr`\n",
    "  - `xr.DataArray` and `xr.Dataset` are respectively the two principle data object constructs\n",
    "  - NetCDF data files translate into `xr.Dataset` objects by means of an `open()` call \n",
    "\n",
    "\n",
    "#### Contents\n",
    "\n",
    "flag re-do this"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## S0. utility code cell\n",
    "\n",
    "The *golive* notebooks make use of a common set of package imports and utility methods. \n",
    "The imports are copied by rote to each notebook and the methods are imported from \n",
    "`golive_utility.py`. The line-magic call `%matplotlib inline` and the two `pip install` calls \n",
    "are also rote. The `gdal` and `pygeotools` imports seem to be problematical (flag May 3 2019).\n",
    "\n",
    "In a perfect world this is cleaned up. In this world the corresponding utility code cell block \n",
    "is painstakingly re-copied from one notebook to the others in the `golive` collection to make \n",
    "sure these notebooks are all working from the same resource playbook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = 'all'  # default is ‘last_expr’\n",
    "\n",
    "# automatically reload external modules such as golive_library if they have changed\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: utm in /srv/conda/envs/notebook/lib/python3.7/site-packages (0.5.0)\n",
      "Requirement already satisfied: pygeotools in /srv/conda/envs/notebook/lib/python3.7/site-packages (0.6.0)\n",
      "Requirement already satisfied: gdal in /srv/conda/envs/notebook/lib/python3.7/site-packages (from pygeotools) (2.4.2)\n",
      "Requirement already satisfied: matplotlib in /srv/conda/envs/notebook/lib/python3.7/site-packages (from pygeotools) (3.1.1)\n",
      "Requirement already satisfied: numpy in /srv/conda/envs/notebook/lib/python3.7/site-packages (from pygeotools) (1.17.1)\n",
      "Requirement already satisfied: scipy in /srv/conda/envs/notebook/lib/python3.7/site-packages (from pygeotools) (1.3.1)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib->pygeotools) (1.1.0)\n",
      "Requirement already satisfied: cycler>=0.10 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib->pygeotools) (0.10.0)\n",
      "Requirement already satisfied: python-dateutil>=2.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib->pygeotools) (2.7.5)\n",
      "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib->pygeotools) (2.4.2)\n",
      "Requirement already satisfied: setuptools in /srv/conda/envs/notebook/lib/python3.7/site-packages (from kiwisolver>=1.0.1->matplotlib->pygeotools) (41.2.0)\n",
      "Requirement already satisfied: six in /srv/conda/envs/notebook/lib/python3.7/site-packages (from cycler>=0.10->matplotlib->pygeotools) (1.12.0)\n",
      "Collecting tqdm\n",
      "\u001b[?25l  Downloading https://files.pythonhosted.org/packages/dc/88/d3213e2f3492daf09d8b41631ad6899f56db17ce83ea9c8a579902bafe5e/tqdm-4.35.0-py2.py3-none-any.whl (50kB)\n",
      "\u001b[K     |████████████████████████████████| 51kB 427kB/s  eta 0:00:01\n",
      "\u001b[?25hInstalling collected packages: tqdm\n",
      "Successfully installed tqdm-4.35.0\n"
     ]
    }
   ],
   "source": [
    "!pip install utm\n",
    "!pip install pygeotools\n",
    "!pip install tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "# Python Standard Library built-in modules\n",
    "import sys\n",
    "import os\n",
    "import timeit\n",
    "from datetime import timedelta, datetime\n",
    "from random import randint as ri\n",
    "\n",
    "# Third-party modules\n",
    "from tqdm import tqdm\n",
    "from matplotlib import pyplot as plt\n",
    "import pandas as pd\n",
    "import netCDF4\n",
    "import xarray as xr\n",
    "import dask\n",
    "import numpy as np\n",
    "from numpy import datetime64 as dt64, timedelta64 as td64\n",
    "from scipy import stats\n",
    "from scipy.signal import medfilt as mf\n",
    "from ftplib import FTP\n",
    "from ipywidgets import *                # interactivity\n",
    "from traitlets import dlink             # interactivity\n",
    "import utm\n",
    "# from osgeo import gdal \n",
    "# from pygeotools.lib import iolib,timelib\n",
    "\n",
    "# Local application/library specific imports\n",
    "from golive_library import GoliveLibrary as g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## S1. Types and data manipulation packages \n",
    "\n",
    "\n",
    "As noted above the [Python Data Science Handbook](https://jakevdp.github.io/PythonDataScienceHandbook/)\n",
    "is an excellent preparatory resource. It is freely available online where chapters 1 and 2 would be (ideally)\n",
    "mastered first before proceeding here. \n",
    "\n",
    "\n",
    "Python has four native data types: float, int, bool and string. A string behaves something like a list but it is \n",
    "immutable whereas a list has the `.append()` method; it can be modified during execution. \n",
    "\n",
    "\n",
    "The *Numerical Python* library `numpy` (abbreviated `np`) introduces arrays. \n",
    "In fact these arrays can be multi-dimensional and of type `ndarray` meaning *n-dimensional array*. \n",
    "\n",
    "\n",
    "Then we need to touch on `pandas` before jumping to `xarray`. flag this section is incomplete May 3 2019.\n",
    "\n",
    "\n",
    "### Data acquisition and storage\n",
    "\n",
    "\n",
    "The `golive` data files are NetCDF files resident on an ftp server. Unfortunately (against principles) we need\n",
    "to copy them to this local Jupyter environment. For this we use the `wget` command. The files are place in:\n",
    "\n",
    "```\n",
    "<home>/data/golive\n",
    "<home>/data/golive_big\n",
    "``` \n",
    "\n",
    "One single file is typically 15MB. It covers a single time window spanned by two source LANDSAT-8 images.\n",
    "This time window is one of: 16, 32, 48, 64, 80, or 96 days. These files are generated automatically and \n",
    "feature good data and plenty of no-data values from obscurement (clouds). LANDSAT frames are designated \n",
    "by path and row coordinates and each covers approximately 80 x 80 kilometers. Path and row coordinates \n",
    "are intrinsic to satellite tracks, not earth coordinates. The good news is that the data are organized\n",
    "by UTM coordinates, i.e. earth coordinates. The slightly bad news is that the UTM zone is a bit \n",
    "tricky to establish (which one?) and without figuring this out we can get frames from different\n",
    "path-row pairs that appear to have the same coordinates (when in fact they are in different UTM zones).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## S2. IPython and Jupyter notebooks\n",
    "\n",
    "\n",
    "### IPython\n",
    "\n",
    "\n",
    "### Jupyter shortcuts\n",
    "\n",
    "\n",
    "[This blog post](https://www.dataquest.io/blog/jupyter-notebook-tips-tricks-shortcuts/) gives useful shortcuts\n",
    "for operating Jupyter notebooks including the idea of an ESC escape mode:\n",
    "\n",
    "- **a** adds a cell above\n",
    "- **b** adds a cell below\n",
    "- **m** changes a selected cell to markdown\n",
    "- **y** changes a cell to code\n",
    "- **d twice** deletes the current cell\n",
    "\n",
    "- **Ctrl Shift minus-sign** will split the current cell into two from where your cursor is\n",
    "- **f** Find and replace on your code but not the outputs\n",
    "- **o** Toggle cell output: Frequently used here!\n",
    "\n",
    "\n",
    "- **Shift + J** or **Shift + down** does multi-select downwards\n",
    "- **Shift + K** or **Shift + up** does upwards\n",
    "- **Shift + M** to merge\n",
    "\n",
    "- **Shift + Tab** will show the Docstring (documentation) for an object or method where your cursor is. If you have just typed the object into a cell you can keep pressing **Shift + Tab** to cycle through the modes of documentation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Here is a very useful tip for accessing object methods and sub-methods...\n",
      "\n",
      "- Type the name of the object with trailing period\n",
      "- place the cursor after this period\n",
      "- hit Shift + Tab multiple times\n",
      "\n",
      "this will display docstring variants\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print('''\n",
    "Here is a very useful tip for accessing object methods and sub-methods...\n",
    "\n",
    "- Type the name of the object with trailing period\n",
    "- place the cursor after this period\n",
    "- hit Shift + Tab multiple times\n",
    "\n",
    "this will display docstring variants\n",
    "\n",
    "''')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## S3. Exploratory session on golive DataSets\n",
    "\n",
    "\n",
    "This section 'deconstructs' a small set of NetCDF files first by pulling them into a list of `xarray Datasets`. \n",
    "This goes against principle! Better is to concatenate them into a *single* such `Dataset`.  This is \n",
    "done properly further below. Here we use list m\\[\\] as a list of six `xarray DataSets`. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "golive_d = g.set_data_directory('golive_big')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# flag this section needs rewrite to import data in a clean manner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.DataArray 'vx_masked' (y: 828, x: 819)>\n",
      "dask.array<shape=(828, 819), dtype=float32, chunksize=(828, 819)>\n",
      "Coordinates:\n",
      "  * x        (x) float64 3.331e+05 3.334e+05 3.337e+05 ... 5.782e+05 5.785e+05\n",
      "  * y        (y) float64 6.786e+06 6.785e+06 6.785e+06 ... 6.538e+06 6.537e+06\n",
      "Attributes:\n",
      "    grid_mapping:   transverse_mercator\n",
      "    standard_name:  x_velocity_masked\n",
      "    long_name:      x component of velocity (masked)\n",
      "    units:          m/d\n",
      "    masking_info:   masked_where(((del_corr_arr<0.100)&(corr_arr<1.000))|(cor...\n",
      "\n",
      "\n",
      "<xarray.DataArray 'x' (x: 819)>\n",
      "array([333142.5, 333442.5, 333742.5, ..., 577942.5, 578242.5, 578542.5])\n",
      "Coordinates:\n",
      "  * x        (x) float64 3.331e+05 3.334e+05 3.337e+05 ... 5.782e+05 5.785e+05\n",
      "Attributes:\n",
      "    standard_name:  projection_x_coordinate\n",
      "    long_name:      x coordinate of projection\n",
      "    units:          m\n",
      "\n",
      "\n",
      "2016-11-30T21:02:58.936365 : /projects/makl5454/bin/sc_pycorr_nc_oc_v5p8.py -img1dir /lustre/janus_scratch/makl5454/Landsat/ak_hp/p063_r018 -img2dir /lustre/janus_scratch/makl5454/Landsat/ak_hp/p063_r018 LC80630182014050LGN00_B8_hp.tif LC80630182014066LGN00_B8_hp.tif -inc 20 -half_source_chip 10 -trackvmax 15.000000 -plotvmax 15.0 -out_name_base ak_v02_S8_063_018_016_2014_050_2014_066 -log10 -dcam 0.1 -cam 1.0 -cam1 0.5 -output_dir /work/scambos_landsat/Landsat/GO_LIVE/ak_out/p063_r018 -VRT_dir /lustre/janus_scratch/makl5454/Landsat -no_speed_ref -offset_correction_lgo_mask -offset_correction_bilinear_fit -lgo_mask_file_dir /lustre/janus_scratch/makl5454/Landsat -lgo_mask_filename AK_ocean_land_glacier_210_mask_100m_Kienholz_v3.tif  ====> sc_pycorr run\n",
      "\n",
      "\n",
      "<xarray.DataArray 'processing_log' (chars: 4284)>\n",
      "dask.array<shape=(4284,), dtype=|S1, chunksize=(4284,)>\n",
      "Coordinates:\n",
      "  * chars    (chars) int32 0 1 2 3 4 5 6 ... 4277 4278 4279 4280 4281 4282 4283\n",
      "Attributes:\n",
      "    standard_name:  sc_pycorr_process_log\n",
      "    long_name:      log from feature tracking processing\n",
      "\n",
      "\n",
      "sc_pycorr_process_log\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(m[0]['vx_masked'])\n",
    "print('\\n')\n",
    "# Identical: print(m[0].vx_masked)\n",
    "\n",
    "print(m[0].vx_masked.x)\n",
    "print('\\n')\n",
    "print(m[0].history)\n",
    "print('\\n')\n",
    "print(m[0].processing_log)\n",
    "print('\\n')\n",
    "print(m[0].processing_log.attrs['standard_name'])\n",
    "print('\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What *is* the ```processing_log```? By ```type``` it is an ```xarray DataArray```. Taking ```.values``` \n",
    "gives a ```numpy ndarray```. \n",
    "While the following is not optimal code it *does* serve to print the content of ```processing_log```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"# log from /projects/makl5454/bin/sc_pycorr_nc_oc_v5p8.py\\n/projects/makl5454/bin/sc_pycorr_nc_oc_v5p8.py -img1dir /lustre/janus_scratch/makl5454/Landsat/ak_hp/p063_r018 -img2dir /lustre/janus_scratch/makl5454/Landsat/ak_hp/p063_r018 LC80630182014050LGN00_B8_hp.tif LC80630182014066LGN00_B8_hp.tif -inc 20 -half_source_chip 10 -trackvmax 15.000000 -plotvmax 15.0 -out_name_base ak_v02_S8_063_018_016_2014_050_2014_066 -log10 -dcam 0.1 -cam 1.0 -cam1 0.5 -output_dir /work/scambos_landsat/Landsat/GO_LIVE/ak_out/p063_r018 -VRT_dir /lustre/janus_scratch/makl5454/Landsat -no_speed_ref -offset_correction_lgo_mask -offset_correction_bilinear_fit -lgo_mask_file_dir /lustre/janus_scratch/makl5454/Landsat -lgo_mask_filename AK_ocean_land_glacier_210_mask_100m_Kienholz_v3.tif\\nNamespace(Greenland=False, VRT_dir='/lustre/janus_scratch/makl5454/Landsat', bbox=None, cam=1.0, cam1=0.5, datestrfmt='%m/%d/%y', dcam=0.1, debug_tifs=False, half_source_chip=10, half_target_chip=40, img1_name='LC80630182014050LGN00_B8_hp.tif', img1datestr=None, img1dir='/lustre/janus_scratch/makl5454/Landsat/ak_hp/p063_r018', img2_name='LC80630182014066LGN00_B8_hp.tif', img2datestr=None, img2dir='/lustre/janus_scratch/makl5454/Landsat/ak_hp/p063_r018', imgdir='.', inc=20, lgo_mask_file_dir='/lustre/janus_scratch/makl5454/Landsat', lgo_mask_filename='AK_ocean_land_glacier_210_mask_100m_Kienholz_v3.tif', log10=True, max_allowable_pixel_offset_correction=3.0, mpy=False, nlf=False, no_speed_ref=True, nodatavmax=0.333, offset_correction_bilinear_fit=True, offset_correction_bilinear_fit_min_pixels=2000, offset_correction_lgo_mask=True, offset_correction_speedref=False, out_name_base='ak_v02_S8_063_018_016_2014_050_2014_066', output_dir='/work/scambos_landsat/Landsat/GO_LIVE/ak_out/p063_r018', plotvmax=15.0, progupdates=False, speed_ref_dir='/projects/makl5454', speed_ref_filename='LISA_750m_vv.tif', trackvmax=15.0, v=False)\\nUTC: 12/01/16_03:58:34  Wall:   0.004  Proc:   0.010 - Start Program\\n>>>> half_target_chip 40 pixels - maximum trackable velocity 28.125000 m/d (unless speed_ref at that location is lower, then it is 1.5 * (speed_ref_vel/pix_size_x_m) * del_t_days + 5 pixels) <<<<time separation in days: 16.000000 years: 0.043806 - using m/d\\n# got lgo mask data and reprojected to original L8 image local utm# lgo mask reference (lgo mask file: /lustre/janus_scratch/makl5454/Landsat/AK_ocean_land_glacier_210_mask_100m_Kienholz_v3.tif) from lgo mask vrt file: tmp3CR5G0lgo.vrt\\nUTC: 12/01/16_04:02:58  Wall: 263.973  Proc: 263.160 - Done with image to image correlation: 416178 correlations of 678132 possible\\nUTC: 12/01/16_04:02:58  Wall: 263.991  Proc: 263.170 - found lgo_mask (land pixel) offset correction (del_i: -0.000000 del_j: -0.000000 std_i 0.436644 std_j 0.553456) using 64801 pixels out of 88626 possible\\nUTC: 12/01/16_04:02:58  Wall: 264.155  Proc: 263.270 - using lgo_bilinear_masked_correction bilinear offset correction - constant offsets would have been -0.000000 pixels in i, -0.000000 pixels in j lgo_bilinear_masked_correction fields for i and j for land pixels, - using bilinear coefs_i [ -2.43355342e-01   7.92498207e-04   1.44520584e-04  -5.32063129e-07] coefs_j [  3.71111151e-01  -6.65638161e-04   5.16895948e-05  -1.82916670e-06] 64801 valid pixels out of 88626 possible for scene (73.117370 %)\\nout image /work/scambos_landsat/Landsat/GO_LIVE/ak_out/p063_r018/ak_v02_S8_063_018_016_2014_050_2014_066_log10.tif GTiff 819, 828, 4\\nadded  applied_bilinear_x_offset_correction_in_pixels to netCDF file 819, 828 added  applied_bilinear_y_offset_correction_in_pixels to netCDF file 819, 828 added  vx to netCDF file 819, 828 added  vy to netCDF file 819, 828 added  vv to netCDF file 819, 828 added  vx_masked to netCDF file 819, 828 added  vy_masked to netCDF file 819, 828 added  vv_masked to netCDF file 819, 828 added  corr to netCDF file 819, 828 added  del_corr to netCDF file 819, 828 added  lgo_mask to netCDF file 819, 828 added  d2idx2 to netCDF file 819, 828 added  d2jdx2 to netCDF file 819, 828 added  del_i to netCDF file 819, 828 added  del_j to netCDF file 819, 828 UTC: 12/01/16_04:02:59  Wall: 265.390  Proc: 264.360 - At end - resource module reports process ak_v02_S8_063_018_016_2014_050_2014_066 using   0.000 local   0.000 stack   1.632 max memory\\n\""
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = m[0].processing_log.values\n",
    "# print(type(b)) gives numpy.ndarray\n",
    "# print(chr(b[0][0])) gives '35'\n",
    "# print(len(b)) gives 4284 characters\n",
    "# s = ''\n",
    "# for i in range(len(b)): s = s + chr(b[i][0])\n",
    "s = ''.join(chr(b[i][0]) for i in range(len(b)))\n",
    "s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.DataArray 'vv_masked' (y: 828, x: 819)>\n",
      "dask.array<shape=(828, 819), dtype=float32, chunksize=(828, 819)>\n",
      "Coordinates:\n",
      "  * x        (x) float64 3.331e+05 3.334e+05 3.337e+05 3.34e+05 3.343e+05 ...\n",
      "  * y        (y) float64 6.786e+06 6.785e+06 6.785e+06 6.785e+06 6.784e+06 ...\n",
      "Attributes:\n",
      "    grid_mapping:   transverse_mercator\n",
      "    standard_name:  speed_masked\n",
      "    long_name:      magnitude of velocity (masked)\n",
      "    units:          m/d\n",
      "    masking_info:   masked_where(((del_corr_arr<0.100)&(corr_arr<1.000))|(cor...\n",
      "\n",
      "\n",
      "OrderedDict([('GDAL_AREA_OR_POINT', 'Area'), ('Conventions', 'CF-1.6'), ('history', '2016-11-30T21:02:58.936365 : /projects/makl5454/bin/sc_pycorr_nc_oc_v5p8.py -img1dir /lustre/janus_scratch/makl5454/Landsat/ak_hp/p063_r018 -img2dir /lustre/janus_scratch/makl5454/Landsat/ak_hp/p063_r018 LC80630182014050LGN00_B8_hp.tif LC80630182014066LGN00_B8_hp.tif -inc 20 -half_source_chip 10 -trackvmax 15.000000 -plotvmax 15.0 -out_name_base ak_v02_S8_063_018_016_2014_050_2014_066 -log10 -dcam 0.1 -cam 1.0 -cam1 0.5 -output_dir /work/scambos_landsat/Landsat/GO_LIVE/ak_out/p063_r018 -VRT_dir /lustre/janus_scratch/makl5454/Landsat -no_speed_ref -offset_correction_lgo_mask -offset_correction_bilinear_fit -lgo_mask_file_dir /lustre/janus_scratch/makl5454/Landsat -lgo_mask_filename AK_ocean_land_glacier_210_mask_100m_Kienholz_v3.tif  ====> sc_pycorr run'), ('GDAL', 'GDAL 2.0.2, released 2016/01/26')])\n",
      "\n",
      "\n",
      "Coordinates:\n",
      "  * x        (x) float64 3.331e+05 3.334e+05 3.337e+05 3.34e+05 3.343e+05 ...\n",
      "  * y        (y) float64 6.786e+06 6.785e+06 6.785e+06 6.785e+06 6.784e+06 ...\n",
      "  * chars    (chars) int32 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ...\n",
      "\n",
      "\n",
      "<xarray.DataArray 'vx_masked' (y: 828, x: 819)>\n",
      "dask.array<shape=(828, 819), dtype=float32, chunksize=(828, 819)>\n",
      "Coordinates:\n",
      "  * x        (x) float64 3.331e+05 3.334e+05 3.337e+05 3.34e+05 3.343e+05 ...\n",
      "  * y        (y) float64 6.786e+06 6.785e+06 6.785e+06 6.785e+06 6.784e+06 ...\n",
      "Attributes:\n",
      "    grid_mapping:   transverse_mercator\n",
      "    standard_name:  x_velocity_masked\n",
      "    long_name:      x component of velocity (masked)\n",
      "    units:          m/d\n",
      "    masking_info:   masked_where(((del_corr_arr<0.100)&(corr_arr<1.000))|(cor...\n",
      "\n",
      "\n",
      "[[nan nan nan ... nan nan nan]\n",
      " [nan nan nan ... nan nan nan]\n",
      " [nan nan nan ... nan nan nan]\n",
      " ...\n",
      " [nan nan nan ... nan nan nan]\n",
      " [nan nan nan ... nan nan nan]\n",
      " [nan nan nan ... nan nan nan]]\n"
     ]
    }
   ],
   "source": [
    "print(m[0].vv_masked)\n",
    "print('\\n')\n",
    "print(m[0].attrs)\n",
    "print('\\n')\n",
    "print(m[0].coords)\n",
    "print('\\n')\n",
    "print(m[0]['vx_masked'])\n",
    "print('\\n')\n",
    "print(m[0]['vv_masked'].values)\n",
    "\n",
    "# - The array is capable of truncating itself to give us a sense of what is there without printing *everything*\n",
    "# - There are a lot of not-a-number values in the speed field (notice the focus is on 'vv' which is speed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('y', 'x')\n",
      "\n",
      " Coordinates:\n",
      "  * x        (x) float64 3.331e+05 3.334e+05 3.337e+05 3.34e+05 3.343e+05 ...\n",
      "  * y        (y) float64 6.786e+06 6.785e+06 6.785e+06 6.785e+06 6.784e+06 ...\n",
      "\n",
      " vv_masked\n",
      "\n",
      " OrderedDict([('grid_mapping', 'transverse_mercator'), ('standard_name', 'speed_masked'), ('long_name', 'magnitude of velocity (masked)'), ('units', 'm/d'), ('masking_info', 'masked_where(((del_corr_arr<0.100)&(corr_arr<1.000))|(corr_arr<0.500))')]) \n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[0.21738425, 0.3744488 ],\n",
       "       [0.2512998 , 0.27024233]], dtype=float32)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Further proof that y comes first...\n",
    "print(m[0]['vv_masked'].dims)\n",
    "\n",
    "# x is listed first by y comes first in referencing locations; and notice x increases west to east\n",
    "#   whereas y decreases north to south...\n",
    "print('\\n', m[0]['vv_masked'].coords)\n",
    "\n",
    "print('\\n', m[0]['vv_masked'].name)\n",
    "\n",
    "print('\\n', m[0]['vv_masked'].attrs, '\\n')\n",
    "\n",
    "# integer indices...\n",
    "#   m[0]['vv_masked'][a:b,c:d].values uses integers a, b, c, d to print a matrix of speeds, as in...\n",
    "a = 249        # the y values, indexing north to south decreasing numerical (see image above)\n",
    "b = 251        #   and this corresponds to increasing index\n",
    "c = 749        # the x values, indexing west to east increasing numerical\n",
    "d = 751        #   and this corresponds to increasing index\n",
    "m[0]['vv_masked'][a:b,c:d].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coordinates:\n",
      "    x        float64 5.281e+05\n",
      "    y        float64 6.638e+06 \n",
      "\n",
      "[[      nan       nan       nan       nan       nan]\n",
      " [      nan       nan       nan       nan       nan]\n",
      " [      nan       nan       nan       nan       nan]\n",
      " [      nan       nan       nan 18.035414       nan]\n",
      " [      nan       nan       nan       nan       nan]\n",
      " [      nan       nan       nan       nan       nan]] \n",
      "\n",
      "[[0.41044605 0.4309364  0.31737426 0.40143636 0.35311106]\n",
      " [0.28281084 0.27785343 0.38837793 0.32612908 0.3373442 ]\n",
      " [0.3475631  0.30343953 0.29954764 0.29412556 0.29397836]\n",
      " [0.47786003 0.41999614 0.2908547  0.25681612 0.38776925]\n",
      " [0.30885625 0.2575876  0.24626212 0.34321824 0.31327507]\n",
      " [0.37430257 0.35534668 0.40519747 0.2731533  0.44645512]] \n",
      "\n",
      "[[0.35512862 0.35785887 0.417958   0.3520321  0.38790464]\n",
      " [0.43297464 0.46064052 0.57191545 0.5078081  0.39370164]\n",
      " [0.7534419         nan 0.3011668  0.38859782 0.5216916 ]\n",
      " [0.4815417  0.37989134 0.39068517 0.43206763        nan]\n",
      " [0.45732138 0.49818754        nan        nan        nan]\n",
      " [       nan        nan        nan 0.2728113  0.27249703]] \n",
      "\n",
      "[[0.09038115 0.123815   0.09873071 0.19618239 0.10258889]\n",
      " [0.14093068 0.21738158 0.19509925 0.2555801  0.2614072 ]\n",
      " [0.12077411 0.18787687 0.21427721 0.15802394 0.11620727]\n",
      " [0.10328206 0.16715579 0.1077848  0.11574709 0.13659836]\n",
      " [0.1872416  0.17457785 0.11892354 0.13936782 0.09043686]\n",
      " [0.11191403 0.17082685 0.20444609        nan 0.09255682]] \n",
      "\n",
      "[[nan nan nan nan nan]\n",
      " [nan nan nan nan nan]\n",
      " [nan nan nan nan nan]\n",
      " [nan nan nan nan nan]\n",
      " [nan nan nan nan nan]\n",
      " [nan nan nan nan nan]] \n",
      "\n",
      "[[0.14655106        nan        nan        nan 0.10225438]\n",
      " [       nan        nan        nan        nan 0.13121285]\n",
      " [       nan 0.13210182 0.10893457 0.1193234  0.1095689 ]\n",
      " [0.15275988 0.14622024 0.11876038 0.12644033 0.12120885]\n",
      " [0.13783072 0.12154137 0.10713261 0.12358696 0.11866077]\n",
      " [0.14726122 0.13806409 0.12625797 0.12898837 0.10601326]] \n",
      "\n",
      "[[0.28845277 0.23994552 0.28734905 0.2431491  0.19211332]\n",
      " [0.1910194  0.16813856 0.23149344 0.23171276 0.20612305]\n",
      " [0.13903503 0.15781759 0.19114526 0.18024541 0.18726926]\n",
      " [0.20102477 0.2098     0.19125398 0.15796323 0.274036  ]\n",
      " [0.16950211 0.16453557 0.20032121 0.18054876 0.2106169 ]\n",
      " [0.18694265 0.2153948  0.1794939  0.06814034 0.1921458 ]] \n",
      "\n",
      "[[0.08426813 0.11589131 0.10924371 0.12409929 0.1140882 ]\n",
      " [0.09596647 0.1092534  0.1066667  0.117033   0.4785427 ]\n",
      " [0.17040925 0.12402933 0.1261201  0.09994841        nan]\n",
      " [0.05855089 0.12894534 0.11201116 0.11156595 0.12929066]\n",
      " [0.08244964 0.09650747 0.14457332 0.17609367 0.19434917]\n",
      " [       nan        nan 0.23721068 0.26291707 0.25492975]] \n",
      "\n",
      "[[0.01771249 0.00893011 0.01787673 0.01809397 0.02767793]\n",
      " [0.00830157 0.00861039 0.00819861 0.00824569 0.02333687]\n",
      " [0.01032019 0.00808949 0.01716333 0.0220554  0.01759464]\n",
      " [0.0175069  0.00378544 0.01245362 0.0310574  0.03644016]\n",
      " [0.02767775 0.02632458 0.03681707 0.016669   0.02691138]\n",
      " [0.04968441 0.03055906 0.02554211 0.01319485 0.01207088]] \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# integer indexing is possible: y is 493, x is 650 in this case\n",
    "a=m[0]['vv_masked'][493,650]\n",
    "print(a.coords, '\\n')\n",
    "\n",
    "# Demonstrating an anomalous (18 meter per day) speed in the midst of some nans in m[0]; and then\n",
    "#   looking into the subsequent 8 DataSets at the same integer indices. \n",
    "#   WARNING! These are not physically the same locations. Each m[] DataSet has its own spatial mapping\n",
    "for i in range(9):\n",
    "    a=m[i]['vv_masked'].isel(y=slice(490,496),x=slice(647,652))\n",
    "    print(a.values, '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "meters per day\n",
      "meters per day\n",
      "m/d\n"
     ]
    }
   ],
   "source": [
    "# An example of overriding existing content using attributes and =\n",
    "# print(type(m[0].vv_masked)) gives same\n",
    "speed = m[2].vv_masked\n",
    "# print(type(speed)) gives DataArray\n",
    "# print(speed.attrs)    # maps units to 'm/d'\n",
    "speed.attrs['units']='meters per day'\n",
    "print(speed.attrs['units'])             # units is now 'meters per day'\n",
    "print(m[2].vv_masked.attrs['units'])    #   notice 'speed' is not a separate version of m[2].vv_masked\n",
    "print(m[3].vv_masked.attrs['units'])    # units unchanged\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## S3. directory() operating on xarray DataArray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'xarray.core.dataarray.DataArray'>\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['T',\n",
       " 'all',\n",
       " 'any',\n",
       " 'argmax',\n",
       " 'argmin',\n",
       " 'argsort',\n",
       " 'assign_attrs',\n",
       " 'assign_coords',\n",
       " 'astype',\n",
       " 'attrs',\n",
       " 'bfill',\n",
       " 'broadcast_equals',\n",
       " 'chunk',\n",
       " 'chunks',\n",
       " 'clip',\n",
       " 'close',\n",
       " 'combine_first',\n",
       " 'compute',\n",
       " 'conj',\n",
       " 'conjugate',\n",
       " 'coords',\n",
       " 'copy',\n",
       " 'count',\n",
       " 'cumprod',\n",
       " 'cumsum',\n",
       " 'data',\n",
       " 'diff',\n",
       " 'dims',\n",
       " 'dot',\n",
       " 'drop',\n",
       " 'dropna',\n",
       " 'dt',\n",
       " 'dtype',\n",
       " 'encoding',\n",
       " 'equals',\n",
       " 'expand_dims',\n",
       " 'ffill',\n",
       " 'fillna',\n",
       " 'from_cdms2',\n",
       " 'from_dict',\n",
       " 'from_iris',\n",
       " 'from_series',\n",
       " 'get_axis_num',\n",
       " 'get_index',\n",
       " 'grid_mapping',\n",
       " 'groupby',\n",
       " 'groupby_bins',\n",
       " 'identical',\n",
       " 'imag',\n",
       " 'indexes',\n",
       " 'interp',\n",
       " 'interp_like',\n",
       " 'interpolate_na',\n",
       " 'isel',\n",
       " 'isel_points',\n",
       " 'isin',\n",
       " 'isnull',\n",
       " 'item',\n",
       " 'load',\n",
       " 'loc',\n",
       " 'long_name',\n",
       " 'masking_info',\n",
       " 'max',\n",
       " 'mean',\n",
       " 'median',\n",
       " 'min',\n",
       " 'name',\n",
       " 'nbytes',\n",
       " 'ndim',\n",
       " 'notnull',\n",
       " 'persist',\n",
       " 'pipe',\n",
       " 'plot',\n",
       " 'prod',\n",
       " 'quantile',\n",
       " 'rank',\n",
       " 'real',\n",
       " 'reduce',\n",
       " 'reindex',\n",
       " 'reindex_like',\n",
       " 'rename',\n",
       " 'reorder_levels',\n",
       " 'resample',\n",
       " 'reset_coords',\n",
       " 'reset_index',\n",
       " 'roll',\n",
       " 'rolling',\n",
       " 'round',\n",
       " 'searchsorted',\n",
       " 'sel',\n",
       " 'sel_points',\n",
       " 'set_index',\n",
       " 'shape',\n",
       " 'shift',\n",
       " 'size',\n",
       " 'sizes',\n",
       " 'sortby',\n",
       " 'squeeze',\n",
       " 'stack',\n",
       " 'standard_name',\n",
       " 'std',\n",
       " 'sum',\n",
       " 'swap_dims',\n",
       " 'to_cdms2',\n",
       " 'to_dataframe',\n",
       " 'to_dataset',\n",
       " 'to_dict',\n",
       " 'to_index',\n",
       " 'to_iris',\n",
       " 'to_masked_array',\n",
       " 'to_netcdf',\n",
       " 'to_pandas',\n",
       " 'to_series',\n",
       " 'transpose',\n",
       " 'units',\n",
       " 'unstack',\n",
       " 'values',\n",
       " 'var',\n",
       " 'variable',\n",
       " 'where',\n",
       " 'x',\n",
       " 'y']"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For a sense of how much functionality is built into a DataArray\n",
    "print(type(m[0].vv_masked))\n",
    "g.dirnou(m[0]['vv_masked'])                 # dirnou() filters out underscore methods ('no-underscore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## S4. ipyleaflet maps\n",
    "\n",
    "### using map on_interaction() with ipywidget 'label'\n",
    "\n",
    "The following code won't do anything of interest; it needs to follow a cell in which an ipyleaflet\n",
    "*map* is instantiated (called 'm'). This code registers a function called handle_interaction() with\n",
    "the interaction response mechanism built into the map, a method called ```.on_interaction(function)```.\n",
    "\n",
    "\n",
    "The code here is all commented out so it can be blindly run out of context without generating error ink."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Experimental: ipywidget 'label' is a live-readout coupled to 'interaction' events on the map\n",
    "#   This works to generate coordinates from left/right mouse clicks lists but deprecated for ipyleaflet map primitives.\n",
    "# label = Label('label info to be viewed here')\n",
    "# display(label)\n",
    "# a, b = [], []\n",
    "\n",
    "# def handle_interaction(**kwargs):\n",
    "#     t = kwargs.get('type')\n",
    "#     if t == 'mousemove': label.value = str(kwargs.get('coordinates'))\n",
    "#     else: \n",
    "#         label.value = t    # down-side: mouse clicks tend to produce multiple copies of coordinates > a[]/b[]\n",
    "#         if t == 'click': a.append(str(kwargs.get('coordinates')))\n",
    "#         elif t == 'contextmenu': b.append(str(kwargs.get('coordinates')))    # rt-click release\n",
    "            \n",
    "# Alternative handle_interaction() to look under the hood at **kwargs\n",
    "# kwargs keys consist of just 'event', 'type', and 'coordinates'\n",
    "# def handle_interaction(**kwargs):\n",
    "#     out = ''\n",
    "#     for key, value in kwargs.items(): out = out + str(key) + ' : ' + str(value) + ';    '\n",
    "#     label.value = out\n",
    "    \n",
    "# gmap.on_interaction(handle_interaction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell shows the results of some mouse-clicks when the prior cell is enabled and run\n",
    "#   Both list a[] and list b[] will often contain multiple copies of a coordinate, presumably\n",
    "#   from down-press and release of the mouse buttons.\n",
    "# print(a)\n",
    "# print('\\n')\n",
    "# print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This breaks the m[2] tif into 256 x 256 tiles in a sub-folder of golive_d called 'tiles'\n",
    "#   This should be demonstrated as a tile layer\n",
    "\n",
    "if False:\n",
    "    import os, gdal\n",
    "    out_path = golive_d + 'tiles/'\n",
    "    output_filename = 'tile_'\n",
    "\n",
    "    tile_size_x = 256  \n",
    "    tile_size_y = 256\n",
    "\n",
    "    ds = gdal.Open(ftif[2])\n",
    "    band = ds.GetRasterBand(1)\n",
    "    xsize = band.XSize\n",
    "    ysize = band.YSize\n",
    "\n",
    "    for i in range(0, xsize, tile_size_x):\n",
    "        for j in range(0, ysize, tile_size_y):\n",
    "            com_string = \"gdal_translate -of GTIFF -srcwin \" + str(i)+ \", \" + str(j) + \\\n",
    "                \", \" + str(tile_size_x) + \", \" + str(tile_size_y) + \" \" + ftif[2] + \" \" + \\\n",
    "                str(out_path) + str(output_filename) + str(i) + \"_\" + str(j) + \".tif\"\n",
    "            os.system(com_string)\n",
    "            \n",
    "            \n",
    "    g.show_local_image('', '/home/jovyan/data/golive/tiles/tile_512_256.tif', 256, 256)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### On ipyleaflet documentation\n",
    "\n",
    "ipyleaflet maps can be infuriating by doing nothing and reporting nothing. Two ideas to bear in mind...\n",
    "\n",
    "1. Search engines often produce good solution links from simple technical questions... but not always\n",
    "2. The [ipyleaflet documentation](https://ipyleaflet.readthedocs.io) was written for some reason in a very sparse manner\n",
    " \n",
    "It could be uncharitably described as *Pinteresque*, meaning\n",
    "(OED) \"typically characterized by implications of threat and strong feeling produced through colloquial \n",
    "language, apparent triviality, and long pauses.\"\n",
    "\n",
    "\n",
    "Here is an example of code that fails to overlay a single image on a map in Jupyter notebooks:\n",
    "\n",
    "\n",
    "```\n",
    "imageoverlay = ImageOverlay(url='file:/'+fpng[2], bounds=gmap.bounds)\n",
    "gmap.add_layer(imageoverlay)\n",
    "```\n",
    "\n",
    "The problem is the path: It must be relative to the current directory. \n",
    "\n",
    "\n",
    "### On Marker and CircleMarker features\n",
    "\n",
    "\n",
    "A **Marker** is a fixed-color pushpin that can be dragged about.\n",
    "\n",
    "\n",
    "```\n",
    "m1 = Marker(location=m.center, color='red')\n",
    "gmap += m1\n",
    "```\n",
    "\n",
    "\n",
    "A **CircleMarker** is more configurable than a Marker but is not draggable, alas.\n",
    "The color scheme shown here uses RGB in 12 bits via three hex characters: #00F will be blue. \n",
    "\n",
    "\n",
    "```\n",
    "loc1 = [m.center[0], m.center[1] + 0.06]\n",
    "cm0 = CircleMarker(radius=5, weight=1, color='#000', opacity=1.0, fill_opacity=1.0, fill_color='#F00')\n",
    "cm1 = CircleMarker(radius=5, weight=1, color='#000', opacity=1.0, fill_opacity=1.0, fill_color='#0F0')\n",
    "cm0.location = m.center          # to show separat modification of attributes\n",
    "cm1.location = loc1\n",
    "gmap.add_layer(cm0)\n",
    "gmap.add_layer(cm1)\n",
    "# get a CirlceMarker identifier: use \n",
    "cm0.model_id\n",
    "# delete marker from map: use \n",
    "gmap.remove_layer(cm0)\n",
    "# delete marker: use \n",
    "cm0.close()\n",
    "```\n",
    "\n",
    "## S5. Building custom colorbars \n",
    "\n",
    "Use as a model for building custom color maps. I added comments ad nauseum to the original. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAJpCAYAAAAXLM1yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X/UHWV97/33hyRoCWrEBOSXYivHqj2CNYsDi55zFNCDqUc8fbSFtjYqbWxXbfWpriPKOrT1tF3a1l/n0YrxR02f+rMoJY+lKlIouhYgAUH5ZYkslJiUECVAggp38n3+2HPrZjL3npm9Z2Zfs/fntda97j17Zu/53jt8ufZ1fee6RhGBmZlZig6adgBmZmZLcSNlZmbJciNlZmbJciNlZmbJciNlZmbJciNlZmbJciM1AUnHSQpJy5fY/1ZJH+46LrOUOW+sDjdSQyTdJekeSSuHnvttSVeO834R8RcR8duNBVhA0rMkbZF0X/bzZUnPavOcZsP6mDcAkg6R9DeSdkm6X9JVbZ/T6nMjdaDlwOunHUQN24GXA4cBq4HNwKemGpHNo77lDcBGBnnzzOz3/z3dcKyIG6kD/RXwJkmrarzmNZK2S9oh6Y2LT0r6E0l/nz1eHOJYL+m72be384eOPSnrET2QfSt9V5UTR8TuiLgrBkuHCNgHPL1G7GZN6FXeSHoG8FJgQ0TcGxH7IuL6GrFbR9xIHWgLcCXwphqveQFwPPAi4DxJZ4w49peAZwCnAxdIemb2/HuB90bE44GfAz6z+AJJ35D066MCkLQb+BHw/wB/USN2syb0LW/+E/Ad4E+zhu+bkv6vGrFbR9xIFbsA+ANJayoe/6cRsTcivgn8LXBOybE/jIibgJuAE7LnHwGeLml1ROyJiGsWXxARz4mIT4wKICJWAU8AXgd8vWLcZk3qU94cA/wCcD9wFIO82TTU+Fki3EgViIibgc8D51V8yd1Dj7/D4D/6pfz70OOHgEOzx+cC/wG4XdJ1kl5S8dw/ERF7gQuBv5N0eN3Xm02iZ3nzQwYN3J9FxMMR8a/AFQx6dZYQN1JL+2Pgd4CjKxx77NDjpzC4mKGWiLgjIs4BDgfeAVw0fLVUDQcBh1AtbrOm9SVvvlH3XDYdbqSWEBFbgU8Df1jh8P+VXc76bODV2etqkfSbktZExH5gd/b0vgqve6Gk50paJunxwLuA+4Db6sZgNqm+5A1wFfBd4C2Slks6FXg+8MW6MVi73EiN9jagyreyfwW2ApcDfx0RXxrjXGcCt0jaw6AYfHZE/AhA0i2SfmOJ160CPslgbP3bDK7sO3PxtWZTkHzeRMQjwFnAOga58yHgtyLi9jFisBbJNz00M7NUuSdlZmbJciNlZmbJciNlZmbJciNlZmbJciNlZmbJKryfS1u0WsFxXZ7RknUXxK5Q2WHSmQG7ar759V+MiDPHCyw9y1cpHjNqLYY58RCHTDuE6dv+Y2L3wlzlTaeNFMcxWIbSbG3VA3dR/z8ara75gqQ95ij4+f932lFM3w3hZfX4rapz9GcnbzzcZ2ZmyXIjZWZmyXIjZWZmyXIjZWZmyarUSElaJekiSbdLuk3SKZIOk3SZpDuy309sO1izPnHemE2uak/qvcAXIuLnGdwR8zYGNza7PCKOZ7CKcdUbnZnNC+eN2YRKG6nsHkX/BfgIQHYXy90MlrnflB22CXhZW0Ga9Y3zxqwZVXpSPwvcC/ytpK9L+nB258sjImIHQPbbtys3+ynnjVkDqjRSy4FfBD4QEc8F9lJjiELSBklbJG3h3jGjNOufxvJm4b62QjRLX5UVJ7YB2yLi2mz7IgbJdo+kIyNih6QjgZ1FL46IjcBGgLVSbCld0KMe+Z6N3Wj8g6685ERfNZY3T12jeO3GZoP74IZm38+W8KGGP+h7/7zZ9+uB0p5URPw7cLekZ2RPnQ7cCmwG1mfPrQcuaSVCsx5y3pg1o+rafX8AfFzSwcCdwKsZNHCfkXQu8F3gFe2EaNZbzhuzCVVqpCLiRorHZ05vNhyz2eG8MZucV5wwM7NkdXurjhaU35Ho0XyhRcYfxFyreyGGL7TINH0hhJVyT8rMzJLlRsrMzJLlRsrMzJLV+5pUXWU1rJkt1czsH2ZdKKthzWzNyjWoqXNPyszMkuVGyszMkuVGyszMkjV3Naky+ZpVb0s5vQ3c+ihfs+ptjco1qOS4J2VmZslyI2VmZslyI2VmZslyTWpWuAY1s9bsgg0T3PRwo8ssS0uhBlXrH6jhu1/2gHtSNnckPUPSjUM/D0h6Q+6Y50u6f+iYC6YVr1kKppU37knZ3ImIbwEnAkhaBnwPuLjg0K9ExEu6jM0sVdPKG/ekbN6dDnw7Ir4z7UDMeqSzvHFPqsTMzJuaL6slbRna3hgRSw3mnw18col9p0i6CdgOvCkibmkyyK5UqWc1XbeamXlTk+pXQTDJvHEjZbNoV0QU3bb9USQdDLwUeEvB7huAp0bEHknrgH8Ejm82TLOkJJk3Hu6zefZi4IaIuCe/IyIeiIg92eNLgRWSVncdoFmCOs0bN1I2z85hiSELSU+WpOzxSQxy5fsdxmaWqk7zxsN9feTC2MQkHQK8EHjt0HO/CxARFwIvB35P0gLwQ+DsiJjZDz5ft+pXKaWiNuZEzeQHtbRp5I0bKZtLEfEQ8KTccxcOPX4f8L6u4zJL2TTyxsN9ZmaWLDdSZmaWLA/3mdkB5qJGZb3gnpSZmSXLjZSZmSXLw31mZrNm9b3wKzVv65HoXUDcSNWUX8sPPG3JrEx+LT/o6Xp+Ls51zsN9ZmaWLDdSZmaWLDdSZmaWrEo1KUl3AQ8C+4CFiFgr6TDg08BxwF3Ar0bEfe2EadY/zhuzydW5cOIFEbFraPs84PKIeLuk87LtNzcanVn/zUTeeHKvTcskw31nAZuyx5uAl00ejtnMc96Y1VC1kQrgS5Kul7T4HeqIiNgBkP0+vOiFkjZI2iJpy72Tx2vWJ84bswlVHe47NSK2SzocuEzS7VVPEBEbyaaJrZVnFNlccd6YTahSTyoitme/dwIXAycB90g6EiD7vbOtIM36yHljNrnSRkrSSkmPW3wMvAi4GdgMrM8OWw9c0laQZn3jvDFrRpXhviOAi7Pb1i8HPhERX5B0HfAZSecC3wVe0V6YZr3jvDFrQGkjFRF3AicUPP994PQ2gjLrO+eNWTO84oSZmSXLjZSZmSXLjZSZmSXLjZSZmSXLjZTNJUl3SfqmpBslbSnYL0n/R9JWSd+Q9IvTiNMsJdPIG9+Z1+ZZfvHXYS8Gjs9+/hPwgey32bzrNG/ckzIrdhbwdzFwDbBqcaUIM1tS43njRspm0erFxVmzn6IbSxQt/jrsaODuoe1t2XNmsyrJvPFwn82iXRGxtuSYAxZ/jYirhvar4DVe6NVmWZJ5456UzaUlFn8dtg04dmj7GGB7N9GZpWkaeeNGyubOiMVfh20Gfiu7Wulk4P7F+0CZzaNp5Y2H+2weLbX46+8CRMSFwKXAOmAr8BDw6inFapaKqeSNGymbOyMWf71w6HEAv99lXGYpm1beeLjPzMyS5UbKzMyS5UbKzMyS5ZqUmZXaWDRt06wD7kmZmVmy3EiZmVmy3EiZmVmyXJOqSV69zay2D7qmZWNyT8rMzJLlRsrMzJLlRsrMzJLlRsrMzJLlCyfM7ACevLuEDRsfve0PqnXuSZmZWbLcSJmZWbLcSJmZWbJck7K0Pe962KJ6r6l5uNnYUq1Rrdl1YGxlah7eFfekzMwsWW6kzMwsWZUbKUnLJH1d0uez7adJulbSHZI+Leng9sI06yfnjdlk6tSkXg/cBjw+234H8O6I+JSkC4FzgQ80HJ8ViYKii1e+TVWSeZNK6aRTv1NQdPlQwx9EqjWqHqvUk5J0DPDLwIezbQGnARdlh2wCXtZGgGZ95bwxm1zV4b73AP8T2J9tPwnYHREL2fY24OiiF0raIGmLpC33ThSqWe84b8wmVNpISXoJsDMirh9+uuDQwvGmiNgYEWsjYu2aMYM0a5KkYyVdIek2SbdIen3BMc+XdL+kG7OfC2qew3ljM6WLvClSpSZ1KvBSSeuAxzIYW38PsErS8uxb4THA9kmDSZFLPTNpAXhjRNwg6XHA9ZIui4hbc8d9JSJeMuY5ppY3KZRBfJPDTJW5Sin8g1XTRd4coLQnFRFviYhjIuI44GzgXyLiN4ArgJdnh60HLmkqKLM2RcSOiLghe/wggwsbCofdJjiH88ZmShd5U2SSeVJvBv5I0lYGY+0faSYks+5IOg54LnBtwe5TJN0k6Z8lPbuhUzpvrPe6zJtayyJFxJXAldnjO4GTJg3ArAWrJW0Z2t4YEQeMu0g6FPgs8IaIeCC3+wbgqRGxJxuy+0fg+HGCcd5YTySVN4u8dl9Ob2tQ+blTvf1DGrErItaOOkDSCgaJ9vGI+Fx+/3DyRcSlkv5G0uqI2NV8uKPduxo2/krXZ62ntzWo/NyppudNVVFnjb3PtXqtZ5J542WRbO5k85U+AtwWEe9a4pgnZ8ch6SQGufL97qI0S8u08sY9KZtHpwKvBL4p6cbsubcCTwGIiAsZXNzwe5IWgB8CZ0fEXHdPbe5NJW/cSNnciYivUnJDj4h4H/C+biIyS9+08mbuG6mZLd24RmUt6m0NqkwKNSp7FNekzMwsWW6kzMwsWW6kzMwsWXNXk5rb0kzRPaiGze0HY1XMbA2qTNE9qIa5ZtU696TMzCxZbqTMzCxZbqTMzCxZvapJuWzSorKaVRX+B0rS3NaTulBWs6rCda2R3JMyM7NkuZEyM7NkuZEyM7NkuZEyM7NkdXrhxPXPg0fd99FmS52LL0beWs2G7Vrjix9mWp2LL77S+T03p849KTMzS5YbKTMzS5YbKTMzS5YbKTMzS5YbKTMzS5YbKTMzS5YbKTMzS5YbKTMzS5YbKTMzS5YbKTMzS5YbKZtLks6U9C1JWyWdV7D/MZI+ne2/VtJx3UdplpZp5E2vbnpo1gRJy4D3Ay8EtgHXSdocEbcOHXYucF9EPF3S2cA7gF/rPlp4iEO4IZ45jVMn5Rd1/bRDmLrbp3juaeWNe1I2j04CtkbEnRHxMPAp4KzcMWcBm7LHFwGnS2rg9sVmvTWVvHEjZbNotaQtQz/5NcSPBu4e2t6WPVd4TEQsAPcDT2orYLMEJJk3pcN9kh4LXAU8Jjv+ooj4Y0lPY9CSHgbcALwya13Npm1XRIy6GUjRN7sY45ilT+C8sSk6hId4JvWGR69PIG+KVOlJ/Rg4LSJOAE4EzpR0MoOxxndHxPHAfQzGIs36YBtw7ND2McD2pY6RtBx4AvCDGudw3tis6SJvDlDaSMXAnmxzRfYTwGkMxhxhMAb5skkCMevQdcDxkp4m6WDgbGBz7pjNwPrs8cuBf4mIyt8InTc2g1rPmyKValKSlkm6EdgJXAZ8G9idjTlC8djk4ms3LI5xcu8koZo1I/vv9nXAF4HbgM9ExC2S3ibppdlhHwGeJGkr8EfAAZfblmksb+5bKDrErFNd5U1epUvQI2IfcKKkVcDFQNH1sIWtZURsBDYCaK0malHNmhIRlwKX5p67YOjxj4BXTHiOZvLmWSudN5aELvImr9bVfRGxG7gSOBlYlY05QvHYpJnhvDGbRGkjJWlN9k0QST8DnMGgq3cFgzFHGIxBXtJWkGZ947wxa0aV4b4jgU3ZbOODGIxDfl7SrcCnJP0Z8HUGY5FmNuC8MWtAaSMVEd8Anlvw/J0MZiCbWY7zxqwZ3a7dd/3zQFuW3h9edWaWDf/zjpoxaDn3roYP5Sf/D/mdjd3FYp177dA/75/P4RXSXhbJzMyS5UbKzMyS5UbKzMySldb9pNqY6+s6VyX+mHpsVL1qXK5zVfJaf0ytc0/KzMyS5UbKzMyS5UbKzMyS5UbKzMySldaFE20ouxhjTq4YmJM/05pSdjHGnFxY4Qsjps89KTMzS5YbKTMzS5YbKTMzS9bs16TK5GtWM1K8mZE/w1KVr1nNSI3KNaj0uCdlZmbJciNlZmbJciNlZmbJck0qr6c1qp6EaePYtQY2jpi3tCGBQkpPa1Qp1qBG/XMmGG7r3JMyM7NkuSdlNkTSXwH/HXgY+Dbw6ojYXXDcXcCDwD5gISLWdhmnWUrazBv3pMwe7TLgFyLiOcC/AW8ZcewLIuJEN1Bm7eWNe1JlelqjsvFExJeGNq8BXj6tWCobVa9aStt1rJ7WqJqWQrmwC23mjXtSZkt7DfDPS+wL4EuSrpfUwq1xzXqr0bxxT8pm0WpJW4a2N0bET77TSvoy8OSC150fEZdkx5wPLAAfX+Icp0bEdkmHA5dJuj0irmoofrNpSDJv3EjZLNo1arw7Is4Y9WJJ64GXAKdHROG9XiJie/Z7p6SLgZMAN1LWZ0nmjRupHnJZrD2SzgTeDPzXiHhoiWNWAgdFxIPZ4xcBb+swzMmV1bFmsJjSxZyoGfzYKmkzb1yTMnu09wGPYzAUcaOkCwEkHSXp0uyYI4CvSroJ+BrwTxHxhemEa5aE1vLGPSmzIRHx9CWe3w6syx7fCZzQZVxmKWszb9yTMjOzZLknZWY2pnmtQXXJPSkzM0uWGykzM0tWaSMl6VhJV0i6TdItkl6fPX+YpMsk3ZH9fmL74Zr1g/PGrBlValILwBsj4gZJjwOul3QZ8Crg8oh4u6TzgPMYXCdvZn3Pm/w8KhdfAH8M01DaSEXEDmBH9vhBSbcBRwNnAc/PDtsEXEmKyWa99rzrYUvNycspzHV23tg0rb63foP62nZCmVitmpSk44DnAtcCR2SJuJiQhy/xmg2StgzWhLp3smjNesh5Yza+yo2UpEOBzwJviIgHqr4uIjZGxNrBmlBrxonRrLecN2aTqTRPStIKBon28Yj4XPb0PZKOjIgdko4EdrYVZFLy95cCL6ZnhZw3Q/L3l4K5vceU1VPl6j4BHwFui4h3De3aDKzPHq8HLmk+PLN+ct6YNaNKT+pU4JXANyXdmD33VuDtwGcknQt8F3hFOyGa9ZLzxqwBVa7u+ypLXzB1erPhmM0G541ZM7zihJmZJcsLzJpZOU/utSlxT8rMzJLlRsrMzJLlRsrMzJLlRsrMzJLlRsrMzJLlRsrMzJLlRsrMzJLlRsrMzJLlRsrMzJLlRspsiKQ/kfQ9STdmP+uWOO5MSd+StDW7DbzZ3Gozb7wsktmB3h0Rf73UTknLgPcDLwS2AddJ2hwRt3YVoFmCWskb96TM6jsJ2BoRd0bEw8CngLOmHJNZ6sbKGzdSNotWS9oy9FNwW9iRXifpG5I+KumJBfuPBu4e2t6WPWfWZ0nmjYf7bBbtioi1S+2U9GXgyQW7zgc+APxvILLf7wRek3+LgtfGeKGaJSPJvHEjZXMnIs6ocpykDwGfL9i1DTh2aPsYYHsDoZkla1p54+E+syGSjhza/B/AzQWHXQccL+lpkg4GzgY2dxGfWYrazBv3pMwe7S8lnchgGOIu4LUAko4CPhwR6yJiQdLrgC8Cy4CPRsQt0wrYLAGt5Y0bKbMhEfHKJZ7fDqwb2r4UuLSruMxS1mbeeLjPzMyS5UbKzMyS5UbKzMyS5UbKzMyS5UbKzMyS5av7zKzcho3TjsDmlHtSZmaWLDdSZmaWLDdSZmaWLNek6oqihXzNbKTfmY2a1sbczStcqmufe1JmZpYsN1JmZpas0kYqu8viTkk3Dz13mKTLJN2R/S66C6PZ3HLemDWjSk/qY8CZuefOAy6PiOOBy7NtM/upj+G8MZtY6YUTEXGVpONyT58FPD97vAm4Enhzg3GZ9Vrv88ZXBFgixq1JHREROwCy34c3F5LZzHLemNXU+iXokjYA2YWbT2n7dGYzwXljk1izq35n+LXthDKxcXtS9yze0z77vXOpAyNiY0SsjYi1sGbM05nNBOeNWU3j9qQ2A+uBt2e/L2ksIiulOPA5zzHuhXTzZg5qUB/ccOBzr53wz/bk3vZVuQT9k8DVwDMkbZN0LoMke6GkO4AXZttmlnHemDWjytV95yyx6/SGYzGbGc4bs2Z4xQkzM0uWF5gt42KP9d00CiUzsqBsXfkaVZ5rVvW5kTIbIunTwDOyzVXA7og4seC4u4AHgX3AwuAqPLP51GbeuJEyGxIRv7b4WNI7gftHHP6CiNjVflRmaWszb9xImRWQJOBXgdOmHYtZX7SRN26kZkR+7tScl9JWS9oytL0xIupWA/4zcE9E3LHE/gC+JCmAD47x/tWtvhd+xcWMNuTnTk06b6pMWc2qzL2fayaOJSSZN26kbBbtGjXWLenLwJMLdp0fEYsTbM8BPjniHKdGxHZJhwOXSbo9Iq4aP2SzqUsyb9xI2dyJiDNG7Ze0HPgV4Hkj3mN79nunpIuBkwA3UjazppU3nidldqAzgNsjYlvRTkkrJT1u8THwIuDmomPN5kgreeOeVN6MFHNco5rI2eSGLCQdBXw4ItYBRwAXD2rELAc+ERFf6DzKlMzIvKiua1QzppW8cSNllhMRryp4bjuwLnt8J3BCx2GZJa2tvPFwn5mZJcuNlJmZJcvDfXNSrHGNyho1IzWoMq5RTZ97UmZmliw3UmZmliw3UmZmlqzZq0m52FJJvkaV549xzsxJjWlS+RpVnmtWzXNPyszMkuVGyszMkuVGyszMkuVGyszMktXthRPPux62uCLfB2UXVkxsybvW2AHW7PKFDT1RdmHFpHZ9pd33T5F7UmZmliw3UmZmliw3UmZmliw3UmZmliw3UmZmliw3UmZmliw3UmZmliw3UmZmliw3UmZmlqyJGilJZ0r6lqStks5rKiizWea8Matu7EZK0jLg/cCLgWcB50h6VlOBmc0i541ZPZOs3XcSsDUi7gSQ9CngLODWJgIzm1HOmzHdEM+bdggJuG3aAXRukuG+o4G7h7a3Zc+Z2dKcN2Y1TNJIFS1nfsDa2ZI2SNoiaQv3TnA2s9lQP2/uW+ggLLM0TdJIbQOOHdo+BtiePygiNkbE2ohYy5oJzmY2G+rnzRO7vaOOWUoUMd6NgyQtB/4NOB34HnAd8OsRccuI19wLfAdYDewa68TdcYzNWCrGp0ZE6dcWSV/I3qOOXRFxZs3XdMJ5k4Q+xzh3eTN2IwUgaR3wHmAZ8NGI+POKr9sSEUnf9s4xNqMPMXbNeTNdjrFfJhpHiIhLgUsbisVsLjhvzKrzihNmZpasaTVSG6d03jocYzP6EGNf9OGzdIzN6EOMnZioJmVmZtYmD/eZmVmyOm2kUl1YU9JHJe2UdPPQc4dJukzSHdnvJ045xmMlXSHpNkm3SHp9anFKeqykr0m6KYvxT7Pnnybp2izGT0s6eFox9pHzZuz4ks+ZLB7nzQidNVKJL6z5MSA/P+A84PKIOB64PNuepgXgjRHxTOBk4Pezzy+lOH8MnBYRJwAnAmdKOhl4B/DuLMb7gHOnGGOvOG8m0oecAefNSF32pH6ysGZEPAwsLqw5dRFxFfCD3NNnAZuyx5uAl3UaVE5E7IiIG7LHDzJYafJoEoozBvZkmyuynwBOAy7Knp/6Z9kzzpsx9SFnwHlTpstGqm8Lax4RETtg8B87cPiU4/kJSccBzwWuJbE4JS2TdCOwE7gM+DawOyIWF6BL/d89Nc6bBqScM+C8GaXLRqrSwpo2mqRDgc8Cb4iIB6YdT15E7IuIExmsSXcS8Myiw7qNqtecNxNKPWfAeTNKl41UpYU1E3KPpCMBst87pxwPklYwSLaPR8TnsqeTixMgInYDVzKoBazK1qyD9P/dU+O8mUCfcgacN0W6bKSuA47Prlg5GDgb2Nzh+evaDKzPHq8HLpliLEgS8BHgtoh419CuZOKUtEbSquzxzwBnMKgDXAG8PDts6p9lzzhvxtSHnAHnTamI6OwHWMdgBehvA+d3ee6SuD4J7AAeYfDN9VzgSQyu/Lkj+33YlGP8JQbd/W8AN2Y/61KKE3gO8PUsxpuBC7Lnfxb4GrAV+AfgMdP+N+/Tj/Nm7PiSz5ksTufNiB+vOGFmZsnyihNmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1JmZpYsN1I1SDpOUkhavsT+t0r6cNdxmaXEeWJNmutGStJdku6RtHLoud+WdOU47xcRfxERv91YgAUkHSzpoiz2kPT83H5Jeoek72c/fylJbcZks20e80TSiZKul/RQ9vvENuO1pc11I5VZDrx+2kHU9FXgN4F/L9i3AXgZcALwHOAlwGu7C81m1NzkiaSDgUuAvweeCGwCLsmet465kYK/At4kaVWN17xG0nZJOyS9cfFJSX8i6e+zx4tDHuslfVfSLknnDx17kqQtkh7IvqW+q8qJI+LhiHhPRHwV2FdwyHrgnRGxLSK+B7wTeFWNv82syDzlyfMZNMrviYgfR8T/AQScVuNvt4a4kYItwJXAm2q85gXA8cCLgPMknTHi2F8CngGcDlwg6ZnZ8+8F3hsRjwd+DvjM4gskfUPSr9eIZ9izgZuGtm/KnjObxDzlybOBb0REDO3/Bs6jqXAjNXAB8AeS1lQ8/k8jYm9EfBP4W+CckmN/GBE3MUiEE7LnHwGeLml1ROyJiGsWXxARz4mIT4zxdwAcCtw/tH0/cKjrUtaAecmT/L7F/Y8b81w2ATdSQETcDHweOK/iS+4eevwd4KgRxw6Phz/EIAEAzgX+A3C7pOskvaTiucvsAR4/tP14YE/uW6FZbXOUJ/l9i/sfbOjcVoMbqZ/6Y+B3gKMrHHvs0OOnANvrniwi7oiIc4DDgXcAFw1fPTWBW/jpt1Cyx7c08L5mMB95cgvwnNzow3NwHk2FG6lMRGwFPg38YYXD/5ekQyQ9G3h19rpaJP2mpDURsR/YnT1dVOAteu1jJD022zxY0mOHEurvgD+SdLSko4A3Ah+rG59ZkTnJkyuzc/xh9h6vy57/l7rx2+TcSD3a24Aq39L+FdgKXA78dUR8aYxznQncImkPg+Lw2RHxIwBJt0j6jRE3NsV+AAAgAElEQVSv/RbwQwbfZr+YPX5qtu+DwP8HfBO4Gfin7Dmzpsx0nkTEwwwuT/8tBg3ja4CXZc9bx+RShZmZpco9KTMzS5YbKTMzS5YbKTMzS5YbKTMzS5YbKTMzS1bh/V7actBhiuXHdHnGND3CimmHMH3bFogf7C9dqkk6M2BXzTe//osRceaYkSXnoMMUy44tP25RtLAAVtT8OjtODI3HXeUN6560/D/Z+urEcPe+ucubThup5cfAEf/U5RnTtI0jph3C9P3yPRUP3MVgbdM6tLrmC5K27Fh40herH/9IC9+BHq55k4pxYih7zf66Az9VgqgbaN0Pooo6Mfy371c8cHbyxsN9ZmaWrE57UmZW34pH4MgdS+9vo9dSd/8471fWKSn9u+oGUXTCsvdo+oOqYtRrHnmg2nsoYPkjNc9b7/CuuCdlZmbJck/KzGzWKAZd8DrckzIzM6unUk9K0irgw8AvAMFgVeBvMVh6/zjgLuBXI+K+Ue9z6F445ZpRR7Tv6pOne/5kXHPKdM+/98vTPX8Hmsqb5Quw5t7x4xjnIrf8dr6cU7eU00bd7AB1gwR46JB6ryk7R14TH8SwhW9P9voeqtqTei/whYj4eQY3B7uNwd05L4+I4xksxV/1bp1m88J5Y9OxONxX5ydRpY2UpMcD/wX4CAzutRIRu4GzgE3ZYZsY3H/FzHDemDWlynDfzwL3An8r6QTgeuD1wBERsQMgInZIOrzoxZI2ABsADklyqphZKxrLm8cWHmE2ggIOno17NFZppJYDvwj8QURcK+m91BiiiIiNwEaAw35OU7/DYt2aWG9rWNOuOVljebP6aYqjtlc/cZWyx6Q1prJSzt7cfXvzx0P5CFP+PX+c/7smDbLomHzgdT+Y/OvzimKoU6d6pIUVLxJXpSa1DdgWEddm2xcxSL57JB0JkP3e2U6IZr3kvDFrQGkjFRH/Dtwt6RnZU6cDtwKbgfXZc+uBS1qJ0KyHnDc2VTN04UTVybx/AHxc0sHAncCrGTRwn5F0LvBd4BXthGjWW84bswlVaqQi4kZgbcGu05sNx2x2OG/MJtfpskgr98DJV7d7jqavF8hfaJHshRRtXyjR9B++p+X/EGbIiodh1IUT4yzUWveag/w58qNDRRdG1HX/E2q+oO7VH0VBll3hMemFFVU+mLIJw4/aV/Eii3GWRUqUl0UyM7NkeYFZM7NZ456UmZlZ+2auJ1VW85qJOa5N/BHJFtcsb/kCrN619P5x7rtXNgf1kIcevZ0vrdRdzKCozJL/oj/xF/9xJtrm/7Ddq0bvL9seZ5HbOovQLszc/7JLzd9fbGY26zzcZ2Zm1j73pMzMZs0M9aQ6baQO3QunTDA95uoGSjH5mlXd8k7RArWtl3f6UoOq8w+090PtxTFjym56WFaTKqoH5UspZTWoSVUpxZTVqGovMFtllduyGlN+8lbdeVVV6k11bqQ4hzUpD/eZmVmy5q9Ztn5ZuReeU7P77cUsbN7N0P2k3JMyM7Nk9aonVaWeVbduNWmNKlmT1qCaKABaI+rOk6pSBsl/yS5b/6+ufOmmqIZfu65fd0JYlQ+irMaU35+vUdVd26+oLlbn73JNymz2SXqGpBuHfh6Q9IbcMc+XdP/QMRdMK16z2lq4n9S08mb+mmWbexHxLeBEAEnLgO8BFxcc+pWIeEmXsZmlalp540bK5t3pwLcj4jvTDsSsMe3Pk+osb2aukcrXrbqoUfXmnlOjzFYNarWkLUPbGyNi4xLHng18col9p0i6CdgOvCkibmkyyKqW74NVu5feX1Z6qbJu3qTzovLnXLm3/P3zdbHG/586ztp9dbfLalZVYqizdt++Vv+XnWTezFwjZQbsioiiO+I+SnZb95cCbynYfQPw1IjYI2kd8I/A8c2GaZaUJPPGF07YPHsxcENE3JPfEREPRMSe7PGlwApJq7sO0GwsLVw4MaTTvHEjZfPsHJYYspD0ZEnKHp/EIFe+32FsZqnqNG86He5bubd47bulNFHbmbRG1Rt1P6w2Pohak8wKxuY7JOkQ4IXAa4ee+12AiLgQeDnwe5IWgB8CZ0dETCPWZfvgCfcvvb+JeVJlX6TL3rNs7b/8/qJjJp43Nc4HUVZjanreVJV5UiNrUsuW3jespQsnppE3rknZXIqIh4An5Z67cOjx+4D3dR2XWcqmkTce7jMzs2S5J2VmNmt8P6luVKlf1S3F1K1RTWVtvz4sINiHGGfEstw8qbr3jxpnDlTd8k5ZDapoelCZAxbxLvm7S+/LVOWmVk3Pmyo7viiGUX9H1ZrUDPFwn5mZJSvpnpSZmY1hhob73JMyM7Nk9b4nNRPr5nVh0nlRrkFNzUH7i+cZLcqXNJr4Al1W3smXXvLx7V5Vfo7Ov+gX1aTK/tCyAl/dGtQ4a/cNb++v2K/Qft+Z18zMrG1upMzMLFm9H+4zM7OcGbpwolIjJeku4EFgH7AQEWslHQZ8GjgOuAv41Yi4r50wzfrHeWM2uTo9qRdExK6h7fOAyyPi7ZLOy7bf3Gh0Nj2+UKIpE+fNpF+K8zcgLJKv1ZctQFsWT35/lRr+xAvM1t1fdMyk220vMBtaet+wGepJTVKTOgvYlD3eBLxs8nDMZp7zxqyGqo1UAF+SdL2kDdlzR0TEDoDs9+FFL5S0QdIWSVvunTxesz5pJG92TeUGIWZpqDrcd2pEbJd0OHCZpNurniAiNgIbAdZKTjerZ+UeOLnGTcgAri4/pCON5M0vLnPeWE0zNNxXqZGKiO3Z752SLgZOAu6RdGRE7JB0JLCzxTgrqzu5t5c3RfSM5V5oKm8O2g8rhybLPlwyebfKgrJ1a05lNaa6Nauqx9TSRI2q7DV1F7EdZ5HbA7aHzzl/s4ZK/2JJKyU9bvEx8CLgZmAzsD47bD1wSVtBmvWN88asGVV6UkcAF2e3rV8OfCIiviDpOuAzks4Fvgu8or0wzXrHeWPTM0/DfRFxJ3BCwfPfB05vIyizvnPemDXDK06Y9czBuS/IZTWqKqWYcWpKdUzlS32+flRknBpSHWPN1aoQd5kZ6knNXxXOzMx6w42UmZkly8N9ZmazZoaG+9xI2cCcrdVXtPhrbr+A9wLrgIeAV0XEDV3HWUW+RtXGTRDz2q5h9UbdmtWkNa4pm0beuJGyeZZf/HXYi4Hjs5//BHwg+22WvoOizTvzdpo3rkmZFTsL+LsYuAZYla0QYWZLazxv3EjZLFq9uDhr9rOh4JiixV+HHQ3cPbS9LXvObFYlmTce7rOBk3OLGPa7RrUrP1Ze4IDFXyPiqqH9RTfuSXKh1/w8qbw2yiB1pxP1vBSztLo31io6vo0PZ7wLJ5LMG/ekbC4NL/4KLC7+OmwbcOzQ9jHA9m6iM0vTNPLGjZTNnRGLvw7bDPyWBk4G7l+8D5RZ8hZ7UnV+yt5ySnnj4T6bR0st/vq7ABFxIXApg8totzK4lPbVU4rVLBVTyRs3UjZ3Riz+euHQ4wB+v8u4lhKCh0dk6jjLzdW9LVIvVbkEu+6NtOq+X9OvWZheWXRaeeNGysxs1szQihOuSZmZWbLcSJmZWbJmbrjv6pNrHt/H6UCnXHPgc3X/cOuN0OgaUVk96aFD6r+mbLuJmlbjda+6c5bGec98natuTeuQhw48R9E/0FL2VaxJebjPzMysfTPXkzIzm3vuSZmZmbXPPSkrNltr+fXa/oNGly3q1o+Kjsm/f/41+f2T1rCK1K5RNVGDKqsZ7V05en9ZcW7l3tHvV/Seoz6Ih/cvvW9GuZGytB26F065uvy4Ye9uJxSz3lCr95PqlIf7zMwsWe5JmZnNGl84YWZm1r7e96TansM6lesF2rhoIV/XqTuL2RdSTE3ZhRPjTKzNv1++nj/phRJVYmh8Mm9ZDaaoZ1E2Obdsf/7CiHH+qDpXmTzgCyfMzKzvPNxnZmbWPjdSZmaWrKSH+9qoN/VyQdkU5WtURVy3asS+ZbB7VfXjqywwW1aDqluzGqcu1sQE4JGq3MBw0sm7eWU1qaIY6qzOe+++0e+/yMN9ZmZm7Uu6J2VmZmOYx56UpGWSvi7p89n20yRdK+kOSZ+WNGnn3GzmOG/MJlOnJ/V64Dbg8dn2O4B3R8SnJF0InAt8YNQb7F0JV//HseIc29zUoPI3Qiwr6E06b6qKUXWrb+5det9smThv9i2D+59Q/YT5kkbRmqZldat8DaysZlV3XlUj6t5wMD+nCQ4MvIl5T6NiKHq/soLesGUVa1IzpFJPStIxwC8DH862BZwGXJQdsgl4WRsBmvWV88amZnG4r85PoqoO970H+J/A4nTnJwG7I2Ih294GHF30QkkbJG2RtOW+dD8HszY0kje7nTc2x0obKUkvAXZGxPXDTxccGkWvj4iNEbE2ItY+sY0uv1lNko6VdIWk2yTdIun1Bcc8X9L9km7Mfi6oeY7G8maV88bqaqEn1UXeFKlSkzoVeKmkdcBjGYytvwdYJWl59q3wGGD7pME0YdLSyjhTe9peP7ATZfdsmq3i3gLwxoi4QdLjgOslXRYRt+aO+0pEvGTMczSWNwsl86SaWDevbB5UviZWdnx+/6SlmErq1qig3g0Hqyhb62+cCWPDli0sva99XeTNAUp7UhHxlog4JiKOA84G/iUifgO4Anh5dth64JKmgjJrU0TsiIgbsscPMriwoXDYbYJzOG9spnSRN0Ummcz7ZuCPJG1lMNb+kWZCMuuOpOOA5wLXFuw+RdJNkv5Z0rMbOqXzxtqn/YNeXZ2fOm/fYd7UmswbEVcCV2aP7wROmjQAsxaslrRlaHtjRGzMHyTpUOCzwBsi4oHc7huAp0bEnmzI7h+B48cJxnljPZFU3ixKasWJaZQ9erG83Dj3cqo7b6r0/Sqs1Zc3vTrWrohYO+oASSsYJNrHI+Jz+f3DyRcRl0r6G0mrI2JX8+GOtrAcdq1eev84y8mVzWuqu3Zf3ftRLfXcSJPOi5q0HjROTFVqXnUWLVzeak0qybxJqpEy60I2X+kjwG0R8a4ljnkycE9EhKSTGAyNf7/DMM3G18KySNPKGzdSNo9OBV4JfFPSjdlzbwWeAhARFzK4uOH3JC0APwTOjojCy8XN5sRU8saNlM2diPgqxXOWho95H/C+biIya1gLPalp5U2njdSeldOfbjNpDaq3c6LyNaq8Nv6wUXWsO/c0f74ZtbAc7l2z9P6yMso486QmrVGNs35g4+v9jTNPqu57lNXBqkwGq/OHt1uTSpLvJ2VmZsnycJ+lbeWe8a4sNJtn83g/KTMzs671vifVi3lOTSu6T9OkH8Q0alZWycKK9mtSk9aoxrmfVJW5VCOVrZNX5X5SZSadB5X/I4vU+QdcUbEm5Z6UmZlZ+9xImZlZsno/3GdmZjkzNNzXaSO199D0a0i9Lb2Ms75fHWU1q7ouGaM+MKceWQE7jqx+fNlafovvOcl2vtSSP+c4a/dViXuksjlPRcpqTvkg8+co+0NX7X70dtEfWacYNyMNTx0e7jMzs2R5uM/MbNYoat8jKlXuSZmZWbLckzIzmzW+cKI/enshxKSKJvyOkvoVLXPs4RWw/ajxX1+lLl92EcOkF1qMc9PD0rjL/idcdtFD0TFlF05Uec86r69rxWwM4dXh4T4zM0vWzPekzMzmzgwN97knZWZmyer+pofzWiNKXd0a1qRW+qaHVT0yYU1qqfec5Pi69aRxJhjXVvcGhFWCmLRwNvEflVO1d+SelJmZWfvcSJmZWbJ84YSZ2aw5aHaG+9xImSXukYObr0lNqolSS+M1qLL9TfxPu+kaU10Hz0bDU4eH+8zMLFnuSZmZzRpf3WfWb5LOlPQtSVslnVew/zGSPp3tv1bScd1HaZaWaeRNpz2pR1jBNo7o8pRJOoZt0w5h6u6Z4rklLQPeD7wQ2AZcJ2lzRNw6dNi5wH0R8XRJZwPvAH6t+2hhPwexh0PKD1yUwjfoFmI4iP2TvcGKHx/4VN0wC95jUnViqHxbxxZ6UtPKG/ekbB6dBGyNiDsj4mHgU8BZuWPOAjZljy8CTpekDmM0S81U8saNlM2i1ZK2DP1syO0/Grh7aHtb9lzhMRGxANwPPKmtgM0SkGTelA73SXoscBXwmOz4iyLijyU9jUFLehhwA/DKrHU1a8wKHuGImsOj22BXRKwdcUjRN7sY45ilT+C8sSlbxr5ax+9LIG+KVKlJ/Rg4LSL2SFoBfFXSPwN/BLw7Ij4l6UIGY5EfmCQYs45sA44d2j4G2L7EMdskLQeeAPygxjmay5uyG0o1UXuoO+eo7v6iW5lP+p4lh1d5eT6stj+GcQy/53em+3Wmi7w5QOlwXwwsrga6IvsJ4DQGY44wGIN82SSBmHXoOuB4SU+TdDBwNrA5d8xmYH32+OXAv0RE5W+EzhubpsXrJur8VNB63hSpVJOStEzSjcBO4DLg28DubMwRiscmF1+7YXGMkx9MeHWOWQOy/25fB3wRuA34TETcIultkl6aHfYR4EmStjLo/RxwuW2ZxvLmvnrDNmZt6Cpv8ipdgh4R+4ATJa0CLgaeWXTYEq/dCGwE0HMOnqhFNWtKRFwKXJp77oKhxz8CXjHhOZrJm194rPPGktBF3uTVmicVEbslXQmcDKyStDxrXYvGJg+091C45pRx4qyu6/sizahTrmn3/b+8t/yYWTFx3iwsh3vXNBtUfnynqGY06vi620Vr3pXFUDIGVXbKKm93yEPlx0xyzrLX1/W9hfJjYLxpUj+qH04nSof7JK3Jvgki6WeAMxh09a5gMOYIgzHIS9oK0qxvnDdmzajSkzoS2JTNNj6IwTjk5yXdCnxK0p8BX2cwFmlmA84bmxrtT2PhkSaUNlIR8Q3guQXP38lgBrKZ5ThvzJoxe6ugT1rzmtGaVts1JmvRpDeUqjJHadLt/DmqFHvyx+TfI1/Hyq2b10YI+WNW7h29v+52kTo9nltnpHdUx+w1UmZmc26G7tThtfvMzCxdbqTMzCxZ3Q737VkJV588/uu7KKzka1o9qVF1XXOa9GO5ek/5MZZ5ZAXsOHLp/eMsKFe3mFJWnMkXd/KK9udrTg8fPPocJcr+hKIQ8n9W/pi2a1ZFRh1T9SPxcJ+ZmVkHfOGEmdmMcU/KzMysA26kzMwsWf0a7hvnootJryhI8EKKNi6SSODPsqU8UnOB2SrV+rLVV/NXDBwwsbaFGy0WLUI7Qt0/IX+RA8AT7q/3mrILKcrmJ1e5hmWU5TUunChb7LYv3JMyM7Nk9asnZWZmpXzhhJmZWQdmvydVVseqW+ApWsC2hwWdHoY8vxaWw67V1Y+vUpOqW4PKb5dN3i07H8BDh4w+pmSB2bJTVJnMW7cGla9hlS1qW2WR21o1qYo3PZwls99ImZnNGQ/3mZmZdcCNlJmZJavb4b69h8LVE96UMO+USVc6zdWsErw7YBMhTVqDmvRjzvtQwZyVIofurX/uf6gfTtoWWpgnlS++lNWg8ou/ls1pqjJBaMKJPJPehxHq16Dyx+f3112Atui5UcN0VWtSHu4zMzPrgC+cMBsi6a+A/w48DHwbeHVE7C447i7gQWAfsBARa7uM02yUrntSbeaNe1Jmj3YZ8AsR8Rzg34C3jDj2BRFxohsos/bypv89qbIaV92CRg9qVFWkVoPqi4j40tDmNcDLpxXLTyysGF2TKivGFH2lnrTmVBZD2byrIiXH1L234zhzlMrmTU261l9RXaysNDcc44pE50m1mTfuSZkt7TXAPy+xL4AvSbpe0oYOYzIrtTjcV+enQY3mTf97UmYHWi1py9D2xojYuLgh6cvAkwted35EXJIdcz6wAHx8iXOcGhHbJR0OXCbp9oi4qqH4zaYhybxxI2WzaNeo8e6IOGPUiyWtB14CnB4RscR7bM9+75R0MXAS4EbKknDQeLfqSDJvOp4ntbJ47btFbSwol69ZtVFsSeCeUynWoEaV81Y2f7pGSDoTeDPwXyOicIE6SSuBgyLiwezxi4C3tRbUvmVw/xOW3l820aZKTapuDSovX4zJr8tXdDOnfB2sprp/dtHHMGlNadXu0fvL3q9qnIuW7Vt63zS1mTeuSZk92vuAxzEYirhR0oUAko6SdGl2zBHAVyXdBHwN+KeI+MJ0wjVLQmt54+E+syER8fQlnt8OrMse3wmc0GVcZnV0PU+qzbxxT8rMzJKVVk9qVL0Kmqn1dFGj6oFJ/+yeTh/rp33LYPeqn27XmVgD9e/9VEX+PfeuHL2/qOZVd73ACctmRR9b2UdVt8ZUd22/oudGTXtLtSbVprQaKTMzm5gXmDUzM+tAaSMl6VhJV0i6TdItkl6fPX+YpMsk3ZH9fmL74Zr1g/PGpmnKK040qspw3wLwxoi4QdLjgOslXQa8Crg8It4u6TzgPAbXybenqGbV9pykKazl18YpXIPqXHN5s/+gR887ys9BqjPRpqr8e+TrRfkYympQRXOiJp2blVO2ZOE49aC6NaiymlaVW3uN+uc7aP/S+2ZVaU8qInZExA3Z4weB24CjgbOATdlhm4CXtRWkWd84b8yaUevCCUnHAc8FrgWOiIgdMEjIbC2motdsALKFBJ8yfqRmPTVx3hx0TCdx2uyYywsnJB0KfBZ4Q0Q8UPV1EbExItYO1oSqcQtssxnQSN4c9KT2AjRLXKWelKQVDBLt4xHxuezpeyQdmX0bPBLY2VaQIyWwbt4scg1qco3lTb4mlddFTaqs5lR3u0jNGlXd+0sVHV+3jlW2XVaDKlrCsM45qtak5qonJUnAR4DbIuJdQ7s2A+uzx+uBS5oPz6yfnDdmzajSkzoVeCXwTUk3Zs+9FXg78BlJ5wLfBV7RTohmveS8MWtAaSMVEV8FtMTu05sNx2w2OG9smuZquM/MzGxaZm/tvroXUvRwwVlfGzJnQpNNfB1ncde6i792cKFEmboXUhQ91/aFElVuvDjq71Dh/W5n2+w1UjZTVu6p/73hH9oJxaw3NN7t45Pk4T4zM0uWe1JmZjNmli6ccCNllroQPFKwQGtVRbWffE2p7MaIDdePulDlf9J1a1B17zdZZf8B5xhVkxr99jPJw31mZpYs96TMzGbMLA33uSdlZmbJck/KAC8o22v5elXZDQsL36PmPKdJj58RTSxy2wb3pMzMzDrgRspsiKQ/kfQ9STdmP+uWOO5MSd+StDW7DbzZ3GozbzzcZ3agd0fEXy+1U9Iy4P3AC4FtwHWSNkfErV0FaDbKlIb7WskbN1IGwNUnP3rbNaqRTgK2RsSdAJI+BZwFTKeRWjHh5J1GYpiRAkhNdUttRcfP0Uc3Vt54uM9m0WpJW4Z+NtR8/eskfUPSRyU9sWD/0cDdQ9vbsufMkrDYk6rzQ6J5456UzaJdEbF2qZ2Svgw8uWDX+cAHgP8NRPb7ncBr8m9R8No5XJ/aZkySeeNGyuZORJxR5ThJHwI+X7BrG3Ds0PYxwPYGQjNL1rTyxo2U2RBJR0bEjmzzfwA3Fxx2HXC8pKcB3wPOBn69oxDbUXdCT91CyhQKL01MD5t0elf+9UUfQ9k5Rq3ltxTt7/ZWHW3mjRsps0f7S0knMhiGuAt4LYCko4APR8S6iFiQ9Drgi8Ay4KMRccu0AjZLQGt540bKbEhEvHKJ57cD64a2LwUu7Sous5S1mTdupMzMZswsLYs0e43UyTXvNV733uRmXVPA8qECQxMLxuULFl3fWKnqMTV0sTxg/jZcdWtak8YYC5O9vo9mr5EyM5tzs9ST8mReMzNLlhspMzNLlof7euiaUw58rm4prozX8kvIQfvhkIeW3l9WP1q598DX5N+vrOZUt2ZVFuO4xwxpYt28SWtK+e2HDnn0dpU/adQ/bf4c+/eUvx94uM/MzKwT7klZ0g7dCye7F2dWi3tSZmZmHeh/T6rpYsyMujpXx/L0sB6pW5MqqycVPVdWo5p0u0jLX/WrzFEqm/e0d+Wjt+v+2fkaVZU62qiPZX9J/WoW9b+RMjOzR/Fwn5mZWQdKG6nsLos7Jd089Nxhki6TdEf2u+gujGZzy3lj1owqPamPAWfmnjsPuDwijgcuz7bN7Kc+hvPGpmTM28cnqbQmFRFXSTou9/RZwPOzx5uAK4E3NxhXsSYukqh7xcAUZrGmOJE2H1ORFOJMRaN5c9D+R0/IrTtxtuiii/xz+Qm/ZfvLJgzXXZC2SM3X1F38tei5/IUO+Y+h7Phx1t0tuzHi8Pb+H4x+/1k07oUTRyzehTEidkg6vMGYzGaV88Y6cVB0e2feNrV+4YSkDZK2SNoC97Z9OrOZ8Ki82TeHX5/NMuM2UvdIOhIG97YHdi51YERsjIi1EbEW1ox5OrOZMF7eLDusswDNUjPucN9mYD3w9uz3JZVetXIv/MeOZ5HOyazV/KKzZeW7Nib3VqlbLdr7zcnP10Pj5c2yffCE+3+63cRND8tqUKt2j95fNvm3bH/VOGsoq0Hla1ZwYE0p/7Hk99eVP2fREFzZArPDH8u+7dXOO1fzpCR9ErgaeIakbZLOZZBkL5R0B/DCbNvMMs4bs2ZUubrvnCV2nd5wLGYzw3lj0zRXPSkzM7Npmb21+yYtrowz2WcGFrn1ArQJW/4IrKlxZWyVWk/deVJ1j68Sw4TXSJct1lo2pwkO/DPyC8pOGlPZHKiic47qAT1yx3hx9dnsNVJmE5D0aeAZ2eYqYHdEnFhw3F3Ag8A+YGFw9apZGroe7mszb9xImQ2JiF9bfCzpncD9Iw5/QUTsaj8qs7S1mTdupMwKSBLwq8Bp047FrC/ayJt+NVIulFRWd95UXr5GlZf4P8XqwQonP7ExIjbWfI//DNwTEUtVAQL4kqQAPjjG+1eXnyeVV1b/Kar9lM1rartmVRRDyfhU3RpUlbX72q5BVZknVWfa275l1eIac7gvybzpVyNlVs2uUWPdkr4MPLlg1/kRsTjB9hzgkyPOcWpEbM/W37tM0rERTw4AAAZ5SURBVO0RcdX4IZtNXZJ540bK5k5EnDFqv6TlwK8AzxvxHtuz3zslXQycBLiRsiS0ceHEtPLG86TMDnQGcHtEbCvaKWmlpMctPgZeBNxcdKzZHGklb7rtSa3cM/1iRg9vetTE/aUmrVHlldWsyuy5c7LXt+xsckMWko4CPhwR64AjgIsHNWKWA5+IiC+0Fs3yhdHzpMZZA69sXlPb86aKnqv51b9uPahonlTd3kbZOcrmZlX5GEbtX0h77KuVvEn7Tzabgoh4VcFz24F12eM7gRM6Dsussmksi9RW3ni4z8zMkuWelCVt5d7pjxCb9Y1m6M68s9dItV1zSmCdvqL7NtX9s/M1qrwE/kxbtOIROHJH9eOr/N+p7v2fmp43VXSOuvtzqqzVV/cUdedeldWgJl3CsGyu2CzycJ+ZmSVr9npSZmZzzveTMjMz60C3PalD9/ZvnlJPijNNzKUaVlazmtTeS8qPscyKknlSBxw/xvyjsnlSZe9ZVoNqYZ5Ufo5S/uVV1u7bvWr0MWXvWXYvqHy9aZx/imELc1iT8nCfmdmM8XCfmZlZB9yTMjObMe5JmZmZdcA9qZ5cGFFX0YTfYX27fmWurXgYjto+wesnXNW0aDt/RUDZZOCi8+VfU3OJhLoTa4uUXRiR3193su6kF0rkPewLJ8zMrO883GdmZtYB96TMzGbMLPWkur/p4YzWgPqmrGbVtj0ry4+xzIpHJqtJLfWedfZPenyVlVUbntybN85NDyf9GOruL+MFZs3MzBLi4T4zsxkzS/eTck/KzMyS5Z6UWeqWL8DqXT/d7qIiXvccZV/bm54wRHl9pmzOUxV1eyNt/9PMY03KjZSZ2YyZpav7Jhruk3SmpG9J2irpvKaCMptlzhuz6sbuSUlaBrwfeCGwDbhO0uaIuLWp4MxmjfPGujBLPalJhvtOArZGxJ0Akj4FnAU42Ups45hph5CAe6YdwLTUz5uDW5gn1UP7SwZ+flxSrynbX6jC+n+dWrF/2hF0bpLhvqOBu4e2t2XPmdnSnDdmNUzSk1LBc3HAQdIGYAMARy+b4HRmM6F+3jyl5Yhs5szScN8kPaltwLFD28cAB4xJRMTGiFgbEWs5zNOybO7Vz5s1ncVmlhxFHPAlrtoLpeXAvwGnA98DrgN+PSJuGfGae4HvAKuBXUsdlwjH2IylYnxqRJT+71fSF7L3qGNXRJxZ8zWdcN4koc8xzl3ejN1IAUhaB7wHWAZ8NCL+vOLrtkTE2rFP3AHH2Iw+xNg15810OcZ+mWgyb0RcClzaUCxmc8F5Y1adi0RmZpasaTVSG6d03jocYzP6EGNf9OGzdIzN6EOMnZioJmVmZtYmD/eZmVmyOm2kUl1YU9JHJe2UdPPQc4dJukzSHdnvJ045xmMlXSHpNkm3SHp9anFKeqykr0m6KYvxT7Pnnybp2izGT0squdG3DXPejB1f8jmTxeO8GaGzRmpoYc0XA88CzpH0rK7OX+JjQH5+wHnA5RFxPHB5tj1NC8AbI+KZwMnA72efX0px/hg4LSJOAE4EzpR0MvAO4N1ZjPcB504xxl5x3kykDzkDzpuRuuxJ/WRhzYh4GFhcWHPqIuIq4Ae5p88CNmWPNwEv6zSonIjYERE3ZI8fBG5jsOZbMnHGwJ5sc0X2E8BpwEXZ81P/LHvGeTOmPuQMOG/KdNlI9W1hzSMiYgcM/mMHDp9yPD8h6TjgucC1JBanpGWSbgR2ApcB3wZ2R8RCdkjq/+6pcd40IOWcAefNKF02UpUW1rTRJB0KfBZ4Q0Q8MO148iJiX0ScyGBNupOAZxYd1m1Uvea8mVDqOQPOm1G6bKQqLayZkHskHQmQ/d455XiQtIJBsn08Ij6XPZ1cnAARsRu4kkEtYFW2Zh2k/++eGufNBPqUM+C8KdJlI3UdcHx2xcrBwNnA5g7PX9dmYH32eD1wyRRjQZKAjwC3RcS7hnYlE6ekNZJWZY9/BjiDQR3gCuDl2WFT/yx7xnkzpj7kDDhvSkVEZz/AOgYrQH8bOL/Lc5fE9UlgB/AIg2+u5wJPYnDlzx3Z78OmHOMvMejufwO4MftZl1KcwHOAr2cx3gxckD3/s8DXgK3APwCPmfa/eZ9+nDdjx5d8zmRxOm9G/HjFCTMzS5ZXnDAzs2S5kTIzs2S5kTIzs2S5kTIzs2S5kTIzs2S5kTIzs2S5kTIzs2S5kTIzs2T9/2dtN//K//W/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x648 with 8 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# For plotting see http://xarray.pydata.org/en/stable/plotting.html\n",
    "\n",
    "from matplotlib.colors import LinearSegmentedColormap\n",
    "\n",
    "\n",
    "x = np.arange(0, np.pi, 0.1)       # 0., 0.1, 0.2, ..., 3.1\n",
    "y = np.arange(0, 2*np.pi, 0.1)     # 0., 0.1, 0.2, .............., 6.2\n",
    "X, Y = np.meshgrid(x, y)           # X and Y as a tuple form a numpy meshgrid\n",
    "Z = np.cos(X) * np.sin(Y) * 10     # Z is defined as the product of cosine and sine over this grid; so Z is a scalar field\n",
    "                                   #   It is also multiplied by 10 so it will have a range of -10 to 10 a prior\n",
    "    \n",
    "colors = [(1, 0, 0), (0, 1, 0), (0, 0, 1)]  # R -> G -> B\n",
    "                                            # A list of three tuples defining points in color space\n",
    "    \n",
    "n_bins = [3, 6, 10, 100]  # Discretizes the interpolation into bins\n",
    "                          # A list of 4 elements that will correspond to the 4 plots that are pending below\n",
    "                          #   indicating discretization. So the handwriting on the wall is that LinearSegmentColormap\n",
    "                          #   will be doing some interpolation on our behalf\n",
    "\n",
    "cmap_name = 'my_list'\n",
    "\n",
    "fig, axs = plt.subplots(2, 2, figsize=(6, 9))                                    # establishes 2 x 2 matrix of sub-plots\n",
    "fig.subplots_adjust(left=0.02, bottom=0.06, right=0.95, top=0.94, wspace=0.05)   #   with some interstices\n",
    "\n",
    "for n_bin, ax in zip(n_bins, axs.ravel()):             # note the use of both zip and .ravel() methods...\n",
    "                                                       #   clearly this iterates over the sub-plots with assignment\n",
    "                                                       #   of individual discretization values and axes\n",
    "                                                       #\n",
    "                                                       # The zip() function returns a zip object, which is an \n",
    "                                                       # iterator of tuples where the first item in each passed \n",
    "                                                       # iterator is paired together, and then the second item \n",
    "                                                       # in each passed iterator are paired together etc.\n",
    "                                                       #\n",
    "                                                       # numpy.ravel(array, order = ‘C’) : returns contiguous \n",
    "                                                       # flattened array(1D array with all the input-array elements \n",
    "                                                       # and with the same type as it). A copy is made only if needed.\n",
    "    \n",
    "    # Create the colormap\n",
    "    cm = LinearSegmentedColormap.from_list(cmap_name, colors, N=n_bin)   # cmap_name is a string attribute (not functional)\n",
    "                                                                         # colors is fixed: Our 3 stakes in the ground\n",
    "                                                                         #   N = 3 is no linear interpolation\n",
    "                                                                         #   N = 6 suggests trying 5... of interest\n",
    "                                                                         #   N = 10 and 100 follow intuitively\n",
    "    \n",
    "    # Fewer bins will result in \"coarser\" colomap interpolation\n",
    "    im = ax.imshow(Z, interpolation='nearest', origin='lower', cmap=cm)     # ax.imshow displays the Z scalar field\n",
    "    ax.set_title(\"N bins: %s\" % n_bin)                                      #   with a nice title\n",
    "    fig.colorbar(im, ax=ax)                                                 #   and the helpful color bars\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## S6. remarks on pyplot and DataArray native plotting\n",
    "\n",
    "### with print(object.__doc__) to show a helpful docstring\n",
    "\n",
    "Types and features\n",
    "\n",
    "* raster images\n",
    "* scatter plots\n",
    "* line plots\n",
    "* histograms\n",
    "* velocity vector fields (supported on ipyleaflet as well)\n",
    "* titles, annotations, overlays\n",
    "\n",
    "\n",
    "The idea is to take advantage of plotting capabilities built into xarray DataArrays and\n",
    "to integrate that with the matplotlib subplots that give us multiple plots for comparative \n",
    "viewing. For example a GOLIVE NetCDF file yields a vv_masked DataArray: Glacier speed \n",
    "across a UTM grid. This object has a .plot() method and even sub-plot methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['contour', 'contourf', 'hist', 'imshow', 'line', 'pcolormesh']"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# dirobj(m[0].vv_masked) produces a long list of methods including 'plot'.\n",
    "# To examine .plot: \n",
    "g.dirnou(m[0].vv_masked.plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "for .plot()...\n",
      "\n",
      "    Enables use of xarray.plot functions as attributes on a DataArray.\n",
      "    For example, DataArray.plot.imshow\n",
      "    \n",
      "\n",
      "\n",
      "for .plot.contour()...\n",
      "\n",
      "    Contour plot of 2d DataArray\n",
      "\n",
      "    Wraps :func:`matplotlib:matplotlib.pyplot.contour`\n",
      "    \n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    darray : DataArray\n",
      "        Must be 2 dimensional, unless creating faceted plots\n",
      "    x : string, optional\n",
      "        Coordinate for x axis. If None use darray.dims[1]\n",
      "    y : string, optional\n",
      "        Coordinate for y axis. If None use darray.dims[0]\n",
      "    figsize : tuple, optional\n",
      "        A tuple (width, height) of the figure in inches.\n",
      "        Mutually exclusive with ``size`` and ``ax``.\n",
      "    aspect : scalar, optional\n",
      "        Aspect ratio of plot, so that ``aspect * size`` gives the width in\n",
      "        inches. Only used if a ``size`` is provided.\n",
      "    size : scalar, optional\n",
      "        If provided, create a new figure for the plot with the given size.\n",
      "        Height (in inches) of each plot. See also: ``aspect``.\n",
      "    ax : matplotlib axes object, optional\n",
      "        Axis on which to plot this figure. By default, use the current axis.\n",
      "        Mutually exclusive with ``size`` and ``figsize``.\n",
      "    row : string, optional\n",
      "        If passed, make row faceted plots on this dimension name\n",
      "    col : string, optional\n",
      "        If passed, make column faceted plots on this dimension name\n",
      "    col_wrap : integer, optional\n",
      "        Use together with ``col`` to wrap faceted plots\n",
      "    xincrease : None, True, or False, optional\n",
      "        Should the values on the x axes be increasing from left to right?\n",
      "        if None, use the default for the matplotlib function.\n",
      "    yincrease : None, True, or False, optional\n",
      "        Should the values on the y axes be increasing from top to bottom?\n",
      "        if None, use the default for the matplotlib function.\n",
      "    add_colorbar : Boolean, optional\n",
      "        Adds colorbar to axis\n",
      "    add_labels : Boolean, optional\n",
      "        Use xarray metadata to label axes\n",
      "    vmin, vmax : floats, optional\n",
      "        Values to anchor the colormap, otherwise they are inferred from the\n",
      "        data and other keyword arguments. When a diverging dataset is inferred,\n",
      "        setting one of these values will fix the other by symmetry around\n",
      "        ``center``. Setting both values prevents use of a diverging colormap.\n",
      "        If discrete levels are provided as an explicit list, both of these\n",
      "        values are ignored.\n",
      "    cmap : matplotlib colormap name or object, optional\n",
      "        The mapping from data values to color space. If not provided, this\n",
      "        will be either be ``viridis`` (if the function infers a sequential\n",
      "        dataset) or ``RdBu_r`` (if the function infers a diverging dataset).\n",
      "        When `Seaborn` is installed, ``cmap`` may also be a `seaborn`\n",
      "        color palette. If ``cmap`` is seaborn color palette and the plot type\n",
      "        is not ``contour`` or ``contourf``, ``levels`` must also be specified.\n",
      "    colors : discrete colors to plot, optional\n",
      "        A single color or a list of colors. If the plot type is not ``contour``\n",
      "        or ``contourf``, the ``levels`` argument is required.\n",
      "    center : float, optional\n",
      "        The value at which to center the colormap. Passing this value implies\n",
      "        use of a diverging colormap. Setting it to ``False`` prevents use of a\n",
      "        diverging colormap.\n",
      "    robust : bool, optional\n",
      "        If True and ``vmin`` or ``vmax`` are absent, the colormap range is\n",
      "        computed with 2nd and 98th percentiles instead of the extreme values.\n",
      "    extend : {'neither', 'both', 'min', 'max'}, optional\n",
      "        How to draw arrows extending the colorbar beyond its limits. If not\n",
      "        provided, extend is inferred from vmin, vmax and the data limits.\n",
      "    levels : int or list-like object, optional\n",
      "        Split the colormap (cmap) into discrete color intervals. If an integer\n",
      "        is provided, \"nice\" levels are chosen based on the data range: this can\n",
      "        imply that the final number of levels is not exactly the expected one.\n",
      "        Setting ``vmin`` and/or ``vmax`` with ``levels=N`` is equivalent to\n",
      "        setting ``levels=np.linspace(vmin, vmax, N)``.\n",
      "    infer_intervals : bool, optional\n",
      "        Only applies to pcolormesh. If True, the coordinate intervals are\n",
      "        passed to pcolormesh. If False, the original coordinates are used\n",
      "        (this can be useful for certain map projections). The default is to\n",
      "        always infer intervals, unless the mesh is irregular and plotted on\n",
      "        a map projection.\n",
      "    subplot_kws : dict, optional\n",
      "        Dictionary of keyword arguments for matplotlib subplots. Only applies\n",
      "        to FacetGrid plotting.\n",
      "    cbar_ax : matplotlib Axes, optional\n",
      "        Axes in which to draw the colorbar.\n",
      "    cbar_kwargs : dict, optional\n",
      "        Dictionary of keyword arguments to pass to the colorbar.\n",
      "    **kwargs : optional\n",
      "        Additional arguments to wrapped matplotlib function\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    artist :\n",
      "        The same type of primitive artist that the wrapped matplotlib\n",
      "        function returns\n",
      "    \n",
      "\n",
      "\n",
      "for .plot.contourf()...\n",
      "\n",
      "    Filled contour plot of 2d DataArray\n",
      "\n",
      "    Wraps :func:`matplotlib:matplotlib.pyplot.contourf`\n",
      "    \n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    darray : DataArray\n",
      "        Must be 2 dimensional, unless creating faceted plots\n",
      "    x : string, optional\n",
      "        Coordinate for x axis. If None use darray.dims[1]\n",
      "    y : string, optional\n",
      "        Coordinate for y axis. If None use darray.dims[0]\n",
      "    figsize : tuple, optional\n",
      "        A tuple (width, height) of the figure in inches.\n",
      "        Mutually exclusive with ``size`` and ``ax``.\n",
      "    aspect : scalar, optional\n",
      "        Aspect ratio of plot, so that ``aspect * size`` gives the width in\n",
      "        inches. Only used if a ``size`` is provided.\n",
      "    size : scalar, optional\n",
      "        If provided, create a new figure for the plot with the given size.\n",
      "        Height (in inches) of each plot. See also: ``aspect``.\n",
      "    ax : matplotlib axes object, optional\n",
      "        Axis on which to plot this figure. By default, use the current axis.\n",
      "        Mutually exclusive with ``size`` and ``figsize``.\n",
      "    row : string, optional\n",
      "        If passed, make row faceted plots on this dimension name\n",
      "    col : string, optional\n",
      "        If passed, make column faceted plots on this dimension name\n",
      "    col_wrap : integer, optional\n",
      "        Use together with ``col`` to wrap faceted plots\n",
      "    xincrease : None, True, or False, optional\n",
      "        Should the values on the x axes be increasing from left to right?\n",
      "        if None, use the default for the matplotlib function.\n",
      "    yincrease : None, True, or False, optional\n",
      "        Should the values on the y axes be increasing from top to bottom?\n",
      "        if None, use the default for the matplotlib function.\n",
      "    add_colorbar : Boolean, optional\n",
      "        Adds colorbar to axis\n",
      "    add_labels : Boolean, optional\n",
      "        Use xarray metadata to label axes\n",
      "    vmin, vmax : floats, optional\n",
      "        Values to anchor the colormap, otherwise they are inferred from the\n",
      "        data and other keyword arguments. When a diverging dataset is inferred,\n",
      "        setting one of these values will fix the other by symmetry around\n",
      "        ``center``. Setting both values prevents use of a diverging colormap.\n",
      "        If discrete levels are provided as an explicit list, both of these\n",
      "        values are ignored.\n",
      "    cmap : matplotlib colormap name or object, optional\n",
      "        The mapping from data values to color space. If not provided, this\n",
      "        will be either be ``viridis`` (if the function infers a sequential\n",
      "        dataset) or ``RdBu_r`` (if the function infers a diverging dataset).\n",
      "        When `Seaborn` is installed, ``cmap`` may also be a `seaborn`\n",
      "        color palette. If ``cmap`` is seaborn color palette and the plot type\n",
      "        is not ``contour`` or ``contourf``, ``levels`` must also be specified.\n",
      "    colors : discrete colors to plot, optional\n",
      "        A single color or a list of colors. If the plot type is not ``contour``\n",
      "        or ``contourf``, the ``levels`` argument is required.\n",
      "    center : float, optional\n",
      "        The value at which to center the colormap. Passing this value implies\n",
      "        use of a diverging colormap. Setting it to ``False`` prevents use of a\n",
      "        diverging colormap.\n",
      "    robust : bool, optional\n",
      "        If True and ``vmin`` or ``vmax`` are absent, the colormap range is\n",
      "        computed with 2nd and 98th percentiles instead of the extreme values.\n",
      "    extend : {'neither', 'both', 'min', 'max'}, optional\n",
      "        How to draw arrows extending the colorbar beyond its limits. If not\n",
      "        provided, extend is inferred from vmin, vmax and the data limits.\n",
      "    levels : int or list-like object, optional\n",
      "        Split the colormap (cmap) into discrete color intervals. If an integer\n",
      "        is provided, \"nice\" levels are chosen based on the data range: this can\n",
      "        imply that the final number of levels is not exactly the expected one.\n",
      "        Setting ``vmin`` and/or ``vmax`` with ``levels=N`` is equivalent to\n",
      "        setting ``levels=np.linspace(vmin, vmax, N)``.\n",
      "    infer_intervals : bool, optional\n",
      "        Only applies to pcolormesh. If True, the coordinate intervals are\n",
      "        passed to pcolormesh. If False, the original coordinates are used\n",
      "        (this can be useful for certain map projections). The default is to\n",
      "        always infer intervals, unless the mesh is irregular and plotted on\n",
      "        a map projection.\n",
      "    subplot_kws : dict, optional\n",
      "        Dictionary of keyword arguments for matplotlib subplots. Only applies\n",
      "        to FacetGrid plotting.\n",
      "    cbar_ax : matplotlib Axes, optional\n",
      "        Axes in which to draw the colorbar.\n",
      "    cbar_kwargs : dict, optional\n",
      "        Dictionary of keyword arguments to pass to the colorbar.\n",
      "    **kwargs : optional\n",
      "        Additional arguments to wrapped matplotlib function\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    artist :\n",
      "        The same type of primitive artist that the wrapped matplotlib\n",
      "        function returns\n",
      "    \n",
      "\n",
      "\n",
      "for .plot.hist()...\n",
      "\n",
      "    Histogram of DataArray\n",
      "\n",
      "    Wraps :func:`matplotlib:matplotlib.pyplot.hist`\n",
      "\n",
      "    Plots N dimensional arrays by first flattening the array.\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    darray : DataArray\n",
      "        Can be any dimension\n",
      "    figsize : tuple, optional\n",
      "        A tuple (width, height) of the figure in inches.\n",
      "        Mutually exclusive with ``size`` and ``ax``.\n",
      "    aspect : scalar, optional\n",
      "        Aspect ratio of plot, so that ``aspect * size`` gives the width in\n",
      "        inches. Only used if a ``size`` is provided.\n",
      "    size : scalar, optional\n",
      "        If provided, create a new figure for the plot with the given size.\n",
      "        Height (in inches) of each plot. See also: ``aspect``.\n",
      "    ax : matplotlib axes object, optional\n",
      "        Axis on which to plot this figure. By default, use the current axis.\n",
      "        Mutually exclusive with ``size`` and ``figsize``.\n",
      "    **kwargs : optional\n",
      "        Additional keyword arguments to matplotlib.pyplot.hist\n",
      "\n",
      "    \n",
      "\n",
      "\n",
      "for .plot.imshow()...\n",
      "\n",
      "    Image plot of 2d DataArray using matplotlib.pyplot\n",
      "\n",
      "    Wraps :func:`matplotlib:matplotlib.pyplot.imshow`\n",
      "\n",
      "    While other plot methods require the DataArray to be strictly\n",
      "    two-dimensional, ``imshow`` also accepts a 3D array where some\n",
      "    dimension can be interpreted as RGB or RGBA color channels and\n",
      "    allows this dimension to be specified via the kwarg ``rgb=``.\n",
      "\n",
      "    Unlike matplotlib, Xarray can apply ``vmin`` and ``vmax`` to RGB or RGBA\n",
      "    data, by applying a single scaling factor and offset to all bands.\n",
      "    Passing  ``robust=True`` infers ``vmin`` and ``vmax``\n",
      "    :ref:`in the usual way <robust-plotting>`.\n",
      "\n",
      "    .. note::\n",
      "        This function needs uniformly spaced coordinates to\n",
      "        properly label the axes. Call DataArray.plot() to check.\n",
      "\n",
      "    The pixels are centered on the coordinates values. Ie, if the coordinate\n",
      "    value is 3.2 then the pixels for those coordinates will be centered on 3.2.\n",
      "    \n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    darray : DataArray\n",
      "        Must be 2 dimensional, unless creating faceted plots\n",
      "    x : string, optional\n",
      "        Coordinate for x axis. If None use darray.dims[1]\n",
      "    y : string, optional\n",
      "        Coordinate for y axis. If None use darray.dims[0]\n",
      "    figsize : tuple, optional\n",
      "        A tuple (width, height) of the figure in inches.\n",
      "        Mutually exclusive with ``size`` and ``ax``.\n",
      "    aspect : scalar, optional\n",
      "        Aspect ratio of plot, so that ``aspect * size`` gives the width in\n",
      "        inches. Only used if a ``size`` is provided.\n",
      "    size : scalar, optional\n",
      "        If provided, create a new figure for the plot with the given size.\n",
      "        Height (in inches) of each plot. See also: ``aspect``.\n",
      "    ax : matplotlib axes object, optional\n",
      "        Axis on which to plot this figure. By default, use the current axis.\n",
      "        Mutually exclusive with ``size`` and ``figsize``.\n",
      "    row : string, optional\n",
      "        If passed, make row faceted plots on this dimension name\n",
      "    col : string, optional\n",
      "        If passed, make column faceted plots on this dimension name\n",
      "    col_wrap : integer, optional\n",
      "        Use together with ``col`` to wrap faceted plots\n",
      "    xincrease : None, True, or False, optional\n",
      "        Should the values on the x axes be increasing from left to right?\n",
      "        if None, use the default for the matplotlib function.\n",
      "    yincrease : None, True, or False, optional\n",
      "        Should the values on the y axes be increasing from top to bottom?\n",
      "        if None, use the default for the matplotlib function.\n",
      "    add_colorbar : Boolean, optional\n",
      "        Adds colorbar to axis\n",
      "    add_labels : Boolean, optional\n",
      "        Use xarray metadata to label axes\n",
      "    vmin, vmax : floats, optional\n",
      "        Values to anchor the colormap, otherwise they are inferred from the\n",
      "        data and other keyword arguments. When a diverging dataset is inferred,\n",
      "        setting one of these values will fix the other by symmetry around\n",
      "        ``center``. Setting both values prevents use of a diverging colormap.\n",
      "        If discrete levels are provided as an explicit list, both of these\n",
      "        values are ignored.\n",
      "    cmap : matplotlib colormap name or object, optional\n",
      "        The mapping from data values to color space. If not provided, this\n",
      "        will be either be ``viridis`` (if the function infers a sequential\n",
      "        dataset) or ``RdBu_r`` (if the function infers a diverging dataset).\n",
      "        When `Seaborn` is installed, ``cmap`` may also be a `seaborn`\n",
      "        color palette. If ``cmap`` is seaborn color palette and the plot type\n",
      "        is not ``contour`` or ``contourf``, ``levels`` must also be specified.\n",
      "    colors : discrete colors to plot, optional\n",
      "        A single color or a list of colors. If the plot type is not ``contour``\n",
      "        or ``contourf``, the ``levels`` argument is required.\n",
      "    center : float, optional\n",
      "        The value at which to center the colormap. Passing this value implies\n",
      "        use of a diverging colormap. Setting it to ``False`` prevents use of a\n",
      "        diverging colormap.\n",
      "    robust : bool, optional\n",
      "        If True and ``vmin`` or ``vmax`` are absent, the colormap range is\n",
      "        computed with 2nd and 98th percentiles instead of the extreme values.\n",
      "    extend : {'neither', 'both', 'min', 'max'}, optional\n",
      "        How to draw arrows extending the colorbar beyond its limits. If not\n",
      "        provided, extend is inferred from vmin, vmax and the data limits.\n",
      "    levels : int or list-like object, optional\n",
      "        Split the colormap (cmap) into discrete color intervals. If an integer\n",
      "        is provided, \"nice\" levels are chosen based on the data range: this can\n",
      "        imply that the final number of levels is not exactly the expected one.\n",
      "        Setting ``vmin`` and/or ``vmax`` with ``levels=N`` is equivalent to\n",
      "        setting ``levels=np.linspace(vmin, vmax, N)``.\n",
      "    infer_intervals : bool, optional\n",
      "        Only applies to pcolormesh. If True, the coordinate intervals are\n",
      "        passed to pcolormesh. If False, the original coordinates are used\n",
      "        (this can be useful for certain map projections). The default is to\n",
      "        always infer intervals, unless the mesh is irregular and plotted on\n",
      "        a map projection.\n",
      "    subplot_kws : dict, optional\n",
      "        Dictionary of keyword arguments for matplotlib subplots. Only applies\n",
      "        to FacetGrid plotting.\n",
      "    cbar_ax : matplotlib Axes, optional\n",
      "        Axes in which to draw the colorbar.\n",
      "    cbar_kwargs : dict, optional\n",
      "        Dictionary of keyword arguments to pass to the colorbar.\n",
      "    **kwargs : optional\n",
      "        Additional arguments to wrapped matplotlib function\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    artist :\n",
      "        The same type of primitive artist that the wrapped matplotlib\n",
      "        function returns\n",
      "    \n",
      "\n",
      "\n",
      "for .plot.line()...\n",
      "\n",
      "    Line plot of DataArray index against values\n",
      "\n",
      "    Wraps :func:`matplotlib:matplotlib.pyplot.plot`\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    darray : DataArray\n",
      "        Must be 1 dimensional\n",
      "    figsize : tuple, optional\n",
      "        A tuple (width, height) of the figure in inches.\n",
      "        Mutually exclusive with ``size`` and ``ax``.\n",
      "    aspect : scalar, optional\n",
      "        Aspect ratio of plot, so that ``aspect * size`` gives the width in\n",
      "        inches. Only used if a ``size`` is provided.\n",
      "    size : scalar, optional\n",
      "        If provided, create a new figure for the plot with the given size.\n",
      "        Height (in inches) of each plot. See also: ``aspect``.\n",
      "    ax : matplotlib axes object, optional\n",
      "        Axis on which to plot this figure. By default, use the current axis.\n",
      "        Mutually exclusive with ``size`` and ``figsize``.\n",
      "    hue : string, optional\n",
      "        Coordinate for which you want multiple lines plotted.\n",
      "    x, y : string, optional\n",
      "        Coordinates for x, y axis. Only one of these may be specified.\n",
      "        The other coordinate plots values from the DataArray on which this\n",
      "        plot method is called.\n",
      "    xincrease : None, True, or False, optional\n",
      "        Should the values on the x axes be increasing from left to right?\n",
      "        if None, use the default for the matplotlib function.\n",
      "    yincrease : None, True, or False, optional\n",
      "        Should the values on the y axes be increasing from top to bottom?\n",
      "        if None, use the default for the matplotlib function.\n",
      "    add_legend : boolean, optional\n",
      "        Add legend with y axis coordinates (2D inputs only).\n",
      "    *args, **kwargs : optional\n",
      "        Additional arguments to matplotlib.pyplot.plot\n",
      "\n",
      "    \n",
      "\n",
      "\n",
      "for .plot.pcolormesh()...\n",
      "\n",
      "    Pseudocolor plot of 2d DataArray\n",
      "\n",
      "    Wraps :func:`matplotlib:matplotlib.pyplot.pcolormesh`\n",
      "    \n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    darray : DataArray\n",
      "        Must be 2 dimensional, unless creating faceted plots\n",
      "    x : string, optional\n",
      "        Coordinate for x axis. If None use darray.dims[1]\n",
      "    y : string, optional\n",
      "        Coordinate for y axis. If None use darray.dims[0]\n",
      "    figsize : tuple, optional\n",
      "        A tuple (width, height) of the figure in inches.\n",
      "        Mutually exclusive with ``size`` and ``ax``.\n",
      "    aspect : scalar, optional\n",
      "        Aspect ratio of plot, so that ``aspect * size`` gives the width in\n",
      "        inches. Only used if a ``size`` is provided.\n",
      "    size : scalar, optional\n",
      "        If provided, create a new figure for the plot with the given size.\n",
      "        Height (in inches) of each plot. See also: ``aspect``.\n",
      "    ax : matplotlib axes object, optional\n",
      "        Axis on which to plot this figure. By default, use the current axis.\n",
      "        Mutually exclusive with ``size`` and ``figsize``.\n",
      "    row : string, optional\n",
      "        If passed, make row faceted plots on this dimension name\n",
      "    col : string, optional\n",
      "        If passed, make column faceted plots on this dimension name\n",
      "    col_wrap : integer, optional\n",
      "        Use together with ``col`` to wrap faceted plots\n",
      "    xincrease : None, True, or False, optional\n",
      "        Should the values on the x axes be increasing from left to right?\n",
      "        if None, use the default for the matplotlib function.\n",
      "    yincrease : None, True, or False, optional\n",
      "        Should the values on the y axes be increasing from top to bottom?\n",
      "        if None, use the default for the matplotlib function.\n",
      "    add_colorbar : Boolean, optional\n",
      "        Adds colorbar to axis\n",
      "    add_labels : Boolean, optional\n",
      "        Use xarray metadata to label axes\n",
      "    vmin, vmax : floats, optional\n",
      "        Values to anchor the colormap, otherwise they are inferred from the\n",
      "        data and other keyword arguments. When a diverging dataset is inferred,\n",
      "        setting one of these values will fix the other by symmetry around\n",
      "        ``center``. Setting both values prevents use of a diverging colormap.\n",
      "        If discrete levels are provided as an explicit list, both of these\n",
      "        values are ignored.\n",
      "    cmap : matplotlib colormap name or object, optional\n",
      "        The mapping from data values to color space. If not provided, this\n",
      "        will be either be ``viridis`` (if the function infers a sequential\n",
      "        dataset) or ``RdBu_r`` (if the function infers a diverging dataset).\n",
      "        When `Seaborn` is installed, ``cmap`` may also be a `seaborn`\n",
      "        color palette. If ``cmap`` is seaborn color palette and the plot type\n",
      "        is not ``contour`` or ``contourf``, ``levels`` must also be specified.\n",
      "    colors : discrete colors to plot, optional\n",
      "        A single color or a list of colors. If the plot type is not ``contour``\n",
      "        or ``contourf``, the ``levels`` argument is required.\n",
      "    center : float, optional\n",
      "        The value at which to center the colormap. Passing this value implies\n",
      "        use of a diverging colormap. Setting it to ``False`` prevents use of a\n",
      "        diverging colormap.\n",
      "    robust : bool, optional\n",
      "        If True and ``vmin`` or ``vmax`` are absent, the colormap range is\n",
      "        computed with 2nd and 98th percentiles instead of the extreme values.\n",
      "    extend : {'neither', 'both', 'min', 'max'}, optional\n",
      "        How to draw arrows extending the colorbar beyond its limits. If not\n",
      "        provided, extend is inferred from vmin, vmax and the data limits.\n",
      "    levels : int or list-like object, optional\n",
      "        Split the colormap (cmap) into discrete color intervals. If an integer\n",
      "        is provided, \"nice\" levels are chosen based on the data range: this can\n",
      "        imply that the final number of levels is not exactly the expected one.\n",
      "        Setting ``vmin`` and/or ``vmax`` with ``levels=N`` is equivalent to\n",
      "        setting ``levels=np.linspace(vmin, vmax, N)``.\n",
      "    infer_intervals : bool, optional\n",
      "        Only applies to pcolormesh. If True, the coordinate intervals are\n",
      "        passed to pcolormesh. If False, the original coordinates are used\n",
      "        (this can be useful for certain map projections). The default is to\n",
      "        always infer intervals, unless the mesh is irregular and plotted on\n",
      "        a map projection.\n",
      "    subplot_kws : dict, optional\n",
      "        Dictionary of keyword arguments for matplotlib subplots. Only applies\n",
      "        to FacetGrid plotting.\n",
      "    cbar_ax : matplotlib Axes, optional\n",
      "        Axes in which to draw the colorbar.\n",
      "    cbar_kwargs : dict, optional\n",
      "        Dictionary of keyword arguments to pass to the colorbar.\n",
      "    **kwargs : optional\n",
      "        Additional arguments to wrapped matplotlib function\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    artist :\n",
      "        The same type of primitive artist that the wrapped matplotlib\n",
      "        function returns\n",
      "    \n"
     ]
    }
   ],
   "source": [
    "print('for .plot()...')\n",
    "print(m[0].vv_masked.plot.__doc__)\n",
    "print('\\n\\nfor .plot.contour()...')\n",
    "print(m[0].vv_masked.plot.contour.__doc__)\n",
    "print('\\n\\nfor .plot.contourf()...')\n",
    "print(m[0].vv_masked.plot.contourf.__doc__)\n",
    "print('\\n\\nfor .plot.hist()...')\n",
    "print(m[0].vv_masked.plot.hist.__doc__)\n",
    "print('\\n\\nfor .plot.imshow()...')\n",
    "print(m[0].vv_masked.plot.imshow.__doc__)\n",
    "print('\\n\\nfor .plot.line()...')\n",
    "print(m[0].vv_masked.plot.line.__doc__)\n",
    "print('\\n\\nfor .plot.pcolormesh()...')\n",
    "print(m[0].vv_masked.plot.pcolormesh.__doc__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## S7. Showing vector fields\n",
    "\n",
    "\n",
    "It would be pleasant to display localized vector fields using pyplot's plt.quiver() method. This is not yet \n",
    "working here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Does not work... also interesting thought: Does a .velocity() method exist?\n",
    "# X, Y = np.meshgrid(m[2]['y'], m[2]['x'])\n",
    "# U = m[2]['vy_masked']\n",
    "# V = m[2]['vx_masked']\n",
    "# \n",
    "# a = plt.quiver(Y, X, V, U)         # what does 'units' keyword do here if added?\n",
    "# plt.xlim(535000,540000)\n",
    "# plt.ylim(6675000,6680000)\n",
    "# plt.title('ice velocities')\n",
    "# print('velocity vector field')\n",
    "# ax = a.axes\n",
    "# plt.xlim(535000,540000)\n",
    "# plt.ylim(6675000,6680000)\n",
    "# , figsize=(20,20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## David Shean's golive_proc.py program"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This code is golive_proc.py from https://github.com/dshean/vmap/tree/master/vmap\n",
    "# It includes this line at the top: #! /usr/bin/env python\n",
    "\n",
    "if False:    # So the cell runs without doing anything\n",
    "\n",
    "    # Next section (9 lines) commented out as in the source: Just indicates preparation steps\n",
    "    \n",
    "    # Get the time-series dataset for path 46 row 27 of Landsat 8\n",
    "    # wget -r ftp://sidads.colorado.edu/pub/DATASETS/nsidc0710_landsat8_golive_ice_velocity_v1/p046_r027/*nc\n",
    "    # parallel 'gdal_translate NETCDF:\"{}\":vv_masked {.}_vv_masked.tif' ::: *nc\n",
    "    # parallel './golive_rename.py {}' ::: L8*vv.tif\n",
    "    # parallel './golive_rename.py {}' ::: L8*vv_masked.tif\n",
    "\n",
    "    # These lines are running some code while providing a UTM bounding box for Mt. Rainier\n",
    "    # make_stack.py -te '577880 5175800 611800 5205050' 2*vv.tif\n",
    "    # make_stack.py -te '577880 5175800 611800 5205050' 2*vv_masked.tif\n",
    "    # dem_mosaic --median *vv_masked.tif -o vv_masked\n",
    "\n",
    "    in_fn = 'L8_046_027_048_2015_350_2016_033_v1_vv_masked.tif'        # golive result file path 27 row 48 (Rainier)\n",
    "                                                                       #   this will be a 48-day interval\n",
    "    a = in_fn.split('_')\n",
    "    dt1 = timelib.doy2dt(int(a[4]), int(a[5]))                         # Get dt1, dt2, dtc for center; interval in days\n",
    "    dt2 = timelib.doy2dt(int(a[6]), int(a[7]))\n",
    "    dtc = timelib.center_date(dt1, dt2)\n",
    "    ndays = (dt2 - dt1).days\n",
    "    \n",
    "    # Build up an output filename\n",
    "    out_fn = '%s_%iday_%s-%s_%s' % (dtc.strftime('%Y%m%d'), ndays, dt1.strftime('%Y%m%d'), dt2.strftime('%Y%m%d'), in_fn) \n",
    "    print(out_fn)\n",
    "    \n",
    "    # clobber any pre-existing output file \n",
    "    if os.path.exists(out_fn): os.remove(out_fn)\n",
    "        \n",
    "    # set up symlink from out_fn to in_fn\n",
    "    os.symlink(in_fn, out_fn)\n",
    "\n",
    "    # b_ds will now be a GDAL dataset\n",
    "    b_ds = gdal.Open(out_fn)\n",
    "    \n",
    "    # get 'masked array' from b_ds\n",
    "    b = iolib.ds_getma(b_ds)\n",
    "\n",
    "    # threshold for Rainier speed filter\n",
    "    b_thresh = 1.5            # most likely 1.5 meters per day to ignore slow ice noise\n",
    "    \n",
    "    b[b > b_thresh] = np.ma.masked        # does the masking based on speed\n",
    "    \n",
    "    out_fn_masked = os.path.splitext(out_fn)[0]+'_filt.tif'\n",
    "    iolib.writeGTiff(b,out_fn_masked,b_ds)                       # write an additional output file reflecting the speed filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### lgo_mask\n",
    "\n",
    "- 0 is moving ice\n",
    "- 1 is non-moving land\n",
    "- 2 is ocean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'timeit' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<timed exec>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'timeit' is not defined"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# prints time used when the cell is done running\n",
    "# use the timer to ingest this time into a program...\n",
    "\n",
    "start_time = timeit.default_timer()\n",
    "\n",
    "# Here we average time series mask values (0 for ice, 1 for land) across time-series using values.mean()\n",
    "if True:\n",
    "    for j in range(401,403):\n",
    "        for i in range(701,703):\n",
    "            print(lgo.isel(y=j,x=i).data.mean())\n",
    "\n",
    "# ...as here the time interval comes programmatically off the timer\n",
    "print('\\n', timeit.default_timer() - start_time, '\\n\\ncompare with...\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "if True:\n",
    "    # Operating on a single pixel...\n",
    "    #   choosing the mid-point of the lower transect, closest to Malaspina\n",
    "    xMid = utm0[0][0] + (number_of_postings[0]/2.)*relative_postings[0][0]    # These must start at utm0 and add a certain \n",
    "    yMid = utm0[0][1] + (number_of_postings[0]/2.)*relative_postings[0][1]    #   number of the relative_postings vectors\n",
    "    print('UTM of point of interest is', xMid, yMid, 'as easting - northing\\n')\n",
    "    s = vv.sel(x=xMid,y=yMid,method='nearest')           # selects a time series at one location (xMid, yMid)\n",
    "    ss = s.where(s.data >= 0.)                           # substitutes nan for -9999. no data value\n",
    "    # each of the following does something relevant: \n",
    "    #   print(s.values)                           # prints speeds as an ndarray (not a list; no commas!)\n",
    "    #   print(s)                                  # describes DataArray (since it was derived from m['vv_masked'])        \n",
    "    #   s.plot()                                  # uses .plot() native to DataArray to show a plot (-9999 to 30 m/d)\n",
    "    #   print(ss.values)                          # shows the nan value substitution in ss\n",
    "    #   print(ss.mean(skipna=True).values)        # use DataArray .mean() to ignore nans and print a mean of 4.127929 (sd 1.0109311)\n",
    "    #   print(np.std([1., 2., 3.]))               # sidebar numpy operating on a list\n",
    "    #   print(s.x, '\\n')\n",
    "    #   print(s, '\\n')\n",
    "    #   print(s.masking_info, '\\n')\n",
    "    #   print(s['time'])\n",
    "    #   dirobj(s)\n",
    "    #   dirobj(s.x)\n",
    "    #   print(ss.std(), ss.mean(), ss.std()/ss.mean())\n",
    "    #   print(ss.data)                                    # prints all 255 values\n",
    "    print(ss.mean().data)                                 # about 1                    \n",
    "    print(ss.std().data)                                  # about 4\n",
    "    print(ss.std().data/ss.mean().data) \n",
    "    # print(dirobj(ss))\n",
    "    print(ss.count('time').data)                          # 78 values (i.e. not nan)\n",
    "    print(ss.sum().data)                                  # sum of all 78 values"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}