diff --git a/docs/_toc.yml b/docs/_toc.yml index f93a10a..a517deb 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -7,3 +7,4 @@ chapters: - file: radklim_precipitation - file: dmi_precipitation - file: it_dpc_precipitation +- file: radar_precipitation_overview diff --git a/docs/radar_precipitation_overview.ipynb b/docs/radar_precipitation_overview.ipynb new file mode 100644 index 0000000..0b7aa41 --- /dev/null +++ b/docs/radar_precipitation_overview.ipynb @@ -0,0 +1,329 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Overview of Sub-Hourly Radar-Precipitation Datasets\n", + "\n", + "This notebook gives a quick overview of all sub-hourly datasets in the `precipitation` catalog.\n", + "\n", + "It creates two figures:\n", + "1. A spatial overview at one timestep that exists in all selected datasets, with a red bounding box for each dataset domain.\n", + "2. A timeline plot of each dataset's temporal coverage." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "from matplotlib.colors import to_rgba\n", + "from matplotlib.patches import Patch\n", + "import numpy as np\n", + "import pandas as pd\n", + "import cartopy.crs as ccrs\n", + "\n", + "import mlcast_datasets" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "cat = mlcast_datasets.open_catalog()\n", + "source_names = list(cat.precipitation)\n", + "source_names" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "def infer_time_step_minutes(ds, sample_size=128):\n", + " sample = pd.DatetimeIndex(ds[\"time\"].isel(time=slice(0, sample_size)).values)\n", + " if len(sample) < 2:\n", + " return np.nan\n", + "\n", + " deltas = sample.to_series().diff().dropna()\n", + " if deltas.empty:\n", + " return np.nan\n", + "\n", + " return deltas.median().total_seconds() / 60.0\n", + "\n", + "\n", + "def select_plot_variable(ds):\n", + " for var_name, da in ds.data_vars.items():\n", + " if \"time\" in da.dims and da.ndim >= 3:\n", + " return var_name\n", + " raise ValueError(\"Could not find a 2D spatial variable with a time dimension.\")\n", + "\n", + "\n", + "datasets = {}\n", + "summary = []\n", + "\n", + "for name in source_names:\n", + " ds = cat.precipitation[name].to_dask()\n", + " step_min = infer_time_step_minutes(ds)\n", + " datasets[name] = ds\n", + "\n", + " t = pd.DatetimeIndex(ds[\"time\"].values)\n", + " summary.append(\n", + " {\n", + " \"dataset\": name,\n", + " \"time_step_minutes\": step_min,\n", + " \"start\": t.min(),\n", + " \"end\": t.max(),\n", + " \"n_times\": len(t),\n", + " }\n", + " )\n", + "\n", + "summary_df = pd.DataFrame(summary).sort_values(\"dataset\").reset_index(drop=True)\n", + "summary_df" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "subhourly = summary_df[summary_df[\"time_step_minutes\"] < 60].copy()\n", + "subhourly_names = subhourly[\"dataset\"].tolist()\n", + "\n", + "if not subhourly_names:\n", + " raise ValueError(\"No sub-hourly datasets found in cat.precipitation.\")\n", + "\n", + "subhourly" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "time_indexes = {\n", + " name: pd.DatetimeIndex(datasets[name][\"time\"].values) for name in subhourly_names\n", + "}\n", + "\n", + "common_times = None\n", + "for name in subhourly_names:\n", + " if common_times is None:\n", + " common_times = time_indexes[name]\n", + " else:\n", + " common_times = common_times.intersection(time_indexes[name])\n", + "\n", + "if common_times is None or len(common_times) == 0:\n", + " raise ValueError(\"No common timestamp found across all sub-hourly datasets.\")\n", + "\n", + "common_time = common_times[-1]\n", + "common_time" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "def get_data_crs(ds, var_name):\n", + " grid_mapping_name = ds[var_name].attrs.get(\"grid_mapping\")\n", + " if grid_mapping_name and grid_mapping_name in ds:\n", + " crs_wkt = ds[grid_mapping_name].attrs.get(\"crs_wkt\")\n", + " if crs_wkt:\n", + " return ccrs.Projection(crs_wkt)\n", + " return ccrs.PlateCarree()\n", + "\n", + "\n", + "def get_domain_bounds(da):\n", + " spatial_dims = [d for d in da.dims if d != \"time\"]\n", + " if len(spatial_dims) != 2:\n", + " raise ValueError(f\"Expected two spatial dims, got {spatial_dims}\")\n", + "\n", + " y_dim, x_dim = spatial_dims\n", + " x = da.coords.get(x_dim, da[x_dim])\n", + " y = da.coords.get(y_dim, da[y_dim])\n", + "\n", + " if x.ndim != 1 or y.ndim != 1:\n", + " raise ValueError(\n", + " \"Expected 1D coordinates for spatial dims; this helper currently supports rectilinear grids.\"\n", + " )\n", + "\n", + " xmin = float(x.min().values)\n", + " xmax = float(x.max().values)\n", + " ymin = float(y.min().values)\n", + " ymax = float(y.max().values)\n", + " return xmin, xmax, ymin, ymax" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## 1) Spatial overview at one common timestep" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "# create equal area projection over Europe for the plot\n", + "crs = ccrs.LambertAzimuthalEqualArea(central_longitude=10, central_latitude=50)\n", + "fig, ax = plt.subplots(\n", + " 1,\n", + " 1,\n", + " figsize=(12, 8),\n", + " subplot_kw={\"projection\": crs},\n", + ")\n", + "\n", + "palette = plt.cm.tab10(np.linspace(0, 1, len(subhourly_names)))\n", + "legend_handles = []\n", + "\n", + "for i, name in enumerate(subhourly_names):\n", + " ds = datasets[name]\n", + " var_name = select_plot_variable(ds)\n", + " data_crs = get_data_crs(ds, var_name)\n", + "\n", + " # Use first timestep to define the dataset spatial domain via non-NaN points.\n", + " da0 = ds[var_name].isel(time=0)\n", + " domain_mask = da0.notnull().astype(float).where(da0.notnull())\n", + " # print number of valid points in domain_mask\n", + " print(f\"{name}: {domain_mask.sum().values} valid points in domain mask\")\n", + "\n", + " color = palette[i]\n", + " spatial_dims = [d for d in domain_mask.dims if d != \"time\"]\n", + " y_dim, x_dim = spatial_dims\n", + " x = domain_mask[x_dim].values\n", + " y = domain_mask[y_dim].values\n", + " z = domain_mask.values\n", + "\n", + " ax.contourf(\n", + " x,\n", + " y,\n", + " z,\n", + " levels=[0.5, 1.5],\n", + " colors=[to_rgba(color, alpha=0.35)],\n", + " transform=data_crs,\n", + " antialiased=True,\n", + " )\n", + "\n", + " legend_handles.append(\n", + " Patch(facecolor=to_rgba(color, alpha=0.45), edgecolor=\"none\", label=name)\n", + " )\n", + "\n", + "ax.coastlines(resolution=\"50m\", color=\"black\", linewidth=0.8)\n", + "ax.gridlines(draw_labels=[\"left\", \"bottom\"], linestyle=\":\", alpha=0.5)\n", + "ax.set_title(\n", + " f\"Sub-hourly radar-precipitation dataset domains (mask from first timestep)\"\n", + ")\n", + "ax.legend(handles=legend_handles, title=\"Datasets\", loc=\"upper right\")\n", + "# set extent to cover Europe, with adjustable padding\n", + "lon_width = 45\n", + "lat_height = 35\n", + "center = [10, 53]\n", + "ax.set_extent(\n", + " [\n", + " center[0] - lon_width / 2,\n", + " center[0] + lon_width / 2,\n", + " center[1] - lat_height / 2,\n", + " center[1] + lat_height / 2,\n", + " ],\n", + " crs=ccrs.PlateCarree(),\n", + ")\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "## 2) Temporal extent overview" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(12, 1.4 * len(subhourly_names) + 2))\n", + "\n", + "starts = []\n", + "ends = []\n", + "\n", + "for i, name in enumerate(subhourly_names):\n", + " t = time_indexes[name]\n", + " start = t.min()\n", + " end = t.max()\n", + " starts.append(start)\n", + " ends.append(end)\n", + "\n", + " ax.hlines(i, start, end, linewidth=8, alpha=0.9)\n", + " ax.plot([start, end], [i, i], \"|\", color=\"black\", markersize=12)\n", + "\n", + "overlap_start = max(starts)\n", + "overlap_end = min(ends)\n", + "if overlap_start <= overlap_end:\n", + " ax.axvspan(\n", + " overlap_start,\n", + " overlap_end,\n", + " color=\"0.85\",\n", + " alpha=0.5,\n", + " label=\"Common coverage window\",\n", + " )\n", + "\n", + "ax.set_yticks(range(len(subhourly_names)))\n", + "ax.set_yticklabels(subhourly_names)\n", + "ax.set_xlabel(\"Time\")\n", + "ax.set_ylabel(\"Dataset\")\n", + "ax.set_title(\"Temporal extent of sub-hourly radar-precipitation datasets\")\n", + "ax.grid(axis=\"x\", linestyle=\":\", alpha=0.5)\n", + "if overlap_start <= overlap_end:\n", + " ax.legend(loc=\"upper left\")\n", + "\n", + "plt.tight_layout()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "mlcast-datasets", + "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.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}