{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Copyright 2013-2016 The Salish Sea NEMO Project and\n", "# The University of British Columbia\n", "\n", "# Licensed under the Apache License, Version 2.0 (the \"License\");\n", "# you may not use this file except in compliance with the License.\n", "# You may obtain a copy of the License at\n", "\n", "# http://www.apache.org/licenses/LICENSE-2.0\n", "\n", "# Unless required by applicable law or agreed to in writing, software\n", "# distributed under the License is distributed on an \"AS IS\" BASIS,\n", "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n", "# See the License for the specific language governing permissions and\n", "# limitations under the License." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "'Functions for working with geographical data and model results.\\n'" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"Functions for working with geographical data and model results.\n", "\"\"\"" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from salishsea_tools import tidetools\n", "import netCDF4 as nc\n", "import pdb\n", "import pandas as pd\n", "import xarray as xr\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def distance_along_curve(lons, lats):\n", " \"\"\"Calculate cumulative distance in km between points in lons, lats\n", "\n", " :arg lons: 1D array of longitude points.\n", " :type lons: :py:class:`numpy.ndarray`\n", "\n", " :arg lats: 1D array of latitude points.\n", " :type lats: :py:class:`numpy.ndarray`\n", "\n", " :returns: Cummulative point-by-point distance along track in km.\n", " :rtype: :py:class:`numpy.ndarray`\n", " \"\"\"\n", " dist = np.cumsum(haversine(lons[1:], lats[1:], lons[:-1], lats[:-1]))\n", " return np.insert(dist, 0, 0)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def haversine(lon1, lat1, lon2, lat2):\n", " \"\"\"Calculate the great-circle distance in kilometers between two points\n", " on a sphere from their longitudes and latitudes.\n", "\n", " Reference: http://www.movable-type.co.uk/scripts/latlong.html\n", "\n", " :arg lon1: Longitude of point 1.\n", " :type lon1: float or :py:class:`numpy.ndarray`\n", "\n", " :arg lat1: Latitude of point 1.\n", " :type lat1: float or :py:class:`numpy.ndarray`\n", "\n", " :arg lon2: Longitude of point 2.\n", " :type lon2: float or :py:class:`numpy.ndarray`\n", "\n", " :arg lat2: Latitude of point 2.\n", " :type lat2: float or :py:class:`numpy.ndarray`\n", "\n", " :returns: Great-circle distance between two points in km\n", " :rtype: float or :py:class:`numpy.ndarray`\n", " \"\"\"\n", " lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])\n", " dlon = lon2 - lon1\n", " dlat = lat2 - lat1\n", " a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2\n", " c = 2 * np.arcsin(np.sqrt(a))\n", " km = 6367 * c\n", " return km" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def _spiral_search_for_closest_water_point(\n", " j, i, land_mask, lon, lat, model_lons, model_lats\n", "):\n", " # Searches in a spiral pattern around grid element (j,i)\n", " # for the closest water point to the the coordinate (lat,lon)\n", "\n", " jmax, imax = land_mask.shape\n", " # Limit on size of grid search\n", " max_search_dist = int(model_lats.shape[1]/4)\n", " closest_point = None\n", " j_s, i_s = j, i # starting points is j, i\n", " dj, di = 0, -1\n", " # move j_s, i_s in a square spiral centred at j, i\n", " while (i_s-i) <= max_search_dist:\n", " if any([(j_s-j) == (i_s-i),\n", " ((j_s-j) < 0 and (j_s-j) == -(i_s-i)),\n", " ((j_s-j) > 0 and (j_s-j) == 1-(i_s-i))]):\n", " # Hit the corner of the spiral- change direction\n", " dj, di = -di, dj\n", " i_s, j_s = i_s+di, j_s+dj # Take a step to next square\n", " if all([i_s >= 0,\n", " i_s < imax, j_s >= 0,\n", " j_s < jmax,\n", " not land_mask[j_s, i_s]]\n", " ):\n", " # Found a water point, how close is it?\n", " actual_dist = haversine(\n", " lon, lat, model_lons[j_s, i_s], model_lats[j_s, i_s])\n", " if closest_point is None:\n", " min_dist = actual_dist\n", " closest_point = (j_s, i_s)\n", " elif actual_dist < min_dist:\n", " # Keep record of closest point\n", " min_dist = actual_dist\n", " closest_point = (j_s, i_s)\n", " # Assumes grids are square- reduces search radius to only\n", " # check grids that could potentially be closer than this\n", " grid_dist = int(((i_s-i)**2 + (j_s-j)**2)**0.5)\n", " if (grid_dist + 1) < max_search_dist:\n", " # Reduce stopping distance for spiral-\n", " # just need to check that no points closer than this one\n", " max_search_dist = grid_dist + 1\n", " if closest_point is not None:\n", " return closest_point\n", " else:\n", " raise ValueError('lat/lon on land and no nearby water point found')\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "def find_closest_model_point(\n", " lon, lat, model_lons, model_lats, grid='NEMO', land_mask=None,\n", " tols={\n", " 'NEMO': {'tol_lon': 0.0104, 'tol_lat': 0.00388},\n", " 'GEM2.5': {'tol_lon': 0.016, 'tol_lat': 0.011},\n", " }\n", "):\n", " \"\"\"Returns the grid coordinates of the closest model point\n", " to a specified lon/lat. If land_mask is provided, returns the closest\n", " water point.\n", "\n", " Example:\n", "\n", " .. code-block:: python\n", " j, i = find_closest_model_point(\n", " -125.5,49.2,model_lons,model_lats,land_mask=bathy.mask)\n", "\n", " where bathy, model_lons and model_lats are returned from\n", " :py:func:`salishsea_tools.tidetools.get_bathy_data`.\n", "\n", " j is the y-index(latitude), i is the x-index(longitude)\n", "\n", " :arg float lon: longitude to find closest grid point to\n", "\n", " :arg float lat: latitude to find closest grid point to\n", "\n", " :arg model_lons: specified model longitude grid\n", " :type model_lons: :py:obj:`numpy.ndarray`\n", "\n", " :arg model_lats: specified model latitude grid\n", " :type model_lats: :py:obj:`numpy.ndarray`\n", "\n", " :arg grid: specify which default lon/lat tolerances\n", " :type grid: string\n", "\n", " :arg land_mask: describes which grid coordinates are land\n", " :type land_mask: numpy array\n", "\n", " :arg tols: stored default tols for different grid types\n", " :type tols: dict\n", "\n", " :returns: yind, xind\n", " \"\"\"\n", "\n", " if grid not in tols:\n", " raise KeyError('The provided grid type is not in tols.\\\n", " Use another grid type or add your grid type to tols.')\n", "\n", " # Search for a grid point with longitude and latitude within\n", " # tolerance of measured location\n", " j_list, i_list = np.where(\n", " np.logical_and(\n", " (np.logical_and(model_lons > lon - tols[grid]['tol_lon'],\n", " model_lons < lon + tols[grid]['tol_lon'])),\n", " (np.logical_and(model_lats > lat - tols[grid]['tol_lat'],\n", " model_lats < lat + tols[grid]['tol_lat']))\n", " )\n", " )\n", "\n", " if len(j_list) == 0:\n", " raise ValueError(\n", " 'No model point found. tol_lon/tol_lat too small or '\n", " 'lon/lat outside of domain.')\n", " try:\n", " j, i = map(np.asscalar, (j_list, i_list))\n", " except ValueError:\n", " # Several points within tolerance\n", " # Calculate distances for all and choose the closest\n", " test_lons = [model_lons[j_list[n], i_list[n]] for n in range(len(j_list))]\n", " test_lats = [model_lats[j_list[n], i_list[n]] for n in range(len(j_list))]\n", " dists = haversine(\n", " np.array([lon] * i_list.size), np.array([lat] * j_list.size),\n", " test_lons, model_lats[j_list, i_list])\n", " n = dists.argmin()\n", " j, i = map(np.asscalar, (j_list[n], i_list[n]))\n", "\n", " # If point is on land and land mask is provided\n", " # try to find closest water point\n", " if land_mask is None or not land_mask[j, i]:\n", " return j, i\n", " try:\n", " return _spiral_search_for_closest_water_point(\n", " j, i, land_mask, lon, lat, model_lons, model_lats)\n", " except ValueError:\n", " raise ValueError(\n", " 'lat/lon on land and no nearby water point found')\n" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def find_model_point(lon, lat, X, Y, tol_lon=0.016, tol_lat=0.011):\n", " \"\"\"Finds a model grid point close to a specified latitude and longitude.\n", " Should be used for non-NEMO grids like the atmospheric forcing grid.\n", "\n", " :arg lon: The longitude we are trying to match.\n", " :type lon: float\n", "\n", " :arg lat: The latitude we are trying to match.\n", " :type lat: float\n", "\n", " :arg X: The model longitude grid.\n", " :type X: numpy array\n", "\n", " :arg Y: The model latitude grid.\n", " :type Y: numpy array\n", "\n", " :arg tol_lon: tolerance on grid spacing for longitude\n", " :type tol_lon: float\n", "\n", " :arg tol_lat: tolerance on grid spacing for latitude\n", " :type tol_lat: float\n", "\n", " :returns: j-index and i-index of the closest model grid point.\n", " \"\"\"\n", "\n", " # Search for a grid point with longitude or latitude within\n", " # tolerance of measured location\n", " j, i = np.where(\n", " np.logical_and(\n", " (np.logical_and(X > lon - tol_lon, X < lon + tol_lon)),\n", " (np.logical_and(Y > lat - tol_lat, Y < lat + tol_lat))))\n", "\n", " if j.size > 1 or i.size > 1:\n", " raise ValueError(\n", " 'Multiple model points found. tol_lon/tol_lat too big.'\n", " )\n", " elif not j or not i:\n", " raise ValueError(\n", " 'No model point found. tol_lon/tol_lat too small or '\n", " 'lon/lat outside of domain.'\n", " )\n", " return j, i" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def find_closest_model_point(lon, lat, X, Y, bathy, lon_tol=0.0052,\n", " lat_tol=0.00189, allow_land=False):\n", " \"\"\"Returns the grid co-ordinates of the closest non-land model point\n", " to a specified lon/lat.\n", "\n", " e.g. yind, xind = find_closest_model_point(-125.5,49.2,X,Y,bathy)\n", " where bathy, X and Y are returned from get_bathy_data(grid).\n", " yind is the y-index(latitude), xind is the x-index(longitude)\n", "\n", " :arg lon: specified longitude\n", " :type lon: float\n", "\n", " :arg lat: specified latitude\n", " :type lat: float\n", "\n", " :arg X: specified model longitude\n", " :type X: numpy array\n", "\n", " :arg Y: specified model latitude\n", " :type Y: numpy array\n", "\n", " :arg bathy: model bathymetry\n", " :type bathy: numpy array\n", "\n", " :arg lon_tol: tolerance value for seaching in longitude\n", " :type lon_tol: float\n", "\n", " :arg lat_tol: tolerance value for searching in latitude\n", " :type lat_tol: float\n", "\n", " :arg allow_land: whether code should return a land point or closest\n", " water point\n", " :type allow_land: boolean\n", "\n", " :returns: yind, xind\n", " \"\"\"\n", " # Tolerance for searching for grid points\n", " # (default is approx. distances between adjacent grid points)\n", " tol1 = lon_tol # lon\n", " tol2 = lat_tol # lat\n", "\n", " # Search for a grid point with lon/lat within tolerance of\n", " # measured location\n", " x1, y1 = np.where(\n", " np.logical_and(\n", " (np.logical_and(X > lon-tol1, X < lon+tol1)),\n", " (np.logical_and(Y > lat-tol2, Y < lat+tol2))))\n", "\n", " # check size of arrays so we don't go out of bounds in our search\n", " xmax, ymax = bathy.shape\n", " if np.size(x1) != 0:\n", " x1 = x1[0]\n", " y1 = y1[0]\n", " # What if more than one point is returned from this search?\n", " # Just take the first one...\n", " #\n", " # If x1,y1 is masked, search 3x3 grid around.\n", " # If all those points are masked, search 4x4 grid around, etc.\n", " for ii in np.arange(1, 100):\n", " if bathy.mask[x1, y1] and not allow_land:\n", " for i in np.arange(max(0, x1-ii), min(xmax, x1+ii+1)):\n", " for j in np.arange(max(0, y1-ii), min(ymax, y1+ii+1)):\n", " if not bathy.mask[i, j]:\n", " break\n", " if not bathy.mask[i, j]:\n", " break\n", " if not bathy.mask[i, j]:\n", " break\n", " else:\n", " i = x1\n", " j = y1\n", " else:\n", " i = []\n", " j = []\n", " return i, j" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "grid_B = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')\n", "bathy, X, Y = tidetools.get_bathy_data(grid_B)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(403, 46)\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA20AAAN5CAYAAACbtdNHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xm8XXV97//3Jwnz6HAVCArKWHBAaiFSez3WCXBAa1vB\noi1FRQWrYilotabq9aeoFb1WrQJWah2utFhr+Sl6NVyVCghSFBIGlYCBRBoCJGEMfO8fOddHGhIS\nsw9Z35PzfD4ePMhee52c93l4JHllrb1TrbUAAADQp2lDDwAAAGDtRBsAAEDHRBsAAEDHRBsAAEDH\nRBsAAEDHRBsAAEDHJl20VdVbqur+qnr4Wp5/c1X9pKour6p/rKrNx4+/s6p+UVWXjv9z6Hp8rv+/\nqpZU1Vcn+usAAABYH11GW1U9o6o+s4bjuyZ5TpL5a/m4XZK8IcmBrbUnJZmR5MhVTvmb1tqB4/98\nfT2mnJrk6F/7CwAAAJggXUbbuDX9rd8fTnLSOj5uepJtqmpGkq2T3LjKc7X6yVU1rapOraoLq+qy\nqnr1rwa09p0ky3796QAAABOj52j7L4FVVS9KckNr7cdr+4DW2o1JPpTk+iQLktzaWvvWKqecMB5m\np1fVDuPHjh0/7+AkByV5TVXtNpFfCAAAwIbqKtqq6gdVdWmS05O8cJXXn70oyduSvHPV09fw8Tsm\nOSLJbkl2SbJtVb18/OmPJ3l8a+2AJAuzMu6S5LlJXllVP0pyYZKHJ9lr4r86AACAX9+MoQesqrU2\nK1n5mrYkf9xa+9Pxx09IsnuS/6iqSrJrkkuq6qDW2i9X+SmeneRnrbVbxj/un5MckuTzrbWbVznv\n00n+dfzHleQNrbVvPnRfGQAAwIZZryttVXVoVc2rqqur6uS1nPPRqrpm/PbDA8aP7V1VPxq/Wvaj\nqrqtqv7s1x3ZWvtJa22n1trjW2uPS/KLJE9ZLdiSlbdFzqqqLcfj7llJ5o5v2WmV834vyU/Gf/yN\nJK8ffw1cqmqvqtpq1S8ta7iqBwAAsDGs80pbVU1L8rGsDKAbk1xcVf/SWpu3yjmHJdmjtbZXVR2c\n5JNJZrXWrk7ylFV+nl8kOWcCdreMh1RV7Zzk0621F7TWLqqqs5P8KMm94//+1PjHnDoek/cnuS7J\ncePHT8/Kq3iXjofeL5O8ePzn/j9J9snK2yyvT3KsK3IAAMDGVK2t6U0aVzmhalaSd7bWDht/fEqS\n1lp7/yrnfDLJd1prXxp/PDfJWGtt0SrnPDfJO1prvzPxXwYAAMCmaX1uj5yZ5IZVHv9i/NiDnbNg\nDee8LMkXft2BAAAAU9lGeffIqtosyYuSfHljfD4AAIBNxfq8e+SCJI9d5fGu48dWP+cxD3LOYUku\nWe0dHP+Lqnrw+zQBAAA2ca21B7wJ4vpE28VJ9hz/C6dvSnJkkqNWO+erSY5P8qXx18Dduurr2cbP\nX+etket6fR3rZ/bs2Zk9e/bQM2CNfH/SK9+b9Mr3Jj3z/TmxVr4v4gOtM9paa/dV1QlJzsvK2ynP\naK3NrarjVj7dPtVaO7eqDq+qa5MsT3LMKp9466z8+9NeMwFfBwAAwJSyXn+5dmvt61n51verHvu7\n1R6fsJaPvSPJf9vQgQAAAFPZRnkjEjausbGxoSfAWvn+pFe+N+mV70165vtz41jn39O2sVRV62UL\nAADAxlZVa3wjElfaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAA\nOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOiba\nAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAA\nOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOiba\nAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAA\nOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOiba\nAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAA\nOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOiba\nAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAA\nOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOiba\nAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAA\nOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOrZe0VZVh1bVvKq6uqpOXss5H62q\na6rqsqo6YJXjO1TVl6tqblVdUVUHT9R4AACATd06o62qpiX5WJLnJdk/yVFVte9q5xyWZI/W2l5J\njkvyyVWe/kiSc1trv5HkyUnmTtB2AACATd76XGk7KMk1rbX5rbV7k3wxyRGrnXNEkrOSpLV2YZId\nqurRVbV9kt9prX1m/LkVrbXbJ24+AADApm19om1mkhtWefyL8WMPds6C8WOPS/KfVfWZqrq0qj5V\nVVuNMhgAAGAqmbERfv4DkxzfWvthVZ2W5JQk71zTybNnz/7Vj8fGxjI2NvYQzwMAABjGnDlzMmfO\nnHWeV621Bz+halaS2a21Q8cfn5Kktdbev8o5n0zyndbal8Yfz0vyjPGn/7219vjx409PcnJr7YVr\n+DxtXVsAAAA2VVWV1lqtfnx9bo+8OMmeVbVbVW2e5MgkX13tnK8meeX4J5qV5NbW2qLW2qIkN1TV\n3uPnPSvJlRv6RQAAAEw167w9srV2X1WdkOS8rIy8M1prc6vquJVPt0+11s6tqsOr6toky5Mcs8pP\n8WdJ/rGqNkvys9WeAwAA4EGs8/bIjcXtkQAAwFQ2yu2RAAAADES0AQAAdEy0AQAAdEy0AQAAdEy0\nAQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAA\ndEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0\nAQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAA\ndEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0\nAQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAA\ndEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0\nAQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAA\ndEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0\nAQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAA\ndEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0\nAQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAA\ndEy0AQAAdGy9oq2qDq2qeVV1dVWdvJZzPlpV11TVZVX1lFWOX1dV/1FVP6qqiyZqOAAAwFQwY10n\nVNW0JB9L8qwkNya5uKr+pbU2b5VzDkuyR2ttr6o6OMknkswaf/r+JGOttSUTvh4AAGATtz5X2g5K\nck1rbX5r7d4kX0xyxGrnHJHkrCRprV2YZIeqevT4c7WenwcAAIDVrE9MzUxywyqPfzF+7MHOWbDK\nOS3JN6vq4qp69YYOBQAAmIrWeXvkBPjt1tpNVfXfsjLe5rbWvrcRPi8AAMCktz7RtiDJY1d5vOv4\nsdXPecyazmmt3TT+75ur6pysvN1yjdE2e/bsX/14bGwsY2Nj6zEPAABg8pkzZ07mzJmzzvOqtfbg\nJ1RNT3JVVr4RyU1JLkpyVGtt7irnHJ7k+Nba86tqVpLTWmuzqmrrJNNaa8uqapsk5yX569baeWv4\nPG1dWwAAADZVVZXWWq1+fJ1X2lpr91XVCVkZXNOSnNFam1tVx618un2qtXZuVR1eVdcmWZ7kmPEP\nf3SSc6qqjX+uf1xTsAEAALBm67zStrG40gYAAExla7vS5q34AQAAOibaAAAAOibaAAAAOibaAAAA\nOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOiba\nAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAA\nOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOiba\nAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAAOibaAAAA\nOibaYD0tXrw4F1100dAzAACYYkQbrKeXvezYHHzwwfnxj3889BQAAKaQGUMPgMniNa85Og9/+I55\n/OMfP/QUAACmkGqtDb0hSVJVrZctAAAAG1tVpbVWqx93eyQAAEDHRBsAAEDHRBsAAEDHRBsAAEDH\nRBsAAEDHRBsAAEDHRBsAAEDHRBsAAEDHRBsAAEDHRBsAAEDHRBsAAEDHRBsAAEDHRBsAAEDHRBsA\nAEDHRBsAAEDHRBsAAEDHRBsAAEDHRBsAAEDHRBsAAEDHRBuDuv/++9NaG3oGAAB0S7QxmOXLl+fh\nD98pz372i4eeAgAA3RJtDGbatGnZYYeHZ7vttht6CiP67ne/m6222j5nnfW5oacAAGxyZgw9gKlr\nq622yvz584aewQRYsWJF7rprae69996hpwAAbHKql9cTVVXrZQvw62utpaqGngEAMGlVVVprD/gN\nldsjgQkh2AAAHhqiDQAAoGOiDQAAoGOiDQAAoGOiDQAAoGOiDQAAoGOiDQAAoGOiDQAAoGOiDQAA\noGOiDQAAoGOiDQAAoGOiDQAAoGOiDQAAoGOiDQAAoGOiDQAAoGOiDQAAoGOiDQAAoGOiDQAAoGOi\nDQAAoGOibRNzySWX5MlP/u18//vfH3oKAAAwAUTbJuYnP/lJLr/8glx22WVDT2FEn/jE6Tn88D/I\n8uXLh54CAMCAqrU29IYkSVW1XrZMZq21XHvttdlzzz1TVUPPYQRjYy/I+ef/W6699trsscceQ89h\nA91///057rg3Ze+998hJJ71x6DkAQMeqKq21B/wmXrRBp2677bYsXLgw++yzz9BTGMHtt9+eHXbY\nIbvttl+uu+6KoecAAB0TbQADufLKK7P99ttn1113HXoKANAx0QYAANCxtUWbNyIBAADomGgDAADo\nmGgDAADomGgDAADomGgDAADomGgDAADomGgDAADomGgDAADomGgDAADomGgDAADomGgDAADomGgD\nAADomGgDAADomGgDAADomGgDAADomGgDAADo2HpFW1UdWlXzqurqqjp5Led8tKquqarLquqA1Z6b\nVlWXVtVXJ2I0AADAVLHOaKuqaUk+luR5SfZPclRV7bvaOYcl2aO1tleS45J8crWf5o1JrpyQxQAA\nAFPI+lxpOyjJNa21+a21e5N8MckRq51zRJKzkqS1dmGSHarq0UlSVbsmOTzJ6RO2GjbAfffdl1NP\n/WC+973vDT0FAADW2/pE28wkN6zy+Bfjxx7snAWrnPPhJCclaRu4ESbEvHnzcvLJJ+XEE/9q6CmM\n6NZbb80HP/jBLFy4cOgpAAAPuYf0jUiq6vlJFrXWLktS4//AIPbbb7/8wz/8Q84886NDT2FEn//8\n53PSSSflE5/4u6GnMKIrr7wyn/vc59KaP9cDgLWZsR7nLEjy2FUe7zp+bPVzHrOGc34/yYuq6vAk\nWyXZrqrOaq29ck2faPbs2b/68djYWMbGxtZjHqyfqsrRRx899AwmwFFHHZVly5blFa94xdBTGNGr\nX/3mXHDBeXniE5+YJz/5yUPPYQTf+c53suOOO+YpT3nK0FMAJo05c+Zkzpw56zyv1vWnm1U1PclV\nSZ6V5KYkFyU5qrU2d5VzDk9yfGvt+VU1K8lprbVZq/08z0jyltbai9byeZo/aQWYWi644IKcf/53\nc9JJb8mMGevz54j0aPny5dl2223zqEftnkWLfj70HIBJq6rSWnvA3Ynr/BWytXZfVZ2Q5LysvJ3y\njNba3Ko6buXT7VOttXOr6vCqujbJ8iTHTPQXAMCm55BDDskhhxwy9AxGtM022+T97/9QZs7ceegp\nAJukdV5p21hcaQMAAKaytV1pe0jfiAQAAIDRiDYAAICOiTYAAICOiTYAAICOiTYAAICOiTYAAICO\niTYAAICOiTYAAICOiTYAAICOiTYAAICOiTYAAICOiTYAAICOiTYAAICOiTYAAICOzRh6AAAMaenS\npTnnnHOycOHC7LTTTnnJS16S7bbbbuhZAPAr1VobekOSpKpaL1sA2PS11vKh970v7333u/M706dn\nr7vuyjVbbpnv3ndf3vaOd+Qtp5ySqhp6JgBTSFWltfaAX3xcaWMky5cvz3e/+908+9nPzowZvp2A\nyeND73tfPvue9+TSO+/M7v/v4LJluS7JC9/zniTJn7/1rcOMA4BVuNLGSP7yL/8q733vu/PZz342\nr3zlK4eewwhuvPHG3Hrrrdlvv/2GngIPuaVLl2a3Rz/6vwbbKq5L8ptbb535ixZl22233bjjAJiy\n1nalzRuRMJI/+IPfy8te9sd55jOfOfQURjQ29vzsv//+WbJkydBTGNHPf/7z3HLLLUPP6No555yT\n35k+fY3BliS7J3n6tGk555xzNt4oAFgL97MxkgMOOCBf/OLfDz2DCfDa1/5pfvzjudlhhx2GnsII\nbrnlljz+8Y/P/vsfnJ/85AdDz+nWwoULs9dddz3oOXvddVduuummjbTogVpr+dnPfpbHPe5xmTbN\nn7ECTGV+FQCSJCee+IZ85jMf95vDSW777bfPi1/8RznyyJcOPaVrO+20U67ZcssHPeeaLbfMzjvv\nvJEWPdCXv/zl7LnnnvngB/9msA0A9MFr2gCYcibDa9ouv/zyHH30a3Paae/J7/7u7w6yAYCNy2va\nAGDcdtttl7e94x154dZb57rVnrsuyQu33jpvffvbB30Tkic96Um5/PILBBsAXtMGwNT0llNOSZIc\nuKa/p+3tb//V8wAwNLdHAjClLV26NF/5yldy0003Zeedd85LXvISb/MPwCDWdnukaAMAAOiA17QB\nAABMQqINAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACg\nY6INAACgY6INAACgY6INAACgY6INAACgY6INHkKttaEnAKtYtmxZrr766ixbtmzoKQCw3kQbPERu\nuOGGbLXVtnnta9889BSY8lasWJE3vOUtedTMmfnN5zwnj5o5M294y1uyYsWKoacBwDqJNniITJs2\nLVtttW222GKLoacwoi996X9l220fnvPPP3/oKWygN598cs783vdy55lnZtlnP5s7zzwzZ37ve3nz\nyScPPQ0A1km0wUNk5syZWbJkUT7ykfcNPYUR3XHHnVm+fEnuvvvuoaewAZYtW5YzTj89d/zFXyRf\n+EryR69ONt88d/zFX+SMM85wqyQA3RNtAOtwzDF/nBUrVuS5z33u0FPYADfeeGOm77hj8ohHJAt/\nmSy8LrnnnuQRj8j0HXbIjTfeOPREfk133313nvSkp+Xoo48begrARjFj6AEAk8H06dOHnsAG2mWX\nXXLfrbcmixcn7zo5uffEZIstksWLc99tt2WXXXYZeiK/phUrVuTaa3+SrbfeZugpABuFK20AbNK2\n3XbbHPuqV2XrU09Nliz5VbBtfeqpOfbYY7PtttsOPZFf0zbbbJPFi3+Z733v60NPAdgoqpe3JK+q\n1ssWADYtK1asyJtPPjlnnHFGpu+wQ+677bYce+yx+fD7358ZM9x0AkAfqiqttXrA8V5CSbQB8FBb\ntmxZbrzxxuyyyy6usAHQHdEGAADQsbVFm9e0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0\nAQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0Mal8\n4xvfyFOfOpZrrrlm6CkAALBRiDYmlfPPvyCXXHJ+5s2bN/QURvTOd/5/ecUrjsv9998/9BQAgK5V\na23oDUmSqmq9bKFfK1asyM9+9rPsvffeQ09hRLvvvn/mz78yy5YtyzbbbDP0HDbQnXfemT/90xNy\nxBHPy5FH/uHQcwBgUquqtNZq9eOutDGpzJgxQ7BtIn7wg/+da6+9VrBNcvPnz88Xv3hmPv7xvx96\nCiOaO3dujj76VZk/f/7QUwBYzYyhBwBT00477TT0BCbAvvvumwsvvDCPe9zjhp7CiP75n8/JP/7j\nGXn605+a1772tUPPAWAVbo8EAHLHHXfkW9/6Vg499NBsvvnmQ88BmJLWdnukaAMAAOiA17QBAABM\nQqINAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACgY6IN\nAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INNhF33nln/sf/eG+uvPLK\noacAADCBRBtsIr797W/n7W//y7z73R8YegojWrBgQU477bQsW7Zs6CkAQAdEG2winvOc5+QTn/hE\n3vvevxp6CiP60Ic+kje/+c05++yzh57CiC6++OJ87WtfG3oGAJNctdaG3pAkqarWyxaAIV1//fX5\n/Oc/n9e//vXZfvvth57DCHbddZ8sWHB1lixZkh133HHoOYzg3HPPzT777JM99thj6CnAJqyq0lqr\nBxzvJZREGwCbmn/5l6/muuuuzxvfeMLQUxjBVVddlX333Te/9VvPzEUXfXvoOYzg3nvvzXnnnZdn\nPvOZ2XrrrYeeAw8g2gAANsCKFSvytrfNztjYb+fwww8beg4j+PSnP53XvOY1mT37XXnnO98x9Bx4\nANEGAMCUNn/+/Lz1rX+dt73txDzhCU8Yeg48gGgDAADo2NqizbtHAgAAdEy0AQAAdEy0AQAAdEy0\nAQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAA\ndEy0AQAAdGy9oq2qDq2qeVV1dVWdvJZzPlpV11TVZVV1wPixLarqwqr6UVX9uKreOZHjAQAANnXr\njLaqmpbkY0mel2T/JEdV1b6rnXNYkj1aa3slOS7JJ5OktXZ3kme21p6S5IAkh1XVQRP7JcBDZ+nS\npbn44ouHngEAwBS2PlfaDkpyTWttfmvt3iRfTHLEaucckeSsJGmtXZhkh6p69PjjO8bP2SLJjCRt\nIobDxvCqV/1ZDjrooHz/+98fegojuummm7Jw4cKhZwAA/NpmrMc5M5PcsMrjX2RlyD3YOQvGjy0a\nv1J3SZI9kvxta81lCyaNP/qj38vy5Xdk3333XffJdKu1lr32+o1Mn75Zbrvt5qHnMKJrrrkmu+yy\nS7bZZpuhpwDARvGQvxFJa+3+8dsjd01ycFXt91B/TpgoL3rRC/O1r30pj3jEI4aewgiqKkce+coc\neeTRQ09hRFdccUX23nvvHHXUq4aewojuv//+XHfddUPPAJgU1udK24Ikj13l8a7jx1Y/5zEPdk5r\n7faq+k6SQ5NcuaZPNHv27F/9eGxsLGNjY+sxD2DdTj/9o0NPYALMnDkzY2MvyBFHHDr0FEb0gQ/8\nTU455aScffbZeelLXzr0HEZw99135+677872228/9BSYdObMmZM5c+as87xq7cFfYlZV05NcleRZ\nSW5KclGSo1prc1c55/Akx7fWnl9Vs5Kc1lqbVVWPTHJva+22qtoqyTeSvK+1du4aPk9b1xYAYNPw\nzW9+M29609vzhS98Ok960pOGnsMIDjrod/OjH30/ixffLNxgRFWV1lqtfnydV9paa/dV1QlJzsvK\n2ynPaK3NrarjVj7dPtVaO7eqDq+qa5MsT3LM+IfvnOSz469rm5bkS2sKNgBgannOc56TK654ztAz\nmAAHHnhAWmvZYosthp4Cm6x1XmnbWFxpAwAAprK1XWl7yN+IBAAAgA0n2gAAADom2gAAADom2gAA\nADom2gAAADom2gAAADom2gAAADom2gAAADom2gAAADom2gAAADom2gAAADom2gAAADom2gAAADom\n2gAAADom2gAAADom2gAAADom2gAAADom2oANcskll2SzzbbIhz700aGnAABs0kQbsEGqKtOnz8j0\n6f4zMtmdeuppefjDd8lVV1019BQAYA1mDD0AmJwOPPDA3HXX8qFnMAFuvfW2LFlyU+6+++6hpzCi\n5z7393LPPffkO9/511TV0HMAmCDVWht6Q5KkqlovWwCmmhUrVmTGDH+ON9k95jF755577srChfNF\n2yR2881wcKsTAAAbTklEQVQ352lPe3Ze/vKX5V3vetvQc4CNqKrSWnvAf8Dd1wSAYNtEXHvtj3P9\n9dcItklu+fLl+elPL8/VV1879BRG9M1vfjN77nlAfvjDHw49hUnOr9IAsInYYosthp7ABNh9992z\ndOnSbL311kNPYURXXjkvP/3pf+TnP/95nvrUpw49h0nM7ZEAAPAQaK1l4cKF2XnnnYeewiSxttsj\nRRsAAEAHvKYNAABgEhJtAAAAHRNtAAAAHRNtAAAAHRNtAAAAHRNtAAAAHRNtAAAAHRNtAAAAHRNt\nAAAAHRNtAAAAHRNtAAAAHRNtAAAAHRNtAAAAHRNtAAAAHRNtAAAAHRNt0JHTT//7HHHEy3PHHXcM\nPQUAgE6INujIWWf9r3z1q1/IjTfeOPQURtBay+ted2JOOWX20FMAgE1AtdaG3pAkqarWyxYYypIl\nS7JgwYI84QlPGHoKI1ixYkW23HLr7Ljjo/Kf//mLoecwgptvvjmve91b8oY3HJtnPOMZQ88BYBNX\nVWmt1erHXWmDjjzsYQ8TbJuAGTNm5Oqr5+VHP/r3oacwoksuuST/9E//kM9+9ktDT2FE//7v/55X\nver43HrrrUNPAfi1udIGAGvRWsv555+fAw88MNtvv/3QcxjBK17x6nzuc6fn3HPPzWGHHTb0HEZw\nxhmfzXXXXZ93vevtqXrABQmY1NZ2pU20AQCbvMWLF+eCCy7I85///Eyb5kajyWzXXffOggXX5Pbb\nb89222039ByYUKINAIBJb+7cuVmyZEkOOeSQoafAhBNtAAAAHfNGJAAAAJOQaAMAAOiYaAMAAOiY\naAMAAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMA\nAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMmhdtuuy0f+chH8stf/nLoKQAAG5VoAyaFz33uc3nTm96U\nj3/8E0NPYURXX311/u7v/i4rVqwYegoATAozhh4AsD6OPPLI3HrrrTnmmGOGnsKITjzx7fm3f/ty\n9tlnn4yNjQ09hxHMmTMnM2bMyNOf/vShpwBs0qq1NvSGJElVtV62APDQufzyy/P1r38jb3rTG7P5\n5psPPYcN1FrL9OnTs8UW2+bOO28feg4juP/++/PVr341hxxySB71qEcNPQemtKpKa60ecLyXUBJt\nADC5nH76mdlss83yx3/8iqGnMIJvfOMbOfTQQ3PkkcfkC184c+g5jOC2227LRRddlGc961mZNs2r\noCYj0QYAwAPcfvvtOeWUv8orX3lkZs2aNfQcRvDqV5+Q00//2/zrv/5rXvCCFww9hw2wtmjzmjYA\ngCls++23z8c/ftrQM5gAf/InR2XZsjty8MEHDz2FCeZKGwAAQAfWdqXNza4AAAAdE20AAAAdE20A\nAAAdE20AAAAdE20AAAAdE20AAAAdE20AAAAdE20AAAAdE20AAAAdE20AAAAdE20AAAAdE20AAAAd\nE20AAAAdE20AAAAdE20AAAAdE20whS1evDhXXnnl0DMAAHgQog2msMMO+4Psv//+ueGGG4aewoiu\nv/763HrrrUPPAAAeAqINprA/+ZOj8tKXviKPetSjhp7CCG655Zbstttu+e///fChpzABrr/++rTW\nhp4BQEdEG0xhr3/9q3P22Wdliy22GHoKI9huu+3y/Of/YV7ykhcOPYURnX322dltt93ywQ9+eOgp\njGjFihW58cYbh54BbCJEG8Akt9lmm+VrX/tS/vqv3zr0FEa05557Zt99n5oDDnji0FMY0eted2Jm\nzpyZH/7wh0NPYUR33nln7rzzzqFnMMXNGHoAALDSAQcckLlzLx56BhNg1qzfzA9+MCs77bTT0FMY\n0eMet2+SZOHC+QMvYSqrXu6br6rWyxYAAEiS5z3vpUmSb3zjnwZewlRQVWmt1QOO9xJKog0AAJjK\n1hZtXtMGAADQMdEGAADQMdEGAADQMdEGAADQMdEGAADQMdEGAADQMdEGAADQMdEGAADQMdEGAADQ\nMdEGAADQMdEGAADQMdEGAADQMdEGAADQMdEGAADQMdEGAADQMdEGAADQsfWKtqo6tKrmVdXVVXXy\nWs75aFVdU1WXVdUB48d2rapvV9UVVfXjqvqziRwPAACwqVtntFXVtCQfS/K8JPsnOaqq9l3tnMOS\n7NFa2yvJcUk+Of7UiiQnttb2T/K0JMev/rEAPTrvvPOy446Pzrnnnjv0FABgilufK20HJbmmtTa/\ntXZvki8mOWK1c45IclaStNYuTLJDVT26tbawtXbZ+PFlSeYmmTlh6wEeIkuXLs1tt/0yS5cuHXoK\nIzr++JPymMfsk1tuuWXoKQCwQWasxzkzk9ywyuNfZGXIPdg5C8aPLfp/B6pq9yQHJLlwA3YCbFQv\nfelLc/fdd2fzzTcfegojWrTol1m06Prcc889Q09hBPfcc09mzXp2nvzkJ+Uzn/nY0HMANqqN8kYk\nVbVtkrOTvHH8ihtA9wTbpuHLX/77LFt2W3baaaehpzCCe+65J1dc8cP85CdXDD2FEV166aXZddd9\ncs45Xxl6Ckwa63OlbUGSx67yeNfxY6uf85g1nVNVM7Iy2P6htfYvD/aJZs+e/asfj42NZWxsbD3m\nAcDaVZUA3wRsu+22ueWWm/1vuQlYtGhRFiy4Oj/72XVDT4HBzZkzJ3PmzFnnedVae/ATqqYnuSrJ\ns5LclOSiJEe11uaucs7hSY5vrT2/qmYlOa21Nmv8ubOS/Gdr7cR1fJ62ri0AAEx+S5YsyY477piq\nGnoKdKWq0lp7wP8x1nmlrbV2X1WdkOS8rLyd8ozW2tyqOm7l0+1TrbVzq+rwqro2yfIkfzL+SX87\nyR8l+XFV/ShJS/K21trXJ+wrAwBgUnnYwx429ASYVNZ5pW1jcaUNAACYytZ2pW2jvBEJAAAAG0a0\nAQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAA\ndEy0AQAAdEy0AQAAdEy0AQAAdEy0AQAAdEy0ARvNe97zgRx77AlprQ09BQBg0qhefvNUVa2XLcBD\nY7fdfiPXXz8vy5YtyzbbbDP0HDbQihUrcvzxb8msWU/NMce8Yug5ALDJqKq01mr14660ARvNBRd8\nK1dddZVgm+QWLVqUT33qo/mbv/n40FMY0U9/+tP8/u+/IldcccXQUwB4EDOGHgBMHTNnzhx6AhNg\n5syZufDCC7PLLrsMPYURnXfeefmnf/pcnvrUJ2T//fcfeg4j+Ld/Ozff+tb5ef/7353NN9986DnA\nBHN7JABMUffee2++/e1v5xnPeEa23HLLoecwgqc97Tn5wQ++lSuuuCL77bff0HOADbS22yNFGwDA\nJDd//vzMmzcvz3ve84aeAoxAtAEAAHTMG5EAAABMQqINAACgY6INAACgY6INAACgY6INAACgY6IN\nAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACgY6INAACg\nY6INAACgY6INAACgY6INoAM33XRT/vZv/zbLly8fegoA0BnRBtCBD3zgwznhhBNy9tlnDz2FEV16\n6aX5+te/PvQMADYhM4YeAEDyhje8LjvuuF2OOOKIoacwohe96KgsWHB1lixZkh133HHoOYxgzpw5\nmTlzZvbaa6+hpwBTXLXWht6QJKmq1ssWANhQZ5/9z/n5z+fnz//8TamqoeewgRYuXJidd945e+31\nlFx99aVDzwGmiKpKa+0Bv3iINgCA1dx///1529tm54ADnpAjj/zDoecAU4RoAwAA6Njaos0bkQAA\nAHRMtAEAAHRMtAEAAHRMtAEAAHRMtAEAAHRMtAEAAHRMtAEAAHRMtAEAAHRMtAEAAHRMtAEAAHRM\ntAEAAHRMtAEAAHRMtAEAAHRMtAEAAHRMtAEAAHRMtAEAAHRMtAEAAHRMtAGbrMWLF2fhwoVDzwAA\nGIloAzZZBxzwtDz2sbvnnnvuGXoKI1q0aFHuvffeoWcAwCBEG7DJevGLX5wjjjgym2222dBTGMG8\nefOy00475eUvf9XQUxhRay0LFiwYegbApCPagE3W//yfp+bLX/77VNXQUxjBIx/5yDz5yb+dpz3t\nt4aewog++MEPZ9ddd81XvvKVoacATCozhh4AAA/mkY98ZC677HtDz2ACPPGJ+2WvvZ6S3Xfffegp\nAJNKtdaG3pAkqarWyxYAAICNrarSWnvALUJujwQAAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMAAOiY\naAMAAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMA\nAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAMAAOiYaAOYgs4887N55CMfk8suu2zoKQDAOog2gCno5psX\nZ/HiX2Tp0qVDT2FEL3/5q3PIIc/NfffdN/QUAB4i1VobekOSpKpaL1sApoK77rorW2655dAzGNH+\n+x+U6667Krfc8stsscUWQ88BYARVldZaPeB4L6Ek2uD/tnf/MZadZR3Avw/dQrBEBAyl6UKBAtYQ\nSov0R0oDK0jcEqU0ElM0ChhjU0ArIIKthvVH0lIrSCFawEqKhrQEGygEYkG6Jhtt2QBbCunSJZD+\nUtagLEIrm2738Y+5hWG6szOz8+udmc8nmeSe977vPc9Jnr0733vPnAOwcPv378+BAwdyzDHHrHYp\nACyS0AYAADCw2UKbv2kDAAAYmNAGAAAwMKENAABgYEIbAADAwIQ2AACAgQltAAAAAxPaAAAABia0\nAQAADExoAwAAGJjQBgAAMDChDQAAYGBCGwAAwMCENgAAgIEJbQAAAAMT2gAAAAYmtAEAAAxMaAMA\nABiY0AYAADAwoQ1giVxyyZ/l8sv/erXLAADWmXmFtqraWlW7q+qOqnrrLHOurKo9VbWrqk6dNn51\nVe2tqi8vVdEAozlw4EAuvXRbLr30stUuhUXat29fLrzwouzcuXO1SwGAJMmmuSZU1SOSvDfJS5L8\nR5KdVfXx7t49bc45SU7s7mdW1RlJ/jbJmZOnP5jkPUk+tNTFA4xi06ZN2bVrVx71qEetdiks0o4d\nO3LVVVfm/vt/kGuuOW21ywGAuUNbktOT7OnuO5Okqq5Ncm6S3dPmnJtJKOvuW6rqsVV1bHfv7e4d\nVXXCUhcOMJqTTz55tUtgCZxzzjm5/vrrc/bZZ692KQCQZH6h7fgkd0/bvidTQe5wc+6djO1dVHUA\nsMKOOuqonHfeeatdBgD80HxC24rZtm3bDx9v2bIlW7ZsWbVaAAAAltP27duzffv2OedVdx9+QtWZ\nSbZ199bJ9tuSdHe/Y9qcq5Lc1N3XTbZ3J3lRd++dbJ+Q5BPdPeu5Q1XVc9UCAACwXlVVurtmjs/n\n6pE7kzyjqk6oqkcmOT/JDTPm3JDkNyc7OjPJvocC20P7n/wAAACwAHOGtu5+MMkbktyY5KtJru3u\n26vqgqr6ncmcTyX5ZlV9Pcn7krzuofVV9eEk/5bkWVV1V1W9dhmOAwAAYF2a8/TIleL0SAAAYCNb\nzOmRAAAArBKhDQAAYGBCGwAAwMCENgAAgIEJbQAAAAMT2gAAAAYmtAEAAAxMaAMAABiY0AYAADAw\noQ0AAGBgQhsAAMDAhDYAAICBCW0AAAADE9oAAAAGJrQBsOzuuuuuXHfddTl48OBqlwIAa47QBsCy\ne/3r35Lzzz8/O3bsWO1SAGDN2bTaBQCw/l188e/npJNOzGmnnbbapQDAmlPdvdo1JEmqqkepBQAA\nYKVVVbq7Zo47PRIAAGBgQhsAAMDAhDYAAICBCW0AAAADE9oAAAAGJrQBAAAMTGgDAAAYmNAGAAAw\nMKENAABgYEIbAADAwIQ2AACAgQltAAAAAxPaAAAABia0AQAADExoAwAAGJjQBgAAMDChDQAAYGBC\nGwAAwMCENgAAgIEJbQBr1L59+3LfffetdhkAwDIT2gDWoP379+e44zbn1FNfsNqlAADLbNNqFwDA\nwh199NE566yX5MlP3rzapQAAy6y6e7VrSJJUVY9SCwAAwEqrqnR3zRx3eiQAAMDAhDYAAICBCW0A\nAAADE9oAAAAGJrQBAAAMTGgDAAAYmNAGAAAwMKENAABgYEIbAADAwIQ2AACAgQltAAAAAxPaAAAA\nBia0AQAADExoAwAAGJjQBgAAMDChDQAAYGBCGwAAwMCENgAAgIEJbQAAAAMT2gAAAAYmtAEAAAxM\naAMAABiY0AYAADAwoQ0AAGBgQhsAAMDAhDYAAICBCW0AAAADE9oAAAAGJrQBAAAMTGgDAAAYmNAG\nAAAwMKENAABgYEIbAADAwIQ2AACAgQltAAAAAxPaAAAABia0AQAADExoAwAAGJjQBgAAMDChDQAA\nYGBCGwAAwMCENgAAgIEJbQAAAAMT2gAAAAYmtAEAAAxMaAMAABiY0AYAADAwoQ0AAGBgQhsAAMDA\nhDYAAICBCW0AAAADE9oAAAAGJrQBAAAMTGgDAAAYmNAGAAAwMKENAABgYEIbAADAwIQ2AACAgQlt\nAAAAAxPaAAAABia0AQAADExoAwAAGJjQBgAAMDChDQAAYGBCGwAAwMCENgAAgIEJbQAAAAMT2gAA\nAAYmtAEAAAxsXqGtqrZW1e6quqOq3jrLnCurak9V7aqqUxayFgAAgEObM7RV1SOSvDfJLyZ5dpJX\nVdVJM+ack+TE7n5mkguSXDXftSy97du3r3YJMCv9yaj0JqPSm4xMf66M+XzTdnqSPd19Z3c/kOTa\nJOfOmHNukg8lSXffkuSxVXXsPNeyxPzjYWT6k1HpTUalNxmZ/lwZ8wltxye5e9r2PZOx+cyZz1oA\nAABmsVwXIqllel0AAIANpbr78BOqzkyyrbu3TrbflqS7+x3T5lyV5Kbuvm6yvTvJi5I8ba61017j\n8IUAAACsc939sC/ANs1j3c4kz6iqE5L8Z5Lzk7xqxpwbkrw+yXWTkLevu/dW1bfnsXbW4gAAADa6\nOUNbdz9YVW9IcmOmTqe8urtvr6oLpp7u93f3p6rqZVX19ST3JXnt4dYu29EAAACsM3OeHgkAAMDq\nWa4LkbDMqupxVXVjVX2tqv65qh47y7yrq2pvVX35SNbDQi2gN7dW1e6quqOq3jpt/LlV9e9V9aWq\n+nxVPX/lqme9W2x/Tp773aq6vapuq6rLVqZy1rul6M3J82+uqoNV9fjlr5qNYgn+b7988r65q6r+\nqap+cuWqXx+EtrXrbUk+290/k+RzSf5olnkfzNTNzY90PSzUnL1VVY9I8t5M9eazk7yqqk6aPH15\nkrd396lJ3p7kL1ekajaKRfVnVW1J8stJntPdz0lyxQrVzfq32PfOVNXmJC9NcueKVMxGstj+vDHJ\ns7v7lCR7DrWewxPa1q5zk1wzeXxNklccalJ370jynSNdD0dgPr11epI93X1ndz+Q5NrJuiQ5mOSh\nT/B+Ksm9y1grG89i+/PCJJd194Ek6e5vL3O9bByL7c0keVeStyxrlWxUi+rP7v5sdx+czLs5yeZl\nrnfdEdrWrid2994k6e5vJXniCq+H2cynt45Pcve07XsmY0nyxiRXVNVdmfrWzadxLKXF9uezkryw\nqm6uqpucvssSWlRvVtXLk9zd3bctd6FsSIt975zut5J8eskrXOfmc8l/VklVfSbJsdOHknSSPz7E\n9MVeUcYVaZi3Ze7NC5Nc1N0fq6pXJvn7TJ3uA/OyzP25KcnjuvvMqjotyUeSPP2ICmXDWa7erKpH\nJ7k4P/5e6VZKLMhK/N5ZVZckeaC7P3wk6zcyoW1g3T3rL6qTi4scO7kf3pOS/NcCX36x69nAlqA3\n703ylGnbm/Oj0yBf3d0XTfbz0aq6eqnqZmNY5v68J8n1k/3snFzw4Qnd/d9LVD7r2DL25olJnprk\n1qqqyfgXqur07vb/O/OyzO+dqarXJHlZkhcvTcUbi9Mj164bkrxm8vjVST5+mLmVh3/itpD1sBDz\n6a2dSZ5RVSdU1SOTnD9t3r1V9aIkqaqXJLljectlgznS/rxh8tzHMvmFo6qeleRogY0lcsS92d1f\n6e4ndffTu/tpmfpw4VSBjSW0qPfOqtqaqb+3fHl371/+ctcf92lboyaX8v1Ikidn6ipRv9rd+6rq\nuCQf6O5fmsz7cJItSZ6QZG+mrsr3wdnWr/yRsN4soDe3Jnl3pj48urq7L5uMn5XkyiRHJflBktd1\n95dW/khYj5agP4/O1Cm7pyTZn+TN3f2vK38krDeL7c0Zr/WNJM/v7v9ZsQNgXVuC9849SR6Z5KEP\nuW7u7tet8GGsaUIbAADAwJweCQAAMDChDQAAYGBCGwAAwMCENgAAgIEJbQAAwLpQVa+sqq9U1YNV\n9bxZ5myuqs9V1Ver6raq+r1DzHnz5F6cj58x/pSq+l5VvekQa26oqi/Po8Zfq6pbJz87quo5c61x\nc20AAGC9uC3JeUned5g5B5K8qbt3VdVjMnUz+hu7e3cyFeqSvDRTtzeY6a+SfGrmYFWdl+R/51nj\nN5K8sLu/O7lNwgeSnHm4Bb5pAwAA1oXu/lp370lSh5nzre7eNXn8/SS3Jzl+2pR3Zepm4D+mqs7N\nVOD66ozxY5K8MclfzBj/6ar6aFXdMvk5a7LPm7v7u5NpN8/Y9yEJbQAAwIZUVU9NckqSWybbL09y\nd3ffNmPeMUn+MMmf5uGB8M+TXJHk/2aMvzvJO7v7jCSvTPJ3hyjht5N8eq46nR4JAACsGVX1mSTH\nTh9K0kku6e5PLOB1HpPko0ku6u7vV9Wjk1ycqVMjZ9qW5F3dfX9VTX+N5yY5sbvfNAmA0wPdLyT5\n2frRgsdU1U909/2TtT+f5LVJzp6rVqENAABYM7r7UKFqQapqU6YC2z9098cnwycmeWqSWydBa3OS\nL1bV6UnOSPIrVXV5ksclebCqfpDkYJKfq6pvJDk6yROr6nPd/eJMBbgzuvuBQ+z/5CTvT7K1u78z\nZ73dvbgjBgAAGEhV3ZTkD7r7C7M8/6Ek3+7uh10FctqcbyZ53sxQVVVvT/K97n7njPETknyiu0+e\nbP9jkl3dfcVk+7ndfWtVPSXJvyT5je6+eT7H42/aAACAdaGqXlFVd2fqaoyfrKpPT8aPq6pPTh6/\nIMmvJ3lxVX2pqr44uYrjTJ3DXNBkHi5K8vzJpf2/kuSCyfifJHl8kr+Z7P/zcx6Xb9oAAADG5Zs2\nAACAgQltAAAAAxPaAAAABia0AQAADExoAwAAGJjQBgAAMDChDQAAYGBCGwAAwMD+Hy+c8MgsV1Ze\nAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#### print(find_closest_model_point(-124, 48, model_lons = X, model_lats = Y, grid = 'Test', tols = {\"Test\":{'tol_lat':0.002658843994140625, 'tol_lon':0.0009403228759765625}}, land_mask = bathy.mask))\n", "lon, lat = -124.5, 48.538\n", "#closest_point = find_model_point(lon, lat, X, Y, tol_lon = 0.0052, tol_lat = 0.00189)\n", "closest_point = find_closest_model_point(lon, lat, model_lons = X, model_lats = Y, grid = 'NEMO', land_mask = bathy.mask)\n", "print(closest_point)\n", "\n", "r1 = max(0, closest_point[0] - 5)\n", "r2 = min(897, closest_point[0] + 5)\n", "c1 = max(0, closest_point[1] - 5)\n", "c2 = min(397, closest_point[1] + 5)\n", "\n", "fig, ax = plt.subplots(1, 1,figsize=(15,15))\n", "#plt.pcolormesh(X[r1:r2,c1:c2],Y[r1:r2,c1:c2], (bathy.mask).astype('int')[r1:r2,c1:c2])\n", "#plt.pcolormesh(X, Y, (bathy.mask).astype('int'))\n", "\n", "plt.scatter([X[closest_point]], [Y[closest_point]], c = \"c\" , s =40)\n", "#plt.scatter([X[109,143]], [Y[109,143]], c = \"c\" , s =40)\n", "plt.scatter(X[r1:r2,c1:c2],Y[r1:r2,c1:c2], np.invert(bathy.mask)[r1:r2,c1:c2])\n", "plt.scatter([lon], [lat], c = \"r\" , s =80)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4MAAAN5CAYAAACynu6FAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XGM53d93/nX2/iISJtzTCpw5RlIIRSyXFMTZcGn3Inp\nWRFro8YR6kW4dyKYE7WuuEFJtsKg03i1ak84cikgX+ryw0Nx1R7kouhiWhIogmkvTWKmBScka4PR\nHWbGwZtExORCdIHA+/6YH/Yw7M7M7v5mfzPzeTykEfP7/j6/3+/zQz+t/NTnN59PdXcAAAAYyxXz\nngAAAACXnxgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAY0JGJwar6uar6ZlU9+zz3X1VV/0dV\nPVxVv1dVr5he/5tV9ZtV9emq+mRVHZ9ev7Kq/kVV/c50/B17mMObqurRqvrG+eYBAABwEByqGKyq\nV1bV+85xfSHJjyV5bIeHvyvJh7v7B5P8zSQPT6//fJI7u/tlSe6c3k6S/z7JM7v7h5L8SJLbqup5\nu0zx15PcsMs8AAAA5u5QxeBUn+PaP03yD8/3gKr6L5P8t939viTp7r/o7j+Z3v3NJFdNf//eJI9v\neZ2/VFXPSPLdSf48yZ9Mn+/Hquo3quo/VdUHq+q7p8/72939xSR1Se8QAABgnx3GGPy20KqqH0+y\n3t2f2eExfy3JH1XV+6rqU1X1nqp61vS+n0lyd1V9MZurgm+dXv+lJH+W5EtJvpDk7u5+sqq+L8n/\nkuSG7v6RJP85yc/N6L0BAABcFlfOewJ7UVW/leSZSb4nydVV9anpXaeSvC2bXxF9avg5nuLKJD+c\n5E3d/Z+q6p1J7sjm10L/5yRv7u7/s6r+TpKV6fO9IslfJLkmyfcl+b+q6mNJXprkWJL/WFWV5L9I\n8pszfLsAAAD7rrrP9a3Lg6mqXpnkp7r7DdPb/1WSj2VzBa+SLGTza54v7+4/2PK45yb5ze5+wfT2\nf5PkLd39t6vqye7+3i1jn+zu762qe6aP+VfT6/cl+dUk/1+SW7r7f9hhnv93kh/p7i/P8v0DAADM\nyp6+JlpVJ6rqkar6XFW95Txj3j3dSfOhqrpu231XTL+e+cCWa3dW1cb0+qeq6sSFTr67f7e7r+nu\nF3T3X0uykeRlW0NwOu5skvWq+uvTSzckOTP9/fFpZKaqbkjy6PT6F5P8d9PrfynJ9UkeSfJbSX60\nql44ve+7q+pF2//viL8bBAAADrBdY7CqrkhyT5JXZfMrkrdU1Uu2jbkxyQu7+0VJbkty77aneXOe\njq+t3tHdPzz9+bWLeQPbdKYRVlV/tar+zZb7fjrJv6qqh7K5m+j/Or3+95L8k6r6dJJ/NL2dJP9b\nku+pqt9N8mCS+6bx+UdJXp/kf6+q307yG0lePH3Nf1BV60muTfLbVfWeGbwnAACAmdv1a6JVdX02\nj164cXr7jiTd3XdtGXNvkk909wentx9OstTdZ6fHPrwvyT9O8rPd/ePTMXcm+dPu/if78L4AAADY\nwV6+JnptkvUttzem13Ya8/iWMd869uFc1Xn79Gul762qq85xPwAAAPtgX4+WqKpXJznb3Q/lO/+O\n7heSvKC7r0vyRJJ37OdcAAAAeNpejpZ4PMnzttz+1o6d28csnmPM30ny41V1U5JnZfNv8O7v7td1\n9x9uGT9J8qFzvXhVHZ7tTgEAAPZBd898g8q9xOBakh+oqudn8wD21ya5ZduYB5K8KckHp39j+OR0\nB8+3TX++dSzEz3X366a3r+nuJ6aPf02S3z3fBA7T8RfMz6lTp3Lq1Kl5T4NDwueFvfJZ4UL4vLBX\nPitciM3jzWdv1xjs7m9U1e1JPprNr5Xe190PV9Vtm3f3e7r7w1V1U1V9PslXk9y6h9f++ekRFN9M\n8oVs7kIKAADAZbCXlcFMj3148bZr/3zb7dt3eY5/n+Tfb7n9ur1PEwAAgFna1w1k4HJaWlqa9xQ4\nRHxe2CufFS6Ezwt75bPCQbDrOYPzVlV90OcIAACwX6pqXzaQsTIIAAAwIDEIAAAwIDEIAAAwIDEI\nAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAw\nIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEI\nAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAw\nIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEI\nAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAw\nIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEI\nAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAw\nIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEI\nAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAw\nIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwoD3FYFWdqKpHqupzVfWW84x5d1U9WlUPVdV1\n2+67oqo+VVUPbLl2dVV9tKo+W1UfqaqrLu2tAAAAsFe7xmBVXZHkniSvSvLSJLdU1Uu2jbkxyQu7\n+0VJbkty77aneXOSM9uu3ZHkY9394iQfT/LWi3oHAMCRMpmsZHHxWCaTlXlPBeBI28vK4MuTPNrd\nj3X315N8IMnN28bcnOT+JOnuB5NcVVXPTZKqWkhyU5L3nuMx75/+/v4kP3FR7wAAOFJOn747GxvH\nc/r03fOeCsCRtpcYvDbJ+pbbG9NrO415fMuYf5rkHybpbY95TnefTZLufiLJc/Y4ZwDgCFtePpmF\nhbUsL5+c91QAjrQr9/PJq+rVSc5290NVtZSkdhi+PRafcurUqad+X1paytLS0oxmCAAcNG984xvy\nxje+Yd7TAJib1dXVrK6u7vvrVPd5G2xzQNX1SU5194np7TuSdHfftWXMvUk+0d0fnN5+JMkrs/m3\ngv9jkr9I8qwk35Pkl7v7dVX1cJKl7j5bVddMH/+D53j93m2OAAAAR1VVpbt3Wli7KHv5muhakh+o\nqudX1TOTvDbJA9vGPJDkdclT8fhkd5/t7rd19/O6+wXTx328u1+35TGvn/7+U0l+5dLeCgAAAHu1\n69dEu/sbVXV7ko9mMx7v6+6Hq+q2zbv7Pd394aq6qao+n+SrSW7dw2vfleQXq+oNSR5L8pMX/zYA\nAAC4ELt+TXTefE0UAAAY2Ty/JgoAAMARIwYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYB\nAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAG\nJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYB\nAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYB4DKaTFayuHgs\nk8nKvKcCwOCqu+c9hx1VVR/0OQLAXi0uHsvGxvEsLKxlff3MvKcDwCFQVenumvXzWhkEgMtoeflk\nFhbWsrx8ct5TAWBwVgYBAAAOMCuDAAAAzIwYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAY\nBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAA\nGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAY\nBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAC4QJPJShYXj2UyWZn3VAAuWnX3\nvOewo6rqgz5HAGAsi4vHsrFxPAsLa1lfPzPv6QBHXFWlu2vWz2tlEADgAi0vn8zCwlqWl0/OeyoA\nF83KIAAAwAFmZRAAAICZEYMAAAADEoMAAAADEoMAzIzt9gHg8LCBDAAzY7t9AJg9G8gAcODZbh8A\nDg8rgwAAAAeYlUEAAABmRgwCAAAMSAwCAAAMSAwCAAAMSAwCAAAMSAwCAAAMSAwCAAAMSAwCAAAM\nSAwCAAAMSAwCAAAMSAwCAAAMSAwCAAAMSAwCAAAMSAwCAAAMSAwCAAAMSAwCAAAMaE8xWFUnquqR\nqvpcVb3lPGPeXVWPVtVDVXXd9Np3VdWDVfXpqvpMVd25ZfydVbVRVZ+a/pyYzVsCAABgN1fuNqCq\nrkhyT5Ibkvx+krWq+pXufmTLmBuTvLC7X1RVr0hyb5Lru/vPq+pvdfefVdUzkvzHqvrV7v7k9KHv\n6O53zPxdAQAAsKO9rAy+PMmj3f1Yd389yQeS3LxtzM1J7k+S7n4wyVVV9dzp7T+bjvmubMZnb3lc\nXcLcAQAAuEh7icFrk6xvub0xvbbTmMe/NaaqrqiqTyd5Ism/6+61LeNun36t9L1VddUFzx4AAICL\nsu8byHT3N7v7ZUkWkryiqo5N7/qFJC/o7uuyGYq+LgoAAHCZ7Po3g9lc5XveltsL02vbxyzuNKa7\n/6SqPpHkRJIz3f2HW+6eJPnQ+SZw6tSpp35fWlrK0tLSHqYNAABw+KyurmZ1dXXfX6e6e+cBmxu/\nfDabG8h8Kcknk9zS3Q9vGXNTkjd196ur6vok7+zu66vqryT5end/paqeleQjSd7e3R+uqmu6+4np\n438myfHu/rvneP3ebY4AAABHVVWlu2e+38quK4Pd/Y2quj3JR7P5tdL7uvvhqrpt8+5+zzTubqqq\nzyf5apJbpw//q0neP92R9IokH+zuD0/v+/npERTfTPKFJLfN9J0BAABwXruuDM6blUEAAGBk+7Uy\nuO8byAAAAHDwiEEAAIABiUEAAIABiUEAAIABiUEAAIABiUEA4NtMJitZXDyWyWRl3lMBYB85WgIA\n+DaLi8eysXE8CwtrWV8/M+/pAAzP0RIAwGWxvHwyCwtrWV4+Oe+pALCPrAwCAAAcYFYGAQAAmBkx\nCAAAMCAxCAAAMCAxCAAAMCAxCAAAMCAxCLBHDuIGAI4SR0sA7JGDuAGAeXC0BMCcOYgbADhKrAwC\nAAAcYFYGAQAAmBkxCAAAMCAxCAAAMCAxCAAAMCAxCAAAMCAxCAAAMCAxCAAAMCAxCAAAMCAxCAAA\nMCAxCAAAMCAxCAAAMCAxCAAAMCAxCAAAMCAxCAAAMCAxCAAAMCAxCAAAMCAxCAAAMCAxCAAAMCAx\nCAAAMCAxCAAAMCAxCAAAMCAxCAAAMCAxCABwBE0mK1lcPJbJZGXeUwEOqOruec9hR1XVB32OAAAH\nzeLisWxsHM/CwlrW18/MezrAJaiqdHfN+nmtDAIAHEHLyyezsLCW5eWT854KcEBZGQQAADjArAwC\nAAAwM2IQAABgQGIQAABgQGIQgKHYbh8ANtlABoCh2G4fgMPGBjIAMAO22weATVYGAQAADjArgwAA\nAMyMGAQAABiQGAQAABiQGAQAABiQGAQAABiQGASe4jBuAIBxOFoCeIrDuAEADh5HSwD7zmHcAADj\nsDIIAABwgFkZBAAAYGaunPcEADh6nv3s5I//eN6z4DC7+urky1+e9ywAjjZfEwVg5qoS/3RzKXyG\nAJ7ma6IAAADMjBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAY\nkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgE\nAA6dyWQli4vHMpmszHsqAIdWdfe857CjquqDPkcAvl1V4p9uLsVun6HFxWPZ2DiehYW1rK+fuXwT\nA5iDqkp316yf18ogAHDoLC+fzMLCWpaXT857KgCHlpVBAGbOyiCXymcI4GlWBgEAAJgZMQgAADAg\nMQgAADAgMQgAADAgMQhwxDh/DQDYC7uJAhwxB+H8NTtBcql8hgCeZjdRAPbE+WsAwF5YGQRg5qzq\ncKl8hgCeZmUQAACAmRGDAAAAAxKDAAAAA9pTDFbViap6pKo+V1VvOc+Yd1fVo1X1UFVdN732XVX1\nYFV9uqo+U1V3bhl/dVV9tKo+W1UfqaqrZvOWAAAA2M2uMVhVVyS5J8mrkrw0yS1V9ZJtY25M8sLu\nflGS25LcmyTd/edJ/lZ3vyzJdUlurKqXTx92R5KPdfeLk3w8yVtn85YAAADYzV5WBl+e5NHufqy7\nv57kA0lu3jbm5iT3J0l3P5jkqqp67vT2n03HfFeSK5P0lse8f/r7+5P8xMW+CQAAAC7MXmLw2iTr\nW25vTK/tNObxb42pqiuq6tNJnkjy77p7bTrmOd19Nkm6+4kkz7nw6cOFmUxWsrh4LJPJyrynAozo\nK19JXvOa5NprN//3K1+Z94wAGNi+byDT3d+cfk10IckrqurY+Ybu91zg9Om7s7FxPKdP3z3vqQAj\nuvXW5N/+2+T3f3/zf2+9dd4zAmBgV+5hzONJnrfl9sL02vYxizuN6e4/qapPJDmR5EySs1X13O4+\nW1XXJPmD803g1KlTT/2+tLSUpaWlPUwbvtPy8smcPn13lpdPznsqwIgefDD52tc2f//a1zZvA8A2\nq6urWV1d3ffXqe6dF+Sq6hlJPpvkhiRfSvLJJLd098NbxtyU5E3d/eqquj7JO7v7+qr6K0m+3t1f\nqapnJflIkrd394er6q4kX+7uu6Y7lF7d3Xec4/V7tzkCcLBUJf7pPofXvGZzRfBrX0ue+czk1a9O\nfvmX5z2rA8lnCOBpVZXurpk/715Cq6pOJHlXNr9Wel93v72qbkvS3f2e6Zh7srnq99Ukt3b3p6rq\nb2Rzc5grpj8f7O5/PB3/7CS/mM0VxceS/GR3P3mO1xaDAIeM/5A/j698ZfOroQ8+mLziFcn73pdc\n5WSlc/EZAnjaXGNwnsQgwOHjP+S5VD5DAE/brxjc9w1kAAAAOHjEIAAAwIDEIAAAwIDEIAAAwIDE\nIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAczOZrGRx8Vgmk5V5\nTwWGU9097znsqKr6oM8RgG9Xlfinm0vhMzSOxcVj2dg4noWFtayvn5n3dOBAqqp0d836ea0MAgAw\nN8vLJ7OwsJbl5ZPzngoMx8ogADNnVYdL5TME8DQrgwAAAMyMGAQAABiQGAQAABiQGASAA8ZW+wBc\nDjaQAWDmbP5xaWy17zMEsJUNZABgELbaB+BysDIIwMxZ1eFS+QwBPM3KIAAAADMjBgEAAAYkBgEA\nAAYkBgEAAAYkBgEAAAYkBoFDxWHcAACz4WgJ4FBxGPfh4FgALpXPEMDTHC0BEIdxAwDMipVBAGbO\nqg6XymcI4GlWBgEAAJgZMQgAADAgMQgAADAgMQgAADCgK+c9AQCOnquv3twABAA4uOwmCgAAcIDZ\nTRQAAICZEYMAAAADEoMAAAADEoMAAAADEoMAAAADEoMAAAADEoMAAAADEoMAAAADEoMAAAADEoMA\nAAADEoMAAAADEoMAAAADEoMAAAADEoMAAAADEoMAAAADEoMAAPtgMlnJ4uKxTCYr854KwDlVd897\nDjuqqj7ocwQA2G5x8Vg2No5nYWEt6+tn5j0d4BCrqnR3zfp5rQwCAOyD5eWTWVhYy/LyyXlPBeCc\nrAwCAAAcYFYGAQAAmBkxCAAAMCAxCAAAMCAxCAAAMCAxCMBl5/w1AJg/u4kCcNk5fw0A9s5uogAc\nGc5fA4D5szIIAABwgFkZBAAAYGbEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIBwwDuMGAOBycLQE\nHDAO4wYAYCtHS8AgHMYNAMDlYGUQAADgALMyCAAAwMyIQQAAgAGJQQAAgAGJQQAAgAGJQQAAgAGJ\nQQAAgAGJQQAAgAGJQQAAgAGJQQAAgAGJQQAAgAGJQQAAgAGJQQAAgAGJQQAAgAGJQQAAgAGJQQAA\n2MFkspLFxWOZTFbmPRWYqeruec9hR1XVB32OAAAcXYuLx7KxcTwLC2tZXz8z7+kwoKpKd9esn9fK\nIAAA7GB5+WQWFtayvHxy3lOBmbIyCAAAcIBZGQQAAGBmxCAAAMCAxCAAAMCAxCAAcMFstQ9w+NlA\nBgC4YLbaB7h85rqBTFWdqKpHqupzVfWW84x5d1U9WlUPVdV102sLVfXxqvq9qvpMVf30lvF3VtVG\nVX1q+nNiNm8JANhvttoHOPx2XRmsqiuSfC7JDUl+P8laktd29yNbxtyY5PbufnVVvSLJu7r7+qq6\nJsk13f1QVf3lJP85yc3d/UhV3Znk/+3ud+zy+lYGAQCAYc1zZfDlSR7t7se6++tJPpDk5m1jbk5y\nf5J094NJrqqq53b3E9390PT6nyZ5OMm1Wx438zcEAADA7vYSg9cmWd9yeyPfHnTnGvP49jFV9f1J\nrkvy4JbLt0+/Vvreqrpqj3MGAADgEl15OV5k+hXRX0ry5ukKYZL8QpLT3d1V9Y+SvCPJ/3Sux586\ndeqp35eWlrK0tLSv8wUAAJiX1dXVrK6u7vvr7OVvBq9Pcqq7T0xv35Gku/uuLWPuTfKJ7v7g9PYj\nSV7Z3Wer6sok/ybJr3b3u87zGs9P8qHu/qFz3OdvBgEAgGHN828G15L8QFU9v6qemeS1SR7YNuaB\nJK9LnorHJ7v77PS+lSRntofgdHOZb3lNkt+9iPkDHDjOXwMADoM9nTM4PfbhXdmMx/u6++1VdVs2\nVwjfMx1zT5ITSb6a5PXd/emq+tEk/yHJZ5L09Odt3f1rVXV/Nv+G8JtJvpDkti0BufW1rQwCh4rz\n1wCAWdqvlUGHzgPM2GSyktOn787y8sm88Y1vmPd0AIBDTgwCAAAMaJ5/MwgAAMARIwYBAAAGJAYB\nAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAG\nJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYB\nAAAGJAYBAAY1maxkcfFYJpOVeU8FmIPq7nnPYUdV1Qd9jgAAh9Hi4rFsbBzPwsJa1tfPzHs6wHlU\nVbq7Zv28VgYBAAa1vHwyCwtrWV4+Oe+pAHNgZRAAAOAAszIIAADAzIhBAACAAYlBAACAAYlBAACA\nAYlBADgH568BcNTZTRQAzsH5awAcFHYTBYDLyPlrABx1VgYBAAAOMCuDAAAAzIwYBAAAGJAYBAAA\nGJAYBAAAGJAYBAAAGJAYBC6Ig7gBAI4GR0sAF8RB3AAAl5ejJYADwUHcAABHg5VBAACAA8zKIAAA\nADMjBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYk\nBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEA4JCbTFayuHgsk8nKvKfCIVLdPe857Kiq+qDPEQAA\n5mlx8Vg2No5nYWEt6+tn5j0dZqyq0t016+e1MggAAIfc8vLJLCysZXn55LynwiFiZRAAAOAAszII\nAADAzIhBAACAAYlBAOBIsrsiwM78zSAAcCTZXRE4KvzNIADABbC7IsDOrAwCAAAcYFYGAQAAmBkx\nCAAAMCAxCAAAMCAxCAAAMCAxCDAg568BAHYTBRiQ89cA4PCwmygAM+P8NQDAyiAAAMABZmUQAACA\nmRGDAAAAAxKDAAAAAxKDAAAAAxKDAAAAAxKDAAAAAxKDAAAAAxKDAAAAAxKDAAAAAxKDAAAAAxKD\nAAAAAxKDAAAAAxKDAAAAAxKDAAAAAxKDAAAAAxKDAAAAAxKDAAAAAxKDAAAAAxKDAAAAA9pTDFbV\niap6pKo+V1VvOc+Yd1fVo1X1UFVdN722UFUfr6rfq6rPVNVPbxl/dVV9tKo+W1UfqaqrZvOWAAAA\n2M2uMVhVVyS5J8mrkrw0yS1V9ZJtY25M8sLuflGS25LcO73rL5L8bHe/NMl/neRNWx57R5KPdfeL\nk3w8yVtn8H4AAADYg72sDL48yaPd/Vh3fz3JB5LcvG3MzUnuT5LufjDJVVX13O5+orsfml7/0yQP\nJ7l2y2PeP/39/Ul+4pLeCQAAAHu2lxi8Nsn6ltsbeTrozjfm8e1jqur7k1yX5Leml57T3WeTpLuf\nSPKcvU4aAACAS3NZNpCpqr+c5JeSvLm7v3qeYX055gIAwOExmaxkcfFYJpOVeU8Fjpwr9zDm8STP\n23J7YXpt+5jFc42pqiuzGYL/srt/ZcuYs9Ovkp6tqmuS/MH5JnDq1Kmnfl9aWsrS0tIepg0AwGF3\n+vTd2dg4ntOn784b3/iGeU8HLovV1dWsrq7u++tU984LclX1jCSfTXJDki8l+WSSW7r74S1jbkry\npu5+dVVdn+Sd3X399L77k/xRd//stue9K8mXu/uu6Q6lV3f3Hed4/d5tjgAAHE2TyUpOn747y8sn\nxSDDqqp0d838efcSWlV1Ism7svm10vu6++1VdVuS7u73TMfck+REkq8meX13f7qqfjTJf0jymWx+\nDbSTvK27f62qnp3kF7O5ovhYkp/s7ifP8dpiEAAAGNZcY3CexCAAADCy/YrBy7KBDAAAAAeLGAQA\nABiQGASAQ8h2+wBcKn8zCACH0OLisWxsHM/CwlrW18/MezoA7CN/MwgAPGV5+WQWFtayvHxy3lMB\n4JCyMgioUikfAAAU00lEQVQAAHCAWRkEAABgZsQgAADAgMQgAADAgMQgAADAgMQgAADAgMQgcOQ4\njBsAYHeOlgCOHIdxAwBHiaMlAPbIYdwAALuzMggAAHCAWRkEAABgZsQgAADAgMQgAADAgMQgAADA\ngMQgAADAgMQgAADAgMQgAADAgMQgAADAgMQgAADAgMQgAADAgMQgAADAgMQgAADAgMQgAADAgMQg\nAADAgMQgAADAgMQgAADAgMQgAADAgMQgAADAgMQgAADAgMQgAADAgMQgAADAgMQgAMCcTCYrWVw8\nlslkZd5TAQZU3T3vOeyoqvqgzxEA4GIsLh7LxsbxLCysZX39zLynAxxQVZXurlk/r5VBAIA5WV4+\nmYWFtSwvn5z3VIABWRkEAAA4wKwMAgAAMDNiEAAAYEBiEAAAYEBiEAAAYEBiEIADyflrALC/7CYK\nwIHk/DUA2GQ3UQCG4vw1ANhfVgYBAAAOMCuDAAAAzIwYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAY\nBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAA\nGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAY\nBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAAAGJAYBAD2xWSyksXFY5lMVuY9FQDOobp73nPYUVX1\nQZ8jAPCdFhePZWPjeBYW1rK+fmbe0wE4tKoq3V2zfl4rgwDAvlhePpmFhbUsL5+c91QAOAcrgwAA\nAAeYlUEAAABmRgwCAAAMSAwCAAAMSAwCAAAMSAwCzIHz1wCAebObKMAcOH8NANgru4kCHCHOXwMA\n5s3KIAAAwAFmZRAAAICZEYMAAAADEoMAAAADEoMAAAADEoMAAAADEoMAAAADEoMAAAAD2lMMVtWJ\nqnqkqj5XVW85z5h3V9WjVfVQVb1sy/X7qupsVf3OtvF3VtVGVX1q+nPi0t4KAAAAe7VrDFbVFUnu\nSfKqJC9NcktVvWTbmBuTvLC7X5TktiT/bMvd75s+9lze0d0/PP35tYt5AwAAAFy4vawMvjzJo939\nWHd/PckHkty8bczNSe5Pku5+MMlVVfXc6e1fT/LH53nuuqhZAwAAcEn2EoPXJlnfcntjem2nMY+f\nY8y53D79Wul7q+qqPYwHAABgBua5gcwvJHlBd1+X5Ikk75jjXAAAAIZy5R7GPJ7keVtuL0yvbR+z\nuMuYb9Pdf7jl5iTJh8439tSpU0/9vrS0lKWlpZ2eGgAA4NBaXV3N6urqvr9OdffOA6qekeSzSW5I\n8qUkn0xyS3c/vGXMTUne1N2vrqrrk7yzu6/fcv/3J/lQd/+NLdeu6e4npr//TJLj3f13z/H6vdsc\nAQAAjqqqSnfPfL+VXVcGu/sbVXV7ko9m82ul93X3w1V12+bd/Z7u/nBV3VRVn0/y1SS3bpn4v06y\nlOT7quqLSe7s7vcl+fmqui7JN5N8IZu7kAIAAHAZ7LoyOG9WBgEAgJHt18rgPDeQAQAAYE7EIAAA\nwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDEIAAAwIDE\nIAAAwIDEIAAAwIDEIABchMlkJYuLxzKZrMx7KgBwUaq75z2HHVVVH/Q5AjCexcVj2dg4noWFtayv\nn5n3dAA4wqoq3V2zfl4rgwBwEZaXT2ZhYS3LyyfnPRUAuChWBgEAAA4wK4MAAADMjBgEAAAYkBgE\nAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAY\nkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgE\nAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAY\nkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgEAAAYkBgE4DtMJitZXDyWyWRl\n3lMBAPZJdfe857CjquqDPkeAo2Zx8Vg2No5nYWEt6+tn5j0dABhaVaW7a9bPa2UQgO+wvHwyCwtr\nWV4+Oe+pAAD7xMogAADAAWZlEAAAgJkRgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAA\nAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMS\ngwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAA\nAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMSgwAAAAMS\ngwAAAAMSgwAAAAPaUwxW1YmqeqSqPldVbznPmHdX1aNV9VBVvWzL9fuq6mxV/c628VdX1Uer6rNV\n9ZGquurS3goAAAB7tWsMVtUVSe5J8qokL01yS1W9ZNuYG5O8sLtflOS2JP9sy93vmz52uzuSfKy7\nX5zk40neelHvAAAAgAu2l5XBlyd5tLsf6+6vJ/lAkpu3jbk5yf1J0t0PJrmqqp47vf3rSf74HM97\nc5L3T39/f5KfuPDpAwAAcDH2EoPXJlnfcntjem2nMY+fY8x2z+nus0nS3U8kec4e5gIAAMAMXDnv\nCWzR57vj1KlTT/2+tLSUpaWlyzAdAACAy291dTWrq6v7/jrVfd4G2xxQdX2SU919Ynr7jiTd3Xdt\nGXNvkk909wentx9J8spvrfxV1fOTfKi7f2jLYx5OstTdZ6vqmunjf/Acr9+7zREAAOCoqqp0d836\neffyNdG1JD9QVc+vqmcmeW2SB7aNeSDJ65Kn4vHJb4XgVE1/tj/m9dPffyrJr1zY1AEAALhYu8Zg\nd38jye1JPprk95J8oLsfrqrbqurvTcd8OMn/U1WfT/LPk/z9bz2+qv51kt9I8ter6otVdev0rruS\n/FhVfTbJDUnePsP3BQAA/3979x/re13XAfz5CmwxWSaYXCYiIoFpILoktjKzjXVlJcRcibOyH8ss\nyDVrkrIg4g9j7TbLmIrVxA2lQQsoQGRCrhwOvUBgF7xakhiQtmyjXxC8+uN8zvbt7Nx7z73nfM/3\n3Pt+PLbPzvl+Pu/P5/P+nL322fd5vu/v5w3sxT6HiS6aYaIAAMDIFjlMFAAAgEOMMAgAADAgYRAA\nAGBAwiAAAMCAhEEAAIABCYMAAAADEgYBAAAGJAwCAAAMSBgEAAAYkDAIAAAwIGEQAABgQMIgAADA\ngIRBAACAAQmDAAAAAxIGAQAABiQMAgAADEgYBAAAGJAwCAAAMCBhEAAAYEDCIAAAwICEQQAAgAEJ\ngwAAAAMSBgEAAAYkDAIAAAxIGAQAABiQMAgAADAgYRAAAGBAwiAAAMCAhEEAAIABCYMAAAADEgYB\nAAAGJAwCAAAMSBgEAAAYkDAIAAAwIGEQAABgQMIgAADAgIRBAACAAQmDAAAAAxIGAQAABiQMAgAA\nDEgYBAAAGJAwCAAAMCBhEAAAYEDCIAAAwICEQQAAgAEJgwAAAAMSBgEAAAYkDAIAAAxIGAQAABiQ\nMAgAADAgYRAAAGBAwiAAAMCAhEEAAIABCYMAAAADEgYBAAAGJAwCAAAMSBgEAAAYkDAIAAAwIGEQ\nAABgQMIgAADAgIRBAACAAQmDAAAAAxIGAQAABiQMAgAADEgYBAAAGJAwCAAAMCBhEAAAYEDCIAAA\nwICEQQAAgAEJgwAAAAMSBgEAAAYkDAIAAAxIGAQAABiQMAgAADAgYRAAAGBAwiAAAMCAhEEAAIAB\nCYMAAAADEgYBAAAGJAwCAAAMSBgEAAAYkDAIAAAwIGEQAABgQMIgAADAgIRBAACAAQmDAAAAAxIG\nAQAABiQMAgAADEgYBAAAGJAwCAAAMCBhEAAAYEDCIAAAwICEQQAAgAEJgwAAAAMSBgEAAAa0pjBY\nVdur6sGq+mJVvWsPbf6gqnZX1b1Vdfq+9q2qS6rqkaraOS3b1385AAAArMU+w2BVfUuS9yf5kSQv\nT3J+Vb10RZvXJ3lJd39Xkrcl+cAa993R3a+alls34oIY15133rnoLnAQUS+slVphf6gX1kqtsBWs\n5ZPBM5Ls7u6Hu/upJB9Pcs6KNuckuTpJuvuzSZ5TVcesYd9a7wXAMjdV9od6Ya3UCvtDvbBWaoWt\nYC1h8AVJvjrz+pFp3Vra7GvfC6ZhpR+uquesudcAAACsy7weILOWT/yuTHJid5+e5LEkO+bUFwAA\nAFao7t57g6ozk1za3dun1xcl6e7+3Zk2H0hyR3dfO71+MMlrk7x4X/tO61+U5KbuPm2V8++9gwAA\nAIe47t7wr9gdvoY2dyc5aQpsjyZ5U5LzV7S5McmvJLl2Co/f7O7Hq+obe9q3qrZ192PT/ucleWC1\nk8/jogEAAEa3zzDY3U9X1QVJbsvSsNI/7u5dVfW2pc39oe6+uarOrqovJfmPJD+7t32nQ18xTUHx\nTJKvZOkppAAAAGyCfQ4TBQAA4NAzrwfI7Jeqem5V3VZVD1XVJ/b0ZNG9TGD/8ZnJ6/+xqnZuXu/Z\nTOutlWnbhVW1q6rur6r3bk7PWYQNuLdcUlWPzNxftm9e79lMG3Fvmba/s6qeqaqj5t9rFmUD7i2X\nVdV9VXVPVd1aVds2r/dspg2olSum9yz3VtX1VfXtm9d7NtsG1Msbq+qBqnq6ql61lnNuiTCY5KIk\nt3f3KUk+leQ3VzbY2wT23f2m5cnrk1yf5M83redstnXVSlX9UJIfS3Jqd5+a5Pc2qd8sxrrqZbJj\n+f7S3bduRqdZiHXXSlUdl+SsJA9vSo9ZpPXWyxXd/YrufmWSv0pyyeZ0mwVYb63cluTl09P3d6+2\nP4eU9dbL/Ul+PMlfr/WEWyUMnpPkI9PvH0ly7ipt9jWB/bKfSPKxufSSrWC9tfL2JO/t7v9Nku7+\nxpz7y2JtxL3FQ6zGsBG18vtJfmOuvWSrWFe9dPcTM+2enaXnJ3BoWm+t3N7dy/VxV5Lj5txfFmu9\n9fJQd+/Ofrx32Sph8Pnd/XiSTE8Yff4qbfY1gX2q6jVJHuvuL8+royzcemvl5CQ/WFV3VdUdVfW9\nc+0ti7YR95YLpuE5H97TcA0OCeuqlap6Q5Kvdvf98+4oW8K67y1VdXlV/VOSNyf5rTn2lcXakPe4\nk59LcsuG95CtZCPrZU3WMrXEhqiqTyY5ZnZVkk5y8SrND/SpNufHp4IHvTnXyuFJntvdZ1bVq5P8\nWZITD6ijbAlzrpcrk1zW3V1VlyfZkeTnD6ijLNy8aqWqjkjy7iwNEZ09Ngexeb9v6e6Lk1w8fd/n\nwiSXHkA32QI24z1uVb0nyVPdfc2B7M/WsUmZaM02LQx291l72lZVj1fVMdPchNuS/Msqzb6W5PiZ\n18dN65aPcViW5itc05cl2brmXCuPZPpOaXffPT3o4eju/tcN6j6bbJ710t1fn1l/VZKbNqDLLMgc\na+UlSU5Icl9V1bT+81V1RnevdhwOAvN+3zLjmiQ3Rxg8aG3Ce9y3Jjk7yQ9vTI9ZpE28t6zJVhkm\nemOSt06//0ySG1Zpc3emCeyr6luzNIH9jTPbz0qyq7v/eZ4dZeHWWyt/kelmWlUnJ3mWIHhIW1e9\nrHjC33lJHphfV1mwA66V7n6gu7d194nd/eIs/dPplYLgIW2995aTZtqdm2TXKvtzaFhvrWzP0neR\n39Dd/zP/7rJgG5GJlq1thEp3L3xJclSS25M8lKWnJn3HtP7YJH8502771GZ3kotWHONPk/zioq/F\nsrVrJcmzknw0S09b+lyS1y76mixbul6uTvJ3Se7N0j8Sjln0NVm2Zq2sONY/JDlq0ddk2br1kuS6\nmXvLDUmOXfQ1WbZsrezO0hOKd07LlYu+JsuWrpdzs/R9wv9K8miSW/Z1TpPOAwAADGirDBMFAABg\nEwmDAAAAAxIGAQAABiQMAgAADEgYBAAADmlV9caqeqCqnq6qVeclr6rjqupTVfWFqrq/qn51Zttl\nVXVfVd1TVbcuTz9VVa+e1i0v507rj5xe75x+fr2qduyjj2+eznFfVf1NVZ26kX+DVc/paaIAAMCh\nrKpOSfJMkg8m+fXu3rlKm21JtnX3vVV1ZJLPJzmnux+sqiO7+4mp3YVJXtbdb6+qb0vyZHc/M+1/\nX5ami3lmxbE/l+Qd3f23e+njmVmaN/3fpzkmL+3uMzfkD7AHh8/z4AAAAIvW3Q8lSVXtcTL27n4s\nyWPT709U1a4kL0jy4HIQnDw7S8Ey3f3fM+uPWF4/q6pOTvKdy0Gwqp6X5ANJXjg1+bXu/kx33zWz\n213TuedKGAQAAJhRVSckOT3JZ2fWXZ7kp5N8M8nrZtafkeRPkhyf5KdWfiqY5CeTXDvz+n1JdnT3\nZ6rqhUk+keRlK/b5hSS3bMS17I1hogAAwEGvqj6Z5JjZVUk6yXu6+6apzR1J3rnaMNGZ4xyZ5M4k\nv9PdN6yy/V1JjujuS1esPyXJ1Ule091Pzqz/QpK3dPc90+vHk3xt6l+SHJ3kpd39n9P21yV5f5If\n6O5/W/Mf4AD4ZBAAADjodfdZ6z1GVR2e5LokH10tCE6uSXJzkktXnP+hqnoiyfck2Tkd77Qkhy0H\nweXTJPm+7n5qlfOfluRDSbbPOwgmniYKAACMZY/fG8zScM+/7+73/b8dqk6aeXlukl3T+hOq6rDp\n9xclOSXJV2banp/kYyvOcVuSd8wc+xXTz+OTXJ+loaZf3o/rOWCGiQIAAIe0acqHP0zyvCx95+/e\n7n59VR2b5Kru/tGq+v4kn05yf5aGl3aSd3f3rVV1XZKTs/SAmIeT/FJ3P1pVb0lyUZInp22/vTwk\ndTrvl5Kc3d1fnFl3dJI/SvLdSQ5L8unu/uWquirJedPxK8lT3X3GHP8swiAAAMCIDBMFAAAYkDAI\nAAAwIGEQAABgQMIgAADAgIRBAACAAQmDAAAAAxIGAQAABiQMAgAADOj/AI/0zGz5KjWaAAAAAElF\nTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "lon, lat = -123.4085, 46.91\n", "tol_lon, tol_lat = 0.0052, 0.00189\n", "fig, ax = plt.subplots(1, 1,figsize=(15,15))\n", "plt.scatter(X[6:12,6:12],Y[6:12,6:12], s =3 )\n", "#plt.scatter([lon], [lat], c = \"w\" , s =80)\n", "#plt.scatter([X[tuple(reversed(closest_point))]], [Y[tuple(reversed(closest_point))]], c = \"c\" , s =40)\n", "#plt.scatter([X[109,143]], [Y[109,143]], c = \"c\" , s =40)\n", "plt.plot([lon + tol_lon, lon + tol_lon, lon - tol_lon, lon - tol_lon,lon + tol_lon],\n", " [lat - tol_lat, lat + tol_lat, lat + tol_lat, lat - tol_lat,lat - tol_lat])\n", "plt.scatter([lon],[lat], s= 14, color = 'r')" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4AAAANmCAYAAACrKeXIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3W+obel9H/bfb+Y6oNSS7FdFsXDGpQQbidaZxs4fadBR\nrRAxUsdgGv1JjN41kDjJYAkhN2/OudCauHEQLsmbWImLSiqBbKURgrZhCEdYIcRJJLXRVI6Doolk\nOxKYxlZNSTuSn76Yc2fOPXfvc9bee/15/nw+l2HuOWfvtZ+9zrlnr+/+PutZWUoJAAAA+vfY1gMA\nAABgHQIgAADAIARAAACAQQiAAAAAgxAAAQAABiEAAgAADGJyAMzMxzLz85n5qauPfzAz//HV534l\nM//Ijvu8PjP/YWY+n5n/IjP/8pyDBwAAYLpDGsBnI+L5ax//TEScl1L+cEScR8Rf23Gfb0XE+0sp\nb4iIPx4RP5GZ33/sYAEAADjepACYma+PiKcj4iPXPv17EfHaq79/V0T8xs37lVK+Xkr5wtXffzci\nvhQR33PKgAEAADjOvYm3+3BEfDBeCXwRET8ZEf9bZv71iMiI+BO3bSAzn4iIH4yIf3LwKAEAADjZ\nnQ1gZr4jIr5x1eTltS/9+Yh4tpTyvfFSGPw7t2zjOyPiF69u/7unDRkAAIBjZCnl9htk/nRE/Hi8\ndD7fqyLi1RHx9yLinaWU7752u98ppbx2x/3vRcSnI+J/KaX83C2Pc/tAAAAAOldKybtvdbw7A+BD\nN858S0R8oJTyTGY+HxF/oZTymcz8kYj4q6WUH9pxn49GxG+VUt5/x7bLIWOBtVxcXMTFxcXWw4BH\n+NmkZn4+qZWfTWqWmYsHwKnnAO7y5yLi5zLz8Yj491cfR2a+LiJ+vpTyzsx8U0T82Yj4F5n5+Ygo\nEfFXSin/64njBgAA4EAHBcBSymci4jNXf/9HEfHItf9KKf82It557TaPnz5MAAAATnXIdQBhSGdn\nZ1sPAXbys0nN/HxSKz+bjO6gcwCX5BxAAABgZGucA6gBBAAAGIQACAAAMAgBEAAAYBACIAAAwCAE\nQAAAgEEIgAAAAIMQAAEAAAYhAAIAAAxCAAQAABiEAAgAADAIARAAAGAQAiAAAMAgBEAAAIBBCIAA\nAACDEAABAAAGIQACAAAMQgAEAAAYhAAIAAAwCAEQAABgEAIgAADAIARAAACAQQiAAAAAgxAAAQAA\nBiEAAgAADEIABAAAGIQACAAAMAgBEAAAYBACIAAAwCAEQAAAgEEIgAAAAIMQAAEAAAYhAAIAAAxC\nAAQAABiEAAgAADAIARAAAGAQAiAAAMAgBEAAAIBBCIAAAACDEAABAAAGIQACAAAMQgAEAAAYhAAI\nAAAwCAEQAABgEAIgAADAIARAAACAQQiAAAAAgxAAAQAABiEAAgAADEIABAAAGIQACAAAMAgBEAAA\nYBACIAAAwCAEQAAAgEEIgAAAAIMQAAEAAAYhAAIAAAxCAAQAABiEAAgAADAIARAAAGAQAiAAAMAg\nBEAAAIBBCIAAAACDEAABAAAGIQACAAAMQgAEAAAYhAAIAAAwCAEQAABgEAIgAADAIARAAACAQQiA\nAAAAgxAAAQAABiEAAgAADEIABAAAGIQACAAAMAgBEAAAYBACIAAAwCAEQAAAgEEIgAAAAIMQAAEA\nAAYhAAIAAAxCAAQAABiEAAgAADAIARAAAGAQAiAAAMAgBEAAAIBBCIAAAACDEAABAAAGIQACAAAM\nQgAEAAAYhAAIAAAwCAEQAABgEAIgAADAIARAAACAQQiAAAAAgxAAAQAABiEAAgAADEIABAAAGIQA\nCAAAMAgBEAAAYBACIAAw2f3MrYcAwAmylLL1GCIiIjNLLWMBAB52M/ide80GmF1mRill0Xfa7i25\ncQCgbRo/gL6YAgoA7HRb+BMMAdqkAQQAHiLcAfRLAwgAvOyQ8CcoArRHAwgACHMAg9AAAsDgTgl/\ngiNAWzSAADAo4Q1gPBpAABjQnOFPkARox+QAmJmPZebnM/NTVx//YGb+46vP/Upm/pE993t7Zv5q\nZv5aZn5oroEDAIe7nymwAQzskAbw2Yh4/trHPxMR56WUPxwR5xHx127eITMfi4i/ERF/KiLeEBHv\nzczvP364AMCxlgx+QiVAGyYFwMx8fUQ8HREfufbp34uI1179/bsi4jd23PWHI+JflVL+TSnlxYj4\neET86PHDBQAOpfUD4IGpDeCHI+KDEVGufe4nI+JnM/OrEfHfRcR/veN+3xMRX7v28a9ffQ4AWMGa\nwU/IBKjfnauAZuY7IuIbpZQvZObZtS/9+Yh4tpTyP2fmfxkRfyci/uQpg7m4uHj572dnZ3F2drb3\ntgDAfsIYQP0uLy/j8vJy1cfMUsrtN8j86Yj48Yj4VkS8KiJeHRF/LyLeWUr57mu3+51Symtv3PeP\nRcRFKeXtVx//VESUUsrP7HicctdYAIC7bRn+zr2WAxwtM6OUsugv8TsD4EM3znxLRHyglPJMZj4f\nEX+hlPKZzPyRiPirpZQfunH7xyPiX0bEj0TEv42IX4mI95ZSvrRj2wIgAJygltZPCAQ4zhoB8JQL\nwf+5iPi5q5D3768+jsx8XUT8fCnlnaWUb2fmX4yIfxAvnW/4t3eFPwDgNLWEPwDqdlADuCQNIAAc\nrtbgpwUEONwaDeAh1wEEACpSa/gDoF4aQABoTCvBTwsIcBgNIADwkFbCHwB10gACQANaDX5aQIDp\nNIAAQLPhD4D6aAABoFK9BD8tIMA0GkAAGFQv4Q+AumgAAaAivQY/LSDA3TSAADCQXsMfAPXQAALA\nxkYJflpAgNtpAAGgc6OEPwDqcG/rAQDAiAQ/ALagAQSAlY0a/kZ93gA10QACwEoEIAC2ZhEYAFiY\n4Pcwi8EA7GYRGABonPAHQE00gACwAMHvdlpAgEdpAAGgQcIfALXSAALATAS/w2gBAR6mAQSARgh/\nALRAAwgAJxD8TqMFBHiFBhAAKib8AdAaDSAAHEjwm5cWEOAlGkAAqIzwB0DLNIAAMIHgtywtIIAG\nEACqIPwB0It7Ww8AAGol+AHQGw0gAOwg/K3L/gZYhwYQAK4RRADomQYQAK4If9uy/wGWpwEEYHiC\nBwCj0AACMDThry6+HwDL0gACMCRBA4ARaQABGI7wVzffH4DlaAABGIZgAcDoNIAADEH4a4vvF8Ay\nNIAAdE2QAIBXaAAB6Jbw1zbfP4D5aQAB6I7gAAC7aQAB6Irw1xffT4B5aQAB6IKgAAB30wAC0Dzh\nDwCmyVLK1mOIiIjMLLWMBYA2CH7jOHeMAAwgM6OUsuiLmwYQgCYJfwBwOA0gAE0R/MalBQR6pwEE\ngGuEPwA4jQYQgOoJfjygBQR6pgEEYHjCHwDMRwMIQJUEP/bRAgK90gACMCThDwCWoQEEoBqCH1Np\nAYEeaQABGIbwBwDL0wACsCnBj2NpAYHeaAAB6JrwBwDr0gACsDrBj7loAYGeaAAB6I7wBwDbubf1\nAAAYg+AHANvTAAKwOOGPpfjZAjiMBhCAxTg4B4C6aAABWITwx1r8rAFMpwEEYFYOxgGgXhpAAGYj\n/LEVP3sA02gAATiZg28AaIMGEICTCH/Uws8iwN00gAAcxcE2ALRHAwjAwYQ/auVnE+B2GkAAJnNw\nDQBt0wACMInwRyv8rALspwEE4FYOpgGgHxpAAPYS/miVn12A3TSAADzCwTMA9EkDCMBDhD8A6FeW\nUrYeQ0REZGapZSwAIxL86NG5YwugIZkZpZRFX5A1gAAIfwAwCA0gwMAEP0agBQRaoQEEYDHCHwCM\nRwMIMBjBjxFpAYEWaAABmJXwBwBjEwABBiH8MTI//wAvcSF4gM458AUAHtAAAnRM+INX+PcAoAEE\n6JIDXQBgFw0gQGeEP9jPvw9gdBpAgE44sAUA7qIBBOiA8AfT+fcCjEwDCNAwB7IAwCE0gACNEv4A\ngENlKWXrMURERGaWWsYCUDPBD+Zx7rgDqExmRill0Rd6DSBAQ4Q/AOAUGkCABgh+sAwtIFATDSAA\nwh8AMBsNIEClBD9YhxYQqIUGEGBQwh8AsAQNIEBFBD/YhhYQqIEGEGAgwh8AsDQNIMDGBD+ogxYQ\n2JoGEKBzwh8AsCYNIMAGBD+okxYQ2JIGEKBDwh8AsBUNIMBKBD9ogxYQ2IoGEKATwh8AUIN7Ww8A\noGeCH7TjIi4iIuJ822EALEoDCLAQ4Q/a8SD8RURk3t9uIAALEwABZnY/U/iDxlwPgAA9EwABZiT4\nQbu0gMAInAMIMBPhD9qlAQRG4TIQACcS/KBddwW/UiwJA6zHZSAAKif8AQAtEQABjmChFxiDcwGB\n3giAAAcS/KAfzv0DRjM5AGbmY5n5+cz81NXHH8/Mz13995XM/Nye+/1kZn4xM/+PzPy7mfn75ho8\nwJq0fjAmLSDQk0MawGcj4vkHH5RS3lNKebKU8mRE/FJEfPLmHTLzD0TEX4qIJ0sp/0m8tOroe04b\nMsD6BD/olxYQGMmkAJiZr4+IpyPiI3tu8q6I+Nierz0eEf9BZt6LiN8fEb956CABtqL1g/5NCYB+\nDwC9mNoAfjgiPhgRj1ynITOfioivl1K+fPNrpZTfjIi/HhFfjYjfiIjfLqU8d/xwAdbjgA/6NyX8\naQiBntwZADPzHRHxjVLKFyIir/677r2xp/3LzO+KiB+NiD8YEX8gIr4zM//MSSMGWJjWD8ZxV7i7\n/nW/F4Ae3JtwmzdFxDOZ+XREvCoiXp2ZHy2lvC8zH4+IH4uIJ/fc920R8a9LKf9XRERmfjIi/kRE\n/E+7bnxxcfHy38/OzuLs7Gzi0wCYhwM8AGAtl5eXcXl5uepjZimPzOrcf+PMt0TEB0opz1x9/PaI\n+FAp5a17bv/DEfG3I+KHIuL/jYhfiIh/Wkr5mztuWw4ZC8CcBD8Y221N4PWvnTtWARaUmVFKWfSg\n5NTrAL47bkz/zMzXZeanIyJKKb8SEb8YEZ+PiP89Xpo++rdOfEyAWQl/wK4AeHH15zq/L4DWHdQA\nLkkDCKzNgRwQsb/92/d5LSCwlBYaQIAmCX/AXfYFQL8/gJZNWQQGoBsO3ICpXP4B6JEGEBiG8Afs\nsu/8v9v4fQK0SgAEhuBgDTiE9g/olUVggK4JfsBUxzSBFoQB5rTGIjDOAQS6JfwBU0wNfrsuCwHQ\nGgEQ6I7gB8zpeui7uPYnIuIi70cp59sMDOAIpoACXRH+gGOc0vgJgMBc1pgCKgACXRD8gFOcOrVT\nCATm4ELwABMIf8CpnNsHjEIABJp1P1P4A6qQeX/rIQBMIgACTRL8gDlpAIFRWAUUaIrgBwBwPA0g\n0AzhD1jKHA2gaaBACzSAQPUEPwCAeWgAgaoJf8DaTmkDtYBA7TSAQJUEP2BtFoIBRqABBKoj/AFb\n0wICvRIAgaoIfwAAy8lSytZjiIiIzCy1jAVYn+AH1OiUJrCU8/kGAgwhM6OUsuhBkQYQ2JzwBwCw\nDgEQ2Mz9TOEP6Jbfb0CNBEBgEw6MAADW5xxAYFWCH9CKY8//u3m/c8c3wERrnAPoOoDAaoQ/oAWu\nBwj0TAMILE7wA1pzagjUAgLHsAoo0DzhD2jRKQFQgwjUzBRQYBGCHzAi4Q+onQYQmJ3wB4xqXwD0\nexGohQAIzMpBDtAizR0wCgEQmIWLugOt2hX+ljgH0O9IoAYCIHAyBzVAiy6u/lz/eM5tA9RIAASO\npvUDWnZbSFsqwPmdCWzNdQCBoziIAXqxVNjbt13XBAT2cR1AoDpaP6AnW0zV9DsU2JLrAAKTOWgB\nerF08HMOIFArDSBwJ60f0JMawpnfqcBWnAMI3MpBCtCLtYPfXY/nXEDgpjXOARQAgZ0EP6BHa4XA\nqY8jBALXCYDAJoQ/oGdLhsBDty0AAtdZBRRYnfAHsJ7M+1sPARiMVUCBiBD8AABGoAEEhD+AmUyd\nAnpx9QdgbQIgDMzlHYDR1Bi6TAMF1iQAwqAEP4BlTAmZNQZRYAwCIAxG6weMbI3gdcxjaAGBtbgM\nBAxE8ANGtnT4O2b71+9Tyvl8gwGa5DIQwCy0fgDLmiNcagGBNWgAoXOCH8Ar5m4Bj93evvtpAWFs\nGkDgJMIfwLLmDoBaQGBpGkDokOAH8LAaz//bdz8tIIxrjQZQAITOCH8AD1sy/M3dAD4gBMKY1giA\n95bcOLAewQ9gXa7lB7TIOYDQAeEPYL8lgtrS4c+5gMBSBEBomMs7ANytxgCoPQS2IgBCowQ/gLsJ\nWgAPswgMNEbwAzhMjYvATL2/xWBgLK4DCDxE+AMA4BQCIDRC+AM4zlIN4BzbvWsbFoMB5uYyEFA5\nwQ+gLnMGSucoAmvTAELFhD+A09UasqaOSwsIzEkDCBUS/ABOt3Twu7j6c8z9ALaiAYTKCH8AdTs2\n+F2//6G0gMBcXAYCKiH4Acxr7abttsd78LXrtzl0fC4JAf1b4zIQAiBUQPgDWEYNIXCuMVzERZw7\nVoKurREAnQMIGxP+APpktVCgRs4BhI3czxT+ABa2VXBaKvx53QBOpQGEDXgBB1hHDVNAAWriHEBY\nkeAHsJ6ewtjN5+JcQOjTGucAmgIKKxH+ANbVUwAEmIsACAtzrh8ApxJmgbmYAgoLEvwAttFbYNr1\nfEwDhf64DAQ0SvAD2M4I4Q/gWAIgzEz4A2AOdwW/+5laQOBgzgGEGQl/ANvroTHr4TkAdXIOIMxA\n8AOoR8vh6ZixawGhHy4DAQ0Q/gDqMVr4AziUBhCOJPgB1Km1IDXHeLWA0AcNIFRK+ANgDq2FVaB9\nGkA4gOAH0Ibag9US49MCQvs0gFAR4Q+gHTUHwJrHBvTPdQBhAuEPgFMtHfxcFxCYQgCEWwh+AJxK\n4wfUxBRQ2EP4A+BUwh9QG4vAwA2CH0Db5gxdD7Z16Da3DH6mgUK71lgERgCEa4Q/gLYtEf4OfYyt\nWz8BENq1RgB0DiCE4AfQmmObuUO2vdb95mYxGOA2GkCGJ/wBtGNKK3dqEJty/5u3qSX8PSAAQps0\ngLAw4Q+gDbVOu6wt+D2gBQT2EQAZkuAHwClqDX4Ad3EZCIYj/AG0Z63A1VOw83oH7KIBZBheCAEA\nGJ1FYBiC8AfQrrVbuZ5awAgLwkBLLAIDJxL8ANom/AHMSwCkW8IfAFP1HPysCApcZxEYuiT8AfRh\njWDWc/gDuMk5gHRF8APoy5LhbLTgpwWE+jkHEAAY0pzh7Pq2Hvx9tPAH8IAACABUZa5wtms7gh8w\nOlNA6Y5poABtWqr14xWmgULdTAEFAIawZOsHwCusAkp3vLsJMCbh725myQACIACwuVPC28XVHwDu\nJgDSJS0gQFuODXCC3+G0gDA25wACAJs5tfkD4DAaQLqlBQSomwC3HS0gjEsABABWNdc5ewIkwOEE\nQLqmBQSoi9BWDy0gjMk5gADA4pYKfg+2K1gCTKMBpHtaQIC+CX/H0wLCeARAAGBxS4Y0ARBguskB\nMDMfy8zPZ+anrj7+eGZ+7uq/r2Tm5/bc77WZ+YnM/FJmPp+Zf3SuwQMAcBotIIwly8TpcZn5kxHx\nn0XEa0opz9z42s9GxG+XUv6bHff7HyLiM6WUX8jMexHx+0sp39xxuzJ1LHAML3AA25u7rdP+zcPp\nElCHzIxSyqIHrZMawMx8fUQ8HREf2XOTd0XEx3bc7zUR8VQp5RciIkop39oV/gAADiX8ARxu6hTQ\nD0fEByPikbeHMvOpiPh6KeXLO+73fRHxW5n5C1dTRf9WZr7q+OHC8by7CbCdua79xzLMkoFx3BkA\nM/MdEfGNUsoXIiKv/rvuvbGj/btyLyKejIi/WUp5MiL+n4j4qeOHCwC0ZulLQAAw3ZTrAL4pIp7J\nzKcj4lUR8erM/Ggp5X2Z+XhE/Fi8FPJ2+fWI+Fop5Z9dffyLEfGhfQ90cXHx8t/Pzs7i7OxswvBg\nuvNSvMsJ0AkBcF73M82WgZVdXl7G5eXlqo85eRGYiIjMfEtEfODBIjCZ+faI+FAp5a233OczEfFf\nlVJ+LTPP46VFYB4JgRaBYS0CIMC6LPzSDgEQtlXNIjC3eHfcmP6Zma/LzE9f+9Rfjoi/m5lfiIj/\nNCJ++sTHhJN4cQNol/C3LG+SQv8OagCXpAFkTV7gAJY3Z1gT/NbjjVLYTgsNIDTJixvAskz7bJc3\nSaFvAiAAMLslApsQCHA6U0AZmnc5AZY1V2gT/tZntgyszxRQAKBJx174fdd9hD+A+WgAGZoGEGB+\nxwa26/c7NkAyLy0grGuNBlAAZHhCIMB85mr9qIMACOsyBRQA6JrwB7AuAZDheXcTYD4CXV/MkoH+\nCIAAwCyOOW9PYARYlwAIoQUEONUcC79QJy0g9EUABAA2IfwBrE8AhCtaQIB1CYDt0AJCPwRAAACA\nQbgOINzgXU6A6Zz7NxazZWBZrgMIAFRL+ANojwAIN3h3EwB2M0sG2icAAgBH0QACtEcABACOIgCO\nSQsIbbu39QCgRueleIED2EPwA2iXBhAAWJzwB1AHARD2sBgMwHwEwL6YJQPtEgABgFUIgQDbcyF4\nuIN3OQFe4fw/rjNbBua1xoXgLQIDANxJ8APogymgcAfvbgKjE+LYxywZaI8ACAAsRngEqIsACBNo\nAQFgNy0gtMU5gADAIrR/APXRAMJEWkBgRBZ/YQotILRDAAQAABiEAAgA7KUBZCotILRBAIQDmAYK\njEYABOiLRWAAgJ2OCXGC39juZ3qzFCqnAYQDeWEDRnFomBP+AOqnAQQATiL4AbQjSyVtRmaWWsYC\nUzjZHejZ1FAn/LGL2TJwnMyMUsqiB5kaQADgIVNCneAH0CbnAMKRvLsJALuZJQP1EgABgIdoAAH6\nJQDCCbSAQI8EQOagBYQ6CYAAwMuEP4C+CYBwIi0g0BMBkDlpAaE+AiAAAMAgXAYCAIgI1/5jGfcz\nzZaBirgQPMzENBegB7eFO8GPYwmAMM0aF4I3BRQAgEV5kxTqIQDCTLy7CQBA7QRAACAiTPEEGIEA\nCDPSAgK9Eg45lWmgUAcBEACIiN0h7+LqDwB9EABhZlpAoCfCH3PSAsL2BEAA4GX7WkAA+iAAwgK0\ngECrhD2WpgWEbd3begAAwPYEP4AxaABhIVpAoBXCH2vTAsJ2BEAAGNxdAVBABOiHKaAAwE6CH0u6\nn2m2DGxAAwgL8sIGtEDQAxhHlkoOUDOz1DIWmJPzHIBaTQ1+AiJL8mYpvCIzo5Sy6MGjKaCwsPNS\nhECgSYIfQH9MAQWAQd0W8IQ/1uJNUliXAAgrML0FaInwB9AvU0ABYEBCHsCYNICwEi0gUAvX/aM2\npoHCejSAADCIKcFO+APomwYQVqQFBGom/LElLSCsQwAEgEFY9RMAARBWpgUEtiIAUjstICxPAASA\nAQh4AERYBAYAumbhF1pzP9NsGViQBhA24IUNWItLPgBwXZZKDkQzs9QyFliD8xyANd0MeoIftfNm\nKSPKzCilLHqQaAoobOS8FCEQWJygB8B1poACQKes+knLvEkKy9AAwoa0gMAShDsA9hEAAaAjVv0E\n4DYWgYEKaAGBUx0a6oRAWmExGEayxiIwzgEEgA4IgABMoQGESmgBgTmYAkqPtICMwmUgAIBJpoY6\n4Q9gbBpAqIgWEDiG1o8RaAEZgXMAAYBbCX8AHEIABICGCYCMwiwZmIcACBUxvQUAgCVZBAYAGmXh\nF0ZzP9ObpXAii8BAhUxzAaa6LdwJfvRIAKRnFoEBAG4l5DEab5LCaQRAqJB3NwEAWIIACACwg3YV\n6JEACJXSAgKnEF5O82D/2Y91Mg0UjicAAkDD9gUUweU4F1d/AHolAELFtIDAPoLKMnbtU/u5TlpA\nOI4ACAAdElqOZ98BPXMdQGiAdzmBfW6GFeHlNBrA9pgtQ0/WuA7gvSU3DgAsQyiZnzANjEAABIDG\nCCbzsj/bdj9TCwgHcA4gNMALG/DAXWFFmDmM/QWMRgMIAJ0QZqazr/qiBYTpNIDQCC9sQITgMocp\n+9B+BnplFVBoiNVAYWxTQ4nwcrvb9o9916YH37dSzrcdCJxojVVABUBojBAIY9JaHe/Bfrm+f1zu\noR/Xv28CIK1zGQgAQPg7wV2hj3b5fsJxnAMIjXEuIIxHADzcxdUfxpJ5f+shQPUEQAAAgEEIgNAg\nLSCMx/X/7nbsPtAW9kULCLcTAAGgEfsWLhFeTgt/tMf3DY5nERgAaIAD3t1u7heXyiDipRbQiqCw\nmwYQGmUaKIzDdet2O+W5j7zfeuD7B8cTAAGgYs79e9Rc015H3He9mPK9c91c2M0UUGjYeSle4GBg\nIwaYuZ7ziPuuB6b4wukmB8DMfCwi/nlEfK2U8kxmfjwi/tDVl787Iv5dKeXJW+77zyLi10spz5w4\nZgAYxq62a8SD2zmf84j7byTXv7/3M50yATcc0gA+GxHPR8RrIiJKKe958IXM/NmI+O077vt/Prgv\nMB8tIPRrX1AZaeVPwQ9gXpPOAczM10fE0xHxkT03eVdEfOzI+wIAO4y++Esr4W+kQN6Cm98Lb5LC\nw6Y2gB+OiA9GxGtvfiEzn4qIr5dSvnzofYF5aAGhXyOGi1aC383tP/j7aN+vtUzZr/Y93O3OBjAz\n3xER3yilfCEi8uq/694b+9u/u+4LAOwh/NWzrZofk7t5kxReMaUBfFNEPJOZT0fEqyLi1Zn50VLK\n+zLz8Yj4sYjYufjLbffddeOLi4uX/352dhZnZ2eTnwiMTgsIfRktSPTyfEcM7WuZ0rDa/7Tm8vIy\nLi8vV33MLAesjJSZb4mIDzxYyTMz3x4RHyqlvPXQ++74ejlkLMCjBEDoQ89T3XYdxC/1XNbYR7se\no9XvTe1OvQSE1UBpQWZGKWXRA7pTrwP47rgx/TMzXxcRP19KeeeJ2wYOpAWE/rUcLtYIfWu73jj1\n8pxqNMf1/1wSAl5yUAO4JA0gzEMAhD70tgLommPe8ty/Fr83rZijGRcAqd0aDeCky0AAAPVoMWSs\nNea1983Nc85a/N705K79701SOH0KKFAZ00ChD/sWsxAwdqthv9Qwhp5Z4AXmoQEEgAr1dKC79IH7\nlvuqp++2E2NBAAAgAElEQVRTC+bY394kZXQaQOiQFhDa1Vug6DX41TQGHqYphNtpAAGgEj1d/mHO\ng/Bap8LWMAaO401SRmYVUOiYFzhoTw8hcKngV2OzM/rqn1s9fyuC0iurgALAYFoOgCO0frykh9VP\nvUnKqARA6Jh3N4G1zBn8WgkTNwNQK+M+VQ3Ps+U3SmBrFoEBgMq0FCaOHeeu9qjl59zK2E/R43O8\nn+nNUoajAYTOeWGD9iw5jXJOc4S/Bx+3HC5aHvtUNT5HLSAcRwMIAJWYerB6yIFvTeGqlnHMraZ9\nPLc5fybnJgDCcTSAMAAtINTvkAPV2267Zhg5JBz0dCDe03O5zVw/k0uYM/xZDIbRCIAAsKFjw9Gu\n+7V8EN66np5nC4Hdzx4cz3UAYRDe4YR6rXmgOudj3dVE9q63xWBOGXvrU0Av4iJKOT9tQDCDNa4D\n6BxAGMR5KUIgVGqraZtLPWbLIWhkx/wc1hr8lrgd9MIUUAAY3FwHwC1MHZxTD8/15ves9uc0tfU7\nJvxl3j9uUNAYARAGYjEY4LqWzz+sVUv7YVfwa2n8wHEEQADoyLEH8KdOlxs1OLT6vHs7f/FYN5+z\nFpARCIAwGC0g1GfOKZhzbWfX1MARA8IxWthPc4yxhecJPMoqoDAgi8FAXU49kK5hNc4Rw8DU/V7r\n9Mql2+K5zbnq5133syIoW1ljFVABEAYlBEJd1jgYX+LAvbZQs7YWVs2ccwxbjn+pVT133V4AZCtr\nBEBTQAFgELumdp66vdG1tg9OeaNh6+e6RADcd1vnAtKze1sPANiG6wJCXWo4wJ6qlXHyita/Z0u2\nf63vGziUBhAABtTa+V+wJj/n9EwAhIFZERTatuZBqqbkNFvvv9a/d0uN/+Z2r3+fzJKhVwIgAFRi\n7YP0XecE7goqrYeHJS25MuWcWl34Ze0x1PBcYWlWAYXBeYcT6tLiqpIjayEA9jDdd8n9vOt+1z9n\ntgxrWmMVUIvAwOAsBgN1udnGTb0t25gytfPB17doeLe4bytGeI6wiymgAFCpXVM0r3+NdrQS/rY+\nV3GfrZrWi7hwSQi6YwooEBGmgkLthMC61TIV9NTHqPFn6pAxLTWF2oXhWYsLwQMAEVHngTkvqeV7\nU8s45nJMG7nUPtAC0hMBEIgIJ7lDS2qdpjeiJdupQ7Y717ZH+7ka7flChAAIAM0Q/Oqz5fej95+H\n2lYv1QLSCwEQeJkWEGBZLQS2WsbYwjmT0CIBEADgSMcEiJpDRy2tYg1j2EULSA+sAgo8woqgANNs\nvUjJnNurLXStcV7jMY9hRVCWZBVQAAAWVUvrd9PcAbDG5whbEAABAI60dajY+vGXMncovbmtU7Zv\nGiitu7f1AID6nJdiGijAApYONq1b+vn0tr/gGBpAAIAV1Br+aglFtYxjCi0gLRMAgZ1cEgLgYbsC\n3NbXAexFT88FaicAAgDc4XpAOeUC5bUGna3HtfXjH0MLSKucAwjs5VxAgNsdssJkjdM/594WUD8N\nIADAgQ5tBGu+/l8tWnxOWkBaJAACt3IuIMDttgiAc26zpuA1dSw1T6eF2gmAAAC32Bc0tl4ApscA\nOEVt49UC0hrnAAJ3ci4gwG6HNFY1qnVc+9y2CmtrzwW2ogEEALhFr8Gi5edV29i1gLREAAQmcS4g\nMKK5plrWFlhqdtu+uu37YR/DNAIgAMAOtQaKUc/9u2nX+Ld6Tq3vS8biHEAAgGtqPZgf6dp/La3y\n+WCc9zPNlqEJGkBgMi9sQO+WDB2nbruVQDSXrS6vMVVLIRWuEwABgOG1cjBv+ufdLq79WZsVs2mB\nKaDAQVwSAujJ2iHhlGDSc2hbistEwKM0gADAkASC+s05DXTu7/e+7XmTlNoJgMDBnAsItGxqC9fK\ntNBj9fzcgP0EQABgGMcEv16DUivPq+awvu98QxeGp2bOAQSO4lxAoEdLX1uuhvP/Wgl+h9r6eW39\n+DBVlkqmcmVmqWUswDQCINCCKY3eCMFvyW2uYa7zAdd8/qWcr/ZY9CEzo5Sy6AGWKaDA0ZwLCNTs\nevC7+f+bt9t3vznHsub9APbRAAIn0QICNTq2CVo6cNUQBFsNlbeN+9DnpAWkVms0gAIgcBIBEKjJ\nHJcEqGnaZS1NZC2W+L4tvU8EQA5hCihQPdNAgVrUGG56v5RETexnmEYABACGsy8s9D7V8rbnXeN4\nb6qpnZ3KJSGojQAInEwLCNTgkHCz9MIvc21r6SmgU/ZNTWoeG7RCAAQAurTr8g9rNV01BpVDFsa5\nbdXUHmkBGYkACMxCCwjUoMfgMnebeOhCOS3sy1amsEINrAIKzMaKoECrapz+udT2bm73kO1vHbJ2\nPX5L5wVaEZS7rLEK6L0lNw6M5bwUIRBoTq3hb+mwdcz2b06pXdPW4RN6YQooADC0GoNFjWO6qZYp\noi01gM4FpAYCIDAr5wICLdo6xNy05XgOOZ9uzXHetsorMJ0poADA0ASIV+w7x+6ufbT0IixbfI/m\nnu76YBv388KbpWxKAARm51xAYDRLTUPcsmGrJRjXMI5Tw2ANzwEesAoosAgBEGjFsatgrnkO3FoB\nYsrjHHqbU/bTXOOZ01yrpmoB2WWNVUAFQGARAiDQkmOCxtorYtZwvt3Ur891n9vut1Wrdujj7ru9\nAMguawRAi8AAi/DCBvSs1umSc5mr5VrK1ovkLHl7WJoACAAM75Sphmudq1dbkLj+vJeYqrn2OZCH\nuLj251hmyrAVARBYjBYQaMHUA/law8hSbnu+S++LKVNQa/l+1DIOmEoABAC4Qy0H+bWMY5c5QvQh\nYbzmfTGVFpAtCIDAorSAQO2OXdVyKzUtBjPnNnsMfq2NlzEIgAAAB9gyDG4RKA55zLlX+mzFrvFP\nfU5aQNYmAAKL0wICPWg9pJxqiRB4yCIyrdD6Ubt7Ww8AAGBrxy4E42D/OD3us1NXBPVmKWvRAAKr\n8MIG9GiL6ZhbWuqcQKuwwno0gAAAcfh0RIFkv5oWzWmFFpC1aACB1XhhA2p3zPlra9s6UI36vKEX\nAiAAwBUhY35z7dMRvjdWBGUNAiAAwIG2Xvxl6zB06OOfMt6t9zX0RgAEVmUaKFAri5EcZo39YF/D\n/ARAAGB4LQeNlsd+m16f111MA2VpAiCwOi0g0LJagompkcAxBEAAYHitB6kWxr/meYOt0wKyJNcB\nBDZxXooXOIAZbL0YzbGXzhg54MGWNIAAANdYCGYd9uHtvEnKUgRAYDPOBQRq0uo5dVuP+ZjHb3Vf\nQw8EQACAIwgwr1hiX9i/WkCWIQACm9ICArVpMXjcNuZdX9PAwbgEQACAI11c+1OjfeHvtq+vwXmW\n03iTlCVYBRQA4JpTA93agbCWMDX38x45AAp+LEkDCGzOCx3Qg62awEOnfx5zmyWMHPD2OS/FayKL\nEwABAG44NZysHW72PV4t4zj1tiMQ/FiLKaBAFVwYHuB4N8PUoeFqzfZS8HuY4MfaslTyQ5eZpZax\nANsQAIHatNAEnhr+7treXNvRDj5K+OOmzIxSyqIHRKaAAtXwQgj0YsuVQWsIgHOH0t44148tCYAA\nAA2rOVwdG4Rrfk6nEPyogSmgQHVMBQVqU2uIues6f6due4sms+fwB3dZYwqoAAhURwAEanJsIFky\nyCwZ/KY+XsuPsybBj0M4BxAYkhdLoCa1rajZY0iK6PN5eT2jRgIgAMCMWl38Zattr/kYa3GuHzUT\nAIEqeeEEanJMC7i2ngJUqwQ/WiAAAgBMUHPAqnlsU2x52Yy5CH60YvIiMJn5WET884j4Winlmcz8\neET8oasvf3dE/LtSypM37vP6iPhoRPyHEfF7EfHzpZT/fs/2LQIDPMRiMEAtalgF9K7ttTpNU/CD\nV9S2CMyzEfH8gw9KKe8ppTx5Ffp+KSI+ueM+34qI95dS3hARfzwifiIzv/+UAQPj8KIKtGzOYNN6\nSOqV1ylaNCkAXjV5T0fER/bc5F0R8bGbnyylfL2U8oWrv/9uRHwpIr7nuKECAIxl6tTIli/T0GK4\nda4fLZvaAH44Ij4YEY/8pGfmUxHx9VLKl2/bQGY+ERE/GBH/5LAhAiPzAgu07NRw02I46pngRw/u\nDICZ+Y6I+MZVk5dX/1333tjR/t3YxndGxC9GxLNXTSAAQBO2DmGHPP7N2659HmLPBD96cW/Cbd4U\nEc9k5tMR8aqIeHVmfrSU8r7MfDwifiwintx358y8Fy+Fv/+xlPL3b3ugi4uLl/9+dnYWZ2dnE4YH\n9O68FAvCAM1aY4XLXdt/8Lk5Hl/wg2VcXl7G5eXlqo85eRXQiIjMfEtEfKCU8szVx2+PiA+VUt56\ny30+GhG/VUp5/x3btgoosJcACGxp6wA1VwCr8XqGNYdL4Y+11bYK6C7vjhvTPzPzdZn56au/vyki\n/mxE/OeZ+fnM/NxVaAQ4iBdhYCs1B5RDXVz7M+W2a6hx/zrXj54d1AAuSQMI3EULCGxhq/Zt7sef\n+hhbBLJaQqDQx9ZaaAABVuOFGWBZwh/0TwAEAFhB7Q1gLUFsbaZ7MhoBEABgj5ZDUctjX4vgx4gE\nQKApXqyBlh2yCMupj3Pbx6PT+jGyKdcBBAAYUutTI2se+xZjE/rAKqBAo6wICiytlmvgLXHu4NbB\nUPiD3dZYBVQDCABww9YB6bpDW8hDQmVNz3Mpgh88zDmAQJO8oANLOSYU1bLC56Hb2yoArvW4Xivg\nUQIgAMA1x4aopULN1EavpTZv6bFa5AX2EwCBZnlxB2pwapjZOrit+fhrBFWvDXA75wACAFzZYirn\n1gFwLYIf1EEDCDTNCz4wp2PP/6vlHMBjLTmOpVs/0z3hMC4DATTPJSGAuR276uacC8hM3dac4arW\nbe0j+NEbl4EAmOC8FCEQmEVLjVzNYc10T6iXKaAAABu62SBu0fwda994l57yCRxPAAQAONJcQefY\n8+RqCIERu8cx99ic6wfzcA4g0A3TQIE5HHv+3zH3P1Wt00Dn3JbQx0jWOAdQAwgAsMPFtT+33WbK\n55Yyx2Od0j461w/aYxEYoBsWgwHmsC/U3fx8LdMveyT4wXI0gAAAC1i6IZtr+6eef7jEuX7AcgRA\noCsOHIClHNIALt0O1nzZhmO3ZZEXWIdFYIDumAYKLOGUKaBzrha6lLkWvzl0jEIfvMIiMABHcDAB\nzM35fw+77dqFh+wbv69hfRpAoEtaQGBOc6z2WWsLuMX0T8EPdtMAAhzJwQUwp54bv7UvAeH3M2xL\nAAQAuMVcl4CotQGcc5u3bcciL1AH1wEEAFhBrS2ii7nDWJwDCHTNuYDAXGpYBXSJbS65HeEPDuMc\nQACAjS19Qfdj1Tam6+Mx3RPqJQACXXMAAtSglrC2xgXq/d6FugmAAABHaC3U7bpQ+5zPoZTzKOV8\ntu0By3AOIDAE5wICx5rjGoCn3u/UbS/5uEIfzMc5gAAAG9oXnGpp/7Ym/EF7BEBgCM5JAea25bUA\nD9323I9ruie0SwAEANhj7lDVQ3Mo+EHbBEBgGFpAYE6Hhrk1wt/S5/oJf9A+ARAAYAVbt3+nPL7g\nB/0QAAEAbrHV6prHcmkH4DYCIDAU00CBLW0dGKc+vuAH/XIdQGA4rgkIHGPO6wGeet9jHueQ8Ads\nY43rAAqAwJCEQOBQ1wPUXOFt60bwOsEPtrdGALy35MYBAHpRU1ibm/AH49AAAsPSAgI1uDlFc82g\nKfhBXUwBBViQAAjUaukQKPhBndYIgFYBBYZlRVBgRMIfjE0ABACozBINoEs7ABECIDA4LSAwAsEP\neEAABACo0BwtoNYPuMkiMABhQRigboeGQaEP2mQRGAAADgqAwh9wGxeCBwDogOAHTGEKKMAV00CB\nFuxqA4U/6MMaU0A1gAAAjbgZ/gQ/4FAaQIBrtIBA7S7iQvCDTlkEBgCAhwh/wCk0gAA3aAGBGp07\nToLuaQABABD+gNloAAF20AICNRD8YCwaQACAQQl/wBI0gAB7aAGBLQh+MC7XAQQAGITgB6zBFFAA\ngI0Jf8BaTAEFuIVpoMCSBD/gOovAAAB0SvgDtqABBLiDFhCYk+AH7GMRGACATgh+QA1MAQW4g4M2\n4FR+jwC10AACACxE8ANqowEEmMBBHHAovzeAGmkAAQBmJPgBNbMKKMABrAgK7CP4AadyHUAAgAYI\nf0ArNIAAB9ICAg8IfsCcNIAAAJUS/oAWWQQGAOAAgh/QMlNAAY5gGiiMR/ADlmYKKABABYQ/oBca\nQIAjaQGhf4IfsCYNIADARoQ/oEcaQIATaAGhP4IfsBUNIADAioQ/oHcaQIATaQGhfYIfUIM1GkDX\nAQQAhiX4AaMxBRTgRA4goU3+7QIj0gACAEMR/ICRaQABgGEIf8DoLAIDMBOLwUC9BD+gBRaBAQA4\ngeAH8DBTQAFm4kAT6uLfJMCjNIAAQFcEP4D9NIAAM3LgCdvybxDgdhpAAKB5gh/ANFYBBViAFUFh\nHYIf0JM1VgE1BRQAaJLwB3A4DSDAQrSAsAzBD+iVBhAA4BrhD+A0FoEBAKon+AHMQwMIsBAHrDAP\n/5YA5qMBBACqJPgBzM8iMAALsxgMHEbwA0ZlERgAYCjCH8CyNIAAK9ACwu0EPwANIAAwAOEPYD0a\nQICVaAHhYYIfwMPWaACtAgpAl74ZEZdX/39NRJxd/Z/tCX4A29EAAqxIC7i834uIT3zHd8RXHn88\nnnrzm+MNb3xjPP/FL8Yvf/az8X3f/nb86RdfdP7DhoQ/gP00gABwoE98x3fE//e93xtfeu65eOKJ\nJ17+/AsvvBBve9vb4hNf/Wq8+8UXtxvgoAQ/gDp4ExSAbnwzIr7y+OPx3I3wFxHxxBNPxHPPPRdf\nefzx+L83Gd24hD+AegiAACtyILysy4h46s1vfiT8PfDEE0/EU29+c1yuOKaRnZfiZx6gMgIgAN34\nZkS84Y1vvPU2b3jjG+N31hnOsAQ/gHoJgAArc2C8nNdExPNf/OKtt3n+i1+M164znCH5+Qao2+QA\nmJmPZebnM/NTVx9/PDM/d/XfVzLzc3vu9/bM/NXM/LXM/NBcAweAm84i4pc/+9l44YUXdn79hRde\niF/+7GfjbMUxjULrB9CGQxrAZyPi+QcflFLeU0p5spTyZET8UkR88uYdMvOxiPgbEfGnIuINEfHe\nzPz+04YM0D4Hyst4TUR837e/HW9729seCYEPVgH9vm9/O169yej65ecZoB2TLgORma+PiKcj4r+N\niPfvuMm7IuKtOz7/wxHxr0op/+ZqOx+PiB+NiF89arQAcIc//eKL8YmvfjV+4Ad+YO91AJmH4AfQ\nnqnXAfxwRHww4tHTJjLzqYj4einlyzvu9z0R8bVrH/96vBQKAYZ3XooLwy/gsYh494svxjdffDE+\n89xz8ennnovXRsRPRGj+ZiL4AbTrzgCYme+IiG+UUr6QmWcRcfNo5b0R8bE5BnNxcfHy38/OzuLs\n7GyOzQIwoNdExH+x9SA6JPwBzOfy8jIuLy9Xfcwsd/wiz8yfjogfj4hvRcSr4qU3UD9ZSnlfZj4e\nEb8REU+WUn5zx33/WERclFLefvXxT0VEKaX8zI7blrvGAtAjLSAtEPwAlpeZUUpZ9MDgzkVgSil/\npZTyvaWU/ygi3hMR/7CU8r6rL//JiPjSrvB35Z9GxH+cmX8wM3/f1f0/NcfAAYB1CH8A/Tj1OoDv\njhvTPzPzdZn56YiIUsq3I+IvRsQ/iJdWEP14KeVLJz4mQFccXFMrl3YA6M+dU0DXYgooMDLTQKmN\n4AewvjWmgE5dBRQAGIDgB9A3DSBAJbSAbEnwA9heFYvAAAB9E/4AxqEBBKiIFpA1CX4AddEAAgCL\nEP4AxqQBBKiMFpAlCX4A9bIKKAAwC8EPgAhTQAGq40CdufmZAuABDSAAdErwA+AmDSBAhRy4cyo/\nQwDsogEEgI4IfgDcRgMIAJ0Q/gC4i8tAAFTMJSGYQvAD6IPLQAAAewl+ABzKFFCAijnAZx8/GwAc\nQwMIAA0R/AA4hQYQoHIO+HnAzwIAp9IAAkDlBD8A5mIVUIBGWBF0PIIfwFjWWAXUFFAAqJDwB8AS\nNIAADdEC9k/wAxiXBhAABiL8AbA0i8AAwMYEPwDWYgooQGNMA+2H4AfAdaaAAkCnhD8AtqABBGiQ\nFrBdgh8A+2gAAaAjwh8AW9MAAjRKC9gOwQ+AKTSAANA44Q+AmmgAARqmBayX4AfAodZoAF0HEABm\nJPgBUDNTQAEaJmzUxfcDgNppAAHgRIIfAK3QAALACYQ/AFpiERiADlgMZn2CHwBzswgMAFRG8AOg\nZaaAAnRAKFmH/QxA6zSAAHAHwQ+AXmgAATohpCzDfgWgJxpAANhB8AOgRxpAgI4ILfOwHwHolQYQ\nAK4IfgD0znUAATrkuoCHEfwAqMEa1wE0BRSAoQl/AIzEFFAAhiT4ATAiDSBAh4Sb29k/AIxKAwjA\nMAQ/AEZnERiAjlkM5iWCHwAtsAgMAJxI+AOAV2gAATo3agso+AHQGg0gABxB+AOA3TSAAAMYpQUU\n/ABo2RoNoFVAAWie4AcA05gCCjCAngNSz88NAOamAQSgSYIfABxOAwgwiJ4CU0/PBQDWpAEEoBmC\nHwCcRgMIQBOEPwA4nctAAAymtUtCCH4AjMJlIAAYluAHAPMzBRRgMC0EqxbGCAAt0gACUA3BDwCW\npQEEGFCNQavGMQFAbzSAAGxK8AOA9VgFFGBgW64IKvgBwMPWWAXUFFAAVif8AcA2NIAAg1uzBRT8\nAGA/DSAA3RD+AGB7FoEBYFGCHwDUQwMIMLglA5rwBwB10QACMDvBDwDqZBEYACJinsVgBD8AOJ5F\nYABohvAHAPXTAALwsmNaQMEPAOahAQSgasIfALRFAwjAQ6a0gIIfAMxvjQbQKqAATCb4AUDbTAEF\n4CH7Qp7wBwDt0wACcCvBDwD6oQEEYC/hDwD6YhEYAACACrgMBAAAALMRAAEAAAYhAAIAAAxCAAQA\nABiEAAgAADAIARAAAGAQAiAAAMAgBEAAAIBBCIAAAACDEAABAAAGIQACAAAMQgAEAAAYhAAIAPD/\nt3e/IZbVdRzH3x+dXbTtj8HSH93S1lxQsWCL2kCMLDOyRQUNkfRBUEYosUVFG0RQlJFREIhJ9MwI\nWSKmJZP8s/vMNEZ3p9ZEE1FclAolIqnd+PbgnpHb3Zm9M+feMzPMeb8e3fu753fuOfDhO/d7zplz\nJKknbAAlSZIkqSdsACVJkiSpJ2wAJUmSJKknbAAlSZIkqSdsACVJkiSpJ2wAJUmSJKknbAAlSZIk\nqSdsACVJkiSpJ2wAJUmSJKknbAAlSZIkqSdsACVJkiSpJ2wAJUmSJKknbAAlSZIkqSdsACVJkiSp\nJ2wAJUmSJKknbAAlSZIkqSeW3QAmOSXJXJLZobFbkjyeZD7JrUvM25Pkj0kOJ7kryeZpbLgkSZIk\naWVWcgbwC8CRhTdJPgTsBi6qqouA20YnJDkTuAXYWVXvAmaA6ybaYmmVHThwYK03QVqU2dR6Zj61\nXplN9d2yGsAk24CPAz8dGv4ccGtVHQeoqr8tMf1UYEuSGeA1wNH2myutPv9QaL0ym1rPzKfWK7Op\nvlvuGcAfAl8GamhsB3BJkoeSPJjkvaOTquoo8APgWeB54OWqum/CbZYkSZIktTC2AUxyBfBiVT0G\nZOijGeCNVbUL+Apw9yJzzwCuBM4GzgRem+T6aWy4JEmSJGllUlUnXyD5DvAp4DhwOvA64JfAVuB7\nVXWwWe4p4P1V9fehudcAl1fVZ5r3NzTL3LzI95x8QyRJkiRpg6uqjF+qvZllbMBeYC9Akg8CX6qq\nG5PcBFwKHEyyA9g03Pw1ngV2JTkN+DfwYeCRJb6n0x2VJEmSpL6b5DmAPwO2J5kHfg7cCJDkrUn2\nA1TVw8A+4FHgEINLSO+caIslSZIkSa2MvQRUkiRJkrQxTHIGcKwk1zQPgf9vkp1D4x9J8ockh5I8\n0jxTcHTubJLDS6z37CT/ah5MP5fk9i73QxtTV/lsPv9akieTPJ7ko13tgzamNtlMck+SR5PMJ7k9\nyQmX1Vs7NamustksZ93URFaazySnJ9nfZG6+ue/FYuu1dmoiXWWzWXbFtXPs/wBOaB64GvjJyPhf\ngU9U1QtJLgTuBbYtfJjkauAfY9b9VFXtHLOMdDKd5DPJ+cAngfObefclOa883a7la5PNa6vqnwBJ\n9gHXssjdmbF2ajKdZNO6qSlpk8/vV9XBDJ5X/UCSy6vq3kXWbe3UJDrJZtva2WkDWFVPNBuXkfFD\nQ6//lOS0JJuq6liSLcAe4LMs/uNlgTeN0UQ6zOeVwC+q6jjwTJIngfcBv+9iP7TxtMnm0A/sTcBm\n/v+5rcOsnWqtw2xaNzWxFvl8BTjYjB9PMsfQAd8R1k611mE2W9XOTi8BXY4MHhUxV1XHmqFvAbcB\nr4yZek5zGv7BJBd3upHqrZb5PAt4buj9882YNDWLZJMkvwVeYHCGet8SU62d6lTLbFo3tSoWy2cz\nfgawG7h/ianWTnWqZTZb1c6JzwAm+R3w5uEhBkf3vl5Vvx4z90Lgu8Blzft3A+dW1ReTnMPSR1uO\nAm+vqpea62h/leSChaOM0oI1yqc01jSzuaCqPpZkM3AXg8f0jP6xsHZqrDXKprQsXeQzyakM7mj/\no6p6ZpGp1k6NtUbZbGXiBrCqLhu/1ImSbGPwQPkbhnboA8B7kjwNbALelOSBqrp05DuPAS81r+eS\n/AXYAcy12wttVGuRTwZHX9429H5bMya9asrZHF7vf5LMMrgs5P6Rz6ydGmstsol1U8vUUT7vBJ6o\nqh8v8Z3WTo21FtmkZe1czUtAXz1bkuQNwH7gq1X10MJ4Vd1RVduqajtwMYMdHv1xTZKtSU5pXm8H\n3gk83fUOaEObWj6BWeC6JJuTvINBPh/udvO1gY3NZpItSd7SvJ4BrgD+fMKKrJ2arqllE+umpm9s\nPnEjo3YAAADeSURBVJvPvg28vqr2LLkia6ema2rZpGXt7PoxEFcleQ7YBexPck/z0c3AucA3Mrg1\n9FySrWPWtTvJN5u3lwCHm3+IvBu4qape7mYvtFF1lc+qOsIgl0eA3wCf9052WokW2dwCzCZ5jMER\n6ReBO5p1WTs1NV1l07qpaVhpPpOcBewFLhga/3SzLmunpqarbLatnT4IXpIkSZJ6Ys3vAipJkiRJ\nWh02gJIkSZLUEzaAkiRJktQTNoCSJEmS1BM2gJIkSZLUEzaAkiRJktQTNoCSJEmS1BM2gJIkSZLU\nE/8DYKZ1Jihh0iUAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(1, 1,figsize=(15,15))\n", "r1,r2,c1,c2 = 0,200,0,200\n", "plt.pcolormesh(X[r1:r2,c1:c2],Y[r1:r2,c1:c2], (bathy.mask).astype('int')[r1:r2,c1:c2])\n", "plt.scatter([X[20,40]], [Y[20,40]], c = \"w\" , s =80)\n", "\n", "lat = 47.5\n", "lon = -123.7\n", "\n", "model_lons = X\n", "model_lats = Y\n", "grid = 'NEMO'\n", "tols={'NEMO': {'tol_lon':0.0104, 'tol_lat':0.00388}}\n", "i_list, j_list = np.where(\n", " np.logical_and(\n", " (np.logical_and(model_lons > lon - tols[grid]['tol_lon'], model_lons < lon + tols[grid]['tol_lon'])),\n", " (np.logical_and(model_lats > lat - tols[grid]['tol_lat'], model_lats < lat + tols[grid]['tol_lat']))))\n", "plt.scatter([X[20,40]], [Y[20,40]], c = \"w\" , s =80)" ] }, { "cell_type": "code", "execution_count": 74, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(array([148]), array([139]))\n" ] } ], "source": [ "atmos_grid = xr.open_dataset('https://salishsea.eos.ubc.ca/erddap/griddap/ubcSSaAtmosphereGridV1')\n", "#print(atmos_grid)\n", "nemo_lon_lat = [-123.7665023803711,49.30873107910156]\n", "\n", "atmos_ji = find_closest_model_point(\n", " nemo_lon_lat[0], nemo_lon_lat[1],\n", " atmos_grid.longitude.values - 360, atmos_grid.latitude.values,\n", " grid = \"GEM2.5\")\n", "print(atmos_ji)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [], "source": [ "model_lons = X\n", "model_lats = Y\n", "\n", "tols={\n", " 'NEMO': {'tol_lon': 0.0104, 'tol_lat': 0.00388},\n", " 'GEM2.5': {'tol_lon': 0.016, 'tol_lat': 0.011},\n", " }\n", "\n", "grid = 'NEMO'\n", "land_mask = bathy.mask\n", "j_list, i_list = np.where(\n", " np.logical_and(\n", " (np.logical_and(model_lons > lon - tols[grid]['tol_lon'], model_lons < lon + tols[grid]['tol_lon'])),\n", " (np.logical_and(model_lats > lat - tols[grid]['tol_lat'], model_lats < lat + tols[grid]['tol_lat']))))\n", "if len(j_list) == 0:\n", " raise ValueError(\n", " 'No model point found. tol_lon/tol_lat too small or '\n", " 'lon/lat outside of domain.')\n", "try:\n", " j, i = map(np.asscalar, (j_list, i_list))\n", "except ValueError:\n", " # Several points within tolerance\n", " # Calculate distances for all and choose the closest\n", " dists = haversine(\n", " np.array([lon] * i_list.size), np.array([lat] * j_list.size),\n", " model_lons[j_list, i_list], model_lats[j_list, i_list])\n", " n = dists.argmin()\n", " i, j = map(np.asscalar, (i_list[n], j_list[n]))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "numpy.ndarray" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(X)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.5.1" } }, "nbformat": 4, "nbformat_minor": 0 }