diff --git a/.flake8 b/.flake8 index 5cc6e1d..9634ee8 100644 --- a/.flake8 +++ b/.flake8 @@ -1,7 +1,7 @@ [flake8] ; max-line-length = 127 -max-complexity = 10 +max-complexity = 12 statistics = true count = true show-source = true -ignore = E501, W503 \ No newline at end of file +ignore = E501, W503, E203 \ No newline at end of file diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml index 13f6eae..371b0c9 100644 --- a/.github/workflows/tests.yaml +++ b/.github/workflows/tests.yaml @@ -49,5 +49,8 @@ jobs: run: | python .github/gee_token.py - name: Test with unittest + env: + S3_ACCESS_KEY: ${{ secrets.S3_ACCESS_KEY }} + S3_SECRET_KEY: ${{ secrets.S3_SECRET_KEY }} run: | python -m unittest discover -s ./tests -p '*.py' \ No newline at end of file diff --git a/README.md b/README.md index b6f742b..92757d1 100644 --- a/README.md +++ b/README.md @@ -27,10 +27,10 @@ You want to download Earth Observation data but don't want to spend hours just f We currently support these data providers: - [Planetary Computer (pc)](https://planetarycomputer.microsoft.com/catalog) - [Google Earth Engine (gee)](https://developers.google.com/earth-engine/datasets) +- [Copernicus Data Space Ecosystem (cdse)](https://dataspace.copernicus.eu/explore-data/data-collections) (not all collections supported) Coming soon: - [Alaska Satellite Facility (asf)](https://asf.alaska.edu/how-to/data-basics/datasets-available-from-asf-sar-daac/) -- [Copernicus Data Space Ecosystem (cdse)](https://dataspace.copernicus.eu/explore-data/data-collections) ## Usage ### Installation @@ -75,15 +75,11 @@ Other data backends work with the same principle, check out the [Demos](https:// ## Limitations Be aware that depending on the collection the data is not mosaicked. -## Limitations -Be aware that depending on the collection the data is not mosaicked. - -## Limitations -Be aware that depending on the collection the data is not mosaicked. - ## Contribute You found a bug or a data source is missing? We encourage you to raise an issue or provide a PR. +We are looking for contributors to add more collections to CDSE. Please get in touch if you are interested. + ## License This work is licensed under the MIT license. diff --git a/docs/demo_files/demo_cdse.ipynb b/docs/demo_files/demo_cdse.ipynb new file mode 100644 index 0000000..0b8ee0d --- /dev/null +++ b/docs/demo_files/demo_cdse.ipynb @@ -0,0 +1,1487 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Downloading with CDSE" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "# install with\n", + "# !pip install terragon-downloader[cdse]" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import terragon\n", + "import geopandas as gpd\n", + "import xarray as xr\n", + "from utils import visualize_sat_images\n", + "import rasterio" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "\n", + "Name: WGS 84 / UTM zone 32N\n", + "Axis Info [cartesian]:\n", + "- E[east]: Easting (metre)\n", + "- N[north]: Northing (metre)\n", + "Area of Use:\n", + "- name: Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.\n", + "- bounds: (6.0, 0.0, 12.0, 84.0)\n", + "Coordinate Operation:\n", + "- name: UTM zone 32N\n", + "- method: Transverse Mercator\n", + "Datum: World Geodetic System 1984 ensemble\n", + "- Ellipsoid: WGS 84\n", + "- Prime Meridian: Greenwich" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# we will use the sample data to download a minicube\n", + "# other datasets can be found here: https://dataspace.copernicus.eu/explore-data/data-collections\n", + "gdf = gpd.read_file(\"data/TUM_OTN.geojson\")\n", + "gdf = gdf.to_crs(\"EPSG:32632\")\n", + "arguments = dict(\n", + " shp=gdf,\n", + " collection=\"SENTINEL-2\",\n", + " start_date=\"2021-01-01\",\n", + " end_date=\"2021-01-05\",\n", + " bands=[\"B02\", \"B03\", \"B04\"],\n", + " resolution=10,\n", + " download_folder=\"tests/download/\",\n", + " filter={\n", + " \"processingLevel\": {\"eq\": \"S2MSI2A\"} # we will use a filter to get only level 2 products\n", + " },\n", + ")\n", + "gdf.crs" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "credentials" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# we use the s3 bucket to download the data: https://documentation.dataspace.copernicus.eu/APIs/S3.html\n", + "# credentials can be created here: https://eodata-s3keysmanager.dataspace.copernicus.eu/panel/s3-credentials\n", + "credentials = {\"aws_access_key_id\": \"\", \"aws_secret_access_key\": \"\"}\n", + "tg = terragon.init(\"cdse\", credentials=credentials)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "download (with eodata AWS S3 backend)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['COP-DEM', 'S2GLC', 'TERRAAQUA', 'SENTINEL-3', 'SENTINEL-5P', 'SENTINEL-1-RTC', 'SENTINEL-1', 'SMOS', 'LANDSAT-7', 'CCM', 'LANDSAT-5', 'ENVISAT', 'LANDSAT-8-ESA', 'SENTINEL-6', 'GLOBAL-MOSAICS', 'SENTINEL-2']\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/mnt/datastorage/home/adrianhohl/Terragon/terragon/copernicus_data_space_ecosystem.py:49: UserWarning: Currently we only support the following collections: ['COP-DEM', 'GLOBAL-MOSAICS', 'LANDSAT-5', 'LANDSAT-7', 'LANDSAT-8-ESA', 'TERRAAQUA', 'S2GLC', 'SENTINEL-1', 'SENTINEL-1-RTC', 'SENTINEL-2']\n", + " warnings.warn(f\"Currently we only support the following collections: {supported_collections}\")\n" + ] + } + ], + "source": [ + "# collections can be viewed using the following command\n", + "# the unsupported collections 'SENTINEL-3','SENTINEL-5P','SENTINEL-6','ENVISAT','SMOS', rely on different file extensions\n", + "# and therefore are not supported with the current workflow, if you are interested in these collections and have already worked with them,\n", + "# please open an issue and describe how you would process them or even better provide a PR\n", + "collections = tg.retrieve_collections()\n", + "print(collections)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 6kB\n",
+       "Dimensions:      (x: 27, y: 18, time: 2)\n",
+       "Coordinates:\n",
+       "  * x            (x) float64 216B 6.976e+05 6.976e+05 ... 6.979e+05 6.979e+05\n",
+       "  * y            (y) float64 144B 5.326e+06 5.326e+06 ... 5.326e+06 5.326e+06\n",
+       "  * time         (time) datetime64[ns] 16B 2021-01-01T10:23:29.024000 2021-01...\n",
+       "    spatial_ref  int64 8B 0\n",
+       "Data variables:\n",
+       "    B02          (time, y, x) uint16 2kB 65535 65535 65535 ... 65535 65535 65535\n",
+       "    B03          (time, y, x) uint16 2kB 65535 65535 65535 ... 65535 65535 65535\n",
+       "    B04          (time, y, x) uint16 2kB 65535 65535 65535 ... 65535 65535 65535\n",
+       "Attributes:\n",
+       "    crs:          EPSG:32632\n",
+       "    data_source:  CDSE\n",
+       "    collection:   SENTINEL-2
" + ], + "text/plain": [ + " Size: 6kB\n", + "Dimensions: (x: 27, y: 18, time: 2)\n", + "Coordinates:\n", + " * x (x) float64 216B 6.976e+05 6.976e+05 ... 6.979e+05 6.979e+05\n", + " * y (y) float64 144B 5.326e+06 5.326e+06 ... 5.326e+06 5.326e+06\n", + " * time (time) datetime64[ns] 16B 2021-01-01T10:23:29.024000 2021-01...\n", + " spatial_ref int64 8B 0\n", + "Data variables:\n", + " B02 (time, y, x) uint16 2kB 65535 65535 65535 ... 65535 65535 65535\n", + " B03 (time, y, x) uint16 2kB 65535 65535 65535 ... 65535 65535 65535\n", + " B04 (time, y, x) uint16 2kB 65535 65535 65535 ... 65535 65535 65535\n", + "Attributes:\n", + " crs: EPSG:32632\n", + " data_source: CDSE\n", + " collection: SENTINEL-2" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# we can now download the data, in the background we search the folder of the corresponding tile for the correct files\n", + "# and use rasterio with its virtual files to open and clip the extend of the files directly,\n", + "# it is possible to download data with up to 4 concurrent workers per account, see here: https://documentation.dataspace.copernicus.eu/Quotas.html\n", + "# make sure to have the right packages installed (depending on the data you want to download), for example jp2 needs libgdal-jp2openjpeg and .nc needs libgdal-hdf5\n", + "ds = tg.create(**arguments, num_workers=4)\n", + "# or\n", + "# items = tg.search(**arguments, num_workers=4)\n", + "# ds = tg.download(items)\n", + "ds" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAADECAYAAAACqSJyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCpUlEQVR4nO3dd3wU17n4/89WrVZdSEJIIEQXvYoOAgQIMFgYbH+dcuMUp13H6YkTx3FLnOT+ch3HThw7107iJHbsJCYGU0XvYJrpRTQ1mnrXatv8/hhAlmGPJJC0kvZ5v168jPfZmTkrNLPPnDnnOQZN0zSEEEIIEbCM/m6AEEIIIfxLkgEhhBAiwEkyIIQQQgQ4SQaEEEKIACfJgBBCCBHgJBkQQgghApwkA0IIIUSAk2RACCGECHCSDAghhBABTpKBFti8eTNf/OIXSUlJISQkhMTERDIzMzl48OAt7z106BCzZ88mNDSUyMhIlixZwoULFxq9Jzs7m+9///uMHTuWyMhIoqOjmTJlCu+9994t+ysoKODb3/42aWlpREZGYjAYePPNN1v8GZrTLoDf/va3LFmyhD59+mAwGJgxY0aLj1VdXc23v/1tEhISsNlsjBo1inffffeW9+3cuZNHHnmEsWPHEhQUhMFgICcnp8XHE6KjkmtHyzT32vHyyy8zceJEYmJiCAoKIikpiYceeogTJ060+JiBTpKBFnj11VfJycnhW9/6FmvWrOGll16isLCQiRMnsnnz5pvvO336NDNmzMDpdPKvf/2LP//5z2RnZzNt2jSKiopuvm/9+vWsXr2apUuX8u9//5u3336bAQMG8MADD/Dcc881Ova5c+d4++23sVqtLFiw4I7a39x2Abz22mvk5uYya9YsYmNj7+h4S5Ys4a9//StPP/00a9euJTU1lU996lP84x//aPS+TZs2sXHjRpKSkpg8efIdHUuIjkyuHS3T3GtHSUkJ8+fP54033mD9+vU8++yzfPTRR0yYMIEzZ87c0bEDliaa7dq1a7e8VlVVpXXv3l1LT0+/+doDDzygxcTEaBUVFTdfy8nJ0SwWi/bDH/7w5mtFRUWa1+u9ZZ/33HOPZrfbNYfDcfM1j8dz8+/79+/XAO0vf/lLi9rf3HZ98nhDhw7V0tLSWnSs1atXa4D2j3/8o9Hrc+bM0RISEjS3233bY/3617/WAO3ixYstOp4QHZlcO5qvJdeO2zl58qQGaD/96U9bdNxAJz0DLRAXF3fLa6GhoQwZMoT8/HwA3G43q1atYunSpYSHh998X+/evZk5cybvv//+zddiYmIwGAy37HP8+PHU1tZSWlp68zWj8e7+qVrSrtY43vvvv09oaCgPPPBAo9e/8IUvcPnyZT788MNWO5YQHZ1cO5qvJdeO27nRG2E2m++qHYFGrsJ3qaKigkOHDjF06FAAzp8/T11dHSNGjLjlvSNGjODcuXM4HA7lPrds2UJsbOxtLyB3qjXa1RLHjx9n8ODBt5yQN45//PjxVjuWEJ2RXDtu706uHR6Ph/r6ek6fPs0jjzxCXFwcX/jCF1qtTYFAkoG79Oijj1JTU8NPfvITQH+GBRAdHX3Le6Ojo9E0jbKyMp/7e+ONN9i6dStPPvkkJpOp1dp5t+26k+P5OtbH2yNEoJJrh+/jtfTaERISgs1mY/DgwZw6dYqtW7fSq1evVmtTIJBk4C789Kc/5e233+bFF19k7NixjWK368JrKrZ27VoeffRR7r//fh577LE7apPX68Xtdt/84/F47rpdvmia1uhYbre7zY4lRFci147WvXbs3r2bPXv28NZbbxEWFsbMmTNlRkELSTJwh5599ll+/vOf8/zzz/ONb3zj5uvdunUDbp+9lpaWYjAYiIyMvCWWlZXFkiVLmDNnDm+//fYdf1l+8YtfxGKx3PyTnp5+V+1S2bZtW6NjWSyWm1MCu3Xr5vNYcPu7DCECgVw7Wv/aMWbMGCZOnMhnPvMZtmzZgqZpPPHEEy1qU6CTERZ34Nlnn+WZZ57hmWeeueUXrl+/fgQHB3Ps2LFbtjt27Bj9+/fHZrM1ej0rK4vFixeTlpbGsmXLsFqtd9y2Z555ptEFJiws7I7b1ZSxY8eyf//+Rq8lJCQAMHz4cN555x3cbnejZ383jj9s2LAWHUuIrkCuHbq2vHaEhYWRkpJCdnZ2i9oU8Pw4k6FTeu655zRAe/LJJ32+58EHH9Ti4uK0ysrKm6/l5uZqVqtVe/zxxxu9NysrS7PZbNrs2bO1urq6ZrXhTqcHtaRdH3cn04PWrFmjAdq7777b6PV58+YppwfJ1ELRVcm1o3nu9NpxQ1FRkRYVFaUtXLiwRccNdNIz0AIvvPACTz31FPPmzeOee+5h7969jeITJ04E9Ow/NTWVhQsX8qMf/QiHw8FTTz1FTEwM3/ve926+f+fOnSxevJj4+HieeOIJDh8+3Gh/Q4YMaTSV50Z1sRtVvw4cOEBoaCgA999/f5Ptb267buz7RrddZWUlmqbdPH5qaiq9e/dWHmv+/PnMmTOHr3/961RWVtK/f3/eeecd1q1bx1tvvdVogFNRURHbtm0DGrL/tWvXEhsbS2xsLGlpaU1+NiE6Mrl2tP61o6Kigjlz5vDpT3+aAQMGEBwcTHZ2Ni+99BL19fU8/fTTTX4u8TH+zkY6k7S0NA3w+efjDhw4oKWnp2t2u10LDw/XFi9erJ07d67Re55++mnl/rZs2dLo/c09tkpz2qVpmvbwww/7PFZz7yqqqqq0b37zm1p8fLxmtVq1ESNGaO+8884t79uyZYvPY7X0rkKIjkiuHa1/7XA4HNojjzyiDR48WAsNDdXMZrPWs2dP7bOf/ax24sSJZn8uoTNomqa1TlohhBBCiM5IZhMIIYQQAU6SASGEECLASTIghBBCBDhJBoQQQogAJ8mAEEIIEeAkGRBCCCECnCQDQgghRICTZEAIIYQIcJIMCCGEEAFOkgEhhBAiwEkyIIQQQgQ4SQaEEEKIACfJgBBCCBHgJBkQQgghApwkA0IIIUSAk2RACCGECHCSDAghhBABTpIBIYQQIsBJMiCEEEIEOEkGhBBCiAAnyYAQQggR4CQZEEIIIQKc2d8NEMJf0oeP9xnbdGxfO7ZECNGZxKZu9Rkr2j+j3drRmqRnQIiPsXi9jKmugB/+EF5/3d/NEUJ0EgnOQj5btBq+/GXIy/N3c1pMegZEYNM0ejodpFZXMq6qgpE1VQRrXnj5ZbBa4XOfg6Agf7dSCNHB2Lz1TKo6wszK/cyq2M8gRy4ejPCmEeLi4Pnn/d3EFjFomqb5uxFCtKuKCti0iVWPfI1x1RXEu5w4DQaO20PZHxrBgbAIXl/+Txg2DFauhIUL/d1iIYS/aRqcPAlZWWx59h0mVh0lWHNyyRLLlohUtoSnsj18DGdHvQ979sCpU2Aw+LvVzSbJgOj6vF44eBCysvQ/e/aAx0O+1cb+0HAOhEVwJCQMh9F0c5NNx/bB4MEwcSL85S9+bLwQwm/KymDjRli3Dtavh4ICsNnYYh3G5vBUNkeMJ9vWu9GXftHT1bBoERw/DkOH+rHxLSOPCUTXdOVKw5f/hg1QUgJhYTB7NrzyCmRk8PlFD6r3sXQp/OEP4HKBxdI+7RZC+I/HA/v2NVw79u3TbyaGDIEHHoCMDJg+nQenf+h7H3Pm6NeaZcskGRCi3dXXw86dDSfx0aN6tj52LHzta/pJPHFiy77Uly7Vn/tt3aqf4EKIrqegoPGNQ3k5REbqNw6PPAJz50KvXs3fX1CQ/mhx2TJ46qm2anWrk8cEosP6r8xHfQc1jX79Ixhy9jRDzp5hwIXzBLmcVISFcXJACqcGDGJPTTBVQXafu7hUX+EzVp1zCjSNldlH2BsawfOJfRrFPzq2t8WfRwjRPgyKqX82bz0z6s8wt2wvc8v3MrT2Ih6M7AsbwvrIiayPmsjZXlY8H3tseIviRJ+hqloH95Xv5J+5P2dIyp84F9T4vc7DGS39OO1CegZEpxHsqmdI8SWGF+YxojCf2JVVuEwmzvfuy5r0uZwcMIhL8Qk3n99VHcy9uwMaDGyMiGZRWRG/TEjG24kGAwkhrtM0BjtyyajYz7zKfUy/PvCvwBpHVtREnkv6MpsiUimzRNzcJNJ49q4OuS58HDXGIO4r38mvu/+/u/0E7UKSAdFhGTSN5PJCRhTlM7wwn/5lVzFpGldCIvgoPplz0yaQ3bcfTmvbTf3bFB7N54uvMLq2ioMh4W12HCFE64l0VzG78iAZlfvJqNhPL1cRDoOFbWEj+Unil9gQO52T9r5tNtq/zmgjKyyV+yokGRDizly5oo/azcrilaz3CXM6qDNbOBHTk78Nn8ax2CSKrn8ph6X0b/PmnAgO4arZSnpFqSQDQnRUHg/s36+P+s/KovjwPkx4OWHrzb+j08gKT2V72EgcRv3GwWgLafMmvR8xhb/n/Q9JzmvkWbu3+fHuliQDwr8UA/829x7K0bhenI/qrn5+14Y0g4HNEVGkV5Ty6x690eRRgRAdw8cH/m3cqE8DvD7w76u9p5IVkUqBNc5vzVsTPp56g4XFFbt4OXaJ39rRXJIMiPalaXD2bMNJvGUL1NZCfLw+avfxx/WR+7GxvKcaQNiONoVH8+mSawyvq+aoPczfzREiMNXVwY4dN+/+OXkSjEYYPx6++U19xlBqKpjN/EkxgLC9VJlC2Bg2hiXlOyUZEALQK/5t3tyQAOTk6KV+p06Fp5/WT+IRIzpsta7D9jCKzRbSK0olGRCivWiaXsXvxnVj2zZwOCAxUb9mPP20Pv0vOtrfLfXpPxFT+VP+C/RwlXDF0s3fzVGSZEC0Pq8XDh1qyOCvV/xjwAB9/m1GBsyYAaGh/m5ps3gNBraER5FeWcaL8UkdNmkRotO7UfHvRgJQUKDP209L02t+ZGToBYA6yTm4MnwiLkxkVuzitZh7/d0cJakzIJQGjRzpM+Ysr7v591iPm2l1tUx31DDVUUc3r4cqDOw0WdlqsrLVFET+J577x0xMUx47qYkBe/VurzJecv6IMl7j8j0LIcxd0+j/xzuq+EPhBT4TP4AzVjsHcg8r9y1EoDOk/tF38KQDAKPmJdWbR4b3NBme00zw5mJC46SlJ1m2UWQFj2J70BDqjI3PVWOyTXns6OCByri3Qr29s0+lMh522XfMVNk4+NbZH2HWPDw08NcA5B+crdy3v0jPgLgjFq+XsY5aptfVkuaoYbDLCcBRaxDvhoaz3RbCnmon7k6SwTflYFAo5UYT6bUVnLH6LmQkhFBLdJYx132EeZ7TzPZkE00t5djYYBrEV60Pst6YQn73Ef5uZqtZEzWNX+S9TLSrnFJLpL+b45MkA6J5NI1kp5OpVVVMrapiQnU1dk2jyGhie7Cd18Kj2GmzU2pq+JVy17j82ODW5TEY2B4cTnptBX+IiPd3c4ToNGxeJ9OqzpFReZKMipMMc1zGg4H9xiR+Z55GlimFfcYkPAb/zBhqa1kRk/kFL5NRsZt3Yhb4uzk+STIgfKusJL2igunXE4CeLhdOg4FDdjt/6N6dzVg5ZbF2mud3d2uTPZJ7a8ro53L4uylCdFzXB/59++pGMipPklZ1lmDNRYElkqyIITyXsICNV/pQZmj7uf4dQYklig9Dh7OgfIckA6KTuDHw72NL/f7B7eai1cqW8HB2hoWxLySEWpOewX98zEAg2GcLpdpgZFad7zUNhAhIZWWwaVPDtSM/n18azGwPG8CTifeSFTGEE7aGUuFcDayEek3kVJ4ueI0Id5W/m+KTJAOB7mMV/9iwAYqL9eU3Z82Cl18m/aWXKAhqu3K/nYnLYGR7cDizaiUZEAHuRsW/G1/+H36o30wMHqyv9pmRQfQTZ6kzWf3d0g4hK3IqPy94hdkVe4H7/N2c25JkINDU18OuXQ0n8ZHrI+7HjoWvfEWfujNp0s2lfgtee82Pje14NtsjWVCcoxdOGjDA380Rov1cutR4qd+yMoiI0Of6v/aafu1ISrr59rqf3uVCYV3IVWsMB0KGsKB8h7+b4pMkA12dpsG5c/oJvG4dbN0KNTXQvbte8e8HP9Ar/sX5r2xnZ7LHFkadwUjwsmXwox/5uzlCtB2HA7Zvb0gATpzQu/nHj4fHHtO//MePB7N8jTTH2sip/ODyX6CqSu997WCkzkAHNyJpqDJ+rCD7ltfCNI1ZaGSgMd8EyR4NJ/Ch1cRmq5nNQWZOmI1oBgODY3so91+E70cEJbnq7vKEJmbgTZw4XhkvKVM/V8zasUkZT+41SBl3BikGPpZc8xl6qaqE+cOHwoEDyv0L4U/Gvr9RxjV32Sde0Ehxl5DhOE+G4wJpjjzsuLhkjCDLOoisoBQ2WgdRarw+8C+hiTFDIVG+YxVNfBl6IpXhHvFNVB0M76kMVznV1xZLE82z1RcrokW3fbVXzVU+3PAVeOcdeOgh9QH8QFK6LsCgaYwBMtDIwMskwAKcBTYGWdhkNbPLaqbGGBij/ttaljWY+QcP6mWVk5P93Rwh7liE18Fsx8WbCUCSpxIHJrYHJfHT0AVkBaVwwhQfMDOG2lJ+SDxHIvszctkySQZE6+nmcTPZUc3jmoc5aMQBVcAmDHwTA1kYuGgwEB0e7O+mdjnbrDa9ROp//gPf/a6/myNEsxk1L+Pq88moO8Pc2hNMdF7ChMYpczf+E5zCOls/tgclUWe0gDPG383tclYnTGLkmmX64mz2jlW8TJKBTsKsaYyur2WKo5rJjmpSXPUAHATewMA6jOyBLlPxryOrMRj156XLlkkyIDq8BHcFGXWnmVt3htl1Z+nmraXcaGNjUDJfi5pPlq0f+eYIfzczIKxJmMwTJ/+uj8G4r2PNKpBkoKO6PvDvoapSpjiqSa2vwa5pFBtN7LGF8mZYDHttIWy7csHfLQ1MS5fCww/rI6wTE/3dGiEaOBywYwe/LvmAuXVnGO66ihcD+4N68Ur4FLKCB7EvKAm3R11/X7S+C6GJMGyYfiMhyYDwqbKy8VK/Fy/yA+CjIDt/DI9lty2UbEsQmtz9+9+iRfoo6vffh298w9+tEYFM0+D06cZL/dbV8ZApnPXBg3g+cjYbgwdSagqMin8d3tKl8OKL+jTvDlTDRZIBf/J64aOPGk7i3bvB7Yb+/WHBAsjIYOrXf0id0ejvlopPiorS51cvWybJgGh/5eWNK/7l5elfLNOmwc9+BhkZ9FqUJQP/OqL774dnn9WXar7nHn+35iZJBtrb1auNK/4VFUFoKKSnw8sv68+i+/a9+fa6R2Uue4e1dCl89av6v2FsrL9bI7oyj0efyvrxin8eD6Sk6N3NGRmQltZ4UJphvf/aK3wbOhQGDtRvJDpQMiB1BlqBLWSIz5hF00g35jDT7WGm28MIjxeAIyYjm80mtEH9uBAeisfH3f/RslrlsfcdO66Mx8WoR6z2TRqjjF8uyfEZi65Q55JffO5byvg9qaOV8X+8809l/MdvvKmM9zEnKeNRA33Pg3aVKBYsBy5fddFN83DSUcj3LRH83dz451xcq95eCIDYuX/3GYt3lDIuZyvzKs8xu/o83Tx1lBttbArry7qw/uyOH0e+rZvP7bXIQuWxq6vrlXFLsfraYbH0UsbDwmt8xpwGdY0Rq12975S4HGX8Url6FlV5tbq31RSuvu4a8b3GQJVFuSlxNvjeoT/xqezVTHzwX7iNja+j59+crd5BG5GegTbQ1+tkrqeaOZ4a0rw1hKJxzWBgq9nEH4IsbDWbKL7+5f9AZLifWyvuVInBxG6jlYUexy3JgBAtFeRxMrE8m5klx5hZcowh1QX6wD97Aq/EjCcrrD8f2nveXOo3TJEIiI4tq/c0/vvYu4y/epTdCeobsvYiyUArCNU8zPDUMserJwB9NRdOYI/Rzq8sMewKqua4ySgD/7qglSYbz7sqidC8VBhkbIdoAU2jf81lZpUcY1bxMSaVncbudXIlKIot3YbxYp97+Y83nFJJNLuc49EDyA+NZ37udkkGOrVPDPy7UpeNBThvsJBlCmWDMYTtJjvV1zP4UHNgLfUbSNaYbPx/rkrmeRz8Uy7aoimfGPi3Jy+PeoOZPVGD+J9+S9jSbTinQnveHPhXWpjj1+aKNmIwkJU0hcUXNvH0hMfwGk3+bpEkA8127VrDwL/16xsG/s2axfcs3dloCuWCUZbrDDRXDSb2Gi0skmRA3I7HAwcPNiwU9omBfw9tN7M7KoU6U8eZYibax7re03nk5DLGFp5gf/wIfzdHkgGfnM7GS/0ePqy/PmYMPPJIw1K/Viv/pxhAKLq+VSYbT7qqCNW8VMujAnH5cuOlfktL9aV+09PhD3/Qrx29ewOwSTGAUHRth2NTuBrcjXl5OyQZ6HBuLPWblQVbtkB1tb6079y58L3v6Uv9du/u71aKDmaV0cbPqWKOp573zbIWRMBxOGDnTv3OPysLjh/Xu/lTU+HRR/Uv/wkTZKlf0YhmMLK+91Qycnfx89Svo/n5RqJL/HZGNVFZyxp8++UuQzQv0zz1THVVku510xd9qd+9BhObjCY2mu0cK6vF8MF6+MD3nN2YUPWCHn3GTvMZ69YzQbltWoFqqUwoyD+vjNcHq7sfzx06poxryb4/26QR8cptw73KMIdz1FOfjufmq3egqU+ea9Xqn018xCyfsbAw9e/UtatHbv79ihE+MphZ5KnlA7ngdyoh92xUxo23W+pW0xhYW0B6yUfMLvqQqRXXB/5ZI9kcOZxNg/6brZFDKbWEYVl1FVbtBnbfdv+R9FMev97hO26O9Si3DS/yPbUPwOWwKeOmevUS5Q7T7ZfqvcHmHOwzFhOk3rfdpH7k6tLcyngk6vPXFqW+rhb1Ve/fnBPqM5ZYr1i6GXBVHbz593URw/hc7QpGnNvLkcj+yu3aWkBduQyaxjCvi5keBzM8DsZ7nfrAPwxsNJrZaDSzw2Ci+hOj/mUOgGjKapOV77prCdY06mTWSJcT7q5hRukRZpccIr3kEEmOIuoNZnZHDOIXSUvYHDWck/aeUvFPtMj+6BRKrOHMv7JPkoG2FqtpzHHVMNNTT5rHQSxeqjGw0xTEk9ZItphsnPKoM2ghmrLKFMST7lpmep2skcFgnZ5R8zC64ixzru5ldskhUivPYNa8ZNsTWRU7kU3dxrAzahj1LpkpJO6c12BkffdxZFzdx69SPuXXZLLLJQMWTWM8XtI1D+l4GIEGzjqOGi28YwlhiymI/cYgXB//oat724Ro0kWjiRMGEws9kgx0VvGOYuYUHSC9+CCzig/SzVVFhcnO1uiRfHfQ19nYbQz5wXGNtjFJMiDu0rr4VD6Vv5khlbmcjEj2Wzu6RDLQR/OSjodZmpepeAgDCoEtmHjFYGJXcDeKDP6fxym6tlWmIL7ursMqFb47B4eDmcWHmF10gNnFBxhWdREvBg5FDOT13veyIWYcB+198HSAOeCi69rbbQjllhDmXd0nyUCLVVXpo/2vj949pDlwAR9i5EWDhY2YOI7hZsU/qyQCoh2sNll53F3LdK/L300Rt6NpcOZMw4yhrVtZVVfHlaBoNsaM49f9PsWWmDGUWCNubnLbAYRCtCK30cymuDHMu7qP3wx8wG/t6BzJgNerz/O/cRLv2qUv9duvH8ybx6dffYMd3DrwT4j2dMZg4qzBxEKPegEY0Y4qKhov9ZubC1arvtTvs88y/oNIToQly8A/4Vfr4sez9NIOBlRf8lsbOm4yUFioV/pbt04v3FFYCCEhMGsWvPSSPne3nz7tZu1rb/q3rUIAGAysNll52O0AlwssTSxfJlqfxwOHDjXM+d+7V39t0CDIzGxY6jdEn3p2Yqt6aqEQ7WFnzDCqTTYyru7zWxvabQnjkQMnKOP1bi9jHNVMq6tgWm0Fw5z6EpLHrXYOxcbzYWQMx8KicN9mqd/BY4Yp970gUb1c5mFPrjLeLUa9Vv24keq5winjUn3GLr3xb+W2+69kKeO/eON9ZTy/Tt1l7SlTF8lJmTDRZ6xn9ybm0+YeVcbrQ5OV8dBgdRdtTbl6mdFzZSeV8V7dfM+D7tFLPc3n8Natt319qNvJuorrieycOcp9iObpN+ofyrjL42RO1WHmVB5iduVhYjxVVBjtbA4bwZ4+k9kWm0pByO1rYkRY1L9DRU0MLg5xqefDOxMqlfGYUvWqpW7FRKeiMPVc+bDz6llSpRZ1nYCKYvXjVWN39WePVZRfiQ9tYvljk7oin9ekHrhZ5FSGcUaq66u4agco456rvn+2LqP6HtsZcvtt/7TnVwyuyGV4+QXl9m3Frz0DPZ0OJtdWMKW2gnG1lYRqXkqMZnbYI/hzRDw7gyMoMVvoHi9LdYrO4YTJQq7RRO9lyyQZaCNWr4txNWeYXnWUaVVHSXHk48XAQXt//hg7nw3hY/gwZCBug5m4pDB/N1eIZlnRcyoP5m3VK+H2b/+aA+2aDNi9HlJrK5lcW8Hk2gqSXPW4MHA4OJRXIxPYYY/gpNUuS/2KzstgYK01mK+9/z688gqYZPDqXdM0+tRfuf7lf4wJNaewe+u5Zo5kR9gIfhb//9gUPooSs/ouW4iObGOPcdSagrAvWwaPP97ux2/bZMDrhSNHICuLNwpOMaquGgsa+ZYgdtkj+F97BPvt4dQaTVTLXH/RRayxBvO1wkK9Xn1amr+b0zl9bODftlP/oaermHqDmQMhg3i5+xK2h43gjK0XGAzkmeTuX3R+tWYbG3qMI/O997pIMnBj4N+NpX6vD/yrMVr4dWwSu+0R5FvVNbGF6MwOm62QmAjLlkky0Fxeb8NSv1lZsGePPvBv4EA2RoxlR9gIPgxJoc4k1w7Rda3oOZXMvf+jz3q5vrJle7n7ZMDp1E/cG6N3P/pIf330aPjiF/XRu5Mn861hvhfrEaIr0QwGWLpETwZ++1u4zaBXAVy50nip35ISCA/Xl/p95RX92pGczM+aGEAoRFexLmG8PvX1P/+B73ynXY99Z8nA+fMNJ/HmzfpSv7Gx+lK/3/mO/l9Z6lcEsqVL4eWXYd8+mOh7RkZAqa/XH53cuHYcParP7x87Fr7+9YalfmVKpghQVZYQ/ftz2bIOngy88AK8+qqeDJjNMGUKPPGEfhKPGiV3QELcMGUKxMXpJ3WgJwNXrsCXvgTbtkFtLcTH69eMH/1In3ERo14CXIiAsnSp3qt+5Qr06NFuh212nYG/PvxfPPy3tyiK6ca+8akU9EzEZW2YZ2oLUycCJ6rUz/qSEvv6jMX3Us/zv2+uuoSjTb05pz46rYx7HOo7lWGT1HUGVDR3njL+79++pI6/u1MZzwnxve42QGyoegT2cEUNBe2KOpf8+46/KePOa+pKfQ/812eU8e5NzDp5L+svyngFvqes9ohW/5ueP6+uYVBachm+9jV93Mz58wFd4e6pyPt5rmIZP4j8FFnBIzhm6dXo5zE4Nkm5fWhUgTJegO9kIila3UMZHFKijGt29flhNainPVeSr4ybnE1cnMp91yk4UbRfuaktRL1g1jVPtfrYB5tYzXWAulZAUIzvn80oi/rnGjc0WRn3VKv/3WrL1T93d2SkMl51Tv2zKwzyPUvIastRbmsPU9dnOPnsZL1n/aWX4L//W/ne1tTsW/nq0BBqg23k9+rFxb59GiUCQojbWLoULl5sGEcToPLM3fBi4G+h0zhmTQroxEiIJkVHw8yZ8N577XrYZicDmtHIxT7J9Ll4sS3bI0TXMWOGfmIvW+bvlvjVquAxaMCi2kP+booQncP99+uP1YrUVSJbU4se8l/s04eo8goiy8raqj1CdB0Wi14Pf9kyfcW8AFViCmNn0CAy6w76uylCdA6LF+v/XbGi3Q7ZomQgv1dPXGYzfS9I74AQzbJ0qb5s7kn1GIOubrl9HHPqjhPilSWBhWhSXJy+smY79iq2KBnwmM3kJfWi74WcNmqOEF3M7Nn63PkAf1SwIngsNlzMdagXiBFCXLd0qV6Fs5164ls8F/Bin2S6X7uGvaaJkaZCCAgKgoUL230wUEdz0RLHMUtPMmvlUYEQzbJkib4U+sqV7XK4FicDOcnJaAYDfS6ql/0VQly3dCkcOwZnz/q7JX613D6OhXUfYdJkIRIhmpSYqNcoaadexWYXHfrG7/7Q8D/nLjDTHsLMj732j1/+VLl9r9hoZdxQ63s+7dG9Z5Tbpo+YrIzbYhOVcbNZXUfgrxvfUsafmfRt9fGJ8BmrvFah3NZQrl63O+3escr4fT3HK+NnL19SxnOzfY9mTYxX//pEedRTyE5Vlyrj2YXlynif4aOVcZNZXWPBcNn3evAnqtX/LnUelzJutzUcO1jTyAN+NXgoL5j1Kbm1jibmeHchWs6n9b8cGAipy3H/tZc+0+K6UfepHx04+qrHGcRdUNQwsajrAJgj1fPdQ2vU89k9wVeU8RCbOvEpt6lHiztrcnzGwuLV1zWv5aoybs9V136xDGji68Gjbnu0M8VnLDZG/XMPcqi7xouD1XUAor2RynheEzPjPaHq1UaDe9h9xpw56t85U3CVMj706w21Yx42jOVbq19n6pfWU2vVj3ni1anK7e/UnZUMzMzUn2VUqT+UEALqDAayjCYWewP8jnjsWP1uZ/lyf7dEiE5hY5/pBHmcTM/b0+bHuvNkoL5ery8uhGjScqOZsZqXXprX303xH4NBv3asWBHQUy2FaK5L4QmciBnInAvb2vxYd5YM9O0Lw4e36xxIITqzdUYTDmCxJ8B7BzIzISdHH0MhhGjShr5pTM/bg83VttNy73xlocxMWL1aH+0ohFCqNhjYZDSx2Ov2d1P8a8YMfaqlPCoQolk29J2B3e1gSsG+Nj3OnScDixfr8x937Gi91gjRhb1vNDFJ89IjkB8VWK2wYIH0KgrRTDmRSWRH92XOha1tepw7TwbGjIGePeWkFqKZ1hjNuIB7A30gYWYmHDoE+eqV5YQQug190piRuxuLx9lmx7jzZODGYKDly2UwkBDNUG4wsNVgItMT4I8K5s/X122QGwkhmmVD3zTCnDVMLGi7ol3NrjNwW5mZ8MorcOQIpaXqioRn8tVFio5m+45frjin3Daoe29l/CsJX1XG1765XBk/dGKvMr7+g7XKePr0aT5j2ReuKbftltJfGY+6Vq+M51ep5/J7StW1AFat+KvP2ENLP6XcNjI4UhmHy8po1UcHlPG9NdnKeG2pOov2hvn+7OZ6dVe+BXXc67l9gvy+wcjvvC59NbLYJtay76oiIvQlWlesgG98A5M7WPn2xJpkZdwR4fvfuTAyRLmtLUR9P1Qdrb7RMbvU1yZ7rHouP+XJ6nhYmM9QfLm67YYmxpv1sqt/NjlNzMUPb2IwbGSC72uPNTFKua07XF37xVatHqvm8AxQxl1NDMYzVBeqj1/i+7prcqhrFHgqE5Rx7Ldunx0+kpzwJOac363e9i7cec8AQFqaPhhIMnwhmuUDw/UTPdAH0GVmwtatUF7u75YI0fEZDGzoM5v03C1tNmj/7pIBqxXuuUcubEI0U7HBwE6DMeAXLuLee8HthjVr/N0SITqF9cnpRNZXwLa2qTlwd8kA6Bn+4cOEVfguJyyEaPAfg6ldVyPrkHr2hHHjpFdRiGY6ETOES6E92uxG4u6TgeuDgfqcP98KzRGi6/vAYNLvittpNbIOKzMT1qzB4lGPexFCoD8qSJ4N778PbVC87O6TgfBwmDWLPuckGRCiOa4YDDB5sjwqWLwYqqsZV/Khv1siRKewvk86XLsGu1t/IOHdJwMAmZkk5hcQ5GjbcolCdBlLl+prewTyYl9Dh0Lfvsy4ttHfLRGiUzgSNwISEuC991p9362TDNx7L0ZNo/eFi62yOyG6vCVL9MW+Vq/2d0v853qtkrSrmzAEclVGIZpJMxjhvvvgP/8Bb+ueM3dXZ+CGxERyY6KJPXGC1XGRt31LQVW5chfXrl7yGatp4u7p3T/+nzK+a8M6ZbykQL12faVdPdBr+4fq+fJmzfd8+FOnTym3tXnV9RtiPLXK+CWHusBNjFU91717P99rppfmVSu3jYxVrznet0ydi3oM6uQy71qkMt49po8yboj3Pc+66JK6PkP+eXV9CFdYE5O0k5P1AXTLlsFDD6nf25UtXkzsiy8yynSYU3EjbvuWxDj1+vA5hWd8xuKs6t5Kp6efMm6x+v79B4i1hqv3b7uqjOMqUYdrfdcpiHSpnxtbQ1KU8ZKwHGW8m0s9KNwSpK4VYE/u6TMWHamuv+AqV8/zN7qGKuMOq7q6ZY/IOmX8fKW6BkPEJd/nd4VDXd+kwqW+ZnuaqLHA/ffr9X3274cJE9TvbYHW6RkAjvVKZPClK5gCfVU2IZpr6VJ9al2t+uLQpU2eTHlQJNNyN/u7JUJ0DtOm6QXLWnnMUaslA8eTErG53Qy8or5jEiLgaRqcOqUPBKqthc0B/EVoNrM7aQbTc2XcgBBNsbkdsHYtREfDBx+06r5b5zEBcCUygqKwUIbnXeJUzybKLQoRYMyaxlR3PQvdDu5x1cGQIRAcDIsWwbBh/m6eX23vnc6Cs8tJrMjlUoS6tLgQgSauppg5uTuZm7uTaZf2wxv10K8ffPnLrXqcVksGMBg4npTImAu5/Fsbh2ZQ17wXoquL8nqZd/3LP8PlIBKNSwYjayzB9F/2LsyaBXa7v5vpd/sTJ1NvCmJa7mbeHfEFfzdHCP/SNIaVZDM3R08ARhWdwmMwsr/7cF4Y+wg//cujkJKiD8BtRa2XDABHkxKZeeIMScWl5MaqB/0I0RUN0Nws0JzMrypnsseJGThksvByUCirLcEcMlnAYODLCxf6u6kdhsNiZ1/iFKbnbpJkQASkIE89U0oPMadwF+k79pBYU0iVxc6WXpN4Y/iDbO41idLrC7/9dPDgNmlDqyYDF2NjqA4KYnhegSQDIiCYNI3JuFig1bNAczIAD3XAZoONbwZHstoSzGWjehUzATt6z+LxnU8RWVdKeXC0v5sjRJuLrS8hvWgPcwp3Mb3kAHaPg9zgHqzpN4P1yVPZ22M0LlMTMwtaUasmA5rRyIleCQzPu8SqsSNbc9dCdBjhHg8PeB0soJ65mpMoNK5iZK3ByhOGULZgpSJUPa1SNLYraQYGTWNy/lbWDFzi7+YI0fo0jZTqs0wt3crswl2MqTiFFwMHI4fy234PsyF2CtmhyXj6NDEtuY0YNE1TL9jdUitW6GVGz5yBgQNvvvylmXOVm71//COfMZdTPV0xOkp9JxFqU1+YXXXqqV1VLvWPKLZbjDJuj/C9XnthlbqGwczBg5TxtNGjlPGqqiJl/GKlejnM48fP+oxFhaj/XeJD1HfE4fXqySyVQeo53tXB6uOfKVTXvK87e9RnbH/2J+o/nD2rryWwciXs2KHXBh81Sh8AuGgRjB0LxlabnBOYpk3TR0l/YvGih55Xz9WvKfdd5yOkh3o++9WSZGU8JExdp6A2TH3tsTdxXbfb1PUsrp73PdffWqWuIVLSREHYkiZuOnuUq+8VY+PVz6xjw7r7jtWXK7ctjlZ/troqdX2Uqlr1da/wirp+SynquL3K9w/v2LLhDf9TXw9btujXjVWrIC8PQkMhI0O/bixYoE8T7ABatWcAgDlz9FHSK1bAD37Q6rsXol243bBnT0MCcPo0BAXpg/5+9ztYuBB69fJ3K7uWzEx46il9uqUMrBSdVWGhXll05UpYvx5qaqB3b33Z7kWLIC1Nv5Z0MK2fDNjtekIgyYDoZEI8HibV1jCtphq6d4fSUoiL07/4f/lLmD1bz+pF28jM1K8Z69frvYtCdAaaxoDKs6Rd2caMq9sg/pj++oQJ8MQTDdOHO/gMu9ZPBkA/kb/0Jb2oSnffXUVC+Fuiy8m06mqm1VQzpq4WM5BtDYJvPqafxKmp0v3fXgYM0Osv3HjUKEQHZfE4GVd8gBlXtjH96jZ61l6m1hTM7u6T4U9/0rv/O9l3X9skAwsX6lnQqlV6UiBEB2HUvAwpK2VK4WUmFVygr9OJ02DgYLCd38R2Z0dIKFctFvb//Of+bmpgysyE//s//TGNuW0uT0LciWhHGTMKdjKrYAfTLu0l1F3DleB4tsVPZ2uPGeyPTcVpCuLYF4Y3vbMOqG3OtthYfb32FSskGRB+Z3e5SC2+ypTCy0wsukqks54yaxA7g4J5rVssH9rt1Mr0v44hM1N/JLN7N0yf7u/WiECmaQwov8Csgu2k5+9gdNFRjGgc6TaUNwd8nq090jgTMajDd/83V9ul3osXw5NP6oMnQtQrQAnR2rpVVbMkJ4cphZcZVVKERfNyITSClb36sDsugVOR0dScO+7vZopPSk2FHj30GwlJBkQ7M3ucjC48SGr+emblb6d39SXqTEHsTJjITyb9hC09p1Jkj1XOJuis2i4ZyMyE739fHwx0331tdhghAAxeL32KSxiZX8CI/AJ6lpfjMhg5HB3LK4NHsieuB1fsMvivwzMa9VHXK1bA//5vl7nrEh1XhKOMyZd2MK1gC5Mu7SDUVc0Vexxbek7juV7T2BOfSr1ZPUW1K2i7ZKB/fxg6VD+p77uP/kPUgyliisN8xpzqKZ+EBJUr48Y6dc+Eu4lll2vd6gvS+ZyDyrhFUabAFhun3LaiWB2vqVevnR0a5bvGAUD/0kvKuCPK93xfb1gT+3ZEKuOMVk/Cfvapv6i3r67Wk82VK2HtOigqgm7dbs79t2RkkBoeTirwbfWeREeyeDH88Y9w4gQMG4bJ5lW+vVvyKJ8xV1y5ctsBRvV89BBDpDJeG64MY1SfnhQ51HUK+jiqfMbON7XydZ36CywuUT29rdqkvvD2tjZxd1zv+9Fbvl197QiPUu/aUqv+wb77y56+g5qm18G5MW141y7wevVaIZ//PixaRI/Ro/m0wcCn1c3oUtp2hE5mJrz2mj4YSIjWkJfXcBJv2QJOJwweDF/4gp4ETJoEJnn+36nNnAlhYfqNRICv6ChaicsFO3c2XDvOnQObTZ8u/Oqr+qD3hMBebbftk4Ff/ELPvIS4E14v7N/fcBIfPaqPMp8+Hf7nf/QEoF8/f7dStKagIJg3T08GfvITf7dGdFZlZbB27fWew7VQUaGPR1m4EH7zG0hPl+JWH9O2ycC4cXq29YnyokKomJ0uep2/os9EWb1ar1cRFaXP3X3iCf2LIiLC380UbWnxYvjMZ+DSJUDGDYjmiSs+y8jT62DGer0nwOOB0aPhW9/SbxzGjJG6IT60bTJwYzDQ8uUwb7IMBhI+hVbU0PtMAX2yC0i8eAWz2wuDrsB//Zd+Ek+eLPPOA8mCBfq/9wcfAJn+bo3ooIweN/3z9jLi9FpGnF5LfPE5XOYgmJsOv/+93gvQUzF+QNzU9lfX6+MGYkuHUNStiVEhInBoGnGXikk+U0DymXxir5bhNRq4nBTH3vQx5AzqyWdfet/frRT+Ehmp13BfvhzmSjIgGgTXlDPs6GZGnF7LsOwNhNSVUxEax9FB81iW8TNO9Z/B75+J93czO522TwauDwYacLFAkoEAZ3a76XntGsmXL5F8+TIhDgcOm5W8AQl8NHUYef0TqA/ueAt4CD9ZvBi++12Cp1VSF9zEsH3RpcVeu8DIw2sZcXgdA7L3YPK4yY8fzpYJX+ZoynxyE8egSff/XWn7ZCAoCObPZ+CWjewe1znLNIo7F1JbS/LlyyRfvkzPwmuYPR7KQ8M4m9Sbi5N7czUpDq9JTmJxG/feC489xsgzW9g7SnoHAonR46bvuX2MPLyOEYfX0eNKNi6zlTODp/Pup3/Fsfj5lEbKqqGtqX0ewmZm0uNf/2J4qJm6iNsXfhkV5bvX4Ii3WLn7YIt6vqu7rok1yYMVhQCApqbThoX6rpEA4KlTjJVwqbc9U1WnjIfs3KOMJ4So6xR4LWXKeFyS74F6T714m4GhmgaHDjWM/j90SJ/qN2UKfPc7sGgRkYMGMQoYpTyyCHhJSTB6NCPObmJduu8Z38neyz5jzmx1oumKU48mL/KqB6o2VYrG4lLH7b7LCOj77+V7rXurQV0fxWq7pow7UW9vMaq/Hurt6k/v9vquD/H+j27zc62ogKws/brx6zUNq4becw8s+hWWOXMYFhqKTDZtG+2TDCxYgNdoJOHURc5PlN6BLqeuDjZt0k/iVavg8mV9tP+8efDd78L8+RCtLq4ixG0tXsy4X72A2e3EbVYXqRKd0IULDTcO27bpNWmGD4evfU0fODx+vIz+byftkwxERlLUJ5GEUxckGegiQmsc8Prr+km8caOeEPTrBw8+qJ/E06ZBEz02QjQpM5OQp59m2JldHB4609+tEXfJ6PUwMP9D+NH1HoCTJ/XrxMyZ8OKL+uj/5GR/NzMgtdtcrUuD+zBq7U7MjnrcNhkk1uloGvFFlQy8cJUBF6+ReK0c/rRRn/L3zDN6ApCSItNHResaMYJr3Xox6dBqSQY6KVt9FaPPbiD19GrGnMkiorYYYmL07v/nnoO5c/WKk8Kv2i0ZuDy4D2NWbadHdi75Iwa212HFXTB7PPS/eO1mAhBR7aDeauZc7zj2jerDfW+t109qIdqKwcCeMfcw5cBKXv3sryXZ7CTiyvKYcHILqadXM/TidiweF3lxQ9g47vPsT7mHX/1hmpQN72DaLRmoiwyjLCGWhJMXJRnowEIdTlKulDD4cjEDrpUR5PFQFm7ndP8eZPeNJzex283R//dJIiDawYejF7B4w2v0zz3COcWiRMJ/jF4PA/IPMf70OlJPZZF87RQuk4UTydP467xfcSBlPtei+zZsIIlAh9OuJd0uDe7LoJ0fYXB70Mzyy9AhaBrxFTUMvlLMkMsl9CytBCCvWzibh/QmZ1QSRdFhckcm/Ob4wMlUhUQy8dAqSQY6EFt9NaPPbiH1dBbjTm8gsqaYSns0BwbN5p/pP+BIv4XU2aQ+RGdh0DRNPa+uNR09CiNH6tNH5s5tFPpy5lSfm524cFW524Q49XKYpTXqKTTHrxUq4x6XetVFm0OdUxVV+Z4/5HI2scYpkcpoYu9uynhBzrlbX6yvh61bG0bx5uVBaChkZOjP/hcsgFjfU5qEaHef+xwcPqxfQz7hoX/43syrnjmLI0a9PHJInnoku7GpCQ5NLKFhKVHHXYrmlXvUyy+X56s/W5hN/dnWvXCba0Benj5jaOVK2Ly5YdXQ68uGy6qhnVf7FnsfPlwfKbp8+S3JgGhjhYWwZo1+Eq9fD9XV0Lu3Xthl0SK99GuQDOwUHVRmJvz97/pUtL59m36/aB1eLxw40HDjcOSIrBraRbVvMmAw6CVG//1vfREJmT/adjSNQS4n/PKX+km8d6/++oQJ8OMf6yfxsGHS/S86h4wMPVldsQK+8x1/t6ZLC3LVMDp3OzyyXV819OrVhlVDf/xj/d8iMtLfzRStrP2XgcvMhN/+Fg4ehNTUdj98V2bRNCY66phTV8vs2hqSPG54/nm9F+ZPf9JP5u7d/d1MIVouNBRmz5ZkoI3EVF1m/IX1TLywnlF5O7B66mHQIPjsZ2XV0ADR/v+6U6fq1ehWrJBkoBVEeTzMqqtlTl0NaXW1hGkal0xmNgTb2WgP4a3cC2BrqmiqEJ1AZqZema64WKa03iWD5mXA1aNMuJDFxAvr6V94DI/BxPHEibw59Qn29s3gz3+e4O9minbU/smA2axXmVqxAn7+83Y/fOenMRgvi3CxCBeTCsoxAR9Zg3g1PIoNdjunLNaG7n9JBERXsWgRfPWretf1ww/7uzWdTpCrlnGXdjIlZz1TcjYQU3uNqqAIDiSn8964RzmQPItqW6S/myn8xD/9PpmZ8Le/wfnzMvikGSxoTMN9PQFw0w8vNcAGzDweHcsmu50ik3ThiS4uPh4mTtRvJCQZaJaYmqtMztnAlJz1pBbsIMjjIC+iLxsGLOajgfM4kTAej0nKhgt/JQNz5zYMBvrud/3ShI4uGpgPLKKGebiIAAowsBILK7GwBTMODCSGyTxeEUAyM/UStnV1EKyeUhyQNI1+l48x/lQW445mkVJ0FI/ByNEe43l9/A/ZlTyXvKj+QNNTC0Vgad86Ax+3aBFUVuorVQEvPPZ5n2/ddfKAcldRvdSD4qrO1ivjmw/cOnf549xB6q52u1W9DOrl4hxlHNCX/j1zpmEKz65d+rSeceMa5vCOGiWj/0VgO3NGXwPjgw/0cwKY/4rvt9tcter9hahrkNgvqAsJhDfxfWppYhak86r68vvaT5rxmK+uTp/zf2PV0EuXIDxcXy100SJZNVQ0i//6ljMz9ed/gTwYyOWCnTsbEoBz5/S7ndmz4bXX9IU8EhL83UohOo5Bg/Q/K1bcTAYC0tWrDcV/Nm6E2lr9kesDD8iqoeKO+C8ZWLQIvvIV/Rf685/3WzPaXVkZrF2rn8Tr1kF5uf6Fv3ChvoTnrFlgV/c0CBHQFi+GP/8ZPJ7AqXanaXrBnxs3Dvv363VaJk2Cp57Sr6eDB0vPobhj/ksGunfXf5FXrOj6yUB2dsNJvHOnfhEbMwa+9S39JB49WgowCdFcmZl69bu9e2HKFH+3ps2YXQ5Yu6Wh+z8/X1/qNyMDvvENvW5IoPaqilbn3yHomZnwzDN6F1cXYtI0Ut31zHXWMcdZq3dr2myQng6vvKL3AiQm+ruZQnROEyboNxMrVnS5ZCCs+hrDz6xjxOnVpJzfDM4avYT74sUNZcOtTS2IIETL+TcZWLwYHn9cf+bVyYV7vcx06V/+M10OojQvhQYjG612+v3zLX0cQEiIv5spROdnNOpraixfrvcQ0Im7xjWNxGvHGXF6DcNPryb5kj5Y+mLP8axNe5z7fnMfDB0q3f+izfk3GRg4UB8ZvGIFdMLH5H28Hua761lQV8kEdz0W4ITJwpu2MDZYgjlitqIZDHw6M9PfTRWia8nMhNdfh9OngcH+bk2LmF319M/eztAjaxh2aDXdyvNwWEM52T+dv4//I8cHZlAVGgfAfcOkaJhoH/6vVJOZCX/+M4YH5qF18OfmJk1jvMfFPLeTeR4ng7we6oFdFhtPhUSz0RLMJSn+I0TbS0/Xe9qWL4fwjp8MhFYVMeToOoYdXUPK8Q3Y6qspje7F0UELOJpyD9l9puM2y6qhwn/8V2fghr17YdIk3v1/n+JSz563fcv2E+o6A2FNVNCKc6rX9f7N+g2+gxUVkJWlD+JZswZKSyEuTp/2t2gRzJmjL6IihGhf998PBQU88Km9Pt9S7yxV7sLEVWXcWZOsjK9+xkeXpqbByZMNA4f37NFfHz++oW7I8OHS/S86DP/fxo4fD/Hx9Dt/zmcy0O4uXtRP4A8+gO3b9XoAw4fri6QsWqS3uYP3YgjR5WVmwuc+R+SCK5RH9PB3a8Dp1K8XNxKAixf1acJz58Ibb+g3ELJqqOig/J8MGI2waBH9//1vtk9P80+m7PHAhx82nMQnTugjdmfMgN/8Rk8Aevdu/3YJIXy75x4wmRh3/AM2Tvmqf9pQUqL3GK5cqfcgVlZCz54Nd/8zZ8piYaJT8H8yALB4MVGvv050aQml3dpn3myQ282gkhKGFBXrC6DcqIR4zz3w7LN6Nh8W1i5tEULcgehomD6d1OMr2i8Z0DR6FZ9hfPYaxmevhZ/t1cuGp6bC97+vJwAjR0r3v+h0OkYyMGsWTouF/ufOsa8Nk4GoujqGFhUzpLiI/qVlmDWNKyEh8Nhj+kk8YULgVDQToivIzGT4936IzVGFw9Y2ybvJ42Jo3i4mnNETgISyCzjMwRzuNwv++Ef9BqJHB3hMIcRd8P8AwuuyBw4irKqKf3zms7fE7nQAoUHT6FtewdQrhQwtKiahuhq3wcD5qChOxsZwIjaW0uBg9QBCIUTHlZMDffrwwuf/xd5RD9wSvtMBhKG15Yw5vZMxx/Yz9twGQusrKA5LYN/A+ewbuIAjfdJwWoJ9DyAUopPpGD0DwLn+/Vmwdg0h1dXU3MXo/CC3m+HFJYy+VsSowmIinE6qLRZOxXRjQ98+nO7WjXpzh/nYQoi7kZxMTsJIUo+tuG0y0BIJhRdJPbWV1JNbScn9CJPXw9keo1k+8RvsGzSf8/GjpPtfdFkd5lvxYp++eA0G+p0/x9GRo1q0bXSdgzHXChldWMSQklIsXo1LoSFs75nAR91jqQgJR5OTWIguaf/wTOZvfxmTx4WniWnGH2fyuBias4/Uk1tJPbWVhOJc6s1BHB0wkT8ufpIDg2dwzTy+DVsuRMfRYZKBgtJrXIiNJenEMTbHRDaKrVi7ufGbvV44cKBh9P+RI2A2w/TpN0fxJvbrRyIQwIucChEQdo1azANZz9E/ZztHB6U3in3wrejGby4r01cLXblSXz20vFx/3n/fQli0iKD0dFLtdlLbr/lCdAgdJhkAONUzgfmHjxLkclH/ybW4a2r0NQxWroTVq/X1vKOi9JW7fvxjfSWvyEi/tFsI4T8Xeo6iMCqJCUdX3JIMAHD2bMONw44d+lTi0aPhm9/Ubx7GjJG6ISLgdahk4GRiAosOHWbglascS+pFeG0tKZev6KN1N28Gh0NfAfCzn9VP4smT9R4BIUTgMhjYN/xeJh5dzuv3v4TR62Hwxd2kHlsJv1sJZ85AUJBewvj3v9dXDe0oBc6E6CA61DdpeUgIlyMjmXniFGmnTpNYVo7HYNCX7Xz+eT0BGDDA380UQnQwe0cuZuH23/OT/1vM4PM7Ca8tpSysOzy4UF/ZUFYNFUKpQyUDAAf7JDPrxEnOxndnR8ogzvSI5+lly/3dLCFEB3ai/3Qux/YntjSXtdP/m33DFnEuaRwrHpPufyGao8MlA7sHDWD3ILn7F0I0n8dk4WtPn/V3M4TotCRtFkIIIQKcJANCCCFEgOsw5YiFEEII4R/SMyCEEEIEOEkGhBBCiAAnyYAQQggR4CQZEEIIIQKcJANCCCFEgJNkQAghhAhwkgwIIYQQAU6SASGEECLASTIghBBCBLj/H8/MqkcTDJNRAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "ds = ds.where(ds != 65535) # remove no data values\n", + "visualize_sat_images(ds, gdf, list(reversed(arguments[\"bands\"])))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### special arguments for cdse" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "# filter_asset_path: filtering asset paths with regex\n", + "# sometimes the same file names exists with a different folder path, for example for Sentinel-2 the image data is in the folder IMG_DATA\n", + "# defaults are already set for COP-DEM and Sentinel-2\n", + "ds = tg.create(**arguments, num_workers=4, filter_asset_path={\"SENTINEL-2\": \".*/IMG_DATA/.*\"})" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "# use_virtual_rasterio_file: downloading the tile files instead of using the virtual rasterio file - default: bool = True\n", + "# in combination with rm_tmp_files=True, the tile files will be deleted after the download\n", + "ds = tg.create(**arguments, use_virtual_rasterio_file=False, rm_tmp_files=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "# resampling: using a different resampling method in the reprojection process - default: rasterio.enums.Resampling = rasterio.enums.Resampling.nearest\n", + "ds = tg.create(**arguments, resampling=rasterio.enums.Resampling.bilinear)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### filtering meta data" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAFeCAYAAAAYIxzjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABDs0lEQVR4nO3deZxU1Z338c+trav3nW66odn3xQVQFgFll0XQ1uTJM6+ZSTJZxzExyySZxBg1cTIZk4xOHJPMk0kyi0km0SgiArLJoqIgouzQrL2w9AK9VVd3Lff54yItigp9q/tWdX3fr5cvY1H87g/S3fd7zzl1jmGapomIiIgkLZfTDYiIiIizFAZERESSnMKAiIhIklMYEBERSXIKAyIiIklOYUBERCTJKQyIiIgkOYUBERGRJKcwICIikuQUBkRiZMOGDXz6059m5MiRpKenU1paytKlS3njjTfe996dO3cyZ84cMjIyyMnJ4Y477uDo0aOXvOfQoUN8/etfZ8KECeTk5JCXl8e0adN46qmn3levqqqKe++9l5kzZ5KTk4NhGPz2t7+96j/DlfQF8Oijj3LHHXcwaNAgDMPg5ptvvuprtbS0cO+991JSUoLf7+faa6/lD3/4w/vet3XrVj7zmc8wYcIEUlJSMAyD48ePX/X1ROSDKQyIxMjPf/5zjh8/zpe//GVeeOEFHnvsMc6ePcvkyZPZsGHDxfcdOHCAm2++mY6ODv74xz/y61//mkOHDjF9+nRqa2svvu/FF19k5cqVlJeX86c//Yknn3ySYcOGcdddd/HQQw9dcu2KigqefPJJfD4fCxcu7FL/V9oXwC9+8QtOnDjBrFmzKCws7NL17rjjDv7zP/+T733ve6xatYpJkybxiU98gt/97neXvG/9+vWsW7eOsrIypk6d2qVrichHMEUkJs6cOfO+15qbm82ioiJz9uzZF1+76667zIKCArOxsfHia8ePHze9Xq/5jW984+JrtbW1ZjQafV/NRYsWmWlpaWYwGLz4WiQSufi/t2/fbgLmb37zm6vq/0r7eu/1xowZY86cOfOqrrVy5UoTMH/3u99d8vrcuXPNkpISMxwOX/ZajzzyiAmYx44du6rriciH08iASIz06dPnfa9lZGQwevRoKisrAQiHwzz//POUl5eTlZV18X0DBgzglltu4Zlnnrn4WkFBAYZhvK/mDTfcQCAQoKGh4eJrLpe9b+Wr6SsW13vmmWfIyMjgrrvuuuT1T33qU9TU1PDaa6/F7Foi8tH0XSbSjRobG9m5cydjxowB4MiRI7S1tTF+/Pj3vXf8+PFUVFQQDAY/tObGjRspLCy8bPjoqlj0dTX27NnDqFGj8Hg877vWO78uIj1HYUCkG9199920trbyne98B4D6+noA8vLy3vfevLw8TNPk3LlzH1jvV7/6FS+99BL33Xcfbrc7Zn3a7asr1/uga727HxHpGZ6PfouIdMV3v/tdnnzySX72s58xYcKES37tcsP/H/Vrq1at4u677+bOO+/knnvu6VJP0WiUaDR6ybXeHSq60tcHMU2TSCRyyWvvHgmI5bVExB6NDIh0gwcffJAf/OAHPPzww/zd3/3dxdfz8/OByz/5NjQ0YBgGOTk57/u1NWvWcMcddzB37lyefPLJLt8sP/3pT+P1ei/+M3v2bFt9fZhNmzZdci2v13vxI4H5+fkfeC24/AiFiHQfjQyIxNiDDz7IAw88wAMPPMC3v/3tS35tyJAhpKamsnv37vf9vt27dzN06FD8fv8lr69Zs4Zly5Yxc+ZMnn76aXw+X5d7e+CBBy4JJ5mZmV3u66NMmDCB7du3X/JaSUkJAOPGjeP3v/894XD4ktGCd64/duzYq7qWiNjk8KcZRHqVhx56yATM++677wPf87GPfczs06eP2dTUdPG1EydOmD6fz/zmN795yXvXrFlj+v1+c86cOWZbW9sV9dDVjxZeTV/v1pWPFr7wwgsmYP7hD3+45PUFCxa876OF76aPFop0D40MiMTIT37yE+6//34WLFjAokWL2LZt2yW/PnnyZMAaOZg0aRKLFy/mW9/6FsFgkPvvv5+CggK+9rWvXXz/1q1bWbZsGcXFxXz7299m165dl9QbPXr0JR8DfGdnwnd2DNyxYwcZGRkA3HnnnR/Z/5X29U7td4b8m5qaME3z4vUnTZrEgAEDPvRat956K3PnzuWLX/wiTU1NDB06lN///vesXr2a//mf/7lkHUNtbS2bNm0COkcOVq1aRWFhIYWFhcycOfMj/2wi8hGcTiMivcXMmTNN4AP/ebcdO3aYs2fPNtPS0sysrCxz2bJlZkVFxSXv+d73vveh9TZu3HjJ+6/02h/mSvoyTdP867/+6w+81pWOSDQ3N5tf+tKXzOLiYtPn85njx483f//737/vfRs3bvzAa13tiISIXJ5hmqbZ7YlDRERE4pY+TSAiIpLkFAZERESSnMKAiIhIklMYEBERSXIKAyIiIklOYUBERCTJKQyIiIgkOe1AKCIiMTd73A22a6zf/XoMOpEroZEBERGRJKcwICIikuQUBkREJD7V18Mf/whvv+10J72e1gyIiEhccJkmowMtTGxpZFJLExQWgmnC9dfDG2843V6vpjAgIiKO6dPRzqSWRia1NHJdSzMZ0QhNbjdvZGQz8t9+Bq2tcM89cOwYDBrkdLu9lk4tFBGRmPugTxP4oxHGtzYzqaWRic1NlHUEiQD70jLYkZHN9owsDqemEzUM69MELS3WCMH3vw9f/3rP/iGSiEYGRESk+5gmg9rbmNTcyMSWJsYFmvGZJme8PrZnZPMfRaW8mZFFq/sDbkcZGbBgATz9tMJAN1IYEBGR2Kqv55bz9UxsaWJiSyMF4RBBw8Vb6Zn8v6J+bM/MptLnB8O4snrl5fCXfwlVVdCvX/f2nqQUBkRExJ5wGLZtgzVrrH927OA+0+RoSirrc/LZkZHF7rRMQq4ufoBt8WLweuGZZ6z1AxJzWjMgIiJX78SJzpv/+vXQ2Ah5eTB3Lsyfz8d/9Bh1Xp+tS1yyA+HChRAIwEsv2etbLksjAyIi8tECAdi0CVavtgLAwYPgdsPkyfC1r8H8+TBhgvUaUPfTn8f2+uXl8NnPwpkzUFQU29qiMCAiIpdhmrBnT+fT/5Yt0N4OZWXWjf/hh2H2bMjJ6Zl+li6Fz38enn3W+rfElKYJRETEUl8Pa9daN/8XX4SaGkhNhZtvtgLA/PkwYsQVLfzrloOK5swBl8vqTWJKIwMiIskqHIbXXut8+t++3RoRGDsWPvEJ6+Y/fTr4/U53aikvhy99CRoarPUJEjMKAyIiyeRyC/9yc2HePPjCF6x/l5Y63eXl3X473H03PPccfPKTTnfTq2iaQEQkDvzl0rtt18icPPR9r3k7Ohh27AijDx9g9KGDFNedJWoYHOs/gH3DRrB/+EhOlPbHfNfH/hreOGG7l+r2Rts1Wo7vf99r/3F0H81uN/cOGHFFNd7cvc12H8lAIwMiIr2JaVJy5hSjDx9k1OGDDD1+FG84TH1OLvuHjeC5ebdycMhw2lJTne60S9Zn5fLlM5WkR8IfvGuhXDX9TYqIJLiMjiBjaiu5/ukdjD58gJymJjq8Xg4PGsKz8xexf9hIThf2ufId/+LYhqw8/v70SaY3n2d1ToHT7fQaCgMiIgnGFY0y5NwZxtVWMv7sSQadP4sLqC4qZsf469k3bAQVAwcT9nqdbjXmTvtS2J2azpymBoWBGFIYEBFJAPmBZsbVVjLu7EnG1FaRHu6gxZvCnsJ+rB84hj2F/Qnfcq3TbfaI9Vl5fP5sNf5ohKDL7XQ7vYLCgIhIHPKFQ4ysr7kYAEpbzhPFoCK3iNVDrmF3nzKO5hRiGp0L/zId7Lcnrc/K494zlUxrbmR9tj5iGAsKAyIi8cA06dfcwLizlYyrPcmI+lP4ohHqUjPYXdifp0feyN6CUgK+OPnMv4OqUvwc9Kcxp6lBYSBGFAZERJxSXw/r1sGaNTy29o/kBVvpcLnZX1DKH0dNZnef/tRk5PaKhX+xtj4rj7+uq8EXjdLR1dMQ5SKFARGRnhIOw+uvd2768/rr1o5/Y8awrWQou/uUcTC/LyF9ZO4jrcvK42/PVjG5pZHNWblOt5Pw9BUnItKdTp68dMe/8+etHf/mzoXPfc7a8a9fP34fg02HkskxfypHU/zMbmpQGIgBhQERkVgKBGDz5s6jfg8csA7XufFG+MpXrP3+J068eNSvdN36rDw+Xn8GTzRKWFMFtigMiIjYYZqwd2/n0//mzdZRv/37Wzf+73/fOuo3V0+vsbYuK4/P1tYwqbWJVzNznG4noSkMiIhcrYYGa+Hf6tXWcbrV1dbJfjffDP/0T1YIGDlSC/+62SF/GpXeFGY3NSgM2KQwICLyUd678G/7dohGYcwY+PjHO4/6TdD9/hOWYbAuO4+l52r5xxKTqMJXlykMiIhcTmWldeNfvfpDF/6Js9Zn5fGpulNc19rMGxlZTreTsBQGREQA2tpg06bOp//9+zsX/t17r/X0P2mSFv7Fmb2p6Zzy+pjT1KAwYIPCgIgkJ9OEffs6V/2/s/CvXz9YsAAeekgL/xKBYbAhK4+5jfX8c98BmJoq6BKFARFJHu8s/Hvn6f+dhX8zZ8IPf2g9/Y8apYV/CWZ9Vi5/UX+a8W0tvJWWLCc0xJZhmqbpdBMiklxGXHON7Rod59s+8j1u02R8R5CZbQFmBAOM7wjiBg55fWz2p7HZn8bW1jBBmzf/gskzbf1+gLJ0+0Pc7eGo7RoA9Ufesl2jNZRiu0ZmuPWK3ucyTV6o3sfq9FwezS255Nd2nNhlu49koJEBEelV+oZDzAgGmNEWYFowQLYZ5bzLxcv+NP6Q0Yct/jROebwX398eaHKwW4mFqGGwMS2bWYHzPJrTVyM7XaAwICIJLSUa5cb2tosBYFi4gwiwy+fn11k5bPan8bbPr4+d9XIb0rK5q6WekaE2DvjSnG4n4SgMiEhiMU2GtrczpekcM4IBbmhvw2+a1Lg9bPGn8WhOHlv9aTS5tOo/mexMyeC8y82sQKPCQBcoDIhI3MsOh5na0sJNzc1Ma2mhbyhE0DB4LSWVR7Lz2ZyaRoXHp+HhJBYxDDalZjM7cJ4nsov1tXCVFAZEJO64TZNxgQA3NTdzU0sL4wMB3MDhlBRWZ2ezNTOTl0Mu2nU4jbzLhrRslrY2MCQU5IhPu0FeDYUBEYkLxR0d3NTczPSWFqa0tJAdidDodvNKRgZP9evH1owMTvt8F99/JZ8mkOTyuj+DZsPF7ECjwsBVUhgQEUekRKNMam21nv6bmxnW3k4EeDstjf/Kz2drZia709KIaLhXrlDIcLElNZtZbY38e06x0+0kFIUBEekZ7+z4t2YNvzp6lEmtrfhNk1NeL1szMni8qIhXMzJo9OjHknTd+rRsFtadY0AoyAmv3+l2Eoa+60Sk+5w7d+mOf1VV1o5/Hg//UlzMlsxMjqSkaLGXxMw2fyYBw8WsQCO/yVYYuFIKAyISO5HIpUf9vv66ddTv6NFw113Wdr8zZvCZyZOd7lR6qXaXi5dTMy+EgSKn20kYCgMiYk9VVefNf+1a66jfnByYMwc+8xnrqN/+/Z3uUpLIhtRsfhg4SUm43elWEobCgIhcnbY264S/dwLAvn3WUb833ABf/nLnUb+a+xeHvJyaRdAwuCXQ6HQrCUPfrSLy4UwT9u/vvPlv2gTBIJSWWjf+Bx6wjvrNy3O6UxEAAi432/zWVIFcGYUBEXm/yy38S0mxjvp9+GErBIwerYV/Erc2pGXzUH2ldUx1aanT7cQ9hQERsRb+bd9u3fhXr+5c+DdqFNx5JyxYADNmQKo2cpHEsDk1mxBVeP/8Z7jnHqfbiXsKAyLJ6t0L/9ats0YD3ln49zd/Yz39a+GfJKgWl5vX/RlMe/pphYErYJimaTrdhEhvNr5sjO0au6sO2a7hN00W5qYwqz3M7I4wI8NRIsBOr5sNPjcbUjy86XV/5I5/owr72u6llhTbNepPxGY+uCQGB9xNnnyD7Rr154K2a6zZst52DYCB/UfYrtGREoMppPoztn77ncFW/rGtCU6dgj597PfTi2lkQKS3Mk1GAfMxmY/JTExSzwWocRlsSPHwz+kpbErxcN6leX/pndb7/BBshmefhc99zul24prCgEgvkmOazL5w85+PSRkQBDZh8B1cbMtP5aDHpYV/khTOudzWotenn1YY+AgKAyIJzGWaTOKdp/8oNwJuYB/wNAZrMNiMQduFm3+e1+1gtyIOKC+39r9oaNDHXz+EDgMXSTClpsmnzCj/a0aoJcI2InyFKNUYfB4XZbgZY3j4quFmjeG6GAREktLtt1uflnnuOac7iWsaGRCJcz4zylwzenHofywQAbYDP8NgDS5eBx31K3I5ffvC1KnWVMEnP+l0N3FLYUAk3pgmg8IdTAu2MC3YwoT2AH5MqoA1GDyEwToMzunmL3Jl7rwTvvlNaGqCrCynu4lLCgMicSAzGmFysJWpwRamBlvoGwnTjsEbKWk8ntWH3zbWsRe08E+kK+64A77yFVi5Ej7xCae7iUsKAyIOcJkmYzramHYhAIzraMMNHPH4WJeaxSv+DN5ISSPospb17G2qd7ZhkURWVmYdnvXUUwoDH0BhQKSH9AmHLjz5tzK5vZWcaIQmw8Vr/nR+kNuXl/0ZnPZ4nW5TpHcqL4cHH4TWVkhPd7qbuKMwINJdgkHYvJmvnTvNtGArQ8PtRIE9vlT+kJHLy/4M9vhStfBPpCeUl8O3vmWdvVFe7nQ3cUdhQCRWTBMOHLj0qN+2Nha4Pbziz+CX2QVsS0mn0a1vO5EeN3QoXHON9akChYH30U8lETvOn4f1662njTVroLLSOup3xgz4/vdh/nzm3voxLfwTiQfl5fDII9Debn2fykUKAyJXIxKBHTs6n/5fe816beRIa8XyO0f9pr3r9BsFAZH4UF4O998Pa9fC4sVOdxNXFAZEPkp19aVH/TY0QHa2ddTvz39uHfVbVuZ0lyLyUUaPtoL7U08pDLyHwoDIewWDsGVLZwDYs8d6up80Ce6+27r533gjePTtI5JwysvhiScgFAKvPr3zDv00E/mAhX+UlFg3/vvus0YB8vOd7lRE7Covh4cfho0bYd48p7uJGwoDEnf86aNjUifDV/mBv5YVNZkRjjArHOaWUIT+pkkQ2OZxs9HjZn1mKgdaG7nryF54Yi888dMu95GZY/88ME9l2HYNAH/Yfh3DW2i7Rrj+uO0ag1yx+fH16W9/2XaNRZOus13jd7//X9s1/rzdtF0DIHAuYLtG7vBc2zVCIfunbBaklVz6gmmy3XCzeeHtfM2XfcV16gI1tnuJZwoDkhRcpsm1kSizwhFuCYWZEIniAQ65DFZ6PWz0unnF49YJfyK9nWHwvNvP/wm38fdmFlF9zwMKA9KLFUej3BKKMCscYWY4TJ4JjcBmr5u/96Ww0eumyqVTvEWSzQq3n3vCrUyJdvCyWx8xBIUB6UVSzCg3RQPMibQyN9LKWNPa8e9Nt4tf+7xs8Hp4w+3Sjn8iSe5Nw0u14WJxJKgwcIHCgCQu02S42cHcSCtzIy3MiAZIw6Ta8LDelc6jKbDJ4+GcSzd/EelkGgbPu/zcFgnybTMLUw8ICgOSWLLNCLdcePKfG22hzLSO+t3qSuMhbyFr3ensNVLAMD50AaGIJLcV7lQ+HwkwIRpih9vndDuOUxiQuOYyTa6PBpkbbWFepJVJ0TY8wEHDx3J3JmvdGWxxpdFmaO5fRK7c6y4vZ3CxJBpUGEBhQOJQ32iIudFW5kRamR1pJZ8I53Gx0Z3Ol33FrHVlcNKlzUJEpOuihsELbj9LIkG+58lM+m3DFQbEecEgbN1qbfizejXHghVEgR0uP7/05LLWnc7rLh31KyKxtcLt51ORANeYYd4ykvsBQ2FAep5pwsGDnTv+vfSSteNf374wfz5/eaiBDe506g19eYpI93nF5aMBg8WRNt5K8tFG/bSVntHYeOlRvydPgs9nnfD30EPWtr9jx4Jh8Kc/xmYHQhGRDxM2DFZdmCp4OMmnChQGpHtEIvDGG51P/9u2Wa+NGAHLllk3/5kzIT3d6U5FJImtcPv5i0gbI80wB5J4qkBhQGKnpgZefNF6+l+3DurrISvLOuTniSesADBggNNdiohctNmVQhMGt0WCHEjiqQKFAem6dy/8W7MGdu+2htkmToQvfrHzqF8dEyoicarDMFjjTmFxJMg/ezOdbscxCgNy5UwTDh26uOr/vQv/+Pa3rVGAggKnOxURuWIr3KncFTnHkGiYIzE6DTPRJOefWq7cOwv/3nn6P3HCWvg3fTo8+KAVAsaNS+qFNyKS2Da6UmjFYHEkyGOuDKfbcYTCgFwqEoGdOztX/b974d/SpVr4JyK9TpthsNadwpJIkMe8CgPigFx3bG6qvtS8Lv/eomiEmyNBpoebmRWNkI9JI7DJ8LDO5WG918/JYzXwi99Y/3wIlz+1y328oyAjNtMMgyZMt10jv1+J7Rozq+ps16iqPGK7BkB7qv0T2ip27rZdwxxo///jKeOLbdcAyIrar7Hr+FnbNfaciMFZGmZstuU+02L/6604e5btGpmZ9n8+njn91hW97wWXl3+PBBkYDVLpctu+bqJRGEhCPtNkcrSdmyNBbgm3M9oMEQV2Gi5+5fKy3uVmu+EmrKF/EUkS69xegiFYFGnnF640p9vpcQoDycA0GWKGuSUS5OZIO9Mi7aRhctpw8ZLbz7+6M9nkTqEq3OJ0pyIijmg1XLzk8rE42sEvUBiQXiLTjDI90n4hAAQpMyO0A6+5UnjEm8VLnhT2GV4t/BMRueB5t4/HQy0UmxFOG8k1VaAw0EsYpsk10dDFm//EaAceoMLw8KI7lQ3uFF51pxDQUb8iIpf1ottHKAQLIx382mN//VMiURhIYEWmySwizDYj3Bw4RT5RmjDY4k7hW74cXnL7qUzSz8yKiFytRsPFZpeXxQoDEs98pslkosw2I8wiwlhMosAuXPynN52Nbj87XT4t/BMR6aKV7hQeCbVQYEapS6KRVIWBeGaaDMFk9oWn/2lESQdOAxtw86jhZiNuGgwDny/b6W5FRBLearePR0Jwa6SD//b4nW6nxygMxJks02Q6EWabUWYRYQAm7cCruPiR4WUDbvZiaOGfiEg3qDdcvOrysijSrjAgPccwTa4lyqwLw/+TiOIBDmOwGjcbDDdbcRHQzV9EpEc87/bx/VAr2WaUxiSZKlAYcMKpU9ZRv2vWcMhsowBoAjbh5u8NL+txU5kkX4AiIvHmBbePfwy1Mj/SwR+TZHRAYaAntLdfetTv229bw/wTJvBbPGww3GzHpYV/IiJx4IzhZrvLw5JIu8KA2GCacPhw581/40YIBKC42Dro51vfgrlzoaCAh2N0NoGIiMTOSncK3wm1kmFGaUmCkVqFgVhparr0qN/jx62jfm+6Cb73PSsEjB+vhX8iIglgpcvHQ7QyNxLiGY/9A77incJAV0Wj1lG/79z8X3nFOup3+HBYssS6+d98s476FRFJQFUuN28aHhZF2hUG5D3etfCPtWuhrg4yM2H2bHj8cSsADBrkdJciIhIDK90+vhoOkGqaTrfS7RQGPkx7O7z8snXzX736koV/fP7z1s1/8mTwep3uVEREYux5dwr3hQPcEu1wupVupzDwbh+28G/ePPjmN62Ff4WFTncqIiLd7JjLzV7DzZJI7w8Dhmkm1vjHNcNvjEmdloj174xohCltTcwInGd6WyP9wx20Y/CGP5MtadlsTs3mgC/1sgv/iorzbfcx6vqxtmsALCy9wXaNXZETtmvkF9gPShOvGWK7BsDIiZNs16j+1Z9s19h+ao3tGv/4q2ds1wCobAvZrhE5Z/8Al5E3TrZdo19Rru0aAKETb9uu0Z4x0HaNjNSg7Rqt5wO2awBUnNtnu0b//FG2a/TtP9R2jV0vvWTr93850MTng81ktrVBSu9dO5B0IwOGaTKqvZWJLU1Mb2vkumALXkyOeVPYmJbD5tRsXkvNIuBKrrOsRUTk/Vb5Uvl6W5O1TmzxYqfb6TZJEQbywx1MCTQxNdDIlEAjeZEwzYaLV1OzebBgAFtSs6ny9t7EJyIiXXPI7aHC5WHo008rDCQabzTKdcEW6+bf2sjIDmvobG9KGk9nFfJKejaveDMIJ8FGEiIiYoNh8EJKKl9avhxCoV67YLx3hAHTpCzUztRAI9NazzOxrZk0M0qd28sraVn8NreYbWnZnPN0/p8YjjjYr4iIJIwXfKl86dxZa1H5vHlOt9MtEjYMpEci3NBmDf1PbW2kX7idEAZvpmbw73klvJyew6EPWPgnIiJypfa6vdYeMk8/rTDguGgU3nyTzzTUMLW1kfEXFv6d8KawJT2bV9Ky2Z6WRZsW/omISCwZBpSXw3/9FzzxBLh7330mvsPA6dOX7vhXW8unXC5eT83iR4UDeCU9m2ot/BMRke5WXg4//rF1Au3MmU53E3PxFQY6Oqwd/1avtgLAW29Zr0+YAJ/9LMyfz8y/+ZoW/omISM+64QYoLbWmChQGYsw0oaLi0h3/WluhqMial/nGN96345+CgIiI9DiXyxodeOopePRR6797kZ4PA01NsGFDZwA4dsz6qMZNN8F3v9t51G8v+4sWEZEEV14O//qv8NprMGWK093EVPeHgQsL/y456jcchqFDYeFCWLDAOuo3I6PbWxEREemyadOgTx9rqkBh4ApcZuEfmZkwa5aVqubPh8GDu+XSIiIi3cLthttvt8LAI4/0qo+uxyYMvLPw752n/127rNfftfCPKVN67c5NIiKSJMrL4Ze/tEa8r7/e6W5iputhoKKic9X/exf+ff3r1sK/Pn1i2KqIiIjDbr4Z8vKs0YGkDgMdHVBQAM3NWvgnIiLJxeuFpUutTxX84Ae9Zqrg6sPAz39uBYGnnrICgBb+iYhIMikvh9/8BvbuhbFjne4mJq44DDx+z98CkF9XzyeA5f/1WypfWn9VF/Nn2h81mDV/uu0aAGWl9hcwFvcv/Og3fYTb591luwaA334rjH7zgO0akaD9dSFjJw+xXSNW8j9j/+vt4KM7bdeY0He07RoA+en2w3thRpbtGuOusf//sXkqNkue/nvvS7ZrdBw5ZbvGXX/5F7ZrFMXoKfX8mj22a5yuP227hmGm267RanbDqXRz5kBWljVV0EvCwFXfnevz82jKzGTw0WPd0Y+IiEh8S0mBxYutMNBLXP2jumFwdPBABh07bu0gKCIikmzKy2H3bjh0yOlOYqJL4/bHBg0io7WVPmdrY92PiIhI/FuwANLSes3oQJfCQE1JX4IpKQw6pqkCERFJQmlp1i66yRwGTJeL4wMHMPjo8Ri3IyIikiDKy+GNN+D4cac7sa3Ly/uPDh5EfkMDWY2NsexHREQkMSxaZC0m/POfne7Eti6Hgcr+/Qm73RodEBGR5JSZae262wumCrocBkI+L1X9+mndgIiIJK8777RO462udroTW2ztAnR08ED6njqNv60tVv2IiIgkjiVLwOOBZ55xuhNbbIWBY4MGYpgmA4+fiFU/IiIiiSM3F2bPTvipAlthoC0tjdPFxdYGRCIiIsmovBw2b4baxN17x/ZhAUcHD6TsZCWeUCgW/YiIiCSWZcusfz/7rJNd2GI7DBwbNBBvOEy/qqpY9CMiIpJYCgthxoyEniqwHQbO5+bSkJujjxiKiEjyuvNOWL8ezp1zupMusX+mMNZZBQOPH8eIRmNRTkREJLHcfjuEw/Dcc0530iUxCgMDSWsLUnz6TCzKiYiIJJaSEpg6NWGnCjxX+sa/+9kTH/yL0SiUlFA+eCg88sgHvu13P/zuVTV3Of0L82zXADACTbZrvL3toO0as8dPtV0DwF9YaruGx+O1XeM/1/2P7RoPTLnXdg0AP9m2azSdsb/dtnHe/j4cM2+bYLsGwO39brBd43CN/c1VThyyv+q6tPiKf3x9qNyIYbvG/pYG2zUOnT1vu8agcdfZrgHg9mTYrmHU1NmusbfF/vdfWyQ2i9vT/B/9d3JPOMSDkQ7KUtJpMd7/dRUItsSkl+4Qk5EBXC647TZrJaVpxqSkiIhIIlnuduMHFkQjTrdy1WITBgCWLoWKCti/P2YlRUREEsVJw8VOw8Xt0bDTrVy12IWB2bMhPR2WL49ZSRERkUTyjMvNvGiE1AQbJY9dGPD7YcGChN50QURExI5nXR7SgXkJNlUQuzAA1lTB669DTU1My4qIiCSCIy4Xuw2DpQk2VRDbMLBoEbjdsGJFTMuKiIgkiuUuDwujEXwJNFUQ2zCQl2dtyaipAhERSVLPujxkAbMTaKogtmEArAMbNmyAJvuf4xcREUk0+wyDg4bBsqQOA0uXQkcHrF4d89IiIiJxzzBY7vKwKBrGkyBTBbEPAwMGwDXX6COGIiKStJ51uckDZpiJcWZP7MMAWFMFK1dCKDbbQIqIiCSSXYaL4xgsiyTGpwq6JwwsXQqNjbBpU7eUFxERiWuGwbNuN0uiYVwJMFXQPWHg2muhrExTBSIikrSedXkoAqYmwFRB94QBw7BGB5Yv18FFIiKSlLYbLqoxWJYAGxB1TxgAKwxUVsKbb3bbJUREROKVaRgsd7tZGolgxPmDcfeFgRkzICdHUwUiIpK0nnV5KMVkYpxPFXi6rbLXa21PvHw5PPggAA0NrbbLHqw8YbsGwNuH7NepaaywXSOlaIDtGgCfK/m87Rqrfvus7Ro7926zXePF51bZrgEwe8Z02zUOHT1ju0b+yKG2a+SeabddA6CyucF2jUiDYbvG88v/03aN/1P+Cds1AHJSc2JQxf55LM1v7rBdY1vrIds1AAINHbZrRDPtf5142u3fQL3E5iYcjXTtyX6rCWeApeH4nirovpEBsKYK3noLjh3r1suIiIjEo6hhsMJwc4cZies1dN0bBhYsAJ8PnnuuWy8jIiISr54x3AzEjOs1dN0bBjIzYfZsHVwkIiJJa5PhogHg6aedbuUDdW8YAGuqYMsWaLA/VykiIpJowobB84bbCgNxOlXQ/WHgttsgErG2JxYREUlCzxhuOHgQ9u1zupXL6v4w0Lcv3HijpgpERCRprTdc1tT5U0853cpldX8YAGuqYM0a3HH+0QoREZHu0GEYsHhx3K4b6JkwsGwZtLbS78TJHrmciIhI3LnzTti9Gw4fdrqT9+mZMDByJAwbxuCKIz1yORERkbizYAGkpcXl6EDPhIELBxcNPHI0bldSioiIdKu0NLj11iQOAwDLlpEeCFB86lSPXVJERCSulJfDjh1wIjZb68dKz4WByZMJpKYy6LCmCkREJEktWmTtzBtnowM9Fwbcbo4PGczgIwoDIiKSpLKyYP78JA4DwNFhQ8htOEdOvXYjFBGRJFVeDq+8AjX2T7uMlR4NA5VlAwh5PPpUgYiIJK/bbgOPB555xulOLurRMBDxejg5cACDNFUgIiLJKjcXZs2Kq6mCHg0DAMeGDqG45hRpra09fWkREZH4UF4OmzZBba3TnQDg6cmLHa47RU1OOrMMg7S33uSt4UOuukZV8/mY9HLmdLXtGq3NzbZr/OGX/267BsDLa1fbrlFf1Wi7RlPaOds1Nr8Wm3k0j3nIdo39B/bbruGP2g++BZGA7RoA1UH7W4IX+Apt1ygaUmq7RsPJFts1AHIKU2zXGHzO/nNVxDhmu8bJMzm2awAUFQyyXcMoTrddo7ba/vqyyiNnbNcACGX6YlLnomXL4ItftM7t+exnY1u7C3p8ZKDVn8KRPgWMO2n/ZiwiIpKQ+vSBGTPiZqqgx8MAwJ6yfgyvOY0vFHLi8iIiIs4rL4f16+Gc/RFVuxwJA7vLSvFGo4yqPu3E5UVERJx3++0QDsOKFU534kwYqM/MoCY3m7GVmioQEZEkVVoKU6bAU0853UnPLiB8t939S7npQAWuaJSoy5FMIiIi4ozaWli5ElpaYOdOa4TA49gt2ZmRAYDdZf1I7+hg8Jn4+FiFiIhItzFN2LMHfvhDmDoViorg05+G9HT40Y/A7Xa0PcdiSGV+LufTUhl/spqKvkVOtSEiItItvKbJjHA7i0NBFoXbYNw46+Y/bx78x3/AwoVWKIgDzo1JGAa7y0oZe7KKP99wHRiGY62IiIjEQn40wq3hIItCQeaFgmRhctJws9Lr54vLn4Kbbwa/3+k238e5MIC1bmD6gQpKzp2nJi/XyVZERESunmkyKhpmUaiNxaEgkyMduIHX3V5+7M9kpdfP2y4vGAZfXLDA6W4/kKNhoKK4D21eL+NPVisMiIhIQvCYJtPfNfw/JBqhFYN1nhS+mJrLC14/Z1zOrgG4Wo6GgYjbzb5+fRl7sprV1451shUREZEPlBeNsODC8P/8UJBsTKouDP/f6/XzksdPMIGnux0NA2BNFXzy2KvktrRyLsP+XtYiIiK2mSbDibDQbGeh2cGUphBuYIfby7/4M1np8bPL7e01690cDwP7+/Ul7HIxtrKaLaOGO92OiIgkKY9pMpXQxQAwlAgBYCM+7k7N4QVvKqcSbPj/SjkeBoI+H4eL+zDupMKAiIj0rBwzyjyzg4W0M9fsIBeTGlysMnx808jgJXy0GQahlBifWhhnHA8DYJ1VUP7aTlLbO2jr5X/hIiLirKFmmIVmBwvNdqYSwgPswsPPjVRWGinswoPZS4b/r1RchIE9/Uv52LY3GF1VwxtDBjrdjoiI9CJu02TKu4b/hxMhCLyEj68amawyfFQbvXP4/0oZpmmaTjcBwKRJMGgQ/PGPH/q2v7llXkwu98yeN23XCHVEbNfIy82zXQMgw59iu0aoLWC7RnPI/pdTYX6B7RoAadmptmucbbZ/tOgto0bYrjHzumtt1wBobra//fexJvtHj+/Zc9h2jdx0+99/AMXp9m8CWe32d3ZvSim1XaMlNTZ/JwfPttuu0Xb4bds1th/a3/XffP48rF5tnQi4apV1THBRESxeDEuWwJw51m6AAsTJyAAAy5bBP/0TtLdDiv0bm4iIJJmKCuvmv2IFbNliHf5zzTVw991WAJg4EXQw3mXFTxhYuhTuuw82bIBbb3W6GxERiXfhMLz6amcAOHAAfD6YNQsee8waBSgrc7rLhBA/YWDMGBg8GJYvVxgQEZHLa2yENWusm/8LL0BDA/TpA4sWwT/+I8ydCxkZTneZcOInDBiGNVXwu9/BE09oKEdERCzHjlk3/+eeg82bIRSyTgD8whes4f8bbtA9w6b4CQNgTRX89KewfTvceKPT3YiIiANcpsnYYBv8wz9YIWDvXvB64ZZbrHvE4sUwcKDTbfYq8RUGpk6F/HxrqkBhQEQkaaRFI9zY2sqM1hamBVrJjUTgV7+yhv8ffBDmzYPMTKfb7LXiKwx4PNaQz7PPWnM/IiLSaxWHQkxvbWF6awsT2gL4TJMjPh/Ls7LZkp7Bfxw/Cu7k/vx/T4mvMADWVMFvfwuHD8OwYU53IyIiMWKYJmOCwYsBYFhHOyFgZ2oa/1pQyNb0DKq979qFVkGgx8RfGJg7F/x+a6rg6193uhsREbEhNRrlxkAr01tbmNbaQn4kwnmXm5fT0/mPvHy2paXTqpu+4+IvDKSnW3NDzz6rMCAikoD6tLUy6fw5pre2MPHC8P9Rn4/nLwz/7/anEk2yvf/jXfyFAbCmCj7zGTh71vr8qIiIxC3DNBnZ2MDUszVMPXOKYc3nCQNvpqbxeH4hm9MzqPbpELp4Fp9hYMkSa9+BFSvgb/7G6W5EROQ9/OEwE+vPMPVMDVNqT5HfHqTJ62NbYTFPDhnJpqZ6WjT8nzDiMwwUFlofM1y+XGFARCROFLYFmHK2hmlnT3F9/RlSolFOpGfyYskAXikqYU9OPpELm/+0tZ53tlm5KvEZBsCaKvjud6G1VSdLiYg4wDBNhjeeY+rZGqadrWF403nChsHbuQX8v+HjeKWohKp0ffa/N4jvMPD3fw8vvgi33+50NyIiScEXDjP1TI01/3+2hoL2IM0eL9sK+/L7wSN5rbCYFq/m/3ub+A0Dw4bB6NHWVIHCgIhIt8kOBBhfWcU1lVWMPHWalEiEyrQM1peU8XKfEnbnFlwc/pfeKX7DAFijA//+79YxlR6r1aGji2JSuqDO/tBWR6v9PtJTztsvArja7E+lhCMR2zUCYfsfFzpy/A3bNQC8pv0a/kL7n2ZprLNfo7W9w3YNgIzcVNs1hjZU264RzA3brhHNtP9nARgazLFf5Dr7T8oP3v8b+31cKdOEN9/sPPr3jTesg35uugm+fC8sWUL/ESP4OPDxnutKHBTfYWDZMvjhD+Hll2HmTKe7ERFJXG1tsGGDdfN//nmoroasLOvI+K98xfp3Xp7TXYpD4jsMTJwIfftaUwUKAyIiV+f0aevGv2IFrFsHgQAMGQJ33WV9hHv6dOs0QEl68R0GXC647TYrDPzkJ9beAyIicnmmCW+91Tn8v3279XN0yhS4/34rAIwapZ+l8j7xHQbAmir45S9hzx4YN87pbkRE4kswCBs3dg7/V1ZaR/3Onw9/93ewcCEUFDjdpcS5+A8Dt9xifWEvX64wICICcOYMrFxpBYC1a639WAYOtB6eliyxplW1/a9chfgPAykpsGCBdXDRffc53Y2ISM8zTdi9u3P4//XXrdcnT4bvfMcKAGPGaPhfuiz+wwBYafcv/gKqqpzuRESkZ7S3w6ZNnQHgxAnIyLBOdf31r63hfx3kJjGSGGFg4UJrn4HnnnO6ExGR7lNbCy+8YN3816yBlhYoK7Oe/JcsgZtvtkZLRWIsMcJATo71TfDsszAsNpsOiYg4zjTJq21k4MFKmDYNXn3Vev2GG+Bb37ICwLhxGv6XbpcYYQCs3Qi/+lVSym6nPUULY0QkMbnCEUpOnGHgwSoGHqoi+1wLIa8HFi2GX/0KFi2CIj30SM9KnDBw221wzz0MPlnD/mEDne5GROSKpQSCDDhczaCDVfSvqCGlPURzVhrHR/Tj+Ij+VA8s5gvf/x+n25QkljhhoKwMrr+e4ccqFQZEJL6ZJrl1jdbT/8EqiitrcZkmZ0rz2TV1DMdH9KOuOFfD/xI3EicMACxdyuCHf4ArEiHqdjvdjYjIRa5IlL4Xhv8HHaoiu6GZkNdN1eC+vLRkMieGlxLITHO6TZHLSrgw4P/e9xhQfYZjZSVOdyMiSS6lvZ0Bp04xcO9pyiqqSQmGaMlM5fiI/my5dRJVg4qJeBPrx6wkp8T6Kh0/nvOZ6Qw/VqUwICKOyGlqYmBNDQNrqulbV4fLNDlbks9bk0dzbEQ/6vrmafhfEk5ihQHD4PCg/ow4cpI1MybpG05Eup0rGqW4rpZB1TUMrKkhp6WZsNtNZVERmyZM4ETfElon5jjdpogtiRUGgEOD+jHp7QMU1zZwuk++0+2ISC+U0tFB2alTDKyppuzUKfyhEK1+P8dLSnj52mupKioi7Em4H58iHyjhvppTJgyn48Ut3FhXy97rhnS5zrW5ubZ7eStaZ7tGaozOEg+3BW3XCKSatmv4YvDHyczItF8EiLTFYOQoZL+Xg81ttmukb33Vdg2AknT729dGveds1+hTlm27xv3/stx2jUscPty59e+WLRCJwHXXwSc/CUuWkH799YxxuRgT26uKxIWECwOm20XNiIGU7j/K3rmTnW5HRBJVOAyvvNIZAA4etLb6nT0bHn8cFi+Gfv2c7lKkRyRcGACoGTWIgbsOkl7fSGu+/ScMEUkSjY2werV181+1ChoarN3+Fi+GH/0I5syB9HSnuxTpcQkZBk4PKyPicVOy/yiHb7rO6XZEJJ4dOdL59L95szUicM018Ld/a+39P3EiuFxOdyniqIQMA5EUH2eG9KN0/zGFARG5VCRiHfjzTgDYvx98Ppg1Cx57zBoFKCtzukuRuJKQYQCgZtRgJix/CV9rGx3pqU63IyJOamqyjvxdscI6Ari+HgoLrUN/fvADmDsXMmOzMFWkN0rgMDCICcs30vfgcU5cP8rpdkSkpx071vn0v2kThEIwdix87nPW8P8NN4C2LRe5IgkbBtoz0qjvX0zpvmMKAyJJwIialJ4+x/Cjpxl+7Aw8Ohi8Xrj5ZvjJT6wAMHCg022KJKSEDQNgTRWM3vA67o4QkVh8wF1E4oqvI8zgE2cZfvQMw46fIb2tg9ZUHxUD+9DnF7+GefMgK8vpNkUSXkKHgepRgxi/5hX6HKni1KhBTrcjIjGQ3RRg+NEzDD96mgHV9XgiUc7mZbJrTBmHBhdRVZyH6TK45s47nW5VpNdI6DDQUphLU0EOpfuPKgyIJCrzneH/Mww/dpqiumYiLoMTpfmsu2k0hwYXcT5bn/0X6U4JHQYAakYPZuAb+yEa1WeFRRKEtyPM4JO1DD96mmHHz5AR6CDg91IxsIgtNwznyIA+tKdo6k+kpyR8GKgeNYiRm3eSf/I09QN1rLFIvMoOBBlVU8/oU3UMqT2PJxKlNjeDt0f159DgYir75mIq0Is4IuHDQEO/YoIZaZTuP6owIBJHDNOk9Fwzo2rqGF1TT0ljCxHD4FhBNuunjeLQ4CLO5WQ43aaI0AvCAC6DmpEDKdl/jLcXTAMjBifViUiXeMMRhp49x+iaOkaeqicr2EHA6+Fg33xeGlnGweI8gj4v6aXaAEgkniR+GACqRw1m8I59ZNaeo7lPntPtiCSVrLZ2RtXUMepUPUPPnMMbjVKbkcqusiL2leRzIj+bqIb/ReJarwgDZ4f0I+zzUrrvKAcUBkS6l2lSer7FCgA1dfQ7bw3/Hy/IZs3YQewvKaAuM83pLkXkKhimaZpONxETd94JlZXw2mtX9PbPLr3J9iX3Hj1tu0ZJn9icq9DQGrFdY8+Zs7ZrREJh2zX8wdhk1NrmZts1Qh0dMegkx3aF0gH59tsAqo5XdO03trXB+vXW1r/PPw81NZCdDbfeau38t2AB5CmIiySqXjEyAMDSpfBXf2X9kCrRQkIR206dsm78K1bAunVWIBg6FD7+cSsA3HSTtR2wiCS83hMGFi2yDiV57jn4whec7kYk8Zgm7NrVefjPjh3W3h3TpsGDD1oBYMQILdIV6YV6TxjIy4MZM2D5coUBkSsVDMKGDZ3D/1VV1l7/CxbAl79sTQPkx2aKQkTiV+8JA2BNFXzjG9bZ5jq8ROTyTp+GlSutALB2LQQCMHgwlJdbT//Tp4PP53SXItKDel8YuPdeWL0aPvYxp7sRiQ+myahQB3MDAbjxRnj9dWv4f8oUuP9+KwCMGqXhf5Ek1rvCwMCBcM011lSBwoAkMZ9pMjXYxpxAK3PaAvSLhGk2DCgrg7vvhoULoaDA6TZFJE70rjAA1ujAY49BKKSVzpJU8iNhZrcFmBsIMCMYIN00qXR7eDEtnbWpabzmT+Xon/7kdJsiEod6XxhYtgweegg2bYI5c5zuRqT7mCYjQx3MaQswN9DKdR3tALzpS+Fn2bmsS03jgNen4X8R+Ui9Lwxce601FLp8ucKA9Do+02RysI25ba3MDgQoi4RpNQw2+9P4en4W61PTqHf3vm9rEeleve+nhmHAbbdZYeBf/1VPRZLwCohyZ0sTc9sCzGwLkGGaVLs9rE1NY11aOq/6/bQb2vtfRLqu94UBsKYKHn8c3nwTrr/e6W5ErpLJaKIsIcQSQkwhgqveGv5/IiuXtWlp7Nfwv4jEUO8MAzNmQE6ONTqgMCAJwIvJDMIXAkCYwURpBdbi4bOksrNfH2o1/C8i3aR3/nTxeq2PTj37rLWNqkgcyiPKwgsBYD4hsoFKDJ7Hywq8bMBDO9bTf6mCgIh0o977E2bZMvjd7+DYMRg0yOluRACTke8a/p9KBDewHTc/wc8KPOzCDWj4X0R6Vu8NAwsWWFuqPvectce6iAM8wPQLQ/9LCDGUKAFgHR6+QCor8XIKLf4TEWf13jCQmQmzZ1tTBQoD0oNygVuBJcACIIdWqi8M/997Yfi/TU//IhJHem8YAGs3wrvvhvr69528NrJsqO3y9S3nbdfI7t/Hdg0A1+F22zWiNVW2a0RS/LZrmL402zUAOtrrYlLnihw82Hn078svQyRiLV5dsgSWLKH0+uv5vGHw+Z7rSETkivXuMLBkiXWc8cqV8Fd/5XQ30puEw7B1a2cAOHwY/H5rNOrf/g0WL4bSUqe7FBG5Ir07DJSUWKe0LV+uMCD2nT8Pq1ZZN/9Vq6z/Li62bvw//rEVBNLTne5SROSq9e4wANZUwcMPQ1sbpKY63Y0kmsOHO5/+t2yxhv+vvRbuuccaeZowwToOWEQkgSVHGPj2t2H9eusJTuTDhMPw6qudAeDAAUhJgVmz4Gc/s76G+vd3uksRkZjq/WFg1CgYNsyaKlAYkMtpbIQ1a6yb/wsvQEMD9Oljfb388IfWgVcZGU53KSLSbXp/GDAMa3Tgv/4LfvELcLud7kjiwdGjnU//mzZZIwLjx8MXv2gN/0+apOF/EUkavT8MgBUGfvxjeO01mDrV6W7ECZEIbNvWGQD27bM2pbrlFnj0UWsUYMAAp7sUEXFEcoSBKVOgsNCaKlAYSB5NTfDii53D/3V11tfBokXw/e/D3LnW5lQiIkkuOcKA220N/S5fDj/6kdPdSDfqFwkzryPA3FAbFBRAKARjx8JnP2t9Ddxwg6aKRETeIznCAFhTBb/+tbU6fORIp7uRGHGZJteFO5gbCjC3o41RkRAdwDavH37yE2v4XwdViYh8qOQJA3PnQlqaNTqgMJDQ0swoMzuCzAsFmN3RRoEZpcFwsd6byr+kZvOSN5UWl4uae+5xulURkYSQPGEgNRXmzbMOLvrmN53uRq5SaSTM3FAbczsCTA0FSQEOub38b0oGa32p7PCkEDV0+I+ISFckTxgAWLYMPvUpOH3a6U7kIximybXvGv4fEwkRArZ5/DyclstaXyon3F6n2xQR6RWSKwwsWmTtO7BihdOdyGWkmSa3hDtY3GEN//cxo5wzXGzwpvKz1Gw2elNp1mf/RURiLrnCQEEB3HSTNVUwuNDpbgQoiUZYEO5gQbiDGZEO/ECFy8PTKRm8eGH4P6LhfxGRbpVcYQCsqYJ/+Ae8/coJ+TTM3NMM0+TaaJhbLwSA8dEwYeBVt5fvp6SzyuPjtD/L6TZFRJKKYZqm6XQTPeroURgyhOeW3Mbh4SNsldq8d4ftdjJjNO/dpyNqu8ZPX1wbg04uIxCAdeus6ZmVK+HUKcjJgVtvtT77v2AB5OZ2z7VFROQjJd/IwODBMHYsQysqbIcB+RDV1fD881YAWL8egkHrwKhPfMIKANOmgVcjMyIi8SD5wgDAsmUM+vGPcUUiRLUbXWyYJuzc2bn3/86d1k5/06ZZW/8uWQIjFL5EROJRcoaBpUtJ/cEPKK2uprKszOluEldbm/XUv2KFNQpQUwPZ2daw/1e/ak0D5OU53aWIiHyE5AwDEybQnJHB0IrDCgNX69SpzuH/deusQDBkCHzsY9bT//TpGv4XEUkwyRkGDIMjQ4Yy5EgFG2+ZZe09IJdnmrBrV+fw/44d4HJZpz8+8IAVAEaO1N+hiEgCS84wAFQMHcq1b+2isPYstX2KnG4nvgSDsHEjPPecNQpQVWUd9btgAXzpS9bwf0GB012KiEiMJG0YqOxfRrvPx9CKCoUBIKO93TrVccUKWLsWWlut0/7uuMN6+p8xA3w+p9sUEZFukLRhIOp2c2zQYIYcqeDVqdOcbqfnmSZ9W1oYU1vHmLpa+jc2wZatMGUK3HefFQBGj9bwv4hIEkjaMADWVMHilc+T2dRIc1a20+10O3c0ytCGc4ypq2V0bR15wSBBt5uD+fm8PKYf/3fjS1CobZpFRJJNUoeBY4MGE3G5GFpxhDevv97pdrpFRkcHo+rqGFNby/D6BvyRCA1+P3sLC9hbWMiR3FwiFw7/+b8KAiIiSSmpw0BHSgqV/csYcuRw7wkDpklxaytjaq2n/wGNjQCczM5m/aCB7Cso5FRGuob/RUTkoqQOA2BNFczasJ6UYJB2v9/pdrrEHY0yvL6eMbV1jK6tJT8YpN3t5mB+Hv87ejT7C/JpSUlxuk0REYlTSR8GjgwZypz16xh07CgHRo12up0rltHRwbVn67juTC3j6+pIDUc4509hX0EhewsLqMjNJaytlkVE5AokfRhoyczkdFExQysq4jsMmCYlLa1cf7aW687UMuzceVzAkewsVg4eyPG8AmoyMjT8LyIiVy3pwwBYUwU3vP4a7nCYiCd+/krc0SgjG85x3ZlarjtbS1GgjXaXiz2F+fx63Gh29SnkvN8a/o/FEcYiIpKc4ufO56CKoUO56eWt9K88yfFBgx3tJaOjg2vO1nHd2VrG19aTFg7T4E/hzT6F/HefQvYW5BHS8L+IiMSQwgBQn1/A+exshlYc7vkwYJr0bQ1w/ZmzXHe2luEN1vD/0ewsVg0awJtFhRzPytTwv4iIdJukDQN11ccv+e89xX0Yf+gQ9aOGY17hjXf5qg1du3goBFu3dh7+U1EBfj/MmWPt/Ld4MYNLShgMlHftCiIiIlcsacPAe+0tLeWmg4fp19BAZX5+7C9w7hysWmXd/FevhvPnoW9fWLwYfvpTmD0b0tJif10REZGPoDBwwcmCfFp9PkZX1cQuDBw61Pn0v3UrRCJw3XXWyX9LlsD111vHAYuIiDhIYeCCqMvFgZK+jKquZs0147pWJByGl1/uDACHDkFKivXU//jj1ihAv36xbVxERMQmhYF32devlAnHT1DQ1ExdVuaV/abz561h/xUrrGmAc+egqMi68f/zP1vrANLTu7VvEREROxQG3uVwcREht4tR1TVsyRrxge/Lb25hZE0NzJoFW7ZYIwLXXAN3320N/0+cqOF/ERFJGAoD7xLyeDhcVMTo6mq2jOoMA65olLK6ekbWnGJUTQ19mpoJu1wwbx489pg1ClBW5mDnIiIiXacw8B77+pVyx+s7yG9uoeTcOUZV1zDi1GnSOjpoTknhYElf1owfR0VREQ8+9Wen2xUREbFNYeA9DpT0BeDrK1cBcCo7m9eGDmZ/SQlV+XlXvAeBiIhIolAYeI9Wv5/nJlyPyzTZX9KX8xla/CciIr2bwsBlvDZsiNMtiIiI9BgteRcREUlyCgMiIiJJTmFAREQkySkMiIiIJDmFARERkSSnMCAiIpLkDNM0TaebEBEREedoZEBERCTJKQyIiIgkOYUBERGRJKcwICIikuQUBkRERJKcwoCIiEiSUxgQERFJcgoDIiIiSU5hQEREJMn9f+bifhTJGkX1AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# it is also possible to filter for more meta data from the tiles,\n", + "# we match them with the meta data from the stac items: https://documentation.dataspace.copernicus.eu/APIs/STAC.html\n", + "# we support four types of filters: eq, ueq, lt, gt, in (equal, unequal, less than, greater than, python \"in\" with list)\n", + "# for example we can use this to filter for cloud cover:\n", + "arguments[\"filter\"].update({\"cloudCover\": {\"lt\": 90}})\n", + "ds = tg.create(**arguments)\n", + "ds = ds.where(ds != 65535) # remove no data values\n", + "visualize_sat_images(ds, gdf, list(reversed(arguments[\"bands\"])))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### filtering items manually\n", + "some properties are hard to filter with the available filter options, hence you can use your own logic to filter them\n", + "\n", + "here an example how you can filter them: some Sentinel-2 tiles are available in different processing baseline versions" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "import re\n", + "import pandas as pd\n", + "\n", + "\n", + "# function to filter for the highest Sentinel-2 processing baseline\n", + "def filter_doubled_processing_baseline_by_name(items):\n", + " \"\"\"use the name of the item to filter the processing highest baseline, except for the 99.99 fill value\n", + " see https://sentinels.copernicus.eu/en/web/sentinel/user-guides/sentinel-2-msi/naming-convention\n", + " \"\"\"\n", + " filtered_items = {}\n", + "\n", + " # re like \"N0301\" or \"N0500\"\n", + " pattern = re.compile(r\"(.*_N)(\\d{4})(_.*)\")\n", + "\n", + " for item in items:\n", + " match = pattern.match(item[\"id\"])\n", + " if match:\n", + " base_id = match.group(1) + match.group(3)\n", + " proc_date = pd.to_datetime(base_id.split(\"_\")[-1].split(\".\")[0])\n", + " base_id = \"_\".join(base_id.split(\"_\")[:-1])\n", + " num = int(match.group(2))\n", + " # If the base_id is not in the dictionary\n", + " # or if the current num is larger (not when it is 99.99, which seems to be a fill value)\n", + " # or if the num is the same but the processing date is newer\n", + " # update the entry\n", + " if (\n", + " base_id not in filtered_items\n", + " or (\n", + " filtered_items[base_id][\"num\"] < num < 9999\n", + " or filtered_items[base_id][\"num\"] == 9999 != num\n", + " )\n", + " or (\n", + " proc_date > filtered_items[base_id][\"proc_date\"]\n", + " and num == filtered_items[base_id][\"num\"]\n", + " )\n", + " ):\n", + " filtered_items[base_id] = {\n", + " \"item\": item,\n", + " \"num\": num,\n", + " \"proc_date\": proc_date,\n", + " }\n", + " else:\n", + " raise ValueError(f\"Could not extract processing baseline from {item['id']}.\")\n", + "\n", + " # Extract the filtered list\n", + " result = [entry[\"item\"] for entry in filtered_items.values()]\n", + "\n", + " return result" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 3kB\n",
+       "Dimensions:      (x: 27, y: 18, time: 1)\n",
+       "Coordinates:\n",
+       "  * x            (x) float64 216B 6.976e+05 6.976e+05 ... 6.979e+05 6.979e+05\n",
+       "  * y            (y) float64 144B 5.326e+06 5.326e+06 ... 5.326e+06 5.326e+06\n",
+       "  * time         (time) datetime64[ns] 8B 2021-01-01T10:23:29.024000\n",
+       "    spatial_ref  int64 8B 0\n",
+       "Data variables:\n",
+       "    B02          (time, y, x) uint16 972B 65535 65535 65535 ... 65535 65535\n",
+       "    B03          (time, y, x) uint16 972B 65535 65535 65535 ... 65535 65535\n",
+       "    B04          (time, y, x) uint16 972B 65535 65535 65535 ... 65535 65535\n",
+       "Attributes:\n",
+       "    crs:          EPSG:32632\n",
+       "    data_source:  CDSE\n",
+       "    collection:   SENTINEL-2
" + ], + "text/plain": [ + " Size: 3kB\n", + "Dimensions: (x: 27, y: 18, time: 1)\n", + "Coordinates:\n", + " * x (x) float64 216B 6.976e+05 6.976e+05 ... 6.979e+05 6.979e+05\n", + " * y (y) float64 144B 5.326e+06 5.326e+06 ... 5.326e+06 5.326e+06\n", + " * time (time) datetime64[ns] 8B 2021-01-01T10:23:29.024000\n", + " spatial_ref int64 8B 0\n", + "Data variables:\n", + " B02 (time, y, x) uint16 972B 65535 65535 65535 ... 65535 65535\n", + " B03 (time, y, x) uint16 972B 65535 65535 65535 ... 65535 65535\n", + " B04 (time, y, x) uint16 972B 65535 65535 65535 ... 65535 65535\n", + "Attributes:\n", + " crs: EPSG:32632\n", + " data_source: CDSE\n", + " collection: SENTINEL-2" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "items = tg.search(num_workers=4, **arguments)\n", + "# filter the items by the highest processing baseline, but discard fill value 99.99 if possible\n", + "items = filter_doubled_processing_baseline_by_name(items)\n", + "tg.download(items)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/mnt/datastorage/home/adrianhohl/Terragon/terragon/copernicus_data_space_ecosystem.py:417: UserWarning: No data found in bounds.\n", + " warnings.warn(\"No data found in bounds.\")\n", + "/mnt/datastorage/home/adrianhohl/Terragon/terragon/copernicus_data_space_ecosystem.py:206: UserWarning: Multiple files found for band VV: [PosixPath('Sentinel-1/SAR/SLC/2022/01/01/S1A_IW_SLC__1SDV_20220101T053510_20220101T053537_041263_04E777_7EED.SAFE/measurement/s1a-iw1-slc-vv-20220101t053512-20220101t053537-041263-04e777-004.tiff'), PosixPath('Sentinel-1/SAR/SLC/2022/01/01/S1A_IW_SLC__1SDV_20220101T053510_20220101T053537_041263_04E777_7EED.SAFE/measurement/s1a-iw2-slc-vv-20220101t053510-20220101t053535-041263-04e777-005.tiff'), PosixPath('Sentinel-1/SAR/SLC/2022/01/01/S1A_IW_SLC__1SDV_20220101T053510_20220101T053537_041263_04E777_7EED.SAFE/measurement/s1a-iw3-slc-vv-20220101t053511-20220101t053536-041263-04e777-006.tiff')].\n", + "Adding them as new bands.\n", + " warnings.warn(f\"Multiple files found for band {band}: {f_paths}.\\nAdding them as new bands.\")\n" + ] + }, + { + "data": { + "text/plain": [ + "Data variables:\n", + " VV_s1a-iw1-slc-vv-20220101t053512-20220101t053537-041263-04e777-004_1 (time, y, x) complex64 616B ...\n", + " VV_s1a-iw2-slc-vv-20220101t053510-20220101t053535-041263-04e777-005_1 (time, y, x) complex64 616B ..." + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# some collection like, Sentinel-1 - IW_SLC__1S, have multiple files (if not filtered)\n", + "# if multiple files are found, the files will be added as separate bands, with the filename as data_var\n", + "from shapely.geometry import Polygon\n", + "\n", + "gdf = gpd.GeoDataFrame(\n", + " geometry=[\n", + " Polygon(\n", + " [\n", + " (9.295, 47.425),\n", + " (9.296, 47.425),\n", + " (9.296, 47.426),\n", + " (9.295, 47.426),\n", + " ]\n", + " )\n", + " ],\n", + " crs=\"EPSG:4326\",\n", + ")\n", + "args = dict(\n", + " shp=gdf,\n", + " collection=\"SENTINEL-1\",\n", + " start_date=\"2022-01-01\",\n", + " end_date=\"2022-01-03\",\n", + " bands=[\"VV\"],\n", + " resolution=10,\n", + " filter={\"productType\": {\"eq\": \"IW_SLC__1S\"}},\n", + ")\n", + "ds = tg.create(**args)\n", + "ds.data_vars" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "sits", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/demo_files/demo_gee.ipynb b/docs/demo_files/demo_gee.ipynb index ad1fa04..7ec2c1d 100644 --- a/docs/demo_files/demo_gee.ipynb +++ b/docs/demo_files/demo_gee.ipynb @@ -703,7 +703,7 @@ "ds = tg.create(**arguments, num_workers=10)\n", "# or\n", "# items = tg.search(**arguments, num_workers=10)\n", - "# ds = tg.download(items, remove_tmp=False) # do not remove the temporary downloaded files\n", + "# ds = tg.download(items)\n", "ds" ] }, diff --git a/docs/demo_files/terragon_workflow.ipynb b/docs/demo_files/terragon_workflow.ipynb index e3df73e..f7674f9 100644 --- a/docs/demo_files/terragon_workflow.ipynb +++ b/docs/demo_files/terragon_workflow.ipynb @@ -92,7 +92,8 @@ "items = tg.search(**args)\n", "print(items)\n", "# here we will filter by taking the first item only\n", - "items = item[:1]" + "items = item[:1]\n", + "tg.download(items)" ] } ], diff --git a/docs/index.rst b/docs/index.rst index de85a63..4bcab25 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -15,6 +15,7 @@ source/base source/gee source/pc + source/cdse .. toctree:: :hidden: @@ -24,6 +25,7 @@ demo_files/terragon_workflow demo_files/demo_gee demo_files/demo_pc + demo_files/demo_cdse .. toctree:: :hidden: diff --git a/docs/source/cdse.rst b/docs/source/cdse.rst new file mode 100644 index 0000000..35a535c --- /dev/null +++ b/docs/source/cdse.rst @@ -0,0 +1,6 @@ +Copernicus Data Space Ecosystem (CDSE) Class +============================================ +.. automodule:: terragon.copernicus_data_space_ecosystem + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 54cd374..713b780 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -1,6 +1,12 @@ Changelog ========= +v0.1.0 +------ +- add support for Copernicus Data Space Ecosystem (CDSE) API +- fix for date extraction + merging of files for GEE API + v0.0.2 +------ - documentation - fix isometric pixels - various bug fixes diff --git a/docs/source/contributing.rst b/docs/source/contributing.rst index c301afb..bf9b6b0 100644 --- a/docs/source/contributing.rst +++ b/docs/source/contributing.rst @@ -12,7 +12,7 @@ Adding a new data source ------------------------ Each data provider has its own class in the package and inherits the structure from the Base class in 'terragon/base.py'. A new data provider will get a new class in the same directory. The workflow of a data provider should match the general workflow described in the demo_files. -Make sure to include tests (and later also documentation) as described in this document. +Make sure to include tests (and later also a demo and a documentation) as described in this document. Linting and Styling ------------------- @@ -40,7 +40,7 @@ The testing uses the unittest framework, each data source has its own file in th Code Documentation ------------------ -The documentation builds on sphinx and readthedocs. Each class/data source has its own .rst file in the folder 'docs/source/' +The documentation builds on sphinx and readthedocs. Each class/data source has its own .rst and demo file in the folders 'docs/source/' and 'docs/demo_files/'. They need to be listed in the 'index.rst'. .. code-block:: console diff --git a/pyproject.toml b/pyproject.toml index 5b3b774..5cdc03e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -36,6 +36,10 @@ pc = [ "odc-stac", "planetary-computer", ] +cdse = [ + "boto3", + "rasterio>=1.4.0", +] [project.urls] Homepage = "https://github.com/drnhhl/terragon" diff --git a/requirements.txt b/requirements.txt index 35a80cd..50cd65a 100644 --- a/requirements.txt +++ b/requirements.txt @@ -11,4 +11,6 @@ odc-stac pystac-client # gee earthengine-api -geedim \ No newline at end of file +geedim +# cdse +boto3 \ No newline at end of file diff --git a/terragon/__init__.py b/terragon/__init__.py index f2eb135..9730187 100644 --- a/terragon/__init__.py +++ b/terragon/__init__.py @@ -1,3 +1,3 @@ from .init import init # noqa: F401 -__version__ = "0.0.2" +__version__ = "0.1.0" diff --git a/terragon/copernicus_data_space_ecosystem.py b/terragon/copernicus_data_space_ecosystem.py new file mode 100644 index 0000000..460a500 --- /dev/null +++ b/terragon/copernicus_data_space_ecosystem.py @@ -0,0 +1,633 @@ +import itertools +import re +import warnings +from pathlib import Path +from urllib.parse import urljoin, urlparse + +import boto3 +import geopandas as gpd +import pandas as pd +import rasterio +import requests +import rioxarray as rxr +import xarray as xr +from joblib import Parallel, delayed +from rasterio.vrt import WarpedVRT +from shapely.geometry import box + +from .base import Base +from .utils import align_coords, meters_to_crs_unit + + +class CDSE(Base): + """Class to interact with the Copernicus Data Space Ecosystem. The images are downloaded from the AWS bucket. + The packages rasterio/rioxarray/boto will be used to download the images. + Currently only these collections are supported: COP-DEM, GLOBAL-MOSAICS, LANDSAT-5, LANDSAT-7, LANDSAT-8-ESA, + TERRAAQUA, S2GLC, SENTINEL-1, SENTINEL-1-RTC, SENTINEL-2. + + :param credentials: credentials to authenticate, expected format: {'aws_access_key_id': , 'aws_secret_access_key': }. If None, it will fallback to the credential handling from boto3/rasterio. + :param base_url: the URL for the STAC catalog, defaults to "https://catalogue.dataspace.copernicus.eu/stac/" + :param end_point_url: the URL for the data endpoint, defaults to "https://eodata.dataspace.copernicus.eu" + """ + + _file_extensions = [ + ".jp2", + ".tif", + ".tiff", + ".nc", + ".dt2", + ".dt1", + ".img", + ".JP2", + ".TIF", + ".TIFF", + ".NC", + ".DT2", + ".DT1", + ".IMG", + ] + + _supported_collections = [ + "COP-DEM", + "GLOBAL-MOSAICS", + "LANDSAT-5", + "LANDSAT-7", + "LANDSAT-8-ESA", + "TERRAAQUA", + "S2GLC", + "SENTINEL-1", + "SENTINEL-1-RTC", + "SENTINEL-2", + ] + + def __init__( + self, + credentials: dict, + base_url: str = "https://catalogue.dataspace.copernicus.eu/stac/", + end_point_url: str = "https://eodata.dataspace.copernicus.eu", + ): + """Initialize class and save the credentials. + + :param credentials: credentials to authenticate, expected format: {'aws_access_key_id': , 'aws_secret_access_key': }. If None, it will fallback to the credential handling from boto3/rasterio. + :param base_url: the URL for the STAC catalog, defaults to "https://catalogue.dataspace.copernicus.eu/stac/" + :param end_point_url: the URL for the data endpoint, defaults to "https://eodata.dataspace.copernicus.eu" + :raises ValueError: when the credentials are in the wrong format + """ + super().__init__() + self.base_url = base_url + self.end_point_url = end_point_url + if credentials: + if "aws_access_key_id" not in credentials or "aws_secret_access_key" not in credentials: + raise ValueError( + "aws_access_key_id or aws_secret_access_key not in credentials, could not initialize." + ) + else: + # fallback to credentials saved in the environment/files from aws + credentials = {"aws_access_key_id": None, "aws_secret_access_key": None} + self.credentials = credentials + + def retrieve_collections(self, filter_by_name: str = None): + """Search the collections provided by Copernicus Data Space Ecosystem. + + :param filter_by_name: name to filter the collections for, defaults to None + :raises RuntimeError: if the request to the collections endpoint fails + :return: a list of collection names + """ + collections_url = urljoin(self.base_url, "collections") + response = requests.get(collections_url) + + if response.status_code == 200: + data = response.json() + collections = [collection["id"] for collection in data["collections"]] + if filter_by_name: + collections = [ + collection for collection in collections if filter_by_name in collection.lower() + ] + warnings.warn( + f"Currently we only support the following collections: {self._supported_collections}" + ) + return collections + else: + raise RuntimeError("Failed to retrieve collections") + + def search( + self, + *args, + resampling=rasterio.enums.Resampling.nearest, + use_virtual_rasterio_file=True, + rm_tmp_files=True, + filter_asset_path={"COP-DEM": ".*/DEM/.*", "SENTINEL-2": ".*/IMG_DATA/.*"}, + **kwargs, + ): + """Search for items in the Copernicus Data Space Ecosystem collections via stac. For a description of the args/kwargs parameters see the Base class function. + + :param resampling: rasterio Resampling method is used to reproject the cubes, defaults to rasterio.enums.Resampling.nearest + :param use_virtual_rasterio_file: use rasterio virtual file function when True, when False whole file is downloaded, defaults to True + :param rm_tmp_files: only used with 'use_virtual_rasterio_file=False' to remove the files after the minicube is created, defaults to True + :param filter_asset_path: manual filtering of the filepath in the AWS bucket, used for collections with ambiguous file names, defaults to {"COP-DEM": ".*/DEM/.*", "SENTINEL-2": ".*/IMG_DATA/.*"} + :raises ValueError: when no items are found or parameters are in the wrong format + :raises RuntimeError: when the corresponding files for the items are not found + + :return: a list of items + """ + super().search(*args, **kwargs) + self._parameters.update( + { + "rm_tmp_files": rm_tmp_files, + "use_virtual_rasterio_file": use_virtual_rasterio_file, + "resampling": resampling, + "filter_asset_path": filter_asset_path, + } + ) + + if self._param("collection") not in self._supported_collections: + warnings.warn(f"Currently we only support collections: {self._supported_collections}") + if self._param("num_workers") > 4: + warnings.warn( + "More than 4 workers are not recommended, because only 4 concurrent connections are allowed: https://documentation.dataspace.copernicus.eu/Quotas.html." + ) + + shp_4326 = self._reproject_shp(self._param("shp")) + bbox = shp_4326.total_bounds + start_date = self._param("start_date") + end_date = self._param("end_date") + # make end_date inclusive + start_date = ( + f"{start_date}T00:00:00.000" if start_date and "T" not in start_date else start_date + ) + end_date = f"{end_date}T23:59:59.999" if end_date and "T" not in end_date else end_date + datetime = f"{start_date}/{end_date}" if start_date and end_date else None + + data = { + "bbox": bbox.tolist(), + "datetime": datetime, + "collections": [self._param("collection")], + "limit": 1000, + } + items = self._get_pages(data) + + if len(items) == 0: + raise ValueError(f"No items found for {self._param('collection')} between {datetime}.") + + # apply filters + return self._filter_items(items, self._param("filter")) + + def _filter_items(self, items, filter): + if filter is not None and len(filter) > 0: + for option in filter: + for k, v in filter[option].items(): + if k == "eq": + items = [item for item in items if item["properties"][option] == v] + elif k == "ueq": + items = [item for item in items if item["properties"][option] != v] + elif k == "in": + # value should be list + if not isinstance(v, list): + raise ValueError(f"Filter option {k} needs a list.") + items = [item for item in items if item["properties"][option] in v] + elif k == "lt": + items = [item for item in items if float(item["properties"][option]) < v] + elif k == "gt": + items = [item for item in items if float(item["properties"][option]) > v] + else: + raise ValueError(f"Filter option {k} not supported.") + + if len(items) == 0: + raise ValueError(f"No items found for {self._param('collection')} after filtering.") + + return items + + def _get_pages(self, data): + """loop through the api pages and return all items.""" + items = [] + for i in range(1, 100): + _data = data.copy() + _data["page"] = i + page = requests.post(urljoin(self.base_url, "search"), json=_data).json() + + if "features" not in page: + raise ValueError(f"There was an error with the request: {page}") + if len(page["features"]) == 0: + break + else: + items.extend(page["features"]) + + if i == 100: + raise ValueError( + "Max number of pages reached. Consider using a smaller time frame." + ) + + return items + + def download(self, items): + """Download the items from Copernicus Data Space Ecosystem as xr.Dataset or download the files. + + :param items: items to download + :return: xarray.Dataset or list of filenames + """ + if self._param("create_minicube"): + ds = self._download_to_minicube( + items, + self._param("shp"), + self._param("collection"), + self._param("bands", raise_error=True), + self._param("resolution"), + self._param("resampling"), + self._param("filter_asset_path"), + self._param("use_virtual_rasterio_file"), + ) + ds = self._prepare_cube(ds) + return ds + else: + return self._download_to_files( + items, + self._param("shp"), + self._param("collection"), + self._param("bands", raise_error=True), + self._param("resolution"), + self._param("resampling"), + self._param("filter_asset_path"), + self._param("use_virtual_rasterio_file"), + ) + + def _download_to_files( + self, + items, + shp, + collection, + bands, + resolution, + resampling, + filter_asset_path, + use_virtual_rasterio_file, + ): + """Download all the items and return the file paths.""" + fns = Parallel(n_jobs=self._param("num_workers"))( + delayed(self._download_to_file)( + item, + shp, + collection, + band, + resolution, + resampling, + filter_asset_path, + use_virtual_rasterio_file, + ) + for item, band in itertools.product(items, bands) + ) + # extract list of lists + fns = [fn for sublist in fns for fn in sublist] + return fns + + def _download_to_file( + self, + item, + shp, + collection, + band, + resolution, + resampling, + filter_asset_path, + use_virtual_rasterio_file, + ): + """Download item to file.""" + f_paths = self._get_asset_path(item, collection, band, resolution, filter_asset_path) + # replace extension to tif and collapse folders to name (there can be multiple files with the same name) + fns = [] + for f_path in f_paths: + fn = self._param("download_folder") / Path("_".join(f_path.with_suffix(".tif").parts)) + if not fn.exists(): + clipped = self._download_file(f_path, shp, resampling, use_virtual_rasterio_file) + fn.parent.mkdir(parents=True, exist_ok=True) + clipped.rio.to_raster(fn) + fns.append(fn) + return fns + + def _download_to_minicube( + self, + items, + shp, + collection, + bands, + resolution, + resampling, + filter_asset_path, + use_virtual_rasterio_file, + ): + """Download all items and merge them into a dataset.""" + # do not rely on threading here, because it will mess up the xarrays + datasets = Parallel(n_jobs=self._param("num_workers"))( + delayed(self._download_item)( + item, + shp, + collection, + band, + resolution, + resampling, + filter_asset_path, + use_virtual_rasterio_file, + ) + for item, band in itertools.product(items, bands) + ) + # extract list of lists [[item1band1, item1band2, ...], [item2band1, item2band2, ...], ...] + datasets = [datasets[i : i + len(bands)] for i in range(0, len(datasets), len(bands))] + # combine them to dataset with time data [ds1, ds2, ...] + time_data = Parallel(n_jobs=self._param("num_workers"))( + delayed(self._combine_bands)(datasets[i], shp, bands, resolution, resampling) + for i in range(len(datasets)) + ) + + # combine the time data + if all([ds is None for ds in time_data]): + raise RuntimeError("No items were found.") + + # extract time information for each dataset + if len(time_data) != len(items): + raise RuntimeError("Lengths of downloaded items and requested items do not match.") + # skip items which were not found + time_data, times = map( + list, + zip( + *[ + (ds, item["properties"]["datetime"]) + for ds, item in zip(time_data, items) + if ds is not None + ] + ), + ) + + time_data = align_coords(time_data, shp, resampling) + + # add time coords (would have been removed by reproject_match in align_coords) + time_data = [ + ds.assign_coords(time=("time", pd.to_datetime([time]).tz_convert(None))) + for ds, time in zip(time_data, times) + ] + + # merge time + data = xr.concat(time_data, dim="time", join="exact").sortby("time") + return data + + def _download_item( + self, + item, + shp, + collection, + band, + resolution, + resampling, + filter_asset_path, + use_virtual_rasterio_file, + ): + """Download all bands of an item and merge them into a dataset.""" + # go through the bands + f_paths = self._get_asset_path(item, collection, band, resolution, filter_asset_path) + if len(f_paths) > 1: + # if there are multiple assets, merge them + datasets = [] + for f_path in f_paths: + ds = self._download_file(f_path, shp, resampling, use_virtual_rasterio_file) + if ds is None: + continue + if "band" not in ds.coords: + # add band dim + ds = ds.expand_dims("band") + # rename band each band coord with filename to make sure they are unique + identifyable later + ds = ds.assign_coords( + band=[f"{f_path.stem}_{old}" for old in ds.coords["band"].values] + ) + datasets.append(ds) + if len(datasets) > 1: + warnings.warn( + f"Multiple files found for band {band}: {f_paths}.\nAdding them as new bands." + ) + if not all([ds.rio.crs == datasets[0].rio.crs for ds in datasets]): + datasets = align_coords(datasets, shp, resampling) + if not all( + [ds.rio.resolution()[0] == datasets[0].rio.resolution()[0] for ds in datasets] + ): + datasets = self._align_resolutions(datasets, shp, resolution, resampling) + # add them as new bands + clipped = xr.concat(datasets, dim="band") + elif len(datasets) == 0: + clipped = None + else: + clipped = datasets[0] + else: + # single asset + clipped = self._download_file(f_paths[0], shp, resampling, use_virtual_rasterio_file) + + return clipped + + def _combine_bands(self, band_data, shp, bands, resolution, resampling): + # check if all bands were found + if len(band_data) != len(bands): + raise RuntimeError("Length of band_data and bands do not match.") + if all([ds is None for ds in band_data]): + warnings.warn("Not all bands were downloaded.") + return None + if not all([ds is not None for ds in band_data]): + raise RuntimeError("Bands were not downloaded.") + + # resample if needed + band_data = self._align_resolutions(band_data, shp, resolution, resampling) + + # pad everything to make sure it has shape of shp + band_data = [ + ds.rio.pad_box(*list(shp.total_bounds)).rio.clip_box(*list(shp.total_bounds)) + for ds in band_data + ] + + # merge bands + if not len(band_data) > 0: + return + # remove band dimension + for i, ds in enumerate(band_data): + if "band" in ds.dims: + if len(ds.band) > 1: + band_data[i] = band_data[i].to_dataset(dim="band") + band_data[i] = band_data[i].rename( + {old: f"{bands[i]}_{old}" for old in band_data[i].data_vars} + ) + else: + band_data[i] = band_data[i].drop_vars("band").squeeze("band") + band_data[i] = band_data[i].to_dataset(name=bands[i]) + ds = xr.combine_by_coords(band_data) + + return ds + + def _align_resolutions(self, datasets, shp, resolution, resampling): + """unify the resolutions of the list of the xr.Datasets.""" + ress = [ds.rio.resolution() for ds in datasets] + resolution = meters_to_crs_unit(resolution, shp) + # round degrees to 8 decimal places for cm resolution + ress = [(round(res[0], 8), round(res[0], 8)) for res in ress] + resolution = [round(res, 8) for res in resolution] + if len(set(ress)) > 1 or (len(ress) > 0 and ress[0] != resolution): + # get closest resolution + idx = sorted( + enumerate(ress), + key=lambda item: abs(resolution[0] - abs(item[1][0])) + + abs(resolution[1] - abs(item[1][1])), + )[0][0] + if ress[idx] != resolution: + datasets[idx] = datasets[idx].rio.reproject( + shp.crs, resolution=resolution, resampling=resampling + ) + # reproject rest to master + datasets = [ + (ds.rio.reproject_match(datasets[idx], resampling=resampling) if i != idx else ds) + for i, ds in enumerate(datasets) + ] + return datasets + + def _download_file(self, f_path, shp, resampling, use_virtual_rasterio_file): + """wrapper to decide whether to use rasterio or patch download.""" + if use_virtual_rasterio_file: + return self._download_file_rasterio(f_path, shp, resampling) + else: + return self._download_file_tile(f_path, shp, resampling) + + def _download_file_rasterio(self, f_path, shp, resampling): + """Download a band of an item and clip it to the shapefile.""" + session = rasterio.session.AWSSession( + aws_unsigned=False, + endpoint_url=( + urlparse(self.end_point_url).netloc + if "://" in self.end_point_url + else self.end_point_url + ), + aws_access_key_id=self.credentials["aws_access_key_id"], + aws_secret_access_key=self.credentials["aws_secret_access_key"], + ) + with rasterio.env.Env(session=session, AWS_VIRTUAL_HOSTING=False): + clipped = self._clip_to_region("s3://eodata/" + str(f_path), shp, resampling) + return clipped + + def _download_file_tile(self, f_path, shp, resampling): + """Download a band of an item and clip it to the shapefile.""" + # download the file locally + download_path = self._param("download_folder") / f_path + download_path.parent.mkdir(parents=True, exist_ok=True) + if not download_path.exists(): + _s3 = boto3.Session( + aws_access_key_id=self.credentials["aws_access_key_id"], + aws_secret_access_key=self.credentials["aws_secret_access_key"], + region_name="default", + ).resource("s3", endpoint_url=self.end_point_url) + + _s3.Bucket("eodata").download_file(str(f_path), download_path) + + # clip to shp + clipped = self._clip_to_region(download_path, shp, resampling) + + # optionally remove the downloaded file + if self._param("rm_tmp_files"): + download_path.unlink(missing_ok=True) + + return clipped + + def _get_asset_path(self, item, collection, band, resolution, filter_asset_path): + # set up session to read file structure + _s3 = boto3.Session( + aws_access_key_id=self.credentials["aws_access_key_id"], + aws_secret_access_key=self.credentials["aws_secret_access_key"], + region_name="default", + ).resource("s3", endpoint_url=self.end_point_url) + + # extract item path + try: + s3_path = item["assets"]["PRODUCT"]["alternate"]["s3"]["href"] + except KeyError: + raise RuntimeError( + "It seems that no s3 path exists for this item. Returned item: ", item + ) + folder_name = "/".join(s3_path.split("/")[2:]) + response = _s3.Bucket("eodata").objects.filter(Prefix=folder_name) + + # filter for extension + paths = [ + obj for obj in response if any([obj.key.endswith(x) for x in self._file_extensions]) + ] + if len(paths) == 0: + raise RuntimeError("No file with valid extension found.") + + # transform to Path object + paths = [Path(obj.key) for obj in paths] + + # filter for band + if band is not None and len(band) > 0: + paths_new = [path for path in paths if band in path.name] + if len(paths_new) > 0: + paths = paths_new + else: + # try again with lower and upper case + paths = [ + path for path in paths if band.lower() in path.name or band.upper() in path.name + ] + if len(paths) == 0: + raise RuntimeError(f"Band {band} not found, is it written correctly?") + + # apply regex path filters + if filter_asset_path and collection in filter_asset_path: + pattern = re.compile(filter_asset_path[collection]) + paths = [path for path in paths if re.search(pattern, str(path))] + if len(paths) == 0: + raise RuntimeError( + "There are no files matching the filter_asset_path: ", + filter_asset_path, + ) + + paths_new = [ + path for path in paths if f"{resolution}m" in str(path) + ] # filter for resolution (currently only applies for S2) + if len(paths_new) > 0: # if it was found take it, otherwise will be resampled later + paths = paths_new + + if len(paths) == 0: + raise RuntimeError("No file found.") + + return paths + + def _clip_to_region(self, file, shp, resampling): + try: + ds = rxr.open_rasterio(file) + gcps = ds.rio.get_gcps() + if gcps and not any(c in ds.coords for c in ["x", "y"]): + with rasterio.open(file) as src: + gcps, crs = src.get_gcps() + # use WarpedVRT to get the correct transform/coordinates + with WarpedVRT(src, src_crs=crs, resampling=resampling) as vrt: + ds = rxr.open_rasterio(vrt) + if ds.rio.crs is None: + warnings.warn("No crs found, continuing with EPSG:4326.") + ds = ds.rio.write_crs("EPSG:4326") + src_crs = ds.rio.crs + if src_crs != shp.crs: + # reproject shp and clip with margin + shp_crs = shp.to_crs(src_crs) + # create bbox with 10% margin + margin = ( + (shp_crs.bounds.maxx.item() - shp_crs.bounds.minx.item()) * 0.1, + (shp_crs.bounds.maxy.item() - shp_crs.bounds.miny.item()) * 0.1, + ) + shapely_box = box( + shp_crs.bounds.minx.item() - margin[0], + shp_crs.bounds.miny.item() - margin[1], + shp_crs.bounds.maxx.item() + margin[0], + shp_crs.bounds.maxy.item() + margin[1], + ) + gdf = gpd.GeoDataFrame(geometry=[shapely_box], crs=src_crs) + ds = ds.rio.clip_box(*list(gdf.total_bounds)) + # then reproject to shp crs + ds = ds.rio.reproject(shp.crs, resampling=resampling) + ds = ds.rio.clip_box(*list(shp.total_bounds)) + else: + ds = ds.rio.clip_box(*list(shp.total_bounds)) + ds.load() + ds.close() + except rxr.exceptions.NoDataInBounds: + warnings.warn("No data found in bounds.") + ds = None + + return ds diff --git a/terragon/google_earth_engine.py b/terragon/google_earth_engine.py index 92a6b21..a0b0101 100644 --- a/terragon/google_earth_engine.py +++ b/terragon/google_earth_engine.py @@ -1,7 +1,6 @@ import hashlib import json import math -import re import warnings from typing import List, Union @@ -11,10 +10,11 @@ import rioxarray as rxr import xarray as xr from joblib import Parallel, delayed +from rasterio.enums import Resampling from rasterio.transform import from_origin from .base import Base -from .utils import meters_to_crs_unit, rm_files +from .utils import align_coords, meters_to_crs_unit, rm_files class GEE(Base): @@ -25,6 +25,10 @@ class GEE(Base): :param credentials: unused, kept for compatibility, defaults to None """ + _GEE_ID_PROP_NAME = "system:id" + _GEE_DATE_PROP_NAME = "system:time_start" + _GEE_ADD_BAND = "FILL_MASK" + def __init__(self, credentials: dict = None) -> None: """Initialize class and GEE. @@ -47,14 +51,14 @@ def retrieve_collections(self, filter_by_name: str = None) -> None: "GEE does not have a collection endpoint. Please, visit https://developers.google.com/earth-engine/datasets/catalog" ) - def search(self, rm_tmp_files=True, **kwargs) -> ee.ImageCollection: - """Search for items in the GEE collections. For a description of the kwargs parameters see the Base class function. + def search(self, *args, rm_tmp_files=True, **kwargs) -> ee.ImageCollection: + """Search for items in the GEE collections. For a description of the args/kwargs parameters see the Base class function. :param rm_tmp_files: remove temporarily downloaded files after creating the minicube, defaults to True :raises ValueError: when parameters are missing or in the wrong format :return: ee.ImageCollection """ - super().search(**kwargs) + super().search(*args, **kwargs) self._parameters.update({"rm_tmp_files": rm_tmp_files}) img_col = ee.ImageCollection(self._param("collection")) @@ -128,25 +132,28 @@ def download(self, img_col: ee.ImageCollection) -> Union[xr.Dataset, List]: def _download_img(self, img_col, i, tmp_dir, shp, region, transform, shape): """Download a single image from the GEE ImageCollection.""" img = ee.Image(img_col.get(i)) - # get the system id - id_prop = next( - (prop for prop in img.propertyNames().getInfo() if "system:id" in prop), - None, - ) - if not id_prop: + id_prop = self._get_img_property(img, self._GEE_ID_PROP_NAME) + date_prop = self._get_img_property(img, self._GEE_DATE_PROP_NAME) + if id_prop is None: warnings.warn( f"Could not find system:id property in image {i}. \ - Using consecutive numbers of images, but this can lead to problems wiht overwriting files." + Using consecutive numbers of images, but this can lead to problems with overwriting files." ) - img_id = i + id_prop = i else: - img_id = img.get(id_prop).getInfo() # replace the / with _ to avoid problems with file paths - img_id = img_id.replace("/", "_") + id_prop = id_prop.replace("/", "_") + if date_prop is None: + warnings.warn( + f"Could not find system:time_start property in image {i}. \ + Using the current date, but this can lead to problems." + ) + # current date in ms + date_prop = int(pd.Timestamp.now().timestamp() * 1000) # create a unique filename through geometry since we are downloading clipped images geom_hash = hashlib.sha256(shp.geometry.iloc[0].wkt.encode("utf-8")).hexdigest() - fileName = tmp_dir.joinpath(f"{img_id}_{geom_hash}.tif") + fileName = tmp_dir.joinpath(f"{date_prop}_{id_prop}_{geom_hash}.tif") if not fileName.exists(): img = geedim.MaskedImage(img) img.download( @@ -158,28 +165,43 @@ def _download_img(self, img_col, i, tmp_dir, shp, region, transform, shape): ) return fileName + def _get_img_property(self, img, prop_name): + """Get a property from an image.""" + prop = img.get(prop_name).getInfo() + if prop is None: + # fallback to in + id_prop = next( + (prop for prop in img.propertyNames().getInfo() if prop_name in prop), + None, + ) + if id_prop is not None: + prop = img.get(id_prop).getInfo() + return prop + def _merge_gee_tifs(self, fns) -> xr.Dataset: """merge the tifs and crop them to the shp""" if len(fns) < 1: raise ValueError("No files provided to merge.") - date_pattern = r"\d{8}" def load_tif(fn): da = rxr.open_rasterio(fn) - time_str = re.findall(date_pattern, str(fn))[0] - da = da.assign_coords(time=pd.to_datetime(time_str, format="%Y%m%d")) + # first string is date, see _download_img + time_str = fn.name.split("_")[0] + da = da.assign_coords(time=pd.to_datetime(int(time_str), unit="ms")) return da out = Parallel(n_jobs=self._param("num_workers"), backend="threading")( delayed(load_tif)(fn) for fn in fns ) - ds = xr.concat(out, dim="time").compute() + out = align_coords(out, self._param("shp"), Resampling.nearest) + + ds = xr.concat(out, dim="time") ds = ds.sortby("time") ds = ds.to_dataset(dim="band") ds = ds.rename_vars( {dim: name for dim, name in zip(ds.data_vars.keys(), ds.attrs["long_name"])} ) - if "FILL_MASK" in ds.data_vars: - ds = ds.drop_vars("FILL_MASK") + if self._GEE_ADD_BAND in ds.data_vars: + ds = ds.drop_vars(self._GEE_ADD_BAND) return ds diff --git a/terragon/init.py b/terragon/init.py index 0b65122..0003a3b 100644 --- a/terragon/init.py +++ b/terragon/init.py @@ -15,5 +15,9 @@ def init(api: str, credentials: dict = None, **kwargs) -> object: from .google_earth_engine import GEE return GEE(credentials, **kwargs) + elif api == "cdse" or api == "copernicus_data_space_ecosystem": + from .copernicus_data_space_ecosystem import CDSE + + return CDSE(credentials, **kwargs) else: - raise ValueError(f'API {api} not supported. Please use "pc" or "gee".') + raise ValueError(f'API {api} not supported. Please use "pc", "gee", or "cdse".') diff --git a/terragon/microsoft_planetary_computer.py b/terragon/microsoft_planetary_computer.py index d88d62a..c86a8b2 100644 --- a/terragon/microsoft_planetary_computer.py +++ b/terragon/microsoft_planetary_computer.py @@ -60,13 +60,13 @@ def retrieve_collections(self, filter_by_name: str = None) -> list: else: raise RuntimeError("Failed to retrieve collections") - def search(self, **kwargs): - """Search for items in the Planetary Computer collections. For a description of the kwargs parameters see the Base class function. + def search(self, *args, **kwargs): + """Search for items in the Planetary Computer collections. For a description of the args/kwargs parameters see the Base class function. :raises ValueError: when no items are found or parameters are in the wrong format :return: a list of items """ - super().search(**kwargs) + super().search(*args, **kwargs) bounds_4326 = self._reproject_shp(self._param("shp")).total_bounds catalog = pystac_client.Client.open( diff --git a/terragon/utils.py b/terragon/utils.py index 3913172..698014e 100644 --- a/terragon/utils.py +++ b/terragon/utils.py @@ -1,7 +1,56 @@ +import warnings + import pyproj from shapely.geometry import Point +def indices_are_identical(datasets: list) -> bool: + """Check if the indices from the xarray datasets are equal.""" + if len(datasets) == 1: + return True + + # get all indices + coords = list(set(key for ds in datasets for key in list(ds.indexes.keys()))) + + # check if all indices exist in all datasets + if not all(all(coord in ds.indexes for ds in datasets) for coord in coords): + return False + + # check if all indices are equal + if not all( + all(datasets[0].indexes[coord].equals(ds.indexes[coord]) for ds in datasets) + for coord in coords + ): + return False + + return True + + +def align_coords(datasets, shp, resampling): + """unify the crs and indices of the datasets. Often the coordinates differ by a very small amount + due to rounding errors, xarray will merge these into one dataframe, here the unify them to prevent + empty pixels in the dataset.""" + # make sure the coordinates are the same and match (e.g. if there are other crs) + if not all([ds.rio.crs == datasets[0].rio.crs for ds in datasets]) or not indices_are_identical( + datasets + ): + crs_list = [ds.rio.crs for ds in datasets] + # match all to master (first shp crs) + idx = [i for i, crs in enumerate(crs_list) if crs == shp.crs] + if len(idx) == 0: + warnings.warn( + f"No matching crs found and all crs are different. Using first crs {crs_list[0]} as master." + ) + idx = 0 + else: + idx = idx[0] + datasets = [ + (ds.rio.reproject_match(datasets[idx], resampling=resampling) if i != idx else ds) + for i, ds in enumerate(datasets) + ] + return datasets + + def rm_files(fns): for fn in fns: if fn.exists(): diff --git a/tests/.env b/tests/.env index 5a31c33..9a116b8 100644 --- a/tests/.env +++ b/tests/.env @@ -1 +1,3 @@ GEE_PROJECT_NAME='' # Store your GEE project name for running tests here +S3_ACCESS_KEY='' # copernicus s3 access key +S3_SECRET_KEY='' # copernicus s3 secret key \ No newline at end of file diff --git a/tests/copernicus_data_space_ecosystem.py b/tests/copernicus_data_space_ecosystem.py new file mode 100644 index 0000000..5957f47 --- /dev/null +++ b/tests/copernicus_data_space_ecosystem.py @@ -0,0 +1,274 @@ +import os +import unittest + +import geopandas as gpd +from base import _TestBase +from shapely.geometry import Polygon +from utils import load_env_variables + +import terragon + + +class TestCDSE(_TestBase, unittest.TestCase): + def setUp(self): + super().setUp() + + load_env_variables() # load the .env vars if running locally + credentials = { + "aws_access_key_id": os.environ.get("S3_ACCESS_KEY"), + "aws_secret_access_key": os.environ.get("S3_SECRET_KEY"), + } + self.tg = terragon.init("cdse", credentials) + self.arguments["collection"] = "SENTINEL-2" + self.arguments["bands"] = ["B02", "B03", "B04"] + self.arguments["filter"] = {"processingLevel": {"eq": "S2MSI2A"}} + + def test_error_on_no_band(self): + """test error if no band is given""" + args = self.arguments.copy() + args["bands"] = [] + self.assertRaises(ValueError, self.tg.create, **args) + + def test_s2_mosaic(self): + """test Sentinel-2 data""" + args = self.arguments.copy() + args["collection"] = "GLOBAL-MOSAICS" + args["bands"] = ["B02"] + args["start_date"] = "2022-01-01" + args["end_date"] = "2022-07-31" + args["resolution"] = 10 + args["filter"] = {"platformShortName": {"eq": "SENTINEL-2"}} + + ds = self.tg.create(**args) + + self.assertTrue( + len(ds.time) == 3 + and self.width - 1 <= len(ds.x) <= self.width + 1 + and self.height - 1 <= len(ds.y) <= self.height + 1 + ) + + def test_s1(self): + """test Sentinel-1 data""" + args = self.arguments.copy() + args["collection"] = "SENTINEL-1" + args["bands"] = ["VV", "VH"] + args["start_date"] = "2022-01-01" + args["end_date"] = "2022-01-04" + args["resolution"] = 10 + args["filter"] = {"productType": {"eq": "IW_GRDH_1S-COG"}} + + ds = self.tg.create(**args) + + # test if the dimensions are correct, with +/- 1 pixel distance + self.assertTrue( + len(ds.time) == self.nr_time_steps + and self.width - 1 <= len(ds.x) <= self.width + 1 + and self.height - 1 <= len(ds.y) <= self.height + 1 + ) + + # not supported file format for raw data + args["filter"] = {"productType": {"eq": "IW_RAW__0S"}} # only dat files + self.assertRaises(RuntimeError, self.tg.create, **args) + + def test_multiple_files_s1(self): + """test Sentinel-1 data with SLC, which have multiple files.""" + args = self.arguments.copy() + args["collection"] = "SENTINEL-1" + args["bands"] = ["VV"] + args["start_date"] = "2022-01-01" + args["end_date"] = "2022-01-03" + args["resolution"] = 10 + args["filter"] = {"productType": {"eq": "IW_SLC__1S"}} + # use other region to include two files -> test only for 2 vars and time + args["shp"] = gpd.GeoDataFrame( + geometry=[ + Polygon( + [ + (9.295, 47.425), + (9.296, 47.425), + (9.296, 47.426), + (9.295, 47.426), + ] + ) + ], + crs="EPSG:4326", + ) + + ds = self.tg.create(**args) + + self.assertTrue(len(ds.data_vars) == 2 and len(ds.time) == 1) + + def test_s1_rtc(self): + """test Sentinel-1-RTC data""" + args = self.arguments.copy() + coords = [ + [ + [39.618673, -14.899018], + [39.618673, -14.891632], + [39.631033, -14.891632], + [39.631033, -14.899018], + [39.618673, -14.899018], + ] + ] + polygon = Polygon(coords[0]) + args["shp"] = gpd.GeoDataFrame(index=[0], crs="EPSG:4326", geometry=[polygon]) + args["collection"] = "SENTINEL-1-RTC" + args["bands"] = ["VV"] + args["start_date"] = "2018-02-25" + args["end_date"] = "2018-02-25" + args["resolution"] = 20 + args["filter"] = {} + + ds = self.tg.create(**args) + + nr_time_steps = 2 + # do not test width and height because the shape was taken arbitrarily from the first example + self.assertTrue(len(ds.time) == nr_time_steps) + + def test_s2_landcover(self): + """test Sentinel-2 landcover data with RGB file.""" + args = self.arguments.copy() + args["collection"] = "S2GLC" + args["bands"] = ["S2GLC"] + args["start_date"] = "2019-01-01" + args["end_date"] = "2019-12-31" + args["resolution"] = 10 + args["filter"] = {} + # args["filter_asset_path"] = {"S2GLC": "^(?!.*RGB).*"} + args["filter_asset_path"] = {"S2GLC": ".*RGB.tif"} + + ds = self.tg.create(**args) + + self.assertTrue( + len(ds.data_vars) == 3 + and len(ds.time) == 1 + and self.width - 1 <= len(ds.x) <= self.width + 1 + and self.height - 1 <= len(ds.y) <= self.height + 1 + ) + + def test_merge_bands_multiband_and_band_s2_landcover(self): + """test Sentinel-2 landcover data with both RGB and normal file + this uses two files, one with 3 bands and one with 1 band""" + args = self.arguments.copy() + args["collection"] = "S2GLC" + args["bands"] = ["S2GLC"] + args["start_date"] = "2019-01-01" + args["end_date"] = "2019-12-31" + args["resolution"] = 10 + args["filter"] = {} + + ds = self.tg.create(**args) + + self.assertTrue(len(ds.data_vars) == 4) + + self.assertTrue( + len(ds.time) == 1 + and self.width - 1 <= len(ds.x) <= self.width + 1 + and self.height - 1 <= len(ds.y) <= self.height + 1 + ) + + def test_dem(self): + args = self.arguments.copy() + args["collection"] = "COP-DEM" + args["bands"] = ["DEM"] + args["start_date"] = "2010-01-01" + args["end_date"] = "2010-12-31" + args["resolution"] = 90 + args["filter"] = { + "spatialResolution": {"eq": 90}, + "productType": {"eq": "DGE_90"}, + } + + ds = self.tg.create(**args) + + self.assertTrue( + len(ds.time) == 1 + and self.width // 9 - 1 <= len(ds.x) <= self.width // 9 + 1 + and self.height // 9 - 1 <= len(ds.y) <= self.height // 9 + 1 + ) + + # does contain two product types: DGE_30 and DTE_30 + args["resolution"] = 30 + args["filter"] = {"spatialResolution": {"eq": 30}} + ds = self.tg.create(**args) + + self.assertTrue( + len(ds.time) == 2 + and self.width // 3 - 1 <= len(ds.x) <= self.width // 3 + 1 + and self.height // 3 - 1 <= len(ds.y) <= self.height // 3 + 1 + ) + + def test_modis(self): + args = self.arguments.copy() + args["collection"] = "TERRAAQUA" + args["bands"] = [""] + args["start_date"] = "2021-01-01" + args["end_date"] = "2021-01-16" + args["resolution"] = 500 + args["filter"] = {} + + # it seems that no item has assets or a download url (tested on 15.11.24) + self.assertRaises(RuntimeError, self.tg.create, **args) + + def test_l5(self): + args = self.arguments.copy() + coords = [ + [ + [25.74646, 53.944763], + [25.808258, 53.944763], + [25.808258, 53.918878], + [25.74646, 53.918878], + [25.74646, 53.944763], + ] + ] + polygon = Polygon(coords[0]) + args["shp"] = gpd.GeoDataFrame(index=[0], crs="EPSG:4326", geometry=[polygon]) + args["collection"] = "LANDSAT-5" + args["bands"] = ["B1"] + args["start_date"] = "1984-12-08" + args["end_date"] = "1984-12-08" + args["resolution"] = 120 + args["filter"] = {} + + nr_time_steps = 2 + + ds = self.tg.create(**args) + + # only test time steps since the shape was taken arbitrarily from the first sample + self.assertTrue(len(ds.time) == nr_time_steps) + + def test_l7(self): + args = self.arguments.copy() + coords = [[[32.9, 62.5], [33.0, 62.5], [33.0, 62.49], [32.9, 62.49], [32.9, 62.5]]] + polygon = Polygon(coords[0]) + args["shp"] = gpd.GeoDataFrame(index=[0], crs="EPSG:4326", geometry=[polygon]) + args["collection"] = "LANDSAT-7" + args["bands"] = ["B1"] + args["start_date"] = "1999-11-01" + args["end_date"] = "1999-11-01" + args["resolution"] = 30 + args["filter"] = {} + + nr_time_steps = 2 + + ds = self.tg.create(**args) + + # only test time steps since the shape was taken arbitrarily from the first sample + self.assertTrue(len(ds.time) == nr_time_steps) + + def test_l8(self): + args = self.arguments.copy() + args["collection"] = "LANDSAT-8-ESA" + args["bands"] = ["B2"] + args["start_date"] = "2021-01-01" + args["end_date"] = "2021-01-16" + args["resolution"] = 30 + args["filter"] = {"processingLevel": {"eq": "LEVEL2"}} + + ds = self.tg.create(**args) + + self.assertTrue(len(ds.time) == 4 and len(ds.x) == 10 and len(ds.y) == 7) + + +if __name__ == "__main__": + unittest.main()