diff --git a/notebooks/2020-05-19_CSV_data_wrangling.ipynb b/notebooks/2020-05-19_CSV_data_wrangling.ipynb new file mode 100644 index 0000000..3793839 --- /dev/null +++ b/notebooks/2020-05-19_CSV_data_wrangling.ipynb @@ -0,0 +1,1421 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Data wrangling: CSV to NetCDF\n", + "A user came to us as she had issues with some model outputs. The model outputs to CSV. That's a land model, so there are outputs only on the land points. The data is organised as so for each point, the data for each successive year is in a new row and the data for each month is in the column of this month (Jan, Feb, etc.) for that row. See below the output after reading in the data to make it clearer.\n", + "\n", + "She wanted to store the data in a (lat, lon, time) netcdf file. She would like the longitude and latitude to cover both the land and ocean (it's a 0.5° regular grid) and we also need to figure out how to combine the year and month information into one time information.\n", + "\n", + "Xarray does not read in CSV, we will need to read in the file with Pandas and then figure out how to convert it to DataArray and write it to NetCDF." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Get the data and dimensions" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "import xarray as xr" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
LonLatYearJanFebMarAprMayJunJulAugSepOctNovDec
039.75-1.2519010.0400.0390.0690.1440.0560.0220.0170.0050.0310.0270.1320.204
139.75-1.2519020.1580.0140.0520.0860.0990.0460.0130.0110.0240.0860.2210.089
239.75-1.2519030.0120.0100.0240.1400.0560.0380.0150.0150.0380.0200.0330.096
339.75-1.2519040.046-0.0020.0200.0510.1110.0320.004-0.0010.0240.0480.1300.049
439.75-1.2519050.0350.0000.1660.1640.1000.0390.0220.0200.0100.0240.1690.110
................................................
680684592.7540.7520110.0000.0000.0040.0050.0010.0190.0010.0010.0010.0030.0030.000
680684692.7540.7520120.0000.0000.0030.0050.0000.0200.0220.0080.0200.0030.0010.000
680684792.7540.7520130.0000.0000.0080.0010.0080.0170.0160.0210.0010.0020.0020.000
680684892.7540.7520140.0000.0000.000-0.000-0.0000.0040.0270.0040.0010.0030.0030.000
680684992.7540.7520150.0000.0000.0020.0060.0080.0000.0150.0050.0020.0020.0040.000
\n", + "

6806850 rows × 15 columns

\n", + "
" + ], + "text/plain": [ + " Lon Lat Year Jan Feb Mar Apr May Jun Jul \\\n", + "0 39.75 -1.25 1901 0.040 0.039 0.069 0.144 0.056 0.022 0.017 \n", + "1 39.75 -1.25 1902 0.158 0.014 0.052 0.086 0.099 0.046 0.013 \n", + "2 39.75 -1.25 1903 0.012 0.010 0.024 0.140 0.056 0.038 0.015 \n", + "3 39.75 -1.25 1904 0.046 -0.002 0.020 0.051 0.111 0.032 0.004 \n", + "4 39.75 -1.25 1905 0.035 0.000 0.166 0.164 0.100 0.039 0.022 \n", + "... ... ... ... ... ... ... ... ... ... ... \n", + "6806845 92.75 40.75 2011 0.000 0.000 0.004 0.005 0.001 0.019 0.001 \n", + "6806846 92.75 40.75 2012 0.000 0.000 0.003 0.005 0.000 0.020 0.022 \n", + "6806847 92.75 40.75 2013 0.000 0.000 0.008 0.001 0.008 0.017 0.016 \n", + "6806848 92.75 40.75 2014 0.000 0.000 0.000 -0.000 -0.000 0.004 0.027 \n", + "6806849 92.75 40.75 2015 0.000 0.000 0.002 0.006 0.008 0.000 0.015 \n", + "\n", + " Aug Sep Oct Nov Dec \n", + "0 0.005 0.031 0.027 0.132 0.204 \n", + "1 0.011 0.024 0.086 0.221 0.089 \n", + "2 0.015 0.038 0.020 0.033 0.096 \n", + "3 -0.001 0.024 0.048 0.130 0.049 \n", + "4 0.020 0.010 0.024 0.169 0.110 \n", + "... ... ... ... ... ... \n", + "6806845 0.001 0.001 0.003 0.003 0.000 \n", + "6806846 0.008 0.020 0.003 0.001 0.000 \n", + "6806847 0.021 0.001 0.002 0.002 0.000 \n", + "6806848 0.004 0.001 0.003 0.003 0.000 \n", + "6806849 0.005 0.002 0.002 0.004 0.000 \n", + "\n", + "[6806850 rows x 15 columns]" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "fname = \"/g/data/w35/lt0205/research/lpj_guess/runs/CRUNCEP/mgpp.out\"\n", + "df = pd.read_csv(fname, header=0, delim_whitespace=True)\n", + "df" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "months=list(df.columns)\n", + "months=months[3:]" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "lons = np.unique(df.Lon)\n", + "lats = np.unique(df.Lat)\n", + "years = np.unique(df.Year)\n", + "nyears = len(years)\n", + "nrows = len(lats)\n", + "ncols = len(lons)\n", + "nmonths = 12\n", + "lons.sort()\n", + "lats.sort()\n", + "years.sort()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Original code.\n", + "The person requesting help sent us the codes she wrote to try and solve the issue. We put one of those below to discuss the good and bad of it.\n", + "\n", + "```python\n", + "out = np.zeros((ncols, nrows, nyears*nmonths))\n", + "for i,lon in enumerate(lons):\n", + " print(i, ncols)\n", + " for j,lat in enumerate(lats):\n", + " #print(j)\n", + " vals = df[(df.Lat == lat) & (df.Lon == lon)].values[:,3:]\n", + " if len(vals)> 0:\n", + " print(\"Reshape\")\n", + " vals = vals.reshape(nyears*nmonths)\n", + " out[i,j,:] = vals\n", + "\n", + "t1 = pd.to_datetime('1/1/1901')\n", + "time = t1 + pd.to_timedelta(np.arange(nmonths*nyears), 'M')\n", + "\n", + "ds = xr.Dataset(data_vars={\"mgpp\":([\"y\", \"x\", \"time\"],out)},\n", + " coords={\"lat\": ([\"y\"], lats),\n", + " \"lon\": ([\"x\"], lons),\n", + " \"time\": time})\n", + "ds.to_netcdf('test.nc')\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The good idea in this code is not to loop over the times but only looping over the spatial points. The bad idea is the original dataset only has the land points. So by looping on all (lon, lat) pairs, we are looping over points that do not have any data in the DataFrame.\n", + "\n", + "The original DataFrame has 59190 spatial points (length of the Frame / # of years). The loop goes through 194878 values, i.e about 3 times more than needed. \n", + "\n", + "So let's estimate the time it needs to run. We'll time the loop for 1 iteration and then multiply by the number of iterations." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "21.3 ms ± 2.33 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" + ] + } + ], + "source": [ + "%%timeit\n", + "out = np.zeros((ncols, nrows, nyears*nmonths))\n", + "for i,lon in enumerate(lons[0:1]):\n", + " #print(i, ncols)\n", + " for j,lat in enumerate(lats[0:1]):\n", + " #print(j)\n", + " vals = df[(df.Lat == lat) & (df.Lon == lon)].values[:,3:]\n", + " if len(vals)> 0:\n", + " vals = vals.reshape(nyears*nmonths)\n", + " out[i,j,:] = vals" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "100.03737333333332" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Time for total loop in minutes:\n", + "(30.8 * nrows*ncols)/1000/60" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "So we need something that runs faster than 100 min." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## First solution: a loop\n", + "Instead of looping through all the spatial points of the output array, we can loop through the rows of the DataFrame and find the correct place in the output array to put this data.\n", + "\n", + "First, let's create the output DataArray so we can easily find the location of points using longitude and latitude. Instead of using only the land points, we'll create the array for the full grid directly and fill with NaN." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "# Create the axes\n", + "time = pd.date_range(start=f'01/{years[0]}',\n", + " end =f'01/{years[-1]+1}', freq='M')\n", + "# We'll use a generic way to create a regular grid from [-180,180] and\n", + "# [-90, 90] when knowing the resolution. Feel free to reuse as needed.\n", + "dx = 0.5\n", + "Lon = xr.DataArray(np.arange(-180.+dx/2., 180., dx), dims=(\"Lon\"),\n", + " attrs={\"long_name\":\"longitude\", \"unit\":\"degrees_east\"})\n", + "nlon = Lon.size\n", + "dy = 0.5\n", + "Lat = xr.DataArray(np.arange(-90.+dy/2., 90., dy), dims=(\"Lat\"),\n", + " attrs={\"long_name\":\"latitude\", \"unit\":\"degrees_north\"})\n", + "nlat = Lat.size" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "Show/Hide data repr\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Show/Hide attributes\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
xarray.DataArray
  • Lat: 360
  • Lon: 720
  • Time: 1380
" + ], + "text/plain": [ + "\n", + "array([[[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", + "\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", + "\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", + "\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", + "\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", + "\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", + "Coordinates:\n", + " * Lat (Lat) float64 -89.75 -89.25 -88.75 -88.25 ... 88.75 89.25 89.75\n", + " * Lon (Lon) float64 -179.8 -179.2 -178.8 -178.2 ... 178.8 179.2 179.8\n", + " * Time (Time) datetime64[ns] 1901-01-31 1901-02-28 ... 2015-12-31" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "out = xr.DataArray(np.zeros((nlat, nlon, nyears*nmonths)),\n", + " dims=(\"Lat\",\"Lon\",\"Time\"),\n", + " coords=({\"Lat\":Lat, \"Lon\":Lon, \"Time\":time}))\n", + "out[:] = np.nan\n", + "out" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First an example with 1 row to see what's happening.\n", + "In that case, we would loop over each row and find the location of the data in the dataarray. That would mean a loop over 6.8M items!\n", + "\n", + "We use the longitude, latitude and year information contain in each row to find the locations to update in the `out` array using `loc`. " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "350 ms ± 5.79 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" + ] + } + ], + "source": [ + "%%timeit\n", + "row = next(df.iterrows())[1]\n", + "out.loc[dict( \n", + " Lon=row[\"Lon\"],\n", + " Lat=row[\"Lat\"],\n", + " Time=out.Time[(out.Time.dt.year==row[\"Year\"])])] = row[3:]" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "array([0.04 , 0.039, 0.069, 0.144, 0.056, 0.022, 0.017, 0.005, 0.031,\n", + " 0.027, 0.132, 0.204, nan, nan, nan, nan, nan, nan,\n", + " nan, nan, nan, nan, nan, nan])\n", + "Coordinates:\n", + " Lat float64 -1.25\n", + " Lon float64 39.75\n", + " * Time (Time) datetime64[ns] 1901-01-31 1901-02-28 ... 1902-12-31\n" + ] + } + ], + "source": [ + "# Just a check\n", + "print(out.sel(Lon=39.75, Lat=-1.25, Time=slice(\"1901\", \"1902\")))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see only the data for the year 1901 has been updated for this point only in `out`.\n", + "\n", + "The different notation for Time and the longitude or latitude is because we only look for one longitude and one latitude but we need several time indexes.\n", + "\n", + "But this is way to slow. It's 308ms times the number of rows (6.8M). \n", + "\n", + "To improve that time, we have to remember all the years for each point are together in order. So we could read nyears at a time and get the data for all years and months at once.\n", + "\n", + "There is still the problem, if we select all years for a point at once, we get a 2D Dataframe with the years in rows and the months in columns. To solve this, we can simply stack the Dataframe over the months columns. Since Pandas doesn't scramble rows when stacking, we can stack the original Dataframe and then use the fact all points have data for the same number of years and months. If we had some missing years for some points, we could modify the following code to allow for it but it would be slower." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "# We stack the months columns of the whole DataFrame at once since we\n", + "# don't have missing years for any point. Doing it once will save time.\n", + "df_stack = df[months].stack()" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.07 ms ± 88.8 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" + ] + } + ], + "source": [ + "%%timeit\n", + "#Example for one spatial point\n", + "rows = df[0:nyears]\n", + "# If we had missing years, we could add the missing years rows and then\n", + "# stack only the rows for the point here.\n", + "#rows_stack = rows[months].stack()\n", + "out.loc[dict( \n", + " Lon=rows[\"Lon\"].min(),\n", + " Lat=rows[\"Lat\"].min())] = df_stack[0:nyears*nmonths]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If each pack is 2.07 ms, then for the whole DataFrame we need 2 min (see below) which is very good!" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2.042055" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(df.index)//nyears*2.07/1000/60" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 1min 5s, sys: 552 ms, total: 1min 5s\n", + "Wall time: 1min 5s\n" + ] + } + ], + "source": [ + "%%time\n", + "#df_stack = df[months].stack()\n", + "for nr in range(0,len(df.index),nyears):\n", + " rows = df[nr:nr+nyears]\n", + " thislon = rows[\"Lon\"].min()\n", + " thislat = rows[\"Lat\"].min()\n", + " out.loc[dict( \n", + " Lon=thislon,\n", + " Lat=thislat)] = df_stack[nr*nmonths:(nr+nyears)*nmonths]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And a plot to check." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEWCAYAAACOv5f1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOy9eXwkZ33g/f3V0YdarWM093jGMx7fJzazmEAAbwDHmMSw4QjgAHEcHJIlECdsgF12IVk2b0iyywv7EoiXhJgsBEgg4ASDyZqwSRxwbGPjazzjY8ZzaUYzGh2tVnfX9bx/PFWlak1Lamla6m6pvp9Pf6Su86nq7t/zq98pSilSUlJSUtYWRrsHkJKSkpKy8qTCPyUlJWUNkgr/lJSUlDVIKvxTUlJS1iCp8E9JSUlZg6TCPyUlJWUNkgr/lDMQkSdE5Lp2jyMlJWX5SIX/GkREphKvQEQqifc3K6UuU0p9v93jnAsR+Q8i8riIlETkgIj8h1nrd4rIP4jItIg8JSKvSqzbIiJ3icgxEVEisnPWvn8uIs6se2TOM5asiPyZiEyKyHER+c1Z6+8QkX3hff7FBa5rvYjcJyKjIjIuIj8QkZcm1l8uIveIyCkRSRN0Us6KVPivQZRSvdELOAT8bGLZF9s9viYQ4B3AIHAD8B4ReUti/V8CDwNDwH8C/lpENoTrAuA7wBvmOf4fJO+RUsqfZ9uPAhcA5wL/FvhtEbkhsf7HwK8BP2riuqaAXwI2hNf2ceBvRcQK17vAV4FbmzhWSsq8pMI/5QxE5GCkLYvIR0Xkr0Tkf4ea9mMicqGIfEhERkTksIhcn9i3X0T+VESGReSoiHxsPs15KSil/kAp9SOllKeU2gd8E3hpeP4LgWuAjyilKkqprwGPEQp7pdQJpdQfAw+0aDjvAP6rUmpMKbUX+F/ALybG+mml1L1AtYnrqiql9imlAvQE56MngXXh+n1KqT8FnmjR2FPWMKnwT2mGnwX+Ai2IHgbuQX93tgG/C/xJYts7AQ84H7gauB745UYHFZG3heaNuV47FhqYiAjwMmYE4mXAc0qpUmKzH4fLm+XXROS0iDwkInM+IYjIILA1PP5Sz9XouI+iJ4u7gM8ppUbO5ngpKY1IhX9KM/yTUuoepZQH/BXaLPH7SikX+DKwU0QGRGQT8BrgN5RS5VBofQJ4S6ODKqW+pJQamOd1qImxfRT9Pf58+L4XmJi1zQRQbPJaP4U242wE/jPw50m7+yx6E8dfyrkaopS6EugD3gb889kcKyVlLlLhn9IMJxL/V4BTCTt4Jfzbi7Z728BwpL2jnwo2LsegROQ9aLPLa5VStXDxFFpwJukDSjRBaE4aDU1KdwNfBH4uPN9nE07g/xieKzr+os81y6lc95QTmoD+EvigiFzVzPFSUhaDtfAmKSlNcxioAevDp4R5EZGbqTcZzebSubR/Efkl4IPAy5VSRxKrngDOE5FiwvRzFfClZi6gAQptf0cp9W7g3bPGMRwe/+8T52rKJh863BfCBs6j3rSUknLWpJp/SstQSg0D3wX+u4j0iYghIrtF5BVzbP/FWVE1s19zCf6bgd8DXq2Uem7WMfcDjwAfEZGciPw74Erga4n9c0A2fJsN30fr3igiveHYrwd+AW17n4svAB8WkUERuRh4F/DnieNlwuMLYIdjavi7E5EXi8hPhvvkReQDwCbg/nC9hMfKRNchItlGx0pJWYhU+Ke0mneghdOTwBjw18CWFp/jY+gwzgcSZpPPJta/BdgTnv/3gTcqpU4m1leYMdk8xYzpCuB9wFFgHPhD4F0L5Dx8BHgWeB74v8AfKqW+k1j/3fD4LwHuCP9/+RzHygKfBkbDMdyINmkdC9efG+4fPVlUgH3zjC0lZU4kbeaSkpKSsvZINf+UlJSUNUgq/FNSUlLWIKnwT0lJSVmDpMI/JSUlZQ2yKuL8169fr3bu3NnuYaSkpHQBDz300Cml1IaFt5yb7ZJXVYKmtj2Fc49S6oaFt1xZVoXw37lzJw8++GC7h5GSktIFiMjzZ3uMKgFvaDKC+U94fv3Znm85WBXCPyUlJWUlEcCUJjfu0Gj6VPinpKSkLBIBMkaT0n++bhBtJBX+KSkpKYtEa/7Nqv6dSRrtk5KSkrJYRJt9mnkteCiRG8JWn8+IyAcbrH+diDwqIo+IyIMi8pOtuIRU809JSUlZJK3S/MMud58GXg0cQderuksp9WRis3uBu5RSSkSuRLfyvPhsz50K/5SUlJRFsiiH7/y8CHgmqk4rIl8GXocujAiAUmoqsX2BFrmQU7NPSkOqlQrVSmXhDRc4RkrK6kQwpbkXsD4010Sv2xIH2obugxFxJFxWfzaRfyciTwHfAn6pFVfQVs1fRG5H93dV6CbbtwA9wFeAncBB4M1KqbE2DbGrmS18c/n8ko8R7duMQI+2Xcr5UlK6gbA5Q7Obn1JK7ZnnULM5Q7NXSv0N8Dci8nLgvwKvavbkc9E2zV9EtgHvBfYopS4HTHQd9g8C9yqlLkDbus5wgKTUE2nprdDWO4XVdC0pqw9pncP3CLA98f4c4Ngc26KU+kdgt4icdeJYu23+FpAXERet8R8DPgRcF66/E/g+8IF2DK7TmUs4LlVoJvebS9NPavPVSqWl2n0q7FO6iRaFej4AXCAiu9ANfN4CvC25gYicDzwbOnyvQTdLGj3bE7dN+CuljorIHwGH0B2JvquU+q6IbArbAaKUGhaRhs2/Q7vZbQA7duxotMmqpVVCdy7B3mgSaERq1klZq7TK4auU8kTkPcA9aOvHnymlnhCRd4frPwu8AXhHqCRXgJ9XLejC1TbhLyKDaK/2LnTLvL8SkV9odn+l1B3otnjs2bOnQxOoF0+zNvbl0JKT5261Vt8M7ZxMZt/3tcJc36O1dh8WSyuTvJRSdwN3z1r22cT/Hwc+3pKTJWin2edVwIGot6qIfB3d5/SEiGwJtf4twEgbx7hitMPkEQn5ThjLfGNIBVFKpyGyiPIOHUo7Qz0PAS8WkR4REeCVwF7gLuCd4TbvBL7ZpvEtK4t10M4WgAu9b5ZcPl+3bycK2uWejDphsmsXsz//iLV8T5qlVRm+7aKdNv/7ReSvgR8BHvAw2ozTC3xVRG5FTxBvatcYl4OFNO1IG59tdpm930LvFzp/Jwr5lMa0+jNb7HclSfq90bQwyatttDXaRyn1EeAjsxbX0E8Ba5LkD65V2tds/0GjH3Cnaf9zmaSW61ztZLn9K628j2c71sWMpd2fy3wI0jKbf7tod6jnmmM+bb6ZfZfyQ04f4TuX6LNpdmJeyrFbydkK5JWc1JebVPNPaYr5vvCNfhDNaOtLoZO1KVhbpoblEoIr4SOZz0+wlM/rbJSidqCTvLpb+qfCvwNoJuKmVT+IVkTRLFckzloX/PNdayPlICmEV0pgzn5SmW8bWPhJJnm8bioLsqhmLh1KKvxXgJW0XS81eqgT6YYxtoKVsqF3GgspOJ38+acO3w5BBcFZm0c6PQqm1bbSbtCy25FothI0EnqL0frnWtapnG25kU4tIZKafbqUpWTOLrUq5lz2+8V8EWdn3y4Hiynr0OjxfzUJ6oXu8dlca6PPsNOVj6XSiu9qJ2Yhi4CRCv/2I4axqC9CKwufNXP8uTS9pQjyZoRSN2mF89GuH3cz9+9shfXs/Vb7ZLpcLEYpmiuhbWkI0uV2n1Uh/COzz0If7HKbTVZy/4WO2Y2TwEo7L5O0QiGA5pSC5GezWHNPygyLzYlp5f0UATNjtux47WBVCP+I1frovBTms+m3QrNdzsmlmz+/uTTRuUIZUwHfpQhdr/mvyjaOnWgj7FZS4bQ0litUN6VDEMEwm3t1KqtK80/SSOtKf4BLox1RN/PVNurU6I+UtYUY3a07r1rhP5u1/iPvpqeehRKJVmsIaEr3IEJHa/XN0N1TV0gU7dPKuiirjaVOfvNFSCzXvV2pz6yZe5J+f1LmQkxp6tWprGrNf3alyrWu/a8Wjflsr6E6XQZpTu9Z69+ZlMaISBrt02nMFT+91Lj6bmS+62wUxdOKe9LqSWW+Yndnw1r4/FNWAAFJa/ssHREZAD4HXA4o4JeAfcBXgJ3AQeDNSqmxpZ5jNWi6i2F2caxmMpnnmiyauXfLeX9bfez4GpvU+lNS5kYwzO7+HrV79J8EvqOUuhi4Ct3G8YPAvUqpC4B7w/ctYbVrfQuVZJiLVuRHVCvNt6RsB9XpctPbtjYTNKVTaOnnKqnNf8mISB/wcuAXAZRSDuCIyOuA68LN7gS+D3zgbM/XyYKpFTTjlG3maSC5zWJ8BNHTQ6eWJxAVoJq180+XyfUUlnlEKUthsWHby/UdlDTJ66w4DzgJfF5EHhaRz4lIAdiklBoGCP9ubOMYgTmanKugjSOqp1OEbCdrzMpoUs+Z53Pt1GtbKzRSZJayTaswTKOpV6fSTpu/BVwD/HrYzP2TLMLEIyK3AbcB7NixY8Ht53IEN0QFDe3CnfrjX4yGvpDGtNhyw41o531q2oQ1x2ecRAKvbtJY7U+P3cJC/qmV+P6JCKbduYK9Gdop/I8AR5RS94fv/xot/E+IyBal1LCIbAFGGu2slLoDuANgz549qqUjk4WrhEZmgeiLOFtQrDSrJYxzqTRT0bE2NYEyM0jgATrCAJiZBEKtPzX5rCxd+b0VkA7W6puhbaNXSh0HDovIReGiVwJPAncB7wyXvRP45nKcPzJRnGHCUQHiO4s6juGUIfAXtd9ysFxJS6tF48329seCP9b8xdD/h98BCTxqpfF4n3j7lI6gk5oQtaq2j4jcICL7ROQZETnD+iEiN4vIo+HrX0TkqlaMv91x/r8OfFFEMsBzwC3oCemrInIrcAh403KdfHayT6QdLpZM/3qciVMo00Z8B2VmWjnMRTGfw3UlGsJ0OtlCkVq5BKA/q/BpTZIKgIj+brTxSW6t0Cpz5YojrYnkERET+DTwarQ15AERuUsp9WRiswPAK5RSYyLyGrTF49qzPXdbv91KqUeAPQ1WvXK5z10X+pcw82R7+5d0vHgCEAMhFCqzbcpN2JlbyVymoKX8kFaVWSn8HJRh6Se9hOCPPv+OEzarjG7/LknrzD4vAp5RSj2njytfBl6HtoIAoJT6l8T2PwTOacWJ15xqU5uamHncB1BB62y8gYdEZgTomAkg4mw1/k7sl7CU4m+xb6ZBZE+sFKSJYF1B276LwmIcvutF5MHE+ztCnyXANuBwYt0R5tfqbwW+3fQ452HNCP86c44KyBWKLT9HZnAz1UoFozqByhS0E3i2EIkmhy4WLp3yFNBM5FJD09esex99Rql9f+WYbYLshO/TYpDFZfieUko1snDoQ51JwwAWEfm3aOH/k82eeD7WhPB3xkfAygJLN+sshiDXj1EZQ2WLM4I+4VCMzA16QfdPBh3DAvexbrJIbhuZgaL3ic8r/VyWl8Ukay1Hracl07okryPA9sT7c4BjZ5xO5Ep0KZzXKKVGW3HiNfPNFq+2IoJ/5oQG4tXq3s9en3KWLJRop4LYjNNUeYf0M2kbXedjCW3+zbwW4AHgAhHZFQa+vAUd8ThzKpEdwNeBtyul9rfqEtaE5p8ZWLkk4VgbCTWVOKIk9DMoqA8zhBUTOl33A2uGuigdo+HyWmkcTHvuJ4Pk01c6AXQF7TcTSUs6eSmlPBF5D3APYAJ/ppR6QkTeHa7/LPBfgCHgj0UEwJvHjNQ0a0L4twvxahiVcZSVIcgPzixXgZ4EZpkdUhbHbEd95NeJQm0Npxy/F9/V/1vZuZ8YGi1PP5uOo/2CP+rk1ZrvhVLqbuDuWcs+m/j/l4FfbsnJEqTCf5nxezdguJUznL+x7T+ZaJQKmZYgXm0mizf09UCD+j6z73nS1h+9T+kIOi7STAQj093is7tH3+FkiwPa9ONMg+8gZgZl2nGM+RlPAB3MbIdbx/0Y0c58b/hpPcmaNsrKIF4NZWbqJoEzJttGE28XfCZriU76nmlaY/ZpJ6nwX2Zy+Tzkd1CdLs+YIbpIsHTej25uvGP74iepINurtf/AA8OKfS9nlHaONP3Z5R5SVpRu+p4BocO3u9s4do8UWgUEmQKoBiG8y2DyWbYuWE0ubxdKDJ16H9VZCgV+3fo5d07DO1OaQ5BWRfu0jc4d2SrDcCuYpZG4nIDM1i6TGuhcLEIjbZVQTjZo6eR6/QDW1ouwN+/WRfY8R2v90T0L/Dnt+5EzGGjuc0hpCa3uI72iCBiG0dSrU+ncka02VECQ758zHLHu/VyCp0M00mRF1E6cDOzNu/U/hqXt/WLoyTZ8woozecPlyrT123gySFRnTVT8TC6LrlsCDwm8RXU868R7ttLMvgeL6rfRIXS75p/a/FcKMTBqUwTZXlAzmYGLaS/YVrpNEw58AJQdCulaSQt9MbSmH2b0muNH8AfOQc3O8AVQAUatRJBtXApE+xT0eWZXiJ2PbhBs7aJb7o2IYNjdLT67QOqsEuYQ8nFdmblizDuAXD7f0Q1Ogv33nbFMAg9xqxi1EoZT1gI88DGmThHYef1/dQJl5zEnwmz68LNIfk7R5DGzYKYQYLZQBMMEw9TLGj0lJFjLpbSTzFVpNqrx0xVPR9L9Nv/unrq6CKM8il/cpIVFxFzZpklW2AmZbJDd8T9AZgS/Mz6CONMz9nvf0c7eagnxqpiei8oWMEojWthbWZzNl2CNHwHAGjuEP6Ar5dYJ/2RIbgOyiQKBuZ4CzthxgtzcZUTSCWCGrr4PaSevlKYRYybWPCr1ENqcG2r90XbJv8vMSvdBbRV+/9YzlqmM1s5FBRjT45g7X4CMHkI5VVQQoJwK1sQx7YcJzUFGOVEvK1HqYTFmuczg5sYrOuQprt100/dqIdIG7mdJ2MnmQeCoUupnRGQd8BVgJ3AQeLNSaqx9I2wNft9mDLeinYuR4zE0Bc0p/KFjnLydStAzgF/ciOGUUXYOv7gJ+9QzBPl+HVorBnZVd+4yL9c9grxj+wAwpsfw+rdiVCbC9+PYG3fqXr+JDl9KDIxaCXztKPZOH0Rli3riAMzSCMb0GMbuF+l9vBriO/W+AjGolUtkC0Wt/c8uNNfoc44c0Cpoa3/os2UxzvBuQaT7k7w6YfTvA/Ym3n8QuFcpdQFwb/i+6zEnjxNkCrrSZxTrP0vrn9fuv5ya4wK26k7GqJUxp06Gb3SpbL+4CXNiGPvgA9jP/gvmrmvqdzpxAE4cwNj9Iqyn78OoljBC81Cw/z6s0YOI72qBrwKM6gTiOXpZNdEXIhTYfu8GAGrlUizUI8EvvlPX3CcydSTDS+NIpOTnEIUEBx4kQ1G7jGYEelfY+GcTlndo5tWptFX4i8g5wGvRdaojXgfcGf5/J/D6lR7XciBOBWvskDYxsMgs3+QE0GohHTowO9mhOx+qNq0197B2kjVxjCBTwO/fgqqUMS9/Jf7j99btY/Stw+hbB4B11fUEE6fwe9fHwlzciha+js7NMCNzkK6oSFDcRJDv109toWB2N144M6ZkPSEzUz9hJJjT5BctVwHiu/VN57t0kl6NiGE09epU2j0t/b/AbwPJWLpNSqlhAKXUsIg0rMcsIrcBtwHs2LFjucd59kyNIv0bMaol/Gyx/gdN6FiMkpIStWck8MKqlE68XUszgrvcrCSWjapMYvkuQSYP0QRw8FHmEpPB5GlgRvOxrr6B4Jkf6gigC1+qF06eBhGUaSO1MuLohu5B74bwM1DhZ2M3bteZPF9+EKOSsFzm87ru03Q5LvAnswV74McTQxD6L7rts+o6bX4RiAhGl5d3aJvwF5GfAUaUUg+JyHWL3T/sgXkHwJ49e+YKxugYzMtfqSNTsr3guygrGwr7xNATTsZI6Me23nBCEK+mncWJxKSmnyKiYyeqXmb617fmAtuECgLEztYv9B2M7RfHbyNbf4R19Q1nHqdaRvl+PCHYJ5+JS0ODfnJTuaJ+crPzjVs/JibluGprSJAfRHynrqGQUSvpdp/hZxpp/PqJwkFqZVR+BRsQLYFuTM5qFd0e7dNOzf+lwE0iciOQA/pE5H8DJ0RkS6j1bwFG2jjGlqKCQGt0gQd+gzLOyRo0CcET5QiIOtP+Pzt7daHaNRJ4M3boLtMkGyGWrW3ipo24tThjN8j1YW/a1XAf78ffRVWnsa99Pd6Pv4tkQwEWJO6rW9M+BDEQr6o/D99B/CwqEtINnLDiVuoriCbuszIzeMf2YW29KF6nxAA7r4v+JRy8+B4qW9DO5uoEyu6pP26HEPswAq8u7DVa10j7r01N1H+/u5E01HPpKKU+pJQ6Rym1E9267HtKqV9AtzB7Z7jZO4FvtmmILce8+GUAWKcPxRpjJIyj7FMgdAonHMG+q23LgY9Rm9Lbh7ZgSdavQWuTMNPIJD5G5Lh0yjoe3q2uaIez5SKYGEUFAd7IYcR3cNefhzo9jBzbh3v/N3AfuAvvoW/hP/l9/IOP4D38nXhf9/5vYF11Pf7oMH5pvO6JwLz4ZVr4ZvIY578Yc9c1GDUdUWSNHZoZgFLYJ/ZhOGUd9VOdxCyNaL9BvI3+fKJkMu/YPtyThxB3RuPXiWfhZxl4YBiI52BOnQy3a0m/2OUh8f2rE/YN/BPOxKn69V3rw5DU5r8M/D7wVRG5FTgEvKnN42kthkmQKYJbReycXhZ4iONpE4DvaDtz8ochMlOuwLR05moYvjgbZdoYvqeFfVjOOMgUdAy7UmCYSNhcJnj2X+PwxG7FeuFrAVAPf4egNE72/HVw9Q04//RlJKPvr9E/RNAzqO/nhu0Ezz9Rl2yXeembGx870tBDjN0vIgN4zz5AdvsVOKeOYNRKiFvBmB7TiWVmRvdvGDsCG3fPJHypYCYfIVE2QtxKwqTnzDy9WbomkbJzcchqJzLblFW3LmkSC4n8F0atBJ5DZv05yzq+5UKMtJlLS1BKfR/4fvj/KPDK+bbvZoJcHwDm1CmCXDH+URvVSbz1xdhso5IFyBJ+AfFDP0G0LnIAx45iNfNU4XsoQht/ZFcOPB22mGwuv4pwfvC1+H9zaDPKyqFME/Ed/MIQZukkGCb+yaNIbuEIJ//J76NqVQj8eKKJ/04cmzH9RPfcmUa8GkF5Evb+M8ZFLyEoDMXHMypjBPlBglw/Rnk0NuUpMXTzGd/DcCv4YTSRPmgHCf7Z5sK5xpbQ6J3Tx7TT2ylDwnmto6C6U/gDHa3VN0NHCP+1hL1pF97RvTpBKLRPi+8S5PowJ47hDcz8GJSZicP9gNjEE5uIwmYl4tZQll1vAoqcx4EfCjwDfF8/Wbg1FDNmqFVB4GMMbEDWbUGZOirHr4V2dMNEDj8B1b0wtAU8F2vLzqYOa1563dynzIaCrKZrB+kJ1kFqU3HNH/XI38NL3jyj3fuezksIfL1eDFS2V+8T+oT8vs0Ni8zpntATBL3rw74Fyyh8kuaY2O+go6mS0RUNfUzhk01mcDPe0b1QGMKcPI54VSQzjXHyOWD+e9vxiCBGGu2TskisbZfgH3wElS3oyB/Txpo8TpDJ10eJJBN+ovwAw4rDDsWfCRcVX9uJIao2GaCsDIZXjQU/0+P6ScCyZ0Iauxz/4CMEJ7UNPpgah6FtiOci1ZI2seUHMcujqG0XY6gAdWw/1tbzULVpyPed1bmtcy7TY3jqn7SN3nf0hG1lwShrAW+Fjs3ws/QLQ9pfU5sCM+rhLGDnkNpUHMUVk9C0jekxZPhp5NwsWLkzi861mNiprQJ9rsAHfC30Al872Wsl3MnjM2W0AfvgA5iXXod/4EeQyWNODOsVgY/KFLpb6CdJhX/KUjB3vgD/wI9QvesRFRBkCzpZ6ekf4O9+EZh2rB3GCT+eEzclMabHdPx/ZPKJ8wH8OJpHfFev8138sZOo0G+Q+Yk3tPnqW4tyqlhbztNmtMDD79OObHFriDOtY/GtDIHdg+wqEJgZzPGjiFvBvf8bKM/B6B3Auur6+Jju/d9Acj0Adctn4z/5fSST059HWDlUGRZGcR2SzSOZPJ5TRlnZmUndMHXYKCR8OTaqZ1C/T+Z4hMI3cu7L+m3Isw8iO19AYJitj5oJtXZ8D5Ur6uACt4I4lcYRYgs4bIN8P2bpJOrUEawXvpbuFpdJJFa2upVU+LcRc9c1+KePgRiYk8dRVo6gWsYaO4xf1CUDxAsbi4gBvot4VQCUlYtt/apS0nZpQPIFAtcBz0V5LkahiD8xU7BsLudmtyK+g2RyOFsuwxo/gvjaDBYLTsNAGRktmFWgbe21EirTowV14CNWAwFqGKjqNMpz8H783XknAOVUoTCkm8WHmry74Xxt0/cdzOP7UJUyMrARd9NF9Yl90UGiTmNRxJfvatNCOKkHVhY1cA4SeFiei3LKiJ1rqfCXwNNPF7USmBkC04oVDqM6iTKtGR9T3Y5GXfBArNmXx5De9QQ9g9gvvKpl4+wIVkEP31T4txmjWiLoGdQx3OUxxLIJsr06eiT6oRlWnTlAPDcM/Qurg9aqBE4VI5PDHxvByBUIpksYPUWCShnJ5AjKJbLX3dy+C10GnB98DVWexN5xIT6A72iHugr0/Qz8OFa/jkQPBSOMCAJdHlp59XV0zKEtqFoF76FvIUPbEKest4ke+QMff2IUY3Cb1tABZediLTnIFjGdCtQqqNPDsOmi2BcTZ/dGzneR0Jnv6yeWTE9irEqXg/Bdgnw/hjuNdeo5nG0tEqpRLolTxnAqgC5CGBSGMKZO6TGBNhtG5sbIzOg0TuwKyiXk2FP6zYYuyMJfDCIzJr0uJRX+bcaYHtO1aApDWG6F4PghrNoUKtOjwzYrEygrLBAmhrZn+w6qUtKhjIZJUC0TlEsEgDm4Ab80jpgmgVPFO/Is5uCGrrdPNkJVyphDWwiKGzEnj+Ot26kjSBLJV8rM6IkzUTJDWVk9OagAyRf0E4Dvaw0+vE9GcRCyBT0xmzbm1CmolvQTweA23QVs6iTe6DBGTx9BtjAjFOOTB2AYBIV1mL7WoO3R5/AGd8TnkcDT2d6zoq90opejbe2Rpu27WkCbNoHZj+G5ZA7/CHfLZUtPAEsklolT1rkMTgVVnsQEnRfie4kJy697qsJzUYaFmjx1RtKQfe2qKMvVEKF1VT1F5HAXBIIAACAASURBVAbgk4AJfE4p9fuz1l8MfB64BvhPSqk/asV5U+HfZozzX4wcfzbMTM1gn7Ob4NQRgvOvjUsJi1MJo0fCbFOninJdgnIJc3AD3rGDGMUBVLVMYBhkX3ULANXv/inm4Aa9vde9lSHnwrzwhYhbI7CyBNlefQ8zhZkGLLMjYhKZtn7vBh1a6TlawIXmNRVqc0G+H/E9gmyvnixyRe2odab1dpkefDGQvs2o0JEZmJmZ5C4xdG4B2kQnYdKdjB+H/m0YtQk9sUSlo8M6QRAgge73LE55pq6TUnFJEGVlwPe0Pd13sJ5/CKwM7o5rFo4AStrtQ19C1OfYmjhOMHEqdFTbqLBsRpx0GDp5I2e5Nlu5qMlT85xwlSK0RKEKS9p/Gng1cAR4QETuUko9mdjsNPBeWlzksrs9FqsEUQFGaQQsG6ysttVPj82UKgatcVo5/cM1TJTn4I8O4x7aT8+bfpvcDbehfL/Ovp+7/tZ4IsAwqfzdp1f4ylpH5RufqHvvjhzEmB4HS5dyIOyTEAnTpkIhDQtlhk3e7Zw214RN3zEzunKnlY1t8VGfX5Xp0SamfP+MkI8mltkaeCggVKagI7R6+rGOPKrXRWWbkwI5qtskBkT2/DB3I+4ZHCbxRdeAldHmvllZ3XNedlhKQiKBHqKyBZQTmhALOhJK3Jp+4nCq+h77rvYv1coYlYk4sixZGmNtIHE474Kv+XkR8IxS6jmllAN8GV3ZOEYpNaKUegBoqQaXav4dgDJtgqKOUAly/agTz+sMXsOE8lj4JeoLs0h1mQcxTCRXIPtT74iP0/OG9zc+gWFiFgdX4lJWBDcM7fQHthH0DCJOmSDXrzXzMCGuGTNIYOchU0DsHu1cDU0vdU1YEs5ZDGumIJ9pa0ds9ONWgTbHGRbKzoc5F9qda9Sm4mMr0yYoTyL9ZcyoUqidBzs3E/sfTWCZQl24L0ppx2s4aYnv4q3boQMGAh/vwW9hX3ANys5p01IDrPEjMPwMxtYL8Ysb4xaV1umD+unHsjH7N+snDqcS103CqWqHeWkc5bkoz8Ec3EgweZpgenLVBRI0wyLMPutF5MHE+zvCwpQA24DDiXVHgGtbMLwFSYV/B2Bv3Ik3/LTWVn0X2f1C1NGndM356GkAtFboVMHKYBQH4kzThTAHN2Js3IF//MCyXsdykn/97fULlELlilprz+tkOEKneNNdrxLaugTeTOnkBtvoDUOfQZh8F5lClJ3HOn0Qv2/LzKZhrf9kEb3IZGT/m5vwDz8WCnE7PrbUpnRZh6gRTNijIArljUw+QX4wruFE4Oks5r4N2LleglNHkKFtGA9+E/OcC/GGdmrl4V+/ibXrcrwDj2s/x3MPY1gZ/Ct/WkcleY5+otl8HoFh6dBVw4JKCaiifB+qZfyJUcTWTyTKc3EP7yd/03ubu9+ribAER5OcUkrtmetIDZatSJXi1OzTIQS5+oqIks2hMnldDyYs9BVkCyjfR+wsKlOINeCFsK99Pd6Bx1HlyeUY+orijB3XppqojHMyo1kFcd2ixRL3U5iHuvViaG3etMnl8zMRW0kMc8ZWTuhPCHyciVOY268Iw3XVTDKVlZlxTIfXFnXzUpkCBAFBfnBmvFEQQODpiUcFSKEPNTGCufMy/GPPII/9H3jw72Kfj3XB1VgbttUNMygM4R/Zp+s/GdZMzkjgaUd4LfQzhU5xyeRQtQruof26qupaJAz1bOa1AEeA7Yn35wDHlm3cCVLNv0PIDG6mOl3GqE5g1MpItkdr6tt1FilBAJaF2nwB/hIcTdbmHbhHntVZl3Bma8MuQfs9RGvgs2LO4/4HSyx7sFDMfF19f9+JTTTVSgXCAm6RGSVpT48qtwa5/roJRKJonnCbyK9gOLpcRHS+6ElAhTWCov2CfL9ODAwzcJWVRUD3N6iVMQc36tDhngGc+7+Nu/d+rE07IK+fcPyJUeznH8QfGwkDCCahPKnDhDM5fKeK0VPUT59Vl2C6RFCeJJgYJfDcOEz21Kd+C4D17/3vS7rv3UnLkrweAC4QkV3AUXSF47e14sALkQr/TkOM+CnAVAGqWkJVyzrGevwEZr6gtUB0lnCzmJdeN5Ny36XUpibO1K7niOFvNbObs8w1UcQTUmJS0pE7FZQY+MVNdfVwlKmPOV957WolkV2rAgIri+G7YZG/IA4ECHJ92r9TGgnLTOQIxo4j+T5yL3g5iIE/egwpDGIWhwj2/ivO/ofDkNnNuIefxgizmr2JUW3SSjwtBtMlaiOn4mqWYhhURxu3qFz1tCjaRynlich7gHvQoZ5/ppR6QkTeHa7/rIhsBh4E+oBARH4DuFQpdVaP8qnw7yByPQXc8ijilHX9n+enoTIJxfXg1XTRMi9hR14C3arxt4wG2alN7TZHAbNo/1w+Hzc2MdxKHLETafvKzNQlhlWny0jvBp3HkcnjnjyEsjJkBjfPOw7xakjkJ/BdxKnoonIqQNlZXcnVtGcK+YWTSpAt6AKATpWoSqyankS5LpIvxFFiKtRmIzOPAdoP4VRRroOZy8TCX4URPqa9Fk0/rSvsppS6G7h71rLPJv4/zjKUP02Ff4dhb9iBc1qb/MxzryLYfx/B+AmCTbu1fXfkSaQ40LTWUfs/n4/ttJGDyjz/BXFRsm6hWqnMhD6eDc0I+0Y1bBY4VrVSiZ8OlJXVgjmM3ZeweF9cqz/xBGE4U/h2Ns4vaEQun4/7/Ubji3IZJPCwTj+vl5s2Xv9W1OlhnfdhZuIkNQDTO6XNhyoATP19cF2U6yCmGdv2o6xwfaoAVRqDyGfQk0fyhdjJ637mg1THS2z5UPeGES+ZtLbP0hCR7cAXgM1AgA5/+qSIrAO+AuwEDgJvVkqNzXWc1Uhm3db4f5XtRey8LvxlWhjrtiwqFCCYGgd0xqo5uAE1ML9m2bEkBHJcPXOBbc/YrkmhHmvqdQsXTp7KFopUK5XY76AkH0fs1EUgRU5d0yawe/SkFib5zTcm8WqIWyXI9SG+q58a3Jq265cnUVsuBCA490qkdJKoeU/QM4hZGiHI9xOUxpB8f1w4MHCqGJaNMz5OdtsO7cCNylFXyqjyJEEo+APHo++W360bV89Fl9Kz4B1dhYjRuCZUF9HOqcsDfkspdQnwYuDfi8ilwAeBe5VSFwD3hu+XnZOfuJ2Tn7h94Q1XGPFqM5EtYuD3b8Y/fgDv+b1U7/5M3bb+gR/hHd1L7XtfoPa9LwA6RFL5PpmXvQVVq2LUyhD2kg2ee3D26TqXeYSv+M7Mmyh5qdE+8x1jCRFC85I8l2HWJXHVlWvI5OPJYb4xiFtBqiXtIwh9Ct7QTiRM/lObL8AcP6q10fhcOnfBcMp1/YyN8WM6qa2mBb9fq2H3FrTjNwhQri5FIVldPsSwbJQfNOxZm/2pd9TlmqwZBH2vm3l1KG3T/JVSw8Bw+H9JRPaiEx5eB1wXbnYnusPXB9owxI7BHzgH6+SzOqHp2L7YGdcIdexpbddNrI+Sv6K8gGD4acTVVUDdB+6i8viDuGVtrx56zx8u56UsnYTWPltDrnsfmVZU4wbrTZ1qKeYlMaiVS3HxtSiG/wyNP/E3jtAJxzrfE11mcDPuyEGdcQsouwfz9CH99JDvR5k23voZAR+Vh7C2XBAvszftQq5+Fcp3MEsnMK/4CYKJUWxAeQ6Vpx5FjBEy51+J8lzsa19PsP8+MO2ub/fZagRJq3q2AhHZCVwN3A9sCicGlFLDItIwDEJEbgNuA9ix4+wrBm64/ROMfWZFHjIWRxi77h/Zj9G3DrFtnOEj+FWH/K7z6jfN5AnKJdwR7TPINToe1AkEf/Su5Rp5a0nUo5m3deASWx/GUTpLjRYKawnlCkWc8RFdg2l2Y5b5zs2s5ucNsDfuBMA7uhfl1qBa0oludjbuJQDoUFjfQ+XO7K2rjuzF2LCD4OQhgrKOJDMKfbGN3955CX5Jmwpr3/8iweQoynXoSYV/PS2K9mknbRf+ItILfA34DaXUpEijhLczCdOj7wDYs2dPSzLi7L7Os14a578YTh/DOO8qXUPmqfvoffnPMPbtv8IdOYb/lx/DGZ8iu3E92Wt+KqxQOX8kUOWuTwGQueoV+KPH6X3FTYzd/RXyG9ZR+sJHGXnoKXLr+tn2kT9ZiUtcVhb0D0ScjeA/m+Mlqo02TbiPv3F3WPU1nOaj307go7K9ZIsDZ+66ZTe4NcQwdSBArkBQnsSwbDIbNp5hpsicf2VdvajZHPnIu8hvGMDKZ+m/9WPNX0PXI10v/NtqkBIRGy34v6iU+nq4+ISIbAnXbwFGVmo8vTf/l5U61aLIrNuKMjMY1RLG9ovx1p9H3xvfjbVuA5WTM75wd9+DBFPjqCA4wzE3F9YVLyPIFem9/AVIvoA1sG65LmP5mMvGT8KEkwzxbBQmO5/wbbR9MpY/cX7v2D4yAxvjTNwFxx3uu5DWnyTIFXWlUbtHJ4mFmcTRE4RRnWwo+CFsPRmVk7ZszOIAmfMuwxrajPWi12IObcU+ZzdimmSvu5naYz9YdZ3fWoIIYtlNvTqVdkb7CPCnwF6l1P9IrLoLeCfw++Hfb7ZheB1H1CPVOX0MY+okft8WvNMnqYyMI6ZBZqCX8vAJVBDglquc+cA/Q6NaLNbWi+L/d3dqrZb52ggupD0nHa2tYFapaHErmJPD8491DnI9DWoKzUPUAAjA790QZwdHju7IPDQnkXko8JGhbfgnDmJs3oXK5PHy/RilE7pC7JPfp+fnPzTvoTb/yvup/sNXmDwwPO93blWyTAmFK0U7zT4vBd4OPCYij4TL/iNa6H9VRG4FDgFvatP4OhZl5RCvir15O1te+fP4T/8I9+izeNNVerdvomdLhkMfuoWeLetwJqfZ+uHPLHzQbmG2UG3GbBKFfkY1/mcfJ4oQirTn2cXcIK7YaY0fIegZrCsCFzVkiQRyPK6FroPFC/742J6jM73DQnPR+MRvouqvYSJD2zANE2VnMYoDYWSKpR3PmQLWxm0LHwcQp0L56EnswlweptWKpMJ/qSil/pnGFe0AXrmSY+kmkjkAzp7XYK3bqjNJp8Zx9j7P+P7DceZlxKlP/dbqqbsyW3jP1uaTy2cJ8YZZuuExVXKfpBkpihwKzUfe4I76+jzJ8xgWZukEXmUCGTinvqdAg6eOJQl+wOvfijVxDN9YH3fWiia3bN/CZrsooUyyOXAqBOu26yglt4qyshjT45y+55v0bB7Cv/9785oQncfvw6s47XcetoE5v09dwlr8zFYN0URgbb0I/6n76d+9jdKhEyg/YODC7dTGS/hVh6nSNOvbPNazRdwKRm0qTHByZmraqACjNoWyMnWOz7g3bqMAAqXibeLjR3X7o/9VgFRLYFozxw2PTeAjUdOVqHibBLoPgOfoCBwxUIYZ9xcAdON4M9OUgJ6PXE8B77SD4dVmEt/cyrz1gZJY2y7R/2y7BOfUEcSrEWR1sxlj6hTeuh3YhTzedHXBY2XOvxLryaeWfC1di9D1mn93jz4lxj5nN4Vrf4p1V19GcccmrL4+3HIVMQ28cuMG292EUZ3UQter6hwF3wMRHe1izrRCjM000f9KzbwCP65rAzMCPxb8Ufy9GGF2rKHPE00GoYadDCedXWdJZXpmWjNG51IzzdzPVvDH9yNqmh5WDzVLJ/B+/N0lHMgIq4564DlxZzExDcQwFm5Yku+jMqrri534g19f/Pm7FtGKQDOvDqUp4S8iF4rIvSLyePj+ShH58PIOLWUxGBe+FPPS68jf9F56b7wZs38Iv+rgVR2swvyRJM7pYzgTnd2HNarTL56jtfxQyAZ57WaMG6wYZthzVsWCOZkBHAvrwA8bknsNTTKI1DVkr5sgDLPuRy1eVdvaw/Mb7qzJVgxy+fySzTyNcDaH2rtSiDONOj08/w4hUfnlCGvsCIihS4lXw6qpvkNm0xb8qi7kNh+H//h/YBgG5eFR3HJ1zUwACsIWoAu/OpVmNf//BXyIsIekUupRdN3plA7EOucyzEtezLrLdmHlMlgL/IC7CWWFdXDC0s7KympzRW0KcauYJR0ZbFTGEc/BqE3pHaN6NVY2zqqNHbyBrztvhdm5hlPGKJ/WRdHyiRiWqLFK1OcXELeqs24DT08AUR2fRLXPZBinM3YcZ+z4Wd+HXE8B69QBxHcQr4Y/ehzrqusX3G+278cbOYxRGsHcfoV+mlAB5sizGD1F3LJ25s7Hzo/fiZnLkBtaY7E+YU+Jpl4dSrMj61FK/eusZS0uhpLSSsRzGHloH16o/TfCmTilK4iGgizYf59O509uc99XOfKRdy37eBfC2naJjqjJFus0ckC3RwwFuN+3WZdNsHMzTwhhw5TYPh5p+qYd99YVTzcqF68aRtLkMcqjmBPH4qYtOgtY1T1J6IbueW0C8hrbyJ3Tx6iVxuOSz60mnuCWgCpPxl2+/JNHMCePo4rrwTBxy1Xc8tx2/8MfvpXDH741TgY0bIva2NLH0nV0ufBv9pnklIjshrCVrMgbCevypHQm5rlX6YJuxfmzlsWtIbUyfnEjQUHbo71TR8C0kGdmz/ftxd6wA3+6DIGPPfJ0LHQj0447dN6MeWZWP944MiNRV8dwK7p5eohRK+mm5YaJVEt4/TORVRKVQUaboAj8ugJy5sSwbqYihm7gHmYWVysVDEg8qbQuJFICT3fyCvxYgC8W7+RRvJNHKVx6Hfa1r6dy16ewtuxEVcsMXLidwlubs+4qP0AFQVzrf/WzdkI9/z26lMLFInIUOAD8wrKNKqUlnPuu2zh85+fxyhWefd9byK3rxw0jOM55y1sxBzejsgWCbAFr7BBKDILiJty7P0vuxTfiDB/E6CmSHeht74UksCaO6Z7GpqkFNYAIfnHDjDBuEOYp4d8ZZ6w2HRnVCVS2iFQmwlh3E5Xp0UX04vr5Sv8/q4cvXlWbgQIDb9259cXnEpFEQa4/jsNvJcZ5e/BPHoIcWFe8nPE7/iMDt/3eoo5hbdjGiXv/kcJb9fv8Te9l6ou/i/IDiu/46Lz7Kj9gx//zeQDENAhcb8HSIquJbg/1bGr0SqnnlFKvAjYAFyulflIpdXBZR5Zy1piXXkdtrNSwFC+A37sefGemPkwYvmjvuFDXgB/azMl/+hcAjv7Or6zk0OfE2nKBFsS+D2bYIN0wEC905EZJTrNaKdb1A5g1SSgxUNleXS7BtPT7aFuvdmYph6QTOfARt4YxPabNT8lqomFRvsgPEPkZWpppjH6qmPjaHUvad/SH+ukuWc68fPQk08fnrucD2uQTJPJJCluG6Nu1BYBDH7plSWPpOlaz2UdEfnOO5QDMKsuQ0oFcdMfXY60fwLQtMn09BDtfgOz/oY5rt2yMgQ34UVvBS34SZVgYO69g6Oqj1EZO4kxOc+xjv9oR2cJSK+NuuQx7+AktWH2DID+AOGXMyZG60saETlhr/DDeup0z4Z1OWbdALKwLHb0OQaaAMvt105RIQHsOEv1KVBBHGwHa9yAGRnUSv18LPsfKkwmfEJRhab9C1P83DKPM9Lcu68Ko6h66PTt3MvbI4zSu6DM32YHeM7JzN/32/1xwv+0f+9O692P7D2P35PDdNeIKlNVf2K0YvvYAv4qut78NeDdw6fIOLaVV7P7kl8n09eBOV8kO9LLuFa/EeO4hXa1x1wsIzr8Wb2intnWHmmqURGXvuoziq98YPz0ceP/b23sxoKtQPn4vAKo2jXhVrPHDmMf2huYZv17LF0GZGcyxQ/ExxHMI8v2YpRHt5LTz2uavAszyaP2TQq2MVEsYtTIq06OXRwLezmrtv1qCwMOOIkATGn6cPRz4Ldf6ze1XYG6/Au/0SQpbhhYdalkbnyJwPYJFCO0jH3nXGd8D5QdUT0/glSurIq+kGaInxIVencq8I1NK/Y5S6neA9cA1SqnfUkr9FvBClqGhcMryseH2T7Dz43ey8f2f5PCXv8rUQ/fpBt6hRmpUdJw3KtAmjP0/1JEu6Gqh08dHqY3rSI4D7387B97/dp7+1Te25Vqsq67HnxjF2f+wdnZaOdTpYSSTQ4mBWTqBdfJZzPEjeofAx+/bjD+gv7LG1CmM8mmsiaO64YqVxZg6hdTKGFM638EaO4R4tZmJAEIzzkz0S5DrA8NCfAezPIr4HmbpxIxDONL4JUykcqfr/QYt5Pl7HkD5AdmB4qL22/rhz2DYFm65uujP8+AH3hn/v+uP/oKezUOL2r+7ka7v5NXsyHYAyXhBB91jN6ULyQ4UyfQXOf7PP8KaOIZx4EcYpRFtl3Yq+Ht/iHtwL6e/83Um/unvkUyO6lgZZ7J8xrGeu/1tbbgCyF1/K+6pE3injuMf3Y+qTochntmwcfkAKlsk07+ebN867GOPY518VrdOzBUJegbAczGmxzGqEwQ9g2GuwDRE9vzYn1CNwz2j0FCVKehlZgZv4ByCydOIW9HhpYZFkJkVZaUUGFZdbaZWcvkXv8XAbb/XlMP32fe9heduf1vcvGjD7Z/AqzoUtm1o6lzV0QmsXIZMsYdn3zeT7rPlQ59m23Uv5ILP/PXSLqKbiMo7dLHNv9mR/QXwryLyURH5CLrj1heWb1gpy4ldyFE+ppOhlJlBeQ7uof1I6ZTWWq0M0ydOofwAM2OjnCrZvjxuuUZh8zqsXIbq6MSC2Z/LTfEdH8U/eRR/9DiS6wHLjgutKSuDvWGmw5tx4Ut1QlSoeUcJWXVmGMNAPFfb8oPQbOO5oUN3fMbJG5eNCEs+2DnU1ot0nX1bC31RisDKokRiE1Cro31awVO//HoAamMlnNL0gpP5c7e/DcO28KpOXS+JiIVKQK8euj/Jq6lvo1Lqv4nIt4GXhYtuUUo9vHzDSllOqqOTHPzePna9+hKO3/kZ+s/XBeGKJw6j/ACvXI0bdme3bsPasI2BC7dTPn4ar+rglMpk+nrwqw61yQqPvfVGrvjLu9tyLT1v+u0zls3lhjPOfzEZdJatTBxHLBtVLaMGt2JOnYQoVt60wfcxprVwE7eiE8XEAN/BLE1rB29k0xUDbJ1p3LB1Y2j2abbw2nKz+5Nf5sD73870yBh2Icdzt78NFQQ6Ezwxoe+95SYu+fxdde8NW4sMK5/BqzhcdMfXzzj+WqETJ/PF0NToRWQHcAr4m+QypdShufdK6VTENNj+st1seMkeRh94hMkDw9rx53gYGQvDNAj8AHdymuKll2MMbKDv6j3kjz7L5IFhMsUCtfESXtXBLmTbfTmLJjO4GefJfwTAyBWQiRMAqP5NeoPjzyH5AtSqqMAn8FxU4GP0FDHyBZRThb5wW9NGfP204Efdu2Y5dY1aaUWuazHs+qO/YO8tN8U5HI3i8xsJfr/qkhsq4parZPpaV6uo65DuT/JqdvTfAv4ufN0LPAd8e7kGBSAiN4jIPhF5RkQ6sLN6d/LMe96MmAaFzes4es//xSlN45SmUb5PbbxEZWQM3/FwJ6cZuHA7GCZ+32aszTvI7LiQTF8Phm1h2BZWLoNbruGWa+y95aZ2X9qi8MdOanNRvoByHVSlrLOdx7ST2xs5ij8xSlAt678To3jDB3EP7Z9pg+hW4yqhkS8giSilK3uK0TFafxLf0eWoA9dDTIPq6ATV0Ym6bY7+zq/w4ze/Rm9fdclvHKA2PrWg4G9XMMCK0qKqngvJOtF8Klz/qIhc04rhN2v2uWLWYK4Bli3rR0RM4NPAq4EjwAMicpdS6snlOuda4MQf/DqGbWEXcvSeu43y8dMo30f5WggMXriD0SeeY3z/EYJQE9y45RykVtJRNbUKmWIhrvdiZCzMTHfGOh///g9Zd8m5SGY/RnGAYGIUKY0jpklQ1uWjlVPFm5zALPQidkYvM0ycZx4lm+2ZMftYmTOzfxO0Mq6/lUyPVlCBwi5kCBwfD6dhgbZMwSY70Mv0yAR+tUZ2oDf+zjz9q2+MvytWXpuMdn/yy/iud4bZaHXRGs2/SVn3GuCC8HUt8Jnw71mxJKOVUupHIvJvzvbk8/Ai4Bml1HMAIvJl4HVAKvxbgO94nLj/MZzJMsUdm2Jt78RDT2HlMgR+gF3IUhsv4Y4cg0f/EckVcE8cxi1XUH5Apq9AdXQCw7bIDi4uvLBTqI2VCJznsMIkJ+UPY/X2Mj2sG+IUdmzD3riViUcfw8xlsAt5crsuxNq4De/4QZTrYu20UNkCRvk0tufgbryw7hyLacy+0hiGkCnmcMs1jDCPozZWb6KaPDiM8hXV0RJW3qYyOkV+qDe2/YN+cjBsCzOXZefH7wTAq7hc/sVvrdzFtIEWxfA3I+teB3xBKaWAH4rIgIhsUUqdVX21Zuv5/2bi9X4R+RIwf63Xs2MbcDjx/ki4LDmm20TkQRF58OTJ5RxK9+D84GvzrnfLVaZHJpg8OEx5+DQAXli5UUwTt6yjWay8jZXPhI1gqgTlEv7YCF6o8ScTgrKDRfxqDb9a4+AH3tkRSWALcfjDt+JXa9TGp3BKOoRV+QFuuYIzNk5tbIrKyDje5CRBaTzWbAGMQh/ukWcB8E8cwh9+FqmVcZ57Ii4FIcsUy99q9tx9L9MjU/hugFv15iz9nR0sYOZs3LKDV/Fwy1WMjM15n/gSPZuHuOTzd2kHcCK5y8prx/fjN7+Wh19/Pftu+zkAHrxxFXVobT7aZ30kq8LXbYmjLCjrmtxm0TSr+SdVOw/tA5hf0pwdjQxlqu6NUnegi82xZ88e1WD7tccFcz8JBs/8EMO2WHfJDiYPDpNfX0T5AZXRCQLXozIyTm6gh+pYGd/VkR9mLsvEwWGy46XYIeiUpuuO606WMTIWYhjxNgc/8E4Cx+W8T3xp+a71LJgeGSM7WIxt3HZfgcLmdSg/oDY+hQoCnFKZqaMnOb33IXrPWU/geJi2zeQD92EX8qijOtah9sgD2M89gfID7B0XYkyPoXL659LJWn8Sr6Inq0xR2/Gj2P3dJdv58AAAIABJREFUn/wy1TE94ZdPlDEzJnZoAvLKFR6/+bXseNULefIdP0t2sICYRmzqsQs5HnnDTzc0Cz721hsB2hYh1goUQjBnC/IzOKWU2jPHugVlXZPbLJpmn1uejLJ9lVL/TSn1ReBnz/bk83AE2J54fw5wbBnPtyrIrJ876VrZeQYu2BGGZ9Zwy1V8x9NleG2L/MYB3LJ2Wmb7sohpoHwf07bwqw5+GNdt5TIEjkth8zoyxQJWIR+3+ov2CZyllRdeKexCji2vuZ7sQJHADyhu30jgemQHiwxeeQlTR0+SKRbwylX6dsw4as1cBqc0Tc8LXsKpR5/h1KPPMHFwmOnhUezN2wnKk5iVMDzUb9xDodO4+hvf5YV/+/dk+zJMHhrh+MPPn7HeymXoWd/Dtp+8BBUopkcmOP//+yoAJx/Zz/Spacae1vWfZuM7Pj0beykdCWsQrZ9Jfoscyd2JIlDNvRagGVm3LPKwWeHfKHNjObM5HgAuEJFdIpJBdw1brZ6jFUG8Gkb/kBZsBRsrl8GvugSOj2EaVEbG9aN/3g4FvKebvx89iRgGblk7eyN7f2T6iUr5AvjVGmKa8csdOUj17s9Q+94XqH2vc3ICz/vElxj+9ky/20xfgcK2DdqfUSnTu20DXrVGdrAXq5CjdOgkTmmaqbCrlT+qw10zxQLZgV56d23HuuwlmIMbdRN3INvbXZ2t3LJLdaxKz1Aet1zFzGXjhK+Ney5m6PKdOKVpTNuMtXmv6jE9Mknv1iK9W4va/JPLsPeWm/T3pJAhU8jglKrYvRkev/m1mDlbO5h9LRT33nLTnJPAY2+9MX5K6ERUk68FaEbW3QW8I4z6eTEwcbb2flhA+IvIa0TkfwLbwlCj6PXnLGMnL6WUB7wHuAfYC3xVKfXEcp1vLeBuOB/jRT9D77YNoYYeIKZ+mnTLNTJ9eQobtSnIsC3ENLAK+TC22yHTV8B3PfywZnttbCo+Tn7DIGJqgRBFgfRsHOjoOOidH78Tr1rDzFgc/8HjiGFgh72Oizs20b97G0GY6FbYPIBfrenM6KMncU8cZuDC7QxddRGDV11O7oqfwPnht3S2tJ3TWdJdhlN2sXszZPvy1CZrTBwYic12/bd+jNHHD2IXcgR+QPmELvOhfIUYQm2yRnVs7msOfIU75cQmJK/iYtoG2b5s7Dh+/ObX8uQ7GhsTIn9BJ6GAQDX3mvc4c8g6EXm3iLw73OxudHj9M+iWur/WimtYyOZ/DHgQuAl4KLG8BNzecI8WoZS6G33RKS1AVIB5fD+5oT5Kh0YQw8CwLZQZYIY/wEjwB66H7+hXdrCI73paEx4o6nWuh1+t6fIPOZ3kFbja1GOHE0ZuqB9lWkgmhzm4EbX1IpzTx5atts1S2Pnrv8mBT/6RzlsoDqJKR6mdnsCr1hDTwDANamNTVEdLBH7AxPcep7CpwOSBYQLXY92VF2HvvAQMk+wle/AOPI61eRLz0uvafWmL5tq//wd+8FOvQPkBucEzfRWXfuFvAdiUWHb1N77Lvtt+jlzGojKqi/7NV9fnwRtfyfTIFEbGpLCxGO9T2FjEmTxz8nCmXHq3dG4kmVrYpNPscc6QdUqpzyb+V+iGWi1lXuGvlPox8GMR+WI4Q6V0KcrKQr5PC7qMiWEauOUadiFL4AcEVSfWwsxcFjOn7f5+tYaZy4bRMSUC19PrDAMzl8UwZjo49Wwewk/0CxbfQ1kZjEIRP0x2qU2eRnyHzODmttyHJMZ5e+KsZrFsfMclcD28qkOm2EPg67aEYgrupEv/uTqD18xY+tU/pEtEOFX8sRGMQh+qNnfP205n4Nx+qmPV2IfTDHYhh191yA8t3O3NzJgYGZPA8XEmK1g5C9/xqY6VUYEi8BUP3vhKcoM53LKLlbdQQcDFn/sGj7zhp3nB1+45m8trKZHm380s1Mzlq0qpNwMPi8gZl6qUunLZRpbSUpRhIYFH6bAu6Oa7no7n7yswcWCEwsYZDcuv1mLTj92TxymVZ2q6FPJ45QpWaCLxqjX8qoOYBut/7T9T+tInsAt58jfeAoGPtety/FyfboASlj/uJHo2D+FOV5FcgdNPHKB/9zbsQo6Rh/ax8YUXxaVshy7ZpidD18PqydGzYztYGV3zJ9+Ht/dBnQgGWFff0L4LOgsu+fxdPH7za+O4/2fe8+ZYI+/Zsi6O4U/SbERXZNe/6qvf5rG33kjgB5imVkLMnE3pyIQW+ol9sn1ZLv7cN+L3URZ5RySOKfBXs/AH3hf+/ZnlHkjK8pLL53H7NtO/cwvF7RsZf/owXsWlMjJGbkAL8kjD96oOFv9/e2ceZ8dV3fnvqaq39evXi3rTvlqSLdvyIiFv2HgFYwhWSDAmzMQQBk1YA5mAITBJYEIwkAkDAyHxBBKTj/ECdmzHjjFesB2DN3mTN1m7tUvdrV5fv62qzvxxq1stqVtqWd39Xnff7+dTn65XVe/Wqfden7p17rm/Y9I6B1d5mnnFhWy784FoQLRIoaMHLxUnKPnEk1U4fR0Ue/qoff9/Q/LdFJ55AG/2ItM77ukEIDY3mgRVAT1/MPH97N528Iv40U2s2N2HBorfl6fQ2UtQDAhKPom0mdwUr0njNg1Ks/ZL9G7dQaq5vnwXMkqcdvN9PPc7V1Azu/aIUMzmP7mWRd+79S21O3iM6fRb/oM31rwfDUOCYoncnh68lDegEzXU5DCtQE87WmGfcnGssE//iPInVfX6wftE5FvA9Ue+y1KpxJrmkl5+NsVt63GTCWLp4oA6o7jOwCCvny8SFH2SDbUUOnqIpZNMf/vZaBgSz1TRuWk3XipmBu7iHlVN9cz648+BGkXQnrv/GTcZJ33xavzNLxH0duKdtxrd9CyA0c4v82fRT3LJaYSlF+l+eR1zLl9JqbsvEi1LsXftJhI1CWact4x8ezexdJJ8e7cZ6M5lceubCXZvQnNmAFRch7A48aOjK/79QV78vXeRrE/h50pUNdcO3Og3ffoaUk31BPkibjJO67ptAJx5xwP03fZNtj/4NMBAj/2VD79nYH8/L13zbhI1CTQI0VBJt6QpZYvmZlAKhwzxnHWXyc565spLeemadxNPx8r6BKDARC9VP9L/wSuG2DaRk3SnLF7UYy1FhVn69nfTt7+XQkeWQkcPhY4eEnUZEvUm86c/jXPvE8/jJKuY+e7LWPLJ64ilk3ipGAs+/RmKPVmT4RIEAxPHvKoUhWfMP7CWSmg8hTdjIUH7XnTH6+W5+CFw5p9O+qxzKXb34cSTpBYtNpOY8kXi6RhhMaBz4w5i6SSp5nqqpk+L3ugSdOxHSyX6Nr5BKZvHq28i1tRy9BNOEM684wFO/qe7OO3m+yhl88TSSRZ89k+pP3nekMf31wWYe8U5hyiEnnbzfUf05NPNmSi7Kkk8k0Qch3gmSTxTNfDe/pvG4bgxl3jazB5+5spLT/g6TwTVkS2VyrFSPT8hIi8DSyM1uf5lK7BufEy0jCZ+y1Kqzr2S7p1dFHvy9OzuxUt59LX14aXiVDXXk283tVi7t7dRzJaYfv6ZpJrr2X7HfWy95W7c+mZybT34uRIbvvUd/HyRDTfcgPh5wjCk5ff/gNAvIekavEVnkDhlBW73PsJEGvFixM5ZXe6PYYCgbjYyczHVs5vo27OPnfc8QGpGC/NWX0F6egNuMkZyWq0J9dQ24LguTrKKoH0PYXc7bY89jrgOVc31hLksWsiR/+WN5b6sUWXpjXey8Ls/w1m0ioZPf4dZf/mPzP3mPxPki9QuaB64UQCUerPH7JGf9IPbWXrjnQMTxTQMB86R68gPzCMY6gaw4t8fpNBdoGNLJ6n65BH7x5PRSPUsJ8eK+f8MI938TWCw1GiPqh4YM6ssY0YYSxEm0lQ1psh35PFSJuMiVm0GK/v2d5jc/qJPLOmRbMiw8ee/pmn5fBJ1GQrb95uZrPH+vH5F3IN68KVu4wD9bB6v3jFF4atqUccBx8Opr7CecRig8TTxpmbC4m6Coo83fS5hNEbhROGwQmcPydlmPcz3ERZ9YvX1A9lNfr5AMpFCvCGKuUwVvBixurrjesvhxWAyM6opdBdJ1ScP0VQaCnFGLK8w6qhCUMnd+hEgxzNoISLNwMDttlKKuaxcuVLXrl1bbjMmDPlcjtjm3/LyX36bQneRqsYUfs7HiR98pJ5/9SVs+NmviCVNaui81Vew9fZfIq6Q78iTrE8SFIMBjSCj75/AzxdoXH4SAFXv/BAaT6PimLKJQKxlQdmuezj83W/Qe+9NdG/dw+yP/jc0nkZCn+6H7qR76x5SzXU4MY+at11A97O/MU8BiQSl3iylbI6qGS2ExTxh0Sd98WrITcxc/9Gm8OjN+Hu24tY345y7esiaBn23fZOO9duY9Zf/eNztv/Lh97wl5VARee4oWjsjYvmZZ+u/P/z4iI6d35g54fONBSNV9fwdEdkIbAUeA7YxxsVcLGOHk+vAn7eC6hm1BMWAzm1d9O4zypZ9bTlmnHcam257EDfmUOguUOgu4FSbHl2+I48Td9FASdSk8HOm5+vEYwQln3kfX0Osvp6q5ecQvP4UaIhTyiGhX7lql/u20rN9H32tHaYwOxBkmkmvXoOXTpqU13SKwqZ1pGa0EGtswZu5gOS8ReTbu+nbYyqBOXGPwvOPQKZhyNP03fG39N3xt+N2WZVAfOFpxGYvAsz1t//gC7T/4AsD+5MrL6PxwgsItr143G2XUzLa5PmPirZP2RjpgO9fA+cCG1R1AXAZ8Jsxs8oypmisCvUS1C2Zw7xLTzEqnimPXEcexxWj9BkobtylekYNC997DtrXDYCX8nBjDkEpoK89S74zdzDtMxln17/ehFNdR9+6p5F0DVLKIfkenGw7Tra9zFd+JMH6/yTs7USDqIat65rqS9GTStPHvzigbeRm6kzIJ9uDv3srYb4PcQ5m+DhejFJnJz3331LOS6oY5IzLkGQV2WcfG7berRb68BZOzOlCo6TtUzZG6vxLqtoOOCLiqOqvgTPH0C7LGJLI1OH2tlL3/j+i6Z3vZMnqM2lff4CwGBIGamQbYg7xmir8fJE3H3iWrT+/n1LeRxwxcXBHiKdjxJKmnOPCr3+H6Z/5Kn6+gBbzaBjizlhkql11txFmu3FOOrfcl34EWjud0s5NhGHIgj/7c0qvP4O6ninirgqOS80HP0MsnaLU3oabziCJpCnl6BdJNtQcUszGiXmk3/exI86z/csfxcnUk92+azwvb9zp+NGXKHQfoNBthgQlkcKrSiGhj/uuj5Nqmkb1Hw7ShKyqI0xPI0xPK5PFb52JPuA7UuffKSLVwOPAzSLyPcZQ2M0y9oSpWjSRQRaeTeM73sGVz95Osj5JobvAgde303CaGeBN1FUTSyfoa8sRS3rMufQsFn/qY7hxlyV/8VfE0omDqX0a4iUTSDJN+uy3EyYzhPE0/r4d+Pt2HN2gMuHkunDrm0k11FBa/wzFXUbSWPwC6sZQL46Kg1dTg1eVwrvg9/DOWw2Og5ZKJOozxOrqDpFEcAq9R5zHzxd48xf3ceD1N4/YN5nRQg63Za7p+avivvvQ6q+FZ39FmKxBxTkkHDQRmNSpnoO4GshhxNx+CWxmbPX8LWOMusapBTXT8VrmIKFPX1sfxUjaN4i0fmLpFG7cI1WfpJT3iTc04GTq6N7ZA14MJ+4Z9cuoglVQ8tFiHkmZQVO8OJrPovlsuS95SNSNIV4MrypJYdd2M9vXjRs5jKgyF44D/kHhAZXodRjgpjME2V7EiyNpo52k3qEVsYIdLxtxuJJPvnPiKX6OlP7efn8FK3XjOKk0+EXEz5vPEsx6RJDrAxF6bv3eCZ//8XPOP+E2RoqqEoxwqVRGWsB98H/ukQIflglHMpWiEGRwu/dSnHM28uTPOff/fIH9v7yf7Y+t582H11M7r4a5l6+ga/Mucu29hMWAN+99DHH+k5rZGQhDcm1RzdfQ581vf43ktFqc2gZIZowT8IsmRFKB+DtfBS+Oe+rbqTntIiQoktbwoJPSEAIfwgBJ15hrUXOji136YZzWLRReftK01duL5Ewxk557bqL+EzcAEG5ZC7EEpbxPqbc4kFI7GXHy3aQ/+FkIikjJOPvWu26lZslC4l17CasbTUhtEL27WmnZv5GaD3ySTV/6E4YeKh8ZFz392xO7gOOkkkM6I+FYk7x6RKR7iKVHRLrHy0jL2JCorkVjCdzeVmTFVfitu9j24GtUt5h5AGAyWDQIiaXjpBqrjBS0a2Znhskalv7v7xOUQjZ/7SvMu/5rpuEwMIOmAK5RzKzU/Hcp5nB69uN270W9JBqrMjctDQ8OUorgtcxBC3nCdb8meP5XZn5AqUjol6g6+0JyrR3k27sodmcJo5oGvTd/HTBPCo4reCmPQnehXJc6poRb1ppZ3mEIoU/rP34TJ9dlqsUFAbnf3msG0R0PN9tO9pa/pvW7n8dxHdrvvpkDN/0d4o696Mdv33HRqLSjTPywz7G0fcZETFtEvoMJGxUxIaSPqmpntO/LwMeAAPisqlaOjuskpF9fv9S6ncSSM1n13Tn4uzaz7of3IK7gpjO4ybip0RuGlLIFEjUp8h1Z1n/uUwTFYOCfVt04YcknzPZA7QykkEVKRjGzEgmrm3B3voxmGgnSDeDFkVKe0I2ZeQkiEAYELzxIbPYinIxJd3Xrm5GdrxF0taNByIEH76X+ggsJWneRWn1omYvi+mdpe/p5Wl9tw0uNtGT2xEOLOfb85EdM//Sfo/EU0z/0Ufb8y98DkDhlJZv+z/fRXz/PSV//FrnH/43cflPuMij6iOOQmTeTzk1jOxj+9BWXDFlT+K0SVnQuz7Epl77Wg8BpkST0BqKSkCKyDFPG7FTgSuDvRaQyYwaTjDBVS5hpRmYuJrFsFWd8/hqW/tcrkaoaMnPNrNx+hc++9iwL3nsefs43OfBJo7u+/W++TCmbI79tswmdaAiqePNPJXHxh8t5eUMSr2uGdD3auQ8JSkgha8I6fgEJfZxCD06+G+/0Cynt3Iz7tvcYhx8G4LiU9mwj395FzVkrKWx+7QjHH48UTFNNdXgpjzDQgfKFk4m+276Jv3ur+Z28+TKlB37C9ht/SOM7LqLp8ssIOoyMeKG7wPov/Cmp89/LrideY8t/vEjzRz6NhiGdG7eT78wNyDaPBf1zVkaLid7zL4vzV9VfDSoO8xSmIDGYgeVbVbWgqlsxZctWlcPGqYa6cTP46RcJM824TbNw6pqReBKvpoZ4TZpUUz1u3KXUWyS7q9XIQiQ93KRR+Cxl89HTQQ7xi+DGkNBHE5XZ8wfQ+KCqVSKmx++4iJ83YZ8oiye21EzQTCw9K8r3N8NgpWyesLeTzg1HZjNpYG4SqRktxNIxNAi5eO2TY39R44j/klHb1GhAXP0iQb5Ibn8H4rpG1K+QN3V8k+bJR0KfYtaUcpRofCUzt2Vgdnnunu+Pia0XPf3bURsXmEqTvMaSP+LgbOFZwOD/op3RtiMQkTUislZE1ra2to6xiZOfZCoFqoTJGsJYitYH7ido34Nz8rnEZp9E/ZI5Uew/SRgq2x95hQVXnU2uwzh8ADfmMe/6r1HKmsG+gVq2FVbAZTBO9gD4JZx8N27PfqSQNU4/DJFSAYIAHG8gD11LJfqe/0/83Vvxs3ny7V3sefg3zPnrHx/Rtrt0FeLF2Hb3Y7hxl6A4vFZNJdapHY5t11/Hrq/994HZym7TLAgDUktOJcz2kDz5DBZe98GBG6f6JeoWzaBx+UJqF0yn9Y5/JVGTwI27hLs30/ihNSTe90nmX30J83/3cgD2ffszI7Zny+f/YNgi8GOFKpQCHdFSqYyZ8xeRh0TklSGWqwcd8xXMfIGb+zcN0dSQn56q3qiqK1V1ZVNT0+hfwBREkxmcXBdOoYdUQy2xC34XdrwK88+g2J0l2VBDviOLG3MpdBfp3raXpddeSMuKk81cgPYsm/7iegqdvTh9nQQb1hohuWRtuS9tWJwlFyANs9ADewjSDYSpWpxSzsyDiKdQL4bT1wFbX0R2vobUNRObPgcnlabmsqOrk6rjoX6JYrZIX1tuyJh/vxTy4QJnEwFv+lyktpnYvJNxahtMSm8xb0Tx/CKFza/R/sQThD0dNJ2/Eg1DMnNbcGIe05bMJFGT4I3v/z+2/M3X0LX3kVi2Cm/Ve+hev/G47OjZ1TFGV3g0Jn6q55g5f1W9XFVPG2K5G0BErsNUCPuwHlSX2wnMGdTMbEwRecsYk8/ljLOKp0AcMudfSu7+m1DfaPGber5F3LiLl/JI1MTp3dVuwj+lEvP/4PdwYw5hMUCDgM57b8ZrnoWbbSfeOPvYBpSRMJnBSWeiazf9DylmUS9BmKzFr5+NM30BYV83wc4NOJk6wlyW/NqHaVm1jJZVy4ZuWEPEi5HrmLh1fQ8n2Po8M1e/j2krlqOlElKMCtn0p/OGAVo015t47xqaP/Jp3IYZxBaeRuMHPjowd6Tl3e+iasY0kvVpo97pxdFENU7bNpyYN6AZNRLOuP1+qhrS7PjqkTOrx4rJEPYpS/qBiFyJqQL2DlXtG7TrHuBnIvJ3wExgMfBMGUyckkjoE1Q34Wbb0Xwfvbtaqb5wFoFfJFGXoWf7Pvy8GarJzK6lZ2cXHRtMOci2dZtZ9MF3sfabt6FvtBNLp6jNZSGXhXlnlPnKjk5s+iLC3lYT5nFcM8ErKCJBCSd7AI0lwHHwmmYR9nSy4Xv/QOPyhVTPmY7bMGNA9+hwnFIfhX07jDxxMeSCJ/7ziGMG16itdPJ9WZza6bhi+ozO3GWEjoc6Ho7jEeayZk5HGCCpNG5Pq5nVm6lDc1kkVUP6kvfjb3oBr2kWM1e/DyeZRuqnI34ev24WXilH5pRTcDJ1ZG/5a9If+uqIbFv8o18cV6johFE4huJ0xVOu3LMfAAngQTE9radU9Y9V9VURuR14DRMO+pSqBmWyccoh+R5IpJFCls0//imphloKrz1D4uQVgKnpO+3kOex5ehNhMaCU93nz8e3Mv3g+Xds7KbW30Xx6C81nmSwXb8XQ1ZgqEb/pJCM9HQ1A4rhIFP7pl2vQRDXksiz8g9/h9X/4BYuuriWeKeJ3dw3TaIl8exde0jwpTQY0VoU/bS5SN5MQcLLtCBDUTic2u4g/bS5Orgsp5pCgZG6kDbOMYF40/8FbupJ8VOUtsfQscF2Cqum4nbtAHNyG6cjsU5C9x6cY3/LF/zvKVzs8/T3/iUy5sn1OUtU5qnpmtPzxoH3fUNVFqrpUVa1s9HgiYgY7+zrJzG2h+R3nUWzdT2n7Bro270Jch/SMBqYtbaFnTy+lbIlUfZK211rxcz4b73ySRH2Gjg3bB2a4ThTU8ZBoYpeKYxYvYW4GYTCQsxf2miIvvfuMln+pbR/F7qGlK7S7jUKHuXFMipm94piwmDgQ+kbCwUsODJCHyQxSKqBufOAmKqFvJsy50fVH2+MnGSVPSabN+6OniaBtN2G2h8KjtxHm++i88c/H/zpHgAKlUEe0VCqVkO1jqRTEwe3eM6BHH/Z0UujoYfNP78RNxgmLJd586AV2/nbbgAgcYDJ+AiVVn2T+t25i4Xd/Vs6reEuoG0PdaBay44IXjxK1Q4KdGwayl8IuI0t9xprL6NmxP8qASg3ZZscj9xOvqSIoGZXUyUC/E1cvGb2Omc/LcYwmkhtDvcQhN1G8GIQ+EhTNul9CahpJnnkRYbrB3Cz8gnlSCAOK2zfiNc0iteIS+vZ3DFvPt6woBKGOaKlUJu+UQ8tx4/R1QLaD4qZ1VM9qwps+l2ogKPnsf2ELhe4Cvbt7iaVjlHpLBMWQ7L4ssbRxCB1bOst7ASeAaBj1bM1r9RLmb2j+SlCktHMzYU8nudYDVJ+8jODFDdR89OvDttn0+e9Gf8fW9vEimUqR78tGg+IO+PnIaRfRRMZo+jgeUupDYwnT64/y//ufCFQciKfMWIpfQqubzM2j0AN+wYjlpaoI2vfS/fI6pq1aybRVFVcEC2V8BnNFZBpwGzAfU0TrGlU9Ir1JRH6CSaDZr6qnjaTtydEdsYwKpaaTTIH1OUuMpk+pSPKsi6ie1USyzpR6rJ5paqwWuguIK6Rb0jiu0LisicZlEzzlVgeN4IXBQCjCm7mQ0vYNOKk0nRu3s/+5N8CLM/Oy8VORrATyudzAZ3JQuylmyl5qaFJ6RUzpzv5wUD/ZDsTPm5ushqZeQlVdFBLyIEqL9aM4v4YBqeZ6Nt10J5tuqsw02EBHtpwgXwIeVtXFwMMcWkt9MP+CUUUYMbbnbzmIOPjT5lF67l/wGqYDUNy0jtYXNlDsyePGXWrn1hEUA/ycT1AKyHXk6esukNpUeVW6jougNNDbhyh10ff7X+BNn0vnI/dR7Okjuy+Lk0pT2PoGQwd8JifJVGogJbj/SUn7bwb9E7ow+6SUNfuj2d1OVa15X6lAWFVPEE+Dqnk6KBWQUp6gfS9Opt6kihYcKOaZfs7J5brcozKOA75XAxdH6zcBj2IyJQ+1R/VxEZl/PA3bnr9lALdrN6IhW+4xU+D9fdvRYp54TZrM3CbcmDOgR+/njSBXX3eBZNKjff0B2tcfKKf5J0SiZohKUpFjU8cj+/QjpJrMMXULG8i+9AyFzp7xNLHyEAdiSTReZcYA+mP3fgHxi1G6bBEpZAmq6sF1CRNp0+sXx6TQaoiU+iDXDY6DeDG0kEOLeZx0DZllp7HjkZfLfaVHcnwx/8Z+NYJoWXMcZ2pR1T0A0d/m0boE2/O3DKDJDFLK4+dMj7d72x7qTq8n395Fdn+PmY5fDCh0FylmS4TFAFeEjmyJGbMzXPLC02W+gtHlYA9XiDU0mm1BOFC5bLiB3klPJHkNUWGbKAtIQ3A0NOGdUs6EcqJMHtHQpNL6RTSWOljwpv/9YYB4cXBcnGINYRgS9HQiXqwixfD6s31GSJtfPL1GAAAW60lEQVSqDjtwISIPAdOH2PWVt2DaiLHO32JiudE/s5Pv4vT/+Sd0PfkojZe/i9KWVwmKPvF0jO6OPIXuAm7MIVGTICwG1NQkKG48QPvuI0sXTjjEMRkp4hys5AWEyRrc+mZKOzfR/tpuMrNr2bd2/YSUZBgVot9KmMgcsk3CwGT/hD4U80g8OfAZ9qd8hlX1OLkuE+cXB6eQRfO9aCGPhgFOPAnpGop7drL/ufW0v9E6IPhWSYxm2EdVLx9un4jsE5EZqrpHRGYA+0flpNiwjyVCghJOroswVoVT10S8Jo04Lq/d9BC5dpPH7ud8gmKAhkopWzTVqfI+01rS/P6+V8t8BSdOMnWwJy/9eexu3Kh9JqsQxyVWHadl5clT1/HDobF+QEUQv4BT6DXS2CWTtaN+KUqbTRC07kQLfUi+h54HbsNf97hRkO3YZ5pMpXHiSSSVNnV/k3GCok/b+nZ2P7e3XFc6PKqE4ciWE+Qe4Lpo/Trg7hNtsB/r/C0kUykz2Bk5PIDklR/BqW8hDJVUgxm0S9TE0UDJdeQJA0Wjx3FnFAtkVARRHQIp9OL27MMp9OBOm07hQBfzLj+L1hc2lNvCspFMpVARNNLyUREkDExd3jCajC+CpDKmAFA+i79nC5Kswt+/C+1uIzl3vknzbNtJ0NUOVXWE9bPR2cvME0AxT9/edrY9ug3HFRx3KL3H8qKMW7bPDcAVIrIRuCJ6jYjMFJH/6D9IRG4BngSWishOETmm0JEN+1gAk8cepmqRYg5/7zYkUricdlIj4jr0tGcJQ6XQXSBRk8BLeXTv6KHuzBrOvGNyFVsbmHTUP0u1mEMLfcTrqkmcvJK6NX9TZgvLi6iaG4AI0h/6EAcJfcJYCi/XZUI6Sy4AwMUUsad9L0H7Xtqffg6AaWGAVNUckmLrpDNIVQ3pubPo6y5QXZ/kitefG+9LHBHjke2jqu3AZUNs3w1cNej1h463bev8LRSyPSY263i43XvJ797Glr+/iennnIwGIeI6xKOJXDNXzaTrzS78nM+0xfWTzvETBjhBFoIi6no4+SKlHRtp/fWjlLJ5FoxQaGyyUtq/Db+6BQ8dKGbfL98QxlJGwXUIFVd3zukQBEjo0wjmKcGLI4kkwe5NaKmEW9tAGD0tiOty8rsXcco/3zPu1zgSjJ7/xFZ2s87fYgh9JPCNeFkYMOO8ZeRaO8l39uElPUqRmqcbcwiKAX1tOVKNky/bJVEzjWLnfjMQKQ6S60YSSaYtX1pu0yqOgQweDY0aarzqGG+IwjeOi5OpQ7w4EosRZntM3ef+uQL5LAdeWn/CJRd/fdY5AGOShdYf9pnIWOdvgaBEmKzF69oN2Q682Ytwu3rI7mknnknS/kYbdfNqSdSnCYs+iZoE2X19nPfIY+W2fExwckalMzZ9EUxfRByYBLJso0KseT7auX9g0Dc51PyIYXDnnUG4ZS2xmfPNrOAoAyjYuRnNZwl7Owl7OgmyPRxYv4NitjQqNt83Zznv2bFuVNoazERX9bTOf4pT2rcVSWYgiP7R0vW4QKK5iarmfaSa69jyq83MPn8R4joUe7Kc/9jjZbV5rPFmLC63CRVNvO6tzzMKkxmcPiPypo5H2LoTt77ZTCjMZSnu2cmeJ1+hbf2Jzxi/5IWneWDJ2UNWUDtRlMqu0jUSbLaPBUqm8pJ6CdSLU9zwAlrIk2quI7P8bKaf2YKbNH1fcZwJVW/WUll4M5finHQuzqJVUMgStO8x+f1VGdQvkWvtYNdTO/GS3qiEa2JJj/qFday96ogx0xPDqnqeGCLyZ8B3gCZVbYu2fRn4GBAAn1XVSTaiWFlI6KPJDIgQJqqRWJLYnCWEfT1U1zZQ2rWZpjPmsXftFgBS9UmW/fTfy2y1ZTIQtO/BqW8m7Gon7Omk/aU32Pf8NqYtrh+YZX6iXLruGX7z9gsBeGT5qlFpE6KYfwU79pFQNucvInMweavbB21bBlwLnIop4/iQiCyx1bzGjjCZQXJdkKodyPWXaTNw0zVIPMmBtS8gjoMbd5hzyXIcd5Ll9FvKRvyCa8jd9V1KHR14VSnaXtlB4+mzmf+tm0b1PI3LmvBScTbe+8aotakKRX9iZ/uUM+zzXeCLmJtoP1cDt6pqQVW3ApuA0btdW45EQzPw1o84BLXTCVoWo4tXUb/8FFrXbQNgx6/XTbgKXZbKJrX682y++ylwXNo3dlDoGF2xvP4Q5ab7NxBPx7hy8wuj0q4yspBPJT8dlKuA+/uAXar6ksghs/dmAU8Ner0z2jZUG2uANQBz584dI0snP06uCykVcGgnTDeYm0ExF4WCTD3VZX/6MeIXXFNuUy2TlLPu+hUAl3/giwA8c+WlrPrlI6PS9tIb7+S25mUsvXDO6M5JURv2GZZjKNX9OfDOod42xLYhP2FVvRG4EWDlypUT+1soIxL6aNd+pLYZ0hDG0zh+0ezUEKe6rrwGWqYED56yAjfmMm1x/ai1+fKHrmLBe89l6YVzmLZk5qi1Czbmf1SGU6oTkdOBBUB/r3828LyIrML09OcMOnw2sHusbLQY2eLSm+uJL0kdnLDjxY2qZejjnTHUPdpiGV28pMeenT1MW1xPLHnibunxc86n+XSTktr6Whutr7Ux95sn3OwAanv+x4+qvsygggQisg1YqaptInIP8DMR+TvMgO9i4JnxtnEqUXziLhKnn4fWNKNewsza7Next1jGiUteeJqfNS0jUZMgKJ54fkdmRjXVs5oIi/6YaQNZ5z+KqOqrInI78BrgA5+ymT5jR+npu0iuvAzN9Qz09MUvgF80zj+0H71l/Oj1Q/z8iaV4fqd6CSlXWLFiOk7Mo3dXKzWjZN9gQlUKEzzbp+zOX1XnH/b6G8A3ymPN1CLs6cCZvoDS7pfxGuYiQcnIOjsO6sRPaCanxXK8rOlYf0Lv3/ftzwAwPx0nLIa0rP593NNGeXLXIGzP3zJx8eJkH7mD+IzZqOOZKkzi4BR6IRidSTYWy3jzts9dTN2qc8fU8U+GmL+Vd5jCFLa+QSmbo+OFl8yG/nqqQKzJps9aJhYv/tNvANh830sUto7ehK7hCFRHtFQq1vlPYZ7/4YPc8pV7qD1pHppIQ1BCCr2oOJT2byu3eRbLcbFlfx9nzsng531qPvr1MT2XneRlmdBUN1fxjsX1uPVNAwO+gNFVn9hjWZYpiCuwdW+WNeNQ+WsyyDtY5z+FceMu3Tt6wIuDG4sKupQgFNSzCvaVRLDjZVMNyzIsJzpgfDyYSV7W+VsmIPfPP4PTr1tFPFNF0LqLWLbdFC0PiqgbJz5tdGdEWk4cf/cbeDNtRbGKQCs7pDMSbMx/ivLubS+RbKihY8MOEqeswOnrxOnrQIq5Q4XeLBWBO+d06/griH55Bxvzt0xIqloa0SCE2mbCeDXim6IuQ1Wy6vrxV9n9hCmFV6lFtS2W8UIV/Ap27CPBOv8pjFPbQHUqDTDg+PtTPQfT/c9/MZ5mWSwVjxV2s0xoSvt3kzrjfMh2IskMGk8NeVysoRGwPX6LpR9Vtdk+lolL5g//iuDNlyDXPaDjD0dmlqTe99lymWixVCzj0fMXkWnAbcB8YBtwjap2HHbMHOCnGAn9ELhRVb93rLbtgO8UR704Wj0N7dyPv+kFiMo0BtteHPVz5XM58rncqLdrsYw3On4F3L8EPKyqi4GHo9eH4wP/Q1VPAc4FPhWVxD0qtuc/xfFmnQJAoCHetBlmslepMCbnSqZS1vlbJg06PjH/q4GLo/WbgEeB6w+xQ3UPsCda7xGR1zEVEF87WsPW+VsAk0oYvPnSwM1grEimhh5XsFgmEqoQjtz5N4rI2kGvb4wqEY6Elsi5o6p7ROSoUrsiMh84C3j6WA1b528ZwJ13RrlNsFgmCIqOXLStTVVXDrfzGCVvR4yIVAN3AJ9T1e5jHV825y8inwE+jYlX3aeqX4y2fxn4GBAAn1XVUay6bLFYLKOAQjBK2T7DlbwFEJF9IjIj6vXPAPYPc1wM4/hvVtU7R3Lesjh/EbkEE8tarqqF/keZaJDiWuBUTBnHh0Rkia3mZbFYKgkFxqnS6T3AdcAN0d+7Dz9ATDH0HwOvq+rfjbThcmX7fAK4QVULAKrafze7GrhVVQuquhXYBKwqk40Wi8UyLKo6ouUEuQG4QkQ2AldErxGRmSLyH9ExFwD/FbhURF6MlquO1XC5wj5LgAtF5BtAHvgzVX0WM0L91KDjdkbbLBaLpXI4vgHft34a1XbgiJJkqrobuCpafwKQ4217zJz/MQYxPKAek5P6NuB2EVnI0Bcw5CcsImuANQBz59qqUxaLZTzR8Ur1HDPGzPkfYxDjE8Cdap6JnhGREGjE9PTnDDp0NrB7mPZvBG4EWLly5cT+FiwWy4RCFYJgYss7lCvmfxdwKYCILAHiQBtmcONaEUmIyAJgMfBMmWy0WCyWYdFQR7RUKuWK+f8E+ImIvAIUgeuip4BXReR2zMw0H/iUzfSxWCyVSCU79pFQFuevqkXgvwyz7xvAN8bXIovFYhk5qjouA75jiZ3ha7FYLG+BUUjjLCvW+VssFstbYJwmeY0Z1vlbLBbLcaKjKO9QLqzzt1gsluNF7YCvxWKxTEGU0Mb8LRaLZWphhN2s87dYLJaphQ37WCwWy9TE5vlbLBbLFENVCSe4to91/haLxfIWsD1/i8VimYJoOLFlx6zzt1gsluNF1Tp/i8VimWoo1vlbLBbL1EOVsFQstxUnhHX+FovFcrzYsI/FYrFMTSa68y9LGUcROVNEnhKRF0VkrYisGrTvyyKySUTeEJF3lcM+i8ViORr9Mf+RLJVKuXr+3wa+pqr3i8hV0euLRWQZcC1wKjATeEhElthSjhaLpaLQid/zL5fzV6AmWq8FdkfrVwO3qmoB2Coim4BVwJPjb6LFYrEMhxJa5/+W+BzwgIj8LSb0dH60fRbw1KDjdkbbjkBE1gBrAObOnTt2llosFsthqCqhb7N9hkREHgKmD7HrK8BlwOdV9Q4RuQb4MXA5IEMcP+QcalW9EbgRYOXKlRN7nrXFYplYqKLB2Pf8RWQacBswH9gGXKOqHYcdkwQeBxIYn/4LVf3LY7U9Zs5fVS8fbp+I/BT4k+jlz4F/itZ3AnMGHTqbgyEhi8ViqRjGKeb/JeBhVb1BRL4Uvb7+sGMKwKWq2isiMeAJEblfVZ86vLHBlCXbB+PQ3xGtXwpsjNbvAa4VkYSILAAWA8+UwT6LxWIZHh23bJ+rgZui9ZuA1UeaoqqqvdHLWLQcMxpSrpj/x4HviYgH5Ili96r6qojcDrwG+MCnbKaPxWKpPI5rklejiKwd9PrGKGw9ElpUdQ+Aqu4RkeahDhIRF3gOOAn4oao+fayGy+L8VfUJYMUw+74BfGN8LbJYLJaRY8o4jljPv01VVw638xjjoyOzx3SSzxSROuDfROQ0VX3laO+xM3wtFovleBnFbJ9jjI/uE5EZUa9/BrD/GG11isijwJXAUZ1/uWL+FovFMnFRk+c/kuUEuQe4Llq/Drj78ANEpCnq8SMiKUzm5PpjNWx7/haLxXKcKIxLqidwA3C7iHwM2A58AEBEZgL/pKpXATOAm6K4vwPcrqr3Hqth6/wtFovleBknVU9VbcfMizp8+27gqmh9HXDW8bZtnb/FYrEcN1bS2WKxWKYek0DeQVQnvjKCiLQCb45R841A2xi1/VaxNo0Ma9PIqUS7xsqmearadCINiMgvMfaNhDZVvfJEzjcWTArnP5aIyNqj5eiWA2vTyLA2jZxKtKsSbZpM2FRPi8VimYJY52+xWCxTEOv8j81INTjGE2vTyLA2jZxKtKsSbZo02Ji/xWKxTEFsz99isVimINb5WywWyxTEOv8IEfmAiLwqIqGIrBy0fb6I5ETkxWj5h0H7VojIyyKySUS+LyJDlaEcdZuifV+OzvuGiLxrvGwawsa/EpFdgz6fq45l43ggIldG590UVUAqCyKyLfo+XuzXdBeRaSLyoIhsjP7Wj7ENPxGR/SLyyqBtw9owHt/bMDZV5G9p0qKqdjHjHqcAS4FHgZWDts8HXhnmPc8A52FqD98PvHucbFoGvISp2bkA2Ay442HTEDb+FfBnQ2wf1sZx+C7d6HwLgXhkx7Iy/a62AY2Hbfs28KVo/UvAt8bYhouAswf/joezYby+t2Fsqrjf0mRebM8/QlVfV9U3Rnp8pK1do6pPqvmF/pQhSqyNkU1XA7eqakFVtwKbgFXjYdNxMKSN43TuVcAmVd2iqkXg1sieSuGYpflGE1V9HDgwQhvG5XsbxqbhKOdvadJinf/IWCAiL4jIYyJyYbRtFqbgfD87o23jwSxgxxDnLpdNnxaRddGjfH/4YDgbx4NynvtwFPiViDwnImuibYeU5gOGLM03xgxnQ7k/u0r7LU1appSw29HKpanqEUUSIvYAc1W1XURWAHeJyKmYsMrhHHfe7Fu0abhzj4pNR5zs6GXmfgT8r+g8/wv438AfjZUtI6Sc5z6cC1R1d1R79UEROWaRjTJTzs+uEn9Lk5Yp5fz1KOXSjvKeAlCI1p8Tkc3AEkzvY/agQ2cDu8fDpujcc4Y496jYdDgjtVFE/h/QX0RiOBvHg3Ke+xDU6K6jqvtF5N8w4YrjKs03RgxnQ9k+O1Xd179eQb+lSYsN+xyDqESaG60vBBYDW6JH5R4ROTfKqPlDhiixNkbcA1wrIgkRWRDZ9Ew5bIocRz+/y8G6oUPaOJa2DOJZYLGILBCROHBtZM+4IiJpEcn0rwPvxHw+xyzNNw4MZ0PZvrcK/S1NXso94lwpC+bHthPTy98HPBBt/z3gVUy2wfPA7wx6z0rMD3Qz8AOiGdNjbVO07yvRed9gUEbPWNs0hI3/CrwMrMP8k844lo3j9H1eBWyIzv+VMv2mFka/m5ei39BXou0NwMPAxujvtDG24xZM+LIU/Z4+djQbxuN7G8amivwtTdbFyjtYLBbLFMSGfSwWi2UKYp2/xWKxTEGs87dYLJYpiHX+FovFMgWxzt9isVimINb5W8qOiPSOQZvv61fzFJHVIrLsLbTx6OFqqhbLZME6f8ukRFXvUdUboperMcqQFoslwjp/S8Ughu+IyCuRBv4Ho+0XR73wX4jIehG5ub9OgYhcFW17Qkz9gnuj7R8RkR+IyPnA+4DvRBrxiwb36EWkUUS2RespEbk1Eha7DUgNsu2dIvKkiDwvIj8Xkerx/XQsltFlSmn7WCqe9wNnAmcAjcCzIvJ4tO8s4FSMpstvgAvEFEf5R+AiVd0qIrcc3qCq/lZE7gHuVdVfAMjw9W0+AfSp6nIRWY6Z0Y2INAJfBS5X1ayIXA/8KfD10bhoi6UcWOdvqSTeDtyiqgFGeOwx4G1AN0a7aCeAiLyIKbLTi9FZ2hq9/xZgzRGtjpyLgO8DqOo6EVkXbT8XEzb6TXTjiANPnsB5LJayY52/pZI4WsnJwqD1APPbfaslKn0OhjyTh+0bSu9EgAdV9UNv8XwWS8VhY/6WSuJx4IMi4opIE6YnfjT1xvXAQhGZH73+4DDH9QCZQa+3ASui9d8/7PwfBhCR04Dl0fanMGGmk6J9VSKyZATXY7FULNb5WyqJf8MoOr4EPAJ8UVX3DnewquaATwK/FJEnMMqnXUMceivwhaga2yLgb4FPiMhvMWML/fwIqI7CPV8kuvGoaivwEeCWaN9TwMkncqEWS7mxqp6WCY2IVKtqb5T980Ngo6p+t9x2WSyVju35WyY6H48GgF8FajHZPxaL5RjYnr/FYrFMQWzP32KxWKYg1vlbLBbLFMQ6f4vFYpmCWOdvsVgsUxDr/C0Wi2UK8v8BgwMhK9zZFScAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "out.sel(Time=\"2015-01\").plot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Another solution: groupby\n", + "In the loop above, we see we handle all the data for each spatial point at once. That's typically what groupby would do if applied to Lon and Lat. Once we group the data per spatial point, we need to apply a function. This function sould return the stacked data for that point (i.e. the timeseries) with the timestamp as an index. This way to overall DataFrame resulting from the groupby will have Lon, Lat and Time as indexes and 1 column of data. This column will be the timeseries at each point stack one after the other. Let's see what it would look like:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "# We reuse the time array we used for defining the out DataArray\n", + "# time = pd.date_range(start=f'01/{years[0]}',\n", + "# end =f'01/{years[-1]+1}', freq='M')\n", + "\n", + "def time_to_date(df):\n", + " df_stack = df[months].stack()\n", + " return(pd.DataFrame(df_stack.values, index=time))" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 2min 8s, sys: 2.63 s, total: 2min 10s\n", + "Wall time: 2min 10s\n" + ] + } + ], + "source": [ + "%%time\n", + "df2 = df.groupby([\"Lon\",\"Lat\"]).apply(time_to_date)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
0
LonLatTime
-179.7566.251901-01-310.0
1901-02-280.0
1901-03-310.0
1901-04-300.0
1901-05-310.0
............
179.7571.252015-08-310.0
2015-09-300.0
2015-10-310.0
2015-11-300.0
2015-12-310.0
\n", + "

81682200 rows × 1 columns

\n", + "
" + ], + "text/plain": [ + " 0\n", + "Lon Lat Time \n", + "-179.75 66.25 1901-01-31 0.0\n", + " 1901-02-28 0.0\n", + " 1901-03-31 0.0\n", + " 1901-04-30 0.0\n", + " 1901-05-31 0.0\n", + "... ...\n", + " 179.75 71.25 2015-08-31 0.0\n", + " 2015-09-30 0.0\n", + " 2015-10-31 0.0\n", + " 2015-11-30 0.0\n", + " 2015-12-31 0.0\n", + "\n", + "[81682200 rows x 1 columns]" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "df2.index.names = [\"Lon\", \"Lat\", \"Time\"]\n", + "df2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see the groupby calculation takes about the same time but we only get a DataFrame as output and not a DataArray. So we still need to convert this DataFrame to a DataArray. Unfortunately, the conversion below crashes the Jupyter Kernel. This is because the DataFrame contains only the landpoints. The DataArray created would be sparse with irregular coordinates. There are ways to solve this issue but we'll see those in another blog on sparse matrices as this one is long enough. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#out = df2.to_xarray()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Conclusion\n", + "Although the code with `groupby` is a lot more compact, it might not be as intuitive to write at first. In this case, it shows a well-thought loop can actually be faster than a `groupby` as the conversion to a spatially full DataArray can be done at the same time as the concatenation of the months and years.\n", + "\n", + "At the opposite, after experience using Pandas and Xarray, the `groupby` solution might be the first one will think about. It isn't a bad solution but then one needs to know how to fill the missing indexes." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Solution summary\n", + "Below is the solution with the loop containing only the necessary code in one cell to clarify what the final code would look like." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "import xarray as xr\n", + "\n", + "fname = \"/g/data/w35/lt0205/research/lpj_guess/runs/CRUNCEP/mgpp.out\"\n", + "df = pd.read_csv(fname, header=0, delim_whitespace=True)\n", + "\n", + "months=list(df.columns)\n", + "months=months[3:]\n", + "years = np.unique(df.Year)\n", + "nyears = len(years)\n", + "nmonths = 12\n", + "years.sort()\n", + "\n", + "# Create the axes\n", + "time = pd.date_range(start=f'01/{years[0]}',\n", + " end =f'01/{years[-1]+1}', freq='M')\n", + "# We'll use a generic way to create a regular grid from [-180,180] and\n", + "# [-90, 90] when knowing the resolution. Feel free to reuse as needed.\n", + "dx = 0.5\n", + "Lon = xr.DataArray(np.arange(-180.+dx/2., 180., dx), dims=(\"Lon\"),\n", + " attrs={\"long_name\":\"longitude\", \"unit\":\"degrees_east\"})\n", + "nlon = Lon.size\n", + "dy = 0.5\n", + "Lat = xr.DataArray(np.arange(-90.+dy/2., 90., dy), dims=(\"Lat\"),\n", + " attrs={\"long_name\":\"latitude\", \"unit\":\"degrees_north\"})\n", + "nlat = Lat.size\n", + "\n", + "# Output array\n", + "out = xr.DataArray(np.zeros((nlat, nlon, nyears*nmonths)),\n", + " dims=(\"Lat\",\"Lon\",\"Time\"),\n", + " coords=({\"Lat\":Lat, \"Lon\":Lon, \"Time\":time}))\n", + "out[:] = np.nan\n", + "\n", + "# We stack the months columns of the whole DataFrame at once since we\n", + "# don't have missing years for any point. Doing it once will save time.\n", + "df_stack = df[months].stack()\n", + "\n", + "for nr in range(0,len(df.index),nyears):\n", + " rows = df[nr:nr+nyears]\n", + " thislon = rows[\"Lon\"].min()\n", + " thislat = rows[\"Lat\"].min()\n", + " out.loc[dict( \n", + " Lon=thislon,\n", + " Lat=thislat)] = df_stack[nr*nmonths:(nr+nyears)*nmonths]\n", + " \n", + "out.sel(Time=\"2015-01\").plot()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}