Skip to content

Instantly share code, notes, and snippets.

@iamgeoknight

iamgeoknight/get_coords.ipynb

Last active January 20, 2024 14:23
Show Gist options
  • Save iamgeoknight/4829aa7a113405a17659f60a15cbc256 to your computer and use it in GitHub Desktop.
Save iamgeoknight/4829aa7a113405a17659f60a15cbc256 to your computer and use it in GitHub Desktop.
Download ZIP
Get the geographical coordinates from NetCDF file using Python
Display the source blob
Display the rendered blob
Raw
get_coords.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Import Libraries"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [],
"source": [
"import xarray as xr\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"import geopandas as gpd\n",
"import math"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Use xarray to read netcdf file. It is recommended to use decode_coords=\"all\" for compatibility with RioXarray. You can rioxarray functionalities like reprojecting, getting CRS, extent etc information "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"ds = xr.open_dataset(\"HLSTimeSeries.nc\", decode_coords=\"all\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Check CRS of netcdf"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"CRS.from_epsg(32756)"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ds.rio.crs"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(217305.0, 6765165.0)"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"float(ds.x[0]), float(ds.y[0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Reprojecting netcdf to 4326"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"ds_4326 = ds.rio.reproject(\"EPSG:4326\")"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(150.09192486031887, -29.210732831128635)"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"float(ds_4326.x[0]), float(ds_4326.y[0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Get XY coordinates which are in 32756 projection. If you want Geographic coordnates, then you can reproject the xarray ataset to 4326"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [],
"source": [
"x, y = np.meshgrid(ds.x, ds.y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Get all variables"
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"['QA', 'Blue', 'Green', 'Red', 'NIR_Narrow', 'SWIR1', 'SWIR2']"
]
},
"execution_count": 78,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"variables = list(ds.var())\n",
"variables"
]
},
{
"cell_type": "code",
"execution_count": 79,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.QuadMesh at 0x7fe450e22f50>"
]
},
"execution_count": 79,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAHUCAYAAAAukLw0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB4vElEQVR4nO3de1xUZf4H8M8Ml+EOSnJTAjItFFFTQzR/6ualUsvcaCu8ppspJWpbRraGruKlNOxmi6FoZbZquuwmRq1KFzVv2Zq13i+EIKXIRZQB5vn94TLbeJ4jc5yBAefzfr3m9ZJnnvOc55w5jA/Pc77foxNCCBARERGRBb2jO0BERETUFHGQRERERCTBQRIRERGRBAdJRERERBIcJBERERFJcJBEREREJMFBEhEREZEEB0lEREREEhwkEREREUlwkNRE7NixA6mpqbh48aLivX79+qFfv36N3qfGcOTIEfzpT39Ct27dEBAQgJYtW6J3795Yv369tH5xcTHGjh2LW265BV5eXoiPj8e//vUvizplZWWYN28e+vXrh5CQEPj4+KBTp05YuHAhrly5YlF33759SEpKQqdOneDr64vg4GAMGDAAW7du1XQc1vSrzqVLlzBr1iy0b98eBoMBgYGB6N+/P44ePVrvfrScr59//hlTp05F3759ERAQAJ1Oh6ysLE3Hda2srCzodDrs3bvXpnbsTafTITU1tVH3uXr1ajz22GO44447oNfrERkZ2aj7t5ezZ88iNTUVBw4cULyXmpoKnU53Q+3asq2tXn75Zdx6661wdXVFQEBAo+23oqICU6dORVhYGDw8PNClSxesXbu20fZP9sdBUhOxY8cOzJ49WzpIeuedd/DOO+80fqcaQW5uLj799FP8/ve/x7p16/Dhhx+iXbt2SEhIwJw5cyzqVlVV4d5778W//vUvLF26FH//+98RHByM++67D3l5eeZ6Z86cQXp6Ou666y5kZGQgOzsbjzzyCFJTUzF06FD89kk8H330EXbv3o0nn3wSf//73/Hee+/BYDDg3nvvxerVq606Bmv7BVz9Eu3Xrx8yMzPx7LPPIjc3FytXrkRcXBwqKyvter6OHTuGDz/8EO7u7njggQesOhay3vvvv49Dhw7h7rvvRtu2bR3dnRt29uxZzJ49WzpImjBhAnbu3Nn4nbLB3//+d8ybNw+jR49GXl4evvjii0bb94gRI7Bq1Sq88soryMnJQY8ePfD4449jzZo1jdYHsjNBTcKrr74qAIiTJ086uiuN6pdffhEmk0lRPmTIEOHl5SWuXLliLnv77bcFALFjxw5zWXV1tejQoYO4++67zWUVFRWioqJC0WbdOf7qq6/MZefOnVPUq6mpEbGxsaJt27ZWHYO1/RJCiOTkZOHt7S2OHz9uVdvX0nK+amtrzf/es2ePACBWrlx5Q/uts3LlSgFA7Nmzx6Z27A2AeOWVVxp1n789v0OGDBERERGNun97sde1ca1XXnlF2OO/mEuXLmmqP3fuXAFA+rvdkD799FMBQKxZs8aifODAgSIsLEzU1NQ0an/IPjiT1ASkpqbi+eefBwBERUVBp9NBp9Nh+/btAJTLbadOnYJOp8Orr76KhQsXIjIyEp6enujXrx+OHDmC6upqvPjiiwgLC4O/vz8efvhhFBcXK/b78ccfIz4+Ht7e3vDx8cHgwYPx3XffNcYhm91yyy3SKfm7774blZWVuHDhgrls48aNuOOOOxAfH28uc3V1xciRI7F7924UFBQAALy9veHt7S1tEwDy8/PNZUFBQYp6Li4u6Natm0W967G2X5WVlXjvvfeQkJCA2267zaq2r6XlfOn1DffrXVJSgnHjxqFly5bw9vbGsGHDcOLECYs6n3/+OR566CG0adMGHh4euP322zFx4kT8+uuvFvXqlmUOHTqExx9/HP7+/ggODsaTTz6J0tJSi7plZWX44x//iMDAQPj4+OC+++7DkSNHGuw4r6chzu/WrVvRr18/BAYGwtPTE7feeit+//vfm2cZ6373Fy1ahHnz5uHWW2+Fh4cHunfvrljePXbsGMaNG4d27drBy8sLrVu3xrBhw3Dw4EFzne3bt6NHjx4AgHHjxpm/e+qWLmVLZh9//DEGDRqE0NBQeHp6Ijo6Gi+++CIuXbpk8/HX7W///v145JFH0KJFC/MsnRAC77zzDrp06QJPT0+0aNECjzzyiMV1FxkZiZdffhkAEBwc3KjLsBs3boSPjw8SEhIsyseNG4ezZ8/i22+/bZR+kH1xkNQETJgwAc8++ywA4JNPPsHOnTuxc+dO3HXXXdfd7u2338Y333yDt99+G++99x7+85//YNiwYRg/fjx++eUXrFixAosWLcIXX3yBCRMmWGyblpaGxx9/HB06dMDf/vY3vP/++ygvL0efPn3w448/1tvnmpoaq17iN0tbWmzbtg2tWrWyGMT88MMPiI2NVdStKzt06NB126y7z6hjx47XrVdTU4Ovvvqq3npa+7Vv3z5cunQJ7dq1w6RJk9CiRQu4u7uje/fu+PTTT63alxrZ+WpI48ePh16vx5o1a5Ceno7du3ejX79+FsvFx48fR3x8PJYtW4bc3FzMmjUL3377Le655x5UV1cr2vz973+P9u3bY8OGDXjxxRexZs0aTJs2zfy+EALDhw/H+++/j+eeew4bN25Ez549cf/991vd74a+bm1x6tQpDBkyBO7u7lixYgW2bNmCBQsWwNvbG0aj0aLuW2+9hS1btiA9PR0ffPAB9Ho97r//foulsbNnzyIwMBALFizAli1b8Pbbb8PV1RVxcXE4fPgwAOCuu+7CypUrAVy9j6fuu+fa74vfOnr0KB544AFkZmZiy5YtmDp1Kv72t79h2LBhdjsXI0aMwO23345169bh3XffBQBMnDgRU6dOxYABA7Bp0ya88847OHToEHr16oVz584BuDpQGT9+PABgy5Yt9R6LEMLqa6I+P/zwA6Kjo+Hq6mpRXvc98MMPP9zQuSAHc+g8Fpldb7mtb9++om/fvuafT548KQCIzp07W0z5p6enCwDiwQcftNh+6tSpAoAoLS0VQghx5swZ4erqKp599lmLeuXl5SIkJEQ8+uij1+1r3f6teW3btk3biRBCLF++XAAQS5cutSh3c3MTEydOVNTfsWOHdJr7t77//nvh6ekpHn744Xr3P3PmTAFAbNq0yar+Wtuvjz76SAAQfn5+onfv3iI7O1v885//FP379xc6nU5s2bLFqv1dS+18/Za9l9uuPY/ffPONACDmzp0r3c5kMonq6mpx+vRpAUD8/e9/N79XtyyzaNEii20mT54sPDw8zMuLOTk50uOcN2+e1ctt1l63Ws+TPZbb1q9fLwCIAwcOqNap+90LCwsTly9fNpeXlZWJli1bigEDBqhuW1NTI4xGo2jXrp2YNm2aufx610Z9S2Z1n2teXp4AIL7//nurt5Wp22bWrFkW5Tt37hQAxOLFiy3K8/Pzhaenp3jhhRcUbfzyyy/17q/uerbmVZ927dqJwYMHK8rPnj0rAIi0tLR626Cmx3LIS83KAw88YDHlHx0dDQAYMmSIRb268jNnziAmJgafffYZampqMHr0aIu/kDw8PNC3b19s27btuvsNCwvDnj17rOrjHXfcYVW9Ojk5OUhKSsIjjzxinl37retFy6i9d+rUKQwdOhTh4eF47733rrv/9957D/PmzcNzzz2Hhx56yFwuhEBtba1F3d/+xWhNv0wmEwDA3d0dOTk58PX1BQD0798f7dq1w1/+8hcMHjwYABR/ubq4uEj3Ud/5aiiJiYkWP/fq1QsRERHYtm0bZs6cCeBqxN+sWbPw6aef4uzZs+bjB4CffvoJDz74oEUb1/4cGxuLK1euoLi4GMHBwebr8tp9P/HEE+Z91sfa6zYqKsqqevbUpUsXuLu746mnnsLkyZPRp08f1WXZESNGwMPDw/yzr68vhg0bho8++gi1tbVwcXFBTU0NFi1ahA8++ADHjh2zmL376aefbrifJ06cwMsvv4ytW7eiuLjYYtbtp59+ks6qavX73//e4ud//vOf0Ol0GDlypMXvRkhICDp37my+NUGrYcOGWX1NWONGvp+oaeMgqRlr2bKlxc/u7u7XLa8Lf6+bmq67F+Fa9d1r4e7uji5duljVRxcXF6vqAcBnn32GESNGYODAgfjwww8VXyqBgYE4f/68Yru6+3CuPW4AOH36NPr37w9XV1f861//ktaps3LlSkycOBFPPfUUXn31VYv3Vq1ahXHjxlmU1f3nYG2/AgMDAVwdUNQNkADAy8sLffv2xaZNm8xlbm5uir6NHTvWoqy+89WQQkJCpGV158FkMmHQoEE4e/Ys/vznP6NTp07w9vaGyWRCz549cfnyZcX2deenjsFgAABz3fPnz8PV1VVRT9YXNQ1x3dpL27Zt8cUXX2DRokVISkrCpUuXcNttt2HKlClITk62qKt2/o1GIyoqKuDv74/p06fj7bffxowZM9C3b1+0aNECer0eEyZMkJ5/a1RUVKBPnz7w8PDA3Llz0b59e3h5eSE/Px8jRoy44XavFRoaavHzuXPnIIRAcHCwtP6N3uPXsmVL+Pv739C217qR7ydq+jhIckK33HILAGD9+vWIiIjQvP2pU6es/kt727ZtVuV4+uyzzzB8+HD07dsXGzZsMA/sfqtTp04WN53WqSuLiYmxKD99+jT69esHIQS2b9+ONm3aqO5/5cqVmDBhAsaMGYN3331XMeC43l+c1vbren9hCyEsBqfX7uva823N+WpIRUVF0rLbb78dwNX7L77//ntkZWVhzJgx5jrHjh274X0GBgaipqYG58+ftxgoyfqi5trBpxrZoLQx9OnTB3369EFtbS327t2LN998E1OnTkVwcDAee+wxcz218+/u7g4fHx8AwAcffIDRo0cjLS3Not6vv/56w7mDtm7dirNnz2L79u3o27evuVyWusQW1/7+1QUsfPXVV+bB82/Jyqwh++NHjajnPrVOnTrho48+Qk1NjcUss9r3EzUPHCQ1Edf+1dyQBg8eDFdXVxw/flwxrW0Ney+35ebmYvjw4bjnnnuwadMm1S+8hx9+GJMnT8a3336LuLg4AFeXpT744APExcUhLCzMXPfMmTPo168famtrsX379usOBrOysjBhwgSMHDkS7733nnRGJjAwUDGDobVfoaGhiI+PxzfffIOysjL4+fkBuBr1lpeXh549e5rb7N69u83nqyF9+OGHFtfOjh07cPr0afNNsnXn8Nq+/fWvf73hffbv3x+LFi3Chx9+iClTppjLteSgacrLbb/l4uKCuLg43Hnnnfjwww+xf/9+i0HSJ598gldffdW85FZeXo5//OMf6NOnj3kWTKfTKc7/p59+ioKCAvNgFtD23dMQn6s1hg4digULFqCgoACPPvqo3dq153Lbww8/jOXLl2PDhg34wx/+YC5ftWoVwsLCzN8N1LxwkNREdOrUCQCwdOlSjBkzBm5ubrjjjjsslmXsJTIyEnPmzMHMmTNx4sQJ3HfffWjRogXOnTuH3bt3w9vbG7Nnz1bdvi4iyx6+/vprDB8+HCEhIXjppZcUCe06dOhgHkw8+eSTePvtt5GQkIAFCxYgKCgI77zzDg4fPmyRMK64uBj9+/dHYWEhMjMzUVxcbJECoU2bNuZZpXXr1mH8+PHo0qULJk6ciN27d1vsv2vXrvUOQqztFwC89tpr6N+/PwYPHowZM2ZAp9Nh8eLF+PXXX/GXv/zFrucLgDkTd12Y9N69e80zDY888oi5XmpqKmbPnm31zN/evXsxYcIEJCQkID8/HzNnzkTr1q0xefJkAMCdd96Jtm3b4sUXX4QQAi1btsQ//vEPfP755/W2rWbQoEH4v//7P7zwwgu4dOkSunfvjm+++Qbvv/++1W3Y67oFgB9//NEcCVpUVITKykrz+e7QoQM6dOhgrqvT6dC3b9/r3jvz7rvvYuvWrRgyZAhuvfVWXLlyBStWrAAADBgwwKKui4sLBg4ciOnTp8NkMmHhwoUoKyuz+L0dOnQosrKycOeddyI2Nhb79u3Dq6++qphRbdu2LTw9PfHhhx8iOjoaPj4+CAsLs/ijo06vXr3QokULPP3003jllVfg5uaGDz/8EN9//722k6dR79698dRTT2HcuHHYu3cv/u///g/e3t4oLCzE119/jU6dOmHSpEma273eHz9a3X///Rg4cCAmTZqEsrIy3H777fjoo4+wZcsWfPDBBw5ZwiU7cNw943StlJQUERYWJvR6vUVkmFp026uvvmqx/bZt2wQAsW7dOotytQSAmzZtEv379xd+fn7CYDCIiIgI8cgjj4gvvviiQY5Ppi4SRe11bXRcUVGRGD16tGjZsqXw8PAQPXv2FJ9//rlFnbrzoPb6bRTUmDFjrlvX2uSe1vSrzldffSX69u0rvLy8hJeXl/jd734nvvnmmwY5X9er+1vPPfec0Ol04qeffrru/uuupdzcXDFq1CgREBAgPD09xQMPPCCOHj1qUffHH38UAwcOFL6+vqJFixYiISFBnDlzRvEZqEUj1e3rt5/BxYsXxZNPPikCAgKEl5eXGDhwoPjPf/7jkGSS1/ssftuX8vJyAUA89thj121v586d4uGHHxYRERHCYDCIwMBA0bdvX5GdnW2uU/e7v3DhQjF79mzRpk0b4e7uLrp27So+++wzi/ZKSkrE+PHjRVBQkPDy8hL33HOP+dr77feJEFcjL++8807h5uZm0X9ZhNqOHTtEfHy88PLyEq1atRITJkwQ+/fvV0TI2RLdphaZtmLFChEXFye8vb2Fp6enaNu2rRg9erTYu3ev1W00pPLycjFlyhQREhIi3N3dRWxsrPjoo48avR9kPzohHJAQhIialLvvvhsRERFYt26do7ty09m8eTOGDh2K77//3jxjfKPq7gd89dVX8ac//clOPSQiNVxuI3JyZWVl+P7777Fq1SpHd+WmtG3bNjz22GM2D5CIqPFxkETk5Pz8/FBVVeXobty0rk0n4YxMJpNFniyZazNVEzUFXG4jIqIGVRcYcD0nT55EZGRk43SIyEocJBERUYM6e/Yszp49e906sbGxjZ7vi6g+HCQRERERSVz/+RNEREREToqDJCIiIiIJDpKIiIiIJDhIaiRffvklhg0bhrCwMOh0OosnvltLCIHXXnsN7du3h8FgQHh4uOLhlURERGQfTEzRSC5duoTOnTtj3LhxN/RQWQBITk5Gbm4uXnvtNXTq1AmlpaX49ddf7dxTIiIiAhjd5hA6nQ4bN27E8OHDzWVGoxEvv/wyPvzwQ1y8eBExMTFYuHCh+WGjP/30E2JjY/HDDz/gjjvucEzHiYiInAiX25qIcePG4ZtvvsHatWvx73//GwkJCbjvvvtw9OhRAMA//vEP3HbbbfjnP/+JqKgoREZGYsKECbhw4YKDe05ERHRz4iCpCTh+/Dg++ugjrFu3Dn369EHbtm3xpz/9Cffccw9WrlwJADhx4gROnz6NdevWYfXq1cjKysK+ffvwyCOPOLj3RERENyfek9QE7N+/H0IItG/f3qK8qqoKgYGBAK4++6iqqgqrV68218vMzES3bt1w+PBhLsERERHZGQdJTYDJZIKLiwv27dsHFxcXi/d8fHwAAKGhoXB1dbUYSEVHRwMAzpw5w0ESERGRnXGQ1AR07doVtbW1KC4uRp8+faR1evfujZqaGhw/fhxt27YFABw5cgQAEBER0Wh9JSIichaMbmskFRUVOHbsGICrg6IlS5agf//+aNmyJW699VaMHDkS33zzDRYvXoyuXbvi119/xdatW9GpUyc88MADMJlM6NGjB3x8fJCeng6TyYSkpCT4+fkhNzfXwUdHRER08+EgqZFs374d/fv3V5SPGTMGWVlZqK6uxty5c7F69WoUFBQgMDAQ8fHxmD17Njp16gTg6pO0n332WeTm5sLb2xv3338/Fi9ejJYtWzb24RAREd30OEgiIiIikmAKACIiIiIJDpKIiIiIJBjd1oBMJhPOnj0LX19f6HQ6R3eHiIiaMCEEysvLERYWBr2+4eYwrly5AqPRaHM77u7u8PDwsEOPmi4OkhrQ2bNnER4e7uhuEBFRM5Kfn482bdo0SNtXrlxBVIQPioprbW4rJCQEJ0+evKkHShwkNSBfX18AVy94Pz8/h/bliV3PScvX9FzcyD0hIiKZsrIyhIeHm//vaAhGoxFFxbU4vS8Sfr43PltVVm5CRLdTMBqNHCTRjalbYvPz83P4IMnN211a7uh+ERGRpca4PcPHVwcf3xvfjwnOcQsJB0lEREROplaYUGtDAqBaYbJfZ5owDpKIiIicjAkCJtz4KMmWbZsTpgAgIiIikuBMEhERkZMxwQRbFsxs27r54CDJSXi5ynNiPLV3jKIso/uqhu4OERE5UK0QqLXhqWS2bNuccLmNiIiISIIzSURERE6GN25bh4MkIiIiJ2OCQC0HSfXichsRERGRBGeSiIiInAyX26zDQRIREZGTYXSbdThIusn8+eDD0vL9Z7tIyx+5/UDDdYaIiKgZ4yCJiIjIyZj++7Jle2fAQRIREZGTqbUxus2WbZsTDpKIiIicTK24+rJle2fAFABEREREEpxJIiIicjK8J8k6HCQRERE5GRN0qIXOpu2dAQdJN5nPCqKl5TU1LtLyvOLbFWWJ3/5RWvfDuOU33jEiIqJmhoMkIiIiJ2MSV1+2bO8MOEgiIiJyMrU2LrfZsm1zwug2IiIiIgnOJBERETkZziRZh4MkIiIiJ2MSOpiEDdFtNmzbnHC5jYiIiEiCM0lEREROhstt1uEg6SZzriBAWm4ocpOWny/zVpT9HBIqbzzuRntFRERNSS30qLVhManWjn1pyjhIIiIicjLCxnuSBO9JIiIiInJenEkiIiJyMrwnyTocJBERETmZWqFHrbDhniQneSwJl9uIiIiIJDiTRERE5GRM0MFkwzyJCc4xlcRBkpNwvSwv9y5SXuj6audYayYicla8J8k6XG4jIiIikuBMEhERkZOx/cZt51huc/hMUkFBAUaOHInAwEB4eXmhS5cu2Ldvn2r9sWPHQqfTKV4dO3a0qHfx4kUkJSUhNDQUHh4eiI6OxubNm83vp6amKtoICQmxaEMIgdTUVISFhcHT0xP9+vXDoUOH7HsCiIiIGtnVe5Jse2kRGRkp/b87KSkJQNP9/9ahg6SSkhL07t0bbm5uyMnJwY8//ojFixcjICBAdZulS5eisLDQ/MrPz0fLli2RkJBgrmM0GjFw4ECcOnUK69evx+HDh7F8+XK0bt3aoq2OHTtatHXw4EGL9xctWoQlS5bgrbfewp49exASEoKBAweivLzcrueBiIjoZrZnzx6L/28///xzADD/391U/7916HLbwoULER4ejpUrV5rLIiMjr7uNv78//P39zT9v2rQJJSUlGDdunLlsxYoVuHDhAnbs2AE3t6vPLIuIiFC05erqqpg9qiOEQHp6OmbOnIkRI0YAAFatWoXg4GCsWbMGEydOtPo4iYiImhKTjc9u0xrd1qpVK4ufFyxYgLZt26Jv375N+v9bh84kZWdno3v37khISEBQUBC6du2K5cuXa2ojMzMTAwYMsBgEZWdnIz4+HklJSQgODkZMTAzS0tJQW2v5SL6jR48iLCwMUVFReOyxx3DixAnzeydPnkRRUREGDRpkLjMYDOjbty927Ngh7UtVVRXKysosXkRERE1N3T1JtrxulNFoxAcffIAnn3wSOp3uhv6/bSwOnUk6ceIEli1bhunTp+Oll17C7t27MWXKFBgMBowePbre7QsLC5GTk4M1a9Yo2t26dSsSExOxefNmHD16FElJSaipqcGsWbMAAHFxcVi9ejXat2+Pc+fOYe7cuejVqxcOHTqEwMBAFBUVAQCCg4Mt2g4ODsbp06el/Zk/fz5mz559I6fCbnRu8tG9T7683Ku4RlHmcVG+1tx5yuvS8u/fmGZl74iIqCkwQW+XPEnXTgYYDAYYDIbrbrtp0yZcvHgRY8eOBYAb+v+2sTh0JslkMuGuu+5CWloaunbtiokTJ+KPf/wjli1bZtX2WVlZCAgIwPDhwxXtBgUFISMjA926dcNjjz2GmTNnWrR7//334/e//z06deqEAQMG4NNPPwVwdYrvt3Q6ywGDEEJRViclJQWlpaXmV35+vlXHQURE1ByFh4ebb4Px9/fH/Pnz690mMzMT999/P8LCwizKtfx/21gcOpMUGhqKDh06WJRFR0djw4YN9W4rhMCKFSswatQouLu7K9p1c3ODi4uLRbtFRUUwGo2K+gDg7e2NTp064ejRowBgvlepqKgIoaGh5nrFxcWK0W4da0bQREREjlYrdKgVNiST/O+2+fn58PPzM5fX93/g6dOn8cUXX+CTTz4xl93I/7eNxaEzSb1798bhw4ctyo4cOSK9yfpaeXl5OHbsGMaPHy9t99ixYzCZTBbthoaGSgdIwNX7iX766SfzBxQVFYWQkBDzHfjA1XXUvLw89OrVy6rjIyIiaopq/3vjti0vAPDz87N41TdIWrlyJYKCgjBkyBBzWVP+/9ahg6Rp06Zh165dSEtLw7Fjx7BmzRpkZGSY8yYAV5ewZPcnZWZmIi4uDjExMYr3Jk2ahPPnzyM5ORlHjhzBp59+irS0NIt2//SnPyEvLw8nT57Et99+i0ceeQRlZWUYM2YMgKvTflOnTkVaWho2btyIH374AWPHjoWXlxeeeOKJBjgbRERENy+TyYSVK1dizJgxcHX930JWU/7/1qHLbT169MDGjRuRkpKCOXPmICoqCunp6UhMTDTXKSwsxJkzZyy2Ky0txYYNG7B06VJpu+Hh4cjNzcW0adMQGxuL1q1bIzk5GTNmzDDX+fnnn/H444/j119/RatWrdCzZ0/s2rXLYhbrhRdewOXLlzF58mSUlJQgLi4Oubm58PX1tfOZICIiajwmoYfJhgg10w1k3P7iiy9w5swZPPnkk4r3mur/tzohnCS3uAOUlZXB398fpaWlFmu2DSly9UJpedC/3KTlsug2k7t8nbo0Qj6mZnQbEZHtGuP/jLp9LN/fDV6+LvVvoKKyvBZ/vGtfo/7/5gh8dttNRlcq/0hlgyEAcC+5oii7Euwpb+NXk7SciIjoZsRBEhERkZMxATZFtznLn8wcJBERETkZ25NJOjTuq9E4x1ESERERacSZJCIiIidj6/PXbNm2OeEgiYiIyMmYoIMJttyT5NjHhTQWDpKIiIicDGeSrMNB0k3GrUx+4Zrc5bEIxhYeirJaN/lfCGrl/e5fpCjbnvOCWheJiIiaBQ6SiIiInMxvn792o9s7Aw6SiIiInIxJ6GCyJU+SDds2J84xFCQiIiLSiDNJRERETsZk43KbsyST5CCJiIjIyZiEHiYbItRs2bY5cY6jJCIiItKIM0k3GWNgrbT8107yj9qlUlnmd0behku1kJZXtFa23XHG69K6hxZOk5YTEVHjqYUOtTYkhLRl2+aEgyQiIiInw+U26zjHURIRERFpxJkkIiIiJ1ML25bM5Ddl3Hw4SCIiInIyXG6zDgdJREREToYPuLWOcxwlERERkUacSWrGIjIXKcq8f3aT1jXJiwEvZVFFmIu0qs9Z+Sq0yV851q72U9kfERE5nIAOJhvuSRJMAUBEREQ3Iy63Wcc5jpKIiIhII84kERERORmT0MEkbnzJzJZtmxMOkoiIiJxMLfSotWExyZZtmxPnOEoiIiIijTiTRERE5GS43GYdDpKaMX2lMlTfqDX03lNZdCVQSKvWeMpTA7heVpbpa+S7u3P269Ly/7wyTb4B2VXnf/5ZWv790L80ck+IyJFM0MNkw2KSLds2J85xlEREREQacSaJiIjIydQKHWptWDKzZdvmhIMkIiIiJ8N7kqzDQRIREZGTEUIPkw1ZswUzbhMRERE5L84kEREROZla6FBrw0Nqbdm2OeEgqRnz+lk5ESjc5HVrJKH+AFDjpQz3r/WTx+9XuspTABguKPtxJcoorevfskJanvjtHxVlt3n9Kq37l04bpeU3u3br5krLQ1uWSssDPSoVZf6eXtK6aqkBDK61irLd96WpdZHIKURkLpKW66qV34X6y/IFG321cpBhunLFto5pYBK23VdkkmeKuek4fLmtoKAAI0eORGBgILy8vNClSxfs27dPtf7YsWOh0+kUr44dO1rUu3jxIpKSkhAaGgoPDw9ER0dj8+bN0jbnz58PnU6HqVOnWpRXVFTgmWeeQZs2beDp6Yno6GgsW7bM5mMmIiKips+hM0klJSXo3bs3+vfvj5ycHAQFBeH48eMICAhQ3Wbp0qVYsGCB+eeamhp07twZCQkJ5jKj0YiBAwciKCgI69evR5s2bZCfnw9fX19Fe3v27EFGRgZiY2MV702bNg3btm3DBx98gMjISOTm5mLy5MkICwvDQw89ZNvBExEROYjJxhu3bdm2OXHoIGnhwoUIDw/HypUrzWWRkZHX3cbf3x/+/v7mnzdt2oSSkhKMGzfOXLZixQpcuHABO3bsgJvb1fWniIgIRVsVFRVITEzE8uXLMXeucilj586dGDNmDPr16wcAeOqpp/DXv/4Ve/fu5SCJiIiaLRN0MNlwX5Et2zYnDh0KZmdno3v37khISEBQUBC6du2K5cuXa2ojMzMTAwYMsBgEZWdnIz4+HklJSQgODkZMTAzS0tJQW2t5f0VSUhKGDBmCAQMGSNu+5557kJ2djYKCAgghsG3bNhw5cgSDBw+W1q+qqkJZWZnFi4iIiJonh84knThxAsuWLcP06dPx0ksvYffu3ZgyZQoMBgNGjx5d7/aFhYXIycnBmjVrFO1u3boViYmJ2Lx5M44ePYqkpCTU1NRg1qxZAIC1a9di//792LNnj2r7b7zxBv74xz+iTZs2cHV1hV6vx3vvvYd77rlHWn/+/PmYPXu2hjNARETU+Jhx2zoOHSSZTCZ0794daWlXo2W6du2KQ4cOYdmyZVYNkrKyshAQEIDhw4cr2g0KCkJGRgZcXFzQrVs3nD17Fq+++ipmzZqF/Px8JCcnIzc3Fx4eHqrtv/HGG9i1axeys7MRERGBL7/8EpMnT0ZoaKh09iklJQXTp083/1xWVobw8HArzwYREVHj4D1J1nHoICk0NBQdOnSwKIuOjsaGDRvq3VYIgRUrVmDUqFFwd3dXtOvm5gYXl/+FrEdHR6OoqAhGoxH79u1DcXExunXrZn6/trYWX375Jd566y1UVVXBaDTipZdewsaNGzFkyBAAQGxsLA4cOIDXXntNOkgyGAwwGAyazgERERE1TQ4dJPXu3RuHDx+2KDty5Ij0Jutr5eXl4dixYxg/fry03TVr1sBkMkGv15vbDQ0Nhbu7O+69914cPHjQYptx48bhzjvvxIwZM+Di4oLq6mpUV1ebt6/j4uICk8mk9VAbhL5aWeZ3XJnXBgAMJfLcRzJVLeSXxeVW8r8cqiU5mK5UynMqXaxWRhgCwO7LysHlo93kS6GmovbScn3IEWl5Q5HldgKAoyWtFGW/HA2U1hVu8mQjnsGXFGUtfKs09E4u3OeitLxbYLm03MdFuc+Hv0mS1t3Y++0b7ldzcTg/TFF2R/hZB/SE7O22Ndbn/9JXyv8YditTfke6qdyaWu2nLBOS7/SGYoKNz25zkhu3HTpImjZtGnr16oW0tDQ8+uij2L17NzIyMpCRkWGuk5KSgoKCAqxevdpi28zMTMTFxSEmJkbR7qRJk/Dmm28iOTkZzz77LI4ePYq0tDRMmTIFAODr66vYztvbG4GBgeZyPz8/9O3bF88//zw8PT0RERGBvLw8rF69GkuWLLH3qSAiImo0wsboNsFBUsPr0aMHNm7ciJSUFMyZMwdRUVFIT09HYmKiuU5hYSHOnDljsV1paSk2bNiApUuXStsNDw9Hbm4upk2bhtjYWLRu3RrJycmYMWOGpv6tXbsWKSkpSExMxIULFxAREYF58+bh6aef1n6wRERETYRJ2DiTxBu3G8fQoUMxdOhQ1fezsrIUZf7+/qisVD5y4bfi4+Oxa9cuq/uxfft2RVlISIhFDiciIiJyHg4fJBEREVHjYnSbdThIIiIicjJcbrOOcwwFiYiIiDTiTNJNptpbPu51qZaH5OuN8hB0Gc9f5OkFAiTpBYL3yusaW7hLy6u9lXkEZm8dI637XJS8f27tlNnO7wwqltY9f8VLWv5LmY+irKZGfu68PEOk5Ym37VWULS/vLa3r/5m8H1V+yvjgi/IsAqj2k6ekyJeUuYZcltb18pSnFzC4Kj9HDzd5nHLnf/5ZWv790L9IyxvbU3uV19MvVcrPGwAqa9yk5Xe3jFMWXnxYWvcvnTZa37kmQnaOACD3xw7Scpmubc9Iy71cjdLybw7frijTn5eff321fPZClg6lxkv+3abWhkeZ8jtIpzEkXxbWX6Xye2uSpP8w6a3/PrYVn91mHQ6SiIiInAyX26zD5TYiIiIiCc4kERERORnOJFmHgyQiIiInw0GSdbjcRkRERCTBmSQiIiInw5kk63CQ1IxdlkSgG1SeOF3rJr+g1cq1qGqhvIxM7vIQXjWGEmWouffPV6R1Wx6St3E5yFtRdiJMGV4MAG6V8lBbX0nIr0nlUMqilPsDANymLHrszn3Squtdu0jLXb9SxhKHfS2PR1ZL+2ByV362FSHyPld7ycurJMdeoRIWrXaeOm19XVFWK898gMr28jDxOyIKFWWXVcL0r1TLy6tqblWUXfqxpbRurac8rUJRlK+0XObDg/Knyp944iWr27CHRT/eLy3/Z6Hy4eCnT3WU1nUpk/9XIQun/65SJUeHCq+flSk2vIrUQuHl5VV+yn4YzqukC1BmLLEbV0mGDa9fVa6leGWZcGm8FAACtoXxN15PHYuDJCIiIifDmSTr8J4kIiIiIgnOJBERETkZziRZh4MkIiIiJ8NBknW43EZEREQNrqCgACNHjkRgYCC8vLzQpUsX7Nv3v6AWIQRSU1MRFhYGT09P9OvXD4cOqUTqNBIOkoiIiJxM3UySLS8tSkpK0Lt3b7i5uSEnJwc//vgjFi9ejICAAHOdRYsWYcmSJXjrrbewZ88ehISEYODAgSgvL7fz0VuPy203GdnTsAGgKkA+HlYLhZe2bZTXlYWaq6UWcKmWtyELY6/2Nmhqw1CijO2VlV2PyV3ZD9nxAYC+Whm6DABv7+yvKIuI/EVa95HbD0jL3y/vqSjzOyPfn2dxlbTc5KasrzfKf+XVrg9ZaLXnefn59yqWn2vZZ1sRJj8W93x3afl/qsMUZTrJU9QBAJXytoWbMhTbWyVlBsrk56O8poWizKSSLkC2PwCIXL1QUebfskJa97YWF6Tl568ocyicKbhFWheV90qL3c8rz5OPyvkQKukdjJKn3nuck59/ryJ5G66S7yC3S/Jzp8at0vr/sKu95HVl5S4q36dq9JLvJrVjMXkp65p0ylQoDUUIHYQNS2Zat124cCHCw8OxcuVKc1lkZORv2hNIT0/HzJkzMWLECADAqlWrEBwcjDVr1mDixIk33FdbcCaJiIiIbkhZWZnFq6pK/kdbdnY2unfvjoSEBAQFBaFr165Yvny5+f2TJ0+iqKgIgwYNMpcZDAb07dsXO3bsaPDjUMNBEhERkZMxQWfzCwDCw8Ph7+9vfs2fP1+6vxMnTmDZsmVo164dPvvsMzz99NOYMmUKVq9eDQAoKro61RgcHGyxXXBwsPk9R+ByGxERkZOxV3Rbfn4+/Pz+t+5qMMhvkzCZTOjevTvS0q5mou/atSsOHTqEZcuWYfTo0eZ6Op1ln4QQirLGxJkkIiIiuiF+fn4WL7VBUmhoKDp06GBRFh0djTNnzgAAQkKuPmfr2lmj4uJixexSY+IgiYiIyMnU3bhty0uL3r174/DhwxZlR44cQUREBAAgKioKISEh+Pzzz83vG41G5OXloVevXrYf8A3ichsREZGTaexkktOmTUOvXr2QlpaGRx99FLt370ZGRgYyMjIAXF1mmzp1KtLS0tCuXTu0a9cOaWlp8PLywhNPPHHD/bQVB0nNmKfkXjatobNqod9SKqGzsrQDamH69kgjALVj9JGHHmsh659qn1XComVPTC+97CGtu/tChLwfrspQ4POd5OHxfic8peUeF60PJ5aF+gNAVaCyzKTyuZjc5F8nsvpqqSe8Dqk85f1n5cmWhVtfr38VbZTXh5tKyLvPWfm5k6UzqPGSX3e1khQMAGCSnKYat5bSuof85OWe55RlQaXy81Gj8nsrC71Xu6bVwuZl508tbF62P7V9qqbdUPldVCuXMajUlaURUP0OUiH7/nUvMUrr+hxRpnGorbp5UwD06NEDGzduREpKCubMmYOoqCikp6cjMTHRXOeFF17A5cuXMXnyZJSUlCAuLg65ubnw9fW94X7aioMkIiIianBDhw7F0KFDVd/X6XRITU1Fampq43WqHhwkERERORlh43KbLbNQzQkHSURERE5GABDWr1JKt3cGjG4jIiIikuBMEhERkZMxQQcdbIhus2Hb5oSDJCIiIifT2NFtzRUHSc1A5ymvS8sNGp6eLQtdVqMWQl0tjzSHXvLwd0OpSl0NIbxqddXSC8jCdbWGEsvaUE1noBLqLFyV9Ssvy7PQHi2Tn9Tay8pfzaqW8s+26rz8s/W4qCxTuw7cLkuLoZeEmldLnvwOABXh1n9pehWphJRfkodAe55Tnmz3onJpXeEu/1rTd/NXlKmFtqt95rL0Ey7V1j9VHpD/fhmK5J+tltB2NWr9k5PXNaikF5B932j5/QTkaUjUzp1eQ0i+WjoDtd9bGdXrQIXs96s6SuV3XJkBALW8AabJ4SCJiIjIyZiEDrpGTCbZXHGQRERE5GSEsDG6zUnC2zi5R0RERCTBmSQiIiInwxu3rePwmaSCggKMHDkSgYGB8PLyQpcuXbBv3z7V+mPHjoVOp1O8OnbsaFHv4sWLSEpKQmhoKDw8PBAdHY3NmzdL25w/f7754XrX+umnn/Dggw/C398fvr6+6NmzJ86cOWPTMRMRETlS3SDJlpczcOhMUklJCXr37o3+/fsjJycHQUFBOH78OAICAlS3Wbp0KRYsWGD+uaamBp07d0ZCQoK5zGg0YuDAgQgKCsL69evRpk0b5OfnSx+St2fPHmRkZCA2Nlbx3vHjx3HPPfdg/PjxmD17Nvz9/fHTTz/Bw0P+sFIiIqLmgDduW8ehg6SFCxciPDwcK1euNJdFRkZedxt/f3/4+/8vlHfTpk0oKSnBuHHjzGUrVqzAhQsXsGPHDri5XY0DjYhQPm29oqICiYmJWL58OebOnat4f+bMmXjggQewaNEic9ltt91m9fHZi+yp4ap1VULe1crVQm210NI/LSG19gh/1toPLftUS7fgc1L59Pfac97yNlTaloUvu6mkVVB7Yr1rhbLcTSWE2uQmn1SWpkSolPcDkpBmANBJQq71GsPE9e7K/pVHt5DWVUtzoOVavxKg/AwB+XWj9XdIduxanzYv26da6g61cy2vKy/XGgqvpQ03SSoT+5CfjxqVz8tVQ0oVtc9LVn45UF73SrDy99N0Wf67TI7j0OW27OxsdO/eHQkJCQgKCkLXrl2xfPlyTW1kZmZiwIABFoOg7OxsxMfHIykpCcHBwYiJiUFaWhpqay0vwKSkJAwZMgQDBgxQtGsymfDpp5+iffv2GDx4MIKCghAXF4dNmzbd0LESERE1FXXRbba8nIFDB0knTpzAsmXL0K5dO3z22Wd4+umnMWXKFKxevdqq7QsLC5GTk4MJEyYo2l2/fj1qa2uxefNmvPzyy1i8eDHmzZtnrrN27Vrs378f8+fPl7ZdXFyMiooKLFiwAPfddx9yc3Px8MMPY8SIEcjLy5NuU1VVhbKyMosXERFRU3N1oGPLPUmOPoLG4dDlNpPJhO7duyMtLQ0A0LVrVxw6dAjLli3D6NGj690+KysLAQEBGD58uKLdoKAgZGRkwMXFBd26dcPZs2fx6quvYtasWcjPz0dycjJyc3NV7y8yma5Osz700EOYNm0aAKBLly7YsWMH3n33XfTt21exzfz58zF79mwtp4CIiIiaKIfOJIWGhqJDhw4WZdHR0VZFjwkhsGLFCowaNQru7u6Kdtu3bw8Xl//dVxAdHY2ioiIYjUbs27cPxcXF6NatG1xdXeHq6oq8vDy88cYbcHV1RW1tLW655Ra4urpq6l9KSgpKS0vNr/z8fGtPBRERUaNhdJt1HDqT1Lt3bxw+fNii7MiRI9KbrK+Vl5eHY8eOYfz48dJ216xZA5PJBL1eb243NDQU7u7uuPfee3Hw4EGLbcaNG4c777wTM2bMgIuLC1xcXNCjRw9N/TMYDDAY5M/oIiIiairEf1+2bO8MHDpImjZtGnr16oW0tDQ8+uij2L17NzIyMpCRkWGuk5KSgoKCAsV9SpmZmYiLi0NMTIyi3UmTJuHNN99EcnIynn32WRw9ehRpaWmYMmUKAMDX11exnbe3NwIDAy3Kn3/+efzhD3/A//3f/6F///7YsmUL/vGPf2D79u12PAtERETUFDl0kNSjRw9s3LgRKSkpmDNnDqKiopCeno7ExERzncLCQsXyVmlpKTZs2IClS5dK2w0PD0dubi6mTZuG2NhYtG7dGsnJyZgxY4am/j388MN49913MX/+fEyZMgV33HEHNmzYgHvuuUf7wTaAlj9UWF85WP4kah9JyLvWp3jL2CPEXo1aOgMt/bBH22rhwa2+r1GU1fjIQ8rVwtVl51ptf2qqWih/vdXOh1r4uCy9g9tl+f7Uyj1/UYY1G0qU5whQP08l7ZQ5ES4Hy/fnqtIPryLlsfsUyPuhdh3IUgOonbtatfwOEi4qofc1KmHs2vZnfWoAtfB4k0r/NByiKtl1rfV7QvZ51aqkcdCSEkHtO0/t98jloiyEX96Py2WSa+mKvG5DYMZt6zj8sSRDhw7F0KFDVd/PyspSlPn7+6OyUi1Zy1Xx8fHYtWuX1f1Qmx168skn8eSTT1rdDhERUZPH9TarOHyQRERERI3M1puvnWQmyeHPbiMiIiJqijiTRERE5GRszZrNZJJERER0U+KN29bhchsRERGRBGeSmoFqf3m5yV0ZLqo3yp8i7V5itGeX6t1njbe7pCagr7b+KdcmN7VwWPnY3vWSSpyyRI23PHhZSwoAvdH6kHzXCm1P99ZLQo/V+qYWpqwltNrnrDwUXta2Wj/U2paFS8vSE1yvbdmT4l1+tv5p7oD8fKilYKgKsP7vR7WQcr3K5SgrN2mMpZe1oVdpQ60fsvOh5frXSi2Fhewz0JoCQJaaQS2dgdrnJfsMVK9HDek4fH6Wf/dW+SsTD9c23Ne0ktDZdvO1k8wkcZBERETkZHhPknW43EZEREQkwZkkIiIiZ8NkklbhIImIiMjJMLrNOlxuIyIiIpLgTBIREZEzcpIlM1twkNQMuKg8y9elXEO8qCRdACAPs1cL01cLyVcP1beNWj9cL8mPW19eZX3bRmX4rSO4XpKXq6VQsJXJXT55rJZGQEZL+DMgT39QGSU/vsuB8n6E7rz+A61/S0t6B7UUAGph87J0BlrD1aUh7xrC9AH5sZg0fIaA/Fi07A+QXzeydhuavN/yz1ZLugW13wu18yFLtaKWluX7N15UlJWVlcH/ry9Z30EbcLnNOhwkERERORveuG0V3pNEREREJMGZJCIiIqej++/Llu1vfhwkERERORsut1mFy21EREREEpxJIiIicjacSbIKB0lERETORuiuvmzZ3glwkNQMBByX59nQGWsUZTWBXtK6WnIZqeZDUsmxI8sZopY7Rq0NGb1RnrNFLe+IzqhMNiPc5UlRdFXyNuTtKs/z1bblvz6y+rK+XZ+Psg0NfVbrh8lXnh9Kb7Q+L5Na/iq1/ESXg5XlVX7yL1hDmdXdaPLU8gWZNORa0pRzyMv2/7TU8v9ozQUlo5aTSpq/SuVOELXzIS3XmM+rofI7XQn2bJB2qXFwkERERORkhLj6smV7Z8BBEhERkbPhPUlWYXQbERERkQRnkoiIiJwNb9y2CgdJRERETkYnrr5s2d4ZcJBERETkbHhPklU4SGoG1EJTZSHoauHxWqilAHC9pBLGfklZpNYPk7u8bbXwcRm1UHi1cH9pGyph/fJ25b8mwmB9WgWttIb7S9uQpB1wOS//DHVVDRemXNVCeZ7cLsvrGkrlYdsu5UZFmdpn6FIu/7yqWyqPUS0s3R6uBKik0pBcpmpXrmoaAZVQfRk3jaHwMmrpOLTc1qrlXGsNx3cvUV4farR816gfNzkLDpKIiIicDe9Jsgqj24iIiJyNsMOribty5YrNbXCQRERERDcFk8mEv/zlL2jdujV8fHxw4sQJAMCf//xnZGZmam6PgyQiIiJnc5POJM2dOxdZWVlYtGgR3N3/97ilTp064b333tPcHgdJREREzuYmHSStXr0aGRkZSExMhIvL/4InYmNj8Z///EdzexwkERER0U2hoKAAt99+u6LcZDKhulrrg8YZ3dYsuFaoPfXe+jB2qITey6g95V1LegG1EHbX8iqVtuVPp5e2rXLcspB3tbQAamH90roqof5q6Qy0/OWh6TPUSHbssnN0vX7UBHopylzPV0rrupTL++F7QnktuKk8Gb2ylfXXqdpnqHYsbhdU8g5oIPvdUEuZ4VItL691U0YFqYX0640qf657K4vcKuV11VIAyPapuj8N1MLm3S7Jj1F2PtT6oda27DNwvSRPC6D2my9LDWByl/82q/VDlqpCryHlQKO6SaPbOnbsiK+++goREREW5evWrUPXrl01t8dBEhERkZO5WTNuv/LKKxg1ahQKCgpgMpnwySef4PDhw1i9ejX++c9/am6Py21ERER0Uxg2bBg+/vhjbN68GTqdDrNmzcJPP/2Ef/zjHxg4cKDm9hw+SCooKMDIkSMRGBgILy8vdOnSBfv27VOtP3bsWOh0OsWrY8eOFvUuXryIpKQkhIaGwsPDA9HR0di8ebO0zfnz50On02Hq1Kmq+504cSJ0Oh3S09Nv5DCJiIiajpv0xm0AGDx4MPLy8lBRUYHKykp8/fXXGDRo0A21pXmQNHbsWHz55Zc3tLNrlZSUoHfv3nBzc0NOTg5+/PFHLF68GAEBAarbLF26FIWFheZXfn4+WrZsiYSEBHMdo9GIgQMH4tSpU1i/fj0OHz6M5cuXo3Xr1or29uzZg4yMDMTGxqruc9OmTfj2228RFhZm0/ESERFR86H5nqTy8nIMGjQI4eHhGDduHMaMGSMdfFhj4cKFCA8Px8qVK81lkZGR193G398f/v7+5p83bdqEkpISjBs3zly2YsUKXLhwATt27ICb29Wb5q69iQsAKioqkJiYiOXLl2Pu3LnS/RUUFOCZZ57BZ599hiFDhmg5PCIioiZJBxvvSbJbT+xLr9dDp1PvXW2ttudiap5J2rBhg3ngsG7dOkRGRuL+++/H+vXrNYfXZWdno3v37khISEBQUBC6du2K5cuXa2ojMzMTAwYMsBgEZWdnIz4+HklJSQgODkZMTAzS0tIUJycpKQlDhgzBgAEDpG2bTCaMGjUKzz//vGI5T6aqqgplZWUWLyIiImeXmpqquE0mJCTE/L4QAqmpqQgLC4Onpyf69euHQ4cOad7Pxo0b8cknn5hfH3/8MV588UWEhoYiIyNDc3s3FN0WGBiI5ORkJCcn47vvvsOKFSswatQo+Pj4YOTIkZg8eTLatWtXbzsnTpzAsmXLMH36dLz00kvYvXs3pkyZAoPBgNGjR9e7fWFhIXJycrBmzRpFu1u3bkViYiI2b96Mo0ePIikpCTU1NZg1axYAYO3atdi/fz/27Nmj2v7ChQvh6uqKKVOm1NsX4Oq9TbNnz7aqrhY1PvJQ4prbA6xuQ0tIrVoKAC3UwubV6FVSA8iohbFL+6ESJl7r6y4tl9GS+gCQpwYweVu/P0Aerq4W2q52jCZfZVoFl/NqKQDk5VqOXX++VN6PQH9FWWmEvM8+Z+X7k11PaikY1FIR6CuU51QtOFvL9SsL+wYA10vWp45QSyOgFsauNyqvJ7XfW7X+VbdUpmHQmv7DJElzoFZXr5ISQRZ6r5UsVN/o7iHvh0r4vusl5e+A2ueilhpAdt1UBjXRIHIHpADo2LEjvvjiC/PPv032uGjRIixZsgRZWVlo37495s6di4EDB+Lw4cPw9fW1eh8PPfSQouyRRx5Bx44d8fHHH2P8+PGa+mzTjduFhYXIzc1Fbm4uXFxc8MADD+DQoUPo0KEDXn/99Xq3N5lMuOuuu5CWloauXbti4sSJ+OMf/4hly5ZZtf+srCwEBARg+PDhinaDgoKQkZGBbt264bHHHsPMmTPN7ebn5yM5ORkffPABPDzkv0j79u3D0qVLkZWVdd2pu99KSUlBaWmp+ZWfn2/VdkRERI3KATduu7q6IiQkxPxq1arV1a4IgfT0dMycORMjRoxATEwMVq1ahcrKSsUkyI2Ki4uzGKBZS/Mgqbq6Ghs2bMDQoUMRERGBdevWYdq0aSgsLMSqVauQm5uL999/H3PmzKm3rdDQUHTo0MGiLDo6GmfOnKl3WyGEeQbrt89nqWu3ffv2FqPU6OhoFBUVwWg0Yt++fSguLka3bt3g6uoKV1dX5OXl4Y033oCrqytqa2vx1Vdfobi4GLfeequ5zunTp/Hcc8+p3jdlMBjg5+dn8SIiIiLg6NGjCAsLQ1RUFB577DHzw2dPnjyJoqIiiwg0g8GAvn37YseOHTbv9/Lly3jzzTfRpk0bzdtqngcMDQ2FyWTC448/jt27d6NLly6KOoMHD75uhFqd3r174/DhwxZlR44ckd5kfa28vDwcO3ZMOnXWu3dvrFmzBiaTCXq93txuaGgo3N3dce+99+LgwYMW24wbNw533nknZsyYARcXF4waNUpxr9LgwYMxatQoi5vEiYiImh1bw/j/u+21994aDAYYDMql/ri4OKxevRrt27fHuXPnMHfuXPTq1QuHDh1CUVERACA4ONhim+DgYJw+fVpTt1q0aGGx+iOEQHl5Oby8vPDBBx9oagu4gUHS66+/joSEBNVlqrpOnjx5st62pk2bhl69eiEtLQ2PPvoodu/ejYyMDIubq1JSUlBQUIDVq1dbbJuZmYm4uDjExMQo2p00aRLefPNNJCcn49lnn8XRo0eRlpZmvrfI19dXsZ23tzcCAwPN5YGBgQgMDLSo4+bmhpCQENxxxx31HhsREVFTZa+M2+Hh4Rblr7zyClJTUxX177//fvO/O3XqhPj4eLRt2xarVq1Cz549r7Z5za0tQgirb3ep8/rrr1tso9fr0apVK8TFxaFFixaa2gJuYJA0atQozTtR06NHD2zcuBEpKSmYM2cOoqKikJ6ejsTERHOdwsJCxfJbaWkpNmzYgKVLl0rbDQ8PR25uLqZNm4bY2Fi0bt0aycnJmDFjht36TkRE5Ozy8/Mtbi2RzSLJeHt7o1OnTjh69Kj5vuKioiKEhoaa6xQXFytml+ozduxYTfXr4/Db7ocOHYqhQ4eqvp+VlaUo8/f3R2Wl/CGbdeLj47Fr1y6r+7F9+/Z665w6dcrq9oiIiJosOy233ej9t1VVVfjpp5/Qp08fREVFISQkBJ9//rn5IbRGoxF5eXlYuHBhvW39+9//tnq/10scLePwQRL9T6fp8ohAtQnCKwHWhyl7FcvDx9VCWW2lFjore3o5oC3kXQu1NtTCxLWEfquFOpffpjxItRQMLtUqT26/oCzTkvrAXmSpGWoCvaR13covSctln4GLxkOpkaRQUAuPV0+VYHuoua7K+pQIqnVVUhdoIft9Ubt21c6HPX7ndJL0E2ptyOqqUfv+sEd6Eru0oZJGQJbe4VJIE027aKdBkrX+9Kc/YdiwYbj11ltRXFyMuXPnoqysDGPGjDE/FiwtLQ3t2rVDu3btkJaWBi8vLzzxxBP1tt2lSxfodDoIcf1O6XQ6zckkOUgiIiKiBvXzzz/j8ccfx6+//opWrVqhZ8+e2LVrlzlQ64UXXsDly5cxefJklJSUIC4uDrm5uVblSLLmHugbxUESERGRk7HXjdvWWrt27fXb0+mQmpoqvem7Pr+NiD9//rw56Co/Px/Lly/H5cuX8eCDD6JPnz6a226YtRYiIiJquuoybtvyakIOHjyIyMhIBAUF4c4778SBAwfQo0cPvP7668jIyED//v2xadMmze1ykERERORsHJBxuyG98MIL6NSpE/Ly8tCvXz8MHToUDzzwAEpLS1FSUoKJEydiwYIFmtvlchsRERE1a3v27MHWrVsRGxuLLl26ICMjA5MnTzYnlH722WfN+Zi04CCJiIjIyTT2PUkN7cKFCwgJCQEA+Pj4wNvbGy1btjS/36JFC5SXq4Q1XwcHSU2I53n5VVfZyvqQYZObfJ1Y7UnUHheV4ZB6eWS15idiy+uqrWNb/0RytaeaawlfVg1T1tBG4b23SMsnJmUryrJOxkvreqfKM9fLwrmFUR7CLrSEVvsoz/P1aEmJIHzl+R1MvsrkcgFHlKkFAPVrSUvYtpbzoXZ8slBuQH5/glqov6Z0EirHZ5eUA2r7rJCkEdCYJkHWhhq1z0X2vaKa3kHlGPWSc632GdqDpuvR9swTDaORUwA0hmuzc2vN1i3DQRIRERE1e2PHjjVn/L5y5QqefvppeHtf/eOtqkr+h1l9OEgiIiJyNjYutzW1maQxY8ZY/Dxy5EhFndGjR2tul4MkIiIiZ3OTLbetXLmyQdplCgAiIiIiCc4kEREROZubbCapoXCQRERE5GRuthQADYWDpCakyl8erlijErUtSxmgV3mqvNrT5mslKQP0GkL6tXKtsP0J3LW+yifCX6Us1xtt359aKHG1n7z+gfJbFWXFp1pKagI+3vKUA7IQaLUQai1pC2Th+Fq5nq9U6Ue1tPxKcICizOOcPHRcX3RjESg3SkuYfkNSS2uhhXpaC/nnIgv313yNVUmuU5VUEGrnWnatq4X6a7nWofJ7q3auZf1T+/5Q+06Qtd3mC3kbd+cvUZTVGq9I65Lj8J4kIiIiIgnOJBERETkb3pNkFQ6SiIiInAzvSbIOl9uIiIiIJDiTRERE5IycZDbIFhwkERERORvek2QVDpKaELVQf0OZ9W3o5dG+cLtk0t6ha9s22t6GatuSp2prDb+V1df6JHDZE8nVnkzvf0J+Pva/21lRdttJedix2wV5KLyW0HS1EG+Tj8oFpYG+XBmSr55GQF5e7a08f+5qIdQqxyKjdnxaQuF1VfKvQC33IainpJCTXWPQ2IbsWtfaDy1cylXekIT7q6URUCO7xtTa0JKiQG+0Pb1Djbf8nMq+rwCVNAKS4wMAr2IP5f5qNKQ4oEbBQRIREZGT4Y3b1uEgiYiIyNlwuc0qjG4jIiIikuBMEhERkZPhcpt1OEgiIiJyNlxuswqX24iIiIgkOJPUDLhWyofssrB+l2p5Xb3R+mG/Wqi/WtirvA1t4ftaqLXdUNTOh+c5ebm3hvOkJdRfrW6twUdabo9zDUm4/6U2Xir700nLPc8pQ+9Vrw8NaQvUzodaCgDpU+9V2lA7d1pSRKidDy2/i1rUeCuPTyvVPtuhbTWuknPdkL/j1S1tT40hTeMA+WcgOz4AqGqhLK+ptsPvrLU4k2QVDpKIiIicDO9Jsg4HSURERM6GM0lW4T1JRERERBKcSSIiInI2nEmyCgdJRERETob3JFmHy21EREREEpxJIiIicjZcbrOKwwdJBQUFmDFjBnJycnD58mW0b98emZmZ6Natm7T+2LFjsWrVKkV5hw4dcOjQIfPPFy9exMyZM/HJJ5+gpKQEUVFRWLx4MR544AHFtvPnz8dLL72E5ORkpKenAwCqq6vx8ssvY/PmzThx4gT8/f0xYMAALFiwAGFhYfY5+GvolelkAAAeF+U5Q2rdlDlNqr3lk4N6N/kVLcur5HpJJYeNSm4Q10tGRZmuSt6G2tSllpw+anW15FZROxZZLiitOVvscSwyjZ0fCgBqvN0VZbL8XABgqlbJsaMhb1Str3J/mklyOwHya1JLPiQ1anm01K52WX3tuZaU9V0vqXyBaHDZW37uKlup/b4oy0wqKZXUv9+Ux+haoe1a17tbX7/GR8tnq5Z3Tu0zl9Vt/N9ba3C5zToOHSSVlJSgd+/e6N+/P3JychAUFITjx48jICBAdZulS5diwYIF5p9ramrQuXNnJCQkmMuMRiMGDhyIoKAgrF+/Hm3atEF+fj58fX0V7e3ZswcZGRmIjY21KK+srMT+/fvx5z//GZ07d0ZJSQmmTp2KBx98EHv37rX94ImIiKhJc+ggaeHChQgPD8fKlSvNZZGRkdfdxt/fH/7+/uafN23ahJKSEowbN85ctmLFCly4cAE7duyAm9vVP2siIiIUbVVUVCAxMRHLly/H3LlzFfv5/PPPLcrefPNN3H333Thz5gxuvfVWq4+TiIioSeFym1UceuN2dnY2unfvjoSEBAQFBaFr165Yvny5pjYyMzMxYMAAi0FQdnY24uPjkZSUhODgYMTExCAtLQ21tZbTnklJSRgyZAgGDBhg1b5KS0uh0+lUZ7qqqqpQVlZm8SIiImpyhB1eTsChg6QTJ05g2bJlaNeuHT777DM8/fTTmDJlClavXm3V9oWFhcjJycGECRMU7a5fvx61tbXYvHkzXn75ZSxevBjz5s0z11m7di3279+P+fPnW7WvK1eu4MUXX8QTTzwBPz8/aZ358+ebZ7r8/f0RHh5uVdtERETU9Dh0uc1kMqF79+5IS0sDAHTt2hWHDh3CsmXLMHr06Hq3z8rKQkBAAIYPH65oNygoCBkZGXBxcUG3bt1w9uxZvPrqq5g1axby8/ORnJyM3NxceHh41Luf6upqPPbYYzCZTHjnnXdU66WkpGD69Onmn8vKyjhQIiKiJkf335ct2zsDhw6SQkND0aFDB4uy6OhobNiwod5thRBYsWIFRo0aBXd3y4iY0NBQuLm5wcXlf1EM0dHRKCoqgtFoxL59+1BcXGwRQVdbW4svv/wSb731FqqqqszbVldX49FHH8XJkyexdetW1VkkADAYDDAY5NEhRERETQbvSbKKQwdJvXv3xuHDhy3Kjhw5Ir3J+lp5eXk4duwYxo8fL213zZo1MJlM0Ov15nZDQ0Ph7u6Oe++9FwcPHrTYZty4cbjzzjsxY8YMxQDp6NGj2LZtGwIDA2/0UK3id0YeKup5pkJafvlWH5v3qRbiKuN24bK0XGesUZQJd/mlpSXkWksYtlZawtK1phzQlIqggdIFaOVSrkzjAAB6g+RzqVb5DI22f7Zqn4usvpbPEACE5FjUPiu1/qmF6jcUWZoP9XJ57L1auLosdYdqegeVVATVXtbPJ1T5y+u6XVKW61XOs9qxyD5HY4v6Vwnqo56CQV5floahIX9vbcEUANZx6D1J06ZNw65du5CWloZjx45hzZo1yMjIQFJSkrlOSkqKdOktMzMTcXFxiImJUbw3adIknD9/HsnJyThy5Ag+/fRTpKWlmdv19fVFTEyMxcvb2xuBgYHm9mpqavDII49g7969+PDDD1FbW4uioiLzbBQRERHd3Bw6k9SjRw9s3LgRKSkpmDNnDqKiopCeno7ExERzncLCQpw5c8Ziu9LSUmzYsAFLly6VthseHo7c3FxMmzYNsbGxaN26NZKTkzFjxgyr+/bzzz8jOzsbANClSxeL97Zt24Z+/fpZ3RYREVGTwuU2qzg84/bQoUMxdOhQ1fezsrIUZf7+/qisrLxuu/Hx8di1a5fV/di+fbvFz5GRkRDCSa4CIiJyPvwvrl58wC0RERGRhMNnkoiIiKhx8cZt63CQRERE5Gx4T5JVOEhqQnyOl0vL9YW/SMu9K5Qh+cJd5RHcGuiMKo/rLpenIpC24astPYEs1FxNQ4bUykK/ZaHSgPyp8oA81NwR7PH0cdm5tkdqBq3h+7LPQDWdhB2uD7VQfy1Pf1d/2ryy7Wpv+f7UymXUwtVr3eRf8y7Vyu8KQ4n8c1ErV2tbvj95eVWA8hjVjsXznFqKAus/cy1pT7R83lfrK89TU00BQNbhIImIiMjJcLnNOhwkERERORsut1mF0W1EREREEpxJIiIicjJcbrMOB0lERETOhsttVuEgiYiIyNlwkGQVDpIc5O4xSxRlgedLtTVSpQyL1knKtBIaQv1VVV2QFuvL3eX13SXlBnldLQG1aikR1NIcmHw8rW5bL0nBAADCqNxn5e0B0rruJQ33sGQtoccmb/m5loXCaw2LlrarkkZALTWApmPR0LY90hmopwuQ/y8iC293qVb5H+eS/FzLUgPUusnD5rVQC71X43FReU7V+uFSrZaiQFmufu7k51qtXAvZda12Paql12C4/82HgyQiIiInw3uSrMNBEhERkbPhcptVmAKAiIiISIIzSURERE5GJwR04sang2zZtjnhIImIiMjZcLnNKlxuIyIiIpLgTJKDuMlCe43ycHChEtavKVhXFmKvsk+dr4/VdVXbVqurRlZfaxsSOrXjViH9q0EtrYJK/2T79P5R5RHoDUiW/kC4q/zK+8rPk1oItBZawuzV6tojxNse7db42B7ire0p9NbX1Rq+LyNLLXA9sv6ppTOQfudB3m+141ZLPyH7XOzRhlrCEdcK+e+F6yVZWhZ53Vo3ZboRk7D9M7QWo9usw5kkIiIiZyPs8LLB/PnzodPpMHXq1P91SQikpqYiLCwMnp6e6NevHw4dOmTbjmzEQRIRERE1mj179iAjIwOxsbEW5YsWLcKSJUvw1ltvYc+ePQgJCcHAgQNRXl7uoJ5ykEREROR06pbbbHndiIqKCiQmJmL58uVo0aKFuVwIgfT0dMycORMjRoxATEwMVq1ahcrKSqxZs8ZOR60dB0lERETOxk7LbWVlZRavqqqq6+42KSkJQ4YMwYABAyzKT548iaKiIgwaNMhcZjAY0LdvX+zYscPmw71RHCQRERE5GXvNJIWHh8Pf39/8mj9/vuo+165di/3790vrFBUVAQCCg4MtyoODg83vOQKj24iIiOiG5Ofnw8/Pz/yzwWBQrZecnIzc3Fx4eHiotqfTWUb4CSEUZY2JgyQHkYbaqoWrq4Wgy2gMeddUXy01QEOxQwoATWkLAPm5biLpDFRTQRjkxyL7WlH7qtFXaPhsVfohfL3l5WppBzQQBmUodo23/LjVw/QbZuJc/an3jRsjrRbyrhbWL+uf2rGYlNkkru5TUt+kdtzyywOGEmWIvFqYvlpKCr1ReYxa21AL99fC5C5JRaC2N8l5EjWNeM3YKZmkn5+fxSBJzb59+1BcXIxu3bqZy2pra/Hll1/irbfewuHDhwFcnVEKDQ011ykuLlbMLjUmLrcRERE5oca8afvee+/FwYMHceDAAfOre/fuSExMxIEDB3DbbbchJCQEn3/+uXkbo9GIvLw89OrVy45HrQ1nkoiIiKhB+fr6IiYmxqLM29sbgYGB5vKpU6ciLS0N7dq1Q7t27ZCWlgYvLy888cQTjugyAA6SiIiInI8QV1+2bG9nL7zwAi5fvozJkyejpKQEcXFxyM3Nha+vr933ZS0OkoiIiJxMU3gsyfbt2y3b1OmQmpqK1NRU2xu3E96TRERERCTBmSQiIiJnY6fotpsdB0kO4nHR9qera6ISJm4PsqfNq2a10JLOoCFpCMlXC71XoxaSr4XWfdrahg4VVrchrsgz6urK5W3o7XHtSVI2uKi06+ajfLq6VrW+8rZNRmWIt9qXqFqoeY23Sjy9rA2VMHYXd+vzxrho+K5R65laagAZtZQDamkEDCXKMvUwfTnXS9VW19Ub5W27Vli/Ty39k6UFuNoP5ShD34hpI3Smqy9btncGXG4jIiIikuBMEhERkbPhcptVOEgiIiJyMk0huq054CCJiIjI2TTBPElNEe9JIiIiIpJw+CCpoKAAI0eORGBgILy8vNClSxfs27dPtf7YsWOh0+kUr44dO1rUu3jxIpKSkhAaGgoPDw9ER0dj8+bN0jbnz58PnU6HqVOnWpQLIZCamoqwsDB4enqiX79+OHTokM3HTERE5Ei2PLfN1qW65sShy20lJSXo3bs3+vfvj5ycHAQFBeH48eMICAhQ3Wbp0qVYsGCB+eeamhp07twZCQkJ5jKj0YiBAwciKCgI69evR5s2bZCfny9Nbb5nzx5kZGQgNjZW8d6iRYuwZMkSZGVloX379pg7dy4GDhyIw4cPN0yadLWnuRvl5bLwbOsDdbX3QxbqDwA6o/Xhtw2ZikBLWL8WaiH9jR2mr8ZUWi4t13kYbG5bC7XUAHYhu9ZVfi/0Ws6pymerr7gsLVf7HdDC9bzNTcDN3fqvbmGw/un2auHqJjeVcnfb/86WpTnQkiZBrQ01aseiNe2AvB/KNmq85dfYN+v/pCgrKyuDv/+fbe6HVXjjtlUcOkhauHAhwsPDsXLlSnNZZGTkdbfx9/eHv7+/+edNmzahpKQE48aNM5etWLECFy5cwI4dO+DmdvWXLSIiQtFWRUUFEhMTsXz5csydO9fiPSEE0tPTMXPmTIwYMQIAsGrVKgQHB2PNmjWYOHGi5uMlIiKi5sOhy23Z2dno3r07EhISEBQUhK5du2L58uWa2sjMzMSAAQMsBkHZ2dmIj49HUlISgoODERMTg7S0NNTWWo7yk5KSMGTIEAwYMEDR7smTJ1FUVIRBgwaZywwGA/r27YsdO3ZI+1JVVYWysjKLFxERUVPD5TbrOHSQdOLECSxbtgzt2rXDZ599hqeffhpTpkzB6tWrrdq+sLAQOTk5mDBhgqLd9evXo7a2Fps3b8bLL7+MxYsXY968eeY6a9euxf79+zF//nxp20VFRQCA4OBgi/Lg4GDze9eaP3++eabL398f4eHhVh0HERFRo6qLbrPl5QQcutxmMpnQvXt3pKWlAQC6du2KQ4cOYdmyZRg9enS922dlZSEgIADDhw9XtBsUFISMjAy4uLigW7duOHv2LF599VXMmjUL+fn5SE5ORm5uLjw8PK67D53O8i4fIYSirE5KSgqmT59u/rmsrIwDJSIiombKoYOk0NBQdOjQwaIsOjoaGzZsqHdbIQRWrFiBUaNGwf2a5zqFhobCzc0NLi7/u0EvOjoaRUVFMBqN2LdvH4qLi9GtWzfz+7W1tfjyyy/x1ltvoaqqCiEhIQCuziiFhoaa6xUXFytml+oYDAYYDI17sywREZFWTCZpHYcut/Xu3RuHDx+2KDty5Ij0Jutr5eXl4dixYxg/fry03WPHjsFk+l/Ew5EjRxAaGgp3d3fce++9OHjwIA4cOGB+de/eHYmJiThw4ABcXFwQFRWFkJAQfP755+Y2jEYj8vLy0KtXLxuOmoiIyMGEHV5OwKEzSdOmTUOvXr2QlpaGRx99FLt370ZGRgYyMjLMdVJSUlBQUKC4TykzMxNxcXGIiYlRtDtp0iS8+eabSE5OxrPPPoujR48iLS0NU6ZMAQD4+voqtvP29kZgYKC5vC5vUlpaGtq1a4d27dohLS0NXl5eeOKJJ2w+dvcSSZiyHZ42rxpSrlIuDW9XCaXXSZ7ErplaCgAtx67SP3uE06uF+2thj35ooRbqLwvJV62r5fpQoZYCwB6pATSlM1BLBaHh+tUS6q+aAsMe14FaWpDySw3StlqyALXzISSpCNw1pBwAAH258vqQtQsAtb7y86EW1i/jekn+uailP9CyP1kKAHukSSDHceggqUePHti4cSNSUlIwZ84cREVFIT09HYmJieY6hYWFOHPmjMV2paWl2LBhA5YuXSptNzw8HLm5uZg2bRpiY2PRunVrJCcnY8aMGZr698ILL+Dy5cuYPHkySkpKEBcXh9zc3IbJkURERNRIuNxmHYc/u23o0KEYOnSo6vtZWVmKMn9/f1RWVl633fj4eOzatcvqfmzfvl1RptPpkJqaitTUVKvbISIiavJM4urLlu2dgMMHSURERNTImHHbKlwsJSIiIpLgTBIREZGT0cHGe5Ls1pOmjYMkIiIiZ2Nr1mwnybjN5TYiIiIiCc4kOYjreWV0niiv0NSGLIeN1jw/svw49sgVpJqrRq1cg4bMQ6TlfKiVN3aeJDWacgupaKjrwx65kzTvU/L7pYOPtK7aUoKW/En2yNdkj/xhmqi1odIPvexasEM+NbVrTF9xWVpu8vG0uq7asbjI9ql2/jX8DrhqyL/UmJgCwDocJBERETkbRrdZhcttRERERBKcSSIiInIyOiGgs+Hma1u2bU44SCIiInI2pv++bNneCXC5jYiIiEiCM0lEREROhstt1uEgyUF050sUZSaVsGi9v6/N+9MSlq5WV0uGVdU21EJnZWHD9ghptgNTabmm+vYIvW8o9kgRofmz1bK/BkwNIO2fxlQV0t8BO4S8q1I7p5J0BlpTT+gast+20ngs0iURrak4tJxTla8Ena8ypYTeWKutH42F0W1W4SCJiIjI2TDjtlV4TxIRERGRBGeSiIiInAwzbluHgyQiIiJnw+U2q3C5jYiIiEiCM0lERERORme6+rJle2fAQVIjGDpsCVxdPSzK3ANbKOo15LRegz6xXhZKrNKulvQC9uibPULe1UL6Gztc3S6flQp7tG2Xz0vlXDf2+dCUqsIOTIH+8n4Ya+QbyFIzaE0hoiHFhpbvDy2pQgDYJ/2HrL7aZ6XWtqS+2rGoXnuSto0tPCQVmwAut1mFy21EREREEpxJIiIicjZMJmkVDpKIiIicDB9LYh0utxERERFJcCaJiIjI2fDGbatwkERERORsBABbwvidY4zEQVJj0FfXQi8snwQt3JWnXu2p3MLX2+p96covyd9ooHBf4AZCfjW03VTbbWgN1W9Tqcrjy1WoheQ3tsb+HFX3pyH0XjWdgaRMX3HZ2q5d3aeWumq/txrTY1hL7RpTvZYaKtWHPdIIaCQ719XeTfOuFt6TZJ2m+ekRERERORhnkoiIiJyNgI33JNmtJ00aB0lERETOhjduW4XLbUREREQSnEkiIiJyNibYFnXDB9wSERHRzYjRbdbhIKkRZG95Hn5+fo2yr0E950jL9eUqYcrGakWZcHezui4AoLzCus5BW7i01hBlWfit1vBsLbQ8sV7z0+Ml4cj2CIPXetz2CBNvqPD9ptw3QP3ak5ZrDZu3sS7ggLQKKuejOdJy7Znc7ZEkhRyFgyQiIiJnwxu3rcJBEhERkbPhIMkqDo9uKygowMiRIxEYGAgvLy906dIF+/btU60/duxY6HQ6xatjx44W9S5evIikpCSEhobCw8MD0dHR2Lx5s/n9ZcuWITY2Fn5+fvDz80N8fDxycnIs2qioqMAzzzyDNm3awNPTE9HR0Vi2bJl9TwARERE1SQ6dSSopKUHv3r3Rv39/5OTkICgoCMePH0dAQIDqNkuXLsWCBQvMP9fU1KBz585ISEgwlxmNRgwcOBBBQUFYv3492rRpg/z8fPj6+prrtGnTBgsWLMDtt98OAFi1ahUeeughfPfdd+YB17Rp07Bt2zZ88MEHiIyMRG5uLiZPnoywsDA89NBDdj4bREREjYQzSVZx6CBp4cKFCA8Px8qVK81lkZGR193G398f/v7+5p83bdqEkpISjBs3zly2YsUKXLhwATt27ICb29WbkCMiIizaGTZsmMXP8+bNw7Jly7Br1y7zIGnnzp0YM2YM+vXrBwB46qmn8Ne//hV79+7lIImIiJovpgCwikOX27Kzs9G9e3ckJCQgKCgIXbt2xfLlyzW1kZmZiQEDBlgMgrKzsxEfH4+kpCQEBwcjJiYGaWlpqK2tlbZRW1uLtWvX4tKlS4iPjzeX33PPPcjOzkZBQQGEENi2bRuOHDmCwYMHS9upqqpCWVmZxYuIiKipqUsBYMvLGTh0JunEiRNYtmwZpk+fjpdeegm7d+/GlClTYDAYMHr06Hq3LywsRE5ODtasWaNod+vWrUhMTMTmzZtx9OhRJCUloaamBrNmzTLXO3jwIOLj43HlyhX4+Phg48aN6NChg/n9N954A3/84x/Rpk0buLq6Qq/X47333sM999wj7c/8+fMxe/bsGzwb9pG7a1b9lZqg+2+dan1lDU/r1vv7Sssb9MnosnB/tXZV+iF7krpd0hY00JPfte5T8/nXkkJB7fqQtKH2h7Qj0i1I21A5H1r6p6UNe4Tpm6quSMv1Bg+b29aUdkNFg6Z9CGyhKKvyd/itv2QDhw6STCYTunfvjrS0NABA165dcejQISxbtsyqQVJWVhYCAgIwfPhwRbtBQUHIyMiAi4sLunXrhrNnz+LVV1+1GCTdcccdOHDgAC5evIgNGzZgzJgxyMvLMw+U3njjDezatQvZ2dmIiIjAl19+icmTJyM0NBQDBgxQ9CclJQXTp083/1xWVobw8PAbOTVEREQNh/ckWcWhQ9zQ0FCLmRsAiI6OxpkzZ+rdVgiBFStWYNSoUXC/5i/E0NBQtG/fHi4uLhbtFhUVwfibvzLd3d1x++23o3v37pg/fz46d+6MpUuXAgAuX76Ml156CUuWLMGwYcMQGxuLZ555Bn/4wx/w2muvSftkMBjM0XJ1LyIioibHJGx/aVBfRLkQAqmpqQgLC4Onpyf69euHQ4cO2fuoNXPoIKl37944fPiwRdmRI0cUN1nL5OXl4dixYxg/fry03WPHjsFk+t+dZUeOHEFoaKhiQPVbQghUVV2dbq6urkZ1dTX0estT5OLiYtEuERERXV9dRPnevXuxd+9e/O53v8NDDz1kHggtWrQIS5YswVtvvYU9e/YgJCQEAwcORHm5PBt9Y3HoIGnatGnYtWsX0tLScOzYMaxZswYZGRlISkoy10lJSZEuvWVmZiIuLg4xMTGK9yZNmoTz588jOTkZR44cwaeffoq0tDSLdl966SV89dVXOHXqFA4ePIiZM2di+/btSExMBAD4+fmhb9++eP7557F9+3acPHkSWVlZWL16NR5++OEGOBtERESNpG65zZaXBsOGDcMDDzyA9u3bo3379pg3bx58fHywa9cuCCGQnp6OmTNnYsSIEYiJicGqVatQWVmpuOe4sTn0nqQePXpg48aNSElJwZw5cxAVFYX09HTzQAW4enP2tctvpaWl2LBhg3lp7Frh4eHIzc3FtGnTEBsbi9atWyM5ORkzZsww1zl37hxGjRqFwsJC+Pv7IzY2Flu2bMHAgQPNddauXYuUlBQkJibiwoULiIiIwLx58/D000/b+UwQERE1JhvvScLVba+N4jYYDDAYrh+0UFtbi3Xr1pkjyk+ePImioiIMGjTIop2+fftix44dmDhxog39tI3DH0sydOhQDB06VPX9rKwsRZm/vz8qKyuv2258fDx27dql+n5mZma9fQsJCbHI4URERET/c21w0iuvvILU1FRpXbWI8h07dgAAgoODLeoHBwfj9OnTDdJvazl8kEQEADln0h3dBc3u6zhTWi7crf+10rm7Scv1kjB2raHL9gj3b6hwac2h/hrSPqgR5RU2tyGjNdTfEWkYZBqqH2r3cDTXFBYyqiksjNWKsv3LpjV0d26MnaLb8vPzLYKUrjeLpBZRXkens0zKIYRQlDU2DpKIiIicjUmgbsnsxreHpkjuuohyAOjevTv27NmDpUuXmm+FKSoqQmhoqLl+cXGxYnapsTHLFRERETW6uojyqKgohISE4PPPPze/ZzQakZeXh169ejmwh5xJIiIicj7CdPVly/YavPTSS7j//vsRHh6O8vJyrF27Ftu3b8eWLVug0+kwdepUpKWloV27dmjXrh3S0tLg5eWFJ5544sb7aAccJBERETmbRs64XV9E+QsvvIDLly9j8uTJKCkpQVxcHHJzc+HrK3+0VGPhIImIiMjZ2OmeJGvVF1Gu0+mQmpqqGhnnKLwniYiIiEiCM0lEN2jLoXmO7oJmg7u+oqm+zlgjKVOGOQMAVELspU+bV0stYIen22uh9an3sjB2tX6oHuN1Ho2koCH1geZjkfTbHudU8/loZA2aGqOB0kw0CD7g1iocJBERETkbARsHSXbrSZPG5TYiIiIiCc4kERERORsut1mFgyQiIiJnYzIBsCFPksmGbZsRLrcRERERSXAmiYiIyNlwuc0qHCQROZHPvpvt6C40a/e3e0FRZvLx1NSGlv9a9OdL5W9IQtBl6QkA9dQAslB41afbN3JqAC191swOKRhU+6GlbUfjIMkqXG4jIiIikuBMEhERkbNp5MeSNFccJBERETkZIUwQ4sYj1GzZtjnhIImIiMjZCGHbbBDvSSIiIiJyXpxJIiIicjbCxnuSnGQmiYMkIiIiZ2MyATob7iviPUlERPRbOUcXOboLdnNfx5lW11WbM9CVX7J+h1pzDsn6oSWPk1rOIrX9qbStRXVkK5vboKaFgyQiIiJnw+U2q3CQRERE5GSEyQRhw3Kbs6QAYHQbERERkQRnkoiIiJwNl9uswkESERGRszEJQMdBUn243EZEREQkwZkkIiIntOXQPEd3QbPBXV+RlsvmNITBRVpXX14lLddp6YhKOgNjC+vTGTicEABsyZPkHDNJHCQRERE5GWESEDYstwkOkoiIiOimJEywbSaJKQCIiIiInBZnkoiIiJwMl9usw0ESERGRs+Fym1U4SGpAdSPtsrIyB/eEiKj5q6mVR6bJiBr53ST6WpWH5JqqlYUmlYfeyuoCqKm+oijT8v1fV7cxZmlqUG1TLskayM/BzUYnnGXOzAF+/vlnhIeHO7obRETUjOTn56NNmzYN0vaVK1cQFRWFoqIim9sKCQnByZMn4eHhYYeeNU0cJDUgk8mEs2fPwtfXFzqdpiwcjaKsrAzh4eHIz8+Hn5+fo7ujCfvuGM2570Dz7j/77hiN2XchBMrLyxEWFga9vuHiqq5cuQKjSq4nLdzd3W/qARLA5bYGpdfrG+yvAXvy8/Nrdl9cddh3x2jOfQead//Zd8dorL77+/s3+D48PDxu+sGNvTAFABEREZEEB0lEREREEhwkOTGDwYBXXnkFBoPB0V3RjH13jObcd6B59599d4zm3HeyHW/cJiIiIpLgTBIRERGRBAdJRERERBIcJBERERFJcJBEREREJMFBUhM1f/589OjRA76+vggKCsLw4cNx+PBhizqffPIJBg8ejFtuuQU6nQ4HDhyweP/UqVPQ6XTS17p168z1IiMjFe+/+OKLFm2dOXMGw4YNg7e3N2655RZMmTJFkbH14MGD6Nu3L9zc3ODu7g6DwdCgfd++fbtqnT179pjbkr3/7rvvSvvu6ekJPz8/tG7d2qZzDwBFRUUYNWoUQkJC4O3tjbvuugvr16+3qFNSUoJRo0bB398f/v7+GDVqFC5evOjQc29N30+dOoXx48cjKioKnp6eaNu2LV555RVFv7Sc+8bqO9A0r3lr+t7Ur/njx4/j4YcfRqtWreDn54dHH30U586ds6jTVK/5+vreENe8p6cnWrdujTlz5jTK89pIOw6Smqi8vDwkJSVh165d+Pzzz1FTU4NBgwbh0qVL5jqXLl1C7969sWDBAmkb4eHhKCwstHjNnj0b3t7euP/++y3qzpkzx6Leyy+/bH6vtrYWQ4YMwaVLl/D1119j7dq12LBhA5577jlznbKyMgwcOBBhYWGIi4vDU089BVdXV4waNarB+t6rVy9FnQkTJiAyMhLdu3e3aG/lypUW9caMGSPt+549e9C2bVtcuHABTz/99A2fewAYNWoUDh8+jOzsbBw8eBAjRozAH/7wB3z33XfmOk888QQOHDiALVu2YMuWLThw4ABGjRrl0HNvTd//85//wGQy4a9//SsOHTqE119/He+++y5eeuklRVvWnvvG6nudpnbNW9P3pnzNX7p0CYMGDYJOp8PWrVvxzTffwGg0YtiwYTCZ/vfE+KZ4zVvT94a45vfs2YM333wTr732GpYsWSLtGzmYoGahuLhYABB5eXmK906ePCkAiO+++67edrp06SKefPJJi7KIiAjx+uuvq26zefNmodfrRUFBgbnso48+EgaDQZSWlgohhHjnnXeEv7+/uHLlirnO/PnzRVhYmDh37lyD9f23jEajCAoKEnPmzLEoByA2btyout31+m4ymW743Ht7e4vVq1dblLVs2VK89957QgghfvzxRwFA7Nq1y/z+zp07BQDxn//8RwjhuHNfX99lFi1aJKKioizKbDn3Ddn3pnrNaz3vTema/+yzz4RerzefHyGEuHDhggAgPv/8cyFE073mrem7jL2veZPJpLodOQZnkpqJ0tJSAEDLli1vuI19+/bhwIEDGD9+vOK9hQsXIjAwEF26dMG8efMsppB37tyJmJgYhIWFmcsGDx6Mqqoq7Nu3z1ynb9++FgnXBg8ejLNnz+LQoUMN2vc62dnZ+PXXXzF27FjFe8888wxuueUW9OjRA++++67FX7bX6/upU6du+Nzfc889+Pjjj3HhwgWYTCasXbsWVVVV6Nevn3m//v7+iIuLM2/Ts2dP+Pv7Y8eOHeY6jjj39fVdprS0VLqfGz33Dd33pnjNaz3vTemar6qqgk6ns2jTw8MDer0eX3/9tXm/TfGat6bvMva+5k+dOmV1n6lx8AG3zYAQAtOnT8c999yDmJiYG24nMzMT0dHR6NWrl0V5cnIy7rrrLrRo0QK7d+9GSkoKTp48iffeew/A1fskgoODLbZp0aIF3N3dUVRUZK4TGRlpUadum1mzZjVY36+tM3jwYISHh1uU/+Uvf8G9994LT09P/Otf/8Jzzz2HX3/91by8cr2+FxYWYsGCBTfU/48//hh/+MMfEBgYCFdXV3h5eWHjxo1o27ateb9BQUGK7YKCgizOqyPOfX19v9bx48fx5ptvYvHixRbltpz7hux7U73mtZ73pnTN9+zZE97e3pgxYwbS0tIghMCMGTNgMplQWFho3m9TvOat6fu1GuKaLyoqQlRUlNX9pobHQVIz8Mwzz+Df//73df+iqc/ly5exZs0a/PnPf1a8N23aNPO/Y2Nj0aJFCzzyyCPmv7SBqzcjXksIYVF+bR3x3xsRjx07ZnFTqT37Xufnn3/GZ599hr/97W+K9357r0mXLl0AXL0f5bflan1fvHjxDZ/7l19+GSUlJfjiiy9wyy23YNOmTUhISMBXX32FTp06Sfdbt+/rnVdr6th67q3pe52zZ8/ivvvuQ0JCAiZMmKBop47Wc9+QfW+q17yW897UrvlWrVph3bp1mDRpEt544w3o9Xo8/vjjuOuuu+Di4qK637p9O/Kat7bvdRrqmpcdNzkWB0lN3LPPPovs7Gx8+eWXaNOmzQ23s379elRWVmL06NH11u3ZsyeAq182gYGBCAkJwbfffmtRp6SkBNXV1ea/gEJCQsx/5dWpu9Hy448/bvC+r1y5EoGBgXjwwQfrba9nz54oKyvDuXPnEBwcLO17cXExgKtT4998843m/h8/fhxvvfUWfvjhB3Ts2BEA0LlzZ3z11Vd4++238e677yIkJEQR+QMAv/zyi8V5bexzb03f65w9exb9+/dHfHw8MjIy6m3bmnPfWH2/tl+AY695rX1vatc8AAwaNAjHjx/Hr7/+CldXVwQEBCAkJMQ8O9JUr3lr+l6nIa75unN/7QwaOR7vSWqihBB45pln8Mknn2Dr1q02T8FmZmbiwQcfRKtWreqtWxdJExoaCgCIj4/HDz/8YDHtnJubC4PBgG7dupnrfPnllzAajea+5+TkICgoCH369GnQvgshsHLlSowePRpubm5WHZ+HhwcCAgIUfa9rLykpCXq9Htu3b7+hc19ZWQkA0Ostf8VcXFzM9yfEx8ejtLQUu3fvNr//7bfforS01Lys6Ihzb03fAaCgoAD9+vXDXXfdhZUrVyrqy1zv3Ddm32X9Ahx7zWvpe1O85n/rlltuQUBAALZu3Yri4mLzQK6pXvPW9B2w/zX/2+MLCwtTLMNRE9BIN4iTRpMmTRL+/v5i+/btorCw0PyqrKw01zl//rz47rvvxKeffioAiLVr14rvvvtOFBYWWrR19OhRodPpRE5OjmI/O3bsEEuWLBHfffedOHHihPj4449FWFiYePDBB811ampqRExMjLj33nvF/v37xRdffCHatGkjnnnmGXOdixcviuDgYPH444+LP/zhD8LLy0t4eXmJWbNmNVjf63zxxRcCgPjxxx8V72VnZ4uMjAxx8OBBcezYMbF8+XLh5+cnpkyZIu37wYMHxeDBgwUAMWnSpBs+90ajUdx+++2iT58+4ttvvxXHjh0Tr732mtDpdOLTTz81t3PfffeJ2NhYsXPnTrFz507RqVMnMXToUIeee2v6XlBQIG6//Xbxu9/9Tvz8888W5+lGz31j9b2pXvPWXjNN9ZoXQogVK1aInTt3imPHjon3339ftGzZUkyfPt2if03xmrem7w1xzR88eFB88sknws/PT7z22muKz5Icj4OkJgqA9LVy5UpznZUrV0rrvPLKKxZtpaSkiDZt2oja2lrFfvbt2yfi4uKEv7+/8PDwEHfccYd45ZVXxKVLlyzqnT59WgwZMkR4enqKli1bimeeecYihFUIIf7973+LPn36NFrf6zz++OOiV69e0vdycnJEly5dhI+Pj/Dy8hIxMTEiPT1dVFdXS/tuMBjs1v8jR46IESNGiKCgIOHl5SViY2MV4d3nz58XiYmJwtfXV/j6+orExERRUlLi8HNfX9/V2vjt311az31j9b0pX/PWXDNCNN1rfsaMGSI4OFi4ubmJdu3aicWLFyvC2pvqNV9f3xvimjcYDCIkJESkpqYy/L+J0gnBNJ9ERERE1+I9SUREREQSHCQRERERSXCQRERERCTBQRIRERGRBAdJRERERBIcJBERERFJcJBEREREJMFBEhEREZEEB0lEREREEhwkEREREUlwkEREjeKXX35BSEgI0tLSzGXffvst3N3dkZub68CeERHJ8dltRNRoNm/ejOHDh2PHjh2488470bVrVwwZMgTp6emO7hoRkQIHSUTUqJKSkvDFF1+gR48e+P7777Fnzx54eHg4ultERAocJBFRo7p8+TJiYmKQn5+PvXv3IjY21tFdIiKS4j1JRNSoTpw4gbNnz8JkMuH06dOO7g4RkSrOJBFRozEajbj77rvRpUsX3HnnnViyZAkOHjyI4OBgR3eNiEiBgyQiajTPP/881q9fj++//x4+Pj7o378/fH198c9//tPRXSMiUuByGxE1iu3btyM9PR3vv/8+/Pz8oNfr8f777+Prr7/GsmXLHN09IiIFziQRERERSXAmiYiIiEiCgyQiIiIiCQ6SiIiIiCQ4SCIiIiKS4CCJiIiISIKDJCIiIiIJDpKIiIiIJDhIIiIiIpLgIImIiIhIgoMkIiIiIgkOkoiIiIgkOEgiIiIikvh/lYR76gNy0ZEAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 640x480 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"ds[variables[1]][0].plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Converting a variables and its band into dataframe of xy coordinates and its raster values"
]
},
{
"cell_type": "code",
"execution_count": 80,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[nan, nan, nan, ..., nan, nan, nan],\n",
" [nan, nan, nan, ..., nan, nan, nan],\n",
" [nan, nan, nan, ..., nan, nan, nan],\n",
" ...,\n",
" [nan, nan, nan, ..., 29., nan, nan],\n",
" [nan, nan, nan, ..., 34., nan, nan],\n",
" [nan, nan, nan, ..., 37., nan, nan]], dtype=float32)"
]
},
"execution_count": 80,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"raster_values = ds[variables[1]][0].values\n",
"raster_values"
]
},
{
"cell_type": "code",
"execution_count": 81,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>x</th>\n",
" <th>y</th>\n",
" <th>raster_value</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>217305.0</td>\n",
" <td>6765165.0</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>217335.0</td>\n",
" <td>6765165.0</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>217365.0</td>\n",
" <td>6765165.0</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>217395.0</td>\n",
" <td>6765165.0</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>217425.0</td>\n",
" <td>6765165.0</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>...</th>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3635</th>\n",
" <td>219255.0</td>\n",
" <td>6763635.0</td>\n",
" <td>29.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3636</th>\n",
" <td>219285.0</td>\n",
" <td>6763635.0</td>\n",
" <td>28.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3637</th>\n",
" <td>219315.0</td>\n",
" <td>6763635.0</td>\n",
" <td>37.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3638</th>\n",
" <td>219345.0</td>\n",
" <td>6763635.0</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3639</th>\n",
" <td>219375.0</td>\n",
" <td>6763635.0</td>\n",
" <td>NaN</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>3640 rows × 3 columns</p>\n",
"</div>"
],
"text/plain": [
" x y raster_value\n",
"0 217305.0 6765165.0 NaN\n",
"1 217335.0 6765165.0 NaN\n",
"2 217365.0 6765165.0 NaN\n",
"3 217395.0 6765165.0 NaN\n",
"4 217425.0 6765165.0 NaN\n",
"... ... ... ...\n",
"3635 219255.0 6763635.0 29.0\n",
"3636 219285.0 6763635.0 28.0\n",
"3637 219315.0 6763635.0 37.0\n",
"3638 219345.0 6763635.0 NaN\n",
"3639 219375.0 6763635.0 NaN\n",
"\n",
"[3640 rows x 3 columns]"
]
},
"execution_count": 81,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df = pd.DataFrame({\"x\": x.flatten(), \"y\": y.flatten(), \"raster_value\": raster_values.flatten()})\n",
"df"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Converting to GeoDataFrame"
]
},
{
"cell_type": "code",
"execution_count": 86,
"metadata": {},
"outputs": [],
"source": [
"\n",
"gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['x'], df['y']), crs=ds.rio.crs)"
]
},
{
"cell_type": "code",
"execution_count": 99,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/tmp/ipykernel_16680/4220489362.py:2: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.\n",
" ].to_file(\"shp/points.shp\")\n"
]
}
],
"source": [
"gdf[~gdf['raster_value'].apply(lambda x: math.isnan(x))].to_file(\"shp/points.shp\")"
]
},
{
"cell_type": "code",
"execution_count": 98,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>x</th>\n",
" <th>y</th>\n",
" <th>raster_value</th>\n",
" <th>geometry</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>217485.0</td>\n",
" <td>6765165.0</td>\n",
" <td>62.0</td>\n",
" <td>POINT (217485.000 6765165.000)</td>\n",
" </tr>\n",
" <tr>\n",
" <th>76</th>\n",
" <td>217485.0</td>\n",
" <td>6765135.0</td>\n",
" <td>63.0</td>\n",
" <td>POINT (217485.000 6765135.000)</td>\n",
" </tr>\n",
" <tr>\n",
" <th>77</th>\n",
" <td>217515.0</td>\n",
" <td>6765135.0</td>\n",
" <td>64.0</td>\n",
" <td>POINT (217515.000 6765135.000)</td>\n",
" </tr>\n",
" <tr>\n",
" <th>145</th>\n",
" <td>217455.0</td>\n",
" <td>6765105.0</td>\n",
" <td>66.0</td>\n",
" <td>POINT (217455.000 6765105.000)</td>\n",
" </tr>\n",
" <tr>\n",
" <th>146</th>\n",
" <td>217485.0</td>\n",
" <td>6765105.0</td>\n",
" <td>57.0</td>\n",
" <td>POINT (217485.000 6765105.000)</td>\n",
" </tr>\n",
" <tr>\n",
" <th>...</th>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3633</th>\n",
" <td>219195.0</td>\n",
" <td>6763635.0</td>\n",
" <td>31.0</td>\n",
" <td>POINT (219195.000 6763635.000)</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3634</th>\n",
" <td>219225.0</td>\n",
" <td>6763635.0</td>\n",
" <td>30.0</td>\n",
" <td>POINT (219225.000 6763635.000)</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3635</th>\n",
" <td>219255.0</td>\n",
" <td>6763635.0</td>\n",
" <td>29.0</td>\n",
" <td>POINT (219255.000 6763635.000)</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3636</th>\n",
" <td>219285.0</td>\n",
" <td>6763635.0</td>\n",
" <td>28.0</td>\n",
" <td>POINT (219285.000 6763635.000)</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3637</th>\n",
" <td>219315.0</td>\n",
" <td>6763635.0</td>\n",
" <td>37.0</td>\n",
" <td>POINT (219315.000 6763635.000)</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>2437 rows × 4 columns</p>\n",
"</div>"
],
"text/plain": [
" x y raster_value geometry\n",
"6 217485.0 6765165.0 62.0 POINT (217485.000 6765165.000)\n",
"76 217485.0 6765135.0 63.0 POINT (217485.000 6765135.000)\n",
"77 217515.0 6765135.0 64.0 POINT (217515.000 6765135.000)\n",
"145 217455.0 6765105.0 66.0 POINT (217455.000 6765105.000)\n",
"146 217485.0 6765105.0 57.0 POINT (217485.000 6765105.000)\n",
"... ... ... ... ...\n",
"3633 219195.0 6763635.0 31.0 POINT (219195.000 6763635.000)\n",
"3634 219225.0 6763635.0 30.0 POINT (219225.000 6763635.000)\n",
"3635 219255.0 6763635.0 29.0 POINT (219255.000 6763635.000)\n",
"3636 219285.0 6763635.0 28.0 POINT (219285.000 6763635.000)\n",
"3637 219315.0 6763635.0 37.0 POINT (219315.000 6763635.000)\n",
"\n",
"[2437 rows x 4 columns]"
]
},
"execution_count": 98,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gdf[~gdf['raster_value'].apply(lambda x: math.isnan(x))\n",
" ]"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "datascience-dev.guru",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment