diff --git a/notebooks/wp4/extreme_indices.ipynb b/notebooks/wp4/extreme_indices.ipynb new file mode 100644 index 0000000..6b2102a --- /dev/null +++ b/notebooks/wp4/extreme_indices.ipynb @@ -0,0 +1,944 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "## Extreme indices (temperature/precipitation)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Import packages" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "import math\n", + "import tempfile\n", + "\n", + "import cartopy.crs as ccrs\n", + "import icclim\n", + "import matplotlib.pyplot as plt\n", + "import xarray as xr\n", + "from c3s_eqc_automatic_quality_control import diagnostics, download, plot\n", + "from xarrayMannKendall import Mann_Kendall_test\n", + "\n", + "plt.style.use(\"seaborn-v0_8-notebook\")\n", + "plt.rcParams[\"hatch.linewidth\"] = 0.5" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "## Define Parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "parameters" + ] + }, + "outputs": [], + "source": [ + "# Time period\n", + "year_start = 1971\n", + "year_stop = 2000\n", + "\n", + "# Choose annual or seasonal timeseries\n", + "timeseries = \"JJA\"\n", + "assert timeseries in (\"annual\", \"DJF\", \"MAM\", \"JJA\", \"SON\")\n", + "\n", + "# Variable\n", + "variable = \"precipitation\"\n", + "assert variable in (\"temperature\", \"precipitation\")\n", + "\n", + "# Choose CORDEX or CMIP6\n", + "collection_id = \"CMIP6\"\n", + "assert collection_id in (\"CORDEX\", \"CMIP6\")\n", + "\n", + "# Define region for analysis\n", + "area = [72, -22, 27, 45]\n", + "\n", + "# Define region for request\n", + "cordex_domain = \"europe\"\n", + "\n", + "# Interpolation method\n", + "interpolation_method = \"bilinear\"\n", + "\n", + "# Chunks for download\n", + "chunks = {\"year\": 1}" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Define models and variable-dependent parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "models_cordex = (\n", + " \"clmcom_clm_cclm4_8_17\",\n", + " \"clmcom_eth_cosmo_crclim\",\n", + " \"cnrm_aladin63\",\n", + " \"dmi_hirham5\",\n", + " \"knmi_racmo22e\",\n", + " \"mohc_hadrem3_ga7_05\",\n", + " \"mpi_csc_remo2009\",\n", + " \"smhi_rca4\",\n", + " \"uhoh_wrf361h\",\n", + ")\n", + "\n", + "match variable:\n", + " case \"temperature\":\n", + " resample_reduction = \"max\"\n", + " index_names = (\"SU\", \"TX90p\")\n", + " era5_variable = \"2m_temperature\"\n", + " cordex_variable = \"maximum_2m_temperature_in_the_last_24_hours\"\n", + " cmip6_variable = \"daily_maximum_near_surface_air_temperature\"\n", + " models_cmip6 = (\n", + " \"access_cm2\",\n", + " \"awi_cm_1_1_mr\",\n", + " \"cmcc_esm2\",\n", + " \"cnrm_cm6_1_hr\",\n", + " \"ec_earth3_cc\",\n", + " \"gfdl_esm4\",\n", + " \"inm_cm5_0\",\n", + " \"miroc6\",\n", + " \"mpi_esm1_2_lr\",\n", + " )\n", + " case \"precipitation\":\n", + " resample_reduction = \"sum\"\n", + " index_names = (\"CWD\", \"R20mm\", \"RR1\", \"RX1day\", \"RX5day\")\n", + " era5_variable = \"total_precipitation\"\n", + " cordex_variable = \"mean_precipitation_flux\"\n", + " cmip6_variable = \"precipitation\"\n", + " models_cmip6 = (\n", + " \"access_cm2\",\n", + " \"bcc_csm2_mr\",\n", + " \"cmcc_esm2\",\n", + " \"cnrm_cm6_1_hr\",\n", + " \"ec_earth3_cc\",\n", + " \"gfdl_esm4\",\n", + " \"inm_cm5_0\",\n", + " \"miroc6\",\n", + " \"mpi_esm1_2_lr\",\n", + " )\n", + " case _:\n", + " raise NotImplementedError(f\"{variable=}\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Define ERA5 request" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "request_era = (\n", + " \"reanalysis-era5-single-levels\",\n", + " {\n", + " \"product_type\": \"reanalysis\",\n", + " \"format\": \"netcdf\",\n", + " \"time\": [f\"{hour:02d}:00\" for hour in range(24)],\n", + " \"variable\": era5_variable,\n", + " \"year\": [\n", + " str(year)\n", + " for year in range(year_start - int(timeseries == \"DJF\"), year_stop + 1)\n", + " ], # Include D(year-1)\n", + " \"month\": [f\"{month:02d}\" for month in range(1, 13)],\n", + " \"day\": [f\"{day:02d}\" for day in range(1, 32)],\n", + " \"area\": area,\n", + " },\n", + ")\n", + "\n", + "request_lsm = (\n", + " request_era[0],\n", + " request_era[1]\n", + " | {\n", + " \"year\": \"1940\",\n", + " \"month\": \"01\",\n", + " \"day\": \"01\",\n", + " \"time\": \"00:00\",\n", + " \"variable\": \"land_sea_mask\",\n", + " },\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Define model requests" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "request_cordex = {\n", + " \"format\": \"zip\",\n", + " \"domain\": cordex_domain,\n", + " \"experiment\": \"historical\",\n", + " \"horizontal_resolution\": \"0_11_degree_x_0_11_degree\",\n", + " \"temporal_resolution\": \"daily_mean\",\n", + " \"variable\": cordex_variable,\n", + " \"gcm_model\": \"mpi_m_mpi_esm_lr\",\n", + " \"ensemble_member\": \"r1i1p1\",\n", + " \"area\": area,\n", + "}\n", + "\n", + "request_cmip6 = {\n", + " \"format\": \"zip\",\n", + " \"temporal_resolution\": \"daily\",\n", + " \"experiment\": \"historical\",\n", + " \"variable\": cmip6_variable,\n", + " \"year\": [\n", + " str(year) for year in range(year_start - 1, year_stop + 1)\n", + " ], # Include D(year-1)\n", + " \"month\": [f\"{month:02d}\" for month in range(1, 13)],\n", + " \"day\": [f\"{day:02d}\" for day in range(1, 32)],\n", + " \"area\": area,\n", + "}\n", + "\n", + "\n", + "def get_cordex_years(\n", + " year_start,\n", + " year_stop,\n", + " timeseries,\n", + " start_years=list(range(1951, 2097, 5)),\n", + " end_years=list(range(1955, 2101, 5)),\n", + "):\n", + " start_year = []\n", + " end_year = []\n", + " years = set(range(year_start - int(timeseries == \"DJF\"), year_stop + 1))\n", + " for start, end in zip(start_years, end_years):\n", + " if years & set(range(start, end + 1)):\n", + " start_year.append(start)\n", + " end_year.append(end)\n", + " return start_year, end_year\n", + "\n", + "\n", + "if collection_id == \"CORDEX\":\n", + " models = models_cordex\n", + " model_key = \"rcm_model\"\n", + " request_sim = (\n", + " \"projections-cordex-domains-single-levels\",\n", + " [\n", + " {\n", + " **request_cordex,\n", + " \"start_year\": start_year,\n", + " \"end_year\": end_year,\n", + " }\n", + " for start_year, end_year in zip(\n", + " *get_cordex_years(year_start, year_stop, timeseries)\n", + " )\n", + " ],\n", + " )\n", + "elif collection_id == \"CMIP6\":\n", + " models = models_cmip6\n", + " model_key = \"model\"\n", + " request_sim = (\n", + " \"projections-cmip6\",\n", + " download.split_request(request_cmip6, chunks=chunks),\n", + " )\n", + "else:\n", + " raise ValueError(f\"{collection_id=}\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Functions to cache" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "def select_timeseries(ds, timeseries, year_start, year_stop):\n", + " if timeseries == \"annual\":\n", + " return ds.sel(time=slice(str(year_start), str(year_stop)))\n", + " return ds.sel(time=slice(f\"{year_start-1}-12\", f\"{year_stop}-11\"))\n", + "\n", + "\n", + "def compute_indices(ds, index_names, timeseries, tmpdir):\n", + " labels, datasets = zip(*ds.groupby(\"time.year\"))\n", + " paths = [f\"{tmpdir}/{label}.nc\" for label in labels]\n", + " datasets = [ds.chunk(-1) for ds in datasets]\n", + " xr.save_mfdataset(datasets, paths)\n", + "\n", + " ds = xr.open_mfdataset(paths)\n", + " in_files = f\"{tmpdir}/rechunked.zarr\"\n", + " chunks = {dim: -1 if dim == \"time\" else \"auto\" for dim in ds.dims}\n", + " ds.chunk(chunks).to_zarr(in_files)\n", + "\n", + " datasets = [\n", + " icclim.index(\n", + " index_name=index_name,\n", + " in_files=in_files,\n", + " out_file=f\"{tmpdir}/{index_name}.nc\",\n", + " slice_mode=\"year\" if timeseries == \"annual\" else timeseries,\n", + " )\n", + " for index_name in index_names\n", + " ]\n", + "\n", + " return xr.merge(datasets).drop_dims(\"bounds\")\n", + "\n", + "\n", + "def compute_trends(ds):\n", + " datasets = []\n", + " (lat,) = set(ds.dims) & set(ds.cf.axes[\"Y\"])\n", + " (lon,) = set(ds.dims) & set(ds.cf.axes[\"X\"])\n", + " coords_name = {\n", + " \"time\": \"time\",\n", + " \"y\": lat,\n", + " \"x\": lon,\n", + " }\n", + " for index, da in ds.data_vars.items():\n", + " ds = Mann_Kendall_test(\n", + " da - da.mean(\"time\"),\n", + " alpha=0.05,\n", + " method=\"theilslopes\",\n", + " coords_name=coords_name,\n", + " ).compute()\n", + " ds = ds.rename({k: v for k, v in coords_name.items() if k in ds.dims})\n", + " ds = ds.assign_coords({dim: da[dim] for dim in ds.dims})\n", + " datasets.append(ds.expand_dims(index=[index]))\n", + " ds = xr.concat(datasets, \"index\")\n", + " return ds\n", + "\n", + "\n", + "def add_bounds(ds):\n", + " for coord in {\"latitude\", \"longitude\"} - set(ds.cf.bounds):\n", + " ds = ds.cf.add_bounds(coord)\n", + " return ds\n", + "\n", + "\n", + "def get_grid_out(request_grid_out, method):\n", + " ds_regrid = download.download_and_transform(*request_grid_out)\n", + " coords = [\"latitude\", \"longitude\"]\n", + " if method == \"conservative\":\n", + " ds_regrid = add_bounds(ds_regrid)\n", + " for coord in list(coords):\n", + " coords.extend(ds_regrid.cf.bounds[coord])\n", + " grid_out = ds_regrid[coords]\n", + " coords_to_drop = set(grid_out.coords) - set(coords) - set(grid_out.dims)\n", + " grid_out = ds_regrid[coords].reset_coords(coords_to_drop, drop=True)\n", + " grid_out.attrs = {}\n", + " return grid_out\n", + "\n", + "\n", + "def compute_indices_and_trends(\n", + " ds,\n", + " index_names,\n", + " timeseries,\n", + " year_start,\n", + " year_stop,\n", + " resample_reduction=None,\n", + " request_grid_out=None,\n", + " **regrid_kwargs,\n", + "):\n", + " assert (request_grid_out and regrid_kwargs) or not (\n", + " request_grid_out or regrid_kwargs\n", + " )\n", + " ds = ds.drop_vars([var for var, da in ds.data_vars.items() if len(da.dims) != 3])\n", + " ds = ds[list(ds.data_vars)]\n", + "\n", + " # Original bounds for conservative interpolation\n", + " if regrid_kwargs.get(\"method\") == \"conservative\":\n", + " ds = add_bounds(ds)\n", + " bounds = [\n", + " ds.cf.get_bounds(coord).reset_coords(drop=True)\n", + " for coord in (\"latitude\", \"longitude\")\n", + " ]\n", + " else:\n", + " bounds = []\n", + "\n", + " ds = select_timeseries(ds, timeseries, year_start, year_stop)\n", + " if resample_reduction:\n", + " resampled = ds.resample(time=\"1D\")\n", + " ds = getattr(resampled, resample_reduction)(keep_attrs=True)\n", + " if resample_reduction == \"sum\":\n", + " for da in ds.data_vars.values():\n", + " da.attrs[\"units\"] = f\"{da.attrs['units']} / day\"\n", + " with tempfile.TemporaryDirectory() as tmpdir:\n", + " ds_indices = compute_indices(ds, index_names, timeseries, tmpdir).compute()\n", + " ds_trends = compute_trends(ds_indices)\n", + " ds = ds_indices.mean(\"time\", keep_attrs=True)\n", + " ds = ds.merge(ds_trends)\n", + " if request_grid_out:\n", + " ds = diagnostics.regrid(\n", + " ds.merge({da.name: da for da in bounds}),\n", + " grid_out=get_grid_out(request_grid_out, regrid_kwargs[\"method\"]),\n", + " **regrid_kwargs,\n", + " )\n", + " return ds" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "## Download and transform ERA5" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "transform_func_kwargs = {\n", + " \"index_names\": sorted(index_names),\n", + " \"timeseries\": timeseries,\n", + " \"year_start\": year_start,\n", + " \"year_stop\": year_stop,\n", + "}\n", + "ds_era5 = download.download_and_transform(\n", + " *request_era,\n", + " chunks=chunks,\n", + " transform_chunks=False,\n", + " transform_func=compute_indices_and_trends,\n", + " transform_func_kwargs=transform_func_kwargs\n", + " | {\"resample_reduction\": resample_reduction},\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Download and transform models" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "interpolated_datasets = []\n", + "model_datasets = {}\n", + "for model in models:\n", + " print(f\"{model=}\")\n", + " model_kwargs = {\n", + " \"collection_id\": request_sim[0],\n", + " \"requests\": [request | {model_key: model} for request in request_sim[1]],\n", + " \"chunks\": chunks if collection_id == \"CMIP6\" else None,\n", + " \"transform_chunks\": False,\n", + " \"transform_func\": compute_indices_and_trends,\n", + " }\n", + " # Original model\n", + " model_datasets[model] = download.download_and_transform(\n", + " **model_kwargs,\n", + " transform_func_kwargs=transform_func_kwargs,\n", + " )\n", + "\n", + " # Interpolated model\n", + " ds = download.download_and_transform(\n", + " **model_kwargs,\n", + " transform_func_kwargs=transform_func_kwargs\n", + " | {\n", + " \"request_grid_out\": request_lsm,\n", + " \"method\": interpolation_method,\n", + " \"skipna\": True,\n", + " },\n", + " )\n", + " interpolated_datasets.append(ds.expand_dims(model=[model]))\n", + "\n", + "ds_interpolated = xr.concat(interpolated_datasets, \"model\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Mask land and change attrs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Mask\n", + "lsm = download.download_and_transform(*request_lsm)[\"lsm\"].squeeze(drop=True)\n", + "ds_era5 = ds_era5.where(lsm)\n", + "ds_interpolated = ds_interpolated.where(lsm)\n", + "model_datasets = {\n", + " model: ds.where(diagnostics.regrid(lsm, ds, method=\"bilinear\"))\n", + " for model, ds in model_datasets.items()\n", + "}\n", + "\n", + "# Edit attributes\n", + "for ds in (ds_era5, ds_interpolated, *model_datasets.values()):\n", + " ds[\"trend\"] *= 10\n", + " ds[\"trend\"].attrs = {\"long_name\": \"trend\"}\n", + " for index in index_names:\n", + " ds[index].attrs = {\"long_name\": \"\", \"units\": ds[index].attrs[\"units\"]}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plotting functions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def hatch_p_value(da, ax, **kwargs):\n", + " default_kwargs = {\n", + " \"plot_func\": \"contourf\",\n", + " \"show_stats\": False,\n", + " \"cmap\": \"none\",\n", + " \"add_colorbar\": False,\n", + " \"levels\": [0, 0.05, 1],\n", + " \"hatches\": [\"\", \"/\" * 3],\n", + " }\n", + " kwargs = default_kwargs | kwargs\n", + "\n", + " title = ax.get_title()\n", + " plot_obj = plot.projected_map(da, ax=ax, **kwargs)\n", + " ax.set_title(title)\n", + " return plot_obj\n", + "\n", + "\n", + "def hatch_p_value_ensemble(trend, p_value, ax):\n", + " n_models = trend.sizes[\"model\"]\n", + " robust_ratio = (p_value <= 0.05).sum(\"model\") / n_models\n", + " robust_ratio = robust_ratio.where(p_value.notnull().any(\"model\"))\n", + " signs = xr.concat([(trend > 0).sum(\"model\"), (trend < 0).sum(\"model\")], \"sign\")\n", + " sign_ratio = signs.max(\"sign\") / n_models\n", + " robust_threshold = 0.66\n", + " sign_ratio = sign_ratio.where(robust_ratio > robust_threshold)\n", + " for da, threshold, character in zip(\n", + " [robust_ratio, sign_ratio], [robust_threshold, 0.8], [\"/\", \"\\\\\"]\n", + " ):\n", + " hatch_p_value(da, ax=ax, levels=[0, threshold, 1], hatches=[character * 3, \"\"])\n", + "\n", + "\n", + "def set_extent(da, axs, area):\n", + " extent = [area[i] for i in (1, 3, 2, 0)]\n", + " for i, coord in enumerate(extent):\n", + " extent[i] += -1 if i % 2 else +1\n", + " for ax in axs:\n", + " ax.set_extent(extent)\n", + "\n", + "\n", + "def plot_models(\n", + " data,\n", + " da_for_kwargs=None,\n", + " p_values=None,\n", + " col_wrap=3,\n", + " subplot_kw={\"projection\": ccrs.PlateCarree()},\n", + " figsize=None,\n", + " layout=\"constrained\",\n", + " area=area,\n", + " **kwargs,\n", + "):\n", + " if isinstance(data, dict):\n", + " assert da_for_kwargs is not None\n", + " model_dataarrays = data\n", + " else:\n", + " da_for_kwargs = da_for_kwargs or data\n", + " model_dataarrays = dict(data.groupby(\"model\"))\n", + "\n", + " if p_values is not None:\n", + " model_p_dataarrays = (\n", + " p_values if isinstance(p_values, dict) else dict(p_values.groupby(\"model\"))\n", + " )\n", + " else:\n", + " model_p_dataarrays = None\n", + "\n", + " # Get kwargs\n", + " default_kwargs = {\"robust\": True, \"extend\": \"both\"}\n", + " kwargs = default_kwargs | kwargs\n", + " kwargs = xr.plot.utils._determine_cmap_params(da_for_kwargs.values, **kwargs)\n", + "\n", + " fig, axs = plt.subplots(\n", + " *(col_wrap, math.ceil(len(model_dataarrays) / col_wrap)),\n", + " subplot_kw=subplot_kw,\n", + " figsize=figsize,\n", + " layout=layout,\n", + " )\n", + " axs = axs.flatten()\n", + " for (model, da), ax in zip(model_dataarrays.items(), axs):\n", + " pcm = plot.projected_map(\n", + " da, ax=ax, show_stats=False, add_colorbar=False, **kwargs\n", + " )\n", + " ax.set_title(model)\n", + " if model_p_dataarrays is not None:\n", + " hatch_p_value(model_p_dataarrays[model], ax)\n", + " set_extent(da_for_kwargs, axs, area)\n", + " fig.colorbar(\n", + " pcm,\n", + " ax=axs.flatten(),\n", + " extend=kwargs[\"extend\"],\n", + " location=\"right\",\n", + " label=f\"{da_for_kwargs.attrs.get('long_name', '')} [{da_for_kwargs.attrs.get('units', '')}]\",\n", + " )\n", + " return fig\n", + "\n", + "\n", + "def plot_ensemble(\n", + " da_models,\n", + " da_era5=None,\n", + " p_value_era5=None,\n", + " p_value_models=None,\n", + " subplot_kw={\"projection\": ccrs.PlateCarree()},\n", + " figsize=None,\n", + " layout=\"constrained\",\n", + " cbar_kwargs=None,\n", + " area=area,\n", + " **kwargs,\n", + "):\n", + " # Get kwargs\n", + " default_kwargs = {\"robust\": True, \"extend\": \"both\"}\n", + " kwargs = default_kwargs | kwargs\n", + " kwargs = xr.plot.utils._determine_cmap_params(\n", + " da_models.values if da_era5 is None else da_era5.values, **kwargs\n", + " )\n", + " if da_era5 is None and cbar_kwargs is None:\n", + " cbar_kwargs = {\"orientation\": \"horizontal\"}\n", + "\n", + " # Figure\n", + " fig, axs = plt.subplots(\n", + " *(1 if da_era5 is None else 2, 2),\n", + " subplot_kw=subplot_kw,\n", + " figsize=figsize,\n", + " layout=layout,\n", + " )\n", + " axs = axs.flatten()\n", + " axs_iter = iter(axs)\n", + "\n", + " # ERA5\n", + " if da_era5 is not None:\n", + " ax = next(axs_iter)\n", + " plot.projected_map(\n", + " da_era5, ax=ax, show_stats=False, cbar_kwargs=cbar_kwargs, **kwargs\n", + " )\n", + " if p_value_era5 is not None:\n", + " hatch_p_value(p_value_era5, ax=ax)\n", + " ax.set_title(\"ERA5\")\n", + "\n", + " # Median\n", + " ax = next(axs_iter)\n", + " median = da_models.median(\"model\", keep_attrs=True)\n", + " plot.projected_map(\n", + " median, ax=ax, show_stats=False, cbar_kwargs=cbar_kwargs, **kwargs\n", + " )\n", + " if p_value_models is not None:\n", + " hatch_p_value_ensemble(trend=da_models, p_value=p_value_models, ax=ax)\n", + " ax.set_title(\"Ensemble Median\")\n", + "\n", + " # Bias\n", + " if da_era5 is not None:\n", + " ax = next(axs_iter)\n", + " with xr.set_options(keep_attrs=True):\n", + " bias = median - da_era5\n", + " plot.projected_map(\n", + " bias,\n", + " ax=ax,\n", + " show_stats=False,\n", + " center=0,\n", + " cbar_kwargs=cbar_kwargs,\n", + " **default_kwargs,\n", + " )\n", + " ax.set_title(\"Ensemble Median Bias\")\n", + "\n", + " # Std\n", + " ax = next(axs_iter)\n", + " std = da_models.std(\"model\", keep_attrs=True)\n", + " plot.projected_map(\n", + " std, ax=ax, show_stats=False, cbar_kwargs=cbar_kwargs, **default_kwargs\n", + " )\n", + " ax.set_title(\"Ensemble Standard Deviation\")\n", + "\n", + " set_extent(da_models, axs, area)\n", + " return fig\n", + "\n", + "\n", + "common_title = f\"{year_start=} {year_stop=} {collection_id=} {timeseries=}\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot ensembles" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for index in index_names:\n", + " # Index\n", + " da = ds_interpolated[index]\n", + " fig = plot_ensemble(da_models=da, da_era5=ds_era5[index])\n", + " fig.suptitle(f\"{index}\\n{common_title}\")\n", + " plt.show()\n", + "\n", + " # Trend\n", + " da_trend = ds_interpolated[\"trend\"].sel(index=index)\n", + " da_trend.attrs[\"units\"] = f\"{da.attrs['units']} / decade\"\n", + " fig = plot_ensemble(\n", + " da_models=da_trend,\n", + " da_era5=ds_era5[\"trend\"].sel(index=index),\n", + " p_value_era5=ds_era5[\"p\"].sel(index=index),\n", + " p_value_models=ds_interpolated[\"p\"].sel(index=index),\n", + " center=0,\n", + " )\n", + " fig.suptitle(f\"Trend of {index}\\n{common_title}\")\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot models" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for index in index_names:\n", + " # Index\n", + " da_for_kwargs = ds_era5[index]\n", + " fig = plot_models(\n", + " data={model: ds[index] for model, ds in model_datasets.items()},\n", + " da_for_kwargs=da_for_kwargs,\n", + " )\n", + " fig.suptitle(f\"{index}\\n{common_title}\")\n", + " plt.show()\n", + "\n", + " # Trend\n", + " da_for_kwargs_trends = ds_era5[\"trend\"].sel(index=index)\n", + " da_for_kwargs_trends.attrs[\"units\"] = f\"{da_for_kwargs.attrs['units']} / decade\"\n", + " fig = plot_models(\n", + " data={\n", + " model: ds[\"trend\"].sel(index=index) for model, ds in model_datasets.items()\n", + " },\n", + " da_for_kwargs=da_for_kwargs_trends,\n", + " p_values={\n", + " model: ds[\"p\"].sel(index=index) for model, ds in model_datasets.items()\n", + " },\n", + " center=0,\n", + " )\n", + " fig.suptitle(f\"Trend of {index}\\n{common_title}\")\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot bias" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "with xr.set_options(keep_attrs=True):\n", + " bias = ds_interpolated - ds_era5\n", + "\n", + "for index in index_names:\n", + " # Index bias\n", + " da = bias[index]\n", + " fig = plot_models(data=da, center=0)\n", + " fig.suptitle(f\"Bias of {index}\\n{common_title}\")\n", + " plt.show()\n", + "\n", + " # Trend bias\n", + " da_trend = bias[\"trend\"].sel(index=index)\n", + " da_trend.attrs[\"units\"] = f\"{da.attrs['units']} / decade\"\n", + " fig = plot_models(data=da_trend, center=0)\n", + " fig.suptitle(f\"Trend bias of {index}\\n{common_title}\")\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Boxplot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "weights = collection_id == \"CMIP6\"\n", + "mean_datasets = [\n", + " diagnostics.spatial_weighted_mean(ds.expand_dims(model=[model]), weights=weights)\n", + " for model, ds in model_datasets.items()\n", + "]\n", + "mean_ds = xr.concat(mean_datasets, \"model\")\n", + "mean_bias_ds = diagnostics.spatial_weighted_mean(bias, weights=weights)\n", + "for is_bias, ds in zip((False, True), (mean_ds, mean_bias_ds)):\n", + " for index, da in ds[\"trend\"].groupby(\"index\"):\n", + " df_slope = da.to_dataframe()[[\"trend\"]]\n", + " ax = df_slope.boxplot()\n", + " ax.scatter(\n", + " x=[1] * len(df_slope),\n", + " y=df_slope,\n", + " color=\"grey\",\n", + " marker=\".\",\n", + " label=\"models\",\n", + " )\n", + "\n", + " # Ensemble mean\n", + " ax.scatter(\n", + " x=1,\n", + " y=da.mean(\"model\"),\n", + " marker=\"o\",\n", + " label=f\"{collection_id} Ensemble Mean\",\n", + " )\n", + "\n", + " # ERA5\n", + " labels = [f\"{collection_id} Ensemble\"]\n", + " if not is_bias:\n", + " da = ds_era5[\"trend\"].sel(index=index)\n", + " da = diagnostics.spatial_weighted_mean(da)\n", + " ax.scatter(\n", + " x=2,\n", + " y=da.values,\n", + " marker=\"o\",\n", + " label=\"ERA5\",\n", + " )\n", + " labels.append(\"ERA5\")\n", + "\n", + " ax.set_xticks(range(1, len(labels) + 1), labels)\n", + " ax.set_ylabel(f\"{ds[index].attrs['units']} / decade\")\n", + " plt.suptitle(f\"Trend{' bias ' if is_bias else ' '}of {index}\")\n", + " plt.legend()\n", + " plt.show()" + ] + } + ], + "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.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}