diff --git a/notebooks/contributions/Vessel_Density_from_Sentinel-1.ipynb b/notebooks/contributions/Vessel_Density_from_Sentinel-1.ipynb new file mode 100644 index 00000000..da1ff565 --- /dev/null +++ b/notebooks/contributions/Vessel_Density_from_Sentinel-1.ipynb @@ -0,0 +1,873 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "32bdafe2", + "metadata": { + "tags": [ + "papermill-error-cell-tag" + ] + }, + "source": [ + "An Exception was encountered at 'In [5]'." + ] + }, + { + "cell_type": "markdown", + "id": "fd59ebb0", + "metadata": { + "papermill": { + "duration": 0.004526, + "end_time": "2023-01-31T19:15:29.601478", + "exception": false, + "start_time": "2023-01-31T19:15:29.596952", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "# Vessel Density based on Copernicus Sentinel-1\n", + "---\n", + "_Author: Alessandro Cimbelli, Oct 2022_\n", + "\n", + "In this updated proposal, the marine area of interest is divided and analysed into areas of 100 km x 100 km\n", + "#### Methodology \n", + "1. Definition of the boxes ID to be analysed.\n", + "2. Extraction of the footprints of the Sentinel-1 (VV band) images acquired in the time period (one month).\n", + "3. Extraction of vessel centroids in each image\n", + "4. Map of vessel density for each 1-km grid tile of the selected area\n", + "5. Production of geojson and csv files\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "269da3e9", + "metadata": { + "execution": { + "iopub.execute_input": "2023-01-31T19:15:29.609997Z", + "iopub.status.busy": "2023-01-31T19:15:29.609339Z", + "iopub.status.idle": "2023-01-31T19:15:29.618067Z", + "shell.execute_reply": "2023-01-31T19:15:29.617302Z" + }, + "papermill": { + "duration": 0.014794, + "end_time": "2023-01-31T19:15:29.620017", + "exception": false, + "start_time": "2023-01-31T19:15:29.605223", + "status": "completed" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "year = 2020\n", + "months = [1,2,3,4,5,6,7,8,9,10,11,12]" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7b1a4422", + "metadata": { + "execution": { + "iopub.execute_input": "2023-01-31T19:15:29.628053Z", + "iopub.status.busy": "2023-01-31T19:15:29.627232Z", + "iopub.status.idle": "2023-01-31T19:15:29.630878Z", + "shell.execute_reply": "2023-01-31T19:15:29.630241Z" + }, + "papermill": { + "duration": 0.009528, + "end_time": "2023-01-31T19:15:29.632458", + "exception": false, + "start_time": "2023-01-31T19:15:29.622930", + "status": "completed" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "#from edc import check_compatibility\n", + "#check_compatibility(\"user-2022.10-14\", dependencies=[\"SH\", \"XCUBE_GEN\"])" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "4b5a4898", + "metadata": { + "execution": { + "iopub.execute_input": "2023-01-31T19:15:29.640369Z", + "iopub.status.busy": "2023-01-31T19:15:29.639775Z", + "iopub.status.idle": "2023-01-31T19:15:29.643502Z", + "shell.execute_reply": "2023-01-31T19:15:29.642836Z" + }, + "papermill": { + "duration": 0.009586, + "end_time": "2023-01-31T19:15:29.645133", + "exception": false, + "start_time": "2023-01-31T19:15:29.635547", + "status": "completed" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "#from edc import setup_environment_variables\n", + "#setup_environment_variables()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "74013f50", + "metadata": { + "execution": { + "iopub.execute_input": "2023-01-31T19:15:29.653012Z", + "iopub.status.busy": "2023-01-31T19:15:29.652114Z", + "iopub.status.idle": "2023-01-31T19:15:32.248482Z", + "shell.execute_reply": "2023-01-31T19:15:32.247593Z" + }, + "papermill": { + "duration": 2.602875, + "end_time": "2023-01-31T19:15:32.250652", + "exception": false, + "start_time": "2023-01-31T19:15:29.647777", + "status": "completed" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "from xcube_sh.config import CubeConfig\n", + "from xcube_sh.cube import open_cube\n", + "from xcube_sh.sentinelhub import SentinelHub\n", + "from xcube.core.geom import (\n", + " clip_dataset_by_geometry, \n", + " mask_dataset_by_geometry,\n", + " convert_geometry,\n", + " rasterize_features\n", + ")\n", + "import xarray as xr\n", + "import numpy as np\n", + "import shapely.geometry\n", + "import IPython.display\n", + "import calendar\n", + "\n", + "import datetime\n", + "import os, sys\n", + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "import rasterio, rioxarray\n", + "import geopandas as gpd\n", + "import pandas as pd\n", + "import warnings\n", + "import scipy\n", + "import pyproj\n", + "from shapely.ops import transform, cascaded_union, unary_union\n", + "from shapely.geometry import shape, Point, Polygon, LineString, MultiPoint, MultiPolygon, mapping \n", + "from rasterio.features import sieve, shapes\n", + "from rasterio.plot import show, plotting_extent\n", + "from rasterio.enums import Resampling\n", + "import json\n", + "from oauthlib.oauth2 import BackendApplicationClient\n", + "from requests_oauthlib import OAuth2Session\n", + "import requests\n", + "import csv\n", + "import shutil" + ] + }, + { + "cell_type": "markdown", + "id": "411a5f3c", + "metadata": { + "tags": [ + "papermill-error-cell-tag" + ] + }, + "source": [ + "Execution using papermill encountered an exception here and stopped:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "05f95f2c", + "metadata": { + "execution": { + "iopub.execute_input": "2023-01-31T19:15:32.258779Z", + "iopub.status.busy": "2023-01-31T19:15:32.257905Z", + "iopub.status.idle": "2023-01-31T19:15:33.148608Z", + "shell.execute_reply": "2023-01-31T19:15:33.147383Z" + }, + "papermill": { + "duration": 0.896904, + "end_time": "2023-01-31T19:15:33.150613", + "exception": true, + "start_time": "2023-01-31T19:15:32.253709", + "status": "failed" + }, + "tags": [] + }, + "outputs": [ + { + "ename": "DriverError", + "evalue": "grid100k_geo.geojson: No such file or directory", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mCPLE_OpenFailedError\u001b[0m Traceback (most recent call last)", + "File \u001b[0;32mfiona/_shim.pyx:83\u001b[0m, in \u001b[0;36mfiona._shim.gdal_open_vector\u001b[0;34m()\u001b[0m\n", + "File \u001b[0;32mfiona/_err.pyx:291\u001b[0m, in \u001b[0;36mfiona._err.exc_wrap_pointer\u001b[0;34m()\u001b[0m\n", + "\u001b[0;31mCPLE_OpenFailedError\u001b[0m: grid100k_geo.geojson: No such file or directory", + "\nDuring handling of the above exception, another exception occurred:\n", + "\u001b[0;31mDriverError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn [5], line 6\u001b[0m\n\u001b[1;32m 3\u001b[0m client_secret \u001b[38;5;241m=\u001b[39m os\u001b[38;5;241m.\u001b[39menviron[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mSH_CLIENT_SECRET\u001b[39m\u001b[38;5;124m'\u001b[39m]\n\u001b[1;32m 5\u001b[0m grid_name \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mgrid100k_geo.geojson\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m----> 6\u001b[0m grid100k \u001b[38;5;241m=\u001b[39m \u001b[43mgpd\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mread_file\u001b[49m\u001b[43m(\u001b[49m\u001b[43mgrid_name\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 8\u001b[0m res1 \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0.0054\u001b[39m \u001b[38;5;66;03m# 0.00054 # lower resolution for satellite footprints\u001b[39;00m\n\u001b[1;32m 9\u001b[0m res2 \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0.00018\u001b[39m \u001b[38;5;66;03m# Higher resolution for vessel shape and number\u001b[39;00m\n", + "File \u001b[0;32m/opt/conda/envs/edc-default-2022.10-14/lib/python3.9/site-packages/geopandas/io/file.py:253\u001b[0m, in \u001b[0;36m_read_file\u001b[0;34m(filename, bbox, mask, rows, engine, **kwargs)\u001b[0m\n\u001b[1;32m 250\u001b[0m path_or_bytes \u001b[38;5;241m=\u001b[39m filename\n\u001b[1;32m 252\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m engine \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfiona\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[0;32m--> 253\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_read_file_fiona\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 254\u001b[0m \u001b[43m \u001b[49m\u001b[43mpath_or_bytes\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfrom_bytes\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mbbox\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mbbox\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmask\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmask\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrows\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mrows\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\n\u001b[1;32m 255\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 256\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m engine \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpyogrio\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[1;32m 257\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m _read_file_pyogrio(\n\u001b[1;32m 258\u001b[0m path_or_bytes, bbox\u001b[38;5;241m=\u001b[39mbbox, mask\u001b[38;5;241m=\u001b[39mmask, rows\u001b[38;5;241m=\u001b[39mrows, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs\n\u001b[1;32m 259\u001b[0m )\n", + "File \u001b[0;32m/opt/conda/envs/edc-default-2022.10-14/lib/python3.9/site-packages/geopandas/io/file.py:294\u001b[0m, in \u001b[0;36m_read_file_fiona\u001b[0;34m(path_or_bytes, from_bytes, bbox, mask, rows, **kwargs)\u001b[0m\n\u001b[1;32m 291\u001b[0m reader \u001b[38;5;241m=\u001b[39m fiona\u001b[38;5;241m.\u001b[39mopen\n\u001b[1;32m 293\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m fiona_env():\n\u001b[0;32m--> 294\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m \u001b[43mreader\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpath_or_bytes\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mas\u001b[39;00m features:\n\u001b[1;32m 295\u001b[0m \n\u001b[1;32m 296\u001b[0m \u001b[38;5;66;03m# In a future Fiona release the crs attribute of features will\u001b[39;00m\n\u001b[1;32m 297\u001b[0m \u001b[38;5;66;03m# no longer be a dict, but will behave like a dict. So this should\u001b[39;00m\n\u001b[1;32m 298\u001b[0m \u001b[38;5;66;03m# be forwards compatible\u001b[39;00m\n\u001b[1;32m 299\u001b[0m crs \u001b[38;5;241m=\u001b[39m (\n\u001b[1;32m 300\u001b[0m features\u001b[38;5;241m.\u001b[39mcrs[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124minit\u001b[39m\u001b[38;5;124m\"\u001b[39m]\n\u001b[1;32m 301\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m features\u001b[38;5;241m.\u001b[39mcrs \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124minit\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;129;01min\u001b[39;00m features\u001b[38;5;241m.\u001b[39mcrs\n\u001b[1;32m 302\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m features\u001b[38;5;241m.\u001b[39mcrs_wkt\n\u001b[1;32m 303\u001b[0m )\n\u001b[1;32m 305\u001b[0m \u001b[38;5;66;03m# handle loading the bounding box\u001b[39;00m\n", + "File \u001b[0;32m/opt/conda/envs/edc-default-2022.10-14/lib/python3.9/site-packages/fiona/env.py:408\u001b[0m, in \u001b[0;36mensure_env_with_credentials..wrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 405\u001b[0m \u001b[38;5;129m@wraps\u001b[39m(f)\n\u001b[1;32m 406\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mwrapper\u001b[39m(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[1;32m 407\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m local\u001b[38;5;241m.\u001b[39m_env:\n\u001b[0;32m--> 408\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mf\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 409\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 410\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(args[\u001b[38;5;241m0\u001b[39m], \u001b[38;5;28mstr\u001b[39m):\n", + "File \u001b[0;32m/opt/conda/envs/edc-default-2022.10-14/lib/python3.9/site-packages/fiona/__init__.py:264\u001b[0m, in \u001b[0;36mopen\u001b[0;34m(fp, mode, driver, schema, crs, encoding, layer, vfs, enabled_drivers, crs_wkt, **kwargs)\u001b[0m\n\u001b[1;32m 261\u001b[0m path \u001b[38;5;241m=\u001b[39m parse_path(fp)\n\u001b[1;32m 263\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m mode \u001b[38;5;129;01min\u001b[39;00m (\u001b[38;5;124m'\u001b[39m\u001b[38;5;124ma\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mr\u001b[39m\u001b[38;5;124m'\u001b[39m):\n\u001b[0;32m--> 264\u001b[0m c \u001b[38;5;241m=\u001b[39m \u001b[43mCollection\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpath\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdriver\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdriver\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mencoding\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mencoding\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 265\u001b[0m \u001b[43m \u001b[49m\u001b[43mlayer\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mlayer\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43menabled_drivers\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43menabled_drivers\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 266\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m mode \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mw\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[1;32m 267\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m schema:\n\u001b[1;32m 268\u001b[0m \u001b[38;5;66;03m# Make an ordered dict of schema properties.\u001b[39;00m\n", + "File \u001b[0;32m/opt/conda/envs/edc-default-2022.10-14/lib/python3.9/site-packages/fiona/collection.py:162\u001b[0m, in \u001b[0;36mCollection.__init__\u001b[0;34m(self, path, mode, driver, schema, crs, encoding, layer, vsi, archive, enabled_drivers, crs_wkt, ignore_fields, ignore_geometry, **kwargs)\u001b[0m\n\u001b[1;32m 160\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mmode \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mr\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[1;32m 161\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39msession \u001b[38;5;241m=\u001b[39m Session()\n\u001b[0;32m--> 162\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msession\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mstart\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 163\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mmode \u001b[38;5;129;01min\u001b[39;00m (\u001b[38;5;124m'\u001b[39m\u001b[38;5;124ma\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mw\u001b[39m\u001b[38;5;124m'\u001b[39m):\n\u001b[1;32m 164\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39msession \u001b[38;5;241m=\u001b[39m WritingSession()\n", + "File \u001b[0;32mfiona/ogrext.pyx:540\u001b[0m, in \u001b[0;36mfiona.ogrext.Session.start\u001b[0;34m()\u001b[0m\n", + "File \u001b[0;32mfiona/_shim.pyx:90\u001b[0m, in \u001b[0;36mfiona._shim.gdal_open_vector\u001b[0;34m()\u001b[0m\n", + "\u001b[0;31mDriverError\u001b[0m: grid100k_geo.geojson: No such file or directory" + ] + } + ], + "source": [ + "# Your client credentials\n", + "client_id = os.environ['SH_CLIENT_ID']\n", + "client_secret = os.environ['SH_CLIENT_SECRET']\n", + "\n", + "grid_name = \"grid100k_geo.geojson\"\n", + "grid100k = gpd.read_file(grid_name)\n", + "\n", + "res1 = 0.0054 # 0.00054 # lower resolution for satellite footprints\n", + "res2 = 0.00018 # Higher resolution for vessel shape and number\n", + "#dres = (res1-res2)/2\n", + "###############################\n", + "# csv parameters\n", + "eo_sensor = 'Sentinel-1'\n", + "indicator_code = 'E1b'\n", + "y_axis = 'Vessel density'\n", + "indicator_name = 'Vessel density'\n", + "data_provider = \"cimbelli\"\n", + "description = 'Number of vessels in the AOI'\n", + "method = 'Object detection of vessels from Sentinel-1 images'\n", + "update_freq = 'Monthly'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2030021e", + "metadata": { + "papermill": { + "duration": null, + "end_time": null, + "exception": null, + "start_time": null, + "status": "pending" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "# This function extracts the centroids of all the vessels in an image\n", + "def vesselShape(imgtif):\n", + "\n", + " img1tmp = imgtif.astype(rasterio.uint8)\n", + " img2tmp = rasterio.features.sieve(img1tmp.to_numpy(), 15, connectivity=8)\n", + " \n", + " with rasterio.Env():\n", + " image = img2tmp #img_vessels.read(1) # first band\n", + " results = (\n", + " {'properties': {'raster_val': v}, 'geometry': s}\n", + " for i, (s, v) \n", + " in enumerate(\n", + " shapes(image, mask=None, transform=imgtif.rio.transform())))\n", + " geoms = list(results)\n", + " npoly = len(geoms)\n", + " npoly_1 = 0\n", + " \n", + "# for each vessel get shape information\n", + " \n", + " ppc = MultiPoint([])\n", + " if npoly > 0:\n", + " for p in range(npoly):\n", + " val = int(geoms[p]['properties']['raster_val'])\n", + "\n", + " if val==1:\n", + " npoly_1 += 1\n", + " pp = shape(geoms[p]['geometry'])\n", + " pp_centroids = pp.centroid #.coords \n", + " ppc = ppc.union(pp_centroids) \n", + " return ppc, npoly_1\n", + " else:\n", + " return None, 0\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "250f5a60", + "metadata": { + "papermill": { + "duration": null, + "end_time": null, + "exception": null, + "start_time": null, + "status": "pending" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "def sentinelInfo(data, box, n):\n", + "\n", + " # Create a session\n", + " client = BackendApplicationClient(client_id=client_id)\n", + " oauth = OAuth2Session(client=client)\n", + "\n", + " # Get token for the session\n", + " token = oauth.fetch_token(token_url='https://services.sentinel-hub.com/oauth/token', client_secret=client_secret)\n", + "\n", + " # All requests using this session will have an access token automatically added\n", + " resp = oauth.get(\"https://services.sentinel-hub.com/oauth/tokeninfo\") \n", + " token1 = token['access_token']\n", + " \n", + " data = '{\"collections\":[\"sentinel-1-grd\"], \"datetime\":\"' + data + '\",\"bbox\": [' + box + '] ,\"limit\": ' + str(n) + '}' # 2020-01-01T00:00:00Z/2020-01-31T23:59:59Z - 10,54.27,11,54.6\n", + " headers = {\"Content-Type\": \"application/json\" , \"Authorization\": \"Bearer \" + token1 }\n", + " endpoint = \"https://services.sentinel-hub.com/api/v1/catalog/search\"\n", + " \n", + " ggg = requests.post(endpoint, data = data, headers=headers)\n", + " \n", + " res = json.loads(ggg.content.decode('UTF-8'))\n", + " id = res['features'][0]['id']\n", + " sat = id[2:3]\n", + " \n", + " ddatetime = res['features'][0]['properties']['datetime'][:-1]\n", + " \n", + " return sat, ddatetime \n", + " " + ] + }, + { + "cell_type": "markdown", + "id": "10c296d4", + "metadata": { + "papermill": { + "duration": null, + "end_time": null, + "exception": null, + "start_time": null, + "status": "pending" + }, + "tags": [] + }, + "source": [ + "### Definition of the boxes ID to be analysed. Please edit them in the list below the map.\n", + "\n", + "_The mask of all European water areas has been produced from the coastline shapefile available at the page:_\n", + "\n", + "https://www.eea.europa.eu/data-and-maps/data/eea-coastline-for-analysis-1/gis-data/europe-coastline-shapefile\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3df7aecf", + "metadata": { + "papermill": { + "duration": null, + "end_time": null, + "exception": null, + "start_time": null, + "status": "pending" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "%%html\n", + "" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fc1c1582", + "metadata": { + "papermill": { + "duration": null, + "end_time": null, + "exception": null, + "start_time": null, + "status": "pending" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "with open('grid100k_geo.geojson', 'r') as f:\n", + " data = json.load(f)\n", + "\n", + "from ipyleaflet import Map, basemaps, basemap_to_tiles, GeoJSON, Marker, DivIcon, WidgetControl\n", + "m = Map(\n", + " basemap = basemap_to_tiles(basemaps.OpenStreetMap.Mapnik),\n", + " center = (49.0, 9.0),\n", + " zoom = 5\n", + " )\n", + "geo_json = GeoJSON(\n", + " data = data,\n", + " style = {'opacity':1, 'fillOpacity':0.1, 'weight':1, 'fillColor':'green'},\n", + " hover_style = {'color':'blue', 'fillOpacity':0.5, 'weight':1.2}\n", + " )\n", + "for i in range(0,len(grid100k)):\n", + " extx = grid100k.iloc[i][\"geometry\"].bounds[2]-grid100k.iloc[i][\"geometry\"].bounds[0]\n", + " lo = grid100k.iloc[i][\"geometry\"].bounds[0] + extx*0.25 #grid100k.iloc[i][\"geometry\"].centroid.x\n", + " la = grid100k.iloc[i][\"geometry\"].centroid.y\n", + " marker = Marker(\n", + " location = [la, lo], \n", + " icon = DivIcon(html=f\"\"\"
{grid100k.iloc[i]['Aoi_Id']}
\"\"\",),\n", + " draggable=False\n", + " )\n", + " m.add_layer(marker)\n", + "m.add_layer(geo_json)\n", + "m.layout.height = '800px'\n", + "m" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c9b8b293", + "metadata": { + "papermill": { + "duration": null, + "end_time": null, + "exception": null, + "start_time": null, + "status": "pending" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "# TILE LIST\n", + "\n", + "# 'BE1','DE1','DE2','DK2','EE1','GR7',\n", + "# 'ES1','ES18','ES3','ES4','ES5','FR2',\n", + "# 'FR3','FR4','IE2','IT3','IT16','IT17',\n", + "# 'IT18','IT19','LT1','LV1','NL3','NL4',\n", + "# 'NO3','NO4','PL1','PT2','RO1','RU1',\n", + "# 'TR1','TR2'\n", + "\n", + "bbox_list = ('NL3',)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a3a295f5", + "metadata": { + "papermill": { + "duration": null, + "end_time": null, + "exception": null, + "start_time": null, + "status": "pending" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "######################################################\n", + "# Density of vessels for each 1km cell\n", + "######################################################\n", + "\n", + "warnings.filterwarnings('ignore')\n", + "# generate csv\n", + "\n", + "gen_date = datetime.datetime.today().strftime('%Y%m%dT%H%M%S')\n", + "csv_file = indicator_code + '_' + data_provider + '_' + gen_date + '.csv'\n", + "header = ['AOI','Country','Region','City','Site Name','Description','Method','EO Sensor','Input Data', \n", + " 'Indicator code','Time','Measurement Value','Reference Description','Reference time','Reference value','Rule', \n", + " 'Indicator Value','Sub-AOI','Y axis','Indicator Name','Color code','Data Provider','AOI_ID','Update Frequency']\n", + "\n", + "if os.path.exists('out/' + csv_file):\n", + " os.remove('out/' + csv_file) \n", + "with open('out/' + csv_file, 'w', encoding='UTF8') as f:\n", + " writer = csv.writer(f)\n", + " writer.writerow(header)\n", + "\n", + "# start processing\n", + "\n", + " for b in range(len(bbox_list)):\n", + " \n", + " for month in months:\n", + " \n", + " \n", + " d1 = datetime.datetime(year, month, 1,12,0,0)\n", + " d2 = datetime.datetime(year, month, calendar.monthrange(year, month)[1],12,0,0)\n", + " time_interval= d1.strftime(\"%Y-%m-%d\"),d2.strftime(\"%Y-%m-%d\")\n", + " month = str(month) if (month)>9 else \"0\"+ str(month)\n", + " \n", + " arr_s1 = []\n", + " poly = grid100k[grid100k.Aoi_Id == bbox_list[b]]\n", + "\n", + " # --------------------------\n", + " # csv file\n", + " id100 = poly.id_100k.values[0]\n", + " aoi_center = str(round(poly.geometry.centroid.values[0].y,6)) + \", \" + str(round(poly.geometry.centroid.values[0].x,6))\n", + " co = poly.Country.values[0]\n", + " port = poly.Port.values[0]\n", + " city = poly.City.values[0]\n", + " bbox = str(round(poly.geometry.bounds.values[0][0],3)) + ',' + str(round(poly.geometry.bounds.values[0][1],3)) + ',' + str(round(poly.geometry.bounds.values[0][2],3)) + ',' + str(round(poly.geometry.bounds.values[0][3],3)) \n", + "\n", + " # ---------------------------------------------\n", + " tmp_masks = []\n", + " tmp_masks.clear() \n", + " # ---------------------------------------------\n", + " config_s1 = CubeConfig(dataset_name = 'S1GRD',\n", + " band_names = ['VV'],\n", + " tile_size = [512, 512],\n", + " crs = \"http://www.opengis.net/def/crs/EPSG/0/4326\",\n", + " spatial_res = res2, \n", + " bbox = bbox,\n", + " time_range = time_interval, \n", + " ) \n", + " s1 = open_cube(config_s1)\n", + " n_img = len(s1.VV.time)\n", + " print(\"aoi\", bbox_list[b], \" - \",n_img, 'images in the month of', month, '-',year)\n", + " print(\" day\") \n", + " # ---------------------------------------------\n", + " points = MultiPoint([]) \n", + " \n", + " for j in range(n_img): \n", + "\n", + " pp1 = pd.DataFrame()\n", + " img0 = s1.VV[j]\n", + "\n", + " # ---------------------------------------------\n", + " # for footprints\n", + " tmp_mask = img0.notnull() \n", + " tmp_mask = tmp_mask.astype(rasterio.uint8)\n", + " tmp_mask.rio.write_crs(\"EPSG:4326\", inplace=True)\n", + "\n", + " tmp_masks.append(tmp_mask)\n", + "\n", + " # ---------------------------------------------\n", + " # define geojson filename\n", + " daytime = str(s1.VV.time[j].values).replace(\"-\",\"\").replace(\":\",\"\")\n", + " daytime = daytime[:daytime.find(\".\")]\n", + " daytime1 = str(s1.VV.time[j].values).replace(\"-\",\"/\").replace(\"T\",\" \")\n", + " daytime1 = daytime1[:daytime1.find(\".\")]\n", + "\n", + " day = daytime[6:11]\n", + " hour = int(day[-2:])\n", + " hour1 = str(hour) if (hour)>9 else \"0\"+ str(hour)\n", + " print(\" \", day, end=\"..\")\n", + "\n", + "\n", + " # mask image with land and fixed semarks (windfarms, ...)\n", + " img1 = img0 > 0.7\n", + " img2 = img1.astype(rasterio.uint8)\n", + " img2.rio.write_crs(\"EPSG:4326\", inplace=True)\n", + "\n", + " land = gpd.read_file('sea_4326.geojson')\n", + " clipped = img2.rio.clip(land.geometry, land.crs, drop=True, invert=False)\n", + "\n", + " # assess the number of vessels\n", + " pp, nvess = vesselShape(clipped) \n", + " print(nvess, \"vessels\")\n", + "\n", + " # ---------------------------\n", + " # geojson\n", + " # save vessels' centroids from a single image to a geojson file \n", + "\n", + " if not(pp.is_empty):\n", + " if type(pp) == shapely.geometry.point.Point: # only 1 vessel\n", + " points1 = [Point(pp.x, pp.y)]\n", + " pp1 = gpd.GeoDataFrame(geometry=points1, crs=4326)\n", + " pp1.rename(columns ={'geometry':'0'}) \n", + " elif len(pp)>0:\n", + " pp1 = gpd.GeoDataFrame(pp, geometry=0, crs=4326) #'EPSG:4326')\n", + "\n", + "\n", + " # only if there are vessels\n", + " if len(pp1)>0:\n", + "\n", + " # ===================== csv ========================== #\n", + "\n", + " daytime2 = daytime1[:10]\n", + " daytime2 = daytime2.replace('/', '-') + \"T\" + hour1 + \":00:00Z/\" + daytime2.replace('/', '-') + \"T\" + hour1 + \":59:59Z\"\n", + "\n", + " # it could happen 2 images in the same day....\n", + "\n", + " ssat, ddate = sentinelInfo(daytime2, bbox, 1)\n", + " input_data = 'S1' + ssat +' - GRD'\n", + " measurement_value = len(pp1)\n", + " indicator_value = pp.wkt\n", + "\n", + " #['AOI','Country','Region','City','Site Name','Description','Method','EO Sensor','Input Data', \n", + " #'Indicator code','Time','Measurement Value','Reference Description','Reference time','Reference value','Rule', \n", + " #'Indicator Value','Sub-AOI','Y axis','Indicator Name','Color code','Data Provider','AOI_ID','Update Frequency']\n", + "\n", + " info = [aoi_center,co,'',city,port,description,method,eo_sensor,input_data,\n", + " indicator_code,ddate,measurement_value,'',ddate,'','','',\n", + " poly.geometry.values[0],'',indicator_name,'RED',data_provider,bbox_list[b],update_freq]\n", + " writer.writerow(info)\n", + "\n", + " # ---------------------------\n", + "\n", + " # geojson \n", + "\n", + " date_geojson = ddate.replace(\"-\",\"/\").replace(\"T\",\" \")\n", + " date_geojson1 = ddate.replace(\"-\",\"\").replace(\":\",\"\")\n", + " pp1['TIMESTAMP UTC'] = date_geojson\n", + " pp1['CLASSIFICATION'] = \"Vessel\"\n", + "\n", + " geojson_file = indicator_code + \"_\" + str(bbox_list[b]) + \"_\" + date_geojson1 + \".geojson\"\n", + " #print(daytime, geojson_file)\n", + "\n", + " if os.path.exists('out/' + geojson_file):\n", + " os.remove('out/' + geojson_file) \n", + " #pp1.to_crs('EPSG:3035').to_file('out/' + geojson_file, driver='GeoJSON')\n", + " pp1.to_crs('EPSG:4326').to_file('out/' + geojson_file, driver='GeoJSON')\n", + "\n", + " with open('out/' + geojson_file, 'r') as file:\n", + " filedata = file.read()\n", + " filedata = filedata.replace('/', '\\/')\n", + " #remove info on crs\n", + " filedata = filedata.replace('\"crs\": { \"type\": \"name\", \"properties\": { \"name\": \"urn:ogc:def:crs:OGC:1.3:CRS84\" } },\\n', '')\n", + "\n", + " with open('out/' + geojson_file, 'w') as file:\n", + " file.write(filedata)\n", + " # ---------------------------\n", + "\n", + " points = points.union(pp) \n", + "\n", + " else:\n", + " #nv = 0\n", + " #print(\"n\", end=\"..\")\n", + " continue\n", + "\n", + " else:\n", + " #nv = 0\n", + " #print(\"n\", end=\"..\")\n", + " continue\n", + "\n", + " # ---------------------------------------------\n", + " # for footprints\n", + " print(\" footprints\")\n", + " mask = xr.concat(tmp_masks, dim='time')\n", + " mask2 = mask.sum(dim='time')\n", + " mask2 = mask2.astype(rasterio.uint8)\n", + " mask2.rio.write_crs(\"EPSG:4326\", inplace=True)\n", + "\n", + " #vectorize raster of footprints\n", + " with rasterio.Env(): \n", + " results = (\n", + " {'properties': {'raster_val': v}, 'geometry': s}\n", + " for i, (s, v) \n", + " in enumerate(\n", + " shapes(mask2.values, mask=None, connectivity=4, transform=mask2.rio.transform())))\n", + " ftprints = list(results)\n", + "\n", + " #---------------------------------------------------------\n", + " # create geojson from footprints\n", + "\n", + " tmp_list = []\n", + "\n", + " for p2 in range(len(ftprints)):\n", + " val = int(ftprints[p2]['properties']['raster_val'])\n", + " ppp2 = shape(ftprints[p2]['geometry']) \n", + " tmp_list.append({\n", + " 'geometry' : Polygon(ppp2),\n", + " 'val': val\n", + " })\n", + " mpl = gpd.GeoDataFrame(tmp_list, crs='epsg:4326')\n", + " mpl.to_file(\"mpl.geojson\", driver='GeoJSON') \n", + "\n", + " # ---------------------------------------------\n", + "\n", + " ppp = gpd.GeoDataFrame(points,geometry=0, crs='EPSG:4326')\n", + " ppp.insert(1, 't', 0)\n", + " ppp.rename_geometry('geometry').to_file(\"vvv.geojson\", driver='GeoJSON')\n", + " points3 = ppp.to_crs('EPSG:3035')\n", + "\n", + " # count vessels for each vector grid with 1-km pixel\n", + " \n", + " print(\" vessel count\")\n", + " grid1k = gpd.read_file('grid1k/'+ bbox_list[b] + '.geojson')\n", + " if os.path.exists('fin.geojson'):\n", + " os.remove('fin.geojson')\n", + " dfsjoin = gpd.sjoin(grid1k, points3,how=\"inner\") #,predicate=\"within\"\n", + " dfpivot = pd.pivot_table(dfsjoin,index='id',columns='t',aggfunc={'t':len})\n", + " dfpivot.columns = dfpivot.columns.droplevel()\n", + " dfpolynew = grid1k.merge(dfpivot, how='left', on='id')\n", + " df = dfpolynew.rename(columns ={dfpolynew.columns[-1]:'vessels'})\n", + " df.columns = df.columns.astype(str)\n", + " df.to_file(\"fin.geojson\", driver='GeoJSON')\n", + " \n", + " # convert vector grid to raster and calc density\n", + " print(\" density\")\n", + "\n", + " # create rasters of number of vessels and acquisitions\n", + " raster1k = 'grid1k/'+ bbox_list[b] + '.tif'\n", + " !gdal_calc.py --quiet --overwrite -A $raster1k --outfile=n_vessels.tif --calc=\"0\"\n", + " shutil.copy('n_vessels.tif', 'n_acquisitions.tif')\n", + " !gdal_rasterize -q -add -a vessels fin.geojson n_vessels.tif\n", + " !gdal_rasterize -q -add -a val mpl.geojson n_acquisitions.tif\n", + "\n", + "\n", + " # raster of the vessel density\n", + " dens_temp = \"out/density_vessels_temp.tif\"\n", + "\n", + " dens = \"out/density_vessels_\" + bbox_list[b] + \"_\" + str(year) + str(month) + '01_cog.tiff'\n", + " !gdal_calc.py --quiet --overwrite -A n_vessels.tif -B n_acquisitions.tif --outfile=$dens_temp --calc=\"(1000*A)/B\"\n", + " !gdal_translate -q -of COG -co COMPRESS=DEFLATE -co BLOCKSIZE=1024 -co RESAMPLING=AVERAGE -co OVERVIEWS=IGNORE_EXISTING -ot Int16 -a_nodata -9999 $dens_temp $dens\n", + "\n", + "\n", + " if os.path.exists(dens_temp):\n", + " os.remove(dens_temp)\n", + " if os.path.exists('fin.geojson'):\n", + " os.remove('fin.geojson')\n", + " if os.path.exists('vvv.geojson'):\n", + " os.remove('vvv.geojson') \n", + " print(\" ok\")\n", + "\n", + "\n", + "print(\"ok\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "57d04b87", + "metadata": { + "papermill": { + "duration": null, + "end_time": null, + "exception": null, + "start_time": null, + "status": "pending" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "shutil.make_archive(\"output\",'zip',\"out\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c3af2d12", + "metadata": { + "papermill": { + "duration": null, + "end_time": null, + "exception": null, + "start_time": null, + "status": "pending" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "shutil.make_archive(\"grid1k\",'zip',\"grid1k\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "864604ac", + "metadata": { + "papermill": { + "duration": null, + "end_time": null, + "exception": null, + "start_time": null, + "status": "pending" + }, + "tags": [] + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:edc-default-2022.10-14]", + "language": "python", + "name": "conda-env-edc-default-2022.10-14-py" + }, + "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.9.13" + }, + "papermill": { + "default_parameters": {}, + "duration": 5.539753, + "end_time": "2023-01-31T19:15:34.073256", + "environment_variables": {}, + "exception": true, + "input_path": "/tmp/tmpimosmukp", + "output_path": "/tmp/notebook_output.ipynb", + "parameters": {}, + "start_time": "2023-01-31T19:15:28.533503", + "version": "2.3.4" + }, + "properties": { + "authors": [ + { + "id": "d8a2b943-31ef-402a-a1c9-315163cc36f1", + "name": "ryy83337@xcoxc.com" + } + ], + "description": "Vessel Density based on Copernicus Sentinel-1", + "id": "9dc6586e-9aea-485e-acac-755f4eb5ef17", + "license": "Commercial", + "license_notes": null, + "license_price": "5.00", + "name": "Vessel Density based on Copernicus Sentinel-1", + "requirements": [ + "eurodatacube", + "eurodatacube-xcube-gen" + ], + "tags": [ + "EO Data", + "Download Service", + "Jupyter", + "Sentinel Hub", + "Vector Data", + "xcube" + ], + "tosAgree": true, + "type": "Jupyter Notebook", + "version": "0.0.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file