add era5 land spain

notebooks/wp5/era5_land_temperature_spain.ipynb

@@ -0,0 +1,292 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0",
+   "metadata": {},
+   "source": [
+    "# Monitor climate change over Europe with land reanalysis data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1",
+   "metadata": {},
+   "source": [
+    "## Import libraries"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import cartopy.crs as ccrs\n",
+    "import fsspec\n",
+    "import geopandas as gpd\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import pymannkendall as mk\n",
+    "import shapely.geometry\n",
+    "import xarray as xr\n",
+    "from c3s_eqc_automatic_quality_control import diagnostics, download, plot\n",
+    "\n",
+    "plt.style.use(\"seaborn-v0_8-notebook\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3",
+   "metadata": {},
+   "source": [
+    "## Set parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Time\n",
+    "year_start = 1997\n",
+    "year_stop = 2022\n",
+    "\n",
+    "# External files\n",
+    "shapefile_url = \"https://www.eea.europa.eu/data-and-maps/data/eea-reference-grids-2/gis-files/spain-shapefile/at_download/file\"\n",
+    "observed_csv = (\n",
+    "    \"observed-annual-average-mean-surface-air-temperature-of-spain-for-1901-2022.csv\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5",
+   "metadata": {},
+   "source": [
+    "## Set request"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "collection_id = \"reanalysis-era5-land-monthly-means\"\n",
+    "request = {\n",
+    "    \"product_type\": \"monthly_averaged_reanalysis\",\n",
+    "    \"variable\": \"2m_temperature\",\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",
+    "    \"time\": \"00:00\",\n",
+    "    \"area\": [44, -10, 36, 0],\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7",
+   "metadata": {},
+   "source": [
+    "## Download data and convert to Celsius"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds = download.download_and_transform(collection_id, request, chunks={\"year\": 1})\n",
+    "da = ds[\"t2m\"]\n",
+    "with xr.set_options(keep_attrs=True):\n",
+    "    da -= 273.15\n",
+    "da.attrs[\"units\"] = \"°C\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9",
+   "metadata": {},
+   "source": [
+    "## Select and cut Spain map"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "10",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def clip_shapefile(data, shapefile_url):\n",
+    "    shapefile_crs = \"EPSG:4326\"\n",
+    "    with fsspec.open(f\"simplecache::{shapefile_url}\") as file:\n",
+    "        gdf = gpd.read_file(file, layer=\"es_100km\").to_crs(shapefile_crs)\n",
+    "\n",
+    "    data = data.rio.set_spatial_dims(x_dim=\"longitude\", y_dim=\"latitude\")\n",
+    "    data = data.rio.write_crs(shapefile_crs)\n",
+    "    data_clip = data.rio.clip(\n",
+    "        gdf.geometry.apply(shapely.geometry.mapping), gdf.crs, drop=False\n",
+    "    )\n",
+    "    return data_clip\n",
+    "\n",
+    "\n",
+    "da = clip_shapefile(da, shapefile_url)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "11",
+   "metadata": {},
+   "source": [
+    "## Plot annual mean"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "12",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "da_time_mean = diagnostics.time_weighted_mean(da)\n",
+    "plot.projected_map(\n",
+    "    da_time_mean.where(da_time_mean), projection=ccrs.PlateCarree(), cmap=\"YlOrRd\"\n",
+    ")\n",
+    "_ = plt.title(\"\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "13",
+   "metadata": {},
+   "source": [
+    "## Plot annual spatial mean"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "14",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "da_spatial_mean = diagnostics.spatial_weighted_mean(da)\n",
+    "da_annual_mean = diagnostics.annual_weighted_mean(da_spatial_mean)\n",
+    "trend, h, p, z, tau, s, var_s, slope, intercept = mk.original_test(da_annual_mean)\n",
+    "\n",
+    "# Plot bars\n",
+    "ax = da_annual_mean.to_pandas().plot.bar()\n",
+    "ax.set_ylabel(f\"{da_annual_mean.attrs['long_name']} [{da_annual_mean.attrs['units']}]\")\n",
+    "ax.bar_label(ax.containers[0], rotation=90, fmt=\"%.2f\", padding=2.5)\n",
+    "plt.show()\n",
+    "\n",
+    "# Plot lines\n",
+    "da_annual_mean.plot(label=\"Data\")\n",
+    "plt.plot(\n",
+    "    da_annual_mean[\"year\"],\n",
+    "    np.arange(da_annual_mean.sizes[\"year\"]) * slope + intercept,\n",
+    "    label=\"Trend Line\",\n",
+    ")\n",
+    "plt.legend()\n",
+    "plt.grid()\n",
+    "plt.title(\"Annual mean\")\n",
+    "plt.show()\n",
+    "\n",
+    "# Print significance\n",
+    "is_significant = p < 0.05\n",
+    "print(f\"The trend is{'' if is_significant else ' NOT'} significant.\")\n",
+    "print(f\"Trend: {slope:f} {da_annual_mean.attrs['units']}/year\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "15",
+   "metadata": {},
+   "source": [
+    "## Comparison with in-situ data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "16",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Read the CSV file into a DataFrame\n",
+    "observed = pd.read_csv(observed_csv)\n",
+    "mask = (observed[\"Category\"] >= year_start) & (observed[\"Category\"] <= year_stop)\n",
+    "observed = observed[mask]\n",
+    "\n",
+    "# Trend and significance\n",
+    "trend, h, p, z, tau, s, var_s, slope, intercept = mk.original_test(\n",
+    "    observed[\"Annual Mean\"]\n",
+    ")\n",
+    "is_significant = p < 0.05\n",
+    "print(f\"The observed trend is{'' if is_significant else ' NOT'} significant.\")\n",
+    "print(f\"Trend: {slope:f} {da_annual_mean.attrs['units']}/year\")\n",
+    "\n",
+    "# bias\n",
+    "bias = np.mean(np.array(da_annual_mean - observed[\"Annual Mean\"]))\n",
+    "print(f\"Bias: {bias} {da_annual_mean.attrs['units']}/year\")\n",
+    "\n",
+    "# Plot the first line\n",
+    "plt.plot(\n",
+    "    observed[\"Category\"],\n",
+    "    observed[\"Annual Mean\"],\n",
+    "    label=\"In-situ temperature\",\n",
+    "    color=\"blue\",\n",
+    "    marker=\"o\",\n",
+    ")\n",
+    "\n",
+    "# Plot the second line\n",
+    "plt.plot(\n",
+    "    da_annual_mean[\"year\"],\n",
+    "    da_annual_mean,\n",
+    "    label=\"ERA5 land temperature \",\n",
+    "    color=\"red\",\n",
+    "    marker=\"x\",\n",
+    ")\n",
+    "\n",
+    "# Add labels and title\n",
+    "plt.xlabel(\"Year\")\n",
+    "plt.ylabel(\"Temperature(°C)\")\n",
+    "plt.title(\" Annual temperature from ERA5-Land and in-situ temperature\")\n",
+    "\n",
+    "# Add legend\n",
+    "plt.legend()\n",
+    "plt.grid()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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
+}