Skip to content

Add radar precipitation datasets overview notebook #29

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
leifdenby wants to merge 4 commits into
base: main
Choose a base branch
from
Draft

Add radar precipitation datasets overview notebook #29

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
c314945
wip radar overview notebook
leifdenby Feb 13, 2026
9d7bd0b
Merge branch 'main' of https://github.com/mlcast-community/mlcast-dat…
leifdenby Feb 13, 2026
5b49344
minor tweak to overview plot
leifdenby Feb 13, 2026
be7a5cc
tweak extent
leifdenby Feb 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions docs/_toc.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,3 +7,4 @@ chapters:
- file: radklim_precipitation
- file: dmi_precipitation
- file: it_dpc_precipitation
- file: radar_precipitation_overview
329 changes: 329 additions & 0 deletions docs/radar_precipitation_overview.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
}
Loading