{
"cells": [
{
"cell_type": "markdown",
"id": "3200aa72-5629-4976-86d0-f1f0793d0038",
"metadata": {},
"source": [
"Complex ICESat-2 Queries in the cloud\n",
"===============================\n",
"\n",
"Although the data referenced in here is stored as columnar parquet files, there is a [compatible raster rotation of the tessalation](https://ui.adsabs.harvard.edu/link_gateway/2007MNRAS.381..865C/doi:10.1111/j.1365-2966.2007.12297.x) that keeps the same indexing and nesting structure. Tiling 2 and 3 dimensional arrays at defined chunk intervals may be an elegent way for data fusion across gridded and non-gridded large environmental data sets. \n",
"\n",
"Reproducing this notebook\n",
"-------------------------\n",
"\n",
"The data in this notebook lives in the cloud, and is 'free' to access, in that it does not require payment to me. But it does require an AWS account, and amazon will charge you for compute. The data transfer cost will likely not be noticeable if you are in the US West-2 region where we are hosting it... but might be pricey if you cross regions or try and export it at volume outside of the cloud.\n",
"\n",
"I use an instance with 256GB of ram to run this notebook; the first data pull is 'small' at ~20GB in memory. The second pull is over 100GB of memory, and will crash the kernel and dask if you are on too small of a machine. Everything here is under development, but particularly the indexing scheme, which has already changed. The most recent version of [healpix morton indexing](https://doi.org/10.1016/j.heliyon.2017.e00332) is at the [mortie library](https://github.com/espg/mortie)... but it produces breaking changes will not work with the data reference in this notebook (i.e., the data already stored in S3), which is why the old functions are replicated instead of imported.\n",
"\n",
"One obvious way to run this as a smaller notebook would be to modify to only grab a single orbital period; as written, this grabs a two year time series for each of the basins.\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "ee405995-b870-4e25-881e-060def8c1109",
"metadata": {},
"outputs": [],
"source": [
"import vaex\n",
"\n",
"import numpy as np\n",
"from numpy import r_\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"\n",
"import mhealpy as mh\n",
"from tqdm.notebook import tqdm\n",
"\n",
"from numba import jit, int64, vectorize\n",
"\n",
"from dask_gateway import Gateway"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "17e31546-8fdd-4158-975c-3abbf5a430b3",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
""
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gateway = Gateway()\n",
"options = gateway.cluster_options()\n",
"options.worker_specification = '4CPU, 16GB'\n",
"cluster = gateway.new_cluster(options)\n",
"cluster.scale(32)\n",
"client = cluster.get_client(set_as_default=True)\n",
"client"
]
},
{
"cell_type": "markdown",
"id": "8ee649e7-6532-4334-8acc-b10ed91e2483",
"metadata": {
"tags": []
},
"source": [
"Tiling\n",
"--------\n",
"\n",
"The [mhealpy](https://mhealpy.readthedocs.io/) library is an superset of the [healpy](https://healpy.readthedocs.io/) library that adds support for multiresolution healpix maps. The astronomy community uses a different indexing scheme then the 'morton' ordering used in this notebook. There's actually three different schemes that are provided for: `ring`, `nested` and `uniq`, with the first two being fixed resolution, and the last being multiresolution.\n",
"\n",
"There are facilities for converting between to representations, and certain index methods are optimized for different types of plotting methods. We use the optimized healpix code, and then convert from nested format to morton indexing.\n",
"\n",
"A few things to note here--\n",
"\n",
"1. The data has already been tiled and uploaded; what is being shown is traversing the morton/healpix embedding space tree to retrieve that data. How to write the data shards is a post for another day.\n",
"2. All of the functions in the next cell are out of date; the mortie library has an updated number scheme that fixes several bugs, and marks southern hemisphere data as negative integers. However, the S3 shards were uploading using this previous scheme, so we're using old functions to index to that data structure.\n",
"3. The shard and chunk sizes add 2 levels to the tree. An updated data upload would add an addition tree level for the S3 partitions. The reason will be apparent later on; we're downloading quite a lot of data for fine scale queries. This would be an even more pronounced issue if this was ATL03 data.\n",
"\n",
"In other words, the next update to this will break the numbering scheme here, change tree depth, and be more amenable to ATL03 data."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "321d5b83-4a40-4ebb-820c-852549e47350",
"metadata": {},
"outputs": [],
"source": [
"# These are old and outdated\n",
"def unique2parent(unique):\n",
" '''\n",
" Assumes input is UNIQ\n",
" Currently only works on single resolution\n",
" Returns parent base cell\n",
" '''\n",
" nside = mh.uniq2nside(unique)\n",
" if nside.all():\n",
" order = int(np.log2(nside.max()))\n",
" unique = unique // 4**(order-1)\n",
" parent = (unique - 16) // 4\n",
" return parent\n",
"\n",
"def heal_norm(base, order, addr_nest):\n",
" N_pix = mh.hp.order2nside(order)**2\n",
" addr_norm = addr_nest - (base * N_pix)\n",
" return addr_norm\n",
"\n",
"@vectorize([int64(int64, int64, int64)])\n",
"def fastNorm2Mort(order, normed, parents):\n",
" # General version, for arbitary order\n",
" if order > 18:\n",
" raise ValueError(\"Max order is 18 (to output to 64-bit int).\")\n",
" mask = np.int64(3*4**(order-1))\n",
" res = np.zeros(order,dtype=np.int64)\n",
" num = 0\n",
" for j, i in enumerate(range(order, 0, -1)):\n",
" nextBit = (normed & mask) >> ((2*i) - 2)\n",
" num += (nextBit+1) * 10**(i-1)\n",
" mask = mask >> 2\n",
" if parents is not None:\n",
" parents = parents- 6\n",
" parents = parents * 10**(order)\n",
" num = num + parents\n",
" return num"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "e7bc29b9-64d0-489b-97ab-2ba25d99da03",
"metadata": {},
"outputs": [],
"source": [
"# Read in drainage basin polylines\n",
"# Source: http://icesat4.gsfc.nasa.gov/cryo_data/drainage_divides/Ant_Grounded_DrainageSystem_Polygons.txt\n",
"basins = pd.read_csv('./Ant_Grounded_DrainageSystem_Polygons.txt', \n",
" delim_whitespace=True, \n",
" names = ['Lat', 'Lon', 'basin'])"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "b2783518-4028-4a3e-b45a-f0c3ab73c6f7",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Lat
\n",
"
Lon
\n",
"
basin
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
-74.525273
\n",
"
-61.122545
\n",
"
1
\n",
"
\n",
"
\n",
"
1
\n",
"
-74.525435
\n",
"
-61.123664
\n",
"
1
\n",
"
\n",
"
\n",
"
2
\n",
"
-74.525467
\n",
"
-61.123826
\n",
"
1
\n",
"
\n",
"
\n",
"
3
\n",
"
-74.525576
\n",
"
-61.124567
\n",
"
1
\n",
"
\n",
"
\n",
"
4
\n",
"
-74.525609
\n",
"
-61.124498
\n",
"
1
\n",
"
\n",
"
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
\n",
"
\n",
"
1238996
\n",
"
-74.523392
\n",
"
-61.109736
\n",
"
27
\n",
"
\n",
"
\n",
"
1238997
\n",
"
-74.523938
\n",
"
-61.113444
\n",
"
27
\n",
"
\n",
"
\n",
"
1238998
\n",
"
-74.524484
\n",
"
-61.117151
\n",
"
27
\n",
"
\n",
"
\n",
"
1238999
\n",
"
-74.525030
\n",
"
-61.120859
\n",
"
27
\n",
"
\n",
"
\n",
"
1239000
\n",
"
-74.525273
\n",
"
-61.122545
\n",
"
27
\n",
"
\n",
" \n",
"
\n",
"
1239001 rows × 3 columns
\n",
"
"
],
"text/plain": [
" Lat Lon basin\n",
"0 -74.525273 -61.122545 1\n",
"1 -74.525435 -61.123664 1\n",
"2 -74.525467 -61.123826 1\n",
"3 -74.525576 -61.124567 1\n",
"4 -74.525609 -61.124498 1\n",
"... ... ... ...\n",
"1238996 -74.523392 -61.109736 27\n",
"1238997 -74.523938 -61.113444 27\n",
"1238998 -74.524484 -61.117151 27\n",
"1238999 -74.525030 -61.120859 27\n",
"1239000 -74.525273 -61.122545 27\n",
"\n",
"[1239001 rows x 3 columns]"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"basins"
]
},
{
"cell_type": "markdown",
"id": "fe933122-e532-4613-ab34-9b639f0bdcf2",
"metadata": {},
"source": [
"Querying Data\n",
"-------------\n",
"\n",
"\n",
"Healpix spatial queries occur seperate from the data; all that we need to do is define the grid cell resolution:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "ae91d141-253c-47f7-8900-341aca928874",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"6508.1853 km at tesselation order 0\n",
"3254.09265 km at tesselation order 1\n",
"1627.046325 km at tesselation order 2\n",
"813.5231625 km at tesselation order 3\n",
"406.76158125 km at tesselation order 4\n",
"203.380790625 km at tesselation order 5\n",
"101.6903953125 km at tesselation order 6\n",
"50.84519765625 km at tesselation order 7\n",
"25.422598828125 km at tesselation order 8\n",
"12.7112994140625 km at tesselation order 9\n",
"6.35564970703125 km at tesselation order 10\n",
"3.177824853515625 km at tesselation order 11\n",
"1.5889124267578125 km at tesselation order 12\n",
"0.7944562133789063 km at tesselation order 13\n",
"0.39722810668945313 km at tesselation order 14\n",
"0.19861405334472657 km at tesselation order 15\n",
"0.09930702667236328 km at tesselation order 16\n",
"0.04965351333618164 km at tesselation order 17\n",
"0.02482675666809082 km at tesselation order 18\n",
"0.01241337833404541 km at tesselation order 19\n"
]
}
],
"source": [
"def order2res(order):\n",
" res = 111* 58.6323*.5**order\n",
" return res\n",
"\n",
"for res in range(20):\n",
" print(str(order2res(res))+ ' km at tesselation order ' + str(res))"
]
},
{
"cell_type": "markdown",
"id": "6b3a37b1-16dc-4ccc-8440-72f0a2524b0f",
"metadata": {},
"source": [
"Requesting an 'order' of 11 will select grid cells that roughly 3km per side. However, the healpix query only guarantees that the observation is *somewhere* in that cell; we don't know if it is at the edge, or at center. Going up one order of resolution doesn't solve this problem for us-- if we want to execute a 'buffer', we need to calculate a neighborhood function. Fortunately, the base healpix library provides functions to do this, and we can cast the results to the morton numbering scheme. The neighborhood function will return all of the cells touching cell that the observation is present in, so calling it at order 11 will return a 9km area about the observation location, with a guarantee of at least 3km of buffer in every direction. This is quite fast:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "7509d8a0-07fb-4edd-9d2e-9ea0eb8dbe94",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 24.6 ms, sys: 0 ns, total: 24.6 ms\n",
"Wall time: 20.6 ms\n"
]
}
],
"source": [
"%%time\n",
"# Set order\n",
"order = 11\n",
"# Subset basin 2\n",
"b2 = basins[basins.basin == 2]\n",
"# Get all nieghbors\n",
"B2_11 = mh.hp.pixelfunc.get_all_neighbours(2**order, b2.Lon.values, \n",
" b2.Lat.values, nest=True, \n",
" lonlat=True)\n",
"# Drop repeats\n",
"Basin2_11 = np.unique(B2_11.ravel())"
]
},
{
"cell_type": "markdown",
"id": "0eedfb87-22db-4d7e-aa5d-2203ac7e3344",
"metadata": {},
"source": [
"Again, with healpix the *query* is separate from the data, and we can visualize what our cells will look like, even if we don't know *a priori* if those cells actually have anything inside them:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "040e3350-2a1a-4c2f-9556-86d375ab4f70",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAJQCAYAAAANCFe9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA2oklEQVR4nO3deXhU9b3H8c9JQkgggISwJSyCCCiLhEUKWEooEUGwQEEFqiJapbbXW29xL4IoFsXlohZLb6vUKyKgoBZQCATEyo4E2VSgIEtYA1cEWUxy7h+QMZP1N5M5ObO8X8+Th8yZM2e+CZk5n/ltx7Jt2xYAAEA5otwuAAAAhAZCAwAAMEJoAAAARggNAADACKEBAAAYITQAAAAjhAYAAGCE0AAAAIwQGgAAgBFCA1AJvvjiC91111264oorFB8fr/j4eF155ZW69957tWHDBrfLc4VlWZowYULYPycQTmLcLgAId9OnT9fvfvc7tWrVSv/5n/+pNm3ayLIs7dixQ7NmzVKXLl20a9cuXXHFFW6XWqlWr16tRo0ahf1zAuHE4toTgHM+++wz9ezZUzfeeKPeffddxcbGFttn7ty56tGjh5KTk12oEADM0T0BOOiZZ55RdHS0pk+fXmJgkKRhw4Z5BYZRo0YpISFBu3btUv/+/ZWQkKDGjRvrD3/4g86fP+/12BMnTui+++5TSkqKYmNj1bx5cz3++OPF9rMsS7/73e/0xhtvqFWrVoqPj1fnzp21Zs0a2batKVOmqFmzZkpISFDv3r21a9euYnW+/vrruuaaaxQXF6fExEQNHjxYO3bs8NrHl9pL6io4ePCg7rnnHjVu3FixsbFKTk7W0KFDdeTIkVJ/x6mpqfrpT39abHteXp5SUlI0ZMiQMp/z8OHDuvfee9WoUSPFxsaqWbNmevLJJ5Wbm+vZp0uXLrrxxhu9HteuXTtZlqX169d7ts2bN0+WZWnLli2l1guENBuAI3Jzc+34+Hi7W7duPj3ujjvusGNjY+2rrrrKfv755+2lS5faTzzxhG1Zlv3kk0969jt79qzdvn17u3r16vbzzz9vL1myxB43bpwdExNj9+/f3+uYkuymTZva3bt3t+fNm2fPnz/fbtmypZ2YmGg/8MAD9i9+8Qt7wYIF9syZM+369evb7du3t/Pz8z2Pf+aZZ2xJ9vDhw+2FCxfab775pt28eXO7Vq1a9tdff+1z7QU1jR8/3nP7wIEDdsOGDe2kpCT7xRdftJcuXWrPnj3bHj16tL1jx45Sf19Tp061JXnVYdu2vWjRIluS/eGHH5b6nIcOHbIbN25sN23a1J4+fbq9dOlS+6mnnrKrVq1qjxo1yrPfI488YickJNgXLlywbdu2Dx8+bEuy4+Pj7UmTJnn2+81vfmPXr1+/1FqBUEdoABxScGK59dZbi92Xm5tr//DDD56vwifoO+64w5Zkz5kzx+sx/fv3t1u1auW5/Ze//KXE/Z599llbkr1kyRLPNkl2gwYN7NOnT3u2vf/++7Yku0OHDl7P/9///d+2JPuLL76wbdu2T548acfHxxcLIvv27bOrVq1qjxgxwufaC2oqfAIfPXq0XaVKFXv79u3Ffl9lOX78uB0bG2s/9thjXttvvvlmu379+vYPP/xQ6nPee++9dkJCgv3NN994Pfb555+3Jdnbtm2zbdu2ly5dakuyV65cadu2bb/11lt2jRo17Pvuu89OS0vzPO7KK6/0+n0A4YbuCcAFnTp1UpUqVTxfL7zwgtf9lmVp4MCBXtvat2+vb775xnM7MzNT1atX19ChQ732GzVqlCRp2bJlXtvT0tJUvXp1z+2rrrpKktSvXz9ZllVse8FzrV69WmfPnvUct0Djxo3Vu3fvYs9jUntJPvroI6WlpXme31SdOnU0cOBA/eMf/1B+fr4k6eTJk/rggw90++23Kyam9PHeCxYsUFpampKTk5Wbm+v56tevnyTpk08+kST16NFDcXFxWrp0qSQpIyNDvXr10g033KBVq1bp+++/1/79+7Vz50716dPHp/qBUEJoABySlJSk+Pj4Ek+Wb7/9ttavX68PP/ywxMdWq1ZNcXFxXtuqVq2qc+fOeW7n5OSoQYMGXid8SapXr55iYmKUk5PjtT0xMdHrdsEYi9K2FzxXwXEaNmxYrM7k5ORiz2NSe0mOHTvm98yG0aNH6+DBg8rIyJAkzZo1S+fPny8WdIo6cuSI/vnPf3oFuCpVqqhNmzaSpOPHj0uS4uLi1KNHD09oWLZsmdLT09WrVy/l5eXp008/9Tw3oQHhjCmXgEOio6PVu3dvLVmyRIcOHfI66V599dWSpL179/p9/Dp16mjt2rWybdsrOBw9elS5ublKSkry+9hFn0eSDh06VOy+7OzsgD1P3bp1deDAAb8e27dvXyUnJ+uNN95Q37599cYbb6hr166e33NpkpKS1L59e02aNKnE+wsPUP35z3+uJ554QuvWrdOBAweUnp6uGjVqqEuXLsrIyFB2drZatmypxo0b+/UzAKGAlgbAQY8++qjy8vI0ZswY/fDDDwE99s9//nOdPn1a77//vtf2N99803N/IHTr1k3x8fF66623vLYfOHBAmZmZAXuefv36afny5frqq698fmx0dLRuu+02vf/++/r000+1YcMGjR49utzHDRgwQFu3btUVV1yhzp07F/sqHBr69Omj3NxcjRs3To0aNVLr1q0925cuXarMzExaGRD2CA2Ag3r06KE///nPWrRokTp27KhXXnlFmZmZWrFihWbNmqWxY8dKkmrWrOnzsW+//Xa1b99ed9xxh1566SUtXbpUEyZM0GOPPab+/fsH7AR22WWXady4cfrwww91++2366OPPtJbb72ltLQ0xcXFafz48QF5nokTJyopKUk9e/bU1KlTlZmZqXnz5umee+7Rl19+We7jR48erfPnz2vEiBGKj4/XLbfcYvScVapUUffu3fXaa68pMzNTixYt0rRp0zRgwACvlo9OnTqpdu3aWrJkidLT0z3b+/Tpo82bN+vIkSOEBoQ9uicAh40ZM0bdunXT1KlT9dJLLyk7O1uWZalRo0bq3r27li1bpt69e/t83Li4OC1fvlyPP/64pkyZomPHjiklJUVjx44N2Im8wKOPPqp69erp5Zdf1uzZsxUfH69evXrpmWee0ZVXXhmQ50hJSdG6des0fvx4TZ48WTk5Oapbt66uu+66YuMuStKyZUt1795dq1at0siRI1WrVq1yH9OwYUNt2LBBTz31lKZMmaIDBw6oRo0aatasmW644QbVrl3bs29UVJR69eql+fPne4WDbt26qXr16jp79qzS0tL8++GBEMGKkAAAwAjdEwAAwAihAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARQgMAADASUos7pUcNc7sEAADCUkb+3HL3oaUBAAAYITQAAAAjhAYAAGCE0AAAAIwQGgAAgBFCAwAAMEJoAAAARggNAADACKEBAAAYITQAAAAjhAYAAGCE0AAAAIwQGgAAgBFCAwAAMEJoAAAARggNAADACKEBAAAYITQAAAAjhAYAAGCE0AAAAIwQGgAAgBFCAwAAMEJoAAAARggNAADACKEBAAAYITQAAAAjhAYAAGCE0AAAAIwQGgAAgBFCAwAAMEJoAAAARggNAADACKEBAAAYITQAAAAjhAYAAGCE0AAAAIwQGgAAgBFCAwAAMEJoAAAARggNAADACKEBAAAYITQAAAAjhAYAAGCE0AAAAIwQGgAAgBFCAwAAMEJoAAAARggNAADACKEBAAAYITQAAAAjhAYAAGCE0AAAAIwQGgAAgBFCAwAAMEJoAAAARggNAADACKEBAAAYITQAAAAjhAYAAGCE0AAAAIwQGgAAgBFCAwAAMEJoAAAARggNAADASIzbBQAAIsvi7CzP932TO7hWB3xHSwMAoNIUDgwIPYQGAEClIDCEPronAACOIiyED1oaAACOKRoY+iZ38BrHQKAILYQGAIAjSgoMCG10TwAAAs50hgRBIrQQGgAAAWPS3UCXROgiNAAAAoLuiPDHmAYAQIX5EhgYCBm6aGkAAFQIKzxGDloaAAB+IzBEFkIDAKDCCAyRgdAAAPBLRcYjEDJCE6EBAFAhBIDIQWgAAPisorMemDURmggNAACfMPgxchEaAADGAhEYCB2hi3UaAADlojsBEqEBAFCC8kICrQyRidAAAPBw8voRtFaEPsY0AAAqHa0MoYmWBgCAJGe7DuiWCA+0NAAAACOEBgCAF1oCUBpCAwCg0gYpEkhCG2MaACCCOTlbAuGH0AAAEaakVgXCAkwQGgAgApTV/eB0YGB9hvBBaACAMEfLAgKF0AAAEaQywwLjJcIPsycAIEJw0kZFERoAIIgVfFpfnJ1V6vcmj69stDKEJ8u2bdvtIkylRw1zuwQAcIyTJ/jKPmmzbHToycifW+4+tDQAAAKKwBC+GAgJAC4rqSm/rBPv4uwsz7aC75khgcpA9wQAVLLSuiHC4STvbytD4SAEd5h0T9DSAACVKJKa7ssLAiz6FHoIDQBQCThBeuP3EZoIDQDgsEiZfmgSBMraJ5JaYUIVoQEAHBTJJ8KSBnCWdF9JQYIxDsGJKZcA4IDCCzBJkRcYTJm2PCA4EBoAwGGREhhMBz36OqsCwYPuCQAIoEgZv+CLsgJDSb+for9DuiqCB6EBABzASc7/VoLyFreCe+ieAIAAoSm9bL6c/AkKwYnQAAAIaoSx4EFoAIAAoDn9Iqd+doJDcGBMAwBUACez8gUqSDAg0n2EBgDwE60LJQvU76LwlTwRHOieAAA/cCKrPAXhgWDmPkIDAPiItRjcQVBzH6EBAHxAYKh8hX/nBAd3ERoAwBCBwR2Ff8/8zt1FaAAAP3DyQiRi9gQAlINZEsBFtDQAQBkIDMCPCA0AgJDBQEh3ERoAAIARQgMAIOgV7hqitcE9hAYAKAUnJ8AboQEASlB0ACSDIIMLgc4dhAYAQMghxLmDdRoABAVWW0RZmPoaHGhpAOA6mpqB0EBLAwDXlBUWCu7rm9yhwp8yK9KKwada4EeEBgCuKO1EXnR7ebdLUzRslHSMkgIJrR7Bh/+T4EFoAFCpSjoBOPFp3uREU14goZUh+PB/4i5CAwDHLc7OKvFTfEkngKItDqW1QPjCn2NwcnIfLQzBx7Jt23a7CFPpUcPcLgGAjwL56b1okDB5zqL7FgQYf4+NykPLT+XKyJ9b7j6EBgABU9YnQ97w4S8CXeUwCQ10TwCoMMICEBkIDQD8Uto4hQKEBVQUCzoFH0IDAGP+Tn8EEB4IDQDKxIwDuIFWhuBEaABQjK8tCLypwyn8bQUXQgMQYcobi1Ae3sSByMUFq4AI5G9LAoEBlYm/t+BDSwMQ5vwdk1B0ESTewFFZAjXAlnERgUdoACJUeW+ivMnCDczICW6EBiBM8eaLUBPIlgH+/p1BaADCUGVdSRIIFKdO8vzdBxahAQgzXOQHoSbQYw+4VoVzmD0BAAhLdFEEHi0NQJigSwKhyKlWBjiD0ACEMK4uiVDm9JRIXgOBR2gAQhRXl0Qoo0UgNBEagBBDWEC44W83dBAagBBBVwTChS/dEiazgWi1qDyEBiAEMMgR4cIkMJQXAkxDAlMvA4/QAAQxPkEhXPl7Iuc14S5CAxBkyntT5FMTwpk/oYAui8pDaACCAEEBkaakrgNeB8GPFSEBl5U3wJE3SoQz0xYBX14HXBLbObQ0AC5gJgTwI6dO8ryWAo/QAFQCml2Bi/omd3B0vAEzJpxF9wTgMAZkAWY40Qc/WhoQUZxeTdGJ/lkg3JQ0+JGLVYUGQgP8EmpLGZsuFlPWSG5fRnmXJFh/N4CbnAjsvNacQ2iAT0Ihyfv6ab/w/mU9lqAABCcCQ+UhNMBjcXZWmS84k5OmWwP+fD2hF64jEAOzCgcQ3rSAyhMKH2TCiWXbtu12EabSo4a5XUJYCpUV2PzpHuAEDoQ3WhkCJyN/brn7MHsiwgXyRO/0C7agVgIDAInA4Aa6JyKYv4GhrBenyQvXny6M8gIDbxhAZCEwuIPQEKHcvNSyP89T0rgD3iiAyERgcA+hIcKE8if1UKgRgHP44OA+xjSAFx4AwAgtDSEsEIMYCQwAQgFdEsGB0BCCAjXjgRcegGBHWAguhIYIw4sOAOAvQkOICOUBjADgD1oZgg8DIUMAy6QCiGQEhuBBS0OI4kUEIJwF8pLZCBxCQxBjTjKASETravAiNAQJXiQAwDiGYMeYBgAAYITQECTKS9QkbgCRhPe84ERoCBF0XwAIZ4uzsxj8GAIIDSGE4AAgHPHeFjoIDUGub3IHr9Rd8OIqnMoBIFzQyhDcCA1BpOiLpbQXT+GwQHgAEMqYLRFamHIZZEp70ZTU2lDY4uwsXnAAQhbvX6GBlgYAAGCEloYQVFIiZ+QxgFBD12roITSEib7JHbwGSZZ0PwAEK7pYQwPdEwAAVxASQg+hIYwUnZ4JAEAgERrCUEF4KG/GBQAAviA0hDlaHgCECj7cBD9CQwQgOAAIVoXfn3ivCn6EhgjC6pEAgIogNESIouMbCBAAggHvQ6GF0AAAAIwQGiJISVMySfkAAFOsCBmBSuqqKLodAICiaGmIcKzlAMBNBe9BfGgJDYQGAIBryrpmDoIPoQEAABghNIAuCgCuKPx+Q/dEaCA0wK+gQLgAECgEhtDB7Amob3KHckNASfcz6wJARZi89yC4EBrgpSAA+PJCXpydRXAA4DMCQ+ihewKSvMMCL2QATuN9JjTR0oByldeKUBA0aG0AgPBm2bZtu12EqfSoYW6XgFIwvgGAqZJaGXjvcF9G/txy96GlAQHFFCoACF+MaQAAVKqiHyj4gBE6CA0ICK6gCcAXBIXQRGiAowgOAErCe0NoYkwDAqqkJakZ5wCgLMy+Ch20NMAxJb0J8OkCQFEEhtBBaAAAAEbonoCjyuqu4NMFAIQWWhpQaYrOsGDJagAS7wWhhNAAAHAdLY+hgdCASldSiwOAyMK6LqGJ0ADXFA4PvGEAQPAjNCBoEByAyEKXROghNAAAACOEBriOTxsAEBoIDQgqdFEAQPAiNCAo0NoARB6uSxN6CA0IGsykACIHgSE0ERoQNAq/iRAcgMhAYAgthAYEDd48gMjAh4LQRWgAALiCDwqhh9CAoMIFrYDwxlVuQxuhAUGP4AAAwYHQgKBT0icQggMQ+ngdhz5CA4JSwcWsuBomEB6YYhkeCA0IerzBAEBwiHG7AMAXJQ2iKtoCQcgAAGcQGhCSyuqqWJydRXAAgghdi+GD7gkEPX/ecJiuCQQPQnz4IDQg6BUdEFnS/aUhOADuK/o65HUZuggNCFmFwwTBAQheRV+ftDyELsY0IGSU90ZT1uBIVqEDgIqjpQFhqbRwwFgHAPAfoQFhq7wuC8IDAPiG0ICwVl53BOEBqBys7hoeCA0Ie+XNvpAID0BlYmxR6CI0IGKYhgcAzih4/fE6C13MnkDEKe8qmixLDQAlo6UBUNmtEHRdAIHHayo00dIAFNI3uUOpb2Zc2hfwHyEhPNDSABRBIACAktHSAJTA13EPpT0GiATlvR4YJxQ+CA1AgBAkEKnK6tYzuR+hg9AAGPJncRrGQSASlPR6ICSEJ8Y0AH4omG3haxDgjRThpLSZRQTk8EVLA1BBRd8gTWZfAOGq8OXqGcsQfggNQICVFyK4TDdCmS/hl7/x8ENoABzGIDBUhqLjZ5z4lF/a3zHhIHIQGoBKUNKa+7Q4wFduh08CMAgNgMsWZ2cRHBwQqJkrwdAc78+Juuhsn0CHVP5mIxOhAahEXOXPGeX9Pivr9x2s/6+FWwgIqagIQgPggqLNvLyR+8/tE7U/63cE+nndwN9sZCI0AAg55Z2cy1sG3GT/kh5b3n4VDRBOnoQLt3Ixngb+smzbtt0uwlR61DC3SwAcwXx2c4zgrzjT0OBPOEPoysifW+4+rAgJBKHSVtqLdCWFK39W5sRFFf0b42808tA9ASAkcIIKnLKmAPtzHEQOQgMQBEqbVcFgM64e6jRfwwK/+8jGmAYgCFVGv33hQFLa925jrIfzSvsdc4XWyGMypoGWBiBIlbT6XkVHvZd2vJK+L/z8nDTCV2mrPPJ/jpLQ0gAEOZPm48pc3rcyTia0MACVj5YGIAyYzP03XYPAtOuhrONV9gmdwAAED1oagBBUtPugqECPfSiPk89HaAAqh0lLg8+h4Y033tCdd95Z4n25ubl66KGH9OKLL/pySGOEBsB9FVlZ0dfjEhiAyuNIaIiKitLtt9+u1157TfHx8Z7t33zzjW6++WZlZWXp/PnzvldrgNAABIdAtz4QGAD3ObIi5Ouvv653331XXbp00fbt2yVJ8+fPV2pqqo4cOaJPPvnE90oBhBSnV2EkMADByefQMGrUKK1du1a2bevaa6/V0KFDNXToUPXs2VNZWVn6yU9+4kSdAIJQWeGBFRyB8OP3QMjt27erc+fOOnfunLp27apVq1bJsqxA1+eF7gkgePlycSOnxkUA8J8jYxokacGCBbrjjjuUkJCgm266Sa+99prS0tI0c+ZM1atXz69iTRAagNDjS4sDgQFwjyNjGh566CHddNNN6tatmzZt2qRXXnlFH330kbZs2aIOHTowpgGAXwgMQPDzuaUhNjZWTz/9tB566CGv7YcPH9att96qzz77TD/88ENAiyxASwMAAM5wZEXI5cuXq0ePHsW2N2jQQJmZmRo/fryvhwQAACGAFSEBAIAzYxoAAEBkMgoN0dHRWrdu3cUHREUpOjq61K+YGK6BBQBAODI6wz/xxBNq1KiR53un12MAAADBhzENAACAMQ0AACBw/AoNe/fu1b333quWLVuqTp06atmype69917t2bMn0PUBAIAg4XNoyMrKUmpqqmbMmKGUlBRdf/31SklJ0YwZM5SamqqsrCwHygQAAG7zeUxDr169lJ2draVLl6pJkyae7d98843S09OVkpKi5cuXB7xQiTENAAA4xZExDevWrdOTTz7pFRgkqWnTppowYYLWrl3r6yEBAEAI8Dk01KpVS7Vq1Srxvssuu0w1a9ascFEAACD4+BwaRowYob/97W8l3vc///M/Gj58eIWLAgAAwcfn5Rs7duyod999V9dee62GDx+uBg0a6PDhw5o1a5aOHj2qYcOGad68eZ79hwwZEtCCAQCAO3weCBkVVXbjhGVZKjikZVnKy8vzv7oiGAgJAIAzHLs0NgAAiDw+h4af/exnTtQBAACCHMtIAwAAI4QGAABghNAAAACMEBoAAIARQgMAADBCaAAAAEZ8mnK5ePFizZ8/X1u3blVOTo4sy1JiYqLatm2rX/7yl0pPT3eqTgAA4DKjloYzZ87ohhtuUL9+/fT222/rwoULatq0qZo0aaILFy7o7bff9tz//fffO10zAABwgVFLwx//+EetXbtWM2fO1LBhwxQT4/2wvLw8zZ07V/fdd5/++Mc/6sUXX3SkWAAA4B6jloY5c+boueee0/Dhw4sFBkmKjo7WrbfeqmeffVazZ88OeJEAAMB9RqHhxIkTatmyZbn7tWzZUidOnKhwUQAAIPgYhYbWrVtr1qxZ5e43a9YstW7dusJFAQCA4GM0puHhhx/WiBEjtG/fPt15551q166dEhMTZVmWcnJytGXLFs2YMUMff/yx3n77badrBgAALjAKDbfeeqvy8vL00EMP6ZZbbpFlWV7327athg0b6s0339Qtt9ziSKEAAMBdlm3btunO+fn5WrNmjWedBkmqU6eO2rVrp65duyoqytm1otKjhjl6fAAAIlVG/txy9/EpNLiN0AAAgDNMQoNPK0JKF1sb/v3vf3utCNm8eXPHWxkAAIC7jM/0O3fu1PDhw1WzZk21atVK3bt3V7du3dSqVSvVrFlTI0eO1Ndff+1krQAAwEVGLQ2bNm1Sr169VLVqVf3qV79S+/btlZiYKOniGg5ffPGF5s+fr4ULF+qTTz7RNddc42jRAACg8hmNabj++ut17tw5LVy4UDVq1Chxn++++04DBgxQXFycFi9eHPBCJcY0AADglICNaVi9erXmzp1bamCQpBo1auiRRx7RzTffbF4hAAAIGUZjGmJiYnT+/Ply97tw4UKJ16YAAAChzyg0pKWlady4cTpw4ECp+xw8eFDjx49X7969A1YcAAAIHkbNAi+88IKuu+46tWjRQr179/YMhCy8jHRmZqbq1Kmj+fPnO10zAABwgfHiTsePH9dzzz2n999/X7t371bBwyzLUosWLTR48GCNHTtWSUlJjhXLQEgAAJzh2IqQ586d08mTJyVJtWvXVlxcnO/V+YHQAACAMxxZEVKS4uLi1LBhQ38eCgAAQpRPoWHz5s2Ki4tTq1atJF1cUnrmzJnaunWrGjVqpJEjR3oWfQIAAOHFKDScOHFC6enpysrKkiTddNNNmjNnjgYOHKglS5Z49ps8ebLWrVunlJQUR4oFAADuMZpyOWnSJO3cuVOTJk3Sq6++qk2bNum2227zzJr47rvvtHjxYuXl5empp55yumYAAOACo5aGDz74QBMnTtTvf/97SVKbNm3Uq1cvvfbaa+rVq5ckKT09XY8//rhefvllp2oFAAAuMmppyM7OVmpqqud2p06dJElt27b12q9du3Y6ePBgAMsDAADBwig01K1bV/v37/fc/uabbySp2AqR+/fvV926dQNYHgAACBZGoaFnz56aMGGCVqxYoc8//1y//e1v1a1bN02aNElHjx6VJB06dEjPPvusunTp4mjBAADAHUaLO+3evVtdu3b1LOjUpEkTrV69WjfddJO++OILpaSk6MCBA4qKitKaNWt0zTXXOFIsizsBAOCMgK4Iefz4cX3wwQeqUqWKBg0apJo1a+r//u//NHnyZG3evFnJycn67W9/q44dO1a48NIQGgAAcIZjy0i7hdAAAIAzTEKD0ZiGwo4cOVLm/Rs2bPD1kAAAIAT4HBo6dOigzMzMEu+bOnWqrrvuugoXBQAAgo/PoeHqq69W3759NWHCBM/lsb/99lsNGTJEDzzwgH79618HvEgAAOA+n0PD0qVL9dhjj+npp59Wnz59tGDBAnXo0EHLly/Xu+++q1deecWJOgEAgMv8HgiZmZmpAQMG6Pz582rTpo0+/PBDXX755QEuzxsDIQEAcIYjAyEl6dSpU/rzn/+sc+fOqV69etqzZ49Wrlzpz6EAAECI8Dk0bNq0SZ06ddKyZcs0Z84c7d69W4MGDdKdd96pu+66S+fOnXOiTgAA4DKfQ0P37t1Vo0YNbdy4UUOHDlW1atX0v//7v5o+fbreeecdXXvttU7UCQAAXOZzaBg1apRWr16tK664wmv73XffrTVr1ig3NzdgxQEAgOAR8BUhv//+e1WrVi2Qh/RgICQAAM4wGQgZ48+BL1y4oJMnT8qyLF122WWKjY313OdUYAAAAO4y7p7IycnRo48+qtatW6tatWpKTk5Ww4YNVa1aNbVu3VqPP/64cnJynKwVAAC4yKh7Ys+ePfrpT3+qY8eOKS0tTe3bt1diYqIk6cSJE9qyZYuWL1+uevXq6ZNPPlGzZs0cKZbuCQAAnBGw7omxY8eqdu3aWrVqlZo0aVLiPvv27dOAAQP04IMP6t133/WtUgAAEPSMuicyMzP11FNPlRoYJKlJkyZ68skntWzZsoAVBwAAgodRaMjNzVV8fHy5+8XHxzPlEgCAMGUUGrp27arJkyfrzJkzpe5z5swZTZ48Wd26dQtYcQAAIHgYjWl47rnnlJaWpubNm2vo0KFq166dEhMTZVmWcnJytGXLFs2bN0/ff/+9VqxY4XDJAADADUahoWPHjlq3bp3GjRunGTNm6OzZs173x8fHa+DAgXryySfVqlUrRwoFAADu8nlFyLy8PO3evduzJkOdOnV0xRVXKDo62pECC2PKJQAAznBkRcjo6Gi1bNnSr4IAAEDo8ik02LattWvXauvWrcrJyZFlWUpMTFTbtm3VtWtXWZblVJ0AAMBlxqHhnXfe0YMPPqjs7GwV7dGwLEvJycmaMmWKbr311oAXCQAA3Gc05XL27NkaMWKErr76as2cOVNbt25Vdna2srOztXXrVs2cOVNt27bVyJEjNXdu+X0iAAAg9BgNhExNTdW1116r6dOnl7nfPffco/Xr12vTpk0BK7AwBkICAOAMk4GQRi0NX375pUaMGFHufiNGjNCXX35pckgAABBijEJDYmKidu7cWe5+u3bt8lz9EgAAhBej0DBs2DA9/PDDmjNnjvLz84vdn5+fr7lz5+qRRx7RzTffHPAiAQCA+4zGNJw5c0aDBw/W0qVLVaNGDV111VVey0jv2LFDp0+fVp8+fTR//nxVq1bNkWIZ0wAAgDNMxjT4tCLkokWLNH/+fG3bts1rRch27dpp8ODBuuGGG/yv1gChAQAAZwQ8NLiN0AAAgDMCNnsCAADAODTs379fU6dO1fTp0/Xtt996tt11113q2rWrhg4dqk8//dSxQgEAgLuMuie2b9+u7t2769SpU5Kkyy+/XCtWrFBaWppOnDihFi1a6Ouvv9bZs2f12WefqUuXLo4US/cEAADOCFj3xIQJE1SvXj2tWrVK27ZtU4sWLfSLX/xCSUlJ2rt3r9avX69///vfatu2rZ566qkKFw4AAIKP0QWrVq1apWeffVY/+clPJEkvv/yyrr76as2ZM0e1atWSdHEWxdixY/Xwww87Vy0AAHCNUUvDyZMn1ahRI8/txo0bS5IaNGjgtV9ycrJOnDgRwPIAAECwMAoNjRo10saNGz23169fL0nFLkz1+eefe4ULAAAQPoy6J4YOHarx48fr1KlTqlmzpl566SX9+te/1sSJE9WoUSN17txZa9as0TPPPKORI0c6XTMAAHCB0eyJ7777TgMGDPBMqRw0aJDmzJmjMWPG6PXXX5dlWbJtW02bNtXatWtVr149R4pl9gQAAM4I+IqQu3fvVpUqVdSkSRPPto8//libN29WcnKyhgwZourVq/tXrQFCAwAAzmAZaQAAYIRlpAEAQMAQGgAAgBFCAwAAMEJoAAAARggNAADACKEBAAAYITQAAAAjhAYAAGCE0AAAAIwQGgAAgBFCAwAAMEJoAAAARggNAADACKEBAAAYITQAAAAjhAYAAGCE0AAAAIwENDRERUWpUaNG+stf/qLc3NxAHhoAALgsoKGhZ8+eatiwoe6//361bNkykIcGAAAuiwnkwVasWCFJOn36tFauXBnIQwMAAJc5MqYhISFB/fv3d+LQAADAJT6Hhm3btpV5/4cffuh3MQAAIHj5HBq6du2qGTNmFNuem5urBx54QIMHDw5EXQAAIMj4HBqGDh2q0aNHa9SoUTp79qwkae/everRo4emTZum5557LuBFAgAA9/k8EHLGjBn62c9+pv/4j//Qhg0bNGbMGI0bN061atXSypUr1bVrVyfqBAAALrNs27b9eWBWVpa6d++u8+fPq1OnTsrIyFCtWrUCXZ+X9Khhjh4fAIBIlZE/t9x9/Jo9sX//fv3mN79Rbm6u2rdvr02bNmnq1Kn+HAoAAIQIn0PDwoULlZqaqoMHD+qTTz7Rxo0b9fDDD2vixIlKT0/X0aNHnagTAAC4zOfuiejoaPXr109vvvmmEhMTPdszMjL0q1/9SjExMTp48GDAC5XongAAwCmOdE8888wzWrBggVdgkKT09HRlZWWxfDQAAGHK74GQpcnPz1dUlDMXz6SlAQAAZzg2ELLMAzoUGAAAgLuM1mno3bu3pk2bptatW6t3795l7mtZlpYtWxaQ4gAAQPAwCg2FezDy8/NlWZbRvgAAIHwEfEyDkxjTAACAM1wZ0wAAAMKTz9eekKS8vDzNmTNHy5cvV05OjurUqaO0tDQNGzZMMTF+HRIAAAQ5n7snjh8/rhtuuEGff/65YmJiVKdOHeXk5Cg3N1epqalavHixkpKSHCmW7gkAAJzhSPfEAw88oK+++kozZ87U2bNndejQIZ09e1ZvvfWWdu7cqQceeMCvYgEAQHDzuS/hn//8p55++mkNHz7csy06OlojRozQ0aNHNWHChEDWBwAAgoTPLQ22batNmzYl3te2bVumXAIAEKZ8Dg19+vTR0qVLS7wvIyNDvXr1qmhNAAAgCPncPTFu3DgNGTJEeXl5GjFihBo0aKDDhw9r5syZmjdvnubNm6cTJ0549i96YSsAABCafJ49UfjaEoVXhiw4TNHVIvPy8ipSnxdmTwAA4AyT2RM+tzQ88cQTZS4jDQAAwhPLSAMAgMCt05CamqpJkyZpx44dFS4KAACEJqPQkJaWpr/97W9q27atrr76ao0bN05ZWVkOlwYAAIKJT90T69ev13vvvad58+Zp165datasmX75y19q6NChuvbaa52sUxLdEwAAOMWke8LvMQ2bN2/We++9p/fee087duxQSkqKhgwZoqFDh+q6665zZLAkoQEAAGc4GhoK27FjhydAbN68WfXr19ehQ4cqethiCA0AADij0kJDYbt379a8efP04IMPBvKwkggNAAA4xZXQ4CRCAwAAzgjY4k69e/f2up2ZmelfRQAAIGQZhYamTZs6XQcAAAhyRqHhjTfecLoOAAAQ5Hy+NDYAAIhMhAYAAGCE0AAAAIwQGgAAgBFCAwAAMEJoAAAARggNAADAiFFoOHr0qPLy8ry2bdmyRYMHD1bDhg2VkpKioUOHavv27Y4UCQAA3GcUGho2bKiNGzd6bm/btk3du3fXsmXLlJqaqvbt22vJkiXq0aOHdu3a5VixAADAPUahoeg1rZ544gnVrl1bW7Zs0aJFi/TRRx9p8+bNql69up5++mlHCgUAAO7ya0zDihUr9Oijj3pdk6JZs2Z68MEHtWzZsoAVBwAAgodfoeHbb79Vu3btim1v3769jh49WuGiAABA8DG6YJUkffXVV4qJubh73bp1derUqWL7nDp1SvHx8YGrDgAABA3j0DBq1CjP97Zta82aNerfv7/XPl988YUaN24csOIAAEDw8PvS2A0aNCi2beXKlbr++usrXhUAAAg6ll10akQQS48a5nYJAACEpYz8ueXuw4qQAADASEBDw/79+7Vv375AHhIAAAQJ44GQJpo3by7btpWbmxvIwwIAgCAQ0NBw++23Kz8/P5CHBAAAQSKgoeHvf/97IA8HAACCCAMhAQCAkYCGhnPnzjEQEgCAMBXQ0LBw4UI1a9YskIcEAABBgu4JAABgxGgg5MSJE40Otn379goVAwAAgpfRMtJRUVGyLEsmK05blqW8vLyAFFcUy0gDAOCMgC0jnZSUpLvvvlvHjh0r84splwAAhC+j7onU1FR9/fXXqlOnTpn71axZMyBFAQCA4GPU0nDNNddo8+bN5e5XvXp1NWnSpMJFAQCA4GM0puH06dPKyclR06ZNK6OmUjGmAQAAZ5iMaTDqnkhISFBCQkKFCwIAAKGLdRoAAIARQgMAADBCaAAAAEYIDQAAwAihAQAAGCE0AAAAI4QGAABghNAAAACMGK0ICQAAQEsDAAAwQmgAAABGCA0AAMAIoQEAABghNAAAACOEBgAAYITQAAAAjBAaACgnJ0f16tXT3r173S6lmAULFig1NVX5+flulwJEPEIDAP3pT3/SwIEDdfnll0uS9u3bp4EDB6p69epKSkrS/fffrwsXLvh83M8//1zp6em67LLLVKdOHd1zzz06ffq01z7lPdeAAQNkWZbefvvtCv2MACqO0ABEuLNnz+rvf/+77r77bklSXl6ebrzxRp05c0b/+te/9M477+i9997TH/7wB5+Om52drT59+qhFixZau3atPv74Y23btk2jRo3y7GP6XHfeeadeeeWVCv+sACrIBhDR3nvvPTspKclze9GiRXZUVJR98OBBz7ZZs2bZVatWtb/99lvj406fPt2uV6+enZeX59m2adMmW5K9c+dOn55r7969tiR79+7dfv2MAAKDlgYgwq1cuVKdO3f23F69erXatm2r5ORkz7a+ffvq/Pnz2rhxo/Fxz58/r9jYWEVF/fg2Ex8fL0n617/+5dNzNW3aVPXq1dOnn37q+w8IIGAIDUCE27t3r9dJ+/Dhw6pfv77XPrVr11ZsbKwOHz5sfNzevXvr8OHDmjJlii5cuKCTJ0/qsccekyQdOnTI5+dKSUkJyoGaQCQhNAAR7uzZs4qLi/PaZllWsf1s2y5xuySNGTNGCQkJni9JatOmjf7xj3/ohRdeULVq1dSgQQM1b95c9evXV3R0tM/PFR8fr++//97nnw9A4BAagAiXlJSkkydPem43aNCg2Kf8kydP6ocffijWKlBg4sSJysrK8nwVGDFihA4fPqyDBw8qJydHEyZM0LFjx9SsWTOfn+vEiROqW7duRX5UABVEaAAiXGpqqrZv3+653a1bN23dutXThSBJS5YsUdWqVdWpU6cSj1GvXj21aNHC81VU/fr1lZCQoNmzZysuLk7p6ek+Pde5c+e0e/dupaamVvjnBeA/y7Zt2+0iALhny5Yt6tixo44eParatWsrLy9PHTp0UP369TVlyhSdOHFCo0aN0qBBg3ye9vjqq6+qe/fuSkhIUEZGhh588EFNnjxZ999/vyQZP9eKFSs0cOBAHTlyRNWqVQvozw/AHC0NQIRr166dOnfurDlz5kiSoqOjtXDhQsXFxalHjx66+eabNWjQID3//PNej7MsSzNmzCjz2OvWrVN6erratWunv/71r5o+fbonMPjyXLNmzdLIkSMJDIDLaGkAoEWLFmns2LHaunWr1xTJ0uzdu1dXXnmltm/friuvvNLR2o4dO6bWrVtrw4YNnrEQANwR43YBANzXv39/7dy5UwcPHlTjxo3L3f/jjz/WPffc43hgkKQ9e/Zo2rRpBAYgCNDSAAAAjDCmAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARQgMAADBCaAAAAEYIDQAAwAihAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARQgMAADBCaAAAAEYIDQAAwAihAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARQgMAADBCaAAAAEZi3C4AkeXcuXO6cOGC22UACKDY2FjFxcW5XQYqAaEBlebcuXOqFV9bF3TO7VIABFCDBg20Z88egkMEIDSg0ly4cEEXdE7XWQNUJSpWsqJkRVmSZUmX/rWioi7djrp0O9pzX+HtsuT1OFlRP34v72PalnWxI8768Tje2yS78L6Xjm1fut+2Cu136XG25/alfXTxaRVVsH+h+wrdvvg8unTsH/eRCt8u6V/L63ap+xYcJ0o/1lTCvuUdSyrl+6KPibJLPaYs+2INXj+bXeQ4drH7rKL7yZYVdfFIVpF9LMv22t+y7Ev/ffaPfzKe/Wyv21EqvK+tKMtWdFT+xe2Xbhf+uvjfW2S7Cr7Pv7RPwbZ8RVs/fh9lSdGyZV3ar4qVd/H5dPF2tJUvSxf/jbLyFS1d+vfS4wv2U/7FfQr2Vb6ipEvHyb90v62oS/vFKE/RBce59JjoSz9ftAqObyta9qXvdbFmSdGWFC3r0veWovTj18XbUYqSpTOnpWadvtGFCxcIDRGA0IBKF6MqirGqXAwNBSf8S+/wluf78kKDdzDwCg1F7is/NFiyoysYGgpOcH6GhtLDgn58XpPQUPB9RUNDaQGi6L7lhYZLjw9IaLBKDw0/bvcOAuWGhsK3ywkNPwaE0kND4dveocH+MQBYtqpYUcVCQ5RXaLCL/Ft6aIi27Isnd8tStKxLz3vxxF5FF0/4F0ND/qUgUBAUVGJoiC4jNESXEBqiC/6PEREYCAkAAIwQGgAAgBFCAwAAMEJoAAAARggNAADACKEBAAAYYcolKl2ufpB1cb7epX9/nKtneeYiXvzXyo+WZ46ddfExXus0eO6LUqG5eBf/tQtNo7RVfMqlZ5tkF97XnymXUoWmXEpevwbWaQjgOg12kSmXBbftS+s05BeacmlX0joN+WGxToM86zQgchAaUGls21ZCQoL+dXqBlOd2NQACJSEhQbZtl78jQh6hAZXGsiydPn1a+/fvV82aNd0uB0AAnDp1So0bN764UBvCHqEBla5mzZqEBgAIQQyEBAAARggNAADACKEBlaZq1aoaP368qlat6nYpAAKE13VksWyGvAIAAAO0NAAAACOEBgAAYITQAAAAjBAaAACAEUIDKsW0adPUrFkzxcXFqVOnTvr000/dLglABa1cuVIDBw5UcnKyLMvS+++/73ZJcBihAY6bPXu2fv/73+vxxx/Xpk2b9NOf/lT9+vXTvn373C4NQAWcOXNG11xzjV599VW3S0ElYcolHNe1a1d17NhRr732mmfbVVddpUGDBulPf/qTi5UBCBTLsjR//nwNGjTI7VLgIFoa4KgLFy5o48aNuv766722X3/99Vq1apVLVQEA/EFogKOOHz+uvLw81a9f32t7/fr1dfjwYZeqAgD4g9CASlH0srm2bXMpXQAIMYQGOCopKUnR0dHFWhWOHj1arPUBABDcCA1wVGxsrDp16qSMjAyv7RkZGerevbtLVQEA/BHjdgEIf//1X/+l2267TZ07d1a3bt3017/+Vfv27dOYMWPcLg1ABZw+fVq7du3y3N6zZ4+ysrKUmJioJk2auFgZnMKUS1SKadOm6bnnntOhQ4fUtm1bvfTSS+rZs6fbZQGogBUrVigtLa3Y9jvuuEMzZsyo/ILgOEIDAAAwwpgGAABghNAAAACMEBoAAIARQgMAADBCaAAAAEYIDQAAwAihAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARQgMAADDy/8GNnDL3HjlxAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"nside = 2**11\n",
"m = np.zeros(mh.hp.nside2npix(nside))\n",
"m[Basin2_11.ravel()] += 1\n",
"mh.hp.gnomview(m, rot=[0,-90], xsize=800, ysize=800, nest=True)"
]
},
{
"cell_type": "markdown",
"id": "839174a4-5b8f-46a6-9b7a-6add3ade60a0",
"metadata": {},
"source": [
"That's our 3km buffer for basin 2.\n",
"\n",
"Our parquet data is sharded and stored at a different resolution; we can down-sample our result to figure out the S3 paths that we need to download to get the data:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "9da557c3-3c82-49e7-8399-4248bcb8481e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 2.28 ms, sys: 0 ns, total: 2.28 ms\n",
"Wall time: 1.86 ms\n"
]
}
],
"source": [
"%%time\n",
"# Convert to Morton indexing\n",
"order = 11\n",
"uniq = mh.nest2uniq(2**order, Basin2_11)\n",
"parents = unique2parent(uniq)\n",
"normed = heal_norm(parents,order, uniq)\n",
"# Get 'full' resolution \n",
"mortonsB2_11 = fastNorm2Mort(order,normed.ravel(), parents.ravel())\n",
"# Degrade to shard resolution\n",
"Basin2_6 = mortonsB2_11 // 10**5\n",
"Basin2_6 = np.unique(Basin2_6)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "1bb3fc25-04cc-43d8-98db-e112312913ee",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(6106, 64)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len(mortonsB2_11), len(Basin2_6)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "cee4e1e3-4726-4269-9c2c-1605cf9de148",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"60\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAJQCAYAAAANCFe9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1A0lEQVR4nO3deXwU9f3H8fckARIIICFc4RJEQDkknCUoJZSIIFigQQWqIlql2p+/2uJdBEEsitbigaW/qtSfiICCWkC5AqLlRoJcKlCQI5yBn8gZk8zvD2BhSQLfTWays7uv5+PBg+zsd2c/e8zMe77znVnLtm1bAAAAlxEV7AIAAEBoIDQAAAAjhAYAAGCE0AAAAIwQGgAAgBFCAwAAMEJoAAAARggNAADACKEBAAAYITQApeDrr7/WPffco6uuukpxcXGKi4vT1Vdfrfvvv1+rV68OdnlBYVmWRo4cGfbPCYSTmGAXAIS7iRMn6ne/+52aNGmi//7v/1azZs1kWZY2b96sKVOmqF27dtq6dauuuuqqYJdaqpYtW6Y6deqE/XMC4cTitycA9/z73/9W586ddfPNN+uDDz5Q2bJlC7SZPn26OnXqpKSkpCBUCADmODwBuOi5555TdHS0Jk6cWGhgkKT+/fv7BYbBgwcrPj5eW7duVc+ePRUfH6+6devqj3/8o06fPu332MOHD+uBBx5Q7dq1VbZsWTVs2FBPPfVUgXaWZel3v/ud3n77bTVp0kRxcXFq27atli9fLtu2NW7cODVo0EDx8fHq2rWrtm7dWqDOt956S9ddd51iY2OVkJCgvn37avPmzX5tAqm9sEMFe/bs0X333ae6deuqbNmySkpKUnp6uvbv31/ke5ycnKwbbrihwPS8vDzVrl1b/fr1u+Rz7tu3T/fff7/q1KmjsmXLqkGDBnrmmWeUm5vra9OuXTvdfPPNfo9r0aKFLMvSqlWrfNNmzJghy7K0fv36IusFQpoNwBW5ubl2XFyc3bFjx4Aed9ddd9lly5a1r7nmGvvFF1+0FyxYYD/99NO2ZVn2M88842t38uRJu2XLlnaFChXsF1980Z43b549fPhwOyYmxu7Zs6ffPCXZ9evXt1NSUuwZM2bYM2fOtBs3bmwnJCTYDz/8sP3LX/7SnjVrlj158mS7Ro0adsuWLe38/Hzf45977jlbkj1gwAB79uzZ9jvvvGM3bNjQrly5sv3dd98FXPu5mkaMGOG7vXv3brtWrVp2YmKi/Ze//MVesGCBPXXqVHvIkCH25s2bi3y/xo8fb0vyq8O2bXvOnDm2JPuTTz4p8jn37t1r161b165fv749ceJEe8GCBfbo0aPtcuXK2YMHD/a1e/zxx+34+Hg7JyfHtm3b3rdvny3JjouLs8eMGeNr99vf/tauUaNGkbUCoY7QALjk3Ibl9ttvL3Bfbm6u/dNPP/n+XbiBvuuuu2xJ9rRp0/we07NnT7tJkya+23/7298Kbff888/bkux58+b5pkmya9asaR87dsw37aOPPrIl2a1atfJ7/r/+9a+2JPvrr7+2bdu2jxw5YsfFxRUIIjt37rTLlStnDxw4MODaz9V04QZ8yJAhdpkyZexNmzYVeL8u5dChQ3bZsmXtJ5980m/6rbfeateoUcP+6aefinzO+++/346Pj7e///57v8e++OKLtiR748aNtm3b9oIFC2xJ9pIlS2zbtu13333Xrlixov3AAw/YqampvsddffXVfu8HEG44PAEEQZs2bVSmTBnfv5deesnvfsuy1Lt3b79pLVu21Pfff++7nZGRoQoVKig9Pd2v3eDBgyVJCxcu9JuempqqChUq+G5fc801kqQePXrIsqwC088917Jly3Ty5EnffM+pW7euunbtWuB5TGovzKeffqrU1FTf85uqWrWqevfurX/+85/Kz8+XJB05ckQff/yx7rzzTsXEFD3ee9asWUpNTVVSUpJyc3N9/3r06CFJ+vzzzyVJnTp1UmxsrBYsWCBJmj9/vrp06aKbbrpJS5cu1YkTJ7Rr1y5t2bJF3bp1C6h+IJQQGgCXJCYmKi4urtCN5XvvvadVq1bpk08+KfSx5cuXV2xsrN+0cuXK6dSpU77b2dnZqlmzpt8GX5KqV6+umJgYZWdn+01PSEjwu31ujEVR088917n51KpVq0CdSUlJBZ7HpPbCHDx4sNhnNgwZMkR79uzR/PnzJUlTpkzR6dOnCwSdi+3fv1//+te//AJcmTJl1KxZM0nSoUOHJEmxsbHq1KmTLzQsXLhQaWlp6tKli/Ly8vTFF1/4npvQgHDGKZeAS6Kjo9W1a1fNmzdPe/fu9dvoXnvttZKkHTt2FHv+VatW1YoVK2Tbtl9wOHDggHJzc5WYmFjseV/8PJK0d+/eAvdlZWU59jzVqlXT7t27i/XY7t27KykpSW+//ba6d++ut99+Wx06dPC9z0VJTExUy5YtNWbMmELvv3CA6i9+8Qs9/fTTWrlypXbv3q20tDRVrFhR7dq10/z585WVlaXGjRurbt26xXoNQCigpwFw0RNPPKG8vDwNHTpUP/30k6Pz/sUvfqFjx47po48+8pv+zjvv+O53QseOHRUXF6d3333Xb/ru3buVkZHh2PP06NFDixYt0rfffhvwY6Ojo3XHHXfoo48+0hdffKHVq1dryJAhl31cr169tGHDBl111VVq27ZtgX8XhoZu3bopNzdXw4cPV506ddS0aVPf9AULFigjI4NeBoQ9QgPgok6dOun111/XnDlz1Lp1a7366qvKyMjQ4sWLNWXKFA0bNkySVKlSpYDnfeedd6ply5a666679PLLL2vBggUaOXKknnzySfXs2dOxDdgVV1yh4cOH65NPPtGdd96pTz/9VO+++65SU1MVGxurESNGOPI8o0aNUmJiojp37qzx48crIyNDM2bM0H333advvvnmso8fMmSITp8+rYEDByouLk633Xab0XOWKVNGKSkpeuONN5SRkaE5c+ZowoQJ6tWrl1/PR5s2bVSlShXNmzdPaWlpvundunXTunXrtH//fkIDwh6HJwCXDR06VB07dtT48eP18ssvKysrS5ZlqU6dOkpJSdHChQvVtWvXgOcbGxurRYsW6amnntK4ceN08OBB1a5dW8OGDXNsQ37OE088oerVq+uVV17R1KlTFRcXpy5duui5557T1Vdf7chz1K5dWytXrtSIESM0duxYZWdnq1q1arr++usLjLsoTOPGjZWSkqKlS5dq0KBBqly58mUfU6tWLa1evVqjR4/WuHHjtHv3blWsWFENGjTQTTfdpCpVqvjaRkVFqUuXLpo5c6ZfOOjYsaMqVKigkydPKjU1tXgvHggRXBESAAAY4fAEAAAwQmgAAABGCA0AAMAIoQEAABghNAAAACOEBgAAYITQAAAAjITUxZ3SovoHuwQAAMLS/Pzpl21DTwMAADBCaAAAAEYIDQAAwAihAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARQgMAADBCaAAAAEYIDQAAwAihAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARQgMAADBCaAAAAEYIDQAAwAihAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARQgMAADBCaAAAAEYIDQAAwAihAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARQgMAADBCaAAAAEYIDQAAwAihAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARQgMAADBCaAAAAEYIDQAAwAihAQAQcuZmZQa7hIgUE+wCAAAoDMHAewgNAICgIhyEDkIDAKDUEBBCG6EBAOAqgkL4IDQAAFzhZljontTKtXmjaJw9AQAAjBAaAACAEUIDAAAwQmgAAABGCA0AAMAIoQEAABghNAAAACOEBgAAYITQAAAAjBAaAACAEUIDAAAwQmgAAIQUfncieAgNAADH8cuW4YnQAAAAjBAaAACAEUIDAAAwQmgAAABGCA0AAMAIoQEAABghNAAAHMXpluErJtgFAABCGyEhchAaAADGCAiRjdAAACjAq+GAS0gHF6EBACKYV8MBvInQAAARiLCA4iA0AECYIyDAKYQGAAgjBAS4idAAAGEgEsICgyCDj9AAACEmEgICvInQAAAeR0iAVxAaAMCDCArwIkIDAHgAIQGhgNAAAEFASEAoIjQAQCkiLCCUERoAwEWEBIQTQgMAOIyg4D1FfSZc+yEwhAYAcAhhwRv4HNxDaACAYmLjVHq6J7Xi/fYAQgMABIANV3DwvntDVLALAAAAoYHQAAAAjHB4AgBQLIWdecBhhPBGaAAAFCmcT0kM59fmFkIDAMCHDSkuhdAAABGIcIDiYCAkAEQYtwID4xnCH6EBABBx6GkpHkIDAAAwQmgAgAjCHjZKgtAAAACMEBoAAIARQgMAADBCaAAARBTGdRQfoQEAABghNAAAACOEBgBAiXE1yMhAaAAAAEYIDQBgiL1pRDp+5RIALoGgAJxHaACAQhAWwhOnW5YMoQEAREgoCd67yEFoABCx2NgBgWEgJICIRGAAAkdPA4CIEelBIdKP50f663cCoQFAWIv0oAA4idAAIOwQFApiLxtOIDQACBuEBX8EBTiN0ABEkLlZmWG5IYnksBCOnye8i9AARJgLN7ChvMGJtKAQyp+VF/D+OYPQAESwUO15iJTAEIqfDcIboQGIcOc2wKGwgQrnsBAK739hwvkzQUGEBgCSvHvYItw2Sl56b4FAERoAFOCF3gfCAuA9hAYARQrWmIdQDwwEBIQrQgOASyqtwxahGhQICN7HZ+QcfrAKgDG3NuwEBiA00NMAICBOjncItbBASECkIzQAHlHYBtTLG6mShIdQCwsAziA0AB7m1dMgLxRIeCAshBc+z8jDmAYgRMzNyvT0SvpytXm5dgBm6GkAQoyXex8K63UgLCCYvLaMhDpCAxDCLt4ge2UFSVCAF5X0e+mV5SuYCA1AGPHClRwBLyHAOovQAIQhLx/CABC6CA1AmGNPC4BTOHsCQMijNwUoHfQ0AAhJnKERXJH2nhNMzyA0AAgprLyB4CE0AAgJhAUg+AgNADyNsOA9kXZoAucRGgB4DkEBXsL38TxCAxCAovawWKk4g/cR8DZCA+AALqZUcrxv8CK+l/4IDYDD6I0wx3sChBZCA1BKCBP+IvV1A6GM0AAEWaSNRCcsAKGL0ACg1IRSYAilWktTJIVcvgMFERoAuCpUVryhUicQTIQGABGNsACYIzQAcIVXN8ZerQsIBYQGIADdk1pF1DHdcEBIQHHwvSkcoQFAWGKlDziP0AAgbBAU3EUvG6KCXQCA0hWOG9buSa3C8nUBXkNPAxCBLtzAsvcI+COAFo2eBiDCsYKECcIlJHoaXOXVhYyNRMmE4xkU9DwAZ7B+vDRCQwQq7Y0CC2FoIUAAKAqhAa4LZMNDwPCWcOxVAVB8hAZ4ihMbKIKHs+h5AHAOocFF7KUFR3Hec4KGmXPvE9/ryBIpnzfrgcsjNAAqeqXISqRwBGIgMhEaXMbKNbTx2RWNXgcg8nCdBgAlwtUYgchBaADgCIIDQhnfXzOEhlLAlxGRgl4HILwxpgHwmKI2uqE0doDgED5C6XtXXHxfzdHTAIQI9uIBBBs9DaWEsyjglMKCA98tuIHvFS5GaADCwMVBgpU9YIbeu8AQGoAwRG8EADcwpgGIEOxRASgpQkMpYqWNYOM7CJzH8hA4QkMp40sKIBRwOAuFITQAACIOO3DFQ2gAAABGCA1BQMIFAIQiQgPgMRxLBuBVhAYAgJ9wD6709hYfoSFI+NIC8KJwDwwoGUIDAAAwQmgIInobAKB0sd4tGX57Aq4xXTjpDgVCS6gu2wSGkiM0BFk4/2T23KxMo4XUKwtyuH4OQCCcWh5ZnsIToQE4K5CVJStEAJGI0AAUQ0n2xggcQOnzSo9mqCM0eACHKCIL7wfCXbiuz8DZEwAAwBChwSPCee+TvQ4gcnhxeQ/n9WtpIzQAAAAjhAYAAGCE0OAhdKEBgLNYrzqL0IBS4cXjnACc5bXlnMDgPEIDAAAwQmjwGJIxAMCrCA0eFK7BwWtdlwCc47XlO1zXo8FGaECp8tqKBUDJsVxHDkKDR5GSAQBeQ2gAAIQVdrrcQ2hAqaMrEwgfXlueCQzuIjR4GF9+AICXEBoAAIARQoPHhWtvg9e6NAEEzmvLcbiuL72E0AAAAIwQGkIA6RkA4AWEBgSN17o2AZjz2vLLzlXpIDSECBYIACgc68fSQ2hAUHltbwXA5bHcRq6YYBcAc15M06w8gMjjxLrIqXWHF9eL4YyeBgQdwQMAQgOhAQAAGCE0oEToGgQQKA5NhC5CAzyBQxQAAkFgCA5CA0rMqYWX4ACEP5bz0EZoAACEFHoZgofQAEfQ2wDgcli+Qx+hAQAAGCE0wDF0GQIoCmdMhAdCAzyHLkwAhSEwBB+hAY5ioQZwMXYEwgehAZ7ESgbAhdgh8QZCAxzHmRQAzmE5Di+EBriCvQIADH4MP4QGAABghNAA1zixd0DXJhCa6GUIT4QGeB7BAYhMBAbvITQAAAAjhAa4ijMpgMjjxPJKL4M3ERrgOhZ+IHIQ8MMboQEhg5UREBnY0fAuQgNKBSsBIPxxWCL8ERpQajgFEwhfLJuRgdCAkMPKCQhP9DJ4H6EBpYqVAhB+CPKRg9CAUsdhCiB8cOXHyEJoAAAEFYEhdMQEuwCguOhtKFy4rYD5nMPf3KzMsPvehit6GlDq2Ai4o3tSq7Bc8Ybja0JBrBdCA6EBpYoVAwCELkIDAM+jtyEysFPhfYQGlBpWCO6JhI1qJLxGwOsIDSgVBAYAJlhXeBuhAQhx7IEj3BAcvIvQANexAnBPpAWGSHu9gNcQGuAqAgOA4mDd4U2EBgAAYITQAIQouuoR7uht8B5CA1zDAg83EJYiC+sRbyE0AAAAI4QGuIK9A3ext41IwvrEOwgNcBwLONxGaAKCg9AAAPA8dka8gdAAR7Fgu4+9bADBQmgA4Li5WZm+f4BT+D4FH6EBjmGBdp/XexkKCwqEBziJ71JwERoAOOJyK3NW9kDoiwl2AQgPbBAiU6Cf+4XtS9pr0j2pFd87oJQRGgAYc3IjffG8vH7oBd4xNyuT70uQEBpQYuztlY5grSRL6/N1shcCgDsIDQAkeSv8FVULYQLn0NsQHIQGlIiXNjQoKNw+n3B7PSgZgkPpIzQAIYoNKIDSRmhAWPHaXoebAwfhPK99f5zE9wdOsGzbtoNdhKm0qP7BLgEwxko6NIVzcJC89b0M9/c61MzPn37ZNlzcCXCBl1bMgFexnIQeQgMAXIANGVA0QgMARBAOCaAkCA0AgKChZye0EBoAAIARQgPgMPac4HVeO0TBMhM6CA0AAMAIoQEAIhC9DSiOgEPD22+/XeR9ubm5+sMf/lCiggAg2NiAAYULODTcc889Gjx4sE6ePOk3/fvvv1enTp30+uuvO1YcEGrY2ADFx/LjfQGHhrfeeksffPCB2rVrp02bNkmSZs6cqeTkZO3fv1+ff/6540UCAIDgCzg0DB48WCtWrJBt22rfvr3S09OVnp6uzp07KzMzUz/72c/cqBMA4DCvjWuA9xVrIGSzZs00ffp05efna8aMGWrfvr1mzpypK664wuHyAACRhEMU3las0DBr1izdcMMNqlatmh588EGtWrVKN954ow4cOOB0fQAAwCMCDg2PPvqobrnlFnXs2FFr167Vq6++qk8//VTr169Xq1atGNOAiMUeEkKRFw9RsCx5V8Ch4a9//avGjh2rWbNmKSEhQZKUlpamzMxMNW7cWN26dXO8SAAobWy4gIICDg2LFi3So48+WmB6zZo1lZGRoccff9yRwoBQEiobGC/uVTohXF9XaeH9g6mAQ0OnTp2KnllUlEaPHl2iggA4r3tSq7DfMETCa4wkoRLEI01MsAsA4Aw2mGdc6n1gQwSUjFFoiI6O1rJly9S+fXtFRUXJsqwi21qWpdzcXMcKBHAewaBkAg0Uc7MyI+Y9757UilCFyzIKDU8//bTq1Knj+/tSoQGAcyJlg+UFF7/XbECDL5JCW6gwCg0jRozw/T1y5Ei3agEiAivB0MDnBBTEmAbAAWxgAHfQ2+Atxboi5I4dO3T//fercePGqlq1qho3bqz7779f27dvd7o+AADgEQH3NGRmZio1NVUnTpxQSkqK2rRpo3379mnSpEmaOnWqFi9erFatWrlQKoBQceGeIWMDQofpHj2faeSybNu2A3lAly5dlJWVpQULFqhevXq+6d9//73S0tJUu3ZtLVq0yPFCJSktqr8r8wUiiZsr/GCd7kj3NVBy8/OnX7ZNwD0NK1eu1JtvvukXGCSpfv36GjlypO69995AZwkgxJlstM+1YS8VCF0Bh4bKlSurcuXKhd53xRVXqFKlSiUuCkBoKM4ePuEBCF0BD4QcOHCg/vGPfxR63//8z/9owIABJS4KgLc5cclmLvsMhJ6Aexpat26tDz74QO3bt9eAAQNUs2ZN7du3T1OmTNGBAwfUv39/zZgxw9e+X79+jhYMIHjc2MjT8wCEjoAHQkZFXbpzwrIsnZulZVnKy8srfnUXYSAkUHLF2TiXZo+A1+sDwpUrAyHdOjMCgPcEY2NMzwPgXQGHhp///Odu1AHAI7yy1861HgDv4TLSACR5JywUht4HwBsIDUAE83JQKAy9D0BwERqACBNqQaEo4fI6gFBSrB+sAgAAkYfQAAAAjBAaAACAkYDGNMydO1czZ87Uhg0blJ2dLcuylJCQoObNm+tXv/qV0tLS3KoTAAAEmVFPw/Hjx3XTTTepR48eeu+995STk6P69eurXr16ysnJ0Xvvvee7/8SJE27XDAAAgsCop+FPf/qTVqxYocmTJ6t///6KifF/WF5enqZPn64HHnhAf/rTn/SXv/zFlWIBAEDwGPU0TJs2TS+88IIGDBhQIDBIUnR0tG6//XY9//zzmjp1quNFAgCA4DMKDYcPH1bjxo0v265x48Y6fPhwiYsCAADeYxQamjZtqilTply23ZQpU9S0adMSFwUAALzHaEzDY489poEDB2rnzp26++671aJFCyUkJMiyLGVnZ2v9+vWaNGmSPvvsM7333ntu1wwAAILAKDTcfvvtysvL06OPPqrbbrtNlmX53W/btmrVqqV33nlHt912myuFAgCA4LJs27ZNG+fn52v58uW+6zRIUtWqVdWiRQt16NBBUVHuXisqLaq/q/MHACBSzc+fftk2AYWGYCM0AADgDpPQEPCvXObn5+s///mP3xUhGzZs6HovAwAACC7jLf2WLVs0YMAAVapUSU2aNFFKSoo6duyoJk2aqFKlSho0aJC+++47N2sFAABBZNTTsHbtWnXp0kXlypXTr3/9a7Vs2VIJCQmSzlzD4euvv9bMmTM1e/Zsff7557ruuutcLRoAAJQ+ozENN954o06dOqXZs2erYsWKhbb58ccf1atXL8XGxmru3LmOFyoxpgEAALc4NqZh2bJlmj59epGBQZIqVqyoxx9/XLfeeqt5hQAAIGQYjWmIiYnR6dOnL9suJyen0N+mAAAAoc8oNKSmpmr48OHavXt3kW327NmjESNGqGvXro4VBwAAvMOoW+Cll17S9ddfr0aNGqlr166+gZAXXkY6IyNDVatW1cyZM92uGQAABIHxxZ0OHTqkF154QR999JG2bdumcw+zLEuNGjVS3759NWzYMCUmJrpWLAMhAQBwh2tXhDx16pSOHDkiSapSpYpiY2MDr64YCA0AALjDlStCSlJsbKxq1apVnIcCAIAQFVBoWLdunWJjY9WkSRNJZy4pPXnyZG3YsEF16tTRoEGDfBd9AgAA4cUoNBw+fFhpaWnKzMyUJN1yyy2aNm2aevfurXnz5vnajR07VitXrlTt2rVdKRYAAASP0SmXY8aM0ZYtWzRmzBi99tprWrt2re644w7fWRM//vij5s6dq7y8PI0ePdrtmgEAQBAY9TR8/PHHGjVqlH7/+99Lkpo1a6YuXbrojTfeUJcuXSRJaWlpeuqpp/TKK6+4VSsAAAgio56GrKwsJScn+263adNGktS8eXO/di1atNCePXscLA8AAHiFUWioVq2adu3a5bv9/fffS1KBK0Tu2rVL1apVc7A8AADgFUahoXPnzho5cqQWL16sr776Sg8++KA6duyoMWPG6MCBA5KkvXv36vnnn1e7du1cLRgAAASH0cWdtm3bpg4dOvgu6FSvXj0tW7ZMt9xyi77++mvVrl1bu3fvVlRUlJYvX67rrrvOlWK5uBMAAO5w9IqQhw4d0scff6wyZcqoT58+qlSpkv7v//5PY8eO1bp165SUlKQHH3xQrVu3LnHhRSE0AADgDtcuIx0shAYAANxhEhqMxjRcaP/+/Ze8f/Xq1YHOEgAAhICAQ0OrVq2UkZFR6H3jx4/X9ddfX+KiAACA9wQcGq699lp1795dI0eO9P089g8//KB+/frp4Ycf1m9+8xvHiwQAAMEXcGhYsGCBnnzyST377LPq1q2bZs2apVatWmnRokX64IMP9Oqrr7pRJwAACLJiD4TMyMhQr169dPr0aTVr1kyffPKJrrzySofL88dASAAA3OHKQEhJOnr0qF5//XWdOnVK1atX1/bt27VkyZLizAoAAISIgEPD2rVr1aZNGy1cuFDTpk3Ttm3b1KdPH91999265557dOrUKTfqBAAAQRZwaEhJSVHFihW1Zs0apaenq3z58vrf//1fTZw4Ue+//77at2/vRp0AACDIAg4NgwcP1rJly3TVVVf5Tb/33nu1fPly5ebmOlYcAADwDsevCHnixAmVL1/eyVn6MBASAAB3mAyEjCnOjHNycnTkyBFZlqUrrrhCZcuW9d3nVmAAAADBZXx4Ijs7W0888YSaNm2q8uXLKykpSbVq1VL58uXVtGlTPfXUU8rOznazVgAAEERGhye2b9+uG264QQcPHlRqaqpatmyphIQESdLhw4e1fv16LVq0SNWrV9fnn3+uBg0auFIshycAAHCHY4cnhg0bpipVqmjp0qWqV69eoW127typXr166ZFHHtEHH3wQWKUAAMDzjA5PZGRkaPTo0UUGBkmqV6+ennnmGS1cuNCx4gAAgHcYhYbc3FzFxcVdtl1cXBynXAIAEKaMQkOHDh00duxYHT9+vMg2x48f19ixY9WxY0fHigMAAN5hNKbhhRdeUGpqqho2bKj09HS1aNFCCQkJsixL2dnZWr9+vWbMmKETJ05o8eLFLpcMAACCwSg0tG7dWitXrtTw4cM1adIknTx50u/+uLg49e7dW88884yaNGniSqEAACC4Ar4iZF5enrZt2+a7JkPVqlV11VVXKTo62pUCL8QplwAAuMOVK0JGR0ercePGxSoIAACEroBCg23bWrFihTZs2KDs7GxZlqWEhAQ1b95cHTp0kGVZbtUJAACCzDg0vP/++3rkkUeUlZWli49oWJalpKQkjRs3TrfffrvjRQIAgOAzOuVy6tSpGjhwoK699lpNnjxZGzZsUFZWlrKysrRhwwZNnjxZzZs316BBgzR9+uWPiQAAgNBjNBAyOTlZ7du318SJEy/Z7r777tOqVau0du1axwq8EAMhAQBwh8lASKOehm+++UYDBw68bLuBAwfqm2++MZklAAAIMUahISEhQVu2bLlsu61bt/p+/RIAAIQXo9DQv39/PfbYY5o2bZry8/ML3J+fn6/p06fr8ccf16233up4kQAAIPiMxjQcP35cffv21YIFC1SxYkVdc801fpeR3rx5s44dO6Zu3bpp5syZKl++vCvFMqYBAAB3mIxpCOiKkHPmzNHMmTO1ceNGvytCtmjRQn379tVNN91U/GoNEBoAAHCH46Eh2AgNAAC4w7GzJwAAAIxDw65duzR+/HhNnDhRP/zwg2/aPffcow4dOig9PV1ffPGFa4UCAIDgMjo8sWnTJqWkpOjo0aOSpCuvvFKLFy9WamqqDh8+rEaNGum7777TyZMn9e9//1vt2rVzpVgOTwAA4A7HDk+MHDlS1atX19KlS7Vx40Y1atRIv/zlL5WYmKgdO3Zo1apV+s9//qPmzZtr9OjRJS4cAAB4j9EPVi1dulTPP/+8fvazn0mSXnnlFV177bWaNm2aKleuLOnMWRTDhg3TY4895l61AAAgaIx6Go4cOaI6der4btetW1eSVLNmTb92SUlJOnz4sIPlAQAArzAKDXXq1NGaNWt8t1etWiVJBX6Y6quvvvILFwAAIHwYHZ5IT0/XiBEjdPToUVWqVEkvv/yyfvOb32jUqFGqU6eO2rZtq+XLl+u5557ToEGD3K4ZAAAEgdHZEz/++KN69erlO6WyT58+mjZtmoYOHaq33npLlmXJtm3Vr19fK1asUPXq1V0plrMnAABwh+NXhNy2bZvKlCmjevXq+aZ99tlnWrdunZKSktSvXz9VqFCheNUaIDQAAOAOLiMNAACMcBlpAADgGEIDAAAwQmgAAABGCA0AAMAIoQEAABghNAAAACOEBgAAYITQAAAAjBAaAACAEUIDAAAwQmgAAABGCA0AAMAIoQEAABghNAAAACOEBgAAYITQAAAAjBAaAACAEUdDQ1RUlOrUqaO//e1vys3NdXLWAAAgyBwNDZ07d1atWrX00EMPqXHjxk7OGgAABFmMkzNbvHixJOnYsWNasmSJk7MGAABB5sqYhvj4ePXs2dONWQMAgCAJODRs3Ljxkvd/8sknxS4GAAB4V8ChoUOHDpo0aVKB6bm5uXr44YfVt29fJ+oCAAAeE3BoSE9P15AhQzR48GCdPHlSkrRjxw516tRJEyZM0AsvvOB4kQAAIPgCHgg5adIk/fznP9d//dd/afXq1Ro6dKiGDx+uypUra8mSJerQoYMbdQIAgCCzbNu2i/PAzMxMpaSk6PTp02rTpo3mz5+vypUrO12fn7So/q7OHwCASDU/f/pl2xTr7Ildu3bpt7/9rXJzc9WyZUutXbtW48ePL86sAABAiAg4NMyePVvJycnas2ePPv/8c61Zs0aPPfaYRo0apbS0NB04cMCNOgEAQJAFfHgiOjpaPXr00DvvvKOEhATf9Pnz5+vXv/61YmJitGfPHscLlTg8AQCAW1w5PPHcc89p1qxZfoFBktLS0pSZmcnlowEACFPFHghZlPz8fEVFufPjmfQ0AADgDtcGQl5yhi4FBgAAEFxG12no2rWrJkyYoKZNm6pr166XbGtZlhYuXOhIcQAAwDuMQsOFRzDy8/NlWZZRWwAAED4cH9PgJsY0AADgjqCMaQAAAOEp4N+ekKS8vDxNmzZNixYtUnZ2tqpWrarU1FT1799fMTHFmiUAAPC4gA9PHDp0SDfddJO++uorxcTEqGrVqsrOzlZubq6Sk5M1d+5cJSYmulIshycAAHCHK4cnHn74YX377beaPHmyTp48qb179+rkyZN69913tWXLFj388MPFKhYAAHhbwMcS/vWvf+nZZ5/VgAEDfNOio6M1cOBAHThwQCNHjnSyPgAA4BEB9zTYtq1mzZoVel/z5s055RIAgDAVcGjo1q2bFixYUOh98+fPV5cuXUpaEwAA8KCAD08MHz5c/fr1U15engYOHKiaNWtq3759mjx5smbMmKEZM2bo8OHDvvYX/7AVAAAITQGfPXHhb0tceGXIc7O5+GqReXl5JanPD2dPAADgDpOzJwLuaXj66acveRlpAAAQnriMNAAAcO46DcnJyRozZow2b95c4qIAAEBoMgoNqamp+sc//qHmzZvr2muv1fDhw5WZmelyaQAAwEsCOjyxatUqffjhh5oxY4a2bt2qBg0a6Fe/+pXS09PVvn17N+uUxOEJAADcYnJ4othjGtatW6cPP/xQH374oTZv3qzatWurX79+Sk9P1/XXX+/KYElCAwAA7nA1NFxo8+bNvgCxbt061ahRQ3v37i3pbAsgNAAA4I5SCw0X2rZtm2bMmKFHHnnEydlKIjQAAOCWoIQGNxEaAABwh2MXd+ratavf7YyMjOJVBAAAQpZRaKhfv77bdQAAAI8zCg1vv/2223UAAACPC/insQEAQGQiNAAAACOEBgAAYITQAAAAjBAaAACAEUIDAAAwQmgAAABGjELDgQMHlJeX5zdt/fr16tu3r2rVqqXatWsrPT1dmzZtcqVIAAAQfEahoVatWlqzZo3v9saNG5WSkqKFCxcqOTlZLVu21Lx589SpUydt3brVtWIBAEDwGIWGi3/T6umnn1aVKlW0fv16zZkzR59++qnWrVunChUq6Nlnn3WlUAAAEFzFGtOwePFiPfHEE36/SdGgQQM98sgjWrhwoWPFAQAA7yhWaPjhhx/UokWLAtNbtmypAwcOlLgoAADgPUY/WCVJ3377rWJizjSvVq2ajh49WqDN0aNHFRcX51x1AADAM4xDw+DBg31/27at5cuXq2fPnn5tvv76a9WtW9ex4gAAgHcU+6exa9asWWDakiVLdOONN5a8KgAA4DmWffGpER6WFtU/2CUAABCW5udPv2wbrggJAACMOBoadu3apZ07dzo5SwAA4BHGAyFNNGzYULZtKzc318nZAgAAD3A0NNx5553Kz893cpYAAMAjHA0Nb775ppOzAwAAHsJASAAAYMTR0HDq1CkGQgIAEKYcDQ2zZ89WgwYNnJwlAADwCA5PAAAAI0YDIUeNGmU0s02bNpWoGAAA4F1Gl5GOioqSZVkyueK0ZVnKy8tzpLiLcRlpAADc4dhlpBMTE3Xvvffq4MGDl/zHKZcAAIQvo8MTycnJ+u6771S1atVLtqtUqZIjRQEAAO8x6mm47rrrtG7dusu2q1ChgurVq1fiogAAgPcYjWk4duyYsrOzVb9+/dKoqUiMaQAAwB0mYxqMDk/Ex8crPj6+xAUBAIDQxXUaAACAEUIDAAAwQmgAAABGCA0AAMAIoQEAABghNAAAACOEBgAAYITQAAAAjBhdERIAAICeBgAAYITQAAAAjBAaAACAEUIDAAAwQmgAAABGCA0AAMAIoQEAABghNABQdna2qlevrh07dgS7lAJmzZql5ORk5efnB7sUIOIRGgDoz3/+s3r37q0rr7xSkrRz50717t1bFSpUUGJioh566CHl5OQEPN+vvvpKaWlpuuKKK1S1alXdd999OnbsmF+byz1Xr169ZFmW3nvvvRK9RgAlR2gAItzJkyf15ptv6t5775Uk5eXl6eabb9bx48f15Zdf6v3339eHH36oP/7xjwHNNysrS926dVOjRo20YsUKffbZZ9q4caMGDx7sa2P6XHfffbdeffXVEr9WACVkA4hoH374oZ2YmOi7PWfOHDsqKsres2ePb9qUKVPscuXK2T/88IPxfCdOnGhXr17dzsvL801bu3atLcnesmVLQM+1Y8cOW5K9bdu2Yr1GAM6gpwGIcEuWLFHbtm19t5ctW6bmzZsrKSnJN6179+46ffq01qxZYzzf06dPq2zZsoqKOr+aiYuLkyR9+eWXAT1X/fr1Vb16dX3xxReBv0AAjiE0ABFux44dfhvtffv2qUaNGn5tqlSporJly2rfvn3G8+3atav27duncePGKScnR0eOHNGTTz4pSdq7d2/Az1W7dm1PDtQEIgmhAYhwJ0+eVGxsrN80y7IKtLNtu9DpkjR06FDFx8f7/klSs2bN9M9//lMvvfSSypcvr5o1a6phw4aqUaOGoqOjA36uuLg4nThxIuDXB8A5hAYgwiUmJurIkSO+2zVr1iywl3/kyBH99NNPBXoFzhk1apQyMzN9/84ZOHCg9u3bpz179ig7O1sjR47UwYMH1aBBg4Cf6/Dhw6pWrVpJXiqAEiI0ABEuOTlZmzZt8t3u2LGjNmzY4DuEIEnz5s1TuXLl1KZNm0LnUb16dTVq1Mj372I1atRQfHy8pk6dqtjYWKWlpQX0XKdOndK2bduUnJxc4tcLoPgs27btYBcBIHjWr1+v1q1b68CBA6pSpYry8vLUqlUr1ahRQ+PGjdPhw4c1ePBg9enTJ+DTHl977TWlpKQoPj5e8+fP1yOPPKKxY8fqoYcekiTj51q8eLF69+6t/fv3q3z58o6+fgDm6GkAIlyLFi3Utm1bTZs2TZIUHR2t2bNnKzY2Vp06ddKtt96qPn366MUXX/R7nGVZmjRp0iXnvXLlSqWlpalFixb6+9//rokTJ/oCQyDPNWXKFA0aNIjAAAQZPQ0ANGfOHA0bNkwbNmzwO0WyKDt27NDVV1+tTZs26eqrr3a1toMHD6pp06ZavXq1bywEgOCICXYBAIKvZ8+e2rJli/bs2aO6detetv1nn32m++67z/XAIEnbt2/XhAkTCAyAB9DTAAAAjDCmAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARQgMAADBCaAAAAEYIDQAAwAihAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARQgMAADBCaAAAAEYIDQAAwAihAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARQgMAADBCaAAAAEZigl0AIsupU6eUk5MT7DIAOKhs2bKKjY0NdhkoBYQGlJpTp06pclwV5ehUsEsB4KCaNWtq+/btBIcIQGhAqcnJyVGOTul6q5fKRJWVrChZUZZkWdLZ/62oqLO3o87ejvbdd+F0WfJ7nKyo83/Lf562ZZ05EGedn4//NMm+sO3Zedtn77etC9qdfZztu322jc48raLOtb/gvgtun3kenZ33+TbShbcL+9/yu11k23PzidL5mgppe7l5SUX8ffFjouwi5ynLPlOD32uzL5qPXeA+6+J2smVFnZmTdVEby7L92luWffbjs89/ZXztbL/bUbqwra0oy1Z0VP6Z6WdvX/jvzMd70XSd+zv/bJtz0/IVbZ3/O8qSomXLOtuujJV35vl05na0lS9LZ/6PsvIVLZ39/+zjz7VT/pk259oqX1HS2fnkn73fVtTZdjHKU/S5+Zx9TPTZ1xetc/O3FS377N86U7OkaEuKlnX2b0tROv/vzO0oRcnS8WNSgzbfKycnh9AQAQgNKHUxKqMYq8yZ0HBug392DW/5/r5caPAPBn6h4aL7Lh8aLNnRJQwN5zZwxQwNRYcFnX9ek9Bw7u+ShoaiAsTFbS8XGs4+3pHQYBUdGs5P9w8Clw0NF96+TGg4HxCKDg0X3vYPDfb5AGDZKmNFFQgNUX6hwb7o/6JDQ7Rln9m4W5aiZZ193jMb9jI6s8E/ExryzwaBc0FBhYaG6EuEhuhCQkP0uc8YEYGBkAAAwAihAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARTrlEqcvVT7LOnK939v/z5+pZvnMRz/xv5UfLd46ddeYxftdp8N0XpQvOxTvzv33BaZS2Cp5y6Zsm2Re2Lc4pl1KJTrmU/N4GrtPg4HUa7ItOuTx32z57nYb8C065tEvpOg35YXGdBvmu04DIQWhAqbFtW/Hx8fry2CwpL9jVAHBKfHy8bNu+fEOEPEIDSo1lWTp27Jh27dqlSpUqBbscAA44evSo6tate+ZCbQh7hAaUukqVKhEaACAEMRASAAAYITQAAAAjhAaUmnLlymnEiBEqV65csEsB4BCW68hi2Qx5BQAABuhpAAAARggNAADACKEBAAAYITQAAAAjhAaUigkTJqhBgwaKjY1VmzZt9MUXXwS7JAAltGTJEvXu3VtJSUmyLEsfffRRsEuCywgNcN3UqVP1+9//Xk899ZTWrl2rG264QT169NDOnTuDXRqAEjh+/Liuu+46vfbaa8EuBaWEUy7hug4dOqh169Z64403fNOuueYa9enTR3/+85+DWBkAp1iWpZkzZ6pPnz7BLgUuoqcBrsrJydGaNWt04403+k2/8cYbtXTp0iBVBQAoDkIDXHXo0CHl5eWpRo0aftNr1Kihffv2BakqAEBxEBpQKi7+2VzbtvkpXQAIMYQGuCoxMVHR0dEFehUOHDhQoPcBAOBthAa4qmzZsmrTpo3mz5/vN33+/PlKSUkJUlUAgOKICXYBCH9/+MMfdMcdd6ht27bq2LGj/v73v2vnzp0aOnRosEsDUALHjh3T1q1bfbe3b9+uzMxMJSQkqF69ekGsDG7hlEuUigkTJuiFF17Q3r171bx5c7388svq3LlzsMsCUAKLFy9Wampqgel33XWXJk2aVPoFwXWEBgAAYIQxDQAAwAihAQAAGCE0AAAAI4QGAABghNAAAACMEBoAAIARQgMAADBCaAAAAEYIDQAAwAihAQAAGCE0AAAAI4QGAABg5P8BaWQPQ8BMLUUAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Approximate view of data in S3\n",
"nside = 2**6\n",
"b2idx = mh.hp.ang2pix(nside,b2.Lon.values, b2.Lat.values, lonlat=True) \n",
"m = np.zeros(mh.hp.nside2npix(nside))\n",
"m[np.unique(b2idx).ravel()] += 1\n",
"# Centered on the south pole\n",
"mh.hp.gnomview(m, rot=[0,-90], xsize=800, ysize=800)# , nest=True)\n",
"print(len(np.unique(b2idx).ravel()))"
]
},
{
"cell_type": "markdown",
"id": "662d5a56-c568-4c3e-b1c9-a9ccdf6e6a26",
"metadata": {},
"source": [
"We can see that the basin outline for 3km resolution is split over 6,106 'cells', but those cells are in turn split over 64 parquet files-- at most, the pole hole actually results in some empty shards. \n",
"\n",
"You might notice a discrepancy above-- there's 64 parquet files, but the plot above only shows 60 shards plotted. That's because when we degraded to the shard resolution, we did it off the 3km buffer, which crossed over to a few adjacent cells... 4 of which pushed to adjacent shards. By doing the buffer at the higher resolution we don't miss any potential points, but also don't over subscribe to unneeded shards if we did the buffer at the lower resolution.\n",
"\n",
"We can build the file list in a loop, and get the data in parallel on distributed workers:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "a34daada-0cf5-48ff-be81-d724dd0533d8",
"metadata": {},
"outputs": [],
"source": [
"cycles = ['4','5','6','7','8','9','10','11']\n",
"prefix = 's3://geostacks/icesat_2/shards='\n",
"mid = '/cycle=' \n",
"suffix = '/atl06.parquet'\n",
"\n",
"filelist = []\n",
"\n",
"for i in Basin2_6:\n",
" for j in cycles:\n",
" s3path = prefix + str(i) + mid + j + suffix\n",
" filelist.append(s3path)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "ef97fe63-b133-4193-9623-2c51f79fe0a3",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"512"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# 64 shards repeated for years / 8 cycles\n",
"len(filelist)"
]
},
{
"cell_type": "markdown",
"id": "f4881177-f0fe-4020-b914-6ec7851cf25b",
"metadata": {},
"source": [
"We could request all the data be returned using this:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "e5f0d168-1be7-41fe-98bd-5a77df355a88",
"metadata": {},
"outputs": [],
"source": [
"def get_thing_full(s3file):\n",
" return vaex.open(s3file).to_arrow_table()"
]
},
{
"cell_type": "markdown",
"id": "0d1f6423-6267-4e53-8a63-594a1e9c8426",
"metadata": {},
"source": [
"...but we'll go one step further and do the subsetting on the workers as well:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "9a604000-16ab-4eb6-9f3d-50b8230a1905",
"metadata": {},
"outputs": [],
"source": [
"def get_thing(s3file, cells):\n",
" df = vaex.open(s3file)\n",
" # Get to our query resolution\n",
" df['midx11'] = df.midx // 10**7\n",
" df['basin2'] = df.midx11.isin(cells)\n",
" df.select(df.basin2 == True)\n",
" tmp = df.to_arrow_table(selection=True, parallel=True)\n",
" return tmp"
]
},
{
"cell_type": "markdown",
"id": "42b1f3cd-463d-4a78-9cf9-0f8c0a054796",
"metadata": {},
"source": [
"The next cell will fail if the cluster workers aren't spun up yet, so it's good practice to check if they're online... and to resubmit the scale call if not. The failure is related to the semantics of `scatter`, which expects to distribute the data to workers at execution time. This is different than `client.submit`, which will be more than happy to submit jobs to non-existent workers; the jobs will just queue until workers come online without error. For the latter case, sometimes this will actually speed up worker allocation, as ad hoc workers tend to activate faster when they **know** there is work for them!"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "1b924d66-1fb5-41b1-8bf6-09150b8b3856",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"Future: ndarray\n",
" status: \n",
"\n",
"\n",
"finished,\n",
"\n",
"\n",
"\n",
" type: numpy.ndarray,\n",
"\n",
"\n",
" key: ndarray-56609e93e968cc406bde452fcdc64f7d"
],
"text/plain": [
""
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"client.scatter(mortonsB2_11, broadcast=True)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "37d0fb08-5775-4bca-bf15-3328d26e6019",
"metadata": {},
"outputs": [],
"source": [
"futures = []\n",
"for thing in filelist:\n",
" futures.append(client.submit(get_thing, thing, mortonsB2_11, retries=10))"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "3c80a9c7-d8bc-407b-a155-25fa63859f46",
"metadata": {},
"outputs": [],
"source": [
"# ~22GB return\n",
"quiver = list(filter(None, client.gather(futures, errors='skip', direct=True)))"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "dc102393-fd89-40f8-a632-1dfeb2a5f83e",
"metadata": {},
"outputs": [],
"source": [
"big_list = []\n",
"for arrow in quiver:\n",
" big_list.append(vaex.from_arrow_table(arrow))\n",
"\n",
"df = vaex.concat(big_list)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "b7c9b66a-3ad4-4499-9581-3583cc3464e9",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
#
id
lat
lon
slope_y
slope_x
slope_x_sigma
h_li
s_li
q_flag
s_fg
snr
h_rb
bsnow_conf
cloud_flg_asr
cloud_flg_atm
msw_flag
fbsnow_h
bsnow_od
layer_flag
bckgrd
e_bckgrd
n_fit_photons
w_surface_window_final
t_year
track
x
y
midx
shards
cycle
midx11
basin2
\n",
"\n",
"\n",
"
0
1523093
-86.14954851360491
33.92593579210944
-0.0395953431725502
0.00046188640408217907
0.0006559690809808671
2939.65185546875
0.012169521301984787
0
0
0.0
0.16256438195705414
6
2
0
4
56.94702911376953
0.07961465418338776
1
23297.232421875
0.0
482
3.0
2019-07-10 02:32:02.307540
gt1l
232881.35931420964
346225.37644100154
2111134242242111122
2111134
4
211113424224
True
\n",
"
1
1523094
-86.14939581927959
33.92457007097016
-0.06023023650050163
-0.006562550086528063
0.0006581993075087667
2939.58447265625
0.01336896326392889
0
0
0.0
0.1646553874015808
6
2
0
4
59.05659103393555
0.08079415559768677
1
25988.84375
0.0
465
3.0
2019-07-10 02:32:02.310355
gt1l
232882.34804547214
346244.66751049814
2111134242242111124
2111134
4
211113424224
True
\n",
"
2
1523095
-86.14924304204673
33.92320645695685
-0.06262587010860443
-0.003759956220164895
0.0007193252677097917
2939.484375
0.013828647322952747
0
0
0.0
0.18167147040367126
6
2
0
4
58.752349853515625
0.08070683479309082
1
27742.94921875
0.0
461
3.0
2019-07-10 02:32:02.313168
gt1l
232883.35374213205
346263.95772095054
2111134242242111142
2111134
4
211113424224
True
\n",
"
3
1523096
-86.14909011053149
33.92184672722354
-0.025539278984069824
0.008118819445371628
0.0006527432706207037
2939.517578125
0.012108918279409409
0
0
0.0
0.16191399097442627
6
2
0
4
56.649208068847656
0.07967527210712433
1
28792.603515625
0.0
475
3.0
2019-07-10 02:32:02.315974
gt1l
232884.39146603562
346283.2462720065
2111134242242111142
2111134
4
211113424224
True
\n",
"
4
1523097
-86.14893700260534
33.92049143181697
-0.03492480516433716
0.011287172324955463
0.0006484501645900309
2939.72021484375
0.012137381359934807
0
0
0.0
0.16041499376296997
6
2
0
4
54.551910400390625
0.07864658534526825
1
24003.869140625
0.0
487
3.0
2019-07-10 02:32:02.318772
gt1l
232885.46588380224
346302.53291968664
2111134242242111144
2111134
4
211113424224
True
\n",
"
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
\n",
"
116,343,239
1549043
-81.82003381242077
-49.04448300989648
4.835084837395698e-05
0.0017873976612463593
0.0004322136810515076
88.74310302734375
0.007835650816559792
0
0
0.0
0.10365600883960724
-1
2
1
1
29.979246139526367
0.028110453858971596
1
148184.5625
0.0
424
3.0
2021-06-11 21:24:54.182550
gt3r
-670302.7729052545
581772.2705912268
5111432111414434313
5111432
11
511143211141
True
\n",
"
116,343,240
1549044
-81.81986157860989
-49.04480402632936
0.0001920364156831056
-0.0016964105889201164
0.000489795405883342
88.74688720703125
0.007848388515412807
0
0
0.0
0.11242944747209549
-1
2
1
1
29.979246139526367
0.027972066774964333
1
146144.84375
0.0
426
3.0
2021-06-11 21:24:54.185357
gt3r
-670320.1922462292
581780.8045149728
5111432111414434314
5111432
11
511143211141
True
\n",
"
116,343,241
1549045
-81.81968929358989
-49.04512368181247
0.0005683625931851566
0.0004018144973088056
0.0005003189435228705
88.71788024902344
0.009416015818715096
0
0
0.0
0.11978613585233688
-1
2
1
1
29.979246139526367
0.027833567932248116
1
131252.15625
0.0
426
3.0
2021-06-11 21:24:54.188167
gt3r
-670337.6020975797
581789.3578408199
5111432111414434332
5111432
11
511143211141
True
\n",
"
116,343,242
1549046
-81.81951694599303
-49.04544167739985
0.0005571695510298014
0.001968930708244443
0.00044805006473325193
88.7217788696289
0.007712152320891619
0
0
0.0
0.10411503911018372
-1
2
1
1
29.979246139526367
0.027694901451468468
1
128720.234375
0.0
439
3.0
2021-06-11 21:24:54.190980
gt3r
-670355.0003577591
581797.9348786987
5111432111414434334
5111432
11
511143211141
True
\n",
"
116,343,243
1549047
-81.81934457364784
-49.045759009733075
-7.287902553798631e-05
-0.0009209429263137281
0.00045299253542907536
88.7392578125
0.00943058729171753
0
0
0.0
0.11246149241924286
-1
2
1
1
29.979246139526367
0.02755611203610897
1
135125.71875
0.0
458
3.0
2021-06-11 21:24:54.193796
gt3r
-670372.3940364019
581806.5212699741
5111432111414434334
5111432
11
511143211141
True
\n",
"\n",
"
"
],
"text/plain": [
"# id lat lon slope_y slope_x slope_x_sigma h_li s_li q_flag s_fg snr h_rb bsnow_conf cloud_flg_asr cloud_flg_atm msw_flag fbsnow_h bsnow_od layer_flag bckgrd e_bckgrd n_fit_photons w_surface_window_final t_year track x y midx shards cycle midx11 basin2\n",
"0 1523093 -86.14954851360491 33.92593579210944 -0.0395953431725502 0.00046188640408217907 0.0006559690809808671 2939.65185546875 0.012169521301984787 0 0 0.0 0.16256438195705414 6 2 0 4 56.94702911376953 0.07961465418338776 1 23297.232421875 0.0 482 3.0 2019-07-10 02:32:02.307540 gt1l 232881.35931420964 346225.37644100154 2111134242242111122 2111134 4 211113424224 True\n",
"1 1523094 -86.14939581927959 33.92457007097016 -0.06023023650050163 -0.006562550086528063 0.0006581993075087667 2939.58447265625 0.01336896326392889 0 0 0.0 0.1646553874015808 6 2 0 4 59.05659103393555 0.08079415559768677 1 25988.84375 0.0 465 3.0 2019-07-10 02:32:02.310355 gt1l 232882.34804547214 346244.66751049814 2111134242242111124 2111134 4 211113424224 True\n",
"2 1523095 -86.14924304204673 33.92320645695685 -0.06262587010860443 -0.003759956220164895 0.0007193252677097917 2939.484375 0.013828647322952747 0 0 0.0 0.18167147040367126 6 2 0 4 58.752349853515625 0.08070683479309082 1 27742.94921875 0.0 461 3.0 2019-07-10 02:32:02.313168 gt1l 232883.35374213205 346263.95772095054 2111134242242111142 2111134 4 211113424224 True\n",
"3 1523096 -86.14909011053149 33.92184672722354 -0.025539278984069824 0.008118819445371628 0.0006527432706207037 2939.517578125 0.012108918279409409 0 0 0.0 0.16191399097442627 6 2 0 4 56.649208068847656 0.07967527210712433 1 28792.603515625 0.0 475 3.0 2019-07-10 02:32:02.315974 gt1l 232884.39146603562 346283.2462720065 2111134242242111142 2111134 4 211113424224 True\n",
"4 1523097 -86.14893700260534 33.92049143181697 -0.03492480516433716 0.011287172324955463 0.0006484501645900309 2939.72021484375 0.012137381359934807 0 0 0.0 0.16041499376296997 6 2 0 4 54.551910400390625 0.07864658534526825 1 24003.869140625 0.0 487 3.0 2019-07-10 02:32:02.318772 gt1l 232885.46588380224 346302.53291968664 2111134242242111144 2111134 4 211113424224 True\n",
"... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...\n",
"116,343,239 1549043 -81.82003381242077 -49.04448300989648 4.835084837395698e-05 0.0017873976612463593 0.0004322136810515076 88.74310302734375 0.007835650816559792 0 0 0.0 0.10365600883960724 -1 2 1 1 29.979246139526367 0.028110453858971596 1 148184.5625 0.0 424 3.0 2021-06-11 21:24:54.182550 gt3r -670302.7729052545 581772.2705912268 5111432111414434313 5111432 11 511143211141 True\n",
"116,343,240 1549044 -81.81986157860989 -49.04480402632936 0.0001920364156831056 -0.0016964105889201164 0.000489795405883342 88.74688720703125 0.007848388515412807 0 0 0.0 0.11242944747209549 -1 2 1 1 29.979246139526367 0.027972066774964333 1 146144.84375 0.0 426 3.0 2021-06-11 21:24:54.185357 gt3r -670320.1922462292 581780.8045149728 5111432111414434314 5111432 11 511143211141 True\n",
"116,343,241 1549045 -81.81968929358989 -49.04512368181247 0.0005683625931851566 0.0004018144973088056 0.0005003189435228705 88.71788024902344 0.009416015818715096 0 0 0.0 0.11978613585233688 -1 2 1 1 29.979246139526367 0.027833567932248116 1 131252.15625 0.0 426 3.0 2021-06-11 21:24:54.188167 gt3r -670337.6020975797 581789.3578408199 5111432111414434332 5111432 11 511143211141 True\n",
"116,343,242 1549046 -81.81951694599303 -49.04544167739985 0.0005571695510298014 0.001968930708244443 0.00044805006473325193 88.7217788696289 0.007712152320891619 0 0 0.0 0.10411503911018372 -1 2 1 1 29.979246139526367 0.027694901451468468 1 128720.234375 0.0 439 3.0 2021-06-11 21:24:54.190980 gt3r -670355.0003577591 581797.9348786987 5111432111414434334 5111432 11 511143211141 True\n",
"116,343,243 1549047 -81.81934457364784 -49.045759009733075 -7.287902553798631e-05 -0.0009209429263137281 0.00045299253542907536 88.7392578125 0.00943058729171753 0 0 0.0 0.11246149241924286 -1 2 1 1 29.979246139526367 0.02755611203610897 1 135125.71875 0.0 458 3.0 2021-06-11 21:24:54.193796 gt3r -670372.3940364019 581806.5212699741 5111432111414434334 5111432 11 511143211141 True"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df"
]
},
{
"cell_type": "markdown",
"id": "6525894a-6ca6-463f-b9bc-1237263f12df",
"metadata": {},
"source": [
"We can also do the same for **all** of the basins, and do a 3km buffer for the entire continent:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "ff7c402c-2572-44d7-95c2-4eecf2b3e0e6",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 353 ms, sys: 29.4 ms, total: 382 ms\n",
"Wall time: 381 ms\n"
]
}
],
"source": [
"%%time\n",
"# Set order\n",
"order = 11\n",
"# Get all neighbors\n",
"B11 = mh.hp.pixelfunc.get_all_neighbours(2**order, basins.Lon.values, \n",
" basins.Lat.values, nest=True, \n",
" lonlat=True)\n",
"# Drop repeats\n",
"B_11 = np.unique(B11.ravel())\n",
"\n",
"uniq = mh.nest2uniq(2**order, B_11)\n",
"parents = unique2parent(uniq)\n",
"normed = heal_norm(parents,order, uniq)\n",
"# Get 'full' resolution \n",
"mortons11 = fastNorm2Mort(order,normed.ravel(), parents.ravel())\n",
"# Degrade to shard resolution\n",
"Basin11 = mortons11 // 10**5\n",
"Basin11 = np.unique(Basin11)\n",
"\n",
"cycles = ['4','5','6','7','8','9','10','11']\n",
"prefix = 's3://geostacks/icesat_2/shards='\n",
"mid = '/cycle=' \n",
"suffix = '/atl06.parquet'\n",
"\n",
"filelist = []\n",
"\n",
"for i in Basin11:\n",
" for j in cycles:\n",
" s3path = prefix + str(i) + mid + j + suffix\n",
" filelist.append(s3path)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "2e5cd2ed-a45c-4ed9-b120-2c754d9e4db8",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5432"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len(filelist)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "0cd24517-9a58-4322-be40-80e18932517a",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"Future: ndarray\n",
" status: \n",
"\n",
"\n",
"finished,\n",
"\n",
"\n",
"\n",
" type: numpy.ndarray,\n",
"\n",
"\n",
" key: ndarray-45290d95a6ffdfd8d6e31725309b1ea1"
],
"text/plain": [
""
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"client.scatter(mortons11, broadcast=True)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "ecf3e9ac-4399-4191-898d-9c51452544af",
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "e348f9974c954c908f602b0c2a05c522",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/10 [00:00, ?it/s]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 2min 59s, sys: 1min 2s, total: 4min 2s\n",
"Wall time: 24min 8s\n"
]
}
],
"source": [
"%%time\n",
"results = []\n",
"for thing in filelist:\n",
" results.append(client.submit(get_thing, thing, mortons11, retries=10))\n",
" \n",
"# extra logic for bigger dataset-- start download to host pre-emptively\n",
"lst = np.array_split(results, 10)\n",
"big_list = []\n",
"for thing in tqdm(lst):\n",
" futures = thing.tolist()\n",
" quiver = list(filter(None, client.gather(futures, errors='skip', direct=True)))\n",
" for future in futures:\n",
" future.cancel()\n",
" for arrow in quiver:\n",
" big_list.append(vaex.from_arrow_table(arrow))\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "50b1b8fe-2435-454f-92d0-45edf8502385",
"metadata": {},
"source": [
"Slightly different logic for the larger run-- we're chunking the download into 10 partitions, and asynchronously transferring the data back as it completes. This does three things:\n",
"\n",
"1. Avoids brittle behavior in dask, which likes to crash when you start to exceed a dozen or so data GB to transfer in one go\n",
"2. Speeds up the computation by allowing us to send results back parallel to additional computation occurring\n",
"3. Clears the results off of the workers, which frees memory-- this allows us to save money by sizing smaller workers\n",
"\n",
"Note that not only is the download from S3 data stores occurring in parallel, but also the filtering of the subsets is also occurring on the worker for more efficient horizontal scaling."
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "c1f7e0e7-9c25-45e7-ad12-9a2bd0b2303d",
"metadata": {},
"outputs": [],
"source": [
"df = vaex.concat(big_list)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "aa57dbf0-9601-43d3-b9ef-a6a1b7c2ea9b",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
#
id
lat
lon
slope_y
slope_x
slope_x_sigma
h_li
s_li
q_flag
s_fg
snr
h_rb
bsnow_conf
cloud_flg_asr
cloud_flg_atm
msw_flag
fbsnow_h
bsnow_od
layer_flag
bckgrd
e_bckgrd
n_fit_photons
w_surface_window_final
t_year
track
x
y
midx
shards
cycle
midx11
basin2
\n",
"\n",
"\n",
"
0
1523093
-86.14954851360491
33.92593579210944
-0.0395953431725502
0.00046188640408217907
0.0006559690809808671
2939.65185546875
0.012169521301984787
0
0
0.0
0.16256438195705414
6
2
0
4
56.94702911376953
0.07961465418338776
1
23297.232421875
0.0
482
3.0
2019-07-10 02:32:02.307540
gt1l
232881.35931420964
346225.37644100154
2111134242242111122
2111134
4
211113424224
True
\n",
"
1
1523094
-86.14939581927959
33.92457007097016
-0.06023023650050163
-0.006562550086528063
0.0006581993075087667
2939.58447265625
0.01336896326392889
0
0
0.0
0.1646553874015808
6
2
0
4
59.05659103393555
0.08079415559768677
1
25988.84375
0.0
465
3.0
2019-07-10 02:32:02.310355
gt1l
232882.34804547214
346244.66751049814
2111134242242111124
2111134
4
211113424224
True
\n",
"
2
1523095
-86.14924304204673
33.92320645695685
-0.06262587010860443
-0.003759956220164895
0.0007193252677097917
2939.484375
0.013828647322952747
0
0
0.0
0.18167147040367126
6
2
0
4
58.752349853515625
0.08070683479309082
1
27742.94921875
0.0
461
3.0
2019-07-10 02:32:02.313168
gt1l
232883.35374213205
346263.95772095054
2111134242242111142
2111134
4
211113424224
True
\n",
"
3
1523096
-86.14909011053149
33.92184672722354
-0.025539278984069824
0.008118819445371628
0.0006527432706207037
2939.517578125
0.012108918279409409
0
0
0.0
0.16191399097442627
6
2
0
4
56.649208068847656
0.07967527210712433
1
28792.603515625
0.0
475
3.0
2019-07-10 02:32:02.315974
gt1l
232884.39146603562
346283.2462720065
2111134242242111142
2111134
4
211113424224
True
\n",
"
4
1523097
-86.14893700260534
33.92049143181697
-0.03492480516433716
0.011287172324955463
0.0006484501645900309
2939.72021484375
0.012137381359934807
0
0
0.0
0.16041499376296997
6
2
0
4
54.551910400390625
0.07864658534526825
1
24003.869140625
0.0
487
3.0
2019-07-10 02:32:02.318772
gt1l
232885.46588380224
346302.53291968664
2111134242242111144
2111134
4
211113424224
True
\n",
"
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
\n",
"
502,554,965
1643664
-65.1595352894234
-64.14736624840631
3.4028234663852886e+38
0.4343118965625763
0.08810029178857803
814.0770263671875
2.443772315979004
1
2
0.0
5.0
127
5
1
3
29.979246139526367
-2.748199462890625
1
95096.1171875
0.0
32
30.0
2021-06-02 22:11:39.206742
gt3r
-2458716.895432818
1191377.9171924246
5134112212142221244
5134112
11
513411221214
True
\n",
"
502,554,966
1643665
-65.1593575757191
-64.14740984828781
3.4028234663852886e+38
-1.3402669429779053
0.09628021717071533
794.2894897460938
2.651156425476074
1
2
0.0
5.0
127
5
1
3
29.979246139526367
-2.748805522918701
1
94514.2890625
0.0
30
34.18221664428711
2021-06-02 22:11:39.209570
gt3r
-2458735.934009656
1191384.8320884425
5134112212142221422
5134112
11
513411221214
True
\n",
"
502,554,967
1643666
-65.15917990299222
-64.14745429112405
0.026572395116090775
0.13136883080005646
0.09426511079072952
778.3578491210938
2.0801053047180176
1
2
0.0
5.0
127
5
1
3
29.979246139526367
-2.749412775039673
1
98183.234375
0.0
25
30.0
2021-06-02 22:11:39.212404
gt3r
-2458754.985957636
1191391.7087620741
5134112212142221424
5134112
11
513411221214
True
\n",
"
502,554,968
1643667
-65.15900240514684
-64.14750247559053
3.4028234663852886e+38
0.4388608932495117
0.06979343295097351
782.5272216796875
1.838956356048584
1
2
0.0
5.0
127
5
1
3
29.979246139526367
-2.750018835067749
1
105618.1640625
0.0
39
30.0
2021-06-02 22:11:39.215232
gt3r
-2458774.0978894485
1191398.4161995882
5134112212142222313
5134112
11
513411221214
True
\n",
"
502,554,969
1643668
-65.15882489514628
-64.14755037468672
3.4028234663852886e+38
0.23369531333446503
0.07562870532274246
785.801025390625
2.3799521923065186
1
2
0.0
5.0
127
5
1
3
29.979246139526367
-2.7506251335144043
1
106303.5546875
0.0
39
30.0
2021-06-02 22:11:39.218061
gt3r
-2458793.2051524017
1191405.1364588912
5134112212142222331
5134112
11
513411221214
True
\n",
"\n",
"
"
],
"text/plain": [
"# id lat lon slope_y slope_x slope_x_sigma h_li s_li q_flag s_fg snr h_rb bsnow_conf cloud_flg_asr cloud_flg_atm msw_flag fbsnow_h bsnow_od layer_flag bckgrd e_bckgrd n_fit_photons w_surface_window_final t_year track x y midx shards cycle midx11 basin2\n",
"0 1523093 -86.14954851360491 33.92593579210944 -0.0395953431725502 0.00046188640408217907 0.0006559690809808671 2939.65185546875 0.012169521301984787 0 0 0.0 0.16256438195705414 6 2 0 4 56.94702911376953 0.07961465418338776 1 23297.232421875 0.0 482 3.0 2019-07-10 02:32:02.307540 gt1l 232881.35931420964 346225.37644100154 2111134242242111122 2111134 4 211113424224 True\n",
"1 1523094 -86.14939581927959 33.92457007097016 -0.06023023650050163 -0.006562550086528063 0.0006581993075087667 2939.58447265625 0.01336896326392889 0 0 0.0 0.1646553874015808 6 2 0 4 59.05659103393555 0.08079415559768677 1 25988.84375 0.0 465 3.0 2019-07-10 02:32:02.310355 gt1l 232882.34804547214 346244.66751049814 2111134242242111124 2111134 4 211113424224 True\n",
"2 1523095 -86.14924304204673 33.92320645695685 -0.06262587010860443 -0.003759956220164895 0.0007193252677097917 2939.484375 0.013828647322952747 0 0 0.0 0.18167147040367126 6 2 0 4 58.752349853515625 0.08070683479309082 1 27742.94921875 0.0 461 3.0 2019-07-10 02:32:02.313168 gt1l 232883.35374213205 346263.95772095054 2111134242242111142 2111134 4 211113424224 True\n",
"3 1523096 -86.14909011053149 33.92184672722354 -0.025539278984069824 0.008118819445371628 0.0006527432706207037 2939.517578125 0.012108918279409409 0 0 0.0 0.16191399097442627 6 2 0 4 56.649208068847656 0.07967527210712433 1 28792.603515625 0.0 475 3.0 2019-07-10 02:32:02.315974 gt1l 232884.39146603562 346283.2462720065 2111134242242111142 2111134 4 211113424224 True\n",
"4 1523097 -86.14893700260534 33.92049143181697 -0.03492480516433716 0.011287172324955463 0.0006484501645900309 2939.72021484375 0.012137381359934807 0 0 0.0 0.16041499376296997 6 2 0 4 54.551910400390625 0.07864658534526825 1 24003.869140625 0.0 487 3.0 2019-07-10 02:32:02.318772 gt1l 232885.46588380224 346302.53291968664 2111134242242111144 2111134 4 211113424224 True\n",
"... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...\n",
"502,554,965 1643664 -65.1595352894234 -64.14736624840631 3.4028234663852886e+38 0.4343118965625763 0.08810029178857803 814.0770263671875 2.443772315979004 1 2 0.0 5.0 127 5 1 3 29.979246139526367 -2.748199462890625 1 95096.1171875 0.0 32 30.0 2021-06-02 22:11:39.206742 gt3r -2458716.895432818 1191377.9171924246 5134112212142221244 5134112 11 513411221214 True\n",
"502,554,966 1643665 -65.1593575757191 -64.14740984828781 3.4028234663852886e+38 -1.3402669429779053 0.09628021717071533 794.2894897460938 2.651156425476074 1 2 0.0 5.0 127 5 1 3 29.979246139526367 -2.748805522918701 1 94514.2890625 0.0 30 34.18221664428711 2021-06-02 22:11:39.209570 gt3r -2458735.934009656 1191384.8320884425 5134112212142221422 5134112 11 513411221214 True\n",
"502,554,967 1643666 -65.15917990299222 -64.14745429112405 0.026572395116090775 0.13136883080005646 0.09426511079072952 778.3578491210938 2.0801053047180176 1 2 0.0 5.0 127 5 1 3 29.979246139526367 -2.749412775039673 1 98183.234375 0.0 25 30.0 2021-06-02 22:11:39.212404 gt3r -2458754.985957636 1191391.7087620741 5134112212142221424 5134112 11 513411221214 True\n",
"502,554,968 1643667 -65.15900240514684 -64.14750247559053 3.4028234663852886e+38 0.4388608932495117 0.06979343295097351 782.5272216796875 1.838956356048584 1 2 0.0 5.0 127 5 1 3 29.979246139526367 -2.750018835067749 1 105618.1640625 0.0 39 30.0 2021-06-02 22:11:39.215232 gt3r -2458774.0978894485 1191398.4161995882 5134112212142222313 5134112 11 513411221214 True\n",
"502,554,969 1643668 -65.15882489514628 -64.14755037468672 3.4028234663852886e+38 0.23369531333446503 0.07562870532274246 785.801025390625 2.3799521923065186 1 2 0.0 5.0 127 5 1 3 29.979246139526367 -2.7506251335144043 1 106303.5546875 0.0 39 30.0 2021-06-02 22:11:39.218061 gt3r -2458793.2051524017 1191405.1364588912 5134112212142222331 5134112 11 513411221214 True"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df"
]
},
{
"cell_type": "markdown",
"id": "1b08d25c-73da-492b-875b-e938a3e280bd",
"metadata": {},
"source": [
"It's difficult to describe exactly how computationally complicated the above query is. In plain language, we've taken all of the drainage basins borders and applied a ~3km buffer about them, and extracted that data over 8 orbital cycles-- two full years of ICESat-2 data. The original data files were *time ordered*, such that each file was a roughly north/south transect going to or from the pole. Because the drainage basins span the full continent and are complex polylines, it's hard to tell how many ICESat-2 data granules they touch... but a good estimate would be most (i.e., well over half) of them. There's ~4,000 hdf5 data files per 91-day orbit repeat that covers Antarctica, so for 8 orbit repetitions we can safely assume over 10,000 files having data in them related to this query, although I suspect it's likely even more files. Of course, you'd need a way to execute the metadata query to know which ones to open; or you'd need to inspect all 16,000 of them.\n",
"\n",
"Above we've subset a half billion data points from a set that was over 10 billion, culling through a few terabytes of data to retrieve the ~100GB that we actually care about. Run times vary, but seem to be between 15 and 30 minutes, although one run finished in under 10."
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "f52aa03e-126e-43f1-9d06-707314e088fa",
"metadata": {},
"outputs": [],
"source": [
"nside = 2**18\n",
"\n",
"def convertcoords(lon, lat):\n",
" healp = mh.hp.ang2pix(nside, lon, lat, nest=True, lonlat=True)\n",
" return healp\n",
"\n",
"# Convert locations to healpix\n",
"df['nest'] = df.apply(convertcoords, arguments=[df.lon, df.lat])"
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "1166d6ce-6584-4450-b63d-ef751bd07e59",
"metadata": {},
"outputs": [],
"source": [
"# Takes a bit of time to process the nest coordinates for the next few plots...\n",
"# ...it is a billion numbers to convert!\n",
"nside = 2**18\n",
"m = np.zeros(mh.hp.nside2npix(nside))\n",
"m[df['nest'].to_numpy().ravel()] += 1"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "a8579aeb-4661-4234-9797-34516e2ee23d",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAJQCAYAAAANCFe9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+JUlEQVR4nO3deXgUVaL+8beSgAmbEvawiSiIIBIBkc1LkEUQFxRccEMdl9F5rqOD+yCgwozLdV/nqjgOyCgIioCyL15BRAUEwQUEBMKigd8oCEKS+v0B3aaTTlLdXdW19PfzPDyku6urT1d11Xnr1KlThmmapgAAACqR5nYBAACAPxAaAACAJYQGAABgCaEBAABYQmgAAACWEBoAAIAlhAYAAGAJoQEAAFhCaAAAAJYQGoAk+PLLL3X99derZcuWysrKUlZWlk466STddNNN+uyzz9wunisMw9Do0aMD/5lAkGS4XQAg6F5++WX96U9/UuvWrXXbbbepbdu2MgxD69ev16RJk9S5c2dt2LBBLVu2dLuoSbVs2TI1adIk8J8JBInBvScA53z88cc666yzdO6552rKlCmqWrVqmWkmT56s7t27Kycnx4USAoB1nJ4AHDRu3Dilp6fr5ZdfjhoYJGno0KERgWH48OGqUaOGNmzYoIEDB6pGjRpq2rSp/vKXv+i3336LeO+ePXt0yy23qHHjxqpatapOOOEE3X///WWmMwxDf/rTnzR+/Hi1bt1aWVlZ6tSpkz755BOZpqnHHntMLVq0UI0aNdS7d29t2LChTDlfe+01nXbaacrMzFR2drYGDx6s9evXR0wTS9mjnSrYvn27brzxRjVt2lRVq1ZVTk6OhgwZol27dpW7jHNzc9WzZ88yzxcVFalx48a66KKLKvzMnTt36qabblKTJk1UtWpVtWjRQmPGjFFhYWF4ms6dO+vcc8+NeN+pp54qwzC0YsWK8HNTp06VYRhas2ZNueUFfM0E4IjCwkIzKyvL7Nq1a0zvu+aaa8yqVauabdq0MR9//HFz3rx55gMPPGAahmGOGTMmPN2BAwfM9u3bm9WrVzcff/xxc86cOebIkSPNjIwMc+DAgRHzlGQ2b97c7Natmzl16lRz2rRpZqtWrczs7Gzz9ttvNy+44AJzxowZ5sSJE80GDRqY7du3N4uLi8PvHzdunCnJvPzyy82ZM2eab7zxhnnCCSeYxx57rPntt9/GXPZQmUaNGhV+vG3bNrNRo0Zm3bp1zSeeeMKcN2+e+dZbb5nXXXeduX79+nKX19NPP21KiiiHaZrmrFmzTEnm9OnTy/3MHTt2mE2bNjWbN29uvvzyy+a8efPMhx56yDzmmGPM4cOHh6e75557zBo1apiHDh0yTdM0d+7caUoys7KyzLFjx4an++Mf/2g2aNCg3LICfkdoABwSqlguu+yyMq8VFhaahw8fDv8rWUFfc801piTz7bffjnjPwIEDzdatW4cfv/TSS1Gne+SRR0xJ5pw5c8LPSTIbNmxo7tu3L/zcu+++a0oyO3ToEPH5Tz31lCnJ/PLLL03TNM29e/eaWVlZZYLIDz/8YB5zzDHmsGHDYi57qEwlK/DrrrvOrFKlirlu3boyy6siP/30k1m1alXzvvvui3j+kksuMRs0aGAePny43M+86aabzBo1aphbtmyJeO/jjz9uSjK/+uor0zRNc968eaYkc8mSJaZpmuaECRPMmjVrmrfccouZl5cXft9JJ50UsTyAoOH0BOCCjh07qkqVKuF///M//xPxumEYOu+88yKea9++vbZs2RJ+vGDBAlWvXl1DhgyJmG748OGSpPnz50c8n5eXp+rVq4cft2nTRpI0YMAAGYZR5vnQZy1btkwHDhwIzzekadOm6t27d5nPsVL2aD744APl5eWFP9+qOnXq6LzzztM///lPFRcXS5L27t2r9957T1dffbUyMsrv7z1jxgzl5eUpJydHhYWF4X8DBgyQJC1evFiS1L17d2VmZmrevHmSpLlz56pXr14655xztHTpUv3666/aunWrvvvuO/Xp0yem8gN+QmgAHFK3bl1lZWVFrSzffPNNrVixQtOnT4/63mrVqikzMzPiuWOOOUYHDx4MPy4oKFDDhg0jKnxJql+/vjIyMlRQUBDxfHZ2dsTjUB+L8p4PfVZoPo0aNSpTzpycnDKfY6Xs0fz4449xX9lw3XXXafv27Zo7d64kadKkSfrtt9/KBJ3Sdu3apffffz8iwFWpUkVt27aVJP3000+SpMzMTHXv3j0cGubPn6++ffuqV69eKioq0kcffRT+bEIDgoxLLgGHpKenq3fv3pozZ4527NgRUemecsopkqTNmzfHPf86depo+fLlMk0zIjjs3r1bhYWFqlu3btzzLv05krRjx44yr+Xn59v2OfXq1dO2bdviem///v2Vk5Oj8ePHq3///ho/fry6dOkSXs7lqVu3rtq3b6+xY8dGfb1kB9Wzzz5bDzzwgD799FNt27ZNffv2Vc2aNdW5c2fNnTtX+fn5atWqlZo2bRrXdwD8gJYGwEH33nuvioqKdPPNN+vw4cO2zvvss8/Wvn379O6770Y8/8Ybb4Rft0PXrl2VlZWlCRMmRDy/bds2LViwwLbPGTBggBYuXKhvvvkm5vemp6frqquu0rvvvquPPvpIn332ma677rpK3zdo0CCtXbtWLVu2VKdOncr8Kxka+vTpo8LCQo0cOVJNmjTRySefHH5+3rx5WrBgAa0MCDxCA+Cg7t276/nnn9esWbN0+umn69lnn9WCBQu0aNEiTZo0SSNGjJAk1apVK+Z5X3311Wrfvr2uueYaPfnkk5o3b55Gjx6t++67TwMHDrStAjvuuOM0cuRITZ8+XVdffbU++OADTZgwQXl5ecrMzNSoUaNs+ZwHH3xQdevW1VlnnaWnn35aCxYs0NSpU3XjjTfq66+/rvT91113nX777TcNGzZMWVlZuvTSSy19ZpUqVdStWze9+OKLWrBggWbNmqUXXnhBgwYNimj56Nixo2rXrq05c+aob9++4ef79Omj1atXa9euXYQGBB6nJwCH3XzzzeratauefvppPfnkk8rPz5dhGGrSpIm6deum+fPnq3fv3jHPNzMzUwsXLtT999+vxx57TD/++KMaN26sESNG2FaRh9x7772qX7++nnnmGb311lvKyspSr169NG7cOJ100km2fEbjxo316aefatSoUfr73/+ugoIC1atXTz169CjT7yKaVq1aqVu3blq6dKmuuOIKHXvssZW+p1GjRvrss8/00EMP6bHHHtO2bdtUs2ZNtWjRQuecc45q164dnjYtLU29evXStGnTIsJB165dVb16dR04cEB5eXnxfXnAJxgREgAAWMLpCQAAYAmhAQAAWEJoAAAAlhAaAACAJYQGAABgCaEBAABYQmgAAACW+Gpwp75pQ90uAgAAgTS3eHKl09DSAACI2ez8VVH/RrARGgAAMeuf0yHq3wg2QgMAIC7RWhhodQg2QgMAIC7RWhhodQg2QgMAALCE0AAASAinJFIHoQEAkJDKTkkQKoKD0AAAKJcdFT79HIKD0AAAKBcVPkoiNAAAKjQ7f5UtLQ6cpvA/Xw0jDQBIPrtaG2i18D9aGgAAMYmnxaCi99AC4R+EBgBATGJpMQgFgoreQwuEfxAaAAAxiaVlgEAQLIQGAEBMCAKpi9AAAIgZ/RBSE6EBABCz/jkdbLsUE/7BJZcAgLhwmiL10NIAAAAsITQAAABLCA0AgKSK1g+CvhH+QGgAACRVqC9EyaBA/wh/IDQAAFxBUPAfQgMAIOk4HeFPhAYAQNLRyuBPhAYAAGAJoQEA4CpOVfgHoQEAELNEK3qunPAnQgMApLh4AkCiFT1BwZ8IDQCQwmbnr6ICh2WEBgBIMaGWBQIDYkVoAIAUEwoKBAbEitAAAJDEVQyoHKEBAFJcKCzQ8oDKEBoAIMX1z+lguZWB1ojURmgAAFhuZaA1IrURGgAAYbQkoCKEBgBAWEWnKmbnryrzGiEjtRAaACCFRav0yzsF0T+nQ5nXSj8mRAQboQEAUli0gBBvxc9gUcFHaAAARIi34icwBB+hAQDgC5z6cB+hAQDgODsqfFoy3EdoAAAkrLJQQIUfDIQGAEDMSoeEREIBpx38g9AAAIiZnS0HtEL4B6EBAOAqWhr8g9AAALAsngqe/g7BQWgAAFgWawXPgE/BQmgAADiGwBAshAYAAGAJoQEAAFhCaAAA+AJXWbiP0AAA8AX6R7iP0AAAcAQtA8FDaAAAOCKeyzPhbYQGAIBjYgkCoZBBePAuQgMAwBHlDezECJH+RWgAADiivMqfUOBfhAYAgGdxqsJbCA0AAM+iVcJbCA0AAM+ipcFbCA0AAM+ipcFbCA0AgMCgZcJZhAYAQGDQMuEsQgMAALCE0AAAACwhNAAAkoY+B/5GaAAAJE1FfQ5iCRSED3cQGgAAnhDLfSro8OgOQgMAwDNKhwTCgbcQGgAArguFBUKCtxEaAACuIyz4A6EBAOAJsXZupDNk8hEaAACeEGtrA60TyUdoAAD4Ai0L7iM0AAB8gZYF9xEaAACeRguDdxAaAACusRIIaGHwDkIDAMA1BAJ/ITQAAABLCA0AAMASQgMAALCE0AAAACwhNAAAAEsIDQAAxzHWQjAQGgAAjov10kpChjcRGgAAjom38i8ZMggQ3kFoAAA4xo7BmxgAyjsIDQAAx9BKECyEBgCAY5xoJSCIuIfQAACIyquVM6cr3ENoAABEReWM0ggNAADHebXVArEhNAAAHEerRTAQGgAAgCWEBgAAYAmhAQAAWEJoAAAAlhAaAABJwRUU/kdoAAAkBVdQ+B+hAQDgKlog/IPQAABwFS0Q/kFoAADEhRaC1ENoAADEhRaC1ENoAAAAlhAaAACAJYQGAABgCaEBAABYQmgAACRFoldbcLWG+wgNAJDCklkRJ3q1BVdruI/QAAApLNkVMa0F/kZoAIAUU7ri9lNrA9xFaACAFFO64qYih1WEBgAAYAmhAQAAWEJoAIAUQkdEJILQAAAphP4LSAShAQBSmJdaHkqWxUvlwu8M0zRNtwthVd+0oW4XAQCAQJpbPLnSaWhpAAAAlhAaAABxmZ2/ypXTCJy6cA+hAQBShN2Vbf+cDknrWFmy7KHPJDwkH30aAAAAfRoAAMFAq4I3EBoAAJ7H+BLeQGgAAACWEBoAAL7FaYvkIjQAAHyL0xbJRWgAACQdLQT+RGgAACQdLQT+RGgAAACWEBoAAHHhFEPqITQAAOKS7FMMhBT3ERoAAL5APwj3ERoAAIAlhAYACBia8eEUQgMABEz/nA4EBziC0AAAAcT5fziB0AAA8DxaTryB0AAAARHkipWWE28gNABAQCSjYg1yMEHlCA0AAMs44k9thAYA8KnZ+as48kdSZbhdAABAfDjqR7LR0gAAPlKyZYFWBiQboQEAfKRk6wItDUg2QgMAwFW0mPgHoQEA4CpaTPyD0AAAACwhNAAAAEsIDQAAwBJCAwAAsITQAAA+ELQrDIL2fVIFoQEAfCB0hUFQKluumPAnQgMA+AiVLdxEaAAAAJYQGgDAR4JyegL+RGgAAB/x8+kJAo//ERoAwAdKV7h2VMDJrsT9HHhwhGGapul2IazqmzbU7SIAABBIc4snVzoNLQ0AAMASQgMAeBR9AOA1hAYA8Cj6AMBrCA0A4BO0PMBthAYA8AlaHuA2QgMAALCE0AAAcBynVoKB0AAAcBynVoKB0AAAACwhNAAAAEsIDQAA29GHIZgIDQDgAeVVsn6tfOnDEEyEBgDwgPIqWT9Xvn4NPCgfoQEAfMRPFbGfAw+iIzQAgI9QEcNNhAYAAGAJoQEAPGJ2/ipfnX5A6slwuwAAgCM49QCvo6UBADyElgZ4GaEBADxidv4qR1obnAoiBJzUY5imabpdCKv6pg11uwgAAATS3OLJlU5DSwMAALCE0AAAACwhNABACqD/AexAaACAFNA/pwPBAQkjNABAimAcCCSK0AAAAUcLA+xCaAAAj3CqcqeFAXYhNACAR4Qqd1oG4FWEBgDwGFoG4FWEBgBIAbRewA6EBgBIAZz6gB0IDQCQQjj1gUQQGgAAMaG1InURGgAAMaG1InURGgAAgCWEBgAAYAmhAQACjP4HsBOhAQACjP4HsBOhAQACjJYG2InQAABJlsyKnJYG2InQAABJVllFbmeooKUBdiI0AECSVVaR29k6QEsD7ERoAIAki6Ui90NLgR/KCHsQGgDAQ0pXwH5oKfBDGWEPQgMAuMiPIQGpi9AAAC4iJMBPCA0AANvRzyGYCA0AANvRghJMhAYAAGAJoQEA4DhOVwQDoQEA4DhOVwQDoQEAXMZROPyC0AAALuMoHH5BaAAAJA2tKv5GaAAAJA2tKv5GaAAAuI4WCH8gNAAAXEcLhD8QGgAAgCWEBgAIMJr9YSdCAwAEQHnhgGZ/2InQAAAus6M1IMjhgNYS7yA0AIDLglzh24Hl4x2EBgBwGUfS8AtCAwC4LNlH0oQUxIvQAABJ4pXKOhRSvFIe+AehAQCSpLIWhWRX4vQVQKwIDQDgEVTisaGlJPkIDQAAXyJkJR+hAQBcwFFy4liGyUdoAAAXcJScmNn5q1iGLiA0AIBLOFKOH4HBHYQGAHAJFR/8htAAAC6itcEeLMfkIDQAgItKtzZQ+cWn5HJkGTqH0AAAHsIpi8QRxJxDaAAAD3G7gnP78+0W71UWQVsOdiE0AICHuN3SELT7UsS7PJO1HPy2nAkNAOARXqpA3A4vTol1GTu9HPy2nAkNAOARdOZzXqyVtFPrwa/rl9AAAC6anb+qTAXCaIfuKrk+ElkPFQUDv54GMkzTNN0uhFV904a6XQQAAOLm5UA4t3hypdPQ0gAAQJJ4NTBYRWgAAMAhfjv9UBlCAwB4VMkKx+uVj9fL55agDTRFaAAAjypZ4XihWdtKx75UEk8A8PtyIjQAgAd58YjU7xWe3WJZHl5cn/EgNACAB1FBB0v/nA6BCA6EBgDwOD9XNrH2y/Dzd43GrjEfvILQAAAeFYQKNNZ+GdGm8fNyCEJQKInQAAAeFapwglbxxMrK9/dzsPATQgMAwLdCYSEofQa8jtAAAB5D5Wed1y5LjYUf1zOhAQA8xm+VH1IHoQEAPMiPR6EIPt+HBjYsAEFEa0Pw+XEd+z40+HGhAwDgR74PDVbRIgEAQGJSJjSEWiQIDwD8hv2Wc2bnr2L5xsB3oSHRlVs6PPBjAYDU5tYYD36sfwzTNE23C2FV37Sh4b9n56+iPwMAICahusMrdYhXyiFJc4snVzqNb0MDAACwj5XQ4LvTEwCQqjitGiwVrUevrmNaGgAA8BknTmvQ0gAAgE3cPvov+fludeqnpQEAANDSAAAA7ENoAAAgRcV6WoPQAABAioq1MyWhAQACwO1OekgNMYeG8ePHl/taYWGh7rjjjoQKBACAX6RaWIs5NFx//fUaPny4Dhw4EPH8li1b1L17dz3//PO2FQ4AYA035XOHV4aATpaYQ8Nrr72mKVOmqHPnzlq3bp0kadq0acrNzdWuXbu0ePFi2wuZKDYiAKki1SoxJFfMoWH48OFavny5TNPUGWecoSFDhmjIkCE666yztGrVKp155plOlDMhJRM4AQJAKmGflxypspzjHtxp3bp16tSpkw4ePKguXbpo6dKlMgzD7vJFiGdwJ6fuaOalO5MBQDTspxALxwZ3mjFjhnr27Kl69erp1ltv1YoVK9SvXz/t3r07ntk5KrTBlLfhxJsO2RABeB37Kdgt5tBw11136fzzz1fXrl21cuVKPfvss/rggw+0Zs0adejQwZN9GuySKs1PAABEE/PpiapVq+rhhx/WXXfdFfH8zp07ddlll+njjz/W4cOHbS1kSCL3nqCZDgCA8jlyemLhwoVlAoMkNWzYUAsWLNA999wT6yyTonRgKNlqQAsCAMApQapjuMslAKQAWltRGe5y6bIgpUsA/ha6igxIhKXQkJ6erk8//fTIG9LSlJ6eXu6/jIwMRwvsNDs3qqCnenZAgL/YfRUZUo+lGv6BBx5QkyZNwn87PR4D/CHooQhIFWzLyeP300T0aUgSpwaZclvQvg8ApCr6NHhI6cDg1ebAWK8qITAAQOqIKzRs3rxZN910k1q1aqU6deqoVatWuummm7Rp0ya7yxcoJStYN8NDRZ/pRhm9GqAAAJFiPj2xatUq5eXl6ddff1W3bt3UsGFD7dy5U0uXLlVWVpYWLVqkDh06OFJYP5+eiMZLTftOlMVL3w8AUDErpydiDg29evVSfn6+5s2bp2bNmoWf37Jli/r27avGjRtr4cKFsZfWgqCFhmSIVnF7sTL3YpmAVMH2B8mhPg2ffvqpxowZExEYJKl58+YaPXq0li9fHussEacg3WzLi2UCUgGBAbGIOTQce+yxOvbYY6O+dtxxx6lWrVoJFwrWWNnQ/bwzoK8D4Dw/7yOQfDGHhmHDhumVV16J+tr//u//6vLLL0+4UIDEzgwAvCbm4RtPP/10TZkyRWeccYYuv/zycEfISZMmaffu3Ro6dKimTp0anv6iiy6ytcBACM2qQOJCLXpsS97i1f1bzB0h09IqbpwwDEOhWRqGoaKiovhLV0qqd4S08iPy6g/NKan2fQEEn1v7NUeunli8eHFMhfiv//qvmKavSJBDA5UfAMBNjoQGNwU5NMAehC8gcWxH/hHruqpoeoaR9iGuGADgNgKDd1RWJ8S6rhJdt4QGj2FjTYyV5UcwA+zD9uQsr9UJhAaknNDNwwAkzmuVGpxFaEBKija0NgB7sD0FF6EBEEdLgJ3YnoKL0ABEwZESAJQV04iQs2fP1rRp07R27VoVFBTIMAxlZ2erXbt2uvjii9W3b1+nygkkFUdKgDu43NPbLLU07N+/X+ecc44GDBigN998U4cOHVLz5s3VrFkzHTp0SG+++Wb49V9//dXpMgMAfCxaSx7DWfuDpZaGv/71r1q+fLkmTpyooUOHKiMj8m1FRUWaPHmybrnlFv31r3/VE0884UhhAQD+Fy0YEBb8wVJLw9tvv61HH31Ul19+eZnAIEnp6em67LLL9Mgjj+itt96yvZCA2+jjAAAWQ8OePXvUqlWrSqdr1aqV9uzZk3ChAK/hKAgALIaGk08+WZMmTap0ukmTJunkk09OuFAAAMB7LPVpuPvuuzVs2DD98MMPuvbaa3XqqacqOztbhmGooKBAa9as0euvv64PP/xQb775ptNlBgAALrAUGi677DIVFRXprrvu0qWXXirDMCJeN01TjRo10htvvKFLL73UkYLCXyq7bMqPl1X5scyAX7G9eVNMt8YuLi7WJ598Eh6nQZLq1KmjU089VV26dFFamrNjRXFr7GApvVOIZSeR7B0KOzAgOdjW3GPl1tgxhQa3ERpQ3g6FHQ0QPGzXyWUlNMTcNFBcXKwNGzZo+fLl+vTTT7VhwwYVFxfHVUAESzIuS4x2o6loLRZ245JLIPkIDN5juaXhu+++0wMPPKD3339fBw4ciHgtKytLF1xwgUaNGmXp0sx40dIAr+FICIAT3Ni3WGlpsNQRcuXKlerVq5eOOeYYXXnllWrfvr2ys7MlHRnD4csvv9S0adM0c+ZMLV68WKeddlpiJQds5NTGR2AA4AQv71sstTT069dPBw8e1MyZM1WzZs2o0/zyyy8aNGiQMjMzNXv2bNsLKtHS4CYv/4iTjWUBIIhs69OwbNky3XfffeUGBkmqWbOm7rnnHi1dutR6CYEksNofwep0BAbgd/T3SS2WQkNGRoZ+++23Sqc7dOhQ1HtTwP/8XFFG6zxpN3acCLLQ7zva77x/Tgd+/ynEUmjIy8vTyJEjtW3btnKn2b59u0aNGqXevXvbVjggmRIJRn4OVYBV5f3OKxvIDcFhqU/Dpk2b1KNHDxUUFKh3797hjpAlh5FesGCB6tSpo48++kgtWrRwpLD0aYAf0OcBQcFvOTm8spxtHdzpp59+0qOPPqp3331XGzduVOhthmHoxBNP1ODBgzVixAjVrVs3sVJXgNAAN3llwwYAJzg2IuTBgwe1d+9eSVLt2rWVmZkZe+niQGhAshAQAKQa28ZpKC0zM1ONGjWK562ALxAYAKCsmELD6tWrlZmZqdatW0s6MqT0xIkTtXbtWjVp0kRXXHFFeNAnAEDqoZUu2CxdPbFnzx517NhRp59+uk455RQNHjxYhw8f1sCBA3XNNdfoscce02233aZTTz1V27dvd7rMcBA9nQFI8e8LuAQz2CyFhrFjx+q7777T2LFj9dxzz2nlypW66qqrwldN/PLLL5o9e7aKior00EMPOV1mOIgjBABSYpU/+5HgsnR64r333tODDz6oP//5z5Kktm3bqlevXnrxxRfVq1cvSVLfvn11//3365lnnnGqrAAAwEWWWhry8/OVm5sbftyxY0dJUrt27SKm4/SEv9GkCCCEvgnJ46d9r6XQUK9ePW3dujX8eMuWLZJUZoTIrVu3ql69ejYWD8nEDgJACPuD5PHTsrYUGs466yyNHj1aixYt0hdffKFbb71VXbt21dixY7V7925J0o4dO/TII4+oc+fOjhYY9vNTygUAuMfS4E4bN25Uly5dwgM6NWvWTMuWLdP555+vL7/8Uo0bN9a2bduUlpamTz75RKeddpojhWVwJ/gNTbwArPDCvsK2wZ1atmypr7/+Wu+9956qVKmiCy+8ULVq1dKcOXP097//XatXr1avXr106623OhYYYK9Q64LbP9KgC/VAZzkDCIK4hpF2Cy0N8BsCAwC/sNLSYKlPQ0m7du2q8PXPPvss1lkCgUVgQFDQ9wlSHKGhQ4cOWrBgQdTXnn76afXo0SPhQgEAvCWRAEzgCI6YQ8Mpp5yi/v37a/To0eHbY//nP//RRRddpNtvv1033HCD7YUEAHhXZaGAFrfgiDk0zJs3T/fdd58efvhh9enTRzNmzFCHDh20cOFCTZkyRc8++6wT5QQAeBShIHXE3RFywYIFGjRokH777Te1bdtW06dP1/HHH29z8SLRERIA/IuOwd7mSEdISfr555/1/PPP6+DBg6pfv742bdqkJUuWxDMrIFA4dwuUj8DgfzGHhpUrV6pjx46aP3++3n77bW3cuFEXXnihrr32Wl1//fU6ePCgE+XEUVRK3sZOEYgf+zfvizk0dOvWTTVr1tTnn3+uIUOGqFq1avrXv/6ll19+Wf/+9791xhlnOFHOlFJ6wyn5mHvVAwgqQrf3xRwahg8frmXLlqlly5YRz//hD3/QJ598osLCQtsKl8pCwaDkOcDSozgSHgAAyWT7iJC//vqrqlWrZucsw+gICQD+RCdI73OsI+ShQ4e0a9cu7d69W4cOHYp4zanAECQlWwhoLQCAis3OX5Vy+0qvfl/LoaGgoED33nuvTj75ZFWrVk05OTlq1KiRqlWrppNPPln333+/CgoKnCyrL1hZ0SXTNsnbn7y6QQNeFc++ruQp2VTaV3q5VcbS6YlNmzapZ8+e+vHHH5WXl6f27dsrOztbkrRnzx6tWbNGCxcuVP369bV48WK1aNHCkcL64fREZSvbyz+GVMT6AJKrom2O7dFdVk5PWAoNF198sb799lvNnDlTzZo1izrNDz/8oEGDBqlVq1aaMmVK7KW1wO3QEM8Pmo0gWFifAJLBjX2NbX0aFixYoIceeqjcwCBJzZo105gxYzR//nzrJXRYRU3I8TQvx7MCuUTSm6yuk9LTERgAJINX6w5LoaGwsFBZWVmVTpeVleXqJZelO8uUXOjJXPhUNN5ndZ2w7gC4xYv7H0unJ/r06aOioiLNmDFD1atXjzrN/v37de6556pq1aqaM2eO7QWV3Ds9QZM0ACDorJyeyLAyo0cffVR5eXk64YQTNGTIEJ166qnKzs6WYRgqKCjQmjVrNHXqVP36669atGhRouVOWOnWhpLPx3uKAcFUesAsAED5LA/u9M0332jkyJGaOXOmDhw4EPFaVlaWzjvvPI0ZM0atW7d2pKBS9JaGUBCgNQAA4FdeqMNsu3qipKKiIm3cuDE8JkOdOnXUsmVLpaenx1fKGLh99QQAAE5xu+XTttMTJaWnp6tVq1ZxFQgAADtFq2i9cNQeDz+UOabQYJqmli9frrVr16qgoECGYSg7O1vt2rVTly5dZBiGU+WUxKkIAEDZkSK9eGlivLxev1k+PfHvf/9bd955p/Lz81X6LYZhKCcnR4899pguu+wyRwoqcXoCzvL6xgoATrJtcKe33npLw4YN0ymnnKKJEydq7dq1ys/PV35+vtauXauJEyeqXbt2uuKKKzR5cuUfCnhR6eZNAPFL1jYUtG3V69/HUktDbm6uzjjjDL388ssVTnfjjTdqxYoVWrlypW0FLImWBjiJlgbAXmxT/mJbS8PXX3+tYcOGVTrdsGHD9PXXX1uZJeA5QTs3CrgtntY7tkFvsxQasrOz9d1331U63YYNG8J3vwT8iKMiwD7RBtorHQoYdv8Iv4QlS6Fh6NChuvvuu/X222+ruLi4zOvFxcWaPHmy7rnnHl1yySW2FxJIBr9stICflQ4PtPAdEVoOXl8Wlvo07N+/X4MHD9a8efNUs2ZNtWnTJmIY6fXr12vfvn3q06ePpk2bpmrVqjlSWPo0AECw0O/BO2zr01C9enXNmTNHM2bM0CWXXKK0tDRt3LhRGzZsUFpami677DLNnDlTs2fPdiwwAHbweooHvMyJ7SdVWxr8+p1jHkbaTbQ0IFEc1QCV8+p24tVyxark9/DSd7KtpQEICq9snICXld5Ooh0Vu3Gk7Iftt7zlUt7dl/3wnUqyHBq2bt2qp59+Wi+//LL+85//hJ+7/vrr1aVLFw0ZMkQfffSRYwUFADgvWqUXrWJjMLTorCwXPy8vS6cn1q1bp27duunnn3+WJB1//PFatGiR8vLytGfPHp144on69ttvdeDAAX388cfq3LmzI4Xl9ARi4aVmPwCJYXt2nm2nJ0aPHq369etr6dKl+uqrr3TiiSfqggsuUN26dbV582atWLFC33//vdq1a6eHHnoo4YIDdkjVDlaAm2Ld5qxOT2DwBkuhYenSpRo1apTOPPNMtWnTRs8884xWr16tO++8U8cee6wkqU6dOhoxYoRjQ0gDVjBQDOCuWMJ66Q6BpV/zMq+XzymWQsPevXvVpEmT8OOmTZtKkho2bBgxXU5Ojvbs2WNj8QBrSg4Uk8j7gVRk9++/vCBQ3nSl/4722GtStSXTUmho0qSJPv/88/DjFStWSFKZVoUvvvgiIlwAyRLPDqa83sxAKknVvgJ2VPipuNwshYYhQ4Zo1KhRGj16tJ544gldddVVuuGGG/Tggw9q2rRp2rp1qyZPnqxx48bpnHPOcbrMSHElN/ZENnxaJQDrrQKJzDsaL25HfhjG2W2Wrp745ZdfNGjQoPAllRdeeKHefvtt3XzzzXrttddkGIZM01Tz5s21fPly1a9f35HCcvUESnP6KCnR0x6A1yXzN+6XVg2/lNNuVq6eiGlEyI0bN6pKlSpq1qxZ+LkPP/xQq1evVk5Oji666CJVr149vtJaQGiAHazsEFJ1pwGkArbv6GwPDW4jNMALou1wrD4HeBG/VUgMIw1Icq5neDzPAYCfERoQeHZV3pWFj/I6aIb+ppMVvCpZ24iX+KmsXkJogK9Y2dBLVtJ2qmzHWvqa89IdzBirH15lV6D1U+uan8rqJYQGeF6sO7NolbQbZanolEWqDgwDb+qf0yElK1G2wdgRGuB5XrqNrJ2fb3Ve7NjghCDegTEWdP6MD1dPAD7FGBIA7MTVE0BAhY6SCAzB4+SRfrR5V9YHKFVaHmANoQHwAXbcqcPJIFhZP5tUwfYUP0ID4DHRdmil7w8QbSApBINd91aJZ57RrqII2lU/oVa6IHwXN9CnAQA8wIsd80pXsF4rH+xFnwYghXEk5S2VrY9oFbJb6zDaGCMcnUMiNAC+Y3XHHbRmZb+L55RSokf28a738j7Xzy0NbAP2IDQAKcDPO/ugSka/FNZ7dASI+NGnAfCReM97c04aQGWs9GnISEI5ANgk3kqfsOBvJY+MWZdwE6cngBRAc6y/hToixhoYWO+wG6EB8Jl4KgKOTp1X3u3QnfwcINlsDQ1paWlq0qSJXnrpJRUWFto5awDy5rX8OKLkAFx23s209IBLpQf6slImwC62doTs1auX9u/fr9WrV6tJkyb6/vvv7Zq1JDpCAokgcNjPi8vUi2WCPyS9I+SiRYskSfv27dOSJUvsnDWABEU7EkZivLQcWa9IBkf6NNSoUUMDBw50YtZAyot2f4BYULFEimVZVnSXSLexXpEMMYeGr776qsLXp0+fHndhAFSOW2LbK9TvIN4+Al5bF14rD5LL6WAbc2jo0qWLXn/99TLPFxYW6vbbb9fgwYPtKBcAJFVlla1XWhSQWkr/7qzew8SpsT1iDg1DhgzRddddp+HDh+vAgQOSpM2bN6t79+564YUX9Oijj9pWOADOoiI8wspOlSN4uKH0787q79CpG6DFHBpef/11vfrqq5oyZYo6d+6s5557Trm5udq1a5eWLFmiv/zlLwkXCkByWN0BWd3ZRBurwOvKK6dfyo+KJdoHKEjsCL5xdYS89tpr9X//93/6/vvvddttt+mkk07S6tWr1aVLl4QLBMAaL+4IS+6U/HIrZSt3A/XD90B0fu0D5NXfXFyhYevWrfrjH/+owsJCtW/fXitXrtTTTz9td9kAVCCZO8J4A4DfdtZBvCU0/Mmrv7mYQ8PMmTOVm5ur7du3a/Hixfr88891991368EHH1Tfvn21e/duJ8oJwEWlxwCIJ0B49cipPIl+XyCIYh4RMj09XQMGDNAbb7yh7Ozs8PNz587VlVdeqYyMDG3fvt32gkqMCAkESSiIeG1QIq+VB6nD7d+elREhY25pGDdunGbMmBERGCSpb9++WrVqlVq1ahXrLAHEIZlHv058VrTA4IUj+lju7QDYyQ9h1dZ7T0hScXGx0tKcuXkmLQ2Ae9w+CnJbqn9/JJcbvzdHWhoqnaFDgQGAc7dcrmh+oddStcKkxQFu8Or2ZqmloXfv3nrhhRd08sknq3fv3hXP0DA0f/582wpYEi0NgP1KHtEE9Wg6qN8LsJNtLQ0lc0VxcbFM0yz3X3FxcfwlBpB0pcdW8AK7ju7tbCWhxQFO8svvy9KtsRcuXBj+O3T7awDJ5+SpgsquZkjm0bpXwkuIX3boiM4PLU1eL1+I7R0hncTpCaQyP+z4gozlDzt58fdk5fREXKGhqKhIb7/9thYuXKiCggLVqVNHeXl5Gjp0qDIyLDVexIXQADjHizsxOwT1e6FyrPvYOBIafvrpJ51zzjn64osvlJGRoTp16qigoECFhYXKzc3V7NmzVbdu3bgLXRFCA1KdmzvBoOyAg/I9ALs5csnl7bffrm+++UYTJ07UgQMHtGPHDh04cEATJkzQd999p9tvvz2uwgKonBs3gSrdj8KP5/dDdzos2W8DqYF1ba+YWxqOO+44jRkzRrfddluZ15566imNHj1a/+///T+7yheBlgYAAJzhSEuDaZpq27Zt1NfatWsnH/WrBGBR0I/Won2/oH9nuMfPv62YQ0OfPn00b968qK/NnTtXvXr1SrRMADzG7T4ATu9kS38/+j3ASX4+RRbzpQ4jR47URRddpKKiIg0bNkwNGzbUzp07NXHiRE2dOlVTp07Vnj17wtOXvrEVAMSKChxBULJ/kF9/0zH3aSh5bwnDMMJ/h2ZT8jnpyOWZdqFPAwC7ePXW3LAf69gaK30aYm5peOCBB8oEAwDJQUWXuNCyCy0/lmOwsa3YixEhAQCAfVdP5ObmauzYsVq/fn3ChQKARMXbicyvnc+QGNa7fSyFhry8PL3yyitq166dTjnlFI0cOVKrVq1yuGgAEF28zc00U6cmL613vwcYS6HhiSee0KZNm/TJJ5/o/PPP11tvvaXTTz9dLVu21F133aVPP/3U6XICQFisO97QaJAAEhN3n4bVq1frnXfe0TvvvKP169ercePGuuiiizRkyBD16NHDkc6S9GlAqqNTV2xYXoB1jt3lsrT169eHA8Tq1avVoEED7dixI9HZlkFoAGBVycBAeEhdrHvrkhYaStq4caOmTp2qO++8087ZSiI0AKVvHpWqqAjgJ375vboSGpxEaAAAxMIvFbYX2HbJZe/evSP+AXAPHfpiE0+nSQSHlwJDEH5blkaEbN68udPlAICYWD2CjLXS8FIlg2AJwm+L0xMAAinRZmmatWEHP/2ObDs9AQB+k+iO2i87elTOzdMCQfsdERoA+JbTlUEQzkFD4Zu8IXGEBgC+lIxKIGhHianMjXUZxKAS862xAcALqNDhdUH8jdLSAHhQEI9QAPgfoQHwoCAeoTihdLgqeWMqghdgP0IDAN8qeW+J0s8RvOAlQQmxlkLD7t27VVRUFPHcmjVrNHjwYDVq1EiNGzfWkCFDtG7dOkcKCQAVISCgMtwe3R6WBndKT0/XsmXLdMYZZ0iSvvrqK5155pkyDEM9evSQaZr6+OOPlZ6erhUrVujEE090pLAM7oQg89MgMF7GcgTiY9vgTqVzxQMPPKDatWtrzZo1mjVrlj744AOtXr1a1atX18MPPxxfaYEUR0VnD5YjvCgorRxx9WlYtGiR7r333oh7UrRo0UJ33nmn5s+fb1vhAAAIgqCE2bhCw3/+8x+deuqpZZ5v3769du/enXChgFQUlCORIGGdAJEsD+70zTffKCPjyOT16tXTzz//XGaan3/+WVlZWfaVDkghQTkSAXBEEPvXWA4Nw4cPD/9tmqY++eQTDRw4MGKaL7/8Uk2bNrWtcADgpqDt8JFcQfz9WAoN48ePL/Ncw4YNyzy3ZMkS9evXL/FSASkiiEciQcM6CpbQ+nR6vQb1d2Ppkkuv4JJLBE1QdyxBwjpCvPz227HtkksAzvDTDiVVsY5ST2UdYK12kC392wlCx1pbQ8PWrVv1ww8/2DlLAHAVIwn6W2j9WVmHoWlCpy/KE0uQjDbEeXnT+eF3ZuvpiSpVqsg0TRUWFto1ywicngAABE200xhunNpI+umJq6++WldddZWdswQAV/nh6A/+VjIclGzt8CLLl1xa8eqrr9o5OwAAUopXw0IIHSEBoAKVnd8GUomtoeHgwYN0hAQQSAQHwObQMHPmTLVo0cLOWQKA67zeZIzyEfbsxekJALAgWmc1eB+Bz16WOkI++OCDlma2bt26hAoDAH4Qqoj8NuJfKgj6OnH7+1kapyEtLU2GYcjKkA6GYaioqMiWwpXGOA0ICq9fVoXKub3zBuxmZZwGSy0NdevW1eDBgzVu3LgKp5s+fbr+8Ic/WCsdkMKobIKFEIhUYSk05Obm6ttvv1WdOnUqnK5WrVq2FAoIOo5S/a/k+mNdegfblrMsdYQ87bTTtHr16kqnq169upo1a5ZwoYBUQGc6wH5BDAwV7SuSvR+x1Kdh3759KigoUPPmzZNRpnLRpwFBwdEQACuSua+w7d4TNWrUcD0wAEFBYEgNtCQhUV7cVzBOA5BkXtsJwBms5+QKWkjzYmCQCA0AgADwYgWbiHi+TzKCE6EBcEnQjowAuCsZwYnQALgkaEdGKItg6KzZ+asCuYy9/J0IDYCLvLxzQOJKDjcN+wU1eHv5exEaAMBhXq4E/I5lW76SYdWu4EpoAFzEDi/10OpgL5Zn+ZwYtZTQALiEnV3q8epldH4W5OXpxX2EpREhAQAAaGkAAACWEBoAAIAlhAYAAGAJoQEAAFhCaAAAAJYQGgAAgCWEBgAAYAmhAUgRBQUFql+/vjZv3ux2UWwzYsQI/fd//7fbxQBSBqEBSBF/+9vfdN555+n444+P+vrmzZtlGEaZfx9++GHEdIsXL1bHjh2VmZmpE044QS+99FLMZRk9enTUz6pevXp4mkWLFkWd5uuvvw5Pc9ddd2n8+PHatGlTzGUAELsMtwsAwHkHDhzQq6++qlmzZlU67bx589S2bdvw4+zs7PDfmzZt0sCBA3XDDTdowoQJ+vjjj3XLLbeoXr16uvjiiy2XZ8SIEbr55psjnjv77LPVuXPnMtN+8803qlWrVvhxvXr1wn/Xr19f/fr100svvaRHHnnE8ucDiA+hAUgBH3zwgTIyMtS1a9dKp61Tp44aNmwY9bWXXnpJzZo101NPPSVJatOmjT777DM9/vjjMYWGGjVqqEaNGuHHq1ev1rp166K2WtSvX1/HHXdcufM6//zzNXLkSEIDkAScngBSwJIlS9SpUydL055//vmqX7++unfvrilTpkS8tmzZMvXr1y/iuf79++uzzz7T4cOH4y7fK6+8olatWqlnz55lXsvNzVWjRo109tlna+HChWVeP+OMM7R161Zt2bIl7s8HYA2hAUgBmzdvVk5OToXT1KhRQ0888YSmTJmiWbNm6eyzz9all16qCRMmhKfZuXOnGjRoEPG+Bg0aqLCwUD/99FNcZfvtt980ceJEXX/99RHPN2rUSP/4xz/0zjvvaOrUqWrdurXOPvtsLVmyJGK6xo0bh78jAGdxegJIAQcOHFBmZmb4cdu2bcNH5j179tQHH3ygunXr6vbbbw9P06lTJ+3du1ePPvqorrzyyvDzhmFEzDt0o9zSz1f0WSVNnTpVv/zyi66++uqI51u3bq3WrVuHH3ft2lVbt27V448/rrPOOiv8fFZWliTp119/rWQpAEgUoQFIAXXr1tXevXvDj2fNmhU+nRCqdKM588wz9corr4QfN2zYUDt37oyYZvfu3crIyFCdOnWizqOyz3rllVc0aNCgcvtRlC5PyZYPSdqzZ4+kyA6SAJxBaABSQG5ubkRl27x5c0vvW7lypRo1ahR+3LVrV73//vsR08yZM0edOnVSlSpVos6jos/atGmTFi5cqOnTp8dVHklau3atqlSpEnHFBwBnEBqAFNC/f3/de++92rt3r2rXrh11mn/+85+qUqWKcnNzlZaWpvfff1/PPPNMxFUJN998s5577jndcccduuGGG7Rs2TK9+uqrmjRpUlzleu2119SoUSMNGDCgzGtPPfWUjj/+eLVt21aHDh3ShAkT9M477+idd96JmO6jjz5Sz549K2wxAWATE0BKOPPMM82XXnqp3Ndff/11s02bNma1atXMmjVrmh07djT/9a9/lZlu0aJFZm5urlm1alXz+OOPN1988cWI1xcuXGhKMjdt2lRheYqKiswmTZqY9913X9TXH3nkEbNly5ZmZmamWbt2bbNHjx7mzJkzy0zXqlUrc9KkSRV+FgB7GKZ5tBcTgECbNWuWRowYobVr1yotzbkLp15//XWNHTtW69atK/eUhV1mzpypO++8U19++aUyMmg4BZzGVgakiIEDB+q7777T9u3b1bRpU8c+58MPP9S4ceMcDwyStH//fo0fP57AACQJLQ0AAMASBncCAACWEBoAAIAlhAYAAGAJoQEAAFhCaAAAAJYQGgAAgCWEBgAAYAmhAQAAWEJoAAAAlhAaAACAJYQGAABgCaEBAABYQmgAAACWEBoAAIAlhAYAAGAJoQEAAFhCaAAAAJYQGgAAgCWEBgAAYAmhAQAAWEJoAAAAlhAaAACAJYQGAABgSYbbBUBqOXjwoA4dOuR2MQDYqGrVqsrMzHS7GEgCQgOS5uDBgzo2q7YO6aDbRQFgo4YNG2rTpk0EhxRAaEDSHDp0SId0UD2MQaqSVlUy0mSkGZJhSEf/N9LSjj5OO/o4PfxayedlKOJ9MtJ+/1uR8zQN48iJOOP3+UQ+J5klpz06b/Po66ZRYrqj7zPDj49OoyMfq7TQ9CVeK/H4yOfo6Lx/n0Yq+Tja/0bE43KnDc0nTb+XKcq0lc1LKufv0u9JM8udpwzzSBkivptZaj5mmdeM0tPJlJF2ZE5GqWkMw4yY3jDMo6vP/P0nE57OjHicppLTmkozTKWnFR95/ujjkv+OrN5Szyv0d/HRaULPFSvd+P3vNENKlynj6HRVjKIjn6cjj9ONYhk68n+aUax06ej/R98fmk7FR6YJTatipUlH51N89HVTaUeny1CR0kPzOfqe9KPfL12h+ZtKl3n0bx0ps6R0Q0qXcfRvQ2n6/d+Rx2lKk6H9+6QWHbfo0KFDhIYUQGhA0mWoijKMKkdCQ6jCP7qHN8J/VxYaIoNBRGgo9VrlocGQmZ5gaAhVcHGGhvLDgn7/XCuhIfR3oqGhvABRetrKQsPR99sSGozyQ8Pvz0cGgUpDQ8nHlYSG3wNC+aGh5OPI0GD+HgAMU1WMtDKhIS0iNJil/i8/NKQb5pHK3TCULuPo5x6p2KvoSIV/JDQUHw0CoaCgqKEhvYLQkB4lNKSH1jFSAh0hAQCAJYQGAABgCaEBAABYQmgAAACWEBoAAIAlhAYAAGAJl1wi6Qp1WMaR6/WO/v/7tXpG+FrEI/8bxekKX2NnHHlPxDgN4dfSVOJavCP/myUuozRV9pLL8HOSWXLaeC65lBK65FKKWAyM02DjOA1mqUsuQ4/No+M0FJe45NJM0jgNxYEYp0HhcRqQOggNSBrTNFWjRg39374ZUpHbpQFglxo1asg0zconhO8RGpA0hmFo37592rp1q2rVquV2cQDY4Oeff1bTpk2PDNSGwCM0IOlq1apFaAAAH6IjJAAAsITQAAAALCE0IGmOOeYYjRo1Ssccc4zbRQFgE7br1GKYdHkFAAAW0NIAAAAsITQAAABLCA0AAMASQgMAALCE0ICkeOGFF9SiRQtlZmaqY8eO+uijj9wuEoAELVmyROedd55ycnJkGIbeffddt4sEhxEa4Li33npLf/7zn3X//fdr5cqV6tmzpwYMGKAffvjB7aIBSMD+/ft12mmn6bnnnnO7KEgSLrmE47p06aLTTz9dL774Yvi5Nm3a6MILL9Tf/vY3F0sGwC6GYWjatGm68MIL3S4KHERLAxx16NAhff755+rXr1/E8/369dPSpUtdKhUAIB6EBjjqp59+UlFRkRo0aBDxfIMGDbRz506XSgUAiAehAUlR+ra5pmlyK10A8BlCAxxVt25dpaenl2lV2L17d5nWBwCAtxEa4KiqVauqY8eOmjt3bsTzc+fOVbdu3VwqFQAgHhluFwDBd8cdd+iqq65Sp06d1LVrV/3jH//QDz/8oJtvvtntogFIwL59+7Rhw4bw402bNmnVqlXKzs5Ws2bNXCwZnMIll0iKF154QY8++qh27Nihdu3a6cknn9RZZ53ldrEAJGDRokXKy8sr8/w111yj119/PfkFguMIDQAAwBL6NAAAAEsIDQAAwBJCAwAAsITQAAAALCE0AAAASwgNAADAEkIDAACwhNAAAAAsITQAAABLCA0AAMASQgMAALCE0AAAACz5/7q3eOU7+PeRAAAAAElFTkSuQmCC\n",
"text/plain": [
"