Skip to content

Instantly share code, notes, and snippets.

@apatlpo

apatlpo/demo_drifter_on_swath.ipynb

Created August 5, 2024 13:37
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save apatlpo/4216d25529e7bb6e08545b325c368f07 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save apatlpo/4216d25529e7bb6e08545b325c368f07 to your computer and use it in GitHub Desktop.
Download ZIP
swot drifter swath demo
Display the source blob
Display the rendered blob
Raw
demo_drifter_on_swath.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "e23ab5be-17da-466c-9af6-915ac66f8189",
"metadata": {},
"source": [
"# demo: find drifter under swath"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7e7df255-b27a-4101-92b4-0f98c63c6fff",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import xarray as xr\n",
"\n",
"import matplotlib.pyplot as plt\n",
"\n",
"#from pyproj import Geod, Proj\n",
"import geopandas as gpd\n",
"import shapely as shp\n",
"from shapely.geometry import Polygon"
]
},
{
"cell_type": "markdown",
"id": "68216b2d-ded3-495a-89c7-31283776fae2",
"metadata": {},
"source": [
"### utils"
]
},
{
"cell_type": "code",
"execution_count": 63,
"id": "c4d39372-93f9-47b6-8a9f-3c41d5f6ce25",
"metadata": {},
"outputs": [],
"source": [
"\n",
"def build_swath_polygon(da, x=\"x\", y=\"y\", dix = 50, diy = 50, side=None):\n",
" \"\"\" build a shapely polygon of swot swath\n",
"\n",
" Parameters\n",
" ----------\n",
" da: xr.Dataset, xr.DataArray\n",
" Input dataset where coordinates x/y are found\n",
" x, y: str, optional\n",
" Coordinates to use to define bounds\n",
" dix, diy: int, optional\n",
" Subsampling steps in num_pixels and num_lines directions\n",
" side: str, optional\n",
" take either the polygon for the full swath, or the right/left ones\n",
" \"\"\"\n",
"\n",
" # right swath\n",
" if side==\"left\":\n",
" da = da.isel(num_pixels=slice(0, da.num_pixels.size//2))\n",
" elif side==\"right\":\n",
" da = da.isel(num_pixels=slice(da.num_pixels.size//2, -1))\n",
" \n",
" nx = da.num_pixels.size\n",
" ny = da.num_lines.size\n",
" \n",
" rg = lambda size, di: list(range(0, size, dix))+[size-1]\n",
" \n",
" X, Y = [], []\n",
" \n",
" e = da.isel(num_lines=0, num_pixels=rg(nx, dix))\n",
" X += list(e[x].data)\n",
" Y += list(e[y].data)\n",
" \n",
" e = da.isel(num_lines=-1, num_pixels=rg(nx, dix))\n",
" X += list(e[x].data)\n",
" Y += list(e[y].data)\n",
" \n",
" e = da.isel(num_lines=rg(ny,diy), num_pixels=0)\n",
" X += list(e[x].data)\n",
" Y += list(e[y].data)\n",
" \n",
" e = da.isel(num_lines=rg(ny,diy), num_pixels=-1)\n",
" X += list(e[x].data)\n",
" Y += list(e[y].data)\n",
" \n",
" mpt = shp.MultiPoint([[x, y] for x, y in zip(X, Y)])\n",
" poly = mpt.convex_hull\n",
"\n",
" return poly\n",
"\n",
"#def is_inside(poly, *args):\n",
"# if len(args)==2:\n",
"# pt = shp.Point(*args)\n",
"# else:\n",
"# pt = args[0]\n",
"# return poly.contains(pt)"
]
},
{
"cell_type": "markdown",
"id": "1da75077-eab8-4354-b70c-238ef7bfb3c5",
"metadata": {},
"source": [
"### load data"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1c935367-4077-40b0-8ae7-2b5f89a346b6",
"metadata": {},
"outputs": [],
"source": [
"swot = xr.open_zarr(\n",
" \"/home/datawork-lops-osi/aponte/swot/cswot/swot/l3_v1.0.1/3_500.zarr\"\n",
").compute()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c08ef0d5-5d57-453c-b718-4edec4b5d0e5",
"metadata": {},
"outputs": [],
"source": [
"dr = xr.open_dataset(\n",
" \"/home/datawork-lops-osi/aponte/swot/cswot/drifters_harmonized/L2/all_med_variational_1h_v0.nc\"\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3c1981cc-0382-4ed3-9e76-f30556b40d19",
"metadata": {},
"outputs": [],
"source": [
"# build swath polygons\n",
"both = build_swath_polygon(swot, x=\"longitude\", y=\"latitude\", dix = 50, diy = 50, side=None)\n",
"left = build_swath_polygon(swot, x=\"longitude\", y=\"latitude\", dix = 50, diy = 50, side=\"left\")\n",
"right = build_swath_polygon(swot, x=\"longitude\", y=\"latitude\", dix = 50, diy = 50, side=\"right\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "23f272c9-8ba0-4698-89f5-919d28b37e71",
"metadata": {},
"outputs": [],
"source": [
"# stack dimensions\n",
"drs = dr.stack(points=(\"datetime\", \"drifter_id\"))\n",
"lon, lat = drs.longitude.fillna(0.).values, drs.latitude.fillna(0.).values"
]
},
{
"cell_type": "code",
"execution_count": 67,
"id": "51e29b44-a46c-45ae-9b90-85fd13042c77",
"metadata": {},
"outputs": [],
"source": [
"# test if points inside left swath\n",
"drs[\"inside_left\"] = (\n",
" \"points\", \n",
" np.array([left.contains(shp.Point(_lon, _lat)) for _lon, _lat in zip(lon, lat)]),\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 76,
"id": "67c199e9-4cba-4a1a-a377-b5ecbbe4845a",
"metadata": {},
"outputs": [],
"source": [
"# put mask inside original dataset\n",
"dr[\"inside_left\"] = drs[\"inside_left\"].unstack()\n",
"\n",
"# masked dataset\n",
"dr_left = dr.where(dr[\"inside_left\"])"
]
},
{
"cell_type": "code",
"execution_count": 81,
"id": "f54a63cd-13bd-418a-a9c7-5564da00ff96",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAADlCAYAAADzy5WHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+C0lEQVR4nO3dfXxT9d0//tdJmrSktCVxIqJtEUWcdygkww0VxDtUftul7kKaTVG/1340osNd7pqA48KbH19keA+auOEUnZThpiI4RRGZN0OuVry5lHkzhRZFvGtta0ubNPn8/uhOlptzkpzknJObvp6PRx7Q5PTkk57knHc+n/fn/ZGEEAJEREREJrHkuwFEREQ0tDD4ICIiIlMx+CAiIiJTMfggIiIiUzH4ICIiIlMx+CAiIiJTMfggIiIiUzH4ICIiIlOV5bsBiSKRCPbt24eqqipIkpTv5hAREVEGhBDo7u7G6NGjYbGk7tsouOBj3759qK2tzXcziIiIKAt79+7F4YcfnnKbggs+qqqqAAw2vrq6Os+tISIiokx0dXWhtrY2eh1PpeCCD3mopbq6msEHERFRkckkZYIJp0RERGQqBh9ERERkKgYfREREZCoGH0RERGQqBh9EZDqv1wtJkqI3u92e7yYRkYkKbrYLEZU2u92OUCgUd1/iz0RU2tjzQUS68ng8cb0aiTe1QEOSJHi9XpNbS0T5kHXwsXLlSkiShG3btik+ft1110GSJOzZsyfbpyCiIhIIBGCxWNDS0pL1PpqamlIGLkRUGrIKPvbt24fbbrtN9fE333wTDz/8cNaNIqLiUl9fD5/PByGEoc9jtVoRCAQMfQ4iMl5Wwcc111yDhQsXKj4WiUQwb948LFmyJKeGEVFx8Hq9aGtrM+W5IpEIfD5ftCeEwzRExUlz8LFx40bYbDbMmDFD8fFVq1bhtNNOw/HHH59z44ioMMXOVmlqajLkOYQQ0VtdXZ3iNrHDNAxEiIqHpuCjp6cHN9xwA+68807Fxz/99FM88MAD+O///u+M99nf34+urq64GxEVLrvdrlvAUVdXB5vNlna71tZWNDQ0pMz7aGpqgsvl0qVdRGQsTcHH4sWL0djYiEMPPVTx8WuuuQbLli2Dw+HIeJ/Lli1DTU1N9FZbW6ulSURkEnkWi57TYkeOHJly9kustWvXIhKJRHtD3G530u90dHSwZghREcg4+HjjjTewY8cONDY2Kj7+1FNPoaysDOeff76mBixcuBCdnZ3R2969ezX9PhEZKxAIoLy8PKdZLGpy2WdzczOEEPD7/XGBSigUQn19fdwsGavVqkdziUgnksgwPf2WW27BE088EV3mvq+vDzt27MCECRMwYsQIWK1WdHd3R3s9vvnmG7z11luYPHkyKioqsGnTJgwfPjzt83R1daGmpgadnZ3R5yIicwUCAcyfPx/BYDCv7dAyeybdVNyGhgasXbs21yYRkQot1++Mg49Ee/bswRFHHIEXX3wR06ZNS3p827ZtOOOMM7B7926MGTMm4/0y+CDKD4/HY0jvRi70DD6sVisGBgZybRIRqdBy/WZ5daIhqhCDDSOFw+F8N4GI/imr4OPaa6/Fa6+9Fv3/Mcccg3Xr1kUfnz17Nt57773o/0855RTcddddubeWiHI21IKOWJIkGV4IjYjSy3rYxSgcdiHSXyAQwHXXXYfe3t58N0WTuro6tLa2ZrSty+VCR0dH2u0kSUIkEsm1aUSUQMv1mwvLEZWgQCAAu90ene3h8/mKLvAAgLa2tpS5HPL0X4vFklHgAWjLIyEiYzDng6jIDYVhFKvVmpSzUV9fHy3rni6gsNlsutYnIaLcMPggKlJDIeiQRSKRrFe1dTgc6OnpgcVigRCCq+MSFQAOu9DQE+oD3lgHrLkQ2LwYOPBNvlukmcvlGjKBRy78fj96enoAIFodlfkeRPnH4IOGllAf+m85BNgwF9i9Fdh+D8St9cBjVwA9X+W7dWkFAgFIkpRxfsNQJ6+Ay0XniAoLZ7vQ0PLhC8CjFyk+JAQgzbwPeHo+gDCACAALcOkTwJHTTGxkMq/Xa9jqsZmqq6uL5lgUI6fTifb29nw3g6hkcbYLkZoxU9Cv0usuSYDY5AMQwmDgAQARiId/CNxYA6w+C+j8xKSGDpJnc+Q78ABQ1IEHMLjonMViQSAQyHdTiIY89nzQ0BIeAJrXQjxzdVaJh0IIeJocaPlgvwGNi5dp3QrShnU+iIzB8upEiXq+Al5chm9aHsIIDGQ940GSJCyY0Kn6+3rF8h6Ph4GHQQrs+xbRkMRhFyp9PV9B/OZIoGU1RiD3hcXueFV9pdfYZdwTb5kkPQYCAVgsFs5kIaKSxuCDSt+zC6FXaYfW9gH8z/7suuybmprigpFY8iwWn8/Hb+Ym8Hg8+W4C0ZDGYRcqbV++B/H2H3UrLPVVXwRhnWIDFrvKH/YsEeUXez6odIQHgNbtwLpLgUf/HXjsSgRXTtbtIi+EwEWP9umyL0rP7XZDCAEhBJxOZ877S9xHbC+UxWJhLRAiE3G2C5WE6uGV2OYFJo42rjPvoy+DOOo+Bh9Gc7vdaG5ujrsv1wBSXh033X7kUuxEpB3rfNCQEQgEUFZWhvHD+wwNPIQQmPEwAw8jyb0ciYFHruTAIxO9vb0oLy9nLRAigzH4oKLicrmi3eRygmY4HMbPTzQ28PjxQ934x7eGPcWQZbVa4ff7UybZulyurPff0NAQF3j4/f643o+6urqk3wkGg5g3b17Wz0lE6XHYhYpCqu7yURXAvl9VGZrAecWTvXjordyn6ZLysEoqWo+r1lOa2v619JgQEYddqERUVlYqTktNtO7HFTkFHkrXqtj7whGBP77DwEMPeg2r2Gw21ce0Jo5aLMqnwWIvJ09UyBh8UEHyeDzo7e3NaNuTDs3tuSQJCCaU7lj/PrBqRx/mPd2Lqv/bjQPh3J6DBns89BIKhVQfW7duXdrfr6+vjwa2qUqtcwYMkTEYfFBB8Xq9kCQp4zoMdgtQPUz9W3AmwgKwnzQ77r5LjgF8nnI8+MYAA48cOZ3OnHo8hBCKuRlqhg0blvLxysrKjHs1CmFBP6JSxOCD8i72W6jWk/3UMdaccz2sI04EDnQl32+RcMnxrMOXKz2WsW9tbc241sftt9+e8vFMe9SIyDhMOKW8yyV4sFuAjgVVcNiy24cQIvr8sf8HBnM9OOSSOz1PMZm8V9I9n9EJrERDFRNOqSjI02W1kqdmioEQ+lub4fhBI2CpyaoNsc8vSRLw3Qtx1VY7fraBuR7FSMvwTCb0qKxKRMkYfJDp5MJg2XyjFEJg2bJlKLNIaG6sAVafCTTfD0Q69WncyONxovcWrH6TuR65kAuGmd1rkMnU2IaGhoz25XQ6o0NGHo8HkiTlVIAscZVjoqEs6+Bj5cqVkCQJ27ZtAwAMDAxg9erVOOOMMzB9+nRMmjQJV155Jb744gu92kol4qqrrkI4rP3KLveUtLW1YZzLAs9h+udjiG23oPGSGbrvNx29v7HrqZDbFivTdq5duxZ+vz/tdh0dHdFAQU6ADgaD8Pl8cUFEZWVlVgEJZ9LQUJZVzse+ffvw/e9/H21tbXjxxRcxbdo0fPLJJxg3bhx27NiBE088Ef39/Zg5cyb6+/vx0ksvZbxv5nyUPj2+9Z0/xoKn5wzXoTXJIkJg+FLzh1wkSSqZ/AKjXkd9fb3qTJVsntPj8ei2wq3D4cDtt9+OxsZGxceV3vdcS4ZKieE5H9dccw0WLlwYd5/dbseVV16JE088EQBQXl6OuXPn4uWXX8a+ffuyeRoqUbmOo4+vAjZdVqlTa5JZJAk/NbBcu5pSCTyMpHfhr+bm5ujwUK49PL29vfD5fKq9IEr758wbGqo0Bx8bN26EzWbDjBnxXdMjR47EvffeG3dfRUUFgMGuSiJZe3t7TgWnnrs8t4qmmdjfzUCgmOiRGNra2hqXq6J2S0dtWOarr75K2tbhcOTcbqJipCn46OnpwQ033IA777wzo+23b98Ot9uNMWPGqG7T39+Prq6uuBuVvmwLTrnsQK0zt6Ji6fQGBTZ/zGzTQhIIBFTLoAPx+Rl2u93QtmgNRuQeEaVeDg650FClKfhYvHgxGhsbceih6etZf/XVV1i9ejVWrVqVcrtly5ahpqYmequtrdXSJCpSSj0XmfSG3HVuuaG9Hn0DAoes6E4qt06ZMSJvxePxwOfzZbzfVKXX9SBX4dVj1gpnvdBQlXHw8cYbb2DHjh2qyVSxBgYGMHv2bNx8882YPHlyym0XLlyIzs7O6G3v3r2ZNomKlFqW/0cffZTy96wScO7RxuRiDESAf1/Xg5pl3fiWa8hlLdU6Kdmor69XTAhtaGgwvAaH1+uFxWKB3W6H3W7XVIWX9UGIUsv4TL5p0yYcOHAA06dPBwD09fUBAK699lqMGDECq1evxlFHHYVIJII5c+Zg6tSpmDt3btr9lpeXo7y8PMvmUyFSm0FgsVhSXpw6OjpS7vc6jwUjh1tzbp+SMgtw/enleOIDJgBmy4iEWaUEU7/fH/0SZETPgdfrjQswtPakuN3upGFF9nAQxcu6vPqePXtwxBFHRKfaynw+H6qqqvCb3/wGALBlyxaMHTsWY8eOzWi/nGpbvPSctpjIbgH6fl1l+En86JXf4sN2jrmkI4RAZWVlXB6DEcFH4vFOvLAnPp5NGwKBAK677jocOHBAcR8222COUWwQ4na78frrr0dL8qfr8Un1vuUsJyoVeSuvvmDBAvz973/HrFmz0NLSgpaWFqxfv1736XFUWOTqj0YFHoA+C8il87+fD+DjDgYemZAkKS/TRLNNVFbj9XqjyaBKCaQNDQ0IBoMIBoNxSabNzc2IRCIQQmQ01JQqwHC5XDm/DqJik1Xwce2112L27Nlx/3/33XexfPly/PWvf4XH44nefve73+naYCocgUDA8KBDtuuzsOHfEBv+1Icwv4QWrEzKomupGpo4vALE91A4HA6sXbs28wamIYRQfA3phhuJShFXtSVNAoEA5s2bp3tiYSp2C9C9sAr2MuN6PnqDEVTf+i2DjxwYMcslNrBV2r9SxdNM2qFWKdWs06Eew0VEhUbL9dv8Mo5UtJS+KZrh4vFWQwMPAGjeO8DAIwdqvRKxF9nYhdrSSVVGXeZyubLuNVDat5kzVORcEZnL5cr4b0NUCriqLaUl53TkI/AAgAWnGVtUTAiB2X/qM/Q5Sl1TUxOs1tQzkWILgaXicrkyyhNTCjzkAmOBQABWqzX6fPJwTH19fdLz+/1+CCHyevHn0AsNNQw+SJUZiaTpfM8FnDDKuOBDCIHj7+rGfsYeOYtEIrosG692IU63b5vNhrvvvhvA4MrJsUODTU1N0RWRE/l8PsXni63tIUlS2uBKq8S1XrjKLQ0lDD4oiZmJpOk8+3+MXcclFAphFyv6Fz15Vopc/yOXHAqXywVJkpLqeygFV5nc1GaztLa2xv3c1NSkuigdUalh8DHEVVZWRhe/kql9EzTb8DJgxDBjh1woP2IvzvX19Tntq6GhIWlWSjb5G3LZdL2HQOThJqXXmbikgM/nYw8IDQkMPoYwj8cTrdXQ29ury1oVerr2FJvh7fl76orulCG/3591wmZbW1vaoRq1aao2m01xOmx7e3tGU3OBweAlVU5TqgXttGhra4sL8oHBuiWJf7empqacAzKiQsfgo0Sk6vb1eDzRRLvYWyEMq6Ty1QHjp5/8f28buwjZUOHz+aLf8I2iFGSEQiHVXIy1a9cqBiByoTC/358y6HA4HBBCIBwOJ61km+nN4XDE7bO3tzdpGKa9vT0p/6Otrc3w1XmJ8onBR5FRCiLSnfBbWlqKssrs+v8dQDhiXAAihMBTH4QN238pk0uOJ9K7XoXD4YDf74/+rLTysZyLoZQvkRiwWCyW6FCj0kq5drs9OvtFj+Xue3p6kno2lIZ1Wltbk15bKBRiAEIli8FHgUsMNooxiMhGjQ34ckEVrBbjvkk/834QQVZTz0owGEz6tp4pteEQq9Wa1HPQ09MTt5J2c3NztNcikc/ni/YqxC57HysSiaiWhW9oaEB/f39GK3dr0d7erhg0JWpubk56XVoXtSMqFqxwWqCMXKSt0NktQPeCKthtxgUeQgh859ZutAcNe4qSJp82LBaLLr0dkiRh9uzZmsqZ51JkTGa323H33XfrHnAoiQ2EUv3NEov5KSXUEhWivC0sR/ophMAjX8mnZx9pNTTwAIDO3hADjyzF9lzIi6vF3jL5li+TJAl+vx+RSETzBba9vV11+CeV2KEVI3o6cpX4d1i/fn2eWkJkHJZXLwCFNMMkVuIqn5mUvNbDGaON/3s07zX8KUqOzWZDMJg+Ymtubk7ZIyJJEpxOJ5YuXZrzhV9uT7rPkN/vz2uQkctnfNasWTq2hKgwMPigjJkReBxcDvzn1GGGPocQAj95kiVNtaqqqsp4WzMXHgQGg4urr74a4fC/EogzDZaMlji9Np3ExFm1HqFUXwYKbDSdKAmHXYaguro6xWmAiWK/reldWlrNfTONrWgKALs+C+HLfkOfomg0NDRkNEwiSRKWLl1qQouy09jYiIGBgbjhn0IIPAAkJbimStT1er0ZFflLl3wem6SuV50SIj3xXVkAzFpN0+FwwOFwRIs6qWX9y4QQqKysRCAQMOWbrFUCjjoo+S2p57c4IQTOf5S9HrKmpqaM8ovuu+++gsuNKEZ1dXVJZdVjJdYcUfqCoPWLgLyCbuIie0T5xNkuBaJQ8z7MNM5lwQfXDDf0OfqDQVQsY/ChVYGdJopKulkuibNbZErDRmqz4IQQmmfIORwOXWqZEMk426UIpRsCGQoOsUQMv8ht/9jQ3ZcEdtPrJ1WZ9EAggLKyMsXAw+l0Kg4bJQYXTqcz+pmRa6DIt3Q9qr29vSzjTnnDs4zBEsucy/8mSjcEUuqGlwEvXVVlaA+QEAKzmGiqyu12QwhherJoKUvMywgEAnEVVmMTZGVutxvt7e1J9yeWZbfZbIrbydrb2+OCkUzaR2QWBh8mkr+1tLS0aCqPbiSbzVYQXeq/+L7xi8h19jDRNJWWlhZ+E9ZZYnKpz+dLWWFVCIHm5mbFxxMLqmlNqE0VhBCZjcGHgYrhRB4KhZKm9slduWaeqPr6jX0uIQTOeIC9Hunwm7B+vF5vRn9POehIV2Qtdhgl1YyZQCAAq9WqutAkUSFgwqkBSu0DHpuYZkShsYPLgc+vN3bI5fn/6cM5zxTG1MtiVGCniYJntVpTDl8ZUTJdLXE1FR5X0pOW6zeLjOlMaWXNYtfb22toYHCvCbU9Xu3lSTYXLpcrZX4BxVMLPPS82GcTbMTKpjQ9kV447KKzq6++Ot9NyJtsVzldtLnP0G9gQgjc/ipXB81FR0dHNGE69mZW8blSIUlSzsOx8krX6QIPeThH7VYoRdhoaOKwi45KbbjFDMPLgK/+qwrlduP+dnvagzhiJfM9jGKxWBRnbQxFhXwOcLvdqsmsRHowpc7HypUrIUkStm3bFnf//fffj4kTJ2LKlCm44IIL8Omnn2b7FEUhEAgwkStLdgvw5fXGBh4AMHcDv+EZiVNzBxX6OYCzmaiQZBV87Nu3D7fddlvS/Y8//jiWLFmCZ599Fq+++iomT56MmTNnluzJKdN1GEjZeeOsqCgz9oTdG4zghb2l+f4j0qqtrQ1lZWUlmZtGxSWr4OOaa67BwoULk+5funQp5syZg5EjRwIA5s+fj3feeQd/+ctfcmtlAcpkzJVSO2+s8fkCbR0DCBfUwCINdVarVTUPI3aRP6fTGZ1eGzv9XeutoaEh7vnD4TC/NFHeaQ4+Nm7cCJvNhhkzZsTd39HRgZ07d8ZV76ypqcHRRx+NLVu25N7SAmG32wu+e7VYHK2wiJyehBA4/xHmeuTKrIUPi51SAKFk1qxZqo/Flkhvb2+PVinNZabR2rVrFWe2KFVaJjKLpqm2PT09uOGGG7B582b098eXivz448FFM0aNGhV3/6hRo6KPKenv74/bV1dXl5YmmSrd3H3S5tVPwjjjyNTbCAFIVccCh9Zh4P1nUaYhXnnw9X7s5rpZOUusrEmpNTc3K35BUVoozizBYDCpTVoWoSPSm6avnosXL0ZjYyMOPfTQpMfkksHl5eVx95eXl6dct2TZsmWoqamJ3mpra7U0yTQMPPRllYDHd4UzmGIrgG/3Ah/GBB4TGgDPfwCHnAS1ORZCCPznc0w0NUNitz5Bcfgj31Nb0/XIEJkp4+DjjTfewI4dO9DY2Kj4uLwqa2KPSH9/f8oVWxcuXIjOzs7obe/evZk2yTQWi4WBRwp2ux1+vz/jkuxWCdj+fxzY6RuO/oHU2w9+W+uOu0+8uRaYeCng+yus899Kek4hBKTGlxCyGbdSsN/vz7quSalpamqC1+vNdzMoDU6zpUKS8bDLpk2bcODAAUyfPh0A0Nc3OJZ+7bXXYsSIEVixYgUAYP/+/XG/t3//fpx99tmq+y0vL0/qLSkkzO9QpxRoZHIRGueywHPY4FuvwqY970OSJGD1OQDswPizII3yAJ+/g1c/6kJHCLjyyT58efPJmverBRP24jU1NeleLpz0xXMZFZKMz/yLFy/Gzp07sW3bNmzbtg3r1q0DANx1113Ytm0bPB4PTj755LhxxK6uLnzwwQc466yz9G+5weT6HaRMLQlRbQZQXV0d/H4/AKBCj0ku4X4g3A3segL4vAVAH34w1oZfPN3HlWuJElgsyad6+fNIlA+6Tjf49a9/jTVr1uDLL78EANxzzz04/vjjcf755+v5NKbgN1t1dXV1GWff22w2CCHQ2tqKxsZGCCHwxuN3p/ydbGvuSpKE935RhWGs+K27dOuAcAiqsCX2UjqdTtUhdCIzZLWw3LXXXovXXnst+v9jjjkG69atw0UXXYQvvvgC5557LioqKuB0OrFx40bFqLuQVVZW5rsJBSmTfA5JkuK2U0yy+/LvSfuN7WXKpcPJapFwyfFleOitgex3QklCIfW1cQpshQZKoLT+DhcJpHzLKvi46667VB9rbGws6oja5XKlnJ0zFMhrQMQGBJnOaMgoMbfykLgfBwOW3IIOWTgi8Md3GHgQAYNfpBI/k6zbQoUgq+CjVA31HA+lWSOGOPXnEK/cHhds6PWn3/pxCAe4xhkRACh+kWKvBxWC4hoPodIwbASkhsd0360QApf9mRVNiVKRF8JUGo4hMguDDzKc1+tFWVnZ4DTczk+AP/w78Ld7dO9Z+f32Puxn7GEqJpoWtlTHh7WLKJ847BJjqBcTq6+vR2trq277CwQCuOqqqyCEgN0C2N5aD3HHpujwlt7DXKve4niL2fR8v5D+WltbVT9nxTYRgEoLgw+Kamtrg8fj0a0Sojxd2W4BOq6vgsNuXE5NcCCC//1y6AaORGo4G4kKEUPfGEO510PW0tKS0zcil8sVHVOW/XCc1dDAAwB+/kQvwjzHEhEVBQYfMdgNOUiuu5FpvZPKyspowKG0AurcE1IXqMpVOBLB6r8zcCQiKha82sYIh8PRlR/dbnfGC6WVqt7eXtTX16fcRpKktHVR3Kl3kTP/y/3s9cgDrpJKRNli8JGgubkZQoi4vAchxJA90ba1tSnenzi0ouZXJwE1lcb2fDz/KXs9zOb3+7lKKhFljcFHhuSgROlW6iRJQnl5OQKBQPTnTN1yXoWhxduCIYFnP+IsF7PNnz8/+n4gItKKwQdlJBgMwufzaQokrBLw7lfGvsW27QkhyI4PwzkcjrifY98P6YbmiIgSMfgwSeLJO1vFNAR0zEEWnDza2NncN72osHCdzuT8n2L5uxshVV5PW1sbLBYLe0KIKGMMPnKU6QyZ/v7+rBZ0cjgcSUM8LS0tmveTCb3LLdfWGDu9NhgcwI79+nV7+P1+xWE0ObchduiN4gkhonVdiIjSYfCRo0wvROFwGB0dHdGLVyarxLrdbvT09CTdb7fbNbczE+FwGJIk6TbleOvuML4N/uvvE/u30tJfIYRAv8L9W94J6zbLpa6uTtNqzJyWTUSUPZ5Bc5BLIuXatWsVu/GdTmc0QEmcTeDxeCBJEoJBfYYalIaChBBZFVuz2ZJntAQjwEHLu4EzfwPU1ENqfAm4sRO4sRP2xV8Dl21C6km6/2wTJJR7/iMut0MIgUs3K4Uk2Vm4cKGm7WOnZQ81kiQpHm8iokwx+MhSpgW4UlEaPlFb7trlcuk+3JKuPkem3G43QqFQ3H1WCTjhYAt+NrEM6PsM+J4PcB0Vs0EZMPY0OBZ/Dcx+LGUPhkUC0LwaZRJw1cZ+PLizH9+5tRvtOqZ7XH311aqPeb1exfsznWrq9/vh9/t1y/vJNyEEQqEQHA4H/H5/vptDREVIEgU2gN3V1YWamhp0dnaiuro6381RlWmvh1wtVAu32510YTNyuqrerBLw2n844E5INo0IwPKf7wI1hyf/Uvd+iNvGI93LFELoHniQfgrsdDLkxJ4nGhoa8OGHH6KlpUXxnEKkNy3Xb/Z8mEDrsuNGJZQaTT7xjXVakgIPYLAHQ9xxHPDxy0B4IP7BqlGQFrQCkxuR6volSRJ+c065ns0mHenRI0jZSfzbNzU1Rc8lLS0t0cKAkiQZljdGlCkGHybIZtlxl8uV8nGbzVZwXd7yt97PuiL4uH1AcRtJAvDwTHTdMgY48E38g8NGAOctHwxCJl2BkELqiRACv3pOv1wP0ldvb29R9dKVCovFomkYNRQK8ThRXjH40Mjr9WY800Frj0esxAXaEgONUCiEq666Kuv9G2WYFfhmYRXGusoGezBsLsWZKtXoRuTWeqB7v8JORgD/z12wLf4c+OGqaE+IEAJH3s4hl0Qc6hja6uvrs34PMAChfGHwEcPr9aKsrCyaYOj1euO6KiVJQlNTU9oPujxbRUuPR7p9NjY2JtUJKcSLzhUnlcFqGTyhSRKAUDvUBkksEhC5bTyw52/JwzAAYKsAJl4K6Vcf4bZX++Bc1o3dyTOPhzxeQIYml8sFSZJU11+Kxfo0VGgYfGDwm4McWITDYTQ1NUV/1kqttyPV1MRMTwpVVVUZtyOXk01sUTOt00mPcmp7S1kkAA+dh55bRgJ/WQBsvgF4/mbgk5Z/BSSV38GSv5WhM5RyV5QFzlgpPoFAAJIkJfWOZkLpnJD4BSv2xhweMgqDD6iv3JoNtd6OYDCoOSBILFedTTszKWaWitKCemokScI9b1dkFfRUIgz8jx/Yvgp49XZg9Zn45pYjgHc3AaE+/OhHP8rlZZCK3t5eXSuTqlWJpdwEAgHY7XZIkpTyeKn97SVJgsfjSbmNEr2m4xMlGrJTba1Wa1bFtFJxOp2qdTpiZdpNbrVaMTDwr+EILd3rsYdV6ff8fj9+/vOfJ9XnAAYDlrVr16ruW60ddgtw/lFWXDqxEheNdwCjxkF89lba6bPp9AQFXMu7uYBcgaurq8squbpUeDwexZlquf5d1PZrlgK7RFABM2yq7YYNGzBz5kycffbZOPXUUzFp0iSsX78+bpv29nZcccUVmDhxIqZNm4ZTTz0VL730kvZXYTCtgYckSapl0RsaGiCEyCjwkLfPhFzuXL5loq6uLulkoTR04vP5FAMPYHCKXiAQUMx5USu4ZbcA7ddX4YmGSlw0HgB6IT57E9I5dwIT5wDH/ghwHYVIFuexSruEGUfqu+4M6W/kyJH5boKp5OFa+aYWILS1tWU9jBEIBFT3m0tCO1G+aer5mDFjBrxeLy677DIAwMaNG/Fv//ZvePPNN3HCCScAAC699FJ89NFH2LZtG+x2OzZs2IBLL70UH3zwAUaNGpX2OYzs+aivr89q6MKIAj1erzernBKt5KDJaDOOtOKZnyqfWAcAlB32fWDWaqBiBLDjd8CuDUDXfkR6PoMlg7jq1y8cwNJXmPRR6IZC70e255FE6XoY7XZ70heExN/RI9mYBchIL1qu35qCj9dffx0TJkxAWdlgAanu7m5UV1fj8ccfx4UXXggAOOGEE3D++edj+fLlAAbHDCsrK+O20avxWmX7QTX64p3vbtVcqVU0TSSEgPT/bgUOi+mJ6f8W2BEAXv8jQp0fQCktVwiB6v/bjW+VS4dQgXE4HIoLIpYCl8uVMtFTkiTcd999aGxsRCAQSJtPoxasVVZWJuVb2Gy2pHWdEs9psecqpX0AgzVBwuFwynYRZcOwYZdJkyZFA49QKIQVK1bg2GOPxdlnnx3d5uKLL8YzzzwT/YD+4Q9/AAAccsghml5EvsjDFpkkWOqlubk55awSm82W1+mUbrc7OuQU2450FU0TSZIE8dvpwFPXAfveHJzNUj4cOP2XwC+aYVv8NXDZBmDk8dHhmeDAYDn1Qg08bDZbUU1jdDqdGQ37yUOJQgjNs2HkQmOlNFMiEAjAarUqBh7y50NelFFeHbmxsTHuPKJUH0gekqmvr4/e53K5koIGp9OZdkHJxHNIT09P0vR8m83GwIMKQlYJp/PmzcOjjz6K4447Dn/84x9x+OHxa3XccMMNWLVqFUaOHIndu3fD5/Nh5cqVivvq7+9Hf/+/ylB1dXWhtrbWkJ6P2CTTQuxq9Hq9WLduneqFzOFwmJZ9nunbIhAI4Fe/vA5bvcgoAIkVBGD/wQIg3A3sfAQ4dwVw8o8HF50L9QF7XgXGTAFsFaYNU2VD/lvpHSD6/X40Njbq3rWej0A202TsQqQ2zJLtOUStR0KNUo+HLPZYFksATKXLsGGXWOFwGDfddBPWrFmD1157DYceeigAYNGiRXj66afx/PPPY+TIkdi6dSv27t2LOXPmKO7nxhtvxE033ZR0f6EvLGcktZNdbBdtukAlV5r3Gx4A9r+Nzt+djxocyPp5+wBU2A4CLl8fHZ7Ra4zdKLHj8EO54JeW/KJCDP5jBQIB1dlgVqsVq1ativZwAOrDMWqvM93wTSwGFVQsTAk+gMEZI2PGjMEll1yCFStW4Msvv8Shhx6KNWvW4Cc/+Ul0u6OOOgo333yz4kwJM3s+ionaeHG6w6VH/oj8jTsr4QHg83eBPa8g+Owi2C1AMALYJG2r+woh8N07u/F+d3bNMJtRvR9mc7vdWb9/hBAIBAK46qqrNF8w0yVfGi3Tz438BSCXz1lsQJIqCGFuBhUbw3I+Erv+LBYLxo0bh127dgEAdu/ejXA4jDFjxsRtN2bMGPzpT39S3Gd5eTmqq6vjbvSv8eJE6abdNjc3ZzRGHztOnXjLOvAABodMRk8AfjAP9sWfAz95HPbFn0P65QfA4adnfFGSJAlPz6nIvh0m0zIdutCpvS/SkQtgZfN9Rq4qnGo6t1aBQACVlZUpK3immyYrk3PBRo4cmdH2qbS0tETzYdrb21X/3gw8qJRpCj4mTpyYdN9nn32G0aNHAwAOO+yw6H2J2wwbNizbNg5panP5E0+esRJ7TBwOR9KJzZQub1sFMO7MwX+rRgH/sXEwCDnuIqD2dKS6RgkhcMGaPuPbOAQ0NDSkrAkRm6goXxgTq+uaKTYQyeXm8/lyzpFyu92w2WzRxNBUQUdisnpismes3t7etCtXE5U0oYEkSWLTpk3Rnx955BFhsVjEyy+/HL3vnHPOEaeffro4cOCAEEKIp556SgAQTz75ZEbP0dnZKQCIzs5OLU0raU6nUwDI6maz2fLdfHW9HUJsuVmIZxcJ8fSvhFhzsRB3f198cl21GF+V3evlbfBmt9uF3+8XDQ0NSY+53e6k+xwOh+J+nE5n9HAp7Wuo3/x+f8Zvd7/fr7gPt9ttwIeLyHxart+acj5WrlyJpqYmWK3WaPXNRYsW4YILLohu8/XXX2PBggXYuXMnKioq0NfXh/nz50cLk6VjVnn1YqU1U17D4S04hTzDhZTV1dUVXHKwPLVYr/dSLrkYSoXDYpVyjRQqfaYlnBqBwYc2qRLfinl6YyrpZr/IyYulkoNRKoQQhgeUkiRh9uzZ+PDDD3Ut3KfnZymT92VOSd9EecLgg4YEpV4g+ZtjJtUlyVyxp5rEoNlisei+0GM6qZ7T6NNiJgFIgZ2aidIybLYLUSFRquDY29sLj8eTFHiINLM25IqeapQq3xZaz4pQWfgwsd2FcFFrbm5OmtmR6YKLubBarfD7/dHnlJ/fjOeOVWjHg8hsDD6oqLW3tydNLU7sbnc4HHE/K81CkGdYqM1QkGc7xJo9e3bSdloDEqfTmXEZc5fLBb/fn3IGVOKQhtpFNfHiZ8YKqZIkJc3w8Hg80dkpZuT3hMNh+Hy+pJkxzC0iMheHXagkpEvETXybp0v8S0ePxEqlPIJc25UNuXBWIfXkxC7QJktVdVRvZp4W1QqNFdipmSgtDrvQkJNqhoDSRTUYDEa/9ada1E+NHjM6Ojo6kgprxbYrk2JxFotFccEyLZR6dfJNCJHUQ+Hz+bIKPOx2uyk9O1rZ7XZIkqQYeJg9DERkNgYfVPKUhkdkcuJjXV0dbDabia2KJw/7xFb3XLZsWdrfi0QipidqFpv+/n60trYq5r74/X7Y7fak3zFyRV656qpaIOV2u/Naap7IDAw+qGQo9RSkOpHX19dH80Pa2tpMH+5QElvds9DqZRSi2EBCa29BfX09fD6f4oqxvb29kCQJZWVluld7TTU82NDQUNAL7hHphcEHlYwHHngg7me1FUWtVqvuF3chRMY9J3KSaary24kkSYrO0shkOGYoSAw21q5dq6n3KvH4Kw3NxCao6rHmjMfjUbxfnk3FHg8aKphwSiUjMW9BCKFp6XKt0tUUkRNKY9uVj3oWxS6bU1QgEMCtt96KBQsWqBbrUnq/aC2CprVtSrk1BXYKJsoai4zRkKMUAGitxeFwOHDsscemXTystbUVgPZptZQsn6efxOOXOPsok5lHWkutKwU8RKWCs11oyEkMPDKZwRI7S0TuwlcLPOS8AjnwUEpSzJVc6EyPuht+v181B6IQZn7YbLa8X3gT67/Ezj5KlRAai71YVCxia+roNYyYCwYfVJJaWlrS9kzIFw75YqyWCJh4kQKgS3JqbB5H4nh/7OwMpedX43a7IYRAY2Mj1q5dGxeASJIUDaCy7bWJDZBib5nmr8ivVynJ02w9PT2a/ray2JybXKc5s/eMjGaxWCBJUtIXq6ampqSif6ZKu+6tybQsyUskg0FLptfV1SU9V11dXc77bWhoyPi1+f1+4XK5hCRJqvuz2WyalncXQoiGhoak/TgcDs37kSm1qxg4nc6Uxz/2Z6fTmfPzFePfiIqD2+3WfC7Sk5brN3M+qCQUyzfIXJZjl2WSTJkPcs0UtVlGNCgxP6nATsFUwIxYFVrP9x8TTmlIKrQARJ7ZwosxxWLSKWlh5Iw9IH/BB3M+qGQIhVyETG65JGDG5gwk5kPIq6Yy8CAiLbxebzQxNJfAQ657o5aAr6XWkN7K8vbMRAUi26mz/MZKepDfd+whI3noMp2Ghoa4BHWLxRItLaA2Ayvd1HKzseeD6J+yWWCOSCu1IFeeoSVJkmolVCpNgUBAcUaKTK6AK98SK+FGIhEIIVQDj8S1ihwOR14DD4DBB1FUc3Nz3IJj2UzDJEonEomkfW/JgYjFYtF9bZlCEAgEogvsSZIEq9UKq9Wa99oT+eDxeBQrJAP/GsqVe2ezlVhGINUq4GZhwikRUR5pTShM7HIvVIFAANddd13KhfTScTgcuP322wtqVpee6uvrFdeY0nsILnGJh1xn3KlhwikRUZFob2/XVKwtduVjSZJQXl5eEL0jXq8XFosFdrsddrsdPp8v48BDrVhbb29vdGG/ysrKgnidevF4PIqLGxqdpG5U4KEVgw8iogIRG4hkmoMUDAajF+jYW+LFWg4Ocr2Ix87EiL01NTVBCIFQKJS2ArDVaoXL5YpWvA2Hw6rLAcjkQKQUghCXy5WU3+F2u3MeXikmHHYhIioCehaYkgOGTNem0bIaszy9MxwO45JLLsl6iCjVsI3FYsG9995blMMxSjNajBwKAeKHXYwctmORMSKiIcCIipe5MGq6sFIgYrPZCmKNIC2UVko2I4fHrMJ2huV8bNiwATNnzsTZZ5+NU089FZMmTcL69euTtnv77bcxc+ZMTJ8+HcceeyymTJmCd955R9urICKilNauXatYOE9pCEOSJLjdbrhcLjgcDk2L4sVu63A44hZEjL0ZlavQ2NiInp6euEX9QqFQUQ2/WCwWxeGoYkgeNkSmC8YIIcS5554r1qxZE/35qaeeEhaLRbz99tvR+95//31RV1cXvS8YDIrvf//7YuPGjRk9BxeWIyIiNUhYCLEYKC0K6Xa7TXt+v9+fdsFMPRi2sNzrr7+OCRMmoKxssDBqd3c3qqur8fjjj+PCCy8EAMyaNQtjx47FrbfeGv29999/H06nEyNHjkz7HBx2ISIiNYlDTRouYXmhVrXU7HabMfRiSs5HKBTCLbfcgj//+c/YsWMHhg8fjlAoFA1GzjvvPMMbT0REQ0/shbTQgw+lirb5aHNiTZF8Bx9ZTbWdN28eDj74YLzwwgvYvHkzhg8fDgD48MMP0dfXh46ODlx44YWYMmUKLrjgAmzbtk11X/39/ejq6oq7ERERKSmmPA8l+QqWEqfxJpZcN1tWwce9996Lr7/+GmeeeSamTJmCzz77DACiVfoWLVqE22+/Ha+++irmzp2LM888E9u3b1fc17Jly1BTUxO91dbWZvlSiIio1M2bNy/6f62LQRYCSZIKoox8LpVn9ZB1kTGr1Yobb7wRQgjccccdgzv7Z0b0pZdeirFjxwIAfvjDH2Ly5Mm4++67FfezcOFCdHZ2Rm979+7NtklERFTiYuuNzJ49O48tyYxST4dcpdZut5veFnldoXyvXaUp+EicU22xWDBu3Djs2rULAKK9FocffnjcdvX19di9e7fiPsvLy1FdXR13IyIiSqdYpqmqVW8NhULRgm9m6enpgRAi74vLaQo+Jk6cmHTfZ599htGjRwMYDDqOOOKI6DCM7PPPP0ddXV0OzSQiIipOcj2WVCXztdRdKQWaXu2uXbvw9NNPR3/+wx/+gPfffx9z5syJ3nf99dfj4Ycfxtdffw1gcHruK6+8gquvvlqnJhMRERWf5uZm1SBECFEQuSBm0TTVduXKlWhqaoLVakU4HIYkSVi0aBEuuOCCuO3uuOMOrFmzBtXV1YhEIli4cCFmzpyZ0XNwqi0REamJHaKoq6sr6sXYXC5XdKKGrNCnDqfCtV2IiKgkmbVOiVlK6fUYXueDiIgoH1LlTVDxYPBBRERFI3HxumKs9SFzuVxxP6vNiilFDD6IiKioJAYcHo8nTy3JXn19fVK+R7FMHdYDgw8iIioqsYXGAKClpaWoyq4HAoG4dVaAodXrATD4ICKiIpTY++Hz+QqmdHkqXq8XPp8v7j6/3z+kej0AznYhIqIilS7fw+l0or293aTWpJe4siwA2Gy2pOrhxYqzXYiIqOSlqxra0dEBSZLyPiQTCARgsViSAg+3210ygYdW7PkgIqKiV1lZmXalVrfbnTRbxmhKvR35aovR2PNBRERDirxgmnxT0tLSEl3IzcjcEI/HE30epcCjrq6u5AIPrRh8EBFRyRFCpFzQVF7WPtXS9rFBhCRJqKysjPtZ7dbS0qK4P7fbDSFEUZeE10tZvhtARERkBPkiny4xVV7aPp10wzqplOIwSy7Y80FERCVNHopxOBymPackSfD7/dHnZuARj8EHERENCYl5IX6/P6PfSywA1tDQEL1PkiQ0NDTE7VcIgUgkgsbGRt1fQ6ngbBciIiLKGWe7EBERUcFi8EFERESmYvBBREREpmLwQURERKYquDofcv5rV1dXnltCREREmZKv25nMYym44KO7uxsAUFtbm+eWEBERkVbd3d2oqalJuU3BTbWNRCLYt28fqqqqMqo4Z5auri7U1tZi7969JT8FmK+1NPG1lp6h8joBvtZiIIRAd3c3Ro8eDYsldVZHwfV8WCwWHH744fluhqrq6uqiejPkgq+1NPG1lp6h8joBvtZCl67HQ8aEUyIiIjIVgw8iIiIyFYOPDJWXl2PJkiUoLy/Pd1MMx9damvhaS89QeZ0AX2upKbiEUyIiIipt7PkgIiIiUzH4ICIiIlMx+CAiIiJTFVydj3x64oknsHTpUgwbNgwWiwX33XcfjjvuONXtX3nlFfzyl79EeXk5+vv7sWLFCpx22mkmtjg769evx+rVqxEOh9HV1YW6ujqsWLECY8eOVdz+oYcewq233opRo0bF3f/MM89g2LBhZjQ5KzfeeCOefPJJjBgxInpfTU0NNmzYoPo7xXpMjznmmKTj88knn2D06NF46aWXkrYvtmMaDAaxZMkSrFixAv/4xz8wZsyYuMfvv/9+3H///Rg2bBhGjBiB3/72tzjssMNS7lPr590saq91YGAADz30EB599FFIkoTOzk5MmDABt956K0aOHKm6v2w+B2ZJdVwvv/xyvPfee6ioqIjeN378eNx///0p91mIxzXV6xwxYgROOumkuO3/8Y9/YPr06Xj44YcV91fIxzRjgoQQQuzYsUMMHz5cvPfee0IIIdasWSMOO+ww0dXVpbj9nj17RHV1tXjxxReFEEJs27ZNVFdXiz179pjV5KzZbDaxefNmIYQQ4XBYzJkzR4wbN04cOHBAcfsHH3xQPPjggya2UB9LliyJHp9MFPMxnTp1atJ9F198sVi1apXi9sV0THfv3i1OOeUUcdlllwkAYvfu3XGP//nPfxaHHHKI+Pzzz4UQQtx0003ipJNOEuFwWHWfWj/vZkn1Wvfu3SsqKirEW2+9JYQQoq+vT5x11lnitNNOS7lPrZ8Ds6Q7rnPmzEm6L51CPK7pXqfSZ3fSpEli06ZNqvss1GOqBYdd/mn58uU4//zzMX78eADAT3/6UwwMDGDNmjWK299zzz045phjMG3aNADA1KlTMX78eKxcudKsJmftRz/6Ec455xwAgxVlr776anz44YfYuXNnnluWX8V8TB988MG4n9vb2/H888/D6/XmqUX6+fbbb/HII4/giiuuUHx86dKlmDNnTvTb//z58/HOO+/gL3/5i+o+tX7ezZLqtdrtdlx55ZU48cQTAQxOx5w7dy5efvll7Nu3z+ym5izdcc1GIR7XdK8z8bP77rvvYt++fZgxY4YZzcsbBh//9MILL8Dj8UR/tlgsmDRpErZs2aK4/ZYtW+K2BwCPx6O6fSF57LHH4n6WuzWDwWA+mlMwivmYHnHEEXE/NzU14bzzzoPT6cxTi/Rz/PHH46ijjlJ8rKOjAzt37ow7bjU1NTj66KNTHjetn3ezpHqtI0eOxL333ht3XzF/dlO91mwV4nFN9zoTP7sPPfQQLrvsMlitVqObllcMPgB8/fXX6OzsTBr/HjVqFD7++GPF3/n44481bV/Itm/fjtGjR2PKlCmq22zatAnTp0/HqaeeilmzZuGNN94wsYXZ+/3vf49p06ZhypQpmDNnDj766CPVbUvpmD700ENpv1EW6zGNJR8bLcctm897odq+fTvcbndSDkwiLZ+DQrJs2TJMmzYNp556KubNm4fPP/9cddtSOK7hcBiPPvooLr/88rTbFusxlTH4ANDb2wsASdXkysvLo48p/Y6W7QuVnFR5zz33wGazKW5zyCGHYNy4cXjmmWfwyiuv4LzzzsPkyZML/mJVV1eHk08+GVu2bMHLL7+MI444ApMmTcKnn36quH2pHNNdu3Zh//79OPvss1W3KdZjmijbz67W3ylEX331FVavXo1Vq1al3E7r56BQHH300Tj99NOxdetWbN26Ff39/TjllFPw7bffKm5fCsd18+bNGDNmDI455piU2xXrMY3F4AOAw+EAMHghjtXf3x99TOl3tGxfqObOnYsf//jHuPjii1W3Oe+887Bs2bLoh/qKK67AhAkTcNttt5nVzKxceeWV+MUvfoGysjJYLBYsXrwYFRUVuO+++xS3L5VjKnfbplrSuliPaaJsP7taf6fQDAwMYPbs2bj55psxefLklNtq/RwUikWLFuEnP/kJLBYL7HY77rjjDrS1taGpqUlx+1I4rpn0WALFe0xjMfgAcNBBB6Gmpgb79++Pu3///v2q00/Hjh2raftCtGDBApSVlWHp0qWaf/fII48sum4+q9WKMWPGqLa7FI6p3G2bTRJfMR5T+dhoOW7ZfN4LSSQSwZw5czB16lTMnTtX8++n+xwUqurqahx88MGq7S7249rR0YEtW7bgkksu0fy7xXhMGXz80/Tp09HS0hL9WQiBnTt34qyzzlLc/swzz4zbHgBaWlpUty80y5cvx549e/Db3/4WkiTh9ddfx+uvv6647cKFC5O6LT/99FPU1taa0dSszZ8/P+m+ffv2qba72I8pADz33HM48sgj0ybyFesxTeR0OnHyySfHHbeuri588MEHKY+b1s97IZk3bx4OO+wwLF68GMBgonSqnAatn4NCkdju/v5+fP311ynbXczHdd26dZg5cyaqq6vTblusxzROnqf6FowdO3aIqqoq8f777wshhHjkkUfi5odffvnl4qc//Wl0e7kmxF//+lchhBAvvfSSqKqqKoqaEH6/Xxx33HHib3/7m2hubhbNzc1iyZIl0boPia916tSp4p577on+/NxzzwmLxSK2bNlidtM1GTNmjNiwYUP059/97neivLxc7Nq1SwhRWsdUNmvWLPH73/8+6f5SOKYvvviiap2PUaNGiS+++EIIIcQtt9ySVOfjjDPOEIsWLYr+nO7znm9qr/X6668XU6dOjX5um5ubxc9+9rO4mg+JrzXd5yDf1F6r3W4Xzc3N0Z9//etfi4MOOihaz0WI4jquaq9T9r3vfU9s3bpV8bFiO6aZYIXTf/re976HNWvWwOv1Rivjbd68GVVVVQCAvr4+hEKh6Pb19fXYtGkT/uu//gt2ux39/f14+umnUV9fn6+XkJHu7m7MmzcPkUgEP/jBD+Iek+ebJ77W66+/HqtWrcJjjz2GcDiMSCSCxx9/HGeeeaapbddq6dKluOuuu3DnnXeiv78fdrsdzz//PL773e8CKJ1jKvvmm2/wwgsv4IEHHkh6rJiPaTAYxDnnnINvvvkGADB79mzU1tZGp4xfdNFF+OKLL3DuueeioqICTqcTGzdujMt5OXDgQFwuQLrPe76keq3vvvsuli9fDgBJU8Jj67kkvtZ0n4N8SXdcb7vttmheQ29vL77zne/gxRdfjKvmWgzHNd3rBID33nsPX375ZbTGUKJiOaZaSEIIke9GEBER0dDBnA8iIiIyFYMPIiIiMhWDDyIiIjIVgw8iIiIyFYMPIiIiMhWDDyIiIjIVgw8iIiIyFYMPIiIiMhWDDyIiIjIVgw8iIiIyFYMPIiIiMhWDDyIiIjLV/w/2tPA5xQx3zgAAAABJRU5ErkJggg==",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1,1)\n",
"ax.scatter(dr.longitude, dr.latitude, c=\"k\", s=1)\n",
"ax.scatter(dr_left.longitude, dr_left.latitude, c=\"C1\", s=2)\n",
"ax.set_aspect(np.cos(np.deg2rad(40)))"
]
},
{
"cell_type": "markdown",
"id": "7685f9c0-c009-49fb-b793-167c2eb05c15",
"metadata": {},
"source": [
"### play with geopandas"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0580ae00-b366-47ce-8781-bdde8d7eac65",
"metadata": {},
"outputs": [],
"source": [
"gdf = gpd.GeoDataFrame(dict(i=range(3)), geometry=[both, left, right],index=[\"both\", \"left\", \"right\"])"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f6cd90e9-8800-44c4-83e5-6e814c444164",
"metadata": {},
"outputs": [],
"source": [
"gdf.loc[[\"both\"]].plot()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2130cdcb-73a5-4b62-9c6c-adfa8d2ba917",
"metadata": {},
"outputs": [],
"source": [
"gdf.plot(column=\"i\", alpha=.8)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d162fb95-5573-4b81-aad6-25d641a69bf6",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.13"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment