From dc59cff1bd4b81e87e264154f0045457be9afb35 Mon Sep 17 00:00:00 2001 From: Mattia Almansi Date: Fri, 19 Jul 2024 11:26:26 +0200 Subject: [PATCH] add fapar --- .../wp5/satellite_fapar_completeness.ipynb | 253 ++++++++++++++++++ 1 file changed, 253 insertions(+) create mode 100644 notebooks/wp5/satellite_fapar_completeness.ipynb diff --git a/notebooks/wp5/satellite_fapar_completeness.ipynb b/notebooks/wp5/satellite_fapar_completeness.ipynb new file mode 100644 index 0000000..e586752 --- /dev/null +++ b/notebooks/wp5/satellite_fapar_completeness.ipynb @@ -0,0 +1,253 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# FAPAR satellite data completeness" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## Import packages" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import fsspec\n", + "import geopandas\n", + "import matplotlib.pyplot as plt\n", + "import shapely.geometry\n", + "from c3s_eqc_automatic_quality_control import download, utils\n", + "\n", + "plt.style.use(\"seaborn-v0_8-notebook\")" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "## Define variables" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "# Request parameters\n", + "year_start = 2014\n", + "year_stop = 2014\n", + "nominal_days = [3, 13, 21, 23, 24]\n", + "variables = [\"fapar\"]\n", + "\n", + "# Region\n", + "lon_slice = slice(-13, 35)\n", + "lat_slice = slice(72, 30)\n", + "\n", + "shapefile_url = \"https://figshare.com/ndownloader/files/23392280\"" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "## Set the data request" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "# Define request\n", + "collection_id = \"satellite-lai-fapar\"\n", + "request = {\n", + " \"satellite\": \"proba\",\n", + " \"sensor\": \"vgt\",\n", + " \"horizontal_resolution\": \"1km\",\n", + " \"product_version\": \"V2\",\n", + " \"year\": [str(year) for year in range(year_start, year_stop + 1)],\n", + " \"month\": [f\"{month:02d}\" for month in range(1, 12 + 1)],\n", + " \"nominal_day\": [f\"{day:02d}\" for day in nominal_days],\n", + " \"format\": \"zip\",\n", + " \"variable\": variables,\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## Define function to compute missing values count" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "def compute_missing_values_count(ds):\n", + " ds.rio.set_spatial_dims(x_dim=\"longitude\", y_dim=\"latitude\", inplace=True)\n", + " ds.rio.write_crs(\"epsg:4326\", inplace=True)\n", + "\n", + " da_mvc = ds[\"fAPAR\"].isnull().sum(\"time\") / ds.sizes[\"time\"] * 100\n", + " da_mvc.attrs[\"long_name\"] = \"Missing values\"\n", + " da_mvc.attrs[\"units\"] = \"%\"\n", + "\n", + " return da_mvc.to_dataset(name=\"fAPAR\")" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "## Download and preprocess data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "# Download and cache\n", + "ds = download.download_and_transform(\n", + " collection_id,\n", + " request,\n", + " transform_func=utils.regionalise,\n", + " transform_func_kwargs={\n", + " \"lon_slice\": lon_slice,\n", + " \"lat_slice\": lat_slice,\n", + " },\n", + " chunks={\"year\": 1, \"nominal_day\": 1, \"variable\": 1},\n", + " cached_open_mfdataset_kwargs={\"combine\": \"nested\", \"concat_dim\": \"time\"},\n", + ")\n", + "\n", + "# Shapefile\n", + "with fsspec.open(f\"simplecache::{shapefile_url}\") as file:\n", + " world_shape = geopandas.read_file(file, layer=\"Continents\")" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "## Define plotting function" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "def imshow_and_hist(da, shape):\n", + " \"\"\"Plot map and histogram side-by-side.\n", + "\n", + " Parameters\n", + " ----------\n", + " da: DataArray\n", + " DataArray to plot\n", + " shape: GeoDataFrame\n", + " Geopandas object with polygons\n", + "\n", + " Returns\n", + " -------\n", + " figure, axes\n", + " \"\"\"\n", + " fig, (ax_imshow, ax_hist) = plt.subplots(\n", + " 1, 2, figsize=[10, 5], gridspec_kw={\"width_ratios\": [3, 2]}\n", + " )\n", + "\n", + " da = da.rio.clip(\n", + " shape.geometry.apply(shapely.geometry.mapping),\n", + " shape.crs,\n", + " drop=True,\n", + " )\n", + " da.plot.imshow(ax=ax_imshow)\n", + " ax_imshow.set_title(\"Map\")\n", + "\n", + " da.plot.hist(bins=50, ax=ax_hist)\n", + " ax_hist.set_ylabel(\"Frequency\")\n", + " ax_hist.yaxis.set_label_position(\"right\")\n", + " ax_hist.yaxis.tick_right()\n", + "\n", + " # Compute and show no data percentage\n", + " missing_data_perc = (da == 100).sum() / da.notnull().sum() * 100\n", + " ax_hist.set_title(\n", + " f\"Percentage of pixels with no data respect the total: {float(missing_data_perc):f} %\"\n", + " )\n", + "\n", + " fig.suptitle(\", \".join(list(shape.CONTINENT)))\n", + " return fig, (ax_imshow, ax_hist)" + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": {}, + "source": [ + "## Plot Europe map" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "da_mvc = ds[\"fAPAR\"].isnull().sum(\"time\") / ds.sizes[\"time\"] * 100\n", + "da_mvc.attrs[\"long_name\"] = \"Missing values\"\n", + "da_mvc.attrs[\"units\"] = \"%\"\n", + "\n", + "da_mvc.rio.write_crs(\"epsg:4326\", inplace=True)\n", + "imshow_and_hist(da_mvc, world_shape[world_shape.CONTINENT == \"Europe\"])" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}