diff --git a/README.md b/README.md
index e450921..2b097ba 100644
--- a/README.md
+++ b/README.md
@@ -149,3 +149,22 @@ planet_vi_list: list = ['NDVI','AVI','CIRE','NIRv','NDWI']
## examples
see the "notebook" sub-folder for detailed jupyter notebook examples
+
+## Version history
+
+The current development version is **0.8.0**. Previous releases are listed below in reverse chronological order.
+
+| Version | Release date |
+| --- | --- |
+| 0.8.0 | Current |
+| 0.7.0 | 2026-08-31 |
+| 0.6.0 | 2026-06-05 |
+| 0.5.0 | 2026-01-19 |
+| 0.4.5 | 2026-01-05 |
+| 0.4.1 | 2025-12-04 |
+| 0.3.9 | 2025-12-04 |
+| 0.3.5 | 2025-11-05 |
+| 0.3.3 | 2025-10-19 |
+| 0.2.26 | 2025-09-25 |
+| 0.2.15 | 2025-05-05 |
+| 0.2.9 | 2025-04-03 |
\ No newline at end of file
diff --git a/notebooks/google_drive_access.ipynb b/notebooks/google_drive_access.ipynb
deleted file mode 100644
index beb27ff..0000000
--- a/notebooks/google_drive_access.ipynb
+++ /dev/null
@@ -1,212 +0,0 @@
-{
- "cells": [
- {
- "metadata": {},
- "cell_type": "markdown",
- "source": "# Using the Google Drive in Python and openEO environment",
- "id": "ed30f45ec0b3dea4"
- },
- {
- "metadata": {},
- "cell_type": "markdown",
- "source": [
- "In WEED we want to use files from the GoogleDrive for the processing - mainly in the benchmarking part of the project. MOreover, some of the results of the benchmarking should be directly uploaded to the GDrive and therefore available for all user. The following steps had to be implemented to access the GDrive from Python.\n",
- "\n",
- "1. Sign in to the WEED Google Account and open the Google Cloud Console to create a project.
\n",
- " - go to https://console.cloud.google.com/
\n",
- " - Create a New Project: Click on “Select a project” at the top and then “New Project”. Give it a name and create. Name is \"WEED-2024\"\n",
- "2. enable the Google Drive API for the WEED project
\n",
- " - Navigate to APIs & Services: In the left sidebar, go to “APIs & Services” > “Dashboard”.\n",
- " - Enable APIs: Click on “Enable APIs and Services”. Search for “Google Drive API” and enable it for your project.\n",
- "3. create credentials\n",
- " - Generate Credentials: In the left sidebar, go to “APIs & Services” > “Credentials”.\n",
- " - Create Credentials: Click “Create Credentials” > “Service Account”. Fill in details, choose a role, and create a JSON key. Save this JSON key securely.\n",
- " - Name: gdrive_access\n",
- " - Email: gdrive-access@weed-2024.iam.gserviceaccount.com (automatic generation)\n",
- " - enable the service account status\n",
- "4. JSON access key\n",
- " - Download and Use JSON Key: Store the downloaded JSON key securely on your local machine. This key will be used for authenticating your application to access the Drive API. (The service account key can only be retrieved the first time)\n",
- " - Therefore, if you need a new key just add a new key to the gdrive_access credentials!!!!!\n",
- " - BEST: store the key in a secure location like the VITO valut and access it from there client-side (https://confluence.vito.be/pages/viewpage.action?spaceKey=EP&title=Vault+user+guide)\n",
- " - link to the VITO vault which you can access with your TERRASCOPE credentials: https://vault.vgt.vito.be/ui/vault/auth?with=ldap\n",
- "5. make sure that for each file and folder which should be available in Python the servie account email address is added to the \"access\" of these files and folders. Currently only the folder \"openeo_tests\" under WEED/working/WP4_ToolboxDvlpt is granted access in this way.\n",
- "6. install the pydrive2 python package in the WEED environment"
- ],
- "id": "30ee94f8ff2cdab2"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": [
- "from pydrive2.fs import GDriveFileSystem\n",
- "import pandas as pd\n",
- "import hvac\n",
- "from getpass import getpass\n",
- "from eo_processing.utils.storage import WEED_storage"
- ],
- "id": "910da53275c67ed2"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": [
- "def get_WEED_credentials(username: str = 'buchhornm', key: str = 'gdrive-access') -> str:\n",
- " \"\"\"\n",
- " Retrieves WEED access credentials from Terrascope VAULT using LDAP authentication.\n",
- "\n",
- " This method prompts the user to enter their password for Terrascope VAULT, authenticates\n",
- " with the VAULT using LDAP, and fetches credentials from the WEED KV storage path.\n",
- "\n",
- " :param username: LDAP username used to authenticate with the VAULT, defaults to 'buchhornm'\n",
- " :param key: Key in the KV WEED storage to get value from, defaults to 'gdrive-access'\n",
- " :return: credentials as a string\n",
- " \"\"\"\n",
- " password_prompt = 'Please enter your password for the Terrascope VAULT: '\n",
- " service_account_password = getpass(prompt=password_prompt)\n",
- "\n",
- " client = hvac.Client(url='https://vault.vgt.vito.be')\n",
- "\n",
- " client.auth.ldap.login(\n",
- " username=username,\n",
- " password=service_account_password,\n",
- " mount_point='ldap'\n",
- " )\n",
- "\n",
- " secret_version_response = client.secrets.kv.v2.read_secret_version(mount_point='kv',\n",
- " path='TAP/apps/WEED',\n",
- " raise_on_deleted_version=True)\n",
- "\n",
- " client.logout()\n",
- "\n",
- " return secret_version_response['data']['data'][key]"
- ],
- "id": "f1e1668457583572"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": [
- "# get the credentials for the GDrive service account access from the VITO TERRASCOPE vault\n",
- "gdrive_credentials = get_WEED_credentials(username='deroob', key='gdrive-access')"
- ],
- "id": "9ab997420e6ebd23"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": [
- "# init the fsspec filesystem to access the files & folders available for the service account credentials\n",
- "# gdrive-access@weed-2024.iam.gserviceaccount.com\n",
- "# \"1k27bitdRp41AtHq1xupyqwKaTLzrMUMu\" is the ID of the only folder currently available for this service account\n",
- "# if more folder should be available then add the email_address to the user of the files and/or folders wished\n",
- "# see: https://filesystem-spec.readthedocs.io/en/latest/usage.html#use-a-file-system for more info to interact with file system\n",
- "\n",
- "gdrive = GDriveFileSystem(\"1k27bitdRp41AtHq1xupyqwKaTLzrMUMu\",\n",
- " use_service_account=True,\n",
- " client_json=gdrive_credentials,)"
- ],
- "id": "ced4ef2536d05a02"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": [
- "# list files and folders\n",
- "for root, dnames, fnames in gdrive.walk(gdrive.root):\n",
- " print(root, dnames, fnames)"
- ],
- "id": "df571ce2b8bce11f"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": [
- "# we can interact now quite easy with the data - example for the CSV file read in pandas\n",
- "# Note: you always have to start from the entrance point and then add the sub-folder plus filenames separated by \"/\"\n",
- "with gdrive.open(gdrive.root + \"/\" + 'SK_v5_reference-points_EUNIS2012.csv', 'rb') as f:\n",
- " df = pd.read_csv(f)\n",
- "df.head()\n"
- ],
- "id": "53ea31f7674d07fb"
- },
- {
- "metadata": {},
- "cell_type": "markdown",
- "source": "Now we're going to try the same things but with the WEED_storage object",
- "id": "ff8ed55e5808349"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": "storage = WEED_storage(username='deroob')",
- "id": "c9654e69efd8721f"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": [
- "import json\n",
- "gdrive = GDriveFileSystem(\"1k27bitdRp41AtHq1xupyqwKaTLzrMUMu\",\n",
- " use_service_account=True,\n",
- " client_json=json.dumps(storage.gdrive_credentials),)"
- ],
- "id": "93a1d9795fb3fb2"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": [
- "for root, dnames, fnames in gdrive.walk(gdrive.root):\n",
- " print(root, dnames, fnames)"
- ],
- "id": "390bcf4c3981f702"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": "",
- "id": "3d843e42223f6f7b"
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "weed",
- "language": "python",
- "name": "weed"
- },
- "language_info": {
- "codemirror_mode": {
- "name": "ipython",
- "version": 2
- },
- "file_extension": ".py",
- "mimetype": "text/x-python",
- "name": "python",
- "nbconvert_exporter": "python",
- "pygments_lexer": "ipython2",
- "version": "2.7.6"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/notebooks/storage_access.ipynb b/notebooks/storage_access.ipynb
deleted file mode 100644
index f4aba38..0000000
--- a/notebooks/storage_access.ipynb
+++ /dev/null
@@ -1,48 +0,0 @@
-{
- "cells": [
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": [
- "from eo_processing.utils.storage import WEED_storage\n",
- "from eo_processing.config import generate_storage_options"
- ],
- "id": "7d6499a360e610a5"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": [
- "# define the output folder (make sure you adapt this to your folder structure)\n",
- "#define if S3 workspace is needed\n",
- "workspace_path = f\"tests/Bert/Bert_test_v{str(2)}\"\n",
- "storage = WEED_storage(username='deroob',s3_bucket='test')\n",
- "storage_options = generate_storage_options(workspace_export = True, S3_prefix=workspace_path, local_S3_needed = True, storage=storage)\n"
- ],
- "id": "75033be000cc7f05"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": "sonata_storage = WEED_storage(username='deroob',s3_bucket='blabla',project='sonata')",
- "id": "b40ad56024036845"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": "",
- "id": "4eae5b169c64aa00"
- }
- ],
- "metadata": {},
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/notebooks/test_nobs.ipynb b/notebooks/test_nobs.ipynb
deleted file mode 100644
index eae896e..0000000
--- a/notebooks/test_nobs.ipynb
+++ /dev/null
@@ -1,460 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "id": "d46e36dcd1f7a0ef",
- "metadata": {},
- "source": [
- "# example how to generate the WEED standard feature cube based on EO data (Sentinel-1 and Sentinel-2)\n",
- "this tests include the generation of the feature cube with and without NVBT band. The NVBT band is a measure of the number of valid input data timesteps after cloud masking and the internal temporal binning.\n",
- "The NVBT band can be useful for identifying areas with high or low levels of observation, which can be important for a variety of applications, such as monitoring changes in land use or assessing the quality of data in a given area."
- ]
- },
- {
- "cell_type": "code",
- "id": "initial_id",
- "metadata": {
- "ExecuteTime": {
- "end_time": "2026-07-01T10:10:45.301457300Z",
- "start_time": "2026-07-01T10:10:43.592557200Z"
- }
- },
- "source": [
- "from eo_processing.utils.helper import init_connection\n",
- "from eo_processing.openeo.processing import generate_master_feature_cube\n",
- "from eo_processing.config.settings import get_advanced_options, get_job_options, get_collection_options"
- ],
- "outputs": [],
- "execution_count": 1
- },
- {
- "cell_type": "markdown",
- "id": "c9abe4dec39f3145",
- "metadata": {},
- "source": [
- "#### declare space and time"
- ]
- },
- {
- "cell_type": "code",
- "id": "1c319454e1c292f1",
- "metadata": {
- "ExecuteTime": {
- "end_time": "2026-07-01T10:10:48.565131400Z",
- "start_time": "2026-07-01T10:10:48.548481800Z"
- }
- },
- "source": [
- "# the time context is given by start and end date\n",
- "year = 2024\n",
- "start = f'{year}-01-01'\n",
- "end = f'{year+1}-01-01' # the end is always exclusive\n",
- "\n",
- "# the space context is defined as a bounding box dictionary with south,west,north,east and crs\n",
- "# we take as example an 10x10km tile in EU LAEA grid around Vienna\n",
- "AOI = {'west': 4780000, 'east': 4790000, 'south': 2830000, 'north': 2840000, 'crs': 3035}"
- ],
- "outputs": [],
- "execution_count": 2
- },
- {
- "cell_type": "markdown",
- "id": "39fc1ee260fa1492",
- "metadata": {},
- "source": [
- "### get processing_options for the eo_processing functions, collection_options and job_options"
- ]
- },
- {
- "cell_type": "code",
- "id": "fa16159203be2475",
- "metadata": {
- "ExecuteTime": {
- "end_time": "2026-07-01T10:11:06.858284600Z",
- "start_time": "2026-07-01T10:11:06.842773300Z"
- }
- },
- "source": [
- "processing_options = get_advanced_options(provider='cdse', skip_check_S1=False, skip_check_S2=True)\n",
- "job_options = get_job_options(provider='cdse', task='feature_generation')\n",
- "collection_options = get_collection_options(provider='cdse')\n",
- "processing_options.update({'openeo_chunk_size': 64})"
- ],
- "outputs": [],
- "execution_count": 3
- },
- {
- "cell_type": "code",
- "id": "72616c9c-4583-4132-a38c-23988aeb8434",
- "metadata": {
- "ExecuteTime": {
- "end_time": "2026-07-01T10:11:08.354813700Z",
- "start_time": "2026-07-01T10:11:08.305942900Z"
- }
- },
- "source": [
- "processing_options"
- ],
- "outputs": [
- {
- "data": {
- "text/plain": [
- "{'provider': 'cdse',\n",
- " 's1_orbitdirection': 'DESCENDING',\n",
- " 'target_crs': 3035,\n",
- " 'resolution': 10.0,\n",
- " 'time_interpolation': False,\n",
- " 'ts_interval': 'dekad',\n",
- " 'S2_temporal_reducer': 'median',\n",
- " 'S1_temporal_reducer': 'mean',\n",
- " 'SLC_masking_algo': 'mask_scl_dilation',\n",
- " 'S2_max_cloud_cover': 95,\n",
- " 'S2_bands': ['B02',\n",
- " 'B03',\n",
- " 'B04',\n",
- " 'B05',\n",
- " 'B06',\n",
- " 'B07',\n",
- " 'B08',\n",
- " 'B8A',\n",
- " 'B11',\n",
- " 'B12'],\n",
- " 's2_tileid_list': None,\n",
- " 'skip_check_S1': False,\n",
- " 'skip_check_S2': True,\n",
- " 'apply_cloud_mask': True,\n",
- " 'get_NVBT': False,\n",
- " 'optical_vi_list': ['ABDI1',\n",
- " 'ABDI2',\n",
- " 'AWEInsh',\n",
- " 'AVI',\n",
- " 'BLFEI',\n",
- " 'CIRE',\n",
- " 'EVI',\n",
- " 'IRECI',\n",
- " 'MBWI',\n",
- " 'MNDWI',\n",
- " 'MNDVI',\n",
- " 'NDMI',\n",
- " 'NDVI',\n",
- " 'NDVIMNDWI',\n",
- " 'NDWI',\n",
- " 'NMDI',\n",
- " 'NIRv',\n",
- " 'S2WI',\n",
- " 'S2REP',\n",
- " 'WRI'],\n",
- " 'radar_vi_list': ['VHVVD', 'VHVVR', 'DpRVIVV'],\n",
- " 'S2_scaling': [0, 10000, 0, 1.0],\n",
- " 'S1_db_rescale': True,\n",
- " 'append': True,\n",
- " 'openeo_chunk_size': 64}"
- ]
- },
- "execution_count": 4,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "execution_count": 4
- },
- {
- "cell_type": "markdown",
- "id": "a05d287923cf93eb",
- "metadata": {},
- "source": [
- "### establish connection to openEO"
- ]
- },
- {
- "cell_type": "code",
- "id": "d06f89bea20dce40",
- "metadata": {
- "ExecuteTime": {
- "end_time": "2026-07-01T10:11:15.268379200Z",
- "start_time": "2026-07-01T10:11:14.743118500Z"
- }
- },
- "source": "con = init_connection(provider='cdse')",
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Authenticated using refresh token.\n"
- ]
- }
- ],
- "execution_count": 5
- },
- {
- "cell_type": "markdown",
- "id": "adcab1031cc5e6c2",
- "metadata": {},
- "source": [
- "### run the feature cube generation WITHOUT NVBT band"
- ]
- },
- {
- "cell_type": "code",
- "id": "f32c4461ec4331a5",
- "metadata": {
- "ExecuteTime": {
- "end_time": "2026-07-01T10:11:21.277193500Z",
- "start_time": "2026-07-01T10:11:21.260122500Z"
- }
- },
- "source": [
- "# update job_options due toBerts setting from last inference runs\n",
- "# ToDO: optimize the job settings for feature_cube_generation_with_nobs & feature_cube_generation and put in settings + add correct task to 'get_job_options'\n",
- "job_options.update({\n",
- " \"driver-memory\": \"4G\",\n",
- " \"driver-memoryOverhead\": \"4G\",\n",
- " \"executor-memory\": \"5G\",\n",
- " \"executor-memoryOverhead\": \"3g\",\n",
- " \"max-executors\": 10,\n",
- " \"python-memory\": \"disable\",\n",
- " \"allow_empty_cubes\": True,\n",
- " \"soft-errors\": 0.05})"
- ],
- "outputs": [],
- "execution_count": 6
- },
- {
- "cell_type": "code",
- "id": "fc107ebf279499cf",
- "metadata": {},
- "source": [
- "# get master cube without nobs\n",
- "data = generate_master_feature_cube(con, AOI, start, end, **collection_options, **processing_options)"
- ],
- "outputs": [],
- "execution_count": null
- },
- {
- "cell_type": "code",
- "id": "f0b7c0577d597daf",
- "metadata": {},
- "source": "data.execute_batch(r'C:\\Users\\buchhorm\\Downloads\\test_cube\\features_cube_v5.tif', title='feature without nobs (10x10km)', job_options=job_options)",
- "outputs": [],
- "execution_count": null
- },
- {
- "cell_type": "markdown",
- "id": "3f775e9b75ec1382",
- "metadata": {},
- "source": [
- "## run with activated NVBT generation"
- ]
- },
- {
- "cell_type": "code",
- "id": "95a0e6270cf65cc3",
- "metadata": {},
- "source": [
- "# now we run same with nobs_perc band\n",
- "processing_options.update({'get_NVBT': True})\n",
- "data2 = generate_master_feature_cube(con, AOI, start, end, **collection_options, **processing_options)"
- ],
- "outputs": [],
- "execution_count": null
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": "data2.execute_batch(r'C:\\Users\\buchhorm\\Downloads\\test_cube\\features_cube_with_nobs_v5.tif', title='feature with nobs (10x10km)', job_options=job_options)",
- "id": "b522dd0e9b4adaa5"
- },
- {
- "metadata": {},
- "cell_type": "markdown",
- "source": [
- "## now we run this 10x10km test in a bigger context - we generate the full EO feature cube with NVBT band\n",
- "Note: since the new ECDC version and alphaEarth STACs are not complete THIS loading has to be tested later"
- ],
- "id": "3684f8b64901068a"
- },
- {
- "cell_type": "code",
- "id": "bcb5d1b835c659eb",
- "metadata": {
- "ExecuteTime": {
- "end_time": "2026-07-01T10:12:04.872738900Z",
- "start_time": "2026-07-01T10:11:27.488000800Z"
- }
- },
- "source": [
- "from habitat_mapping.openeo.feature_cubes import create_EOfeature_cube_WEED_V1\n",
- "processing_options.update(target_crs = 3035)\n",
- "processing_options.update({'get_NVBT': True})\n",
- "processing_options.update({'openeo_chunk_size': 16})\n",
- "job_options.update({\"allow_empty_cubes\": True})\n",
- "job_options.update({\n",
- " \"driver-memory\": \"4G\",\n",
- " \"driver-memoryOverhead\": \"4G\",\n",
- " \"executor-memory\": \"5G\",\n",
- " \"executor-memoryOverhead\": \"3G\",\n",
- " \"max-executors\": 10,\n",
- " \"python-memory\": \"disable\",\n",
- " \"soft-errors\": 0.05,\n",
- "})\n",
- "data3 = create_EOfeature_cube_WEED_V1(con, AOI, start, end, collection_options, processing_options)"
- ],
- "outputs": [
- {
- "name": "stderr",
- "output_type": "stream",
- "text": [
- "Deriving band listing from unordered `item_assets`\n",
- "The specified bands ['precipitation-flux', 'temperature-mean'] in `load_stac` are not a subset of the bands [] found in the STAC metadata (unknown bands: ['precipitation-flux', 'temperature-mean']). Working with specified bands as is.\n"
- ]
- }
- ],
- "execution_count": 7
- },
- {
- "metadata": {
- "ExecuteTime": {
- "end_time": "2026-07-01T10:42:53.609735700Z",
- "start_time": "2026-07-01T10:12:16.483423100Z"
- }
- },
- "cell_type": "code",
- "source": "data3.execute_batch(r'C:\\Users\\buchhorm\\Downloads\\test_cube\\eo_feature_cube_beta2_10x10.tif', title='eo_feature cube (10x10km)', job_options=job_options)",
- "id": "641442d062a77b1e",
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "0:00:00 Job 'j-2607011012174934aeec7870e2685f99': send 'start'\n",
- "0:00:07 Job 'j-2607011012174934aeec7870e2685f99': queued (progress 0%)\n",
- "0:00:13 Job 'j-2607011012174934aeec7870e2685f99': queued (progress 0%)\n",
- "0:00:19 Job 'j-2607011012174934aeec7870e2685f99': queued (progress 0%)\n",
- "0:00:27 Job 'j-2607011012174934aeec7870e2685f99': queued (progress 0%)\n",
- "0:00:37 Job 'j-2607011012174934aeec7870e2685f99': queued (progress 0%)\n",
- "0:00:50 Job 'j-2607011012174934aeec7870e2685f99': queued (progress 0%)\n",
- "0:01:05 Job 'j-2607011012174934aeec7870e2685f99': running (progress 8.9%)\n",
- "0:01:24 Job 'j-2607011012174934aeec7870e2685f99': running (progress 11.5%)\n",
- "0:01:48 Job 'j-2607011012174934aeec7870e2685f99': running (progress 14.5%)\n",
- "0:02:19 Job 'j-2607011012174934aeec7870e2685f99': running (progress 18.0%)\n",
- "0:02:56 Job 'j-2607011012174934aeec7870e2685f99': running (progress 22.0%)\n",
- "0:03:43 Job 'j-2607011012174934aeec7870e2685f99': running (progress 26.5%)\n",
- "0:04:41 Job 'j-2607011012174934aeec7870e2685f99': running (progress 31.4%)\n",
- "0:05:41 Job 'j-2607011012174934aeec7870e2685f99': running (progress 35.8%)\n",
- "0:06:41 Job 'j-2607011012174934aeec7870e2685f99': running (progress 39.7%)\n",
- "0:07:42 Job 'j-2607011012174934aeec7870e2685f99': running (progress 43.1%)\n",
- "0:08:42 Job 'j-2607011012174934aeec7870e2685f99': running (progress 46.2%)\n",
- "0:09:42 Job 'j-2607011012174934aeec7870e2685f99': running (progress 49.0%)\n",
- "0:10:43 Job 'j-2607011012174934aeec7870e2685f99': running (progress 51.5%)\n",
- "0:11:43 Job 'j-2607011012174934aeec7870e2685f99': running (progress 53.7%)\n",
- "0:12:43 Job 'j-2607011012174934aeec7870e2685f99': running (progress 55.8%)\n",
- "0:13:44 Job 'j-2607011012174934aeec7870e2685f99': running (progress 57.7%)\n",
- "0:14:44 Job 'j-2607011012174934aeec7870e2685f99': running (progress 59.4%)\n",
- "0:15:44 Job 'j-2607011012174934aeec7870e2685f99': running (progress 61.0%)\n",
- "0:16:44 Job 'j-2607011012174934aeec7870e2685f99': running (progress 62.4%)\n",
- "0:17:45 Job 'j-2607011012174934aeec7870e2685f99': running (progress 63.8%)\n",
- "0:18:45 Job 'j-2607011012174934aeec7870e2685f99': running (progress 65.1%)\n",
- "0:19:45 Job 'j-2607011012174934aeec7870e2685f99': running (progress 66.3%)\n",
- "0:20:45 Job 'j-2607011012174934aeec7870e2685f99': running (progress 67.4%)\n",
- "0:21:45 Job 'j-2607011012174934aeec7870e2685f99': running (progress 68.4%)\n",
- "0:22:45 Job 'j-2607011012174934aeec7870e2685f99': running (progress 69.4%)\n",
- "0:23:46 Job 'j-2607011012174934aeec7870e2685f99': running (progress 70.3%)\n",
- "0:24:46 Job 'j-2607011012174934aeec7870e2685f99': running (progress 71.1%)\n",
- "0:25:46 Job 'j-2607011012174934aeec7870e2685f99': running (progress 72.0%)\n",
- "0:26:47 Job 'j-2607011012174934aeec7870e2685f99': running (progress 72.7%)\n",
- "0:27:47 Job 'j-2607011012174934aeec7870e2685f99': running (progress 73.5%)\n",
- "0:28:49 Job 'j-2607011012174934aeec7870e2685f99': running (progress 74.2%)\n",
- "0:29:50 Job 'j-2607011012174934aeec7870e2685f99': finished (progress 100%)\n"
- ]
- },
- {
- "data": {
- "text/plain": [
- ""
- ],
- "text/html": [
- "\n",
- " \n",
- " \n",
- " \n",
- " \n",
- " "
- ]
- },
- "execution_count": 8,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "execution_count": 8
- },
- {
- "metadata": {},
- "cell_type": "markdown",
- "source": "## now upscale the the final 20x20km tiles",
- "id": "3d9134ce13226561"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": [
- "# we take as example an 20x20km tile in EU LAEA grid around Vienna\n",
- "AOI = {'west': 4780000, 'east': 4800000, 'south': 2820000, 'north': 2840000, 'crs': 3035}\n",
- "# increasing driver memory\n",
- "job_options.update({\n",
- " \"driver-memory\": \"10G\",\n",
- " \"driver-memoryOverhead\": \"4G\",\n",
- "})\n",
- "data4 = create_EOfeature_cube_WEED_V1(con, AOI, start, end, collection_options, processing_options)"
- ],
- "id": "b7fe80d4f70e9b87"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "outputs": [],
- "execution_count": null,
- "source": "data4.execute_batch(r'C:\\Users\\buchhorm\\Downloads\\test_cube\\eo_feature_cube_beta2_20x20.tif', title='eo_feature cube (20x20km)', job_options=job_options)",
- "id": "a2be83f50c1aa7d0"
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "weed",
- "language": "python",
- "name": "weed"
- },
- "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.14.3"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/scripts/generate_openeo_pg.py b/scripts/generate_openeo_pg.py
deleted file mode 100644
index abbe063..0000000
--- a/scripts/generate_openeo_pg.py
+++ /dev/null
@@ -1,157 +0,0 @@
-#%%
-import os
-import json
-import openeo
-import sys
-
-sys.path.append(r'C:\Git_projects\eo_processing\src')
-from eo_processing.openeo.preprocessing import ts_datacube_extraction, extract_S1_datacube, extract_S2_datacube, S2_BANDS
-
-# Establish OpenEO connection (adjust if needed for your environment)
-connection = openeo.connect("openeo.dataspace.copernicus.eu").authenticate_oidc()
-
-# Define the spatial and temporal extent for the tests
-spatial_extent = {
- 'east': 4880000,
- 'south': 2898000,
- 'west': 4878000,
- 'north': 2900000,
- 'crs': 'EPSG:3035'
-}
-temporal_extent = ["2021-01-01", "2022-01-01"]
-
-# Directory to save the process graphs
-output_directory = "resources/unit-tests"
-os.makedirs(output_directory, exist_ok=True)
-
-#%%
-
-def generate_process_graphs(function, scenarios, connection, bbox, start, end, output_dir):
- for filename, processing_options in scenarios:
- try:
- # Create the data cube with the specified processing options
- datacube = function(
- connection=connection,
- bbox=bbox,
- start=start,
- end=end,
- **processing_options
- )
-
- # Convert the process graph to JSON format
- process_graph_json = datacube.to_json()
-
- # Save the process graph JSON to a file
- output_path = os.path.join(output_dir, filename)
- with open(output_path, "w") as json_file:
- json_file.write(process_graph_json)
- print(f"Process graph saved: {output_path}")
-
- except ValueError as e:
- # If an invalid option was provided, log the error (this is expected for certain tests)
- print(f"Error generating process graph for {filename}: {e}")
-
-#%%
-
-# Define the different test scenarios and processing options
-ts_test_scenarios = [
-
- # Sentinel-2 Only, Basic
- ("ts_datacube_extraction_S2.json", {"S1_collection": None}),
-
- # Combined Sentinel-1 and Sentinel-2
- ("ts_datacube_extraction_combined.json", {}),
-
- # Sentinel-2 with Cloud Masking
- ("ts_datacube_extraction_S2_with_masking.json", {"SLC_masking_algo": "mask_scl_dilation"}),
-
- # Temporal Aggregation and Interpolation
- ("ts_datacube_extraction_S1_interpolation.json", {
- "S1_collection": "SENTINEL1_GRD",
- "ts_interval": "P1M",
- "time_interpolation": True
- }),
-
- # Combined with Custom CRS and Resolution
- ("ts_datacube_extraction_combined_custom_crs.json", {
- "target_crs": "EPSG:3857",
- "resolution": 20.0
- })
-]
-
-# Define test scenarios for `extract_S1_datacube`
-s1_test_scenarios = [
- # Basic Sentinel-1 Data Extraction
- ("extract_S1_basic.json", {"S1_collection": "SENTINEL1_GRD"}),
-
- # Sentinel-1 with Temporal Aggregation
- ("extract_S1_temporal_aggregation.json", {
- "S1_collection": "SENTINEL1_GRD",
- "ts_interval": "P1M"
- }),
-
- # Sentinel-1 with Resampling and Custom CRS
- ("extract_S1_custom_crs.json", {
- "S1_collection": "SENTINEL1_GRD",
- "target_crs": "EPSG:3857",
- "resolution": 30.0
- })
-]
-
-# Define test scenarios for `extract_S2_datacube`
-s2_test_scenarios = [
- # Basic Sentinel-2 Data Extraction
- ("extract_S2_basic.json", {"S2_collection": "SENTINEL2_L2A"}),
-
- # Sentinel-2 with Cloud Masking
- ("extract_S2_with_masking.json", {
- "S2_collection": "SENTINEL2_L2A",
- "SLC_masking_algo": "mask_scl_dilation"
- }),
-
- # Sentinel-2 with Temporal Aggregation and Resampling
- ("extract_S2_temporal_resampling.json", {
- "S2_collection": "SENTINEL2_L2A",
- "ts_interval": "P1M",
- "target_crs": "EPSG:3857",
- "resolution": 20.0
- })
-]
-
-#%%
-
-# Generate process graphs for `ts_datacube_extraction`
-generate_process_graphs(
- function=ts_datacube_extraction,
- scenarios=ts_test_scenarios,
- connection=connection,
- bbox=spatial_extent,
- start=temporal_extent[0],
- end=temporal_extent[1],
- output_dir=output_directory
-)
-
-# Generate process graphs for `extract_S1_datacube`
-generate_process_graphs(
- function=extract_S1_datacube,
- scenarios=s1_test_scenarios,
- connection=connection,
- bbox=spatial_extent,
- start=temporal_extent[0],
- end=temporal_extent[1],
- output_dir=output_directory
-)
-
-# Generate process graphs for `extract_S2_datacube`
-generate_process_graphs(
- function=extract_S2_datacube,
- scenarios=s2_test_scenarios,
- connection=connection,
- bbox=spatial_extent,
- start=temporal_extent[0],
- end=temporal_extent[1],
- output_dir=output_directory
-)
-
-print("All process graphs have been generated successfully!")
-# %%
diff --git a/scripts/global_120x120km_UTM_tiling-grid.py b/scripts/global_120x120km_UTM_tiling-grid.py
deleted file mode 100644
index b692e37..0000000
--- a/scripts/global_120x120km_UTM_tiling-grid.py
+++ /dev/null
@@ -1,120 +0,0 @@
-"""
-this script generates a global 120x120km tiling grid in UTM.
-
-"""
-
-
-'''
-STEPS
-1. load the UTM zones gpkg
-2. load the global land masses gpkg
-3. loop over the 120 UTM zones
- a. load the full 120x120 km grid depending on the hemisphere and change the epsg to the UTM zone
- b. load a copy of the UTM zone polygon and warp to correct epsg
- c. select the grid cells intersecting the UTM zone polygon and delete non-needed
- d. load a copy of the land masses polygon and warp to correct epsg
- e. select the grid cells intersecting the land masses and dlete non-needed ones
- f. run over the remaining 120x120km tiles and add the naming as well as
- add the openEO bbox dict using the bounds of the polygon (plus add the bounds in WGS84)
- g. convert the geodataframe back to EPSG:4326 and save it to a new gpkg
-4. combine all UTM zones into one gpkg and clean up the metadata
-5. save the final global processing grid in the package ressources
-
-'''
-
-import geopandas as gpd
-import pandas as pd
-import math
-from eo_processing.utils.mgrs import UTM_2_grid20id
-import os
-
-# standard paths
-path_grid_n = os.path.normpath(r'C:\Users\buchhorm\Downloads\120x120km_grid\basic_120x120km_grid_nothern_no-crs_v2.gpkg')
-path_grid_s = os.path.normpath(r'C:\Users\buchhorm\Downloads\120x120km_grid\basic_120x120km_grid_southern_no-crs_v2.gpkg')
-path_land = os.path.normpath(r'C:\Users\buchhorm\Downloads\120x120km_grid\land_sea_mask_20kmbuffered_EPSG4326_v2.gpkg')
-path_utm = os.path.normpath(r'C:\Users\buchhorm\Downloads\120x120km_grid\UTM_zones_high-res_EPSG4326.gpkg')
-
-# load UTM zones
-gdf_utm = gpd.read_file(path_utm)
-
-# loop over all UTM zones
-lZones = gdf_utm.name.unique().tolist()
-lFiles = []
-
-for UTMzone in lZones:
- print(f'processing zone: {UTMzone}')
- # get the UTM zone polygon and the EPSG
- if UTMzone[-1] == 'N':
- epsg = 32600 + int(UTMzone[:2])
- else:
- epsg = 32700 + int(UTMzone[:2])
-
- gdf_zone = gdf_utm[gdf_utm.name == UTMzone].copy()
-
- path_out = os.path.normpath(
- r'C:\Users\buchhorm\Downloads\120x120km_grid\results\UTM_zone_{0}.gpkg'.format(
- epsg))
- os.makedirs(os.path.dirname(path_out), exist_ok=True)
-
- if not os.path.exists(path_out):
- # get the land masses
- gdf_land = gpd.read_file(path_land)
- gdf_land = gdf_land.clip(gdf_zone)
- gdf_land = gdf_land.dissolve()
-
- # bring all to correct epsg
- gdf_zone = gdf_zone.to_crs(epsg=epsg)
- gdf_land = gdf_land.to_crs(epsg=epsg)
-
- # short cut if the gdf_land is empty
- if gdf_land.empty:
- continue
-
- # now we load the basic 120x120km grid and set the crs to the correct epsg
- if UTMzone[-1] == 'N':
- gdf_grid = gpd.read_file(path_grid_n).set_crs(epsg=epsg, allow_override=True)
- else:
- gdf_grid = gpd.read_file(path_grid_s).set_crs(epsg=epsg, allow_override=True)
-
- # select grid cells which are intersecting with utm zone
- gdf_grid = gdf_grid[gdf_grid.intersects(gdf_zone.union_all(method='coverage'))]
-
- # intersecting
- gdf_grid = gdf_grid[gdf_grid.intersects(gdf_land.union_all(method='coverage'))]
-
- # add the bbox dict for openEO
- gdf_grid['bbox_dict'] = gdf_grid.apply(
- lambda row: {
- 'west': math.floor(row.geometry.bounds[0] / 100) * 100.,
- 'south': math.floor(row.geometry.bounds[1] / 100) * 100.,
- 'east': math.ceil(row.geometry.bounds[2] / 100) * 100.,
- 'north': math.ceil(row.geometry.bounds[3] / 100) * 100.,
- 'crs': epsg
- }, axis=1
- )
-
- # save to disk
-
-
- gdf_grid[['left', 'top', 'right', 'bottom', 'row_index','col_index', 'bbox_dict', 'geometry']].to_file(
- path_out)
- else:
- # load
- gdf_grid = gpd.read_file(path_out)
- gdf_zone = gdf_zone.to_crs(epsg=epsg)
- gdf_grid = gdf_grid.clip(gdf_zone)
- gdf_grid.to_file(path_out)
-
- lFiles.append(path_out)
-
-
-# load all geopackages which were produced and convert to EPSG:4326 and merge them into one GeoDataFrame
-result = gpd.GeoDataFrame(columns=['left', 'top', 'right', 'bottom', 'row_index','col_index', 'bbox_dict', 'geometry'], geometry='geometry', crs='EPSG:4326')
-
-for file in lFiles:
- gdf_tmp = gpd.read_file(file)
- gdf_tmp = gdf_tmp.to_crs(epsg=4326)
- result = pd.concat([result, gdf_tmp], ignore_index=True)
-
-# save final result to disk
-result.to_file(os.path.normpath(r'C:\Users\buchhorm\Downloads\120x120km_grid\ECDC_global_120x120km_grid.gpkg'))
diff --git a/scripts/global_20x20km_UTM_tiling-grid.py b/scripts/global_20x20km_UTM_tiling-grid.py
deleted file mode 100644
index a5a8029..0000000
--- a/scripts/global_20x20km_UTM_tiling-grid.py
+++ /dev/null
@@ -1,114 +0,0 @@
-"""
-this script generates a global 20x20km tiling grid in UTM using the grid20id system for naming. the tiles are
-optimized to the global land masses within the Sentinel-2 sensing area (including a 10km buffer).
-
-"""
-
-
-'''
-STEPS
-1. load the UTM zones gpkg
-2. load the global land masses gpkg
-3. loop over the 60 UTM zones
- a. load the full 20x20 km grid and change the epsg to the UTM zone
- b. load a copy of the UTM zone polygon and warp to correct epsg
- c. clip the grid to the UTM zone polygon
- d. load a copy of the land masses polygon and warp to correct epsg
- e. clip the grid to the land masses
- f. run over the remaining 20x20km tiles and add the grid20 name (using the centeroid of the polygon) as well as
- add the openEO bbox dict using the bounds of the polygon
- g. convert the geodataframe back to EPSG:4326 and save it to a new gpkg
-4. combine all UTM zones into one gpkg and clean up the metadata
-5. save the final global processing grid in the package ressources
-
-'''
-
-import geopandas as gpd
-import pandas as pd
-import math
-from eo_processing.utils.mgrs import UTM_2_grid20id
-import os
-
-# standard paths
-path_utm = os.path.normpath(r'C:\Users\BUCHHORM\Downloads\global_20x20km_opneEO_processing_grid\global_high-res_UTMzones.gpkg')
-path_land = os.path.normpath(r'C:\Users\BUCHHORM\Downloads\global_20x20km_opneEO_processing_grid\ne_10m_land\Land_masses_10km_buffered_Clipped_S2.gpkg')
-path_grid = os.path.normpath(r'C:\Users\BUCHHORM\Downloads\global_20x20km_opneEO_processing_grid\basic_20x20km_grid_no-crs.gpkg')
-
-# load UTM zones
-gdf_utm = gpd.read_file(path_utm)
-
-# loop over all UTM zones
-lZones = gdf_utm.name.unique().tolist()
-lFiles = []
-
-for UTMzone in lZones:
- print(f'processing zone: {UTMzone}')
- # get the UTM zone polygon and the EPSG
- if UTMzone[-1] == 'N':
- epsg = 32600 + int(UTMzone[:2])
- else:
- epsg = 32700 + int(UTMzone[:2])
-
- gdf_zone = gdf_utm[gdf_utm.name == UTMzone].copy()
-
- # get the land masses
- gdf_land = gpd.read_file(path_land)
- gdf_land = gdf_land.clip(gdf_zone)
- gdf_land = gdf_land.dissolve()
-
- # bring all to correct epsg
- gdf_zone = gdf_zone.to_crs(epsg=epsg)
- gdf_land = gdf_land.to_crs(epsg=epsg)
-
- # short cut if the gdf_land is empty
- if gdf_land.empty:
- continue
-
- # now we load the basic 20x20km grid and set the crs to the correct epsg
- gdf_grid = gpd.read_file(path_grid).set_crs(epsg=epsg, allow_override=True)
-
- # clipping
- gdf_grid = gdf_grid.clip(gdf_zone)
-
- # intersecting
- gdf_grid = gdf_grid[gdf_grid.intersects(gdf_land.union_all(method='coverage'))]
-
- # name the 20x20 grid cells with the grid20id
- if str(epsg)[2] == '6':
- zone_letter = 'Z'
- else:
- zone_letter = 'A'
-
- gdf_grid['grid20id'] = gdf_grid.apply(
- lambda row: UTM_2_grid20id(row['geometry'].centroid.x, row['geometry'].centroid.y, int(str(epsg)[-2:]), zone_letter), axis=1, result_type='expand')
-
- # add the bbox dict for openEO
- gdf_grid['bbox_dict'] = gdf_grid.apply(
- lambda row: {
- 'west': math.floor(row.geometry.bounds[0] / 100) * 100.,
- 'south': math.floor(row.geometry.bounds[1] / 100) * 100.,
- 'east': math.ceil(row.geometry.bounds[2] / 100) * 100.,
- 'north': math.ceil(row.geometry.bounds[3] / 100) * 100.,
- 'crs': epsg
- }, axis=1
- )
-
- # save to disk
- gdf_grid[['grid20id', 'bbox_dict', 'geometry']].to_file(
- os.path.normpath(
- r'C:\Users\BUCHHORM\Downloads\global_20x20km_opneEO_processing_grid\results\UTM_zone_{0}.gpkg'.format(
- epsg)))
- lFiles.append(os.path.normpath(
- r'C:\Users\BUCHHORM\Downloads\global_20x20km_opneEO_processing_grid\results\UTM_zone_{0}.gpkg'.format(epsg)))
-
-
-# load all geopackages which were produced and convert to EPSG:4326 and merge them into one GeoDataFrame
-result = gpd.GeoDataFrame(columns=['grid20id', 'bbox_dict', 'geometry'], geometry='geometry', crs='EPSG:4326')
-
-for file in lFiles:
- gdf_tmp = gpd.read_file(file)
- gdf_tmp = gdf_tmp.to_crs(epsg=4326)
- result = pd.concat([result, gdf_tmp], ignore_index=True)
-
-# save final result to disk
-result.to_file(os.path.normpath(r'C:\Users\BUCHHORM\Downloads\global_20x20km_opneEO_processing_grid\final_global_grid.gpkg'))
diff --git a/scripts/laea_tiling_grids_1_create.py b/scripts/laea_tiling_grids_1_create.py
deleted file mode 100644
index f96f6b7..0000000
--- a/scripts/laea_tiling_grids_1_create.py
+++ /dev/null
@@ -1,212 +0,0 @@
-"""
-This script is used to generate the base LAEA tiling grids in 20, 100 and 50km resolution covering the full panEU
-
-1. step create tiling grid in 100km over all areas via function name-to-bbox
-2. filter by buffered panEU layer
-3. write out (in EPSG:3035 and EPSG:4326)
-4. write out the high res layers in EPSG:3035 and EPSG:4326
-5. create also the 50x50 and 20x20 km variants
-
-"""
-
-from eo_processing.utils.geoprocessing import laea100km_id_to_extent, laea50km_id_to_extent, laea20km_id_to_extent
-import pandas as pd
-import geopandas as gpd
-from shapely.geometry import box
-import os
-from tqdm import tqdm
-
-def create_100k_grid():
- ### start with the 100km tiling grid for EU
- # filter N and E
- N_min = 8
- N_max = 75
- E_min = 9
- E_max = 80
-
- print('start creating combinations')
- # create all combinations
- entries = []
- for N in range(N_min, N_max+1):
- for E in range(E_min, E_max+1):
- name = f'E{E:02d}N{N:02d}'
- entries.append(name)
- print(f'done creating combinations. number of tiles: {len(entries)} ')
-
- print('start creating geodataframe')
- frames = [] # List to collect DataFrames
-
- for tile_id in tqdm(entries, desc='creating tiles'):
- aoi = laea100km_id_to_extent(tile_id)
-
- df = gpd.GeoDataFrame(
- {"name": tile_id, "bbox_dict": str(aoi), "geometry": [box(aoi['west'], aoi['south'], aoi['east'], aoi['north'])]})
- df.crs = aoi['crs']
- frames.append(df)
-
- # Concatenate all frames at once
- result = pd.concat(frames, ignore_index=True)
- gdf_100 = gpd.GeoDataFrame(result, geometry='geometry', crs='EPSG:3035')
- print('done creating geodataframe')
-
- print('filter tiles to panEU')
- #load the panEU shapefile
- gdf_panEU = gpd.read_file(r'C:\Users\buchhorm\Downloads\new_grids\EU_biogeographic_buffered_final.gpkg')
- #convert to EPSG:3035
- gdf_panEU = gdf_panEU.to_crs(epsg=3035)
-
- #intersect
- gdf_final = gdf_100[gdf_100.name.isin(gdf_100.clip(gdf_panEU).name.unique().tolist())]
-
- print('write out')
- #write out
- gdf_final.to_file(r'C:\Users\buchhorm\Downloads\new_grids\LAEA_100km_tiling_grid_EU_EPSG3035.gpkg')
-
- print('create variants')
- print('create high res version')
- gdf_high = gdf_final.copy()
- gdf_high['geometry'] = gdf_high['geometry'].apply(lambda x: x.segmentize(250))
- gdf_high.to_file(r'C:\Users\buchhorm\Downloads\new_grids\LAEA_100km_tiling_grid_EU_high_res_EPSG3035.gpkg')
-
- print('convert to 4326 versions')
- # create also the EPSG:4326 versions
- gdf_high4326 = gdf_high.to_crs(epsg=4326)
- gdf_high4326.to_file(r'C:\Users\buchhorm\Downloads\new_grids\LAEA_100km_tiling_grid_EU_high_res_EPSG4326.gpkg')
- gdf_final4326 = gdf_final.to_crs(epsg=4326)
- gdf_final4326.to_file(r'C:\Users\buchhorm\Downloads\new_grids\LAEA_100km_tiling_grid_EU_EPSG4326.gpkg')
-
-path_100k = r'C:\Users\buchhorm\Downloads\new_grids\LAEA_100km_tiling_grid_EU_EPSG3035.gpkg'
-if os.path.isfile(path_100k):
- gdf_100k = gpd.read_file(path_100k)
-else:
- create_100k_grid()
- gdf_100k = gpd.read_file(path_100k)
-
-##### now we create the 50x50km versions
-def create_50k_grid(gdf100):
- # we get a list of unique name100
- ltiles = gdf100['name'].unique().tolist()
- print('create 50km tiles')
- # now we loop over all tiles and create all possible 50x50km identifier to create the 50x50km grid
- l50kmtiles = []
- for tile in ltiles:
- parts = tile.lstrip('E').split('N')
- # create the 4 sub-tiles
- l50kmtiles.append(f'E{parts[0]}0N{parts[1]}0')
- l50kmtiles.append(f'E{parts[0]}5N{parts[1]}0')
- l50kmtiles.append(f'E{parts[0]}0N{parts[1]}5')
- l50kmtiles.append(f'E{parts[0]}5N{parts[1]}5')
-
- # create pandas dataframe with the 50x50km tiles
- gdf_50 = pd.DataFrame(l50kmtiles, columns=['name'])
-
- # now we create the bbox_dict column
- gdf_50['bbox_dict'] = gdf_50.apply(lambda row: laea50km_id_to_extent(row['name']), axis=1)
-
- # create the bbox coordinates from the bbox_dict column
- gdf_50['geometry'] = gdf_50['bbox_dict'].apply(
- lambda row: box(row['west'], row['south'], row['east'], row['north']))
-
- # convert to geopandas
- gdf_50 = gpd.GeoDataFrame(gdf_50, geometry='geometry', crs='EPSG:3035')
-
- #intersetc to further filter
- print('filter tiles to panEU')
- #load the panEU shapefile
- gdf_panEU = gpd.read_file(r'C:\Users\buchhorm\Downloads\new_grids\EU_biogeographic_buffered_final.gpkg')
- #convert to EPSG:3035
- gdf_panEU = gdf_panEU.to_crs(epsg=3035)
-
- #intersect
- gdf_final = gdf_50[gdf_50.name.isin(gdf_50.clip(gdf_panEU).name.unique().tolist())]
-
-
- print('write out')
- #write out
- gdf_final.to_file(r'C:\Users\buchhorm\Downloads\new_grids\LAEA_50km_tiling_grid_EU_EPSG3035.gpkg')
-
- print('create variants')
- print('create high res version')
- gdf_high = gdf_final.copy()
- gdf_high['geometry'] = gdf_high['geometry'].apply(lambda x: x.segmentize(250))
- gdf_high.to_file(r'C:\Users\buchhorm\Downloads\new_grids\LAEA_50km_tiling_grid_EU_high_res_EPSG3035.gpkg')
-
- print('convert to 4326 versions')
- # create also the EPSG:4326 versions
- gdf_high4326 = gdf_high.to_crs(epsg=4326)
- gdf_high4326.to_file(r'C:\Users\buchhorm\Downloads\new_grids\LAEA_50km_tiling_grid_EU_high_res_EPSG4326.gpkg')
- gdf_final4326 = gdf_final.to_crs(epsg=4326)
- gdf_final4326.to_file(r'C:\Users\buchhorm\Downloads\new_grids\LAEA_50km_tiling_grid_EU_EPSG4326.gpkg')
-
-path_50k = r'C:\Users\buchhorm\Downloads\new_grids\LAEA_50km_tiling_grid_EU_EPSG3035.gpkg'
-if os.path.isfile(path_50k):
- gdf_50k = gpd.read_file(path_50k)
-else:
- create_50k_grid(gdf_100k)
- gdf_50k = gpd.read_file(path_50k)
-
-# create the 20x20km tiling grid
-def create_20k_grid(gdf100, gdf50):
- # we get a list of unique name100
- ltiles = gdf100['name'].unique().tolist()
- print('create 20km tiles')
- # now we loop over all tiles and create all possible 50x50km identifier to create the 50x50km grid
- l20kmtiles = []
- for tile in ltiles:
- parts = tile.lstrip('E').split('N')
- # create the 25 sub-tiles
- for x in range(0, 10, 2):
- for y in range(0, 10, 2):
- l20kmtiles.append(f'E{parts[0]}{x}N{parts[1]}{y}')
-
-
- # create pandas dataframe with the 50x50km tiles
- gdf_20 = pd.DataFrame(l20kmtiles, columns=['name'])
-
- # now we create the bbox_dict column
- gdf_20['bbox_dict'] = gdf_20.apply(lambda row: laea20km_id_to_extent(row['name']), axis=1)
-
- # create the bbox coordinates from the bbox_dict column
- gdf_20['geometry'] = gdf_20['bbox_dict'].apply(
- lambda row: box(row['west'], row['south'], row['east'], row['north']))
-
- # convert to geopandas
- gdf_20 = gpd.GeoDataFrame(gdf_20, geometry='geometry', crs='EPSG:3035')
-
- #intersetc to further filter
- print('filter tiles to panEU')
- #load the panEU shapefile
- gdf_panEU = gpd.read_file(r'C:\Users\buchhorm\Downloads\new_grids\EU_biogeographic_buffered_final.gpkg')
- #convert to EPSG:3035
- gdf_panEU = gdf_panEU.to_crs(epsg=3035)
-
- #intersect
- gdf_final = gdf_20[gdf_20.name.isin(gdf_20.clip(gdf_panEU).name.unique().tolist())]
-
- # also do intersect with 50km grid since that is also optimized
- gdf_final = gdf_final[gdf_final.name.isin(gdf_final.clip(gdf_50k).name.unique().tolist())]
-
-
- print('write out')
- #write out
- gdf_final.to_file(r'C:\Users\buchhorm\Downloads\new_grids\LAEA_20km_tiling_grid_EU_EPSG3035.gpkg')
-
- print('create variants')
- print('create high res version')
- gdf_high = gdf_final.copy()
- gdf_high['geometry'] = gdf_high['geometry'].apply(lambda x: x.segmentize(250))
- gdf_high.to_file(r'C:\Users\buchhorm\Downloads\new_grids\LAEA_20km_tiling_grid_EU_high_res_EPSG3035.gpkg')
-
- print('convert to 4326 versions')
- # create also the EPSG:4326 versions
- gdf_high4326 = gdf_high.to_crs(epsg=4326)
- gdf_high4326.to_file(r'C:\Users\buchhorm\Downloads\new_grids\LAEA_20km_tiling_grid_EU_high_res_EPSG4326.gpkg')
- gdf_final4326 = gdf_final.to_crs(epsg=4326)
- gdf_final4326.to_file(r'C:\Users\buchhorm\Downloads\new_grids\LAEA_20km_tiling_grid_EU_EPSG4326.gpkg')
-
-path_20k = r'C:\Users\buchhorm\Downloads\new_grids\LAEA_20km_tiling_grid_EU_EPSG3035.gpkg'
-if os.path.isfile(path_20k):
- gdf_20k = gpd.read_file(path_20k)
-else:
- create_20k_grid(gdf_100k, gdf_50k)
- gdf_20k = gpd.read_file(path_20k)
diff --git a/scripts/laea_tiling_grids_2_high_res_Sentinel2_tiles.py b/scripts/laea_tiling_grids_2_high_res_Sentinel2_tiles.py
deleted file mode 100644
index 4444e26..0000000
--- a/scripts/laea_tiling_grids_2_high_res_Sentinel2_tiles.py
+++ /dev/null
@@ -1,89 +0,0 @@
-"""
-extract from the official Sentinel2 tiling grid KML file the UTM bounding box
-
-we create a high resolution S2 tiling grid in EPSG:3035 optimized for panEU
-
-"""
-
-import bs4 as bs
-import pandas as pd
-import geopandas as gpd
-import shapely.wkt
-
-def get_epsg(tile_id):
- """ get epsg code from the Sentinel-2 tile name"""
- # calculate the UTM zone_number and zone_letter from lon, lat
- zone_number = int(tile_id[:2])
- zone_letter = tile_id[2]
-
- # figure out from zone_letter if N or S
- northern = (zone_letter >= 'N')
-
- # get EPSG number from zone number and northern
- if northern:
- target_EPSG = 32600 + zone_number
- else:
- target_EPSG = 32700 + zone_number
-
- return target_EPSG
-
-# filter to tiles in UTM zones needed for pan EU
-#df2 = pd.read_csv(r'C:\Users\buchhorm\Downloads\new_grids\utm_codes_europe.csv')
-#df2['utm'] = df2['ZONE'].astype(str) + df2['ROW_']
-#lNeeded = df2['utm'].tolist()
-
-
-## Start read out everything we need with beautiful soup from Sentinel-2 KML
-xml_file = r'C:\Users\buchhorm\Downloads\120x120km_grid\S2_tiles.kml'
-soup = bs.BeautifulSoup(open(xml_file), 'html.parser')
-
-#find all VRTRasterBands in the soup
-tiles = soup.find_all('description')
-
-data = []
-
-for element in tiles:
-
- text = bs.BeautifulSoup(element.text, 'html.parser')
- try:
- rows = text.table.find_all('tr')
- except:
- continue
- for row in rows:
- try:
- cols = row.find_all('td')
- except:
- continue
- cols = [ele.text.strip() for ele in cols]
- if cols[0] == 'TILE_ID': data1 = cols[1]
- if cols[0] == 'UTM_WKT': data2 = cols[1]
-
- data.append((data1, data2, get_epsg(data1)))
-
-df = pd.DataFrame(data, columns=['tile_id', 'geometry', 'epsg'])
-data = None
-soup = None
-tiles = None
-
-# convert the geometry string into a shapely geometry
-df['geometry'] = df['geometry'].apply(lambda x: shapely.wkt.loads(x))
-
-# now we convert the UTM into EPSG:3035 and combine all
-lEPSG = df.epsg.unique().tolist()
-
-list_results: list[gpd.GeoDataFrame] = []
-for tile in lEPSG:
- print(f'run EPSG zone: {tile}')
- df3 = df[df.epsg == tile].copy()
- gdf: gpd.GeoDataFrame = gpd.GeoDataFrame(df3, geometry=df3.geometry, crs=f'EPSG:{tile}')
-
- # add extra points in 250m intervall
- gdf['geometry'] = gdf['geometry'].apply(lambda x: x.segmentize(250))
-
- gdf2 = gdf.to_crs(epsg=4326)
- list_results.append(gdf2)
-
-result = gpd.GeoDataFrame(pd.concat(list_results, ignore_index=True), crs=list_results[0].crs)
-
-# write out
-result.to_file(r'C:\Users\buchhorm\Downloads\120x120km_grid\Sentinel2_tiling_grid_high_res_EPSG4326.gpkg', driver='GPKG')
diff --git a/scripts/laea_tiling_grids_3_assignment_S2tiles.py b/scripts/laea_tiling_grids_3_assignment_S2tiles.py
deleted file mode 100644
index 110630f..0000000
--- a/scripts/laea_tiling_grids_3_assignment_S2tiles.py
+++ /dev/null
@@ -1,548 +0,0 @@
-"""
-we assign each LAEA grid cell one or multiple Sentinel-2 tileIDs which are needed to load the LAEA tile
-
-Note: we only run this approach for 50km and 20km tiles.... the 100K is created out of the lists of the
- 20K and 50K grids (use one and check against the other)
-
-"""
-
-import geopandas as gpd
-import pandas as pd
-from shapely.ops import unary_union
-from eo_processing.utils.mgrs import LL_2_MGRSid
-import itertools
-
-### declaration
-# load the files
-gdf_s2 = gpd.read_file(r"C:\Users\buchhorm\Downloads\new_grids\Sentinel2_tiling_grid_EU_high_res_EPSG3035.gpkg")
-
-gdf_laea = gpd.read_file(r"C:\Users\BUCHHORM\Downloads\new_grids\LAEA_20km_tiling_grid_EU_high_res_EPSG3035.gpkg")
-path_out = r'C:\Users\buchhorm\Downloads\new_grids\S2_tile_info_20K_grid.gpkg'
-
-results = []
-
-# run over LAEA grid
-for row in gdf_laea.itertuples():
- print(f'* process LAEA tile {row.name}')
- #filter the s2 tiles to buffered BBOX of LAEA grid
- aoi = row.geometry.buffer(100)
- xmin, ymin, xmax, ymax = aoi.bounds
- gdf_s2_c = gdf_s2.copy()
- s2_aoi = gdf_s2_c.cx[xmin:xmax, ymin:ymax]
-
- # case when no S2 intersecting tiles are found for grid cell - error
- if s2_aoi.empty:
- print(f'--- error: no S2 tile matches the bounds of LAEA grid {row.name}')
- results.append([row.name, 'error', None, None])
- continue
-
- # SINGLE WINNER
- if s2_aoi.shape[0] == 1:
- print(f' - single winner in first attempt.')
- results.append([row.name, 'single', s2_aoi.tile_id.iloc[0], s2_aoi.tile_id.iloc[0]])
- continue
-
- # DECISION FOR MULTIPLE SINGLE WINNER
- single_match = []
- for tile in s2_aoi.itertuples():
- if tile.geometry.contains(row.geometry):
- single_match.append([row.name, 'single', tile.tile_id, tile.tile_id])
- if len(single_match) == 0:
- pass
- elif len(single_match) == 1:
- print(f' - single winner in second attempt.')
- results.append(single_match[0])
- continue
- elif len(single_match) > 1:
- print(f' - multiple single winners - we chose the best match')
- # calculate the MGRS identifier of centroid of LAEA grid cell
- ## get the latlon of centroid of the LAEA grid cell
- lon, lat = gpd.GeoSeries(row.geometry.centroid, crs=3035).to_crs(epsg=4326).get_coordinates().iloc[0]
- MGRSid = LL_2_MGRSid(lon, lat)
-
- iDone = False
- # try first if we have a direct hit to MGRS grid
- for element in single_match:
- if element[2] == MGRSid:
- results.append(element)
- iDone = True
- break
- if iDone:
- continue
-
- # now we try if we get a zone + band hit
- for element in single_match:
- if element[2][:3] == MGRSid[:3]:
- results.append(element)
- iDone = True
- break
- if iDone:
- continue
-
- # now we try just a zone hit
- for element in single_match:
- if element[2][:2] == MGRSid[:2]:
- results.append(element)
- iDone = True
- break
- if iDone:
- continue
-
- # damm - fuck it we just take the first hit from the single_match list
- results.append(single_match[0])
- continue
-
- # more sufisticate approach - we have to check the possibilities
-
- # TWO TILEID COVERING LAEA GRID TOGETHER
- if s2_aoi.shape[0] == 2:
- # we recheck if the combination of both s2 tiles completly cover the LEAE grid tile
- double_aoi = unary_union(s2_aoi.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_aoi.tile_id.unique().tolist())])
- continue
- else:
- print(f'--- error: double S2 tile matches do not cover LAEA grid {row.name} completly')
- results.append([row.name, 'error - double result but no full coverage even when combined',
- None, ",".join(str(element) for element in s2_aoi.tile_id.unique().tolist())])
- continue
-
- # THREE S2 TILEID INTERSECTION DECISION TREE
- if s2_aoi.shape[0] == 3:
- # do we have multiple epsg zones in the result
- epsg_test = s2_aoi.epsg.unique().tolist()
-
- # if we have TWO epsg zones we can check if the double of the one epsg zone covers all
- if len(epsg_test) == 2:
- s2_sub_a = s2_aoi[s2_aoi.epsg == epsg_test[0]].copy()
- s2_sub_b = s2_aoi[s2_aoi.epsg == epsg_test[1]].copy()
-
- if s2_sub_a.shape[0] == 2:
- s2_sub = s2_sub_a
- else:
- s2_sub = s2_sub_b
-
- # now we have only the two tileIDs with the same epsg number
- double_aoi = unary_union(s2_sub.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- continue
-
- # OK, we do not have two epsg zones OR we have two epsg zones and the double of the one zone was not enough
- # we check now all double combination of the triplet - first combi wins
- lS2tiles = s2_aoi.tile_id.unique().tolist()
- lCombi = list(itertools.combinations(lS2tiles, 2))
-
- iDone = False
- for element in lCombi:
- s2_sub = s2_aoi[s2_aoi['tile_id'].isin(list(element))].copy()
- double_aoi = unary_union(s2_sub.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- iDone = True
- break
- if iDone:
- continue
-
- # damm we really need all three Sentinel-2 tiles to cover the area
- results.append([row.name, 'triplet',
- None,
- ",".join(str(element) for element in s2_aoi.tile_id.unique().tolist())])
- continue
-
- # FOUR S2 TILEID INTERSECTION DECISION TREE
- if s2_aoi.shape[0] == 4:
- # do we have multiple epsg zones in the result
- epsg_test = s2_aoi.epsg.unique().tolist()
-
- # if we have TWO epsg zones we check separatly
- # we can have two or three tileids in one epsg zone max
- if len(epsg_test) == 2:
- # all tests for the first epsg zone
- s2_sub_a = s2_aoi[s2_aoi.epsg == epsg_test[0]].copy()
-
- # now we have to check if we have only 2 tiles in the epsg filter or three
- if s2_sub_a.shape[0] == 2:
-
- # now we have only the two tileIDs with the same epsg number
- double_aoi = unary_union(s2_sub_a.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_sub_a.tile_id.unique().tolist())])
- continue
-
- if s2_sub_a.shape[0] == 3:
- # now we check all double combinations of the three results
- lS2tiles = s2_sub_a.tile_id.unique().tolist()
- lCombi = list(itertools.combinations(lS2tiles, 2))
-
- iDone = False
- for element in lCombi:
- s2_sub = s2_aoi[s2_aoi['tile_id'].isin(list(element))].copy()
- double_aoi = unary_union(s2_sub.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- iDone = True
- break
- if iDone:
- continue
-
- # we check if the triplet can cover the areas
- tripple_aoi = unary_union(s2_sub_a.geometry.tolist())
- if tripple_aoi.contains(row.geometry):
- print(f' - tripple winner ')
- results.append([row.name, 'triplet',
- None,
- ",".join(str(element) for element in s2_sub_a.tile_id.unique().tolist())])
- continue
-
- # all tests for the second zone
- s2_sub_b = s2_aoi[s2_aoi.epsg == epsg_test[1]].copy()
-
- # now we have to check if we have only 2 tiles in the epsg filter or three
- if s2_sub_b.shape[0] == 2:
-
- # now we have only the two tileIDs with the same epsg number
- double_aoi = unary_union(s2_sub_b.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_sub_b.tile_id.unique().tolist())])
- continue
-
- if s2_sub_b.shape[0] == 3:
- # now we check all double combinations of the three results
- lS2tiles = s2_sub_b.tile_id.unique().tolist()
- lCombi = list(itertools.combinations(lS2tiles, 2))
-
- iDone = False
- for element in lCombi:
- s2_sub = s2_aoi[s2_aoi['tile_id'].isin(list(element))].copy()
- double_aoi = unary_union(s2_sub.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- iDone = True
- break
- if iDone:
- continue
-
- # we check if the triplet can cover the areas
- tripple_aoi = unary_union(s2_sub_b.geometry.tolist())
- if tripple_aoi.contains(row.geometry):
- print(f' - tripple winner ')
- results.append([row.name, 'triplet',
- None,
- ",".join(str(element) for element in s2_sub_b.tile_id.unique().tolist())])
- continue
-
- # we have more than two epsg zones OR the tiles within one single epsg zone are not enough
- # we check all double combinations - first wins
- lS2tiles = s2_aoi.tile_id.unique().tolist()
- lCombi = list(itertools.combinations(lS2tiles, 2))
-
- iDone = False
- for element in lCombi:
- s2_sub = s2_aoi[s2_aoi['tile_id'].isin(list(element))].copy()
- double_aoi = unary_union(s2_sub.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- iDone = True
- break
- if iDone:
- continue
-
- # we check all tripple combinations - first wins
- lCombi = list(itertools.combinations(lS2tiles, 3))
-
- iDone = False
- for element in lCombi:
- s2_sub = s2_aoi[s2_aoi['tile_id'].isin(list(element))].copy()
- triple_aoi = unary_union(s2_sub.geometry.tolist())
- if triple_aoi.contains(row.geometry):
- print(f' - triple winner ')
- results.append([row.name, 'triplet',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- iDone = True
- break
- if iDone:
- continue
-
- # damm we really need all four S2 tiles - a quadruple
- results.append([row.name, 'quadruple',
- None,
- ",".join(str(element) for element in s2_aoi.tile_id.unique().tolist())])
- continue
-
- # now we deal with the rest - we should have not more than quadruples per epsg zone
-
- # do we have multiple epsg zones in the result
- epsg_test = s2_aoi.epsg.unique().tolist()
-
- # if we have TWO epsg zones we check separatly
- if len(epsg_test) == 2:
- # all tests for the first epsg zone
- s2_sub_a = s2_aoi[s2_aoi.epsg == epsg_test[0]].copy()
-
- # now we have to check if we have only 2 tiles in the epsg filter or three
- if s2_sub_a.shape[0] == 2:
-
- # now we have only the two tileIDs with the same epsg number
- double_aoi = unary_union(s2_sub_a.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_sub_a.tile_id.unique().tolist())])
- continue
-
- if s2_sub_a.shape[0] == 3:
- # we check all double combinations - first wins
- lS2tiles = s2_sub_a.tile_id.unique().tolist()
- lCombi = list(itertools.combinations(lS2tiles, 2))
-
- iDone = False
- for element in lCombi:
- s2_sub = s2_aoi[s2_aoi['tile_id'].isin(list(element))].copy()
- double_aoi = unary_union(s2_sub.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- iDone = True
- break
- if iDone:
- continue
-
- # we check if the triplet can cover the areas
- tripple_aoi = unary_union(s2_sub_a.geometry.tolist())
- if tripple_aoi.contains(row.geometry):
- print(f' - tripple winner ')
- results.append([row.name, 'triplet',
- None,
- ",".join(str(element) for element in s2_sub_a.tile_id.unique().tolist())])
- continue
-
- # check quadruple
- if s2_sub_a.shape[0] == 4:
- # we check all double combinations - first wins
- lS2tiles = s2_sub_a.tile_id.unique().tolist()
- lCombi = list(itertools.combinations(lS2tiles, 2))
-
- iDone = False
- for element in lCombi:
- s2_sub = s2_aoi[s2_aoi['tile_id'].isin(list(element))].copy()
- double_aoi = unary_union(s2_sub.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- iDone = True
- break
- if iDone:
- continue
-
- # we check all tripple combinations - first wins
- lCombi = list(itertools.combinations(lS2tiles, 3))
-
- iDone = False
- for element in lCombi:
- s2_sub = s2_aoi[s2_aoi['tile_id'].isin(list(element))].copy()
- triple_aoi = unary_union(s2_sub.geometry.tolist())
- if triple_aoi.contains(row.geometry):
- print(f' - triple winner ')
- results.append([row.name, 'triplet',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- iDone = True
- break
- if iDone:
- continue
-
- quad_aoi = unary_union(s2_sub_a.geometry.tolist())
- if quad_aoi.contains(row.geometry):
- print(f' - quadruple winner ')
- results.append([row.name, 'quadruple',
- None,
- ",".join(str(element) for element in s2_sub_a.tile_id.unique().tolist())])
- continue
-
- # all tests for the second epsg zone
- s2_sub_b = s2_aoi[s2_aoi.epsg == epsg_test[1]].copy()
-
- # now we have to check if we have only 2 tiles in the epsg filter or three
- if s2_sub_b.shape[0] == 2:
-
- # now we have only the two tileIDs with the same epsg number
- double_aoi = unary_union(s2_sub_b.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_sub_b.tile_id.unique().tolist())])
- continue
-
- if s2_sub_b.shape[0] == 3:
- # we check all double combinations - first wins
- lS2tiles = s2_sub_b.tile_id.unique().tolist()
- lCombi = list(itertools.combinations(lS2tiles, 2))
-
- iDone = False
- for element in lCombi:
- s2_sub = s2_aoi[s2_aoi['tile_id'].isin(list(element))].copy()
- double_aoi = unary_union(s2_sub.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- iDone = True
- break
- if iDone:
- continue
-
- # we check if the triplet can cover the areas
- tripple_aoi = unary_union(s2_sub_b.geometry.tolist())
- if tripple_aoi.contains(row.geometry):
- print(f' - tripple winner ')
- results.append([row.name, 'triplet',
- None,
- ",".join(str(element) for element in s2_sub_b.tile_id.unique().tolist())])
- continue
-
- # check quadruple
- if s2_sub_b.shape[0] == 4:
- # we check all double combinations - first wins
- lS2tiles = s2_sub_b.tile_id.unique().tolist()
- lCombi = list(itertools.combinations(lS2tiles, 2))
-
- iDone = False
- for element in lCombi:
- s2_sub = s2_aoi[s2_aoi['tile_id'].isin(list(element))].copy()
- double_aoi = unary_union(s2_sub.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- iDone = True
- break
- if iDone:
- continue
-
- # we check all tripple combinations - first wins
- lCombi = list(itertools.combinations(lS2tiles, 3))
-
- iDone = False
- for element in lCombi:
- s2_sub = s2_aoi[s2_aoi['tile_id'].isin(list(element))].copy()
- triple_aoi = unary_union(s2_sub.geometry.tolist())
- if triple_aoi.contains(row.geometry):
- print(f' - triple winner ')
- results.append([row.name, 'triplet',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- iDone = True
- break
- if iDone:
- continue
-
- quad_aoi = unary_union(s2_sub_b.geometry.tolist())
- if quad_aoi.contains(row.geometry):
- print(f' - quadruple winner ')
- results.append([row.name, 'quadruple',
- None,
- ",".join(str(element) for element in s2_sub_b.tile_id.unique().tolist())])
- continue
-
-
- # we check all double combinations - first wins
- lS2tiles = s2_aoi.tile_id.unique().tolist()
- lCombi = list(itertools.combinations(lS2tiles, 2))
-
- iDone = False
- for element in lCombi:
- s2_sub = s2_aoi[s2_aoi['tile_id'].isin(list(element))].copy()
- double_aoi = unary_union(s2_sub.geometry.tolist())
- if double_aoi.contains(row.geometry):
- print(f' - double winner ')
- results.append([row.name, 'double',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- iDone = True
- break
- if iDone:
- continue
-
- # we check all tripple combinations - first wins
- lCombi = list(itertools.combinations(lS2tiles, 3))
-
- iDone = False
- for element in lCombi:
- s2_sub = s2_aoi[s2_aoi['tile_id'].isin(list(element))].copy()
- triple_aoi = unary_union(s2_sub.geometry.tolist())
- if triple_aoi.contains(row.geometry):
- print(f' - triple winner ')
- results.append([row.name, 'triplet',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- iDone = True
- break
- if iDone:
- continue
-
- # we check all quadruple combinations - first wins
- lCombi = list(itertools.combinations(lS2tiles, 4))
-
- iDone = False
- for element in lCombi:
- s2_sub = s2_aoi[s2_aoi['tile_id'].isin(list(element))].copy()
- triple_aoi = unary_union(s2_sub.geometry.tolist())
- if triple_aoi.contains(row.geometry):
- print(f' - quadruple winner ')
- results.append([row.name, 'quadruple',
- None,
- ",".join(str(element) for element in s2_sub.tile_id.unique().tolist())])
- iDone = True
- break
- if iDone:
- continue
-
- # fuck - still no winner..... I already tested too much and made it way to complicated.... just put all in
- print(f' - multiple winner - use all tileid')
- results.append([row.name, 'multiple - write filter by hand',
- None, ",".join(str(element) for element in s2_aoi.tile_id.unique().tolist())])
-
-# merge the results and write out
-df_result = pd.DataFrame(results, columns=['laea_tileid', 'match', 's2_tileid', 's2_multi_list'])
-
-if gdf_laea.shape[0] != df_result.shape[0]:
- print(' -- error: we have a mismatch between number of LAEA grids and the number of results')
-
-gdf_result = gdf_laea.merge(df_result, how='left', left_on='name', right_on='laea_tileid')
-
-#write out
-gdf_result.to_file(path_out)
\ No newline at end of file
diff --git a/scripts/laea_tiling_grids_4_finalize.py b/scripts/laea_tiling_grids_4_finalize.py
deleted file mode 100644
index 99d31a7..0000000
--- a/scripts/laea_tiling_grids_4_finalize.py
+++ /dev/null
@@ -1,91 +0,0 @@
-"""
-finalizing the LAEA grids
-- the polygons itself have to be in 4326 for better intersection with geojson or geoparquet files
-- name as str
-- bbox_dict as str BUT the crs has to be an int
-- s2_tileid_list as str seperated by comma
-
-we first run 20k, then 50K and in the end 100k (for 100K the tile list is assembled from the 20K)
-
-"""
-
-try:
- import importlib.resources as importlib_resources
-except:
- import importlib_resources
-import eo_processing.resources
-from os.path import normpath
-import geopandas as gpd
-from eo_processing.utils.helper import convert_to_list
-from eo_processing.utils.geoprocessing import laea100km_id_to_extent
-
-### start with 20K
-path_20k = r'C:\Users\buchhorm\Downloads\new_grids\LAEA_20km_tiling_grid_EU_EPSG4326.gpkg'
-path_20k_meta = r'C:\Users\buchhorm\Downloads\new_grids\S2_tile_info_20K_grid.gpkg'
-
-gdf20 = gpd.read_file(path_20k)
-gdf20_meta = gpd.read_file(path_20k_meta)
-
-out_20k = r'C:\Users\buchhorm\Downloads\new_grids\LAEA-20km_add-info.gpkg'
-
-gdf20 = gdf20.merge(gdf20_meta[['name', 's2_multi_list']], how='left', left_on='name', right_on='name')
-gdf20.rename(columns={'s2_multi_list': 's2_tileid_list'}, inplace=True)
-gdf20['bbox_dict'] = gdf20.bbox_dict.str.replace("'EPSG:3035'", "3035")
-
-gdf20[['name', 's2_tileid_list', 'bbox_dict', 'geometry']].to_file(out_20k, driver='GPKG')
-
-
-### now the 50K
-path_50k = r'C:\Users\buchhorm\Downloads\new_grids\LAEA_50km_tiling_grid_EU_EPSG4326.gpkg'
-path_50k_meta = r'C:\Users\buchhorm\Downloads\new_grids\S2_tile_info_50K_grid.gpkg'
-
-gdf50 = gpd.read_file(path_50k)
-gdf50_meta = gpd.read_file(path_50k_meta)
-
-out_50k = r'C:\Users\buchhorm\Downloads\new_grids\LAEA-50km_add-info.gpkg'
-
-gdf50 = gdf50.merge(gdf50_meta[['name', 's2_multi_list']], how='left', left_on='name', right_on='name')
-gdf50.rename(columns={'s2_multi_list': 's2_tileid_list'}, inplace=True)
-gdf50['bbox_dict'] = gdf50.bbox_dict.str.replace("'EPSG:3035'", "3035")
-
-gdf50[['name', 's2_tileid_list', 'bbox_dict', 'geometry']].to_file(out_50k, driver='GPKG')
-
-### finally the 100K
-path_100k = r'C:\Users\buchhorm\Downloads\new_grids\LAEA_100km_tiling_grid_EU_EPSG4326.gpkg'
-out_100k = r'C:\Users\buchhorm\Downloads\new_grids\LAEA-100km_add-info.gpkg'
-
-#create the metadata from the 20k
-gdf100 = gpd.read_file(path_100k)
-gdf100_meta = gdf20.copy()
-
-gdf100_meta['name100'] = gdf100_meta.name.apply(lambda x: x[:3] + x[4:7])
-# convert s2_tileid_list to list
-gdf100_meta['s2_tileid_list'] = gdf100_meta.apply(lambda row : row['s2_tileid_list'].split(','), axis =1)
-
-# update the 's2_tileid_list' column
-s2_grouped = gdf100_meta.groupby('name100')['s2_tileid_list'].apply(
- lambda x: list(set([item for sublist in x for item in sublist]))
- )
-gdf100_meta['s2_tileid_list100'] = gdf100_meta['name100'].map(s2_grouped)
-
-# now we dissolve by name100
-gdf100_meta = gdf100_meta.dissolve(by='name100')
-
-gdf100_meta.reset_index(inplace=True)
-
-gdf100_meta.drop(columns=['name', 's2_tileid_list'], inplace=True)
-
-gdf100_meta.rename(columns={'name100': 'name', 's2_tileid_list100': 's2_tileid_list'}, inplace=True)
-
-#convert list to string
-gdf100_meta['s2_tileid_list'] = gdf100_meta['s2_tileid_list'].apply(lambda x: str(x))
-gdf100_meta['s2_tileid_list'] = gdf100_meta['s2_tileid_list'].apply(lambda x: x.replace("'", ""))
-gdf100_meta['s2_tileid_list'] = gdf100_meta['s2_tileid_list'].apply(lambda x: x.replace(" ", ""))
-gdf100_meta['s2_tileid_list'] = gdf100_meta['s2_tileid_list'].apply(lambda x: x.replace("[", ""))
-gdf100_meta['s2_tileid_list'] = gdf100_meta['s2_tileid_list'].apply(lambda x: x.replace("]", ""))
-
-#add this info to the 100k grid
-gdf100 = gdf100.merge(gdf100_meta[['name', 's2_tileid_list']], how='left', left_on='name', right_on='name')
-gdf100['bbox_dict'] = gdf100.bbox_dict.str.replace("'EPSG:3035'", "3035")
-
-gdf100[['name', 's2_tileid_list', 'bbox_dict', 'geometry']].to_file(out_100k, driver='GPKG')
diff --git a/scripts/test_catalog_check.ipynb b/scripts/test_catalog_check.ipynb
deleted file mode 100644
index dbdc5fc..0000000
--- a/scripts/test_catalog_check.ipynb
+++ /dev/null
@@ -1,286 +0,0 @@
-{
- "cells": [
- {
- "metadata": {},
- "cell_type": "markdown",
- "source": "# improved catalog check for Sentinel-1 and Sentinel-2",
- "id": "a5ed738ce4acf5d0"
- },
- {
- "metadata": {},
- "cell_type": "markdown",
- "source": "#### old way",
- "id": "d1e5dde20aee7b70"
- },
- {
- "cell_type": "code",
- "id": "initial_id",
- "metadata": {
- "collapsed": true,
- "ExecuteTime": {
- "end_time": "2025-11-27T16:11:54.983Z",
- "start_time": "2025-11-27T16:11:35.997053Z"
- }
- },
- "source": [
- "from __future__ import annotations\n",
- "import requests\n",
- "import json\n",
- "from eo_processing.utils.geoprocessing import reproj_bbox_to_ll\n",
- "\n",
- "# determine S1 or S2\n",
- "sentinel = 'S2'\n",
- "start = '2020-01-01T00:00:00.00Z'\n",
- "end = '2024-12-31T23:59:59.00Z'\n",
- "aoi = {'east': 4840000, 'south': 2800000, 'west': 4820000, 'north': 2820000, 'crs': 3035}\n",
- "latlon_box = reproj_bbox_to_ll(aoi)\n",
- "\n",
- "\n",
- "#-------\n",
- "if sentinel == 'S1': satelite = \"SENTINEL-1\"\n",
- "elif sentinel == 'S2': satelite = \"SENTINEL-2\"\n",
- "else: raise ValueError(f\"{sentinel} is not satellite for which this has been implemented\")\n",
- "\n",
- "url= (f\"https://datahub.creodias.eu/odata/v1/Products?$filter=Collection/Name eq '{satelite}' and \"\n",
- " f\"OData.CSC.Intersects(area=geography'SRID=4326;{latlon_box}') and ContentDate/Start \"\n",
- " f\"gt {start} and ContentDate/Start lt {end}&$top={100}\")\n",
- "results = requests.get(url)\n",
- "json_data = json.loads(results.text)\n",
- "print(len(json_data[\"value\"]))"
- ],
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "100\n"
- ]
- }
- ],
- "execution_count": 1
- },
- {
- "metadata": {},
- "cell_type": "markdown",
- "source": "#### improved version based on pySTAC",
- "id": "b9502fa7044e16b8"
- },
- {
- "metadata": {
- "ExecuteTime": {
- "end_time": "2025-11-27T16:14:45.238415Z",
- "start_time": "2025-11-27T16:14:45.228897Z"
- }
- },
- "cell_type": "code",
- "source": [
- "from typing import TYPE_CHECKING\n",
- "from eo_processing.utils.geoprocessing import reproj_bbox_to_ll\n",
- "import pystac_client\n",
- "import pandas as pd\n",
- "from eo_processing.config.data_formats import openEO_bbox_format\n",
- "\n",
- "\n",
- "def catalog_check_CDSE_S2(start: str, end: str, bbox: openEO_bbox_format) -> None:\n",
- " #quickfix on dates that are in date format\n",
- " if not 'Z' in start:\n",
- " start = start + \"T00:00:00.00Z\"\n",
- " if not 'Z' in end:\n",
- " end = end + \"T00:00:00.00Z\"\n",
- "\n",
- " # set the minimum number of S2 images with two satellites (5 daily observation)\n",
- " MIN_VALUE_S2 = 1./5.\n",
- " #in 2017 S2B started in june/july so than only S2A sattelite\n",
- " if pd.to_datetime(start).year == '2017':\n",
- " MIN_VALUE_S2 = 1./10.\n",
- "\n",
- " # the percentage of observations we want to have at least\n",
- " percentage = 0.8\n",
- " # convert the openEO bbox format to a shapely Polygon\n",
- " latlon_box = reproj_bbox_to_ll(bbox)\n",
- " # number of days in the temporal extent\n",
- " temp_extent_days = (pd.to_datetime(end)-pd.to_datetime(start)).days\n",
- "\n",
- " # run the PySTAC-client search\n",
- " # Connect to the Copernicus Data Space Ecosystem STAC API\n",
- " #catalog_url = \"https://catalogue.dataspace.copernicus.eu/stac\"\n",
- " catalog_url = \"https://stac.dataspace.copernicus.eu/v1/\"\n",
- " client = pystac_client.Client.open(catalog_url)\n",
- "\n",
- " search = client.search(\n",
- " collections=['sentinel-2-l2a'],\n",
- " bbox=list(latlon_box.bounds),\n",
- " datetime=f\"{start}/{end}\",\n",
- " fields=[\"id\", \"properties.datetime\"],\n",
- " #query={\"eo:cloud_cover\": {\"lt\": 95}},\n",
- " )\n",
- "\n",
- " # get the dates of all found matches\n",
- " results = []\n",
- " for item in search.items_as_dicts():\n",
- " results.append(item['properties']['datetime'])\n",
- "\n",
- " # count the number of unique dates on which we have observations (resolved tile overlap)\n",
- " df = pd.DataFrame(results, columns=['date'])\n",
- " df['date'] = pd.to_datetime(df['date'])\n",
- " df['date'] = df['date'].apply(lambda x: x.date())\n",
- " nbr_files = df['date'].nunique()\n",
- "\n",
- " print(f'Found {nbr_files} images.')\n",
- "\n",
- " # run the test\n",
- " if nbr_files < MIN_VALUE_S2*percentage*temp_extent_days:\n",
- " raise ValueError(f'not enough S2 images. Found {nbr_files} images.')"
- ],
- "id": "f49831019d51f465",
- "outputs": [],
- "execution_count": 6
- },
- {
- "metadata": {
- "ExecuteTime": {
- "end_time": "2025-11-27T16:17:49.103914Z",
- "start_time": "2025-11-27T16:15:26.656716Z"
- }
- },
- "cell_type": "code",
- "source": "catalog_check_CDSE_S2(start, end, aoi)",
- "id": "a361f2ce20365c6",
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Found 724 images.\n"
- ]
- }
- ],
- "execution_count": 8
- },
- {
- "metadata": {},
- "cell_type": "markdown",
- "source": "#### now Sentinel-1",
- "id": "33aef142792d22e0"
- },
- {
- "metadata": {},
- "cell_type": "code",
- "source": [
- "def catalogue_check_CDSE_S1(orbit_direction: str, start: str, end: str, bbox: openEO_bbox_format) -> str | None:\n",
- " #quickfix on dates that are in date format\n",
- " if not 'Z' in start:\n",
- " start = start + \"T00:00:00.00Z\"\n",
- " if not 'Z' in end:\n",
- " end = end + \"T00:00:00.00Z\"\n",
- "\n",
- " # set the minimum number of S1 images with two satellites\n",
- " MIN_VALUE_S1 = 1./12.\n",
- " # the percentage of observations we want to have at least\n",
- " percentage = 0.8\n",
- " # convert the openEO bbox format to a shapely Polygon\n",
- " latlon_box = reproj_bbox_to_ll(bbox)\n",
- " # number of days in the temporal extent\n",
- " temp_extent_days = (pd.to_datetime(end)-pd.to_datetime(start)).days\n",
- "\n",
- " # run the PySTAC-client search\n",
- " # Connect to the Copernicus Data Space Ecosystem STAC API\n",
- " #catalog_url = \"https://catalogue.dataspace.copernicus.eu/stac\"\n",
- " catalog_url = \"https://stac.dataspace.copernicus.eu/v1/\"\n",
- " client = pystac_client.Client.open(catalog_url)\n",
- "\n",
- " # if we have an orbit_direction given we have to test that first\n",
- " if orbit_direction is not None:\n",
- " if orbit_direction not in ['ASCENDING', 'DESCENDING']:\n",
- " raise ValueError(\n",
- " f'`orbit_direction` value `{orbit_direction}` not recognized.')\n",
- "\n",
- " search = client.search(\n",
- " collections=['sentinel-1-grd'],\n",
- " bbox=list(latlon_box.bounds),\n",
- " datetime=f\"{start}/{end}\",\n",
- " query={\"sat:orbit_state\": {\"eq\": f\"{orbit_direction.lower()}\"},\n",
- " #\"sar:polarizations\": {\"eq\": \"VV&VH\"},\n",
- " },\n",
- " )\n",
- "\n",
- " # get the dates of all found matches\n",
- " results = []\n",
- " for item in search.items():\n",
- " results.append(item.datetime.date())\n",
- "\n",
- " # count the number of unique dates on which we have observations (resolved tile overlap)\n",
- " df = pd.DataFrame(results, columns=['date'])\n",
- " nbr_files = df['date'].nunique()\n",
- "\n",
- " if nbr_files < MIN_VALUE_S1*percentage*temp_extent_days:\n",
- " print(f'Not enough S1 images with orbit {orbit_direction}. \\n' + \\\n",
- " f'Found {nbr_files} images.')\n",
- " else: return orbit_direction\n",
- " #use both orbits -> check with both directions.\n",
- "\n",
- " search = client.search(\n",
- " collections=['sentinel-1-grd'],\n",
- " bbox=list(latlon_box.bounds),\n",
- " datetime=f\"{start}/{end}\",\n",
- " #query={\"sar:polarizations\": {\"eq\": \"VV&VH\"} },\n",
- " )\n",
- "\n",
- " # get the dates of all found matches\n",
- " results = []\n",
- " for item in search.items():\n",
- " results.append(item.datetime.date())\n",
- "\n",
- " # count the number of unique dates on which we have observations (resolved tile overlap)\n",
- " df = pd.DataFrame(results, columns=['date'])\n",
- " nbr_files = df['date'].nunique()\n",
- "\n",
- " if nbr_files < MIN_VALUE_S1*percentage*temp_extent_days:\n",
- " raise ValueError(f'not enough S1 without orbit direction selection. \\n'+ \\\n",
- " f'Found {nbr_files} images.')\n",
- "\n",
- " return None"
- ],
- "id": "bb7d1e0aae2045c2",
- "outputs": [],
- "execution_count": null
- },
- {
- "metadata": {},
- "cell_type": "code",
- "source": "orbit = catalogue_check_CDSE_S1('DESCENDING', start, '2020-12-31T23:59:59.00Z', aoi)",
- "id": "bb7870f8b4c24364",
- "outputs": [],
- "execution_count": null
- },
- {
- "metadata": {},
- "cell_type": "code",
- "source": "",
- "id": "39c717f2367d642b",
- "outputs": [],
- "execution_count": null
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Python 3",
- "language": "python",
- "name": "python3"
- },
- "language_info": {
- "codemirror_mode": {
- "name": "ipython",
- "version": 2
- },
- "file_extension": ".py",
- "mimetype": "text/x-python",
- "name": "python",
- "nbconvert_exporter": "python",
- "pygments_lexer": "ipython2",
- "version": "2.7.6"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/src/eo_processing/_version.py b/src/eo_processing/_version.py
index 35b620f..3ea5268 100644
--- a/src/eo_processing/_version.py
+++ b/src/eo_processing/_version.py
@@ -1,3 +1,3 @@
#!/usr/bin/env python3
-__version__ = '0.7.0'
+__version__ = '0.8.0'
diff --git a/src/eo_processing/config/settings.py b/src/eo_processing/config/settings.py
index 12f11a1..cc7de74 100644
--- a/src/eo_processing/config/settings.py
+++ b/src/eo_processing/config/settings.py
@@ -93,27 +93,29 @@
}
OPENEO_EXTRACT_CDSE_JOB_OPTIONS: Dict[str, Union[str, int, float, List]] = {
- "driver-memory": "8G",
- "driver-memoryOverhead": "5G",
+ "driver-memory": "4G",
+ "driver-memoryOverhead": "2G",
"driver-cores": 1,
- "executor-memory": "2000m",
- "executor-memoryOverhead": "256m",
- "python-memory": "2500m",
+ "executor-memory": "2G",
+ "executor-memoryOverhead": "1G",
+ "python-memory": "disable",
"executor-cores": 1,
- "max-executors": 25,
+ "max-executors": 20,
"logging-threshold": "info"
}
OPENEO_INFERENCE_CDSE_JOB_OPTIONS: Dict[str, Union[str, int, float, List]] = {
- "driver-memory": "1000m",
- "driver-memoryOverhead": "1000m",
+ "driver-memory": "5G",
+ "driver-memoryOverhead": "4G",
"driver-cores": 1,
- "executor-memory": "1500m",
- "executor-memoryOverhead": "256m",
+ "executor-memory": "5G",
+ "executor-memoryOverhead": "5G",
"executor-cores": 1,
- "max-executors": 20,
- "python-memory": "4000m",
+ "max-executors": 15,
+ "python-memory": "disable",
"logging-threshold": "info",
+ "force-s3proxy": True,
+ "soft-errors": 0.05,
"udf-dependency-archives": [
"https://s3.waw3-1.cloudferro.com/project_dependencies/onnx_deps_python311.zip#onnx_deps"
]
@@ -124,11 +126,12 @@
"driver-memoryOverhead": "2G",
"driver-cores": 1,
"executor-memory": "2G",
- "executor-memoryOverhead": "2G",
+ "executor-memoryOverhead": "3G",
"python-memory": "disable",
"executor-cores": 1,
- "max-executors": 25,
- "logging-threshold": "info"
+ "max-executors": 15,
+ "logging-threshold": "info",
+ "soft-errors": 0.05
}
OPENEO_CUBEEXTRACTION_CDSE_JOB_OPTIONS: Dict[str, Union[str, int, float, List]] = {
diff --git a/src/eo_processing/resources/udf_force_into_proba.py b/src/eo_processing/resources/udf_force_into_proba.py
new file mode 100644
index 0000000..6e5d6bc
--- /dev/null
+++ b/src/eo_processing/resources/udf_force_into_proba.py
@@ -0,0 +1,228 @@
+import xarray as xr
+from typing import Dict
+from openeo.udf import inspect
+from openeo.metadata import CubeMetadata
+
+def apply_metadata(metadata: CubeMetadata, context:Dict) -> CubeMetadata:
+ """ Rename the bands by using openeo apply_metadata function
+ :param metadata: Metadata of the input data cube
+ :param context: Context of the UDF
+ :return: renamed labels
+ """
+ band_names = context.get("band_names")
+ typology_schema = context.get("typology_schema")
+ force_ocean = context.get("force_ocean")
+ force_snow = context.get("force_snow")
+ mask_non_mangrove = context.get("mask_non_mangrove")
+
+ if force_ocean:
+ if typology_schema == "EUNIS2021plus":
+ needed = "Level1_class-0_habitat-MH-20000"
+ elif typology_schema == "IUCNGET":
+ needed = "Level1_class-0_habitat-M-80000"
+ else:
+ raise Exception("unknown typology schema - please adjust the code base of UDF")
+ if needed not in band_names:
+ band_names.append(needed)
+
+ if force_snow:
+ if typology_schema == "EUNIS2021plus":
+ needed = "Level1_class-0_habitat-U-90000"
+ elif typology_schema == "IUCNGET":
+ needed = "Level1_class-0_habitat-T-10000"
+ else:
+ raise Exception("unknown typology schema - please adjust the code base of UDF")
+ if needed not in band_names:
+ band_names.append(needed)
+
+ if mask_non_mangrove:
+ if typology_schema == "IUCNGET":
+ needed = "Level1_class-0_habitat-MFT-100000"
+ else:
+ raise Exception("currently mangrove masking is only possible in IUCNGET. set parameter to False for "
+ "EUNIS or adjust the code base of UDF for new typology.")
+ if needed not in band_names:
+ band_names.append(needed)
+
+ return metadata.rename_labels(dimension="bands", target=band_names)
+
+def apply_datacube(cube: xr.DataArray, context: Dict) -> xr.DataArray:
+ """ imprint or mask external data into the probability cube
+
+ :param cube: data cube
+ :param context: dictionary to provide external data - key class_mapping is used to provide external re-mapping dict
+ :return: data cube with remapped values (new cube)
+ """
+
+ # get all needed data together
+ output_band_names = context.get("band_names")
+ typology_schema = context.get("typology_schema")
+ force_ocean = context.get("force_ocean")
+ force_snow = context.get("force_snow")
+ mask_non_mangrove = context.get("mask_non_mangrove")
+ mask_non_ocean = context.get("mask_non_ocean")
+ inspect(message=f"settings: typology_schema: {typology_schema}, force_ocean: {force_ocean}, force_snow: "
+ f"{force_snow}, mask_non_mangrove: {mask_non_mangrove}, mask_non_ocean: {mask_non_ocean},"
+ f"output_band_names ({len(output_band_names)}): {output_band_names}")
+
+ # get the band names of cube handed over to UDF
+ input_band_names = cube.indexes["bands"].values
+ inspect(message=f"input cube band names ({len(input_band_names)}): {input_band_names}")
+
+ # check that we have to imprint "ocean"
+ if force_ocean:
+ inspect(message="imprinting ocean")
+ # check that we have the needed band
+ if "ocean" not in input_band_names:
+ raise ValueError("no ocean band in cube - please add a ocean band to the input cube")
+
+ cube_ocean = cube.sel(bands="ocean")
+ ocean_value = 255
+
+ if typology_schema == "EUNIS2021plus":
+ imprint_band = "Level1_class-0_habitat-MH-20000"
+ elif typology_schema == "IUCNGET":
+ imprint_band = "Level1_class-0_habitat-M-80000"
+ else:
+ raise Exception("unknown typology schema - please adjust the code base of UDF")
+
+ # check that we have the needed output band, if not create
+ if imprint_band not in input_band_names:
+ inspect(message=f"+ create extra proba layer for ocean ({imprint_band})")
+ output_band_names.append(imprint_band)
+ # Create zero-filled band with same spatial dimensions as cube
+ zero_band = xr.zeros_like(cube.isel(bands=0))
+ zero_band = zero_band.assign_coords(bands=imprint_band)
+ # Add the new band to the cube
+ cube = xr.concat([cube, zero_band], dim="bands")
+
+ # now we imprint the ocean in band "Level1_class-0_habitat-M-80000"
+ ocean_mask = cube_ocean == ocean_value
+ inspect(message=f"+ imprint ocean mask into band {imprint_band}")
+ cube.loc[dict(bands=imprint_band)] = xr.where(ocean_mask, 100, cube.loc[dict(bands=imprint_band)])
+
+ # reset the values of all other bands on level 1 to ZERO
+ other_bands_level1 = [x for x in input_band_names if x.startswith("Level1_class-0")]
+ if imprint_band in other_bands_level1:
+ other_bands_level1.remove(imprint_band)
+ inspect(message=f"+ imprint ocean mask into non-marine bands ({other_bands_level1})")
+ cube.loc[dict(bands=other_bands_level1)] = xr.where(ocean_mask, 0, cube.loc[dict(bands=other_bands_level1)])
+
+ # now mask areas wrongly classified as ocean (e.g. in the case of the ocean mask)
+ # NOTE: currently only needed for IUCN-GET typology
+ if mask_non_ocean and (typology_schema == "IUCNGET"):
+ inspect(message=f"+ reset areas in the {imprint_band} band which are not really ocean. "
+ f"(level1 second winner will be used)")
+ # get the mask of the non-ocean areas
+ non_ocean_mask = ~ocean_mask
+ # now mask the non-ocean areas
+ cube.loc[dict(bands=imprint_band)] = xr.where(non_ocean_mask, 0, cube.loc[dict(bands=imprint_band)])
+
+ # check that we have to imprint "snow"
+ if force_snow:
+ inspect(message="imprinting snow")
+ # check that we have the needed band
+ if "snow" not in input_band_names:
+ raise ValueError("no snow band in cube - please add a snow band to the input cube")
+
+ snow_cube = cube.sel(bands="snow")
+ snow_value = 110 # value in LCFM 10m maps for permanent snow & ice
+
+ if typology_schema == "EUNIS2021plus":
+ imprint_band = "Level1_class-0_habitat-U-90000"
+ elif typology_schema == "IUCNGET":
+ imprint_band = "Level1_class-0_habitat-T-10000"
+ else:
+ raise Exception("unknown typology schema - please adjust the code base of UDF")
+
+ # check that we have the needed output band, if not create
+ if imprint_band not in input_band_names:
+ inspect(message=f"+ create extra proba layer for snow ({imprint_band})")
+ output_band_names.append(imprint_band)
+ # Create zero-filled band with same spatial dimensions as cube
+ zero_band = xr.zeros_like(cube.isel(bands=0))
+ zero_band = zero_band.assign_coords(bands=imprint_band)
+ # Add the new band to the cube
+ cube = xr.concat([cube, zero_band], dim="bands")
+
+ # create mask and imprint snow into band
+ snow_mask = snow_cube == snow_value
+ inspect(message=f"+ imprint snow mask into band {imprint_band}")
+ cube.loc[dict(bands=imprint_band)] = xr.where(snow_mask, 100, cube.loc[dict(bands=imprint_band)])
+
+ # reset the values of all other bands on level 1 to ZERO
+ other_bands_level1 = [x for x in input_band_names if x.startswith("Level1_class-0")]
+ if imprint_band in other_bands_level1:
+ other_bands_level1.remove(imprint_band)
+ inspect(message=f"+ imprint snow mask into non-snow bands ({other_bands_level1})")
+ cube.loc[dict(bands=other_bands_level1)] = xr.where(snow_mask, 0, cube.loc[dict(bands=other_bands_level1)])
+
+ # check that we have to mask mangrove
+ if mask_non_mangrove:
+ inspect(message="imprinting mangroves")
+ # check that we have the needed band
+ if "mangrove" not in input_band_names:
+ raise ValueError("no mangrove band in cube - please add a mangrove band to the input cube")
+
+ mango_cube = cube.sel(bands="mangrove")
+ mango_value = 1 # raster value for potential areas of mangrove
+
+ if typology_schema == "IUCNGET":
+ imprint_band = "Level1_class-0_habitat-MFT-100000"
+ else:
+ raise Exception("currently mangrove masking is only possible in IUCNGET. set parameter to False for "
+ "EUNIS or adjust the code base of UDF for new typology.")
+
+ # check that we have the needed output band, if not create
+ if imprint_band not in input_band_names:
+ inspect(message=f"+ create extra proba layer for mangrove ({imprint_band})")
+ output_band_names.append(imprint_band)
+ # Create zero-filled band with same spatial dimensions as cube
+ zero_band = xr.zeros_like(cube.isel(bands=0))
+ zero_band = zero_band.assign_coords(bands=imprint_band)
+ # Add the new band to the cube
+ cube = xr.concat([cube, zero_band], dim="bands")
+
+ # We have to first determine the level3 winner in the level1 MFT results. When the winner is mangrove (MFT1.2) and
+ # outside the mangrove potential area mask THEN we can set the level1 pixel to ZERO. (Again under the assumption that
+ # these areas can not be MFT1.1 or MFT1.3)
+ # since level2 is a single class model we do not have to check it
+
+ level3_bands = ['Level3_class-MFT1_habitat-MFT1.1-100101',
+ 'Level3_class-MFT1_habitat-MFT1.2-100102',
+ 'Level3_class-MFT1_habitat-MFT1.3-100103']
+
+ if all(band in input_band_names for band in level3_bands):
+ # get the level3 winner
+ level3_data = cube.sel(bands=level3_bands).copy()
+ level3_data = level3_data.fillna(0)
+ level3_winner = level3_data.sel(bands=level3_bands).argmax('bands')
+ # get mask of MFT1.2 winner
+ level3_mangrove_winner = level3_winner == 1
+ # check if value is bigger than 0
+ level3_mangrove_winner_value = level3_data.sel(bands='Level3_class-MFT1_habitat-MFT1.2-100102') > 0
+
+ # create mask where mangrove can exist
+ mango_mask = mango_cube == mango_value
+
+ # create the removal mask -> MFT1.2 outside the mangrove potential area mask
+ removal_mask = ~mango_mask & level3_mangrove_winner & level3_mangrove_winner_value
+
+ n_pixels_to_remove = int(removal_mask.sum().values)
+
+ # reset the level1 MFT values to ZERO in these areas
+ inspect(
+ message=f"+ reset areas in the {imprint_band} band were wrongly mangroves were detected ({n_pixels_to_remove} pixel)"
+ f"(level1 second winner will be chosen)")
+ cube.loc[dict(bands=imprint_band)] = xr.where(removal_mask, 0, cube.loc[dict(bands=imprint_band)])
+ else:
+ inspect(message=f"+ no level3 bands found in cube - cannot mask mangroves")
+
+ # filter the cube to only the bands we need
+ cube = cube.sel(bands=output_band_names)
+ inspect(message=f"output cube band names {len(cube.indexes['bands'].values)}: {cube.indexes['bands'].values}")
+ # final check that we have the right number of bands
+ if len(cube.indexes['bands'].values) != len(output_band_names):
+ raise ValueError("wrong number of bands in output cube - please adjust the code base of UDF")
+
+ return cube
\ No newline at end of file
diff --git a/src/eo_processing/resources/udf_max_occurence_hierarchical_merger.py b/src/eo_processing/resources/udf_max_occurence_hierarchical_merger.py
index 6613fe0..3348772 100644
--- a/src/eo_processing/resources/udf_max_occurence_hierarchical_merger.py
+++ b/src/eo_processing/resources/udf_max_occurence_hierarchical_merger.py
@@ -1,11 +1,9 @@
-import os, sys
import pandas as pd
import numpy as np
import xarray as xr
import re
-from typing import Dict, List, Tuple, Union
+from typing import Dict, List
from openeo.udf import inspect
-from datetime import datetime
from openeo.metadata import CubeMetadata
def apply_metadata(metadata: CubeMetadata, context:Dict) -> CubeMetadata:
@@ -14,7 +12,10 @@ def apply_metadata(metadata: CubeMetadata, context:Dict) -> CubeMetadata:
:param context: Context of the UDF
:return: renamed labels
"""
- return metadata.rename_labels(dimension="bands", target=['EUNIS habitat level3'])
+ band_name = context.get('typology', 'EUNIS')
+ final_band_name = f"{band_name} habitat level3"
+
+ return metadata.rename_labels(dimension="bands", target=[final_band_name])
def _select_highest_prob_class(cube: xr.DataArray, raster_codes) -> xr.DataArray:
""" Select per model the highest probability of occurrence class
@@ -150,8 +151,9 @@ def apply_datacube(cube: xr.DataArray, context:Dict) -> xr.DataArray:
### get the list of classes as output from inference run
# use returned metadata to build up the class dictionary
- inspect(message=cube.indexes["bands"].values)
- df = parse_prob_classes_fromStac(cube.indexes["bands"].values)
+ input_band_names = cube.indexes["bands"].values
+ inspect(message=f"input cube band names ({len(input_band_names)}): {input_band_names}")
+ df = parse_prob_classes_fromStac(input_band_names)
inspect(message=f"## context parameters")
inspect(message=f"{df}")
diff --git a/src/eo_processing/utils/alphaearth_utils.py b/src/eo_processing/utils/alphaearth_utils.py
index cb10ca8..b2cc74c 100644
--- a/src/eo_processing/utils/alphaearth_utils.py
+++ b/src/eo_processing/utils/alphaearth_utils.py
@@ -43,7 +43,8 @@
import subprocess
from pathlib import Path
from typing import List, Optional, Tuple
-
+from concurrent.futures import ProcessPoolExecutor
+from concurrent import futures
import geopandas as gpd
import pandas as pd
from eo_processing.utils.storage import WEED_storage
@@ -129,7 +130,7 @@ def get_bucket_and_key(self, s3_url: str) -> Tuple[str, str]:
split_str = s3_url.replace(S3_BASE_URL, "").split("/", 1)
return split_str[0], split_str[1]
- def download_s3_file(self, s3_url, output_dir):
+ def download_s3_file(self, s3_url, output_dir, overwrite: bool = False):
"""
Download a single file from S3 using unsigned requests.
@@ -149,6 +150,9 @@ def download_s3_file(self, s3_url, output_dir):
try:
output_path = Path(output_dir) / Path(key)
+ if output_path.exists() and not overwrite:
+ logger.info(f"File {output_path} already exists, skipping download.")
+ return True, str(output_path)
output_path.parent.mkdir(parents=True, exist_ok=True)
self.s3_client.download_file(bucket, key, str(output_path))
logger.info(f"Download complete: {output_path}")
@@ -302,8 +306,7 @@ def get_embedding_filenames(
def download_embeddings(
- filenames: List[str], output_path: str, overwrite: bool = False
-) -> List[str]:
+ filenames: List[str], output_path: str, overwrite: bool = False, multit: bool = True) -> List[str]:
"""
Download the embeddings for the intersecting grids from sourcecoop or HTTP and save them to the output directory.
@@ -322,21 +325,39 @@ def download_embeddings(
successful_downloads = 0
if boto3_available:
logger.info("Using S3 for downloading files.")
- client = S3Downloader()
- for fl in filenames:
- success, new_path = client.download_s3_file(fl, output_dir)
- if success:
- new_paths.append(new_path)
- successful_downloads += 1
- logger.info(
- f"Downloaded {len(new_paths)} files to {output_dir} out of {len(filenames)}"
- )
+ if multit:
+ with ProcessPoolExecutor() as executor:
+ future_to_key = {executor.submit(download_wrapper, fl, output_dir): fl for fl in filenames}
+
+ for future in futures.as_completed(future_to_key):
+ key = future_to_key[future]
+ exception = future.exception()
+ if future.result()[0]:
+ new_paths.append(future.result()[1])
+ successful_downloads += 1
+
+ else :
+ client = S3Downloader()
+ for fl in filenames:
+ success, new_path = client.download_s3_file(fl, output_dir)
+ if success:
+ new_paths.append(new_path)
+ successful_downloads += 1
+ logger.info(
+ f"Downloaded {len(new_paths)} files to {output_dir} out of {len(filenames)}"
+ )
else:
logger.info("Using HTTP protocol for downloading files.")
new_paths = http_download(filenames, output_dir, overwrite=overwrite)
logger.info(f"Downloaded {len(new_paths)} files to {output_dir}")
return new_paths
-
+
+def download_wrapper(fl, output_dir):
+ client = S3Downloader()
+
+ success, new_path = client.download_s3_file(fl, output_dir)
+
+ return success, new_path
def patch_vrt_relative_path(path_vrt: Path) -> None:
"""
@@ -391,7 +412,7 @@ def translate_to_cog(path_vrt: Path, version):
return path_out
-def build_files_dataframe(files: list[Path]) -> pd.DataFrame:
+def build_files_dataframe(files: list[Path], version: str) -> pd.DataFrame:
"""
Builds a DataFrame from a list of file paths, extracting relevant metadata from their parent directories.
@@ -411,10 +432,11 @@ def build_files_dataframe(files: list[Path]) -> pd.DataFrame:
# Join extracted parts as new columns
files_df = files_df.join(
parent_parts.rename(
- columns={0: "rest", 1: "version", 2: "_unused", 3: "year", 4: "zone"}
+ columns={0: "rest", 1: "random", 2: "_unused", 3: "year", 4: "zone"}
)
)
- files_df = files_df.drop(columns=["rest", "_unused"])
+ files_df = files_df.drop(columns=["rest", "_unused","random"])
+ files_df["version"] = version
files_df["s3_prefix"] = files_df[["version", "year", "zone"]].apply(
lambda x: "/".join(x), axis=1
)
diff --git a/src/eo_processing/utils/stac_helper.py b/src/eo_processing/utils/stac_helper.py
index e88e2f3..3e43ba0 100644
--- a/src/eo_processing/utils/stac_helper.py
+++ b/src/eo_processing/utils/stac_helper.py
@@ -1,5 +1,12 @@
import pystac_client
-
+import os
+import logging
+from urllib.request import urlopen
+from io import BytesIO
+from typing import Optional, List
+import geopandas as gpd
+import pandas as pd
+import re
def get_stac_collection_url(collection_id: str, catalog_url: str = "https://catalogue.weed.apex.esa.int/") -> str:
"""
@@ -25,7 +32,7 @@ def get_stac_collection_url(collection_id: str, catalog_url: str = "https://cata
def query_modelID_asset_url(model_id: str,
catalog_url: str ="https://catalogue.weed.apex.esa.int/",
- collection_id: str ="model-STAC-v2") -> str:
+ collection_id: str ="model-STAC") -> str:
"""
Queries the URL of a specific asset for a given modelID from a STAC catalog.
@@ -50,4 +57,289 @@ def query_modelID_asset_url(model_id: str,
filter_lang="cql2-json",
)
item_collection = search.item_collection()
- return item_collection.items[0].assets["model_valid_geometry"].href
\ No newline at end of file
+ if not item_collection.items:
+ logging.error(f"No items found for model_id: {model_id}")
+ raise ValueError(f"No items found for model_id: {model_id}")
+
+ return item_collection.items[0].assets["model_valid_geometry"].href
+
+def query_modelID_output_bands(model_id: str,
+ catalog_url: str ="https://catalogue.weed.apex.esa.int/",
+ collection_id: str ="model-STAC") -> List[str]:
+ """
+ Queries the output bands associated with the specified model ID from a STAC catalog.
+
+ This function communicates with a STAC catalog to retrieve the model's output
+ bands by searching for items using their model ID. It uses CQL2 filtering
+ to perform the search and returns the output band names from the first matching
+ item's properties.
+
+ :param model_id: (str) The unique identifier for the model to query.
+ :param catalog_url: (str) The URL of the STAC catalog to connect to. Defaults to
+ "https://catalogue.weed.apex.esa.int/".
+ :param collection_id: (str) The ID of the STAC collection to search within. Defaults
+ to "model-STAC".
+ :return: A list of output band names (List[str]) associated with the given
+ model ID.
+ :raises ValueError: If no items are found in the catalog for the provided model ID.
+ """
+ client = pystac_client.Client.open(catalog_url)
+
+ search = client.search(
+ limit=20,
+ collections=[collection_id],
+ filter={"op": "=", "args": [{"property": "properties.modelID"}, model_id]},
+ filter_lang="cql2-json",
+ )
+ item_collection = search.item_collection()
+ if not item_collection.items:
+ logging.error(f"No items found for model_id: {model_id}")
+ raise ValueError(f"No items found for model_id: {model_id}")
+
+ return item_collection.items[0].properties["output_band_names"]
+
+def query_proba_results(df_AOI: gpd.GeoDataFrame, collection_id:str, processing_year:int,
+ stac_url:str = 'https://catalogue.weed.apex.esa.int',
+ info_debug: bool = True, postprocess: bool = True) -> gpd.GeoDataFrame:
+ """
+ Queries and retrieves PROBA results intersecting with a given Area of Interest (AOI) from a STAC catalog.
+
+ This function searches for PROBA results within the bounding box of the provided AOI and retrieves metadata
+ and asset information from the specified STAC catalog. The retrieved data is reformatted and returned
+ as a GeoDataFrame containing details about the intersecting PROBA tiles.
+
+ Arguments:
+ :param df_AOI: A GeoDataFrame representing the Area of Interest (AOI). The GeoDataFrame should contain geometry
+ information and a coordinate reference system. If the coordinate reference system is not EPSG:4326,
+ the function will reproject it to EPSG:4326.
+ :param collection_id: A string specifying the collection identifier to search within the STAC catalog.
+ :param processing_year: An integer specifying the year for which to retrieve PROBA results.
+ :param stac_url: A string specifying the URL of the STAC catalog to query. Defaults to 'https://catalogue.weed.apex.esa.int'.
+ :param info_debug: A boolean flag to enable or disable debug logging information. Defaults to True.
+ :param postprocess: A boolean flag to determine whether to postprocess the retrieved data. Defaults to True.
+
+ Returns:
+ A GeoDataFrame containing metadata and details about the PROBA results that intersect with the AOI. The GeoDataFrame
+ includes additional columns extracted from the metadata and asset information of the intersecting tiles, such as
+ datetime, bounding box, tile ID, and others.
+
+ Raises:
+ ValueError: If no intersecting PROBA tiles are found in the STAC catalog for the specified collection.
+ """
+ if info_debug: print(f"get_modelID_asset_geometry_from_STAC")
+ # convert AOI into BBOX in 4326
+ if df_AOI.crs != "EPSG:4326":
+ df_AOI_4326 = df_AOI.to_crs("EPSG:4326")
+ else:
+ df_AOI_4326 = df_AOI.copy()
+
+ bbox_4326 = df_AOI_4326.total_bounds
+
+ if info_debug: print(f"- searching for PROBA results in {bbox_4326}")
+ if info_debug: print(f"- using STAC url: {stac_url}")
+ if info_debug: print(f"- using collection id: {collection_id}")
+
+ client = pystac_client.Client.open(stac_url)
+
+ search = client.search(
+ collections=[collection_id],
+ bbox=bbox_4326,
+ fields=["properties", "assets.openEO.href"],
+ )
+
+ results = []
+ for item in search.items_as_dicts():
+ results.append(
+ [item['properties']['datetime'], item['properties']['proj:bbox'], item['properties']['proj:shape'],
+ item['properties']['proj:code'], item['assets']['openEO']['href']])
+
+ # build dataframe
+ df_result = pd.DataFrame(results, columns=['datetime', 'file_bbox', 'file_shape', 'file_epsg', 'file_url'])
+ if info_debug: print(f"- found {len(df_result)} intersecting PROBA tiles")
+ # check if there are any results
+ if df_result.empty:
+ ValueError(f"No intersecting PROBA tiles found in the STAC ({collection_id}).")
+
+ if postprocess:
+ if info_debug: print(f"- postprocessing PROBA tiles")
+ # split out from file_url important parts (file_name, tile_id, etc)
+ df_result['basename'] = df_result['file_url'].apply(lambda x: os.path.basename(x))
+ df_result[['project_typology', 'type', 'processing_year', 'tileID', 'model_short', 'inference_run_version',
+ 'procesisng_start']] = df_result['basename'].str.split('_', expand=True)
+ df_result['processing_year'] = df_result['processing_year'].str[-4:].astype(int)
+
+ # first limit results to processing year
+ df_result = df_result[df_result['processing_year'] == processing_year]
+
+ # check if we have tiles smaller than our standard 20x20km grid - yes then make sure tile name is correct
+ def extract_real_tileid(tileid_variant):
+ """
+ Extract the real tileID by removing trailing letter suffixes.
+ Ensures the tileID ends with a number.
+ Example: '48πXH34a' -> '48πXH34'
+ """
+ # Remove any trailing letters after the last digit
+ return re.sub(r'[a-zA-Z]+$', '', tileid_variant)
+
+ for idx, row in df_result.iterrows():
+ if row.file_shape != [2000, 2000]:
+ df_result.at[idx, 'tileID'] = extract_real_tileid(row.tileID)
+
+ # now we can filter out spatial duplicates for same used modelID_short name
+ # NOTE: that assumes that NEVER different inference runs of smae modelID were saved in same STAC catalog
+ df_result = df_result.drop_duplicates(subset=['tileID', 'model_short'], keep='first')
+
+ # last step. we have to prepare the output file_name.
+ # Step 1: Check if we have duplicate tileIDs with different model_short values
+ duplicate_tiles = df_result.groupby('tileID')['model_short'].apply(lambda x: list(x.unique())).to_dict()
+ tiles_with_multiple_models = {k: v for k, v in duplicate_tiles.items() if len(v) > 1}
+
+ if tiles_with_multiple_models:
+ if info_debug: print(f" -- Found {len(tiles_with_multiple_models)} tiles with multiple model_short values")
+
+ # Step 2: For tiles with multiple models, condense model_short names
+ # Create a condensed model_short by combining unique values
+ for tile, models in tiles_with_multiple_models.items():
+ # Sort models to ensure consistent naming
+ strata = [x.split('-')[0] for x in models]
+ condensed_name = '-'.join(sorted(strata)) + '-' + '-'.join(models[0].split('-')[1:])
+ # Update all rows for this tileID with the condensed name
+ df_result.loc[df_result['tileID'] == tile, 'model_short'] = condensed_name
+
+ # Step 3: Now remove duplicate tileIDs (keeping first occurrence)
+ df_result = df_result.drop_duplicates(subset=['tileID'], keep='first')
+ else:
+ if info_debug: print(" -- No duplicate tileIDs found with different model_short values")
+ # Still remove any exact duplicates
+ df_result = df_result.drop_duplicates(subset=['tileID'], keep='first')
+
+ # Step 4: Create the file_prefix column properly
+ df_result['file_prefix'] = df_result.apply(
+ lambda
+ row: f"{row['project_typology']}_mece-cube_year{row['processing_year']}_{row['tileID']}_{row['model_short']}_{row['inference_run_version']}",
+ axis=1
+ )
+ # filter to final needed
+ df_result = df_result[['tileID', 'file_prefix']]
+
+ return df_result
+
+def get_modelID_asset_geometry_from_STAC(df_AOI: gpd.GeoDataFrame, typology_schema: str = 'IUCNGET',
+ model_version: Optional[str]=None,
+ stac_url:str = 'https://catalogue.weed.apex.esa.int',
+ collection_id:str = 'model-STAC', info_debug: bool = True) -> gpd.GeoDataFrame:
+ """
+ Retrieves geometries associated with model IDs from a STAC (SpatioTemporal Asset Catalog) collection.
+
+ This function queries a STAC catalog to find model records based on an input area of interest (AOI),
+ typology schema, and optionally a specific model version. The geometries associated with the
+ models are retrieved and returned as a GeoDataFrame.
+
+ Parameters
+ ----------
+ df_AOI : gpd.GeoDataFrame
+ The area of interest represented as a GeoDataFrame. Must have a valid coordinate reference
+ system (CRS). If not in EPSG:4326, it will be converted to this CRS.
+ typology_schema : str
+ A typology schema filter (e.g., "IUCNGET") to apply to the STAC search. This helps filter
+ results by their topology type. Default is "IUCNGET".
+ model_version : Optional[str]
+ A specific model version (e.g., "1.0") to filter the search results. If None, all versions
+ are considered. Default is None.
+ stac_url : str
+ The URL of the STAC catalog to query. Default is 'https://catalogue.weed.apex.esa.int'.
+ collection_id : str
+ The ID of the STAC collection to search within. Default is 'model-STAC'.
+ info_debug : bool
+ if additional debug messages should be printed. Default is True.
+
+ Returns
+ -------
+ gpd.GeoDataFrame
+ A GeoDataFrame containing the retrieved model IDs, their properties, and associated geometries.
+ The GeoDataFrame is in EPSG:4326 CRS.
+
+ Raises
+ ------
+ ValueError
+ If no intersecting model IDs are found in the specified STAC collection.
+ """
+ if info_debug: print(f"get_modelID_asset_geometry_from_STAC")
+ # convert AOI into BBOX in 4326
+ if df_AOI.crs != "EPSG:4326":
+ df_AOI_4326 = df_AOI.to_crs("EPSG:4326")
+ else:
+ df_AOI_4326 = df_AOI.copy()
+
+ bbox_4326 = df_AOI_4326.total_bounds
+
+ # init the pySTAC search
+ if info_debug: print(f"- searching for modelIDs in {bbox_4326} with typology schema {typology_schema}")
+ if info_debug: print(f"- using STAC url: {stac_url}")
+ if info_debug: print(f"- using collection id: {collection_id}")
+ client = pystac_client.Client.open(stac_url)
+
+ search = client.search(
+ collections=[collection_id],
+ bbox=bbox_4326,
+ filter={"op": "=", "args": [{"property": "properties.topology"}, typology_schema]},
+ filter_lang="cql2-json",
+ fields=["properties.modelID", "properties.topology", "properties.training_year", "properties.model_version", "properties.name_spatial_region", "properties.name_spatial_zone", "assets.model_valid_geometry"],
+ )
+
+ # run the search
+ results = []
+ for item in search.items_as_dicts():
+ if model_version is None:
+ results.append([item['properties']['modelID'], item['properties']['topology'],
+ item['properties']['training_year'], item['properties']['model_version'],
+ item['properties']['name_spatial_region'], item['properties']['name_spatial_zone'], item['assets']['model_valid_geometry']['href']])
+ else:
+ if float(item['properties']['model_version']) == float(model_version):
+ results.append([item['properties']['modelID'], item['properties']['topology'],
+ item['properties']['training_year'], item['properties']['model_version'],
+ item['properties']['name_spatial_region'], item['properties']['name_spatial_zone'], item['assets']['model_valid_geometry']['href']])
+ # build dataframe
+ df_result = pd.DataFrame(results, columns=['modelID', 'typology', 'training_year', 'model_version', 'spatial_region', 'spatial_zone', 'asset'])
+ if info_debug: print(f"- found {len(df_result)} intersecting modelIDs")
+ # check if there are any results
+ if df_result.empty:
+ ValueError("No intersecting modelIDs found in the modelSTAC.")
+
+ ## load the modelIS asset geometry and add
+ if info_debug: print(f"- loading modelIDS geometries from STAC assets")
+ # Initialize an empty list to store geometries
+ geometries = []
+
+ # Loop over each row in df_result
+ for idx, row in df_result.iterrows():
+ parquet_url = row['asset']
+
+ try:
+ # Download and read parquet file in memory
+ with urlopen(parquet_url) as response:
+ parquet_data = response.read()
+
+ # Read from bytes
+ temp_gdf = gpd.read_parquet(BytesIO(parquet_data))
+
+ # Extract the geometry (assuming there's only one geometry per file)
+ # If there are multiple geometries, you can use .union_all() or take the first one
+ if len(temp_gdf) > 0:
+ geometry = temp_gdf.geometry.iloc[0]
+ else:
+ geometry = None
+
+ geometries.append(geometry)
+
+ if info_debug: print(f"- Successfully loaded geometry for {row['modelID']}")
+
+ except Exception as e:
+ print(f"Error loading {row['modelID']}: {e}")
+ geometries.append(None)
+
+ # Add geometries as a new column to df_result
+ df_result['geometry'] = geometries
+
+ # Convert df_result to a GeoDataFrame
+ return gpd.GeoDataFrame(df_result, geometry='geometry', crs='EPSG:4326')
diff --git a/src/eo_processing/utils/storage.py b/src/eo_processing/utils/storage.py
index c319101..cd8ca98 100644
--- a/src/eo_processing/utils/storage.py
+++ b/src/eo_processing/utils/storage.py
@@ -980,7 +980,7 @@ def QueryItems(self, table: str, lcolumns: List[str]) -> List[tuple]:
return lresults
- def AddColumns(self, table: str, column_names: lst(str)) -> None:
+ def AddColumns(self, table: str, column_names: List[str]) -> None:
"""
Adds a new column to a specified table in the database. The method dynamically constructs
an SQL ALTER TABLE statement to add the column with the specified name, ensuring
@@ -1045,11 +1045,60 @@ def UploadTrainingPointData(self, table: str, df: pd.DataFrame, column_names: Li
data = ReadFaker(bulk_df)
self.BulkInsert(table, data, tuple([column.lower() for column in df.columns]))
- for index, row in line_df.iterrows():
+ if len(line_df) != 0:
#upload line
- self.StatusUpdateTiles(table, (row['MGRSid10'], row['year']),
- tuple([column.lower() for column in column_names]),
- row[column_names].to_list())
+ self.StatusUpdateTilesBulk(table, line_df, column_names)
+
+ def StatusUpdateTilesBulk(self, table:str, df: pd.DataFrame, column_names: List[str]) -> bool:
+ # establish connection to data base
+ # ini connection
+ conn = self.create_connection()
+
+ # set all following in a try loop so if even the pre-processing fails then the connection is closed and rolled back
+ try:
+ # create cursor
+ cur = conn.cursor()
+ # prepare UPDATE statement
+ print('** update the tile status...')
+ for index, row in df.iterrows():
+ PK = (row['MGRSid10'], row['year'])
+ lmsg = row[column_names].to_list()
+
+ # Build the query to update multiple columns in a single SQL statement
+ columns_query = ", ".join([f'"{col.lower()}" = %s' for col in column_names])
+ query_values = tuple(lmsg) + (PK,)
+
+ sql_statement = f"""
+ UPDATE {table}
+ SET {columns_query}
+ WHERE (mgrsid10, year) = %s;
+ """
+ cur.execute(sql_statement, query_values)
+
+ # commit transactions
+ conn.commit()
+ # close cursor
+ cur.close()
+
+ except psycopg.Error as e:
+ print("** Could not update the data in the PostgreSQL database - error...")
+ print(e.pgerror)
+ print(e.pgcode)
+ # excecute a rollback when the error didn't closed the connection
+ try:
+ conn.rollback()
+ except:
+ print('No RollBack possible or not needed!')
+
+ if cur.closed == False: cur.close()
+ return False
+
+ finally:
+ # close connection
+ if conn.closed == 0: conn.close()
+
+ print('** No errors - all data successfully updated.')
+ return True
def StatusUpdateTiles(self, table: str, PK: tuple(str, int), lcolumns: List[str], lmsg: List[str]) -> bool:
"""
@@ -1451,6 +1500,26 @@ def delete_collection(self, collection_name: str) -> None:
else:
print(f"Failed to delete collection: HTTP {resp.status_code}\n{resp.text}")
+ def delete_collection_item(self, collection_name: str, item_id: str) -> None:
+ """
+ Deletes a specified item from a collection in the catalog. This method constructs
+ the URL for the item, uses authentication to validate the request, and sends a
+ delete request to remove the item. If the deletion is successful, a confirmation
+ message is printed. Otherwise, an error message with the appropriate HTTP status
+ code and error details is displayed.
+
+ :param collection_name: The name of the collection containing the item.
+ :param item_id: The ID of the item to delete.
+ """
+ catalog_url = self.get_catalog_url().rstrip('/')
+ auth_token = self.get_bearer_auth()
+ item_url = f"{catalog_url}/collections/{collection_name}/items/{item_id}"
+ resp = delete(item_url, auth=auth_token)
+ if resp.status_code == 204:
+ print(f"Item '{item_id}' deleted successfully from collection '{collection_name}'")
+ else:
+ print(f"Failed to delete item: HTTP {resp.status_code}\n{resp.text}")
+
class ReadFaker:
"""
A class that mimics file reading behavior for data in tabular formats.