{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# problem description\n", "\n", "_The approach might be useful for those working on data collected over sub-linear profiles, e.g. for planning or implementation purposes..._\n", "\n", "for Site D: \n", "\n", "profiles parallel to the main one at distace (perpendicular to the mean direction): \n", "```\n", "meters:\n", "-40, -30, -20, -10, +10, +20, +30, +40\n", " ```\n", " \n", "for Site E\n", " \n", "profiles parallel to the main one at distace (perpendicular to the mean direction): \n", "\n", "```\n", "meters:\n", "-20, -15, -10, -5, +5, +10, +15, +20\n", "```\n", " \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# change of CRS (done with ogr2ogr)\n", "\n", "original projection:\n", "\n", "```\n", "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs\n", "```\n", "\n", "desired projection:\n", "\n", "```\n", "+proj=utm +zone=28 +north +ellps=GRS80 +datum=GRS80 +units=m +no_defs\n", "```" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "+ellps=GRS80 +no_defs +proj=utm +units=m +zone=28\n" ] } ], "source": [ "# files as WGS84 geographic, as exported from RTK processing\n", "\n", "input_raw_D = './0_location-profiles-RAW/AGPA-G_D2.shp'\n", "input_raw_E = './0_location-profiles-RAW/AGPA-E.shp'\n", "output_UTM_D = './1_location-profiles-UTM/ERT-siteD-UTM28-GRS80.shp'\n", "# get_input projection\n", "import fiona\n", "from fiona.crs import to_string\n", "\n", "c = fiona.open(output_UTM_D, 'r')\n", "# print(c.crs)\n", "print(to_string(c.crs))\n", "\n", "# wgs = from_epsg\n", "# wgs_proj4_string = to_string(wgs)\n", "\n", "# print(wgs_proj4_string)\n", "\n", "UTM_D = './0_location-profiles-RAW/AGPA-G_D2.shp'\n", "UTM_E = './0_location-profiles-RAW/AGPA-E.shp'\n", "\n", "# one could reproject via \n", "# https://pcjericks.github.io/py-gdalogr-cookbook/projection.html#reproject-a-geometry" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Input and output" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-40, -30, -20, -10, 10, 20, 30, 40]\n", "[-20, -15, -10, -5, 5, 10, 15, 20]\n", "-40\n" ] } ], "source": [ "# arrays for shifted points\n", "delta_D = [-40, -30, -20, -10, +10, +20, +30, +40]\n", "delta_E = [-20, -15, -10, -5, +5, +10, +15, +20]\n", "\n", "delta_TEST = [10]\n", "\n", "print(delta_D)\n", "print(delta_E)\n", "\n", "print(delta_D[0])" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "D\n", "./1_location-profiles-UTM/ERT-siteD-UTM28-GRS80.shp\n", "./2_output/agpa_D-shifted.csv\n", "./2_output/agpa_D-shifted.geojson\n", "[-40, -30, -20, -10, 10, 20, 30, 40]\n" ] } ], "source": [ "delta_select = 'D'\n", "delta = delta_D\n", "\n", "# input\n", "input_file = './1_location-profiles-UTM/ERT-site' + delta_select + '-UTM28-GRS80.shp'\n", "\n", "\n", "# output CSV/json\n", "\n", "outcsv = \"./2_output/agpa_\" + delta_select + \"-shifted.csv\"\n", "outjson = './2_output/agpa_' + delta_select + '-shifted.geojson'\n", "\n", "# delta\n", "\n", "# outjson = '2_output/test.geojson'\n", "# outjson = '2_output/agpa_D-shifted.geojson'\n", "# outjson = '2_output/agpa_E-shifted.geojson'\n", "\n", "print(delta_select)\n", "print(input_file)\n", "print(outcsv)\n", "print(outjson)\n", "print(delta)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Data load" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "schema is:\n", "{'properties': OrderedDict([('antenna he', 'float:24.15'), ('lateral rm', 'float:24.15'), ('name', 'str:80'), ('projected', 'str:80'), ('rms', 'str:80'), ('sample cou', 'float:24.15'), ('solution s', 'str:80')]), 'geometry': 'Point'} \n", "\n", "geometries are: \n", "{'type': 'Point', 'coordinates': (648601.4745557815, 3229147.9738091887)}\n", "{'type': 'Point', 'coordinates': (648600.5361209628, 3229143.138978075)}\n", "{'type': 'Point', 'coordinates': (648599.7661030253, 3229138.156342755)}\n", "{'type': 'Point', 'coordinates': (648598.7990097381, 3229133.272981949)}\n", "{'type': 'Point', 'coordinates': (648597.9483836357, 3229128.388791427)}\n", "{'type': 'Point', 'coordinates': (648596.9545805457, 3229123.46105742)}\n", "{'type': 'Point', 'coordinates': (648596.0947405221, 3229118.56331744)}\n", "{'type': 'Point', 'coordinates': (648595.1762005028, 3229113.609735632)}\n", "{'type': 'Point', 'coordinates': (648594.2721767295, 3229108.636121245)}\n", "{'type': 'Point', 'coordinates': (648593.2952063341, 3229103.763752384)}\n", "{'type': 'Point', 'coordinates': (648592.6392054026, 3229098.807557541)}\n", "{'type': 'Point', 'coordinates': (648591.786515621, 3229094.036386801)}\n", "{'type': 'Point', 'coordinates': (648590.9748514218, 3229088.9411467803)}\n", "{'type': 'Point', 'coordinates': (648589.9975829801, 3229084.1011719736)}\n", "{'type': 'Point', 'coordinates': (648589.144341087, 3229079.2450881554)}\n", "{'type': 'Point', 'coordinates': (648588.2456953059, 3229074.411880112)}\n", "{'type': 'Point', 'coordinates': (648587.5794628204, 3229069.5876915418)}\n", "{'type': 'Point', 'coordinates': (648586.5790606767, 3229064.5499946526)}\n", "{'type': 'Point', 'coordinates': (648585.8623549817, 3229059.698842795)}\n", "{'type': 'Point', 'coordinates': (648585.0456319387, 3229054.725776453)}\n", "{'type': 'Point', 'coordinates': (648584.3042273485, 3229050.012683716)}\n", "{'type': 'Point', 'coordinates': (648583.4001623433, 3229045.0961384033)}\n", "{'type': 'Point', 'coordinates': (648582.5647751306, 3229039.665844577)}\n", "{'type': 'Point', 'coordinates': (648581.665703492, 3229035.211213602)}\n", "{'type': 'Point', 'coordinates': (648580.9057197726, 3229030.0542368554)}\n", "{'type': 'Point', 'coordinates': (648579.9821064426, 3229025.2335872888)}\n", "{'type': 'Point', 'coordinates': (648579.2077628493, 3229020.400097237)}\n", "{'type': 'Point', 'coordinates': (648578.3236284078, 3229015.495011512)}\n", "{'type': 'Point', 'coordinates': (648577.4638452544, 3229010.754280825)}\n", "{'type': 'Point', 'coordinates': (648576.5390187646, 3229005.727932812)}\n", "{'type': 'Point', 'coordinates': (648575.5946259949, 3229000.8785012243)}\n", "{'type': 'Point', 'coordinates': (648574.9364234685, 3228995.664924923)}\n", "{'type': 'Point', 'coordinates': (648574.0533106204, 3228991.0935741314)}\n", "{'type': 'Point', 'coordinates': (648573.1830404663, 3228986.0321496315)}\n", "{'type': 'Point', 'coordinates': (648572.2467538107, 3228980.9361673873)}\n", "{'type': 'Point', 'coordinates': (648571.4503726437, 3228976.1057025185)}\n", "{'type': 'Point', 'coordinates': (648570.5414256933, 3228970.8126438186)}\n", "{'type': 'Point', 'coordinates': (648569.672187713, 3228966.050441338)}\n", "{'type': 'Point', 'coordinates': (648568.3397000737, 3228961.1769017302)}\n", "{'type': 'Point', 'coordinates': (648567.7542779779, 3228956.3556339275)}\n", "{'type': 'Point', 'coordinates': (648566.9048521187, 3228951.458135281)}\n", "{'type': 'Point', 'coordinates': (648566.109367921, 3228946.6162900357)}\n", "{'type': 'Point', 'coordinates': (648565.0290407431, 3228941.8700283184)}\n", "{'type': 'Point', 'coordinates': (648563.3972566999, 3228936.4794101315)}\n", "{'type': 'Point', 'coordinates': (648562.8736383251, 3228931.9387767552)}\n", "{'type': 'Point', 'coordinates': (648561.503112745, 3228927.599268888)}\n", "{'type': 'Point', 'coordinates': (648560.9811715708, 3228922.554788729)}\n", "{'type': 'Point', 'coordinates': (648559.6440178016, 3228917.816313375)}\n" ] } ], "source": [ "# load data\n", "import fiona\n", "import shapely \n", "\n", "c = fiona.open(input_file, 'r')\n", "c.crs\n", "\n", "# number of points, records\n", "len(c) # should be 48 for AGPA-G\n", "\n", "\n", "# PRINT SCHEMA\n", "print('schema is:')\n", "print(c.schema, '\\n')\n", "\n", "# PRINT GEOMETRIES\n", "print('geometries are: ')\n", "for feat in c:\n", " print(feat['geometry'])" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "648601.4745557815\n", "3229147.9738091887\n", "648600.5361209628\n", "3229143.138978075\n", "648599.7661030253\n", "3229138.156342755\n", "648598.7990097381\n", "3229133.272981949\n", "648597.9483836357\n", "3229128.388791427\n", "648596.9545805457\n", "3229123.46105742\n", "648596.0947405221\n", "3229118.56331744\n", "648595.1762005028\n", "3229113.609735632\n", "648594.2721767295\n", "3229108.636121245\n", "648593.2952063341\n", "3229103.763752384\n", "648592.6392054026\n", "3229098.807557541\n", "648591.786515621\n", "3229094.036386801\n", "648590.9748514218\n", "3229088.9411467803\n", "648589.9975829801\n", "3229084.1011719736\n", "648589.144341087\n", "3229079.2450881554\n", "648588.2456953059\n", "3229074.411880112\n", "648587.5794628204\n", "3229069.5876915418\n", "648586.5790606767\n", "3229064.5499946526\n", "648585.8623549817\n", "3229059.698842795\n", "648585.0456319387\n", "3229054.725776453\n", "648584.3042273485\n", "3229050.012683716\n", "648583.4001623433\n", "3229045.0961384033\n", "648582.5647751306\n", "3229039.665844577\n", "648581.665703492\n", "3229035.211213602\n", "648580.9057197726\n", "3229030.0542368554\n", "648579.9821064426\n", "3229025.2335872888\n", "648579.2077628493\n", "3229020.400097237\n", "648578.3236284078\n", "3229015.495011512\n", "648577.4638452544\n", "3229010.754280825\n", "648576.5390187646\n", "3229005.727932812\n", "648575.5946259949\n", "3229000.8785012243\n", "648574.9364234685\n", "3228995.664924923\n", "648574.0533106204\n", "3228991.0935741314\n", "648573.1830404663\n", "3228986.0321496315\n", "648572.2467538107\n", "3228980.9361673873\n", "648571.4503726437\n", "3228976.1057025185\n", "648570.5414256933\n", "3228970.8126438186\n", "648569.672187713\n", "3228966.050441338\n", "648568.3397000737\n", "3228961.1769017302\n", "648567.7542779779\n", "3228956.3556339275\n", "648566.9048521187\n", "3228951.458135281\n", "648566.109367921\n", "3228946.6162900357\n", "648565.0290407431\n", "3228941.8700283184\n", "648563.3972566999\n", "3228936.4794101315\n", "648562.8736383251\n", "3228931.9387767552\n", "648561.503112745\n", "3228927.599268888\n", "648560.9811715708\n", "3228922.554788729\n", "648559.6440178016\n", "3228917.816313375\n", "[648601.47455578 648600.53612096 648599.76610303 648598.79900974\n", " 648597.94838364 648596.95458055 648596.09474052 648595.1762005\n", " 648594.27217673 648593.29520633 648592.6392054 648591.78651562\n", " 648590.97485142 648589.99758298 648589.14434109 648588.24569531\n", " 648587.57946282 648586.57906068 648585.86235498 648585.04563194\n", " 648584.30422735 648583.40016234 648582.56477513 648581.66570349\n", " 648580.90571977 648579.98210644 648579.20776285 648578.32362841\n", " 648577.46384525 648576.53901876 648575.59462599 648574.93642347\n", " 648574.05331062 648573.18304047 648572.24675381 648571.45037264\n", " 648570.54142569 648569.67218771 648568.33970007 648567.75427798\n", " 648566.90485212 648566.10936792 648565.02904074 648563.3972567\n", " 648562.87363833 648561.50311275 648560.98117157 648559.6440178 ] [3229147.97380919 3229143.13897807 3229138.15634276 3229133.27298195\n", " 3229128.38879143 3229123.46105742 3229118.56331744 3229113.60973563\n", " 3229108.63612125 3229103.76375238 3229098.80755754 3229094.0363868\n", " 3229088.94114678 3229084.10117197 3229079.24508816 3229074.41188011\n", " 3229069.58769154 3229064.54999465 3229059.6988428 3229054.72577645\n", " 3229050.01268372 3229045.0961384 3229039.66584458 3229035.2112136\n", " 3229030.05423686 3229025.23358729 3229020.40009724 3229015.49501151\n", " 3229010.75428083 3229005.72793281 3229000.87850122 3228995.66492492\n", " 3228991.09357413 3228986.03214963 3228980.93616739 3228976.10570252\n", " 3228970.81264382 3228966.05044134 3228961.17690173 3228956.35563393\n", " 3228951.45813528 3228946.61629004 3228941.87002832 3228936.47941013\n", " 3228931.93877676 3228927.59926889 3228922.55478873 3228917.81631337]\n" ] } ], "source": [ "# IMPORT GEOMETRIES IN GEOPANDAS\n", "\n", "import pandas as pd\n", "import geopandas as gpd\n", "gdf = gpd.GeoDataFrame.from_file(input_file)\n", "\n", "# gdf.head()\n", "\n", "# create empty dataframe\n", "columns = ['x_UTM','y_UTM']\n", "dictionary_append = []\n", "geom2fit = pd.DataFrame(columns=columns)\n", "# print all rows x,y\n", "# fill a dataframe for each row\n", "for geometry in gdf.geometry:\n", " print(geometry.x)\n", " print(geometry.y)\n", " dictionary_append.append({'x_UTM':geometry.x, 'y_UTM':geometry.y})\n", "geom2fit = geom2fit.append(dictionary_append)\n", "\n", "print(geom2fit.x_UTM.values, geom2fit.y_UTM.values)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Linear fit & perpendicular" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[648601.47455578 648600.53612096 648599.76610303 648598.79900974\n", " 648597.94838364 648596.95458055 648596.09474052 648595.1762005\n", " 648594.27217673 648593.29520633 648592.6392054 648591.78651562\n", " 648590.97485142 648589.99758298 648589.14434109 648588.24569531\n", " 648587.57946282 648586.57906068 648585.86235498 648585.04563194\n", " 648584.30422735 648583.40016234 648582.56477513 648581.66570349\n", " 648580.90571977 648579.98210644 648579.20776285 648578.32362841\n", " 648577.46384525 648576.53901876 648575.59462599 648574.93642347\n", " 648574.05331062 648573.18304047 648572.24675381 648571.45037264\n", " 648570.54142569 648569.67218771 648568.33970007 648567.75427798\n", " 648566.90485212 648566.10936792 648565.02904074 648563.3972567\n", " 648562.87363833 648561.50311275 648560.98117157 648559.6440178 ]\n", "[648581.0571522126, 3229032.6283563185]\n", "center point x,y are : 648580.5592867916 3229032.8950612815\n", "slope of fit is: 5.643287879511387\n", "intercept of fit is: -431096.99035144487\n", "slope of perpendicular to fit is: -0.17720166352501993\n", "slope is: 1.3954152207670156\n", "perpendicular slope is: -0.17538110602788104\n", "check: 1.5707963267948966\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/anrossi/miniconda/envs/python3test/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:106: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance. In a future version, a new instance will always be created and returned. Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.\n", " warnings.warn(message, mplDeprecation, stacklevel=1)\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAALEAAAJdCAYAAAB9MJYVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJztvXmcVNWZ//9+eqFoQHYoELppChBkdWEQEhdcJi5fotHRCU5rTEYkk2U0k8RfnC8TjZo245jEjFEzIcZoYn/jqNGEII57uwSQESNdbM1STTeLFChr21DQ9Pn9cc+lbzfd9FbVt27V83697qtvPffce58uPpw+93zOOVeMMShKkMnxOwFF6SoqYiXwqIiVwKMiVgKPilgJPCpiJfCoiJXAoyJWAo+KWAk8eX4n0J0MHjzYFBcXN4l9+umn9O7d25+EOkGm5bty5cqPjTFDunQTY0zWbGeffbZpzptvvnlCLJ3JtHyB900X/121OaEEHhWxEnhUxErgURErgUdFrAQeFbESeFTESuBRESuBR0WsBB4VsRJ4VMRK4FERK4FHRawEHhWxEnhUxErgURErgadNEYtITxFZISKrRGSNiNxt42UiUikiq0XkcRHJt/ESEamw21IRmea51m22/BoR+ZYnfp2NNYjIdE+8WEQOiciHdvsvz7GzRSQqIptE5CERkWR9KUqwaE9NnAAuMsZMA84ALhORmUAZMAGYAhQA82z5KuACY8xU4F5gIYCITAZuAWYA04A5IjLOnrMauAZ4u4X7bzbGnGG3f/LEfwHMB8bZ7bL2/cpKptGmiO0sklr7Md9uxhizxDPFZAUw0pZfaozZa8svd+PA6cByY0ydMaYeeAu42p6zzhhT2d6kRWQ40NcYs8ze/7fAF9p7vpJZtKtNLCK5IvIhsAt41RjznudYPnAj8D8tnHoz8JLdXw2cLyKDRKQXcAVQ2I7bjxaRv4rIWyJyno2NALZ5ymyzMSULaddsZ2PMMeAMEekPvCAik40xq+3hR4G3jTHveM8RkQtxRHyuvcY6EbkfeBWoBVYB9W3c+iOgyBjziYicDfxRRCYBLbV/W1xoWUTm4zQ7CIfDlJeXNzleW1t7Qiyd0XxboKMzS4G7gO969v8I5DQrMxXYDJx2kuvcB3y9WawcmH6Sc8qB6cBwYL0nfj3wy7Zy19nO3c+bb75pdu58yixdOsq8+aaYpUtHmZ07nzp+nO6Y7SwiQ2wNjIgUAJcA60VkHnApcL0xpsFTvgh4HrjRGLOh2bWGespcA/y+HffOtfsRnAe4mDHmI+CgiMy0vRJfAv7U1u+i+MFrVFbOJ5GoBgyJRDWVlfOJx8uSdof2NCeGA09aMeUAzxhjFotIPVANLLO9W88bY+4B7gQGAY/aeL0xxu02+4OIDAKOAt8w9gFQRK4Gfg4MAV4UkQ+NMZcC5wP32HsdA/7JGLPHXutrwBM4PSMv0dj2VtKKx2hoqGsSaWioIxZbQDhckpQ7tCliY0wFcGYL8RbPNcbMo7G7rfmx81qJvwC80EL8D8AfWjnnfWByq4kracIuAO67DwYNgq9+1YkmEjVJu4M6dkqKGYoxsGIFHDzYGA2FipJ2BxWxkmLmsW9fAfv3w+jRTiQnpxeRSGnS7qAiVlLMJRw5cisAkQiEQqMYP35h0trDkGWrYir+sG2bs+jljTfuZvDgwUm/vtbESsoZPHgwn//851MiYFARK93ATTfdxKJFi1J2fRWxklIaGho4duxYSu+hIlZSRjxextatf0/v3nk88MDQpLp0XvTBTkkJ8XgZlZXzicXqSCSgf//dVFbOB0hqzwRoTaykiFhsAQ0NdVRVQU4OjBrVaDcnGxWxkhJcW3nzZigshB49msaTiYpYSQmurVxV1ejUeePJREWspIRIpBSRAi6/HC680Ikl22520Qc7JSW4D2833vgdYBehUBGRSGnSH+pARaykEGMuZv/+X3LVVVel9D7anFBSRmlpKf/wD//gTiFLGSpiJWVUVFQwevRoUr2ujYpYSQnGGKLRKJFIJOX3UhErSSceL+PPfy5k7969jB5dnjK72UUf7JSk4trN69Y5k0MjkYMps5tdtCZWkoprN0ci8N3vwtixqbObXbQmVpKKaysPGQL/5/+cGE8FWhMrScW1lZctgx07ToynAhWxklQikVIaGgq4805YvNiJpcpudlERK0klHC4hN/cu6uvdgT/hpM9ubo62iZWks2OH03T44hdXsWfPHsLh2Sm9n9bEStKJRqPk5eUxYcKEbrmfilhJOhUVFUyYMIEe7kj4FKPNCSVpxONlxGIL+PKXq6mrG26dutQv4K81sZIUXKcukahm8GAoKvrIOnWvpfzeKmIlKbhOXU0NPP007NuHXZf4sZTfW0WsJAXXkfvgA/jlL+HIEffIrpTfW0WsJAXvxNDevR3b2WFoyu+tIlaSQiRSSk5OL2IxGDMGRBynrpWXBiQVFbGSFMLhEk477ZdUVQmjRzeuQ+y8pyi1qIiVpNHQcBGJRC6XXfYLZs3aklKr2Yv2EytJ49RTT6W2tjblq2A2R0WsJJVQKNTt99TmhJI07rvvPu67775uv6/WxEqXce3mhQurKSws4OabR3Vbexi0Jla6iGs3HzxYzdatMGrUoaS/9rYtVMRKl/DazQ0Nzmu+Uj0xtDkqYqVLuHZzVZXz2V0rJZUTQ5ujIla6hGs3JxIQDsPIkU3j3YGKWOkSrt08Z44zei0vL/UTQ5ujIla6RDhcwvjxCwmFRgGSktfetoWKWOkyeXmX8Y//WMCBA3/sVrvZRUWsdJloNMr69et9cetARawkgWg0CsCUKVN8ub+KWOky0WiUQYMGMXz4cF/uryJWOk08XsayZcW8++6vGDWqll27/p8veejYCaVTuHZzQ0MdEybAsGGJlK9D3BoqYqVTuHYzwK23OjHXbtbeCSUQuLbykSPgfTlSd9rNLipipVO4tvLTT8OVVzZO0e9Ou9lFRax0Cu/s5r59nReQd7fd7KJtYqVTuO3eqqovE4nUEwqNStlrb9tCa2Kl0/Ttew3btjVw4YXf98VudlERK51m3bp1NDQ0MHXqVF/zUBErnWbAgAEsWLCAGTNm+JqHtomVTjN69Gh++MMf+p2G1sRKx3Ht5ieeEF57rahbJ4W2hIpY6RDexbS//W34yU+2dvvs5uaoiJUO4drNe/fC3r3+zG5ujopY6RCurRyLOZ/9mN3cHBWx0iG8i2mD+8JFf+xmFxWx0iG8dnP//jBwoH92s4t2sSkdwnXl/v7vb+czn/nIV7vZRUWsdJhwuIQvf9k/0TZHmxNKh9m1axfPPfcce/fu9TsVQEWsdIJ33nmH6667js2bN/udCqAiVjqA69T9+c/XIgKDBlX4nRKgbWKlnXgnhsZizsKBW7f+M716hXx9qAOtiZV24p0YGos5/cN+O3UuKmKlXbiO3KFDsGNHejh1LtqcUNpFKFREIlFNKARPPQXusmt+OnUuWhMr7cJ16nJy4NRTYdAg/506FxWx0i7cdYjfeWcwf/4zvqxD3BranFDaTThcQnn549TW1vKTn7zndzrH0ZpYaTfGGCoqKnxbwrU1VMRKu4nH43z88ce+z25ujopYaTd+L6bdGipipU1cu/nllz+HCITDlX6n1AQVsXJSvBND58yBJUvg44+/4/sMZy8qYuWkeO1mgJ4908dudlERKyfFtZWPHYMFC2DZsqbxdEBFrJwU11bevh2WLoUDB5rG0wEVsXJSXLvZHf8eiaSP3eyiIlZOims319T0IycHxo0rShu72UVtZ6VNwuES9ux5hvHjNzJ79lq/0zkBrYmVdnHKKadw3nnn+Z1Gi2hNrLSLp556yu8UWkVrYiXwqIiVVnHt5jvuEMaN68G6df/ld0otoiJWWsRrN2/YANu3HyUe/3Za2c0uKmKlRZrPbnYmhh5KK7vZRUWstIhrKxvTOEXfG08nVMRKi7i28scfw8GDjVP008ludlERKy3i2s319XDxxTBxYvrZzS7aT6y0iGsr5+cv4N/+rYZQqMj3dYhbQ0WstEo4XEKvXlfSp08fRMTvdFqlzeaEiPQUkRUiskpE1ojI3TZeJiKVIrJaRB4XkXwbLxGRCrstFZFpnmvdZsuvEZFveeLX2ViDiExvdv9/FZFN9l6XeuKX2dgmEbkjGV+GciLnnnsuc+fO9TuNk9KeNnECuMgYMw04A7hMRGYCZcAEYApQAMyz5auAC4wxU4F7gYUAIjIZuAWYAUwD5ojIOHvOauAa4G3vjUVkIjAXmARcBjwqIrkikgs8AlwOTASut2WVJHL06FHWrVtHcXGx36mclDZFbBxq7cd8uxljzBJ7zAArgJG2/FJjjLuE+HI3DpwOLDfG1Blj6oG3gKvtOeuMMS3NPrwKeNoYkzDGVAGbcP4TzAA2GWNixpgjwNO2rJJEKisrOXr0aNpN0W9Ou9rEtuZbCYwFHjHGvOc5lg/cCNzWwqk3Ay/Z/dVAqYgMAg4BVwDvt3HrETj/EVy22RjA1mbxc1rJfT4wHyAcDlNeXt7keG1t7QmxdKb78n2N11//OQCHD99Kefk64JIOX6U78m2XiI0xx4AzRKQ/8IKITDbGrLaHHwXeNsa84z1HRC7EEfG59hrrROR+4FWgFlgF1Ldx65aeJgwt/wUxreS+ENukmT59upk9e3aT4+Xl5TSPpTPdka9jOT9ILFZHbi4UFe0hJ+dBxo8/vcO9E92Rb4f6iY0x+4BynPYpInIXMAT4treciEwFHgOuMsZ84jn/18aYs4wx5wN7gI1t3HIbUOj5PBLYcZK4kgRcy3n6dJg3D/Lz02+Gs5f29E4MsTUwIlKA8zdlvYjMAy4FrjfGNHjKFwHPAzcaYzY0u9ZQT5lrgN+3cftFwFwRCYnIaGAcTvv7f4FxIjJaRHrgPPwtas8vrLSNay2feSZ4OybS0XKG9jUnhgNP2nZxDvCMMWaxiNQD1cAy24f4vDHmHuBOYBBOTwJAvTHG7Tb7g20THwW+4T4AisjVwM9xavUXReRDY8ylxpg1IvIMsBan6fEN27RBRL4JvAzkAo8bY9Z0+dtQAMdaPnCgmi1bnDETPXo0xtORNkVsjKkAzmwh3uK5xph5NHa3NT/W4vwWY8wLwAutHCsFTvA6jTFLgCWtJq50mkiklGefvZl//ucE990Hs2alr+UMOnZCaYFwuITaWqcdMWZMei2o3RJqOystUl1dQL9+/bjuur1pbTmD1sRKK0SjUaZMmZL2AgYVsdICxpjjIg4CKmKlCfF4GUuXFnPnnQeYNeuPaTmnrjnaJlaO43317dlnA3xEZeV8gLR9qAOtiRUPrlO3enXjEq7p7NS5aE2sHMd15F54AdascfqHvfF0RWti5TiuI9c4Rb9pPF1RESvHiURKqa8vYOvWxin66ezUuaiIleOEwyXk5Pwbx44Fw6lz0Tax0oTt250RrnPnrmHixGDM+NKaWGnC9ddfz9q1aznttNP8TqXdaE2sNCEvL4/TTz/d7zQ6hNbEShPuuOMO3njjDb/T6BAqYuU4e/bs4f7772flypV+p9IhVMQK4FjOTz/tPMjl5f04EGMmXLRNrBwfM1FZ6axHXFi4KxBjJly0JlaOj5mIxaBvXxg0KBhjJlxUxMrxsRH79jl2szsOPt3HTLhoc0IhFCoikajm3nvh6NGm8SCgNbFyfEFtcBZKgWCMmXBRESuEwyXs2nUrd95ZwO7dwRkz4aLNCQWAdet68c47h1iy5CB9+vTxO50OoTWxAjizmyORSOAEDCpixRKNRtN+HeLWUBFnOfF4GeXlRWzYsJ4BA94MlFPnom3iLMZ16nbvruP002HMmP2BcupctCbOYlynbsgQePhh+Oxng+XUuaiIs5jWHLmgOHUuKuIsxnXkvv99+NGPTowHBRVxFuM6dWvXQo5VQpCcOhcVcRYTDpcwaNCP2bPHmaIfNKfORXsnspx4fDwA11zzKrNmdfwVX+mA1sRZTkVFBUBglnFtCRVxlhOJRLjpppsIh8N+p9JptDmR5Vx55ZVceeWVfqfRJbQmzlLi8TL+8pdRLFokLFtWHEi72UVFnIW4dnNVVQ1XXQWLF1dTWTk/sEJWEWch3omhAEVFwbSbXVTEWYhrK8dizqTQUaOaxoOGijgL8S6mPXIk9OzZNB40VMRZiGs3x2LBWky7NbSLLQtxbeWbbrqNfv0+IRQaRSRSGji72UVFnKWEwyWUlgZTtM3R5kSWEovFiEajNDQ0+J1Kl1ERZykPPfQQM2fO9DuNpKAizlKi0SiTJk0iJyf4Egj+b6B0COfdzaNYufINhg1bF1iXzos+2GURrt388cd17N8PRUW1gZzd3BytibOI5nbzmDHBtptdVMRZhGsrT5jgTAwdP75pPKioiLMI11bu0wdmzoRevZrGg4qKOItw7eb/+R+orHRiQbabXVTEWUQ4XMKYMb/gwQeFN94I7uzm5mjvRJZx8OAMjhwxXHHFE8yadZPf6SQFrYmzjGg0ChDYZVxbQkWcZVRUVJCbmxu49zefDBVxFhGPl/H22z9hxIhj/PWvEzLCrQNtE2cNrlu3YMEh9u6FRKI6I9w60Jo4a3Dduh49wF0nJRPcOlARZw2JRA1btsAjj8CuXU3jQUdFnCWEQkVEo/Dcc3DsWNN40FERZwmRSClVVXkUFDQ2JzLBrQMVcdYQDpewc+dYxowJkZMjGePWgfZOZA3GGCord3HttTcxe/Yv/U4nqWhNnCXs3buXHj16BHod4tbQmjhLGDhwIB999FFGzG5ujtbEWUYmTAxtTub9RkoT4vEyli0rpqREuO66vhljNXvR5kQG41rNDQ11LF8OQ4YczBir2YvWxBmMazXX10NNDUQimWM1e1ERZzCupbx1K9TXN66AmQlWsxcVcQbjWspVVc7nSKRpPFNQEWcw7sTQ3Fxnmn5RUeZYzV70wS6DcR/e8vMXcMEFNYRCRYFeh7g1VMQZTjhcwtCh/4CI+J1KytDmRIazf/9+Bg4cyO9+9zu/U0kZKuIMJxqNsm/fPgYOHOh3KilDRZzhZOIU/eaoiDOYeLyMV1/9/+jdG7ZuPTcjLWfQB7uMxbWcN26sIxKBI0dqMtJyBq2JMxbXcj7nHLjkEieWiZYzaE2csbjW8g03tBzPJLQmzlBCoSJqa6Gu7sR4pqEizlAikVIWLcpnzpxGIWei5Qwq4owlHC5h9+6zCYdz6dUrs2Y3N0fbxBnMxo0HmT79cmbP/rPfqaQUrYkzlEQiQWVlZUabHC4q4gxl/fr11NfXZ+QU/eaoiDOUYcOG8fDDD3Puuef6nUrKURFnIPF4GbHYOUya9M8ZbTe76INdhuHazatW1TFkCAwbljmLabeG1sQZhms333MPPP64E8tUu9lFRZxhJBI1HDgAH3/cODHUjWcqKuIMIxQqOv4Ccq+IM9FudlERZxiRSClbtuTbfSeWqXazi4o4wwiHS9iz51z69s1h0KDMefXtydDeiQyktPRJbrllC+edd57fqXQLKuIMpLCwkMLCQr/T6Da0OZFhxONxHnroIbZv3+53Kt1GmyIWkZ4iskJEVonIGhG528bLRKRSRFaLyOMikm/jJSJSYbelIjLNc63bbPk1IvItT3ygiLwqIhvtzwE2PltE9ovIh3a703POZfb+m0TkjmR+KUElHi/jiSemcNttt7FkyfSMd+pc2lMTJ4CLjDHTgDOAy0RkJlAGTACmAAXAPFu+CrjAGDMVuBdYCCAik4FbgBnANGCOiIyz59wBvG6MGQe8bj+7vGOMOcNu99hr5QKPAJcDE4HrRWRiZ76ATMF16jZs2A3AyJE7qaycnxVCblPExqHWfsy3mzHGLLHHDLACGGnLLzXG7LXll7tx4HRguTGmzhhTD7wFXG2PXQU8afefBL7QRlozgE3GmJgx5gjwtL1G1uI6dZs3w6mnQkFB5jt1Lu16sLM130pgLPCIMeY9z7F84EbgthZOvRl4ye6vBkpFZBBwCLgCeN8eCxtjPgIwxnwkIkM915glIquAHcB3jTFrgBHAVk+ZbcA5reQ+H5gPEA6HKS8vb3K8trb2hFg603q+jiNXVdW4DjE4Tp2fv193fL/tErEx5hhwhoj0B14QkcnGmNX28KPA28aYd7zniMiFOCI+115jnYjcD7wK1AKrgPo2bv0BMMoYUysiVwB/BMYBLa2OZ1rJfSG2STN9+nQze/bsJsfLy8tpHktnWst32bIiDh6sZscOuOCCxngoVMSsWSeW7y664/vtUO+EMWYfUA5cBiAidwFDgG97y4nIVOAx4CpjzCee839tjDnLGHM+sAfYaA/FRWS4PXc4sMuWP+A2ZYwxS4B8ERmMU/N6+5BG4tTUWUskUkrPnr1YtAiuu86JZbpT59Ke3okhtgZGRAqAS4D1IjIPuBS43hjT4ClfBDwP3GiM2dDsWkM9Za4Bfm8PLQJusvs3AX+y5YaJXZNURGbYfD8B/hcYJyKjRaQHMNdeI2sJh0sYP34h/fuP4pRTMntiaHPa05wYDjxp28U5wDPGmMUiUg9UA8uszp63vQd3AoOAR2283hgz3V7rD7ZNfBT4hucB8N+BZ0TkZpzGna1LuBb4mr3XIWCufZCsF5FvAi8DucDjtq2c1bz4YoKami/zgx/8wO9UupU2RWyMqQDObCHe4rnGmHk0drc1P9aiD2qbHBe3EH8YeLiVc5YAS1pNPAt5+umn+eSTT7JOxOrYZRDRaDQrJoY2R0WcIezevZudO3dmxRT95qiIM4B4vIynn54MQG7u/Vnh0nnRUWwBx7WbP/mkjlNOgcLCXRk/MbQ5WhMHHNdunj0b/vQnGDgwe+xmFxVxwPFOAPW+5SuTJ4Y2R0UccEKhIhoa4Otfh1dfbRrPFlTEAScSKWXnzp6sWwdHjzqxbLGbXVTEASccLuHw4a8CzuzmbLKbXbR3IgPYunUAIsKXvlRLr169/E6n29GaOAOoqKhg7NixWSlg0Jo4I4hEIhQXF/udhm+oiDOABx54wO8UfEWbEwHn2LFjOKNTsxcVccB56KGvMXBgHs8+KyxbVpx14yZARRxo4vEyli37DYcPNzBwICQS1VkzTd+LijjAxGIL2LSpnuJiyM11Ytk2bgJUxIEmkag5YYq+G88mVMQB5tNPR7BvX9PFtCG7xk2AijjQFBZ+jy9+MY9p0xpj2TZuAlTEgWbq1G/yn//5BJMnjwKya5q+FzU7AkxNTQ1Dh/4ds2Zll2ibk701ce0uOLTP7yy6xBe+8AWuuiqr11EEsrkmfvvHsOKXzAwNhu1nwdCJEJ7kbIPGQV4PvzM8KfX19axdu5ZvfvObfqfiO9kr4ql/D32Hs7/iDXru3wab34AGu75hTj4MPg3CVthDrbj7ntp0DpCPLF/+IIlEgry8n7Bs2XNEIqVZ1xZ2yV4Rj5wOI6ezrv5MwrNnQ/0R+GQjxNfCrjUQXwPVyyD6bOM5PftBeLKttSc64h56OvTs262px+NlvP76vwFO95rr1EH2zHD2kr0ibk5ej8bmxPGl4HDazbvWQXw17FrriHvV03DkYGOZ/kVNxR2eDAPHQG5qvt5YbAGbNx8hJwfcEZiuU6ciVk6koD+MmuVsLsbAvhor6tW29l4LG14Gc8wpkxuCIaedKO4+4S43SRKJGs4/H4YNgx49msazERVxZxCBAaOcbfzljfH6BOyubKyxd62FzW/Cqt83likYaNvZngfJoadDj97tvn0oVMTYsdWMHXtiPBtRESeTvBAMn+psXur2NIrarbn/+hQc/dQWEBhQ7BG12yQZDTm5xy/zx79u597yOvYcfpjespurxz7JxWPeArLTqXNREXcHvQbC6POczaWhAfZVnyjuyiXgrlmeVwBDxkN4MtH6ESxalQdHR2LoR60Zyu833kpBgXDh6CrtnVB8ICfHqWkHjobT5zTGjx6C3esb29nx1bDxFaZ8uovHc4Fc2G36sr6hiEpTSDz2WWbN/ikMnODbr+I3KuJ0I78ATj3T2TxMv+P/cVrOVibIViZIDeNztlKS8zoF9S/Bwp+D5Dg9Im7XX3iSs9+/2PkPk8GoiANCqP8wlu7rx1ImH4/l0MCMvvt4+qpTGh8md0Zh7SKOv0wqv7fz4Hhc3La93WugP79IClARB4TbLx3Pvz4f5dDRY8djofx85l5+EUwaAZM876888qnt217TKO51i+GD3zaW6TPsREdyyHjn4TRgqIgDwhfOHAHAvX9azieHQ/TL380NU15k1qlXA80e6Hr0Pu5IHscYqI07gvaK+72FcCzhlJFcGDyuWfffRMfMSRO7vSVUxAFi1qnl/GT2zTiv23aorHRe2Npmz4QInDLM2cZ63vFzrB72bG4q7u3vw5rnG8v0OMU2SSY1FXdB/yT+dp1HRRwgYrEFrFyZYPt2+PznHV122W7OzXOaEUPGw+RrGuOHDzhNkl1rGntK1jwPK3/TWKbvyGZNkom+jABUEQeIRKKGl1+GDz+EK69sGk86PftC0TnO5mIMHNhha2xbc8fXnnQE4MBPjsH+cSkdAagiDhChUBGbN1f7NzFUBPqNcLbTPtcYP2EE4FqoWQ7RZ5kKEL3XGQE4+gL44u+SnpaKOEAUFt5NTc2XmTGjMZYWdvNJRgD+9eUyzjy1h9McyeuZmtun5KpKSti/fzr19TB+/CBgD6FQUXrbzQX92d9/EsyYndLbqIgDxIYNzvver7nmjax86WJrZLYfmWFcffXVLFq0iIkTJ/qdSlqhIg4Yp5xyCnl5+gfUi4o4QNxyyy289dZbfqeRdqiIA8LGjb/iscceY/v2H2TtOsStoSIOAPF4Ga+88s+AswJmtq5D3Boq4gDgrEPsjJdwjY5sXIe4NVTEASCRqCEWg969YejQpnFFRRwIQqEijIFJk5oOP8jW2c3N0b6aABCJlPLd786noaHueCwt7OY0QWviABAOlzB+/EJCoexeh7g1VMQB4JVXXuHKKx9i6NDXgTeYNWuLCtiDijgArFy5khUrVjB48GC/U0lLVMQBIBqNUlRURL9+/fxOJS1REac58XgZ7733HCNG1LBsWTHwmt8ppR0q4jQmHi9j9epbqK4+etypgx+rU9fFDNZVAAAcjElEQVQMFXEaE4st4ODBQ3zmMzBlihtNqFPXDO0nTmMSiRr694d77jkxrjSiNXEaEwoVcfRoy3GlERVxGhOJlPL97+dw++3eaEidumaoiNOYcLiEmpp+DBrUG9epg++q0dEMFXEas2fPHj76aC8XX3wXs2c3MGvWFuASv9NKO1TEaUw0GgVgSmPXhNICKuI0xhWxTs8/OSriNGbatGl85zvfYfjw4X6nktaoiNOUeLyMvLwbmTPnpyxfPlpdupOgZkcaEo+XsW7dLVRVHaKwEKDxtbcwwsfM0hOtidOQWGwBH310iH/8R3jJWUNbJ4aeBBVxGuJODAUYM6ZpXDkRFXEaEgoVHRex+wJyN66ciIo4DYlESqmqyuXUU6FXLyemE0NbR0WchoTDJWzbFmbMmAJ0YmjbaO9EmvLQQ7+mV69enH/++X6nkvaoiNOUyy67zO8UAoM2J9KQNWvW8Morr1BfX+93KoFARZxmxONllJZ+hjlzLmXZsjHq1LUDbU6kEfF4GZWV89m4sY7iYjh2rOa4U6cPda2jNXEaEYstsM6cLuHaEVTEaUQiUcP+/bBnD01euKhO3clREacRXqfOK2J16k6OijiNiERKmTKlgF/9ylmLGNSpaw/6YJdGuA9vvXsvIJGoSf83hqYJKuI0o6wszqRJv+TSSy/1O5XAoM2JNKKhoYHvf//7vOQOIlbahYo4jYjFYtTV1enE0A6iIk4jdIp+51ARpwnxeBkvvvgVRODgwb9Tu7kD6INdGuDazdu31zFiBOTkbFW7uQOoiNMA127+3vfg0CEn5trNKuK20eZEGuC1lQsKWo4rraMiTgNCoSKqquDuu6GmpmlcaRsVcRoQiZSyfn0Pysshx/6LqN3cflTEaUA4XMK+fRcRCgnDh6MTQzuIPtilCZs3H2Xq1OlcfPEKv1MJHFoTpwkVFRVqcnQSFXEaUFtbS2FhIdOnT/c7lUCizYk0oE+fPqxcudLvNAKL1sQ+E4+XsWxZMeXlOSxbVqx2cyfQmthHXLv5Zz+r4+OP4Z57qtVu7gRt1sQi0lNEVojIKhFZIyJ323iZiFSKyGoReVxE8m28REQq7LZURKZ5rnWbLb9GRL7liQ8UkVdFZKP9OcDGRUQeEpFN9npnec65yZbfKCI3JfNL6S5cu3n16hPtZqX9tKc5kQAuMsZMA84ALhORmUAZMAGYAhQA82z5KuACY8xU4F5gIYCITAZuAWYA04A5IjLOnnMH8LoxZhzwuv0McDkwzm7zgV/Yaw0E7gLOsde7yxV+kEgkajh2DKqrYfTopnGl/bQpYuNQaz/m280YY5bYYwZYAYy05ZcaY/ba8svdOHA6sNwYU2eMqQfeAq62x64CnrT7TwJf8MR/a2+zHOgvIsOBS4FXjTF77L1eBQK3eFkoVMT27XDkiM5u7grtahOLSC6wEhgLPGKMec9zLB+4EbithVNvBty5NquBUhEZBBwCrgDet8fCxpiPAIwxH4nIUBsfAWz1XG+bjbUWbyn3+Ti1OOFwmPLy8ibHa2trT4h1HzcQi90P1HtEHCKRuKHVnPzNt+N0R77tErEx5hhwhoj0B14QkcnGmNX28KPA28aYd7zniMiFOCI+115jnYjcj1Nr1gKrgLZWzJOW0jlJvKXcF2KbNNOnTzezZ89ucry8vJzmse5jNrt3H+Pcc/+TUaMOEQqNanN2s7/5dpzuyLdDvRPGmH0iUo7zp3u1iNwFDAG+6i0nIlOBx4DLjTGfeM7/NfBrW+Y+nBoUIC4iw20tPBzYZePbgELPpUcCO2x8drN4eUd+l3Thuut+xHXX/cjvNAJNe3onhtgaGBEpwHm58HoRmYfTNr3eGNPgKV8EPA/caIzZ0OxaQz1lrgF+bw8tAtwehpuAP3niX7K9FDOB/bbZ8TLwOREZYB/oPmdjgePTTz/1O4XA056aeDjwpG0X5wDPGGMWi0g9UA0sExGA540x9wB3AoOAR2283hjj+ql/sG3io8A3PA+A/w48IyI3AzXAdTa+BKftvAmoA74CYIzZIyL3Av9ry91jjNnTqW/AR2pra+nXrx8//elPue22lh4plPbQpoiNMRXAmS3EWzzXGDOPxu625sfOayX+CXBxC3EDfKOVcx4HHm818QCwevVqGhoaKPa+IknpMGo7+4hO0U8OKmIfWb78WXr1ErZsGaPjJrqAitgn4vEyPvjgdYqLDTk5kEg44yZUyB1HBwD5RCy2gEsvbSA3tzGm0/Q7h4rYJxKJGlp6y5eOm+g42pzwidraEWzdCg0NTeM6bqLjqIh94r33zuVLXwKv16HT9DuHitgntmzJY/jwAQwePAp9f3PX0DaxT0SjUaZNO4dZs3RB7a6iNbEPHD16lHXr1qnJkSRUxD6wYcMGjhw5oivCJwkVsQ+EQkv50Y+GUFBwozp1SUDbxN1MPF7Gjh3fYubMOqDRqQOd4dxZtCbuZmKxBbz7bh3r1zfGdIZz11ARdzOJRA0PPQTPPntiXOkcKuJu5ujRkcTjTWc3gzp1XUFF3M0cPvwVoKmI1anrGiribmbbtmEATJgwAnXqkoP2TnQz0WiUfv368YUvbMXOQVS6iNbE3cwDDzzA0qVLVcBJREXczfTu3ZuJEyf6nUZGoSLuRnbu3Mkdd9zBhg0b2i6stBsVcTfy2ms/5v777+fll8er3ZxEVMTdRDxexttvPwQ4y7jqxNDkoSLuJmKxBWzefJRwGPr0cWJqNycHFXE3kUjUUFV1olOndnPXURF3E3l5hezb13RFeFC7ORmo2dFNjBt3Hy+8cAuJxKHjMbWbk4PWxN1EOFzChAm/4pRTdGJostGauJv4+c9/zqpVq3jssS1+p5JxqIi7iZdeeont27f7nUZGos2JbiIajers5hShIu4G9u7dy7Zt23R2c4pQEXcDb731U7v3PbWbU4CKOMXE42Vs3vwAxcVqN6cKFXGKicUWcPbZCX7zGxgyxImp3ZxcVMQppjVbWe3m5KEiTjH5+YWUlMAf/tA0rnZz8lARp5j8/G+xYwf06NEYU7s5uaiIU8zOnWMAmDBhGGo3pwZ17FKM+666uXM3cMopp/icTWaiNXGKqaioYPTo0SrgFKI1cYqZMWOGzm5OMSriFPOd73zH7xQyHm1OpJAtW35Defkoystz1G5OISriFBGPl7Fw4T9x8cU1bN9u1G5OISriFOHMbj5Cbi4Mc9YQVLs5RaiIU0QiUUMsBsXFNHl/s9rNyUdFnCJCoaIWp+ir3Zx8VMQpol+/O/jkE11MuzvQLrYUMXz4F7n77neJRN4AdhIKFRGJlKrdnAJUxCliwIAB3HnnU36nkRVocyJFVFRUsHXrVr/TyAq0Jk4Rt9xyC3369OH111/3O5WMR2viFNDQ0MDq1at1in43oSJOAStWPEhdXR35+f+pdnM3oCJOMvF4Ga+99n8BGDNGZzd3ByriJOPazSKOWwdqN6cafbBLMolEDZdeCuPGQc+eTeNKalARJ5lQqIhhw6qPD/rxxpXUoM2JJDNs2J0sXpzPzp2NMbWbU4uKOMl8/PEUfvKTo1RVDUFnN3cP2pxIMu7s5muv/Qvjxo3zOZvsQGviJBONRikoKCDSfAymkjJUxEkmGo0yefJkcr0j4ZWUoiJOMroifPejIk4i8XgZv/tdiCuueFzt5m5EH+ySRDxeRmXlfHr0qGPQoEa7GdCeiRSjNXGSiMUW8O67dfzmN3DsmBNTu7l7UBEniUSihrffhsWLdXZzd6MiThKhUBGxmM5u9gMVcZIoKrqH6uqmLyBXu7l7UBEniYMHZ3DkCJx22iDUbu5etHciSWzbto3evXtzzTWvcNZZZ/mdTlahIk4Sl1xyCQcOHPA7jaxERZxEcnK0deYH+q13kXi8jGXLijn/fOF73xukLp0PqIi7gOvS7dtXzTvvwMcf79FJoT6gIu4CsdgCGhrqqKpyPo8erS6dH6iIu4DrxsVizucxY5rGle5BRdwFXDcuFoOCAgiHm8aV7kFF3AUikVJycnpxyikwcybk5KhL5wfaxdYFXDfun/5pAYlEja5B7BMq4i4ydOg/qGh9RpsTXeTll19m5MiRx2c5K92PiriLRKNRtm/fzogRI/xOJWtREXeRaDTKiBEjGDhwoN+pZC0q4i4Qj5exfPnTjBixXSeG+oiKuJPE42WsWXMLW7YcJRLRdYj9REXcSWKxBdTVHWLOHDj7bCemlrM/aBdbJ0kkaujTB2699cS40r1oTdxJQqEi9u+H+voT40r3oiLuJJFIKfffn8PXvtYYU8vZH1TEnSQcLmHr1oGMHt0LnRjqL9om7iT79+9n27aP+frX72P27H/1O52sRmviTrJ69WoAXQEzDVARdxJ3rMTUqVN9zkRREXeSz372s/zHf/wHhYWFfqeS9aiIO0E8XkZt7ef5m7/5HsuXj1aXzmdUxB0kHi9j/fpb+MtfqjlwwKjdnAaoiDtILLaAePwQt98Ob77pxNRu9pc2RSwiPUVkhYisEpE1InK3jZeJSKWIrBaRx0Uk38ZLRKTCbktFZJrnWv9ir7FaRH4vIj1t/CIR+cDGnxSRPBufLSL7ReRDu93pudZl9v6bROSOZH8xrZFI1LB5s7PvXcZV7Wb/aE9NnAAuMsZMA84ALhORmUAZMAGYAhQA82z5KuACY8xU4F5gIYCIjABuBaYbYyYDucBcEckBngTm2ng1cJPn/u8YY86w2z32WrnAI8DlwETgehGZ2NkvoSOEQkVN1pnwxhV/aFPExqHWfsy3mzHGLLHHDLACGGnLLzXG7LXll7txSx5QYGvaXsAOYBCQMMZssGVeBf6ujbRmAJuMMTFjzBHgaeCqtn6XZBCJlBKL5RIOQ58+TkztZn9pl2Nna76VwFjgEWPMe55j+cCNwG0tnHoz8BKAMWa7iPwYqAEOAa8YY14REQHyRWS6MeZ94FrA2281S0RW4Qj+u8aYNcAIYKunzDbgnFZynw/MBwiHw5SXlzc5Xltbe0Ls5IwgFhvI6NH7gaPAUBoa5rFu3QjWrevIdTpHx/P1l27J1xjT7g3oD7wJTPbEfgX8rIWyFwLrgEH28wDgDWAITm3+R+AGe2wW8A5Ojf5D4K823hfoY/evADba/euAxzz3uhH4eVv5n3322aY5b7755gmxtvjwww/NBx980OHzkkFn8vWTtvIF3jcd0GBLW4d6J4wx+4By4DIAEbnLivLb3nIiMhV4DLjKGPOJDV8CVBljdhtjjgLPA5+x111mjDnPGDMDeBvYaOMHjG3KGGOW4NTYg3FqXm9tPRKnpu4Wpk2bxplnntldt1PaoD29E0NEpL/dL8AR43oRmQdcClxvjGnwlC/CEeiNprGdC04zYqaI9LJNiItxampEZKj9GQK+B/yX/TzMlkVEZth8PwH+FxgnIqNFpAcwF1jU+a+h/axcuZInnniCw4cPd8ftlHbQnjbxcOBJ2y7OAZ4xxiwWkXqcnoRlVmfPG6f34E6ch7VHbbzeGDPdGPOeiDwHfADUA3/F9lwAt4vIHHv9Xxhj3rDxa4Gv2XsdwunBMEC9iHwTeBmnl+Nx47SVU0o8XsaDD36d//7vA4wZcxennXafDr1MA9oUsTGmAjjhb6cxpsVzjTHzaOxua37sLuCuFuK3A7e3EH8YeLiVay0Blpws92TirkW8aVMdRUVw7FiNvjE0TVDHrp24axHHYo1LuKpTlx6oiNtJIlHDwYOwe3dTk0OdOv9REbeTUKiILVucfa/drE6d/+j0pHYSiZRy9Oh8nn++jl69nJg6demBiriduA9vsZiuRZxuqIg7wCOPVDJ+fCklJSrcdELbxO3EGMPPfvYzli9f7ncqSjNUxO2kurqagwcP6uzmNERF3E4qKioAnaKfjqiI20E8XsaLLzrj9D/99Is6ny7N0Ae7NnDt5gMH6hg1CvLytqrdnGZoTdwGrt381a/Cb37jxNRuTi9UxG3gtZWdQXknxhV/URG3gTsx9NZbYcOGpnElPVARt0EkUsrGjT2IRiEUcmJqN6cXKuI2CIdL2LfvEnr0gJEj0XWI0xDtnWgHmzfXM2nSmVx88Qd+p6K0gNbE7SAajeoSrmmMirgNEokEZ511Fueff77fqSitoM2JNgiFQixevNjvNJSToDVxGzQ0NLRdSPEVFXELxONlLFtWTHl5Dldf3Y+zzhrrd0rKSVARN8MdK5FIVAOGDRtqgSod9JPGqIib4Y6VAGhogC1bIBJp0LESaYyKuBneMREffQSHDztT9HWsRPqiIm6Gd0xELOb8jER0rEQ6oyJuRiRSSk6OMyd/yBCYMwcikQIdK5HGaD9xM7xT8ydMqGHaNJ2an+6oiFsgHC4hHC5hy5YtFBYWkpub63dKyknQ5kQr1NXVEYlEKC3VZkS6oyJuhbVr12KMYfLkyX6norSBirgZrlv3zDN/A8Dw4VU+Z6S0hYrYg9etq6pyZnIcPvx9devSHBWxB69bt3kzFBeDyCF169Ic7Z3w4HXl5s6FI0dOjCvph4rYQyhUZAf+wIwZTeNK+qLNCQ+uWxePw4cfOjWxzmxOf1TEHsLhEsaPX8hf/jKAf/kXqK8fqTObA4A2J5oRDpdw8OBrhMMvccUVW9s+QfEdrYlbQGc3BwsVMU2nI7377ijWrKnQdYgDhDYneI3KygeP9w/HYjUcPgyjRtX6nJfSXlTEPHZcwADDhsEjj0Bxcbe9cVfpIipidjX51KMHTJwIsN2XbJSOo21ihjb59MorsHy5GhxBQkXMvOPTkQCefBJefjlXDY4AoSLmEsaPX0goNIpDh2DHDjjnnKvV4AgQKmIcg2PWrC0MHOi8aHHWrBt8zkjpCCpiD/quumCiIvZQWVlJnz59KC4u9jsVpQOoiD088MADbN68mZwc/VqChP5r8dpxy3n58tEY86rfCSkdJKtF7Myd+zGJRDV79hh+8INq/vzneTqnLmBktYiduXMJADZtcoyO2trDOqcuYGS1iL1z57yLB+qcumCR1SJuvgLm4MHQt69azkEjq0XsWMvOa0JjMacW1jl1wSOrRexYy9+lR48icnJgwoS+OqcugOhQTC7hM5/5YZOXjyvBIqtrYiUzUBHjOHWXX345xhi/U1E6gYoYePvtt6mpqUFE/E5F6QRZ2SaOx8uIxRbY/uChfPjhMc499xK/01I6SdaJ2F2+1Z0cWlsbZ9s2GD36mM+ZKZ0l65oT3uVbAarsGtoDB5b7k5DSZbJOxM0t5ZwcmD4diop2+5SR0lWyTsTNLeVJk+CBB6CwcJRPGSldJetE7H3ZIsCxY2o1B52sE7G7fGsoNApj4JprhD/96VK1mgNM1okYGmc3jx1bw4EDhrFj/9bvlJQukJUidolGo4DObg46WS1id4q+vnAx2GS1iKPRKOFwmP79+/uditIFss6x83L55ZergDOArBbxDTfcwMiRI/1OQ+kiWSvi/fv3s3//fh1+mQFkbZt48eLFjBo1ii1btviditJFslbEdXV1jBw5ksLCQr9TUbpI1or4lltuYevWreTlZW2LKmPIWhErmYOKWAk8KmIl8KiIlcCjIlYCj4pYCTwqYiXwqIiVwKMiVgKPilgJPCpiJfCoiJXAoyJWAo+KWAk8KmIl8KiIlcCjIlYCj4pYCTwqYiXwqIiVwKMiVgJPmyIWkZ4iskJEVonIGhG528bLRKRSRFaLyOMikm/jJSJSYbelIjLNc61/sddYLSK/F5GeNn6RiHxg40+KSJ6Ni4g8JCKb7PXO8lzrJhHZaLebkv3FKMGhPTVxArjIGDMNOAO4TERmAmXABGAKUADMs+WrgAuMMVOBe4GFACIyArgVmG6MmQzkAnNFJAd4Ephr49WAK8rLgXF2mw/8wl5rIHAXcA4wA7hLRAZ09ktQgk2bIjYOtfZjvt2MMWaJPWaAFcBIW36pMWavLb/cjVvygAJb0/YCdgCDgIQxxn278qvA39n9q4Df2tssB/qLyHDgUuBVY8wee69Xgcs68wUowaddbWIRyRWRD4FdOOJ5z3MsH7gR+J8WTr0ZeAnAGLMd+DFQA3wE7DfGvAJ8DOSLyHR7zrWAuyzPCGCr53rbbKy1uJKFtGv5G2PMMeAMEekPvCAik40xq+3hR4G3jTHveM8RkQtxRHyu/TwAp2YdDewDnhWRG4wxT4nIXOBBEQkBrwD17mVaSuck8RMQkfk4TRGAWhGpbFZkMM5/pKCQafl2+bVVHVrDyRizT0TKcf50rxaRu4AhwFe95URkKvAYcLkx5hMbvgSoMsbstmWeBz4DPGWMWQacZ+OfA06z52yjsVYGp2myw8ZnN4uXt5LzQmy7vCVE5H1jzPTWjqcbmu+JtKd3YoitgRGRAhwxrheReTht0+uNMQ2e8kXA88CNnnYuOM2ImSLSS5w3gV8MrLPnDLU/Q8D3gP+y5ywCvmR7KWbiNEE+Al4GPiciA2wN/zkbU7KQ9tTEw4EnRSQXR/TPGGMWi0g9Tk/CMvt2+ueNMfcAd+I8rD1q4/XGmOnGmPdE5DngA5zmwl9prCFvF5E59vq/MMa8YeNLgCuATUAd8BUAY8weEbkX+F9b7h5jzJ5OfwtKoJFsX2RaRObbJkcg0HxbuEe2i1gJPmo7K4En0CIWkf4i8pyIrBeRdSIyy3PsuyJiRGSw/dxPRP7ssc+/4il7TEQ+tNsiT1xEpFRENtjr3+qJt2iH+5zvxda+/1BE3hWRsTYeEpH/tvm+JyLF3ZxvkYi8Yq+x1r2/iIy2+Wy0+fXoVL7GmMBuOHb1PLvfA+hv9wtxeiuqgcE29n+B++3+EGAP0MN+rm3l+l8Bfgvk2M9D7c8rcEwcAWYC76VJvhuA0+3+14EnPPv/ZffnAv/dzfmWA39r9/sAvez+MzjDDcDpkfpaZ/L1XYhdEHBfnHEa0sKx54BpwBbPl/yvOMaM4BgumzzibE0UK4CxLcR/idO16H6uBIanQb6VwDme8++z+y8Ds+x+Ho75cEIeqcgXmAi828I1xOaRZz/PAl7uTL5Bbk5EgN3Ab0TkryLymIj0FpErge3GmFXNyj8MnI5jlkSB20xj/3ZPEXlfRJaLyBc854wBvmiPvSQi42y8M7Z3d+Q7D1giIttwhgL8e/N8jTH1wH6cbtDuyPc0YJ+IPG+v84Dtrh0E7LP5QNPvsGP5+l2jdqEmno7T3+zWPP8JPAC8B/SzsS001hTXAg/i1ABjcWqZvvbYqfZnxJ4zxq3xgO/Y/WuAd+z+i8C5nlxeB85Og3yf91z/duAxu78GGOnJZTMwqDvytfH9Ntc84A84wxGGAJs89ysEop3J13cxdkHEw4Atns/nWTHtsl/uFvuPUGPLvgic5yn/BjCjhes+AVxr99cDxXZfcBxD6FxzIqX5WlFs9sSLgLV2vzPNiaTki/PMUO6J3wg8gjYnwBizE9gqIuNt6GLgA2PMUGNMsTGmGOdP1Fm2bI0tg4iEgfFATBzrOmTjg4HPAmvtNf8IXGT3L8B5cILW7XA/890L9BMRd9zJ32JtfZuvO0b7WuANYxWS6nxxXNUBIjLEXucinP9cBnjT5oPN70+dytfvGrWLtfEZwPtABY7gBjQ7voXGP3en4oyQiwKrgRts/DM2tsr+vNlzfn+cGiYKLAOmeWrlR3D+zEVxBvqnQ75Xe46VAxEb7wk8i/OwtcKNd0e+9tjf2mtEcf5yuL0WEZvPJptfqDP5qmOnBJ7ANicUxUVFrAQeFbESeFTESuBRESuBR0WsBB4VsRJ4VMRK4Pn/Aft0I1T8Dnx3AAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "from math import atan, degrees\n", "\n", "\n", "# x of geometry in UTM\n", "x_points = geom2fit.x_UTM.values*1.0\n", "# y of geometry in UTM\n", "y_points = geom2fit.y_UTM.values*1.0\n", "\n", "print(x_points)\n", "# center point\n", "center_point = [np.mean(x_points), np.mean(y_points)]\n", "print(center_point)\n", "\n", "# center_point_x = np.mean(x_points)\n", "center_point_x = (max(x_points) + min(x_points))/2\n", "# center_point_y = np.mean(y_points)\n", "center_point_y = (max(y_points) + min(y_points))/2\n", "\n", "print('center point x,y are :', center_point_x, center_point_y)\n", "# fit of points\n", "fit = np.polyfit(x_points,y_points,1)\n", "m, b = np.polyfit(x_points, y_points, 1)\n", "fit_fn = np.poly1d(fit)\n", "# print(fit)\n", "print('slope of fit is: ', fit[0])\n", "print('intercept of fit is: ', fit[1])\n", "### perpendicular line passing through the center\n", "print('slope of perpendicular to fit is: ', -1/fit[0])\n", "m_p = (-1/fit[0])\n", "\n", "# m_angle = degrees(atan(fit[0]))\n", "# m_p_angle = degrees(atan(m_p))\n", "m_angle = atan(fit[0])\n", "m_p_angle = atan(m_p)\n", "print('slope is: ', m_angle)\n", "print('perpendicular slope is: ', m_p_angle) #in radians\n", "print('check: ', m_angle - m_p_angle) # in radians\n", "\n", "# intercept of perpendicular line\n", "b_c = center_point_y -m_p*center_point_x\n", "y_p_points = m_p*x_points + b_c\n", "\n", "# set figure sixe\n", "fig = plt.gcf()\n", "fig.set_size_inches(4.5, 10.5)\n", "\n", "# perpendicular_fit = ((-1/m)*x_points, y_points)\n", "# points and fit\n", "ax1 = plt.plot(x_points,y_points, 'yo', x_points, fit_fn(x_points), '--k')\n", "# plot center point\n", "ax2 = plt.plot(center_point_x, center_point_y, 'o')\n", "# plot perpendicular line\n", "ax3 = plt.plot(x_points, y_p_points)\n", "# plot perpendicular fit\n", "# plt.plot(x_points,(-1/m)*x_points, 'ro')\n", "plt.grid('on')\n", "plt.axes().set_aspect(1) # it's turned off, but lines are perpendicular, uncomment to check ;)\n", "plt.show" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# line azimuth\n", "# see also for a more general case \n", "# https://github.com/marcoalopez/line_azimuth/blob/master/get_lineazi.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Apply transform" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "slope of perpendicular to fit is: -0.17720166352501993\n", "delta is: [-40, -30, -20, -10, 10, 20, 30, 40]\n", "0\n", "-40\n", "48\n", "-30\n", "96\n", "-20\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/anrossi/miniconda/envs/python3test/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:106: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance. In a future version, a new instance will always be created and returned. Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.\n", " warnings.warn(message, mplDeprecation, stacklevel=1)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "144\n", "-10\n", "192\n", "10\n", "240\n", "20\n", "288\n", "30\n", "336\n", "40\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUIAAAIfCAYAAAAFVdUoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzsnXt4nFW1/z8rtya0xhYopYBp4Iii7bGxRWg8KOWiiHpEETwgrcjVo0jlWG6FH3CkxwoHULkIR5QeUAqI3NUCIrYgpqWAphzCRUCStNICQksNNNMks39/zDswDWkzmXfPzJp3r8/zzJPMnvfds/fMO9+91tpr71eccxiGYYRMVbkbYBiGUW5MCA3DCB4TQsMwgseE0DCM4DEhNAwjeEwIDcMIHhNCwzCCx4TQMIzgMSE0DCN4TAgNwwiemnI3oJRsv/32rrm5Oa9j33jjDUaPHl3cBo0QjW0Ca9dIsXaNjELa9dhjj/3dOTc+7xOcc8E8pk+f7vJlyZIleR9bKjS2yTlr10ixdo2MQtoFPOpGoA3mGhuGETwmhIZhBI8JoWEYwWNCaBhG8JgQGoYRPCaEhmEEjwmhYRjBY0JoGEbwmBAahhE8JoSGYQSPCaFhGMFjQmgYRvCYEBqGETwmhIZhBI8JoWEYwWNCaBhG8JgQGoYRPCaEhmEEjwmhYRjBY0JoeGXZMvje9zJ/DaNSCOoudkZxWbYMDjgANm2Cujq4/35obS13qwxjeMwiNLyxdGlGBAcGMn+XLi13iwwjP0wIDW/MnJmxBKurM39nzoxXn7nZRqkw19jwRmtrxh1eujQjgnHcYnOzjVJiQmh4pbXVj2AN5WabEBrFwlxjQyU+3WxzsY3hGFYIRaReRFaIyEoR6RCR70Tli0TkGRF5QkQWikhtVH6UiDwePdpEZGpOXd+Kju8QkVNyyg+PytIismdOebOIbBSR9ujxPzmvTReR/xOR50TkMhERXx+KUX6ybvb8+fHc4qyLfc45mb8mhsZQ5GMRpoD9nXNTgRbgUyIyA1gE7AH8M9AAHB8d/wKwr3PuQ8B84GoAEZkCnADsBUwFPisiu0fnPAEcCjw4xPs/75xriR7/nlN+FXAisHv0+FR+XTYqhdZWmDcvnktsM9lGPgwrhC5DT/S0Nno459zi6DUHrAB2iY5vc86ti45fni0HPgAsd8696ZzrBx4AvhCd85Rz7pl8Gy0iE4FG59yy6P1/Bnw+3/MN/2h1P33PZBvJJK8YoYhUi0g78DJwn3Pu4ZzXaoHZwD1DnHoccHf0/xPAx0VkOxHZBvg08J483n5XEfmziDwgIh+LynYGVuccszoqM8pAR0ejWvfTl4udRavgG/HIa9bYOTcAtIjIWOB2EZninHsievlK4EHn3B9yzxGR/cgI4T5RHU+JyIXAfUAPsBLoH+at1wBNzrlXRWQ6cIeITAaGige6oSoQkRPJuNBMmDCBpXn6Rj09PXkfWyo0tglgxYoJpFKOdFpIpdIsXNhJKtVd7mZt9nm1tkIqFc817uhoZO7cqfT1VVFbm+aSS1YyefKGWO3SRNDtcs6N6AGcB5ya8/8dQNWgYz4EPA+8byv1LAC+MahsKbDnVs5ZCuwJTASezik/EvjxcG2fPn26y5clS5bkfWyp0Ngm55y74orHXEODc9XVzjU0ONfWFq++tjbnFiyIX4/vz2vBgkwfIfN3wQId7fJFktoFPOpGoGvDWoQiMh7oc86tF5EG4EDgQhE5HjgIOMA5l845vgm4DZjtnPvLoLp2cM69HB1zKLBVRyV679eccwMishuZSZG/OudeE5F/RJM2DwNfAS4fri9GcZg8eUMQidTZeGO2bbHjjc3N0N0NTU3Q2Rm7fUbh5OMaTwSuE5FqMjHFm51zvxaRfqALWBZlrtzmnDsfOBfYDrgyKu93zmVTYm4Vke2APuAkF02qiMgXyAjZeOA3ItLunDsI+DhwfvReA8C/O+dei+r6OnAtmRnru3k7FmmUgRASqX2unOlvaqJ61apMjKerKyOKJoZlY1ghdM49Dnx4iPIhz3XOHc/bqTSDX/vYFspvB24fovxW4NYtnPMoMGWLDTcqEt9W17JlfoQriy/BX71qFbuQ8wPsLn9MNWRsiZ2hCp9WV0dHI6edptPN/glwNhkXqwoy7nGB+Bb7EDEhNNThy+pqbx+r1s3+ngi/d457gHdPmlSwW6w5plpJ2FpjIzYdHY0qc+taWtarTaZ2zrEceLdzsWKDtnLGD2YRGrFYtgzmzp1Kf78+i8T3bLZG91N7TLVSMCE0YrF0KfT1VZFO63M/wY+bXQz3U0Sora2NVwnhxFSLjQmhEYuZM6G2Nk1/f7U699MXxUjpSafTwx+UJyHEVIuNCaERi9ZWuOSSlWzYMM2LO6XRNfPtftZ89rOsf+MNxsaYJCkG2Ziqt4TxCsKE0IjN5MkbvPxotM6A+nQ/aW5mrzfeoApwXV2IokRqnzFV0DmobQkTQkMN2leVeGlLdzfVQHXOc0346qfWQW1LWPqMoQbf2/MvWtSkLqWHpiY2kVkvmn0eB63bglVaWo8JoaEG39vzL1y4q7r9ES/8+tfZn2jdaMwYoebbEPgc1EqRp2qusVEw2RhQY2Ojt8C6D9csa42k06LOxf7BD37AS8DXt92WL8WMDWoPJfiIN5YqT9WE0CiI3BhQTc1Upk3T8yPMWiOpVJq6uipVScaf+9znuOaaazj99NPjVUQRtgXzjK9BrRR5qiaERkHkWiPOiUprZOHCTo49djdV+yNeffXVfPnLX2amB9XyOZudjamOGqXne4TS5amaEBoFkWuN1NQ4ldZIKtVNa+tuserx7X5uv/32TJkyxdvW8z5XzqRSu7Joka4ZXt95qlvChNAoiFxrpLFxJa2t02LVpzXnzKf7mXrPe/j3V1/lmQce8NU8L2iOqYK/PNWtYUJoFEzWGlm6dOQ3MMpFc86ZN/ezuZna1av5T6K7jClKpNYcUy0VJoRG2dE8+wmekoyjxOkaIiFUlEitOaZaKkwIjbLj0/3MtUZU0dTEQFcXA0Q/upiJ1L7RGlMtFSaERtnxmXOWa41cdJG//MbYdHZyqAh7Ad/ecUdGx3SLtbqfxRjUfOapbgkTQkMFPhOps9ZIe/tYL23zxdhZs/jZ8uV87Cc/YWaMejS7n8UY1EqRp2pL7IwRo3V96+BlXS0t62PV57ufP//5z3n22Wdj16N9HW9rK8ybF0+4cvvY1ydF76NZhMaIGMoa0cJgaySVKnw227fVdcghh/CrX/2KM888k09+8pOFV0TxtucvhQuaL6XOUzUhNEaE9mB4rosdx4rw3c/FixfjnOORRx6JLYS+V5RoXCrpO091OEwIjRExlDWSShVeXwhBf4D+/n4g4x4//fTTsdvna99A7UslfeSp5oMJoTEihrJGCrW8Qgj6ZzkI2BfYccYMnr722tjt84VPF1TroJYPJoTGiCmGNaLdzY5FczO3A3UAXV185MgjYc0aDxXHx5cLqnlQywebNTbKhs/NO0HvbPYrXV3U8vb2/A1r15azOe8gO8s7eXLhLqj2mezhMIvQKBvFCvprm80eP2kSrqsLBwiwcccdGV1gXVrdT+17Iw6HCaFRVoJwszs7kebmzPripiYeufbaghKqNbufvmOqpV4qaUJojIhKskjizGb7pKoqE4F666buBfqNqsWe4t0BrxRLJU0IjbypNIskzmy2T7F3zsWvhOIlUmsb1MqxVNKE0MibECwS32K/1uPESCgx1cGCH3epZD6YEBp5U6ztsjSJqW+xnz17NgA1NX5+aiHEVH0ulcwXE0Ijb4q1XZYmF9u3+zntd79jT2CSJ/fYF5pjquBvqWS+mBAaI6IY22VptkZitau5me9G/1YPDKjant9nTBX0Wvj5YkJolBztOWfeVpR0d2/+A1O0PT8Ub5ZXk4WfL7ayxCg5WWtk/nw/PxqtK0resR1/zO35tfaz0leVgFmERpkIwRq589JLufTzn+cs4MBJk2K5xZr7mYRJNBNCIy+0xoA0xxtPPfVUngNWjB5NT8zYoOZ+JmESzYTQGJZQrBHwK/idkfjN9BAEDSGuWk6xNyE0hiUEawT8C352M9brr7++8EoifPdTo3VfTrE3ITSGJQRrBIon+GPH+lkipnHljE98b9wwEkwIjWEp5s4iWn6E4Fnwm5u5GnjIT9O8odm6B4+pSyPEhNDIixBmeb0JfnMzdHVxAnBC9rmSRGrNMdVyYkJolJQgLJLu7rc2Yc0+14LmmGo5MSE0SkoScs6GpamJ/q4uIPMDk5iJ1L7RHlMtByaERklJQs7ZsHR2cpgIHwb+M2YiNegV/CS52SaERsmp9JyzfLgrevynBxHUKvhJcrNtrbGxVbSub9V8B7xOjxMj2tfxZu+AV4yBrZSYRWhskXKP0ltDszVy1FFHAVBXV1d4JRFBxFQpf66qCaGxRbS7n1qD/k1NTbS1tXHAAQfEblsQMVXKm0wNJoTGVkhSMHxr+O7njTfeyI033uijaUAYMVUoXzI1mBAaW0Gz++kT39bILQ0NfKq3lzEeZox9EcqgVigmhMZW0ep++sabNdLczGd6e6kD6OpSs6oklEGtUEwIjZIQikXyelcXo4HqbIGyVSUhDGqFYEJolIRQLJL/rarixHSaKqLctBirSrSKfRIHNRNCo2SEYJHMBX4B3Ag0x4gRahb7JA5qJoTGFtEwUg9FuXPOtkY6nWY5sOPGjVBfX3A9msUekjeomRAaQ6JlpB4Kzbs1f/CDH+T555+nPoYIQjLdz6HQMqiZEBpDomWk3hJad2vu6OiIV0FEEt3PoSh3InUWW2tsDInPtbxa1yv7Xt968cUXM2bMGC677DIfzUvMOt7h8NXPOJhFaAxJCEu7fLtl/zjtNE4GNv3Hf8CcOR5a6IdQ1ivHwYTQ2CJJX9rl1S1rbub/kdmVujqdVpNIDWEManExITSKipZg+JbwtqKku5sqdCZSQ/IHtbiYEBpFJZQ74NHUxKauLmqJflSWSF1RmBAaRSeEO+A9fc89HPuBD3Ai8FVLpM4LTYJvQmi8A00XaC6aXbMdd9yRhv335/6dduKrP/95wfVo7iMkd1AzITQ2Q9sFmovm2c+xY8dy//33x64nye5nLtoE34TQ2AxtF2gummc/x4wZQyqVoq+vL1Y9SXY/c9Em+CaExmZou0AHo3X284033ohXQQ5JdT9z0bKiJIsJobEZmtfx+sK72Dc385/AY7Fb5hfN1j2Ud2v+wZgQGu9A6zpeX/hOpHZdXZwNuOi5lkTqUDZu8IEJoVEUgrFGosTpGiCd81wDoWzc4AMTQqMoBGONNDUx0NWFI1pVEiORuhgkbd/AYmFCaBSFYKyRzk5qmptZ39XFmPe8h6qYbrFWwdecuuQDE0KjaARjjXR2MtZDNZoFX3Pqkg9sP0LjLbTuG+hzb0Tw28/e3l5EhMbGxth1hbBvoNY+mkVoAHpHatDtZs+ePRuAN998s/BKIpLufoLePFUTQgPQ735qdbPvvfdeAHbZZZfYbUu6+wn6EqmzmBAaQDizvL772dPTA8CCBQtitw30rpzxiaZE6iwmhAag2/30ie9k6suc4xbgy1/+sp8GeiCUQc0nJoTGW2h1P33jpZ/NzdDVxdeBE7LPlawoCWVQ84kJoeGdICySaAVJdfTQtKIEwhnUfGFCaHgnCIukqQm6ujZ/HgOVYo/eWV7fmBAaRSHxFklnJ6eKsB3wL+96Fx+P4RarFXvC2R9x2IRqEakXkRUislJEOkTkO1H5IhF5RkSeEJGFIlIblR8lIo9HjzYRmZpT17ei4ztE5JSc8sOjsrSI7Dno/eeJyHPRex2UU/6pqOw5ETnTx4cRMqEkU/vkEuAsoOnxx2PVozXJOIuPRGp4W/DPOSfzV9O1lo9FmAL2d871RGL3kIjcDSwCZkXH3AAcD1wFvADs65xbJyIHA1cDe4vIFDJx5b2ATcA9IvIb59yzwBPAocCPc99YRD4IHAFMBnYCfici74te/hHwCWA18IiI3OWce7KgTyFwQrFIimWNNDc3xzo/iJgqiq178hBC55wDeqKntdHDOecWZ48RkRXALtHxbTmnL8+WAx8Aljvn3ozOeQD4AvDfzrmnorLBb38IcJNzLgW8ICLPkRFSgOecc3+NzrspOtaEsAA0X6Bg+yOOBM391BxvzGutsYhUi0g78DJwn3Pu4ZzXaoHZwD1DnHoccHf0/xPAx0VkOxHZBvg08J5h3npnYFXO89VR2ZbKjQLw6X5qdbF9u58XXnghMOTgXRC+3E/NbnZW8OfP1yXQkOdkiXNuAGgRkbHA7SIyxTn3RPTylcCDzrk/5J4jIvuREcJ9ojqeEpELgfvIWJgrgf5h3nqoq8wxtIC7ISsQORE4EWDChAkszfPK6OnpyfvYUlHMNl10USPt7WNpaVlPKrVhRD+gbLs6OhqZO3cqfX1V1NamueSSlUyevKEo7R1JuwAaGxupqZmKc0JNjaOxcSVLlxbetg+eey6XkYnxjPQ7Keb3GKefg9vV0fH2NeHze2xthVQqf5EuyW/ROTeiB3AecGrO/3cAVYOO+RDwPPC+rdSzAPjGoLKlwJ45z+cB83Ke3wu0Ro97t3Tclh7Tp093+bJkyZK8jy0VGtvk3NvtWrDAuepq5yDzd8ECHe3K0taWaVNbW8yKJ03KdDL7mDQpVrt8U2g/c9vV1uZcQ0Pme2xo8PCZxaCQzwt41I1A14a1CEVkPNDnnFsvIg3AgcCFInI8cBBwgHMunXN8E3AbMNs595dBde3gnHs5OubQSNC2xl3ADSLyfTKTJbsDK8hYiruLyK7A38hMqOhZ4xQo2oP+Prfnd+S4KwlMptYeN/ZNPq7xROA6Eakm45Le7Jz7tYj0A13AsihOcptz7nzgXGA74MqovN85l02JuVVEtgP6gJOcc+sAROQLwOXAeOA3ItLunDvIOdchIjeTmQTpj84ZiM75JhkLsRpY6JzriP1pGLEIJehPUxOpri5qiH5ACUym1jyxUQzymTV+HPjwEOVDnuucO55MKs1Qr31sC+W3A7dv4bXvAt8donwxsPidZxjlJPGJ1MD69nY+M24cxwLHTZoUa42xVsEPJZE6i60sCRytF6hmi2TWrFm0AStqajgu5kYLmgXf16CmVexzMSEMGM0XqGaLJDuDGTeRGnQLvi80i30WE8KA0X6BarVIstvyX3zxxbHbVgkrZ+JSCWJvQhgwlXCB+sC34I8fP55XXnmFQw45xEv7bOVM+TEhDBjN7qdPvAp+czMvvfKKuhu5h2LdFwsTwsDR6n76xJvgNzfjsnsQdnUhinal1p7DqR0TQsMLQVgkUeK0AGlAFCVSB5PDWSRMCA0vBGGRNDUx0NWFI9qeX1kidQg5nMXChNDwQhAWSWcnR4nwXuCsnXdmdAITqSGQQW0QJoSGN0KwSG6O/n539epY9WjuYxCD2iBMCANF8yjt0yLR2k/tqUshDGq5mBAGiPZR2pdF0tHRyGmn+etnb28v48aNY4cddii8kghLXdKFCWGAVMIo7cMiaW8f67Wf9fX1vPbaa/EalYOlLunBhDBAfI/SHR2NLFum70JvaVnvtZ977LEHa9euZe3atdTX13tpow+0D2zak6nBhDBIfAfD586dSn+/Pmtk8uQNXq2RDzzzDP8G1O+xh5pEaggjplpsTAgDxWcwvK+vinQ64dZIczM3krmFI11doGhVidaYaiVhQmjEYuZMqK1N099fneics/5oR+rqbIGiVSWgM6ZaSZgQGrFobYVLLlnJhg3TEp1zdlt9PZ/t7aWK6BaKMVaVhBJT1TqoDYUJoRGbyZM3eEmL0Bz0P6aqig8B/wNMjbE9fygxVc2D2lDkdYN3I1lovQm77xvNL1rU5K2Pf/zjH+nfc09SDz8cKzaYjalqvAE7ZMQq6TeaHwqzCAND80jtK+if7WMqtSuLFvnpY0tLC4888ki8SggnplopidRZTAgDQ7P7CX6C/tk+ptPirY9VVVWMHz+el156KVY9ocRUKyWROosJYWAUI+essbFR1Yif7WMqlaaurip22+68806cc7zyyis+mhdETBUqI5E6iwlhYPh2PzdtgpqaqUybpueiz/Zx4cJOjj12t9jtWvXFL3ImIM55aZ8vQlkhVApMCAPEp/s5MADOiUprJJXqprV1t3gVNTfztYEBhGhmMYGJ1KB7NrsUmBAaBZFrjdTUuMQG/enupopMIrWLnmsilBVCxcaE0CiIXGuksXElra3TCq5Lc9Cfpib6o+35a6Lnccim9YwapaiP+J/NrjRMCI2CyVojS5duiFWP6qB/Zyejmpt5o6uLmhiJ1FCctB5f+J7NVmndbwUTwoDQeoGqzznr7GS0h2qKkdbjEx+z2aqt+61gK0sCIXuBnnNO5q+mVSVZN3v+fD8/nOzKmY6OxthtmzRpEiLC7NmzY9eVFfyqqrSX1CWNq4MqbUVJFrMIA0G1+0lxdmv2kdazatUqAKZMmRK7bb7SejRbXeqt+y1gQhgIxbxFoyZ8p/W4KHfwjDPO8NI+H2k9mge1SltRksWEMBCKeYvGiy7Ss7LEa1pPczMLgAf9NM0b2q2uSlpRksWEMCCKdYvG9vax8Sv1hLe0nuZmXFcXpwFzo+dJTKQGvUslS4kJoTFiBlskLS3rC66rGDPZXtJ6osTpGiCd81wLWmOqlYoJoTFiBlskqVRhgqM56E9TE9LVBcTfkRrCialWKiaERkHkWiSFpkhoDvrT2Zlxh7u7MyLoIZE68TFV9OaqDocJYSBovEC1B/2zFqGLGRsMIqaKcgt/GEwIA0DrBVqsoL8msQe/MdViEMRSyWEwIQwAzRdoMYL+WbGPw7XXXguAiMRum6+YapZKEXxtFv7WMCEMgGImU2v5IfoW+/POOw+AMWPGeGmfj5gq6LXuoXKTqcGEMAiKmUwd1/LyxVBin0oVXt9nurt5N/DRgQFPLfSDZuseKjOZGkwIg6FYydRafohDiX3BlldzM5dBZlfqN99UlUxdjHvOVJr1VgxMCI0R4dvy8ok3a6S7e/MfhqJk6mLcc0aTZV8uTAiNEeHT8tJqkfTvsgvVq1bx1jRJjGTqYq6ciYNWy75cmBAaI8bHD1Fz0P/4/ffnmeuu42zgszF2pdbcR9+WvdZBLV9MCI2yoNkiueOOO3gd+NpOO/G3GLFBzX30bdlrFfx8MSFMOFpHas0pPRs2ZPL84u5BqD2vLukTaCPBhDDBaB6pi53S42Mz1jlz5hReCeGunNEm+PlgQphgtI/UIVgkxVw5o6mPlZpIncWEMMFodj994rWfzc3cCSz30zRvaBZ7qNxE6iwmhAlGs/vpE2/9bG6Gri4+B3wu+zyBidSgd1ArFyaECScE9xM89bO7G0dmRYkDJIGJ1KB7UCsXJoRGXgSxtKupib6uLoTohxFzV2rfhDKolQMTQiMvirm0S82PsLOTz4vwL8DZMRKps2gV/CAGtRFiQmjkTQhLu+6OHmd7EEGtgh/EoDZCqsrdAKN4LFsG3/te5q8WstZIdbW/oL+vPl522WWAn81YhxJ8TbS2wrx58YRLex9HglmECUXraK056H/hhRcCmfv7xiWEWd4kJFJnMSFMKJpdUK1B/zPOOIOzzz6b888/P3bbNAu+L5KQSJ3FhDChhBAQ922RzJkzJ/ayuly0Cr5PKj2ROosJYUIJISDu2yI5VYTZwFQPM8Y+CWFQKzcmhAkmhFleXxZJf1MT5wOjALq6VK0qCWFQKzcmhMZWCSHoD7B61Sp2AaqzBYpWlUAYg1o5MSE0tkoIQX+AnwJnkcknqwJ12/P7IJRBrRBMCI1hCSHov0CE+53jbmBsQrfnD2VQKwQTwgSidaTWnHfmnGM5MDbalLVQNIs9hDGoFYIJYcLQPFJr3rF58eLFXH755fEqIRz3U/OgVggmhAlD+0itdcfmgw8+mIMPPjh2u0JxP5OUTA221jhxaF7L6xOf61w//elPIyIceeSRXtrmYx0v6F/L66ufGjCLMGGEYpH4dM3ef/fdTAUm33abp9b5IWnup2ZMCBNICAFxn9vzX0jGNaretCmRidSgN9aoBRNCY4tot0h8bc9fTbITqTVb9lqwGKGxRbIWyfz58X88WmONNDWxCRjIeR4Hjf3UHmvUgFmExlZJukVyy8UXc+nhh3My8KWYmy1o7WcoKT1xMCFMGBovUs2xxsMOO4zPbtzI2rVrM/HBGGjtZygTaHEwIUwQWi9S7RZJfX09zTFFEHTHVEOYQIuDCWGC0HqRarZIqqqqcM7hYi6tgzBmeTWLfRxMCBOE5g08tVokPgQwl6THVJO2oiSLCWGCCGEDT59iv379ej4KHAuq8ge1WvZZkrI9fy4mhAkj6Rt4+rRIlk6cyH1APajalVp7TDWJmBAa70B7HMiXRfKu3l7qyEmmVZJMrTmmmlRMCI13oHm7LJ/8DvgomVUlArGTqX2iNaaaVEwIjSHRul2WT77nXMYd7u7OiGDMZGqNYm9udn6YECYErReoeovEQ0xQs9ibm50fJoQJQPMFqtUiuf/++znwwAOpqamhr68vVpu0i7252cMz7KYLIlIvIitEZKWIdIjId6LyRSLyjIg8ISILRaQ2Kj9KRB6PHm0iMjWnrm9Fx3eIyCk55duKyH0i8mz0d1xUPlNEXheR9uhxbs45n4re/zkROdPnh1JpaF5U73vjhgMOgHPOyfyNs7HBySefDEBtbW3hlUSEshmu735qIh+LMAXs75zricTuIRG5G1gEzIqOuQE4HrgKeAHY1zm3TkQOBq4G9haRKcAJwF7AJuAeEfmNc+5Z4EzgfufcBZGonQmcEdX9B+fcZ3MbJCLVwI+ATwCrgUdE5C7n3JMFfg4VTSizvD4tkmeffRaAf/mXf4ndrlDcz6QmU0MeQugyqfc90dPa6OGcc4uzx4jICmCX6Pi2nNOXZ8uBDwDLnXNvRuc8AHwB+G/gEGBmdNx1wFLeFsKh2At4zjn316ium6I6ghTCUGZ5fQr+of39vBc466mnvLRNo9gXgyQmU0OeMcLIAnsMeC/wI+fcwzmv1QKzgW8NcepxwN3R/08A3xXvN0e6AAAgAElEQVSR7YCNwKeBR6PXJjjn1gA459aIyA45dbSKyErgReBU51wHsDOwKueY1cDe+fQlqYQwy+tzV+rryFz8NX/7m5pEatAbU006MpK1liIyFrgdONk590RU9hPgDefcKYOO3Q+4EtjHOfdqVHYccBIZC/NJYKNz7j9EZL1zbmzOueucc+NEpBFIR275p4FLnXO7i8jhwEHOueOj42cDeznnTh6izScCJwJMmDBh+k033ZRXX3t6ehgzZkzen00pKEWbFi1qYuHCXUmnhaqqNMce28lRR2090VjjZwVbbtfH99sPyATIXfR4cMmSsrcrS0dHI+3tY2lpWc/kyRsKfp+Ojkbmzp1KX18VtbVpLrlk5Vbrq7TvcWvst99+jznn9sz7hOzOG/k+gPPIWGbZ/+8AqgYd8yHgeeB9W6lnAfCN6P9ngInR/xOBZ7ZwTiewPdAK3JtTPg+YN1zbp0+f7vJlyZIleR9bKkrRprY25xoanKuuzvxtayu8XW1tzi1YkF8dxWCLn9ekSc7B249Jk0rYqtJdWwsWZL5HyPxdsEBHu0ZKIe0CHnUj0LV8Zo3HR5YgItIAHAg8LSLHAwcBRzrn0jnHNwG3AbOdc38ZVNcOOcccCtwYvXQXcHT0/9HAndFxO4qIRP/vRWYQfxV4BNhdRHYVkTrgiKgOIya+Znl9zvD65oYFC+gePx5EIOau1GCzvEkgnxjhROC6KE5YBdzsnPu1iPQDXcCySKtuc86dD5wLbAdcGZX3u7dN1FujGGEfcJJzbl1UfgFwc+Q6dwOHR+WHAV+P3msjcESk9v0i8k3gXjIrpBa6TOzQ8EDSN26YNWsWzjkuuOACzjhja3NywxNETJXkxxrzmTV+HPjwEOVDnusycbvjt/Dax7ZQ/ipwwBDlVwBXbOGcxcDioV4zyo/moL+L4uJxRRB0Cz4kf39EX9jKkgpH60gdSm6d5s1wfaFd7H1gQljBaBYI0Jlbt3bt2vgNysE2w00GJoQVTAgjNfj9Ic6ePRuAmhp/l37SY6pJXlGSxYSwgtEch/OJzx/iRz7yEXb93e84s7/fEqlHQFJXlGQxIaxgQonDgb8f4oIbbnj7iaLt+UP6LjViQljhaIzDaaa9q4upRDtSg5rt+cG+y3IybEK1EQahbCV1Cpmtj94ixvb8WvsYynfpE7MIDSAM1+zOO+/kATJ7wf1MJNb2/Fr7CGF8l74xIaxgtAbEtbpmp556KgC3brMNP3vjjVh1ae1jlqR/l74xIaxQNI/UWvPOOiPrb79o95k4aJ/l9YXW79I3JoQViuaRWusa1/7+fgCuv/76eBURjvsZQg4hmBBWLNqXdmlc43rYYYfx+9//nrFjxw5/cB6E4n4mPYcQTAgrlhCWdvkWiF/+8pe+muYV7YNaCJgQVjBJX9rlUyAuvvhinj/tNM4GdvGwB6FPQhjUtGNCGDiag+E+41MbTz+dS4B6ULWiJEvSBzXtmBAGjvY74PmKT73fOerIWUGgaEWJLzQPatoxITSCuAPeDcC/knPBx1hRAjpjcdoHNc2YEBre0Oya3UnmJtr3xFxRAroFP4RBrRjYWuMKROvaT61rXLObsd4LkE7Hjg0OJfhJI4Q+5mIWYYWhaaR+/cVunlx8M/3v+WdAb5LxjjvuyNVXX82aNWsKb1AOIaS7hBZvNCGsMDS5n3+4agHP/GU1NQP30/0/FzGq2lHdUMcHJ+zADn1HAp8quG7f/TzhhBMKP3kQIaS7hLKiJIsJYYWhaY3rzlP2pOeFLgYGNtGL8Hfq6N9UDate5bkfXUHDD3/AaNdPbZ1QPaaBse99HzOO/hbv3mn4iQqf/fy2CNOAWR7zB0NIdwlhRUkWE8IKQ5P7+eF/O54P/9vxLF26lJkzZ7Kpt5dHf3Ypf/vzI/Sv30A//bxJDa+6WlyPY3X7Mzz5p68zun8TDVUD1IyqpXbcu5nUui9TDzuOuvp6//1sbuZ7ZG5+rS1/MDT3UzMmhBWI1jWudfX1fPTEd94r+JW/PsOjP7+CDZ2dDLzRyyYcr1NHaqAG/v4mXb+6m2W3/5ox6b6Me73NKLbZeSdavvhV5s3bp/AGAXR3U03Oha4of9DSXfRgQhgwpXKzx+/2fg4+7/J3HP/0vbfx1D230vvKqwz099GL8AqjGEhVwV/X8uxFF7BNXx/b0E9dbRVVjaPZ/v0fZMZx32b0uO3zatMz227Le159FSGyCmPmD/rG0l10YEIYMOV2s/c46FD2OOjQzcpS/9jA8ut+yEuP/5mBDT30keZNV83fXS1s6Gf1I4/z+MNfYUx/Hw1VaWrqa6jbdlt22/cTtBx27Dve46AxY5j46qt8F9jfQ4xQq9WlPd6oHRPCwNHmZo96VyP7fvPcd5Svfaqdx274MT2rVjHwRooUsJ5aUv018PI/eOGXt/GHG2/OuNc1ULXNKKq23ZZNr69jOXDVYYexf8zdZzRbXSGk9BQTE8IKoqOjkWXLdF6gxQ787/iBFj4z/6p3lK+87Vqef+C3pP7+KgP9/fRSxUvUku6tghfXM/cT+zC6r4+G/le58YhPUN04hgmTP8SMY+cy6l2NI2qDZqsrhJSeYmJCWCEsWwZz506lv1/nBVquXamnHvpVph761c3K3lj3d1Zcdykv/d9KXM9G+iTNG9Tw93QtvN7HqrbH+PNDRzCmv4/6qjTV9bWMGr8tux/wr/zzvx65xffSbnWFkNJTLEwIK4SlS6Gvr4p0Wu8FqmVX6tHjtme/U+Zzzz33cPDBB1NbW8umTZtY/efl/Pnmhbzx4moG3kyxCVhHHZv6q2HNBl64fhEPXHsdY9J91NVC1egGGpua2POob7DD+6Z4E/uOjkZOO02n1RVqSo8JYYUwcybU1qbp768ueyJ1MfFpkVx66aUADAwMALDLh2ewy4dnbHbMpt5eOu68nuf/cB99r61noL+fjVSxztWSfhP+9nQ3T/+/M9imr4/RMkBNXRX/NPZd9LVPIzVlzojda4D29rFqra7QVpRkMSGsEFpb4ZJLVrJhw7SyJ1IXE58WyUMPPcQM4Ifp9BYTqevq699KDM+l5+U1PHzdpbz69JMM9LzJJnH8gxo2ulpYt4lVS5fzp/vbGDOwiYZql3GvJ4zngwd/kfcf+PmttqulZb2a1UFDEdKKkiwmhBXE5MkbvLgqmuNAPi2Sw3p6+BHQACNeVTJmh4kccNoF7yh/YfkS/u/2n/HmmrUMvLmJFMKr1NHXVw2r1/HXn/yU3115FaNdP3W1maWFjbvuxoyjv8W4pt2AzPeoZXWQkcGEMEC0x4F8WSR7AnWAZAs8rCrZdcZ+7Dpj8/sib+rtpf3mn9C9/EH61r9Of98AG6niNVeHe8PBE8/z1KknM7pvE9tImupRVdSNvYKZ02cw/cPfJLqBQEFoHtQqCRPCAPEdB9Ka1vMQsNmeMzFWlWzN/ayrr2evr5zMXl85ebPy11/sZvn//pD1zz/LwBu99BG51+laeK2XrvuW0vm7eznypvsKbpemTTgqGRPCQPG5tEtrWs/Xliyh7qtfzViCMXalLtT9fPdOTRx09vffUf7rn/yQgWf+zMaXXmZM8z8V1KYsPgc1zbPZxcaEsEJYtgwWLWpi1ChdF6fWtJ7e3t7MPx52mvHtfo7ZvYWZJ5wSu11ZfA1qmmezi41t1V8BZC2ShQt35YADdG3Rn03r8bE9v89bEIwdO5b99tuPww8/PHZdvm9BoJXsbHbS+zkUZhFWAFmLJJ0WdSO1r7Qe37OfqVQKgLPPPrvwSiJCian6nM2uNEwIK4CsRZJKpamrq1IXEPeR1lOs2c+Wlpb4lRBGTBXCzCEEc40rgqxFcuyxnbF/OFnL65xzUOVm+3Y/Pw5cAZncQUVkY6qh3B2uUjCLsEJobYVUqpvW1t1i1aM178yn+/n4tttyN1F2nsLt+UNYKllpmBAGhuZkal9u2Q7r1jGKHHdH2fb8ISyVrDRMCAOjGNtlNTY2qhLUHSdNgq4uHNGqkpjb82uMqYJe674SMSEMEN/bZdXUTGXaND0/wjmf+xzn3nIL265Zg8Tcnl+z1eV7f0SNeaqlwiZLlOMzt84nudZIX5+oCvpffvnljF+zhgeXLIkdGxzK6tJC1rqfPz+eQGvOUy0VZhEqZrA1ctFFelzQXGukpsapCfq/taLEE5pjquB3V2qNeaqlwoRQMYOtkfb2seVu0lvkxhobG1fS2jqt4Lp8up/HHHMMAFVVfpydEGKqvvNUKxETQsUMtkZaWtaXu0mbkbVGli7dEKsen0H/3/zmNwDssssusdqUS9JjqlmxX7iwk2OP3U1Nu0qJCaFiBlsjqVQ8wdGac+bT/ezp6WEGcHd3N7VHHglr1nhqZTxyxd45Ued++spTrVRMCJWTa43ECdRrnv306X7+T1UVxwwMUAu4tWvVJFNrjakaGUwIA0F7zpmvZOoT0+nNC5QkU2uNqRoZLH0mEHyv5dWY1tPe3s53AJdbGDOZ2ietrTBvXiahOg6aU3oqFbMIA8H37KfGtJ5jjjmGduf4B3CxCG9OmMDomMnUGt3PYm3Pr202u5SYECpG620atab1PPnkkwAs++hH4Y9/5JGlS5lZYF2a3c9iDWraZrNLiQmhUob6IWrBZ1qPT7HftGkTANdcc028iggnpqp9NrtUmBAqRfMP0VdaT7Gsrj322CN2HcW8O5wmfM9mVyomhEoZ6ocY7T6vAh9pPZUk9kmMqYLf2exKxoRQKUP9EAsVnFCC/i0tLWyzzTY+mgYkP6aaxdcKoUrGhFAxvpd2JTnoD/Dndetg5Uo1SdRZtMZUjbcxIUw4mt1P8Gd1XVxdzSnpNNWAZLfnv/ba+BV7QHtM1bCE6sQTQiI1wLuiFSWSLVCyoiRLNpk6jnBZInXxMIsw4RQz6K/JIvklcDQ5F7Sy7fl9oH1vxErGhDAAihX01+Jmr1+/nvuB2cAvRTIi2NlZsMmkVfB9x1Q1in25MCFUiNYLtJi5dXH6efTRRwNwR00N9PXFaxR6BR/83mhea8J+OTAhVIZWawT0utlLI8tv0qRJhTcoB983RdI4qGkW+3JgQqgM7ReoRjf79ddf5/zzz+dLX/pS/IbhT/A1D2raE/ZLjQmhMrS6n77x3c9zzz3XR7PewudNkTQOaj4T9pOACaEytLqfvvHZzy+JcALwiZj3MPaN9lleX9Z9EjAhVIhG97MY+OjnqgkTuBYYBZBNpFYihjbLWzmYECaYEIL+615+mYlAdbZAYSJ1sWZ5NX0PlY4JYYIJIej/U+ACMkIokMhEatBv3Vc6JoQJJ+lB/8uBx4AHgJqYMULNgh+CdV9OTAgVofUC1R70bwNqnBv2uOHQLPghWPflxIRQCZovUM1Bf+ccvb298SqJ0J66lHTrvpyYECpB+wWqOehfX18fv2GEkbqk3bovF7YNlxJC2S7L51ZSO++8MyLCTI+/Zh/bZYHeLbOyYj9/vh5x1oBZhEoIwRoBvxbJmjVrmAHc9MADqvIHQffkhiVSvxMTQkWEkEjtU/A/7xyLgAZIbDK15kEtSZgQJhDtcSBfgv8RoC63IIHJ1JoHtSRhMcIE4jMOpDXWuHbtWn4NDOQWxkim1tpP37FjY2jMIkwoSb8D3s0330wbcC5wQe6u1AWguZ+aU5eShAmhsUU0u2Vz5sxhzpw5XurS3E/QnbqUFMw1VoJG10xzSs/atWvjVxIRivupNaVHA2YRKkDrSK05pWfixIlAZmVJXEJxP7VPopUTE0IFaHbNQkjpgTDcT9+CnySGdY1FpF5EVojIShHpEJHvROWLROQZEXlCRBaKSG1UfpSIPB492kRkak5d34qO7xCRU3LKtxWR+0Tk2ejvuKhcROQyEXkuqm9azjlHR8c/KyJH+/xQSo1mF9QXPvt45513MgO4HjK5g4rQ7n76WjmTNPKxCFPA/s65nkjsHhKRu4FFwKzomBuA44GrgBeAfZ1z60TkYOBqYG8RmQKcAOwFbALuEZHfOOeeBc4E7nfOXSAiZ0bPzwAOBnaPHntH9e8tItsC5wF7Ag54TETucs6ti/uBlAPNLqgvfPZx1Re/yP1APahLpNa8osTYMsMKocsEYXqip7XRwznnFmePEZEVwC7R8W05py/PlgMfAJY7596MznkA+ALw38AhwMzouOuApWSE8BDgZ1EblovIWBGZGB17n3Putaiu+4BPATfm33VdhOCC+upj08AAdeS4M4oSqW1FSWWS16yxiFSLSDvwMhkBejjntVpgNnDPEKceB9wd/f8E8HER2U5EtgE+Dbwnem2Cc24NQPR3h6h8Z2BVTn2ro7ItlQdPCDOgz1ZVkSbjCgCxd6X2jQ/3U7uLnTTymixxzg0ALSIyFrhdRKY4556IXr4SeNA594fcc0RkPzJCuE9Ux1MiciFwHxkLcyXQP8xby1DN2Ur5OysQORE4EWDChAkszfOK6unpyfvYUpFvmy66qJH29rG0tKwnldpQ8I+oo+PteiZP3hC7Xb6Yfv/99B15JDVr1/LmjjvyyLXXDqkU+bQr3z76JJ92NTY2UlMzFeeEmhpHY+NKli4tvH359FPjNQ8lapdzbkQPMrG5U3P+vwOoGnTMh4DngfdtpZ4FwDei/58BJkb/TwSeif7/MXBkzjnPRK8fCfw4p3yz47b0mD59usuXJUuW5H1sqShlm9ranGtocK66OvO3rU1HuzZu3OgOPPBAt2bNmmGPHa5dI+mjT/L9vNranFuwIH678u2nxmveucLaBTzqRqBr+cwaj48sQUSkATgQeFpEjgcOigQonXN8E3AbMNs595dBde2Qc8yhvB3TuwvIzvweDdyZU/6VaPZ4BvC6y7jO9wKfFJFx0QzzJ6OyikPjDC/odc1OP/10fve737HTTjvFrktrH7MkfW9ETeTjGk8ErhORajIxxZudc78WkX6gC1gmIgC3OefOJ7P8czvgyqi83zm3Z1TXrSKyHdAHnOTenuW9ALhZRI4DuoHDo/LFZGKJzwFvAscAOOdeE5H5wCPRcee7aOKkktAcENe6bf2iRYsAGD9+fLwGobePvrFE6uHJZ9b4ceDDQ5QPea5z7ngyqTRDvfaxLZS/ChwwRLkDTtrCOQuBhVtseAWgfYZXY0rPunWZsfOkk4a8LEaE1j76xhKph8dWlpQR7RaJxpQeFy2pO/fcc+M3DJ19LAa2K/XWMSEsI6FYJD4F/wAyWfaakqhB/6BmbB0TwjITgkXiTfCbm7mXKHdK2YqSUAa1pGJCmBC0WyReBL+7m+pBzzURwqCWVEwIE0IIFsnNo0bxxd5eqoiswpgrSrS6n7ZeufSYECaIpFsk/9bbywxgCVA/aVIst1ir2IOtVy4HJoRlROtordkiWQ7Ue9iMVavYZ7E74JUWE8IyoXm0DsEiCSHJOIQ++sKEsExoH621WSQ33HADANFqpdiEsD2/JVLnjwlhmdDsfvrCZx/POussAEaPHu2lbRDO9vxa2qIZE8IyEYL76dMiuf7661k5cyZf7ulRlT8I+q17Y3hMCMuINvezGPiySPaZNYt9BgYyT5QlU4dg3ScdE8IKJ5SA+O1dXRyCbc9vFAcTwgonhKB/b28v3wc+A9RlC2MkUxejjyFY90nGhDABJD3of9JJJ/EQcDLwY5GMCBboFmvtI+hfJplkTAiNt9Bqkdx6660A/GriRH784oux6tLaRwhjmaRWTAjLgNaRWqtFsmFD5mZDZ555ZrwGobePWZK+TFIrJoQlRvNIrdUiyW7GOmfOnMIbFKG1j74JZRLNFyaEJUb7SK3RItm4cSNXXHFF/EZFaOyjb2xVycgwISwxoeSc+exn/R57cGp3N1xxhZrcQdBvddmqkvwxISwxoeSc+ernLQ0N/GtvL3WAKEukDiF1KRRMCMtAKDlnPvo5preXaqKNWEFVIjUkP3UpFIa9wbuhk6xbVl3tb/ZT443mrwM2AW/tQBhzV2qtDDWwGaXDLMIKJZTZz5uAl4Dfx0ykzqLV/QwldqwVE8IKJoTZT8hszU86HbsezYIfSuxYK+YalxCt7qdvN9sXK1asAPxtxqrd/WxthXnz4gmX9j5qxSzCEqF5pNY6+7nXXnu9lUztgxDcT+0pPVoxISwR2t3PEGY/Q3A/LZG6MEwIS4T2Na6+8Cn4WZfYp1UYQuqSJVKPHBPCEhHKLK9Pwf9E9NCURA3hDGohYUJYQkKY5fUl+P1NTdwBjAJ1W/OHMqiFhAlhBaI9IO5D8FevWsUuQHW2IKErSjQPaiFh6TMVSNYimT/fjwWhMa3nGjIrSt7KHoy5okRjH0Fv6lJomEVYoSR9lvfqHXbgkZdf5iZg7KRJsdxirX0EvalLoWFCGDhaXbOXXnrJW11a+5gl6YNaJWCucYkIwTXz2cfsqhIfhOJ+2qqSwjGLsARoHqk1JhmvWLGCvffeGxEh7WGNcSjup/ZJNM2YEJaAEFwzn3382te+BkBdXV28RuUQgvtpq0oKx1zjEqDV/fSJzz4++eSTzAD+nEpl8gcVod399LFxQ4iYRVgCfI3UHR2NnHZa8q2RQzdt4hqgAdQlU4ewcUOImBCWCB+uWXv72MS72AAfgcw9SrIFipKpNcZUjfiYEFYQLS3rgwiGv0DO1vzgJZnap+WlLaZqxMeEsIKYPHlDELOflzuXcYe7u2Nvz6/V8rKNG3RhQlhhJH3288ILL+TRRx/l508/TX19fez6tFpetnGDLkwIi4zWkVqrQJx11lmk02nGjx/PlVdeGbs+zbl1tnGDHkwIi4jmkVqra5ZNoP7+978fr0ERvi0vjYOaZrGvFEwIi4jmkdqnQBQjrceHW5zFh+WleVCzROr4mBAWEd8jdUdHI8uW6Zr9BH9pPevXr4/fmCKheVAD254/LiaERcS3WzZ37lT6+/VZJL7Seo4++mhmAPNAVRI16A0lGH4wISwyPgPifX1VpNP6LBJfaT0X/vGPvJdoV2plK0q0hxKMeJgQVggzZ0JtbZr+/mqVAXEfgr/Ha69tXqBoRQnoCyUY/jAhrBBaW+GSS1ayYcO0xM5+/hcZt/it+5TEXFHiO6bqC58rhLR+l5WGCWEFMXnyBi+xKa2zn+c4xzLgNyJeVpRojan6CiVo/i4rDduGqwho3SoL9G4jdf/99wOwGCCdjh0bzMZUtfUzi4/tsrR+l5WIWYSe0T5KF2P2c9GiJkaNitfPOXPmANDQ0BCvQRHaY6o+sERqf5gQeqYS8s18r3FNpXZl0aJ4ov/ss89G7fPzYYUQU7VEan+YEHqmEkZp32tc02mJLfp9fX0ALFq0KH7DIpIeUwVLpPaFxQg9kx2lk3zz9SxZ0a+qSscW/auvvpofV1ez4047qdqe3+JwYWAWYREo1lZZF13UqMrCzIr+woWdHHvsbrH6fMJ3v5tRG1CVTK01pmr4xYRQMYOtkfb2seVu0jtobYVUqpvW1t0KruPiiy9mQlcXs9C3Pb/WmKrhFxNCxQy2Rlpa4m1KoDXof9FFF7Eb8GX8JFNn+9nY6MeC1hhTNfxiQqiYwdZIKrWh4Lo0B/1feeUVXgaWvetd7NPTEyuZOrefNTVTmTZNTz+zA1sqlaaurso2blCECaFninmjoDiBes1pPc5lbtW0z4bChT5Lbj+dE1X99BlT1TywVSImhB7RfHGGcj/e3H7W1DhVk0vgJ6YKuge2SsSE0COaL85i3o83Dr29vfEqGERuPxsbV9LaOi1WfVpFvxLyVSsJE0KPaL84Nd6Pt76+nnXr1rHMY6Jktp9Ll8ZztTVb+LaqxC8mhB7xfXFqtEaGEvtUKkaFzc2M7e7m4Ji7zRQDzRY+2KoSn5gQeqaY9x3WwFBiX/AkTnMz/V1dCFDV1YUoSaLOUqy4qqEPE0KlaLZGvFki3d0IObmDSpKosxQrrqpthZBhQqgW3y6oRjebpiY2dXVRS3QhxtyRuhgUI66qcYVQ6NimC0rxuXlD1iI555zMXy0bONx56aXsD9wEMGlSbLdY6wYV2UGtutrfCiGN/axkzCJUjO+lXdrc7NNPP52/AI9vsw2zPIigxpgqhLNCqJIxi9AjWkfqwRZJ3KC/rz5mb+i+7777xq5L+3ZZPrbmB/39rFTMIvSE5pG6mMnUcfr40ksvFX7yIHzGVFXGUyO056pWKiaEntDqfmbRmEztE19pPZoHNLBE6mJhQuiJENby+uxjb28vDQ0NiAjpdNpL+5Iu9lkskdo/JoSe0Op++sSnNXL66aczA5jrnJrdqKE4O1JrHNSMzTEh9EgIFokva2T8FVdwP1APqrbmL8aO1BoHNWNzTAiVEYpF8n7nqCMnbUHRqpKkpy0Z78SEUBmhWCQ3AP9KzgUYc1WJRsG3Gd7KwYRQISFYJHeJcKxz3CgSa2t+0Cv4IexGlBRMCD2g9QLVbJH4mikG3YJfzN2ItPQxCZgQxkTzBarVIlmxYgWNjY3sscce8RoUEULqkmaxTwImhDHRfoFqtEj23Xdfent7+djHPsaDDz4Yu20hpC5ptu6TgAlhTEKZ5fUp+Nn7lCxYsMBb+5KeumQrSoqLCWFMQpnl9Sn4nwZagX1mzVKRO5hF+6BmK0qKx7BCKCL1wIPAqOj4W5xz54nIImBPoA9YAXzNOdcnIkcBZ0Sn9wBfd86tjOr6D+B4wAH/BxzjnOsVkf2Bi4E64DHgOOdcv4jMBO4EXojqu805d35U16eAS8lscPxT59wF8T6Kwglhlteb4Dc3cwuZL1pTIjWEM6gZ7yQfizAF7O+c6xGRWuAhEbkbWATMio65gYzAXUVGtPZ1zq0TkYOBq4G9RWRnYA7wQefcRhG5GThCRH4GXAcc4Jz7i4icDxwNXBPV/Qfn3GdzGyQi1cCPgE8Aq4FHROQu59yThX4QGgjBInm9q4vR6N6eP+mDmvFOhhVC55wjY9kB1EYP55xbnD1GRFYAu0THt+Wcvjxbnov6pPMAABPCSURBVPN+DSLSB2wDvAhsB6Scc3+JjrkPmMfbQjgUewHPOef+Gr3/TcAhQEULYQgWSde4cXxg3TocIJDIRGqwyY1KI68YYWSBPQa8F/iRc+7hnNdqgdnAt4Y49TjgbgDn3N9E5GKgG9gI/NY591sREaBWRPZ0zj0KHAa8J6eOVhFZSUY0T3XOdQA7A6tyjlkN7J1PX7STdIvkQ6+9lnGHu7sTm0gN/gc1jWKfJPISQufcANAiImOB20VkinPuiejlK4EHnXN/yD1HRPYjI4T7RM/HkbHadgXWA78UkVnOuetF5AjgByIyCvgt0B9V8ydgUuSWfxq4A9idyJgY3Myh2i4iJwInAkyYMIGleW5S19PTM+yxHR2NtLePpaVlPZMnx7uZuK82ZWlsbKSmZirOCTU1jsbGlcS94bmPdp1yyin80/TpnHzttZmCGFssL1rURCq1K+m0kEqlWbiwk1TqbVd7JO0qFq2tmQ1ic5sxknZ1dDQyd+5U+vqqqK1Nc8klK4t2rWn4vIaiJO1yzo3oAZxHxjLL/n8HUDXomA8BzwPvyyk7HLgm5/lXgCuHqP+TwM1beO9OYHsyk4735pTPA+YN1/bp06e7fFmyZMlWX29rc66hwbnq6szftra8qy6Y4do0mLY25xYsiN+24eoZSbvIDFjxGpTTrq19B+X6vIZjJO1asCDTP8j8XbBAR7tKSSHtAh51I9C1fGaNxwN9zrn1ItIAHAhcKCLHAweRmeRI5xzfBNwGzHZvx/0g4xLPEJFtyLjGBwCPRufs4Jx7ObIIzwC+G5XvCLzknHMisheZzUpeJWNR7i4iuwJ/A44AvjxcX3yi1fXMxYeb7dP9zOYP+iKEmKrFGktDPq7xROC6KE5YRcZa+7WI9ANdwLJMmO+t1JZzyUyAXBmV9zvn9nTOPSwit5Bxd/uBP5OZUQY4TUQ+G9V/lXPu91H5YcDXo/faCBwRqX2/iHwTuJfMBORCl4kdlgztM7y+8Cn4J510EgBVVf7uGZb0mKolUpeGfGaNHwc+PET5kOc6544nk0oz1GvnkXGnB5efBpw2RPkVwBVbqGsxsHio10pBCNYI+BX8W2+9lRnAD9NpVfmDoNvyskTq4mMrS2KQdGsE/Ar+oa+/zhVAAyQ+mdosuMrChFAB2ndP8SX455DZmv+tKf8EJlNrtu6NLWNCqIAQdk8B2HXSpIwlmCVGMrVWq0uzdW9sGRNCJSR995Snn36af/7b37hrm204eOPGWMnUmgU/lEm0pGFCmCA0B/yPPvpo+vv7+XxVFamYu1NrFvxQJtGShgnhCNE8QmvdkRpg5cqVAEyZMiVeRegWfAhjEi1pmBCOgEoYoTXuSA2QSqUAuOaare2lkR+aBd8n2gU/SZgQjoCQRuhi9bWlpSV+JegVfJ9YMnXpMCEcASEFwn339ePAl0BV7iDoH9wsmbo0mBCOgJAC4V6tkeZmlmb/V5ZIHdLgZmwZE8IRElIg3Js10t29+b5pihKpQxrcjC3jb/W7MSKylkh1dfID4f9F5sY2b20YGTOR+nvfy/z1RWsrzJtXnDxOozIwi7BMhLS29RznuIfMHcBk0qREJlKD/qWSxpYxISwjIa1t/SNQ5YbcRDxvtIcTQlkqmUTMNR4BxXDL4qLdHfO5GWslhBN8uNnav9MkYhZhnmgdpbXPeh5zzDGAn81YLZHaKBYmhHmi1S3TPuv5wAMPALDzzjvHqyjCEqmNYmBCmCeaLS/NKT0vvvginZ2d1NfXx2+gR7QObFkskbq0mBDmiXbLywfeXbLo/sXNMe9fXAw0D2xG6TEhHAGaLS8f+F5NQlcXA0C6q4taRatJIIyBzcgfE8IyoDnfzOdqEkfmFoPZ59pI+sBm5I8JYRkIIt+sqQnX1UU/0T1KErgtfxbNA5uRHyaEeaD1hkiqLZHOTuaKUA98oK6OryR0NQkEMrAlHBPCYdB8cWoP+P8w+vtUtDt1IagW+xwSP7AlHBPCYdB8cWoP+IsIzjn22GOPgusIKbk4pL5qw4RwGLTHfzQH/NMxb9IE4awmAUukLicmhMMQSvxHszUSwmqSLJZIXR5MCPMghPiPb2vk4upq3kynmVBdzdf6+300MTbavwOjfJgQlgjtExvg0RppbuZb6TQCVA0MqNmavxK+A6M8mBCWCO0TG17p7qaKTDK1i55rIKjvwBgRJoQlRPPEhleamujr6sIRXWCKkqmD+Q6MEWFCOAwa3R/NExsA69vb+cy4cXwVOCGhW/NrzyYwRoYJ4VbQ+kPUfr+TBx54gDbg8TFjOCFGbFCz1RVKNkEomBBuBe0/RK33OznkkENwMe9PAvonN0LIJggFu2fJVqiEe2TEoVj3xrjyyitZu3Zt7HqyVtf8+fFFOiv655yT+avlvjNJv8YqBbMIt4J2FzQuxdqI9R7nOOmkk7xYhUmf3LDVJDowIRwGzS5oXIqxESvAL4DD4jfPK5onmGw1SfkxISwBWq0R8LsRa5Ya4AgPVfokpDXLxsgxISwB2oP+Xmhqgih3sA/4v5jVaV45o9XCNwrHhLAEBLGiobMTmptZ3NXFfwGj9t234KrU9jFCs4VvFIbNGg/BsmWwaFGT15nF1laYN6846RZq6OzkP/fck4dFuPbaawuuRnUf8TvTu2wZfO97emaxQ8UswkFkrZFUalcWLdJnjWhf0fDII4/ErkPzxAZYMnUSMSEcRNYaSadFpduj/Ue48847c/bZZ/ONb3yj4DoqIW3JkqmThQnhILLWSCqVpq6uSp01Anp/hGeddRYvvvgi3/zmN2MJIfjpY0dHI6edptfi0m75hoTFCAeRtUaOPbbTyw9HawyoGCsafvKTnwAwbty4+JV5oL19rOpYo8+VM0Y8zCIcgtZWSKW6aW3dLVY9mmNA3lc0NDdz1d//ziVA69FHx2+gB1pa1qtPW7Jkah2YEBYR7TEgXz/Cjxx5JKxdyxeAzwC1t9wC3/9+/IpjMnnyhuSnLRleMCEsIr5neBctamLUKH0/wIZog4VqoBaoWb26rO3JJelrlQ0/mBAWEd8zvFpTejbuuCOj167FAQKxdqQGnStnglgdFDAmhEXG5wyv1pSeR268kZlf/SrS3U11U1OsGzVpdUF9xlS1z2aHiM0aVwBZa6SqKu3NGvE9k/2R8eMZVVvL/ddcE6sezatKfK0O0j6bHSJmEVYAWWtk4cJOjj12N5VB/0cffRSATZs2xarHpwva0dHIsmX63E/fs9lGfEwIKwRfKT3FDvoffPDBsc73GVedO3cq/f363E+fs9mGH0wIA6MYu1J/vKuLfYEH4jcP8BdX7eurIp3WOctr+YO6MCEMjGLsSi3AYuBrPhroiZkzobY2TX9/tc3yGsNiQhggvnelFjL5gzrWk2RobYVLLlnJhg3TLJnaGBYTQqNwcnalrgaat9++3C3ajMmTN3iZiLBk6uRj6TNG4XR2wqRJOKBq0iTe+8orBVeldXMKsFtuhoBZhEY8OjupFqFuzRpSBVah3fX0vUGF1rSekDEhNGLR29sLxMsfrATX0+eNn7Sm9YSMucZGLM466ywARKTgOny7nlk3u6OjMV5FRSCb1mOrSnRhFqERi+uuuw6AHXbYoeA6inWXv5qaqUybpsvi8p3WY/jBhNCIxbp16wBUbM0Pm7vZzok6N9tnWo/hDxNCIxbOOQDOPffcMrckQ+7KmZoap+4uf+AvrcfwhwmhEQvnHL/4xS/K3Yy3yHWzGxtX0to6raB6tM9kG36xyRIjNhMmTCh3EzYju13W5MkbCq5D83Zghn9MCI2C+fa3v42IMGvWrHI3xTuWRB0WJoRGwXz/+9/HOcdPf/rTcjfFO75vtal55YxhMULDA/X19eVuQlHwmUSdG2+86KJGszCVYRahYRSZwfHG9vax5W6SMQgTQsMoMoPjjS0t68vdJGMQ5hobRpEZvHImlSp8NtsoDiaEhlECcuONloqjD3ONDcMIHhNCwzCCx4TQMIzgMSE0DCN4TAgNwwgeE0LDMILHhNAwjOAxITQMI3hMCA3DCB4TQsMwgseE0DCM4DEhNAwjeIYVQhGpF5EVIrJSRDpE5DtR+SIReUZEnhCRhSJSG5UfJSKPR482EZmaU9d/RHU8ISI3ikh9VL6/iPwpKr9ORGqichGRy0Tkuai+aTl1HS0iz0aPo31/MIZhhEM+FmEK2N85NxVoAT4lIjOARcAewD8DDcDx0fEvAPs65z4EzAeuBhCRnYE5wJ7OuSlANXCEiFQB1wFHROVdQFbYDgZ2jx4nAldFdW0LnAfsDewFnCci4wr9EAzDCJthhdBl6Ime1kYP55xbHL3mgBXALtHxbc65ddHxy7PlETVAQ2TxbQO8CGwHpJxzf4mOuQ/4YvT/IcDPordZDowVkYnAQcB9zrnXove6D/hUIR+AYRhGXjFCEakWkXbgZTIC9HDOa7XAbOCeIU49DrgbwDn3N+BioBtYA7zunPst8HegVkT2jM45DHhP9P/OwKqc+lZHZVsqNwzDGDF5bczqnBsAWkRkLHC7iExxzj0RvXwl8KBz7g+554jIfmSEcJ/o+TgyFt6uwHrglyIyyzl3vYgcAfxAREYBvwX6s9UM1ZytlL8DETmRjFvNhAkTWJrnrpg9PT15H1sqNLYJrF0jxdo1MkrSLufciB5kYnOn5vx/B1A16JgPAc8D78spOxy4Juf5V4Arh6j/k8DN0f8/Bo7Mee0ZYCJwJPDjnPLNjtvSY/r06S5flixZkvexpUJjm5yzdo0Ua9fIKKRdwKNuBLqWz6zx+MgSREQagAOBp0XkeDKxuiOdc+mc45uA24DZ7u24H2Rc4hkiso2ICHAA8FR0zg7R31HAGcD/ROfcBXwlmj2eQcadXgPcC3xSRMZFluYnozLDMIwRk49rPBG4TkSqycQUb3bO/VpE+snM8C7L6Bq3OefOB84lMwFyZVTe75zb0zn3sIjcAvyJjOv7Z6IZZeA0EflsVP9VzrnfR+WLgU8DzwFvAscAOOdeE5H5wCPRcec7514briOPPfbY30WkK48+A2xPJn6pCY1tAmvXSLF2jYxC2jVpJAdLxoo0BiMijzrn9hz+yNKhsU1g7Rop1q6RUYp22coSwzCCx4TQMIzgMSHcMlcPf0jJ0dgmsHaNFGvXyCh6uyxGaBhG8JhFaBiGMZKkw0p4AGOBW4CnyeQptua8diqZFSjbR8/fDfwKWAl0AMfkHDsAtEePu3LKF5FJ7H4CWAjURuUCXEYm1edxYFop25Xz+uVAT87zUcAvonY9DDSX+PMS4LvAX6L652j4vMjksf4pKn8IeO//b+9sQq2qojj+WyAqGn6gz+SR8niaYhMho+xDAuNNDKLigQ4KCZ0YQUQ40EGODMVBBQlCQlKD6Ets8BIH2YOCUqSQS0b0lOsXSYUpSKOHu8Fa17s579wTHLx7X99df9jc/XX2XvzPOmt/nLPXTczXcvQU1a/AuVY/6MmrU8DvJsfMHpErt96XylVX76c8N7kN190OqCeb7RafCSyw+DL0o+uLEfG7gf0WHwCuR4p3q0P7m+zmC/AJsCPKP27564FTKeWyskeAjwsK8SpwyOJbgE8T8/UK8BF2+ghY0gt8oYZ5TcTRkcR8jQMjFr8PmGPxz1BPTKAHC3b0iFy59b5Urrp6P0Ufum2YUgZgHuoGTErKvgDWAs2I+F3oWWlBR+IJ2g9sR4MTtfkGsNfipccBU8mFujX7Fv0APlaIE9gojH5A/zftveEUcp3GZluF/Nx8/QY8Fl3/diq+gIeA70vaEOtvhqUfB07kliu33lfJRQ29LwvTbY9wGPgL+FBEfhaRwyIyV0SeA66GEM4W6r8PrEHdgTWA10P7uOBsETkjIj+KyPPFjkq87lR5xEkh12vo0u+PQlt35AohTAI30ZM/qeRaAWy2suMi8mCP8LUd+FpErqD3cV9CvlYBN0TkqLVzwE5uLQJuWL9FTnLKdQeZ9L5Krjp6PxX/Z/3vpYBOkSdpj/TvAQfQPYL5ltekPQKNAu+gI9BKdPSaZ2WD9jts16wo9PUB8G6UHgOeitLfAOtSyAUMovtcrZlEPDL+AjwQpc8Di1LxBdwC3rT4i8B3ufmy9NGo/Z3A4VR8Wf5Nk2kG8CXqqWkAmIj6XgY0csuVW+8r+Kql96W2I7fxupsBWAo0o/QGuzF/GuFNuzGXrO4YsCGqfxJ4tKTdI8BolN5DwesO1UuErsoFPAtci9q6jT1QVC+pus4Xukk+ZHFBHWfk5msAOB/lLwfOpeIL3Usbj/JfBg5Sf2ncVbly630FX7X0vixMq6VxCOEacFlEVlvWM8BPIYQlIYShEMIQOnV/2OpesjqIyP3AauCCebWZZfmLgSfRN1V08rpDZ085XZcrhDAWQlgatfVvCGFlJNdWi48CJ4NpRwq+0Adno8WfRl9SZOUL+AeYLyKrrP0RzBNSCr5QZyELRWTA2tlo9zGg+12jlr8V+Cq3XFYvm95X8FVL70vRyULeqwH9X5Uz6Kv8Y8DCQnmT9lR8EH0l30A/C3jJ8p+wvLP2uy26fhKdZrc+yXgrtGc7B62sgf43SzK5Cm3FS4TZwOfoxvNpYDgxXwvQkb4B/ACs7QW+gBeisvEWLyn4srIRa6OBzlRbb0eHrd8Jk2NWj8iVTe+r5Kqr98XgJ0scDkffY1otjR0Oh6MO3BA6HI6+hxtCh8PR93BD6HA4+h5uCB0OR9/DDaHD4eh7uCF0OBx9DzeEDoej7/EfIdRBpAjQiVUAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# translate the center point according to the delta_D\n", "from math import sqrt, sin, cos, tan\n", "# from geopandas import GeoDataFrame\n", "import pandas as pd\n", "import geopandas as gpd\n", "# from shapely.geometry import Point\n", "import shapely \n", "\n", "#####\n", "# FROM ABOVE:\n", "### perpendicular line passing through the center\n", "print('slope of perpendicular to fit is: ', -1/fit[0])\n", "m_p = (-1/fit[0])\n", "# intercept of perpendicular line\n", "b_c = center_point_y -m_p*center_point_x\n", "y_p_points = m_p*x_points + b_c\n", "#####\n", "# test, put points at given distance along the perpendicular line\n", "\n", "# delta_D = delta_TEST\n", "\n", "# create a dataframe \n", "# columns = ['x_UTM','y_UTM']\n", "# dictionary_append = []\n", "# geom2shift = pd.DataFrame(columns=columns)\n", "\n", "'''\n", "# create a dataframe \n", "columns = ['x_UTM','y_UTM','sPoints','delta_D']\n", "dictionary_append = []\n", "#geom2shift = pd.DataFrame(columns=columns)\n", "#geom2shift.set_geometry('sPoints')\n", "geom2shift = gpd.GeoDataFrame(columns=columns,geometry='sPoints')\n", "print('delta_D is: ',delta_D)\n", "'''\n", "\n", "\n", "# fix by Mikhail Minin start\n", "# to define the data frame\n", "columns = ['x_UTM','y_UTM','sPoints','delta']\n", "dictionary_append = []\n", "#geom2shift = pd.DataFrame(columns=columns)\n", "#geom2shift.set_geometry('sPoints')\n", "geom2shift = gpd.GeoDataFrame(columns=columns,geometry='sPoints')\n", "# fix by Mikhail Minin stop\n", "\n", "\n", "\n", "\n", "\n", "# delta = delta_E\n", "\n", "print('delta is: ',delta)\n", "j=0\n", "for element in delta:\n", " print(j)\n", " print(element) \n", " for i in range(len(x_points)):\n", " y_shift = 0\n", " x_shift = 0\n", " x_shift = (element * cos(m_p_angle))\n", " y_shift = (element * sin(m_p_angle))\n", " # print(x_points[0])\n", " # print(y_points[0])\n", " # print('x_shift is: ',x_shift)\n", " # print('y_shift is: ',y_shift)\n", " x_shifted = x_points + x_shift\n", " y_shifted = y_points + y_shift\n", " # print('x_shifted is: ',x_shifted)\n", " # print('y_shifted is: ',y_shifted)\n", " geom2shift.loc[j+i,'sPoints']=shapely.geometry.point.Point(x_shifted[i],y_shifted[i])\n", " geom2shift.loc[j+i,'x_UTM']=x_shifted[i]\n", " geom2shift.loc[j+i,'y_UTM']=y_shifted[i]\n", " geom2shift.loc[j+i,'delta']=element\n", " j+=len(x_points)\n", " # calculate x, y for delta_d elements\n", "\n", " \n", " # plot each point shifted\n", " # points and fit\n", " fig.set_size_inches(4.5, 10.5)\n", " ax1 = plt.plot(x_points, fit_fn(x_points), '--k')\n", " fig = plt.gcf()\n", " ax2 = plt.plot(x_points[24],y_points[24])\n", " plt.plot(x_shifted,y_shifted, '.', color='blue')\n", " # plot perpendicular line\n", " ax3 = plt.plot(x_points, y_points, '.', color='red')\n", " plt.grid('on')\n", " plt.axes().set_aspect(1) # it's turned off, but lines are perpendicular, uncomment to check ;)\n", " \n", " # plot perpendicular line\n", " ax3 = plt.plot(x_points, y_p_points)\n", "# for i in range(len(x_points)):\n", "# geom2shift.loc[j+i,'sPoints']=shapely.geometry.point.Point(x_points[i],y_p_points[i])\n", "# geom2shift.loc[j+i,'x_UTM']=x_points[i]\n", "# geom2shift.loc[j+i,'y_UTM']=y_p_points[i]\n", "# geom2shift.loc[j+i,'delta_D']=element\n", "# j+=len(x_points)\n", " plt.show\n", " \n", " # populate the pandas dataframe with new results\n", "# for geometry in geom2shift.geometry:\n", "# print(geometry)\n", "# add geometry\n", " \n", " # add attributes based on element\n", " \n", " # save\n", " \n", "# # create empty dataframe\n", "# columns = ['x_UTM','y_UTM']\n", "# dictionary_append = []\n", "# geom2fit = pd.DataFrame(columns=columns)\n", "# # print all rows x,y\n", "# # fill a dataframe for each row\n", "# for geometry in gdf.geometry:\n", "# print(geometry.x)\n", "# print(geometry.y)\n", "# dictionary_append.append({'x_UTM':geometry.x, 'y_UTM':geometry.y})\n", "# geom2fit = geom2fit.append(dictionary_append)\n", " " ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(['x_UTM', 'y_UTM', 'sPoints', 'delta'], dtype=object)" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "geom2shift.columns.values" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# export data as GeoJson / shapefile / excel" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# export geopandas as json/shapefile\n", "\n", "# outjson = '2_output/test.geojson'\n", "# outjson = '2_output/agpa_D-shifted.geojson'\n", "# outjson = '2_output/agpa_E-shifted.geojson'\n", "\n", "\n", "with open(outjson, 'w') as f:\n", " f.write(geom2shift.to_json())\n", "\n", "# export to csv\n", "\n", "with open(outcsv, 'w') as f:\n", " f.write(geom2shift.to_csv())\n", " \n", "# df.to_file('/dev/sdc/Boundary.shp', driver='ESRI Shapefile', crs_wkt=prj)\n", "\n", "# and to csv as\n", "# geom2shift.to_csv()\n", "\n", "# to convert to shapefile you'd need an appropriate driver installed, but it should look something like\n", "\n", "# geom2shift.to_file(driver ='ESRI Shapefile', filename='abc.shp')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Load raster over the same area" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# TBD" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Pick raster values for each point and add to geopandas" ] }, { "cell_type": "code", "execution_count": 119, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# TO BE ADDED LATER \n", "# SEE https://github.com/mminin/snippets/blob/master/getElevationAtPoint.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# visualize materials on leaflet" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import mplleaflet" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAALAAAAD8CAYAAADXLS5JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAGiBJREFUeJztnX2wXVV5h5834RIuUkkgAUMUEzp8iF8JpBFwoApqwDolIs7EVk0dqNVqa9sxFdoZaeM46qQzdhCVYRSNo1UpiRExNGYElBn5aCAhCcSUGARzoRWJwQ8ykIS3f+y1c/e92efsdc5ae++1zlnPzJl7ztr77L3Ovb+7zrvf9e71E1UlkYiVKW13IJFwIQk4ETVJwImoSQJORE0ScCJqkoATUZMEnIiaJOBE1CQBJ6LmiLY70CQzZ87UuXPntt2NhAX333//r1R1VtV+QyXguXPnsnHjxra7kbBARB6z2S+FEImoSQJORE0ScCJqkoATUZMEnIiaocpChMbaTWOsXL+DJ/bu46TpoyxffDpLFsxpu1tRkQTcEms3jXH1mq3s238QgLG9+7h6zVaAJOIeSAJuiZXrdxwSb86+/QdZuX7HIQGnEbqaJOCWeGLvvq7tNiN0ErjFRZyIHCUi94nIgyLykIj8q2n/hojsEJFtInKjiIyY9j8XkS3m8RMReW3hWB8x+z8kIn9XaH+naXtBRBYW2ueKyD4R2Wwe1xe2nS0iW0Vkp4hcKyLi65fSBCdNH+3a3m2EhnGBj+3dhzIu8LWbxg7tv3bTGK//9O3Mu+r7vP7Tt0/YNijYZCGeAy5U1dcC84GLReQc4BvAGcCrgVHgSrP/o8Afq+prgE8ANwCIyKuAvwQWAa8F3iYip5r3bAMuA35ccv6fqep88/hAof2LwPuBU83jYruPHAZvPKN8mj9vrxqhfQh8EKgUsGb8zrwcMQ9V1XVmmwL3AS81+/9EVX9t9r8nbwdeAdyjqs+q6gHgR8DbzXu2q+oO206LyGzgxap6tzn/14Altu8PgTt++lTX9qoR2lXgg4JVHlhEporIZuCXwAZVvbewbQR4D/BfJW+9ArjNPN8GXCAix4vI0cBbgZdZnH6eiGwSkR+JyPmmbQ6wu7DPbtNW1vf3i8hGEdn41FPlommDKgFWjdCuArchhhDESsCqelBV55ONpotMOJDzBeDHqnpX8T0i8kYyAX/MHGM78BlgA5nYHwQOVJz6SeBkVV0A/APwHyLyYqAs3i1doUVVb1DVhaq6cNasyuq8xjh2dKS0PRdg1QhdJfCq41cRSwjS00ycqu4F7sTEmyJyDTCLTFyHEJHXAF8CLlXVpwvv/7KqnqWqFwB7gEcqzvdc/n5VvR/4GXAa2Yj70sKuLwWe6OWztMnaTWP8/vnD/3dHpgjLF58OVI+g3QRuc/wqbEKQEEboyjSaiMwC9qvqXhEZBd4EfEZErgQWAxep6guF/U8G1gDvUdX/mXSsE1T1l2afy4BzLc69R1UPisgpZBdru1R1j4j81lxM3gu8F/hcD5+7VVau38H+g4d/YRxz1BGH0mDHjo6wd9/+w/axCRFsjl9FLGk+mzzwbGCViEwlG7FvUtVbReQA8Bhwt8lgrVHVFcDHgeOBL5j2A6qap8ZWi8jxwH7gQ/nFnoi8nUyAs4Dvi8hmVV0MXACsMOc6CHxAVfeYY30Q+CpZBuQ2xmPt4Okkjr3PZoK1GUG7Cbzq+DZU/QNVTcQ0NdNYKWBV3QIsKGkvfa+qXsl4Sm3ytvM7tH8H+E5J+2pgdYf3bAReVbYtdGzE0W0ErRL4yvU7GCsRcS/xr2uI09RMY6pGaxgXceQjaJXAqy7wqrANccqwzYL4ukhMAm6YOsWRC7wqg1GFrxCnDNuZRluSgBumCXG45oBtprldQhygNMTppY85ScANYzMyuYrDNQfc7zS3bYizdtNYaSK/lz7mJAE3SBPxr48ccF3T3MXPUDbrJGDdx5wk4AZpIv5tIgdc1zS30nuKLQm4QdqMf3vNAXc7h+s0dyeBz+kxfIAk4EZp4uLIRw1EndPckIUJI1MmRsG9hDhFkoAbpImLIx81EHXngIHDy7H6vB0hCbhBmrg4qiv+9Z0DntzP/Qe1r1rlJOAGaeviyGf8G1IOGJKAG6Wti6Mm4t82csCQBNwYPkYmV4FX0VSaz1cOGJKAG8PHyOR69V9Fm2m+fnLAkATcGD5GpkGogej0GfrJAUMScGO4jkyDUgPhmuabTBJwA/gYmQalBsI1zTeZJOAG8CG+GGogmpjmnkwScAP4EF/oNRBNhDhlJAE3gOtXbww1EG3Ev5AE3AiuExAx1EC0Ef9CEnAj2ExATM4Bj45MDSL+tf0WaGKau4wk4Abo9tW7dtMYq+8fm5ADFuAdZ88JKv51vRPadZq7E0nADdBNHGV35yrjo3Po8a/tndCu09ydSAKumSpx2CwQEnr8W+xvp8/hOs3diSTgmqkSRydh5O0xxL/F/k7G53KvZSQB10yVODoZI+Ttg1IDXEcOGJKAa6fqD9dJaHn7oNQA15EDhiTgWvExMg1KDXAdOWBIAq6VOmsgBqEG2DUHDEnAtdJmDUQMNcCu8S8E7hNntl1tvOB2iMjiQvvFpm2niFzl/JuoAdeRKdUAVxO0T5yInAksBV5J5svxBeOYNBX4PHAJcCbwLrNvMNR1+zmkGuAiofvEXQp8y5i9PArsJPsHWATsVNVdqvo88C2zbzAMSw1wWzUQOaH7xM0BflF4nfvBdWoPhmGoAYb67b6qCN0nrpMfnLVPXFtGh8NQAwz1231VEbRPHNnIWhylcz+4Tu1lfW7F6HAYaoChfruvKmyyELNEZLp5nvvE/bTgE/euXnziCvtcBnyz4vS3AEtFZJqIzCPzibsP+G/gVBGZJyJHkl3o3WLzgZvC9at3EGog6o5/wW4Eng3cISJbyISzQVVvBa4HTiTzidssIh83+xd94jaLyMbCsVaLyMPA95jkEyciu8mMD78vIusBVPUh4CbgYbKw40MmnDkAfBhYD2wn8657yOH34B3XCYhBqIGoqwa4SNA+cWbbJ4FPlrSvA9Z17HjLuDhtQjYSf/2exw/bXgxB6vaC67TUle23wMbH9nT9DD5IM3E1MAg1EC5LXfmy+7IhCbgGBqEGIoSlrmxIAq6BQaiBCGGpKxuSgGtgEGogqoxY2q6ByEkC9syg1EBUGbG0XQORkwTsmUGpgQC6GrG0XQORkwTsmUGpgagyYmlrHYjJJAF7ZlBqINq2O7AlCdgzg1ADEYLdgS1JwJ4ZhBqIWHLAkATsnUGogYglBwxJwN5xHZlCWAciBLsDW5KAPTIoNRAhpPlsSQL2yKDUQISQ5rMlCdgjMdRAxDDN3QtJwB4JXRyxTHP3QhKwJ2IQx6DFv5AE7I0YxDFo8S8kAXsjBnHEMM3dK0nAnohBHDFMc/dKErAnYhBHDNPcvZIE7IkYxBHDNHevJAF7IgZxhD7N3Q9JwJ4IXRwxTHP3QxKwB2IQRwzT3P2QBOyBGMQRwzR3PyQBeyAGccSQ5uuHJGAPhF4DAXGk+fohCdiRGGogII40Xz8kATviQ3yuq0Da4Bpjh1YDkZME7IgPGyrXVSBtcImxQ41/IQnYGR82VK53AFfhKsBQ418Ix+jwOBHZICKPmJ8zTPsbROQZs9J7cRX4YIwO67Kh6uUO4CpcBRhq/AvhGB1eBfxQVU8Ffmhe59ylqvPNY4U5VjBGh64XR66rQNrgKsAQayBygjA6JDMpXGWerwKWVHQrGKND14sj16t/G1wvwEKsgcgJxejwRFV9EsD8PKFwjHNN+HKbiLzStFkbHdbtE+c6AeFa5FNFDNPcLoRudPgA8HITvnwOWJsfvqybHfpem0+cD3G4FvlUUWear80aiJxQjA7/T0Rmm/fOJhvpUdXf5OGLcSUaEZGZ9GB0WCc+xFE1wrZdA+Fjob86CcXo8BZgmXm+DPiu2e8lIiLm+SLT36cJxOjQx9X5G8+YdZhARkemWhf5VBFDms+FSp84MqPDVebKfwqZqeCtInIAeIzM6BBgjckSFI0OAQ6o6kJzrNUicjywn4LRIfBp4CYRuQJ4HHinab8c+KA51z5gqbloPCAiudHhVODGNowO+/WCK9ZArL5/bIJABHjH2XMmFPl0O0cVrn5zTaT5XAjF6PBp4KKS9uuA6zq8p1WjQ9v4t5v4Vq7fwb79BydsU8bDg6ZqIPoRuM80nwtpJq5Pmrg4CqEGook0nwtJwH0yDDUQ4F7kUzdJwH0Sw8VR3Wm+tuNfSALum2Gogaj6HG3Hv5AE3DfDUAOR97dTmq/t+BeSgPvGNXZ0dcK0wcetTt3SfG3WQOQkAfeJl5swHZwwq/B1q1O3NF+bNRA5ScB9UNd9cL04YVYx6DUQOUnAfRCDOAa9BiInCbgPYhBHDGk+HyQB94HrxVET4oghzeeDJOAeqXNywKc4Ykjz+SAJuEd8xL+uTpg2+KiB6KfMs8kcMCQB94yP+NfV7MUG1ylimzLPbsdviiTgHvFxceS60F8VPtaBqLvM0xdJwD1S18VRLwv9VVHXOhA+yzx9kQTcI3VdHNku9GfDIK8DMZkk4B5xnYBI60D4JQm4B3xMQKR1IPySBNwDPiYg0joQfkkC7gEfExCuo1+/fQxpmtsnScA94DoB4cPspd8+hjTN7ZMkYEuacpt3EUcs09w+SQK2pAm3eVdxxDLN7ZMkYEuacBpyFUcs09w+SQK2pAmnIVdxdOpj3t52iFMHScAWxBD/AkiH9EHeHsI0t2+SgC2IIf4tHqtTewjT3L5JArYghvjX5hwhTHP7JgnYgiZsWF3F0YSheEg1EDlJwBa4LkLSRPzbhKF4SDUQOaH7xImIXGu84LaIyFmF9ywz+z8iIsuoG4dFSJqIf5sIc0KqgcgJ3SfuEuBU83g/8EVzrOOAa4DXmeNdk4u+DlwXIWni4qjNNF8bNRA5ofvEXQp8zZzmHmC6MYFZTGb3tcecawPGeKYOQq8BjiXNVweh+8R18oOz9olzJYYa4FjSfHUQuk9cJz84a584V6PDGGqAY0nz1UHQPnF09oOz9olzNTqMoQY4hjRfXQTtE2fa32uyEecAz5gQYz3wFhGZYS7e3mLavBNDDXAMab66CN0nbh1ZrLwTeBZ4H4Cq7hGRT5AZHgKsUNU9ff0GutDkxZGLFxxQmebrxwsu9PgXwveJU+BDHd5zI3Bjx457wIf42loHIk/zLVkwx9nsMNT4F9JMXFdiuTgKPc1XJ0nAXYjh4iiGNF+dJAF3wXVkaqoGIvQ0X50kAXfBdWRqswYipDRfnSQBd8F1ZGqzBiKkNF+dJAF3IfQC8WGugchJAu5ADAXiw1wDkZME3IEYCsRjSfPVSRJwB2IoEI8hzVc3ScAdaFMcqQbCniTgDri69DQmjsBvdaqbJOASfLj0tF0DAdUxduzxLyQBl+LDpaeJHLBLDcQgxL+QBFyKD5eeuo0CXWsgBiH+hSTgUnzEjnUbBbrWQAxC/AtJwKX4mCKu2yjQtQYixnXQykgCLqGuKWKfRoGuZtwx1wAXSQIuwbVA3HWdXhtcY+yYa4CLJAFPwkeBuOs6vTZ9rCvNF0MNcJEk4ElUXRyt3TTGlA4KzcXhuk6vTR/rTvOFXANcJAl4Et0ujgCuXrOVg3q4xPOvbxuB1+UF10SaL4Qa4CJJwJPoViBeNvIBTBXhU5e9GqgWOLTvBddEmq8pkoALVH31dhLGC6osWTCnUuB5eqtqhHbpI9hNEded5muKJOACVV+9/U5w5AJfu2nMaoR26aMPL+RYcsCQBDyBqq/eOiY4YOII7drHOr2QQ8sBQxLwBNq6By4foX30cdC8kKtIAjbEcA/cMHohV5EEbIjhHrhh9EKuIgnYEMM9cMPohVxFErAhhnvghtELuYokYEMM66ANoxdyFUnAhhjWQRtGL+QqmjY6/HtzjG0i8k0ROcq0XygiD5j2VSJyhGl/g4g8IyKbzePjhWNdbM6/U0SuwpGY10ELZamrNmjS6HAO8LfAQlV9FTAVWCoiU8i84Zaa9scY98sAuEtV55vHCnOsqcDnyYwQzwTeJSJn9vtLgPDFEUOarw2aNDqEzNJg1IywR5M5Cx0PPFcwhNkAvKOiW4uAnaq6S1WfB75FZorYFzGII4Y0Xxs0ZnSoqmPAv5GZuDxJ5jj0A+BXwIiI5EYwlzPRQutcE77cJiKvNG3WRoc2PnExiCOGNF8bNGZ0aOywLgXmAScBLxKRd5sRfCnwWRG5D/gt4waIDwAvN+HL54C1+eHLutmh75U+cTGII4Y0Xxs0aXT4JuBRVX1KVfeTecmdZ457t6qer6qLgB9jDBBV9Td5+KKq68hG6pn0YHRoQxNG2VUVYFXEkOZrgyaNDh8HzhGRoyUzkLsI2G7ekxsgTiMbsa83r19i9kVEFpn+Pk3mD3eqiMwTkSPJRvBb+vkFNLFItE0FWBUxpPnawGYEng3cISJbyISzQVVvJRPZiWRGh8UUV9HocLOIbAQwcfPNZGHBVnPuG8x7lovIdmAL8D1Vvd20Xw5sE5EHgWvJMhVq3O4/TObOuZ3MfPGhfn4BNuKrKkD3UQFWRQxpvjZo2ujwGuCakvblwPKS9uuA6zocax2Zk6cT3cRnU4CeC7xsH9sKMBuqDBVt0nzdzA5jzAFDmonrGvf5ugeuahGSKppI88UY/8KQC7ipe+CeLTmH7S1EUH+aL9b4F4ZcwE3dA/frSV/D00dHrG8h6nYeX2m+WONfGHIBt3UP3Ium9SaMunPAsca/MOQCbqv+oZeJAZsQxDXNV1VnHDJDK2Afkw+ui/zZ9LEqBPGR5quqMw6ZoRVwnbef2y7yZ9PHbiFIL2m+MmzXcguZoRVwE7efuwrDZg20OtN8oafQYIgF3G+GwacLfBX93gPXdJqvTYZWwP1mGHy6wFdR1z1wvtN8bTK0Aq77DgsfudUm7oHzkeZrk6EVcN1Trz5yqzGk+dpmKAXsknqynXp1jX9jSPOFwNAJ2EfqqYn4N4Y0XwgMnYB9pJ46jXw+498Y0nwhMHQCdk09QfkNedBs/BtCmi8Ehk7ArqmnJtbXbWKVyZhLKIsMnYCXLz6d0ZGpE9p6mZVqYn3dJlaZjLmEssjQCXjJgjl86rJXM2f6KEI2YhWT9v0K3Of6uk3cYRxzCWWRynviBpElC+Z0FFPevnL9Dp7Yu4+Tpo+yfPHpEwR+9ZqtE+LkXkobbai6/81XgXu3c8TCUAq4CheBnzR9tDRL4Tv+7Sa+YYl/IQm4L7oJvGqErqLORVQGLf6FJGDvVI3QVfhaRKUsUzJo8S8kAddCtxG6iqoQxNcExyDEvzCEWYjQiSHNFxJJwIERQ5ovJFIIESChp/lCIgk4QtpM84VGEvAAUmeaLzSSgIcM1zRfaCQBDyEuab7QCN0nTkTkWuMFt0VEzioca5mIPGIey0gMJ6ra9UGWAz/GPB8B7gXOAd5qtgnwTeCDZp/zgBnm+SXAveb5HDIPuVHz+ibgL8j+iX4BnGbaVwBXmOdvJXM5EnPO/FjHAbvMzxnm+Yyqz3L22WdrIg6AjVrx91TV4H3iLgW+Zk5zDzBdRGaTeXNsUNU95lwbMMYzieEidJ+4Tn5w1j5xicEmdJ+4Tn5w1j5xNkaHiXjpKQuhqntF5E6yr+ttBZ+4vyruV/CJu0RLfOLMPrlP3NdV9W7gfNP+FuA0855OfnC7gTdMar+zQ59vYNyv+SkReayXz+yBmWTfMqESav9ebrVXVZBMJtDp5vkocBfwNjInop9gLsoK+58M7ATOm9T+OuAhsthXyAy+/8ZsO8H8nAb8kMxcHOBPmHgRd5+OX8Q9SnYBN8M8P84m6G/6geXFSOpffw+bEXg2sMo4xE8h82S7VUQOkDnL3228CNdo5iZf9IkDOKCZ1eu9IpL7xB0ANjHRJ+5t5vhf1HGfuHVkmYidwLPA+8w/3R4R+QSZbx3AClXdY/FZEgOGmP/CRE2IyEZVXVi9ZzuE3r8qUjll/dxQvUurhN6/rqQROBE1aQRORE0ScAdEZLqI3CwiPxWR7SJybmHbR0VERWSmeX2siHyvUC/yvsK+B43p+WYRuaXQ/lURebSwbb5p71j/UWMfTxaRH5hjPCwic037PBG519SbfFtEjjTt08zrnWb7XNffd9+0nQYJ9UGW5rvSPD+S8VTiy4D1ZBmYmabtn4DPmOezgD3Akeb17zoc/6vA5SXtpfUfNffxTuDN5vkxwNHm+U3AUvP8esbrXf4auN48Xwp8u62/UxqBSxCRFwMXAF8GUNXnVXWv2fxZ4B+ZOPOnwB9Iljc8hkwch985aUen+o9a+igiZwJHqOoGc5zfqeqzZr8LgZvN+1cBSwp9XGWe3wxcZPZvnCTgck4BngK+IiKbRORLIvIiEflTYExVH5y0/3XAK8hmCbcCH1HVF8y2o8xU9j0ismTS+z5pwoTPisg002Zb5+Grj6cBe0VkjTnOSpPzPx7Yq6r5P2KxH4f6aLY/Y/ZvnCTgco4AziKbVFkA/B74F+CfySZqJrMY2ExW4zEfuM6MkAAna5Zn/TPg30XkD0371cAZwB+RzSx+zLTb1nn46uMRZNP4HzV9OYWszLVbP6xrUeomCbic3cBuHa+6u5lMLPOAB0Xk52T1Fw+IyEvIZgjXmK/9nWRT22cAqOoT5ucuslhzgXn9pNn/OeArwKLCucvqP+rq425gk6ruMqPpWnOcX5GFL/lsbbEfh/poth9LFpI0ThJwCar6v8AvRCS/0/Ei4AFVPUFV56rqXLI/4llm38fNPojIicDpwC4RmZGHBiYb8HrgYfN6tvkpZLHlNnOuW4D3mmzEOWRlp0/W1Uey6fgZIpJbGl0IPKzZFdodZOWtAMuA7xb6mN8Fczlwu9m/edq6egz9QfY1uxHYQjYqzZi0/eeMX+GfBPyALLbcBrzbtJ9n2h40P68ovP/2wv5fZ/yuFwE+D/zMbF9YZx/NtjebY2wly47k2YlTyG5W2An8JzDNtB9lXu80209p6++UZuISUZNCiETUJAEnoiYJOBE1ScCJqEkCTkRNEnAiapKAE1GTBJyImv8HuywkeiaFOSUAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ax = geom2shift.plot()" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# +ellps=GRS80 +no_defs +proj=utm +units=m +zone=28\n", "crs_leaflet = {'ellps': 'GRS80',\n", " 'proj': 'utm',\n", " 'units': 'm',\n", " 'zone': '28'}\n", "\n", "mplleaflet.display(fig=ax.figure, crs=crs_leaflet)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Appendix & TBD\n", "\n", "raster conversion to UTM GRS80\n", "\n", "```\n", "# dtm\n", "\n", "gdalwarp -t_srs \"+proj=utm +zone=28 +north +ellps=GRS80 \\\n", "+datum=GRS80 +units=m +no_defs\" \\\n", "dtm-lanzarote-wgs84-utm28n.tif dtm-lanzarote-grs80-utm28n.tif\n", "\n", "# same for shade relief\n", "```\n", "\n", "conversion of geojson to shapefile (table to contain all info on point locations)\n", "\n", "```\n", "ogr2ogr -nlt POINT -skipfailures D.shp agpa_D-shifted.geojson OGRGeoJSON\n", "\n", "ogr2ogr -nlt POINT -skipfailures E.shp agpa_E-shifted.geojson OGRGeoJSON\n", "\n", "```" ] } ], "metadata": { "anaconda-cloud": {}, "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.6.3" }, "toc": { "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "toc_cell": false, "toc_position": { "height": "535px", "left": "0px", "right": "1068px", "top": "111px", "width": "212px" }, "toc_section_display": "block", "toc_window_display": true }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }