Skip to content

Commit

Permalink
mask wp4
Browse files Browse the repository at this point in the history
  • Loading branch information
@malmans2
malmans2 committed Jul 18, 2023
1 parent 6d9736b commit e3f51ab
Showing 1 changed file with 42 additions and 20 deletions.
62 changes: 42 additions & 20 deletions notebooks/wp4/clima_and_bias_cordex_cmip6_regionalised.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
"import cartopy.crs as ccrs\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"import scipy\n",
"import skill_metrics\n",
"import xarray as xr\n",
Expand Down Expand Up @@ -60,8 +61,8 @@
"assert variable in (\"precipitation\", \"temperature\")\n",
"\n",
"# Choose CORDEX or CMIP6\n",
"collection_id = \"cordex\"\n",
"assert collection_id in (\"cordex\", \"cmip6\")\n",
"collection_id = \"CORDEX\"\n",
"assert collection_id in (\"CORDEX\", \"CMIP6\")\n",
"\n",
"# Define region for analysis\n",
"lon_slice = slice(-4, 20)\n",
Expand All @@ -70,6 +71,10 @@
"# Define region for request\n",
"cordex_domain = \"europe\"\n",
"\n",
"# Mask out sea grid nodes\n",
"mask_sea = True\n",
"assert isinstance(mask_sea, bool)\n",
"\n",
"# Chunks for download\n",
"chunks = {\"year\": 1}\n",
"assert \"month\" not in chunks, \"Do not use chunks smaller than 1y\""
Expand Down Expand Up @@ -102,6 +107,8 @@
"]\n",
"\n",
"models_cmip6 = [\n",
" \"cmcc_cm2_hr4\",\n",
" \"mpi_esm1_2_lr\",\n",
" \"access_cm2\",\n",
" \"awi_esm_1_1_lr\",\n",
" \"bcc_esm1\",\n",
Expand Down Expand Up @@ -142,6 +149,11 @@
" ], # Include D(year-1)\n",
" \"month\": [f\"{month:02d}\" for month in range(1, 12 + 1)],\n",
" },\n",
")\n",
"\n",
"request_lsm = (\n",
" request_era[0],\n",
" request_era[1] | {\"year\": \"1940\", \"month\": \"01\", \"variable\": \"land_sea_mask\"},\n",
")"
]
},
Expand Down Expand Up @@ -205,7 +217,7 @@
" return start_year, end_year\n",
"\n",
"\n",
"if collection_id == \"cordex\":\n",
"if collection_id.lower() == \"cordex\":\n",
" weights = False # Do not weight spatial statistics/errors\n",
" periodic = False\n",
" models = models_cordex\n",
Expand All @@ -221,7 +233,7 @@
" for start_year, end_year in zip(*get_cordex_years(year_start, year_stop))\n",
" ],\n",
" )\n",
"elif collection_id == \"cmip6\":\n",
"elif collection_id.lower() == \"cmip6\":\n",
" weights = True # Weight spatial statistics/errors\n",
" periodic = True\n",
" models = models_cmip6\n",
Expand Down Expand Up @@ -371,7 +383,7 @@
"ds_timeseries = xr.concat(\n",
" [\n",
" ds_sim.drop_vars(\"height\", errors=\"ignore\"),\n",
" ds_sim.mean(\"model\").expand_dims(model=[\"ensemble\"]),\n",
" ds_sim.mean(\"model\").expand_dims(model=[collection_id]),\n",
" ds_era.drop_vars(\"height\", errors=\"ignore\"),\n",
" ],\n",
" \"model\",\n",
Expand All @@ -382,8 +394,18 @@
" ds_timeseries, lon_slice=lon_slice, lat_slice=lat_slice\n",
")\n",
"\n",
"\n",
"ds = ds_timeseries.mean(\"year\", keep_attrs=True)\n",
"if mask_sea:\n",
" ds_lsm = download.download_and_transform(\n",
" *request_lsm,\n",
" transform_func=utils.regionalise,\n",
" transform_func_kwargs={\"lon_slice\": lon_slice, \"lat_slice\": lat_slice},\n",
" )\n",
" da_lsm = ds_lsm[\"lsm\"].squeeze(\"time\", drop=True)\n",
" ds = ds.where(da_lsm > 0.5)\n",
"ds = ds.compute()\n",
"with xr.set_options(keep_attrs=True):\n",
" ds = ds_timeseries.mean(\"year\").compute()\n",
" ds_bias = ds.drop_sel(model=\"ERA5\") - ds.sel(model=\"ERA5\")\n",
"for da in ds_bias.data_vars.values():\n",
" da.attrs[\"long_name\"] += \" Bias\""
Expand Down Expand Up @@ -412,7 +434,7 @@
")\n",
"projection = Projection(central_longitude=(lon_slice.stop + lon_slice.start) / 2)\n",
"\n",
"da_to_plot = ds[variable].sel(model=[\"ensemble\", \"ERA5\"])\n",
"da_to_plot = ds[variable].sel(model=[collection_id, \"ERA5\"])\n",
"plot_kwargs = xr.plot.utils._determine_cmap_params(\n",
" da_to_plot.values,\n",
" robust=True,\n",
Expand Down Expand Up @@ -451,11 +473,11 @@
")\n",
"plot_kwargs[\"projection\"] = projection\n",
"plot.projected_map(\n",
" ds_bias[variable].sel(model=[\"ensemble\"]), stats_weights=weights, **plot_kwargs\n",
" ds_bias[variable].sel(model=[collection_id]), stats_weights=weights, **plot_kwargs\n",
")\n",
"plt.show()\n",
"facet = plot.projected_map(\n",
" ds_bias[variable].drop_sel(model=\"ensemble\"),\n",
" ds_bias[variable].drop_sel(model=collection_id),\n",
" col=\"model\",\n",
" col_wrap=min(3, ds_bias.sizes[\"model\"] - 1),\n",
" **plot_kwargs,\n",
Expand Down Expand Up @@ -485,8 +507,8 @@
"# Plot\n",
"fig, ax = plt.subplots(1, 1)\n",
"x = np.linspace(\n",
" df_stats[\"ensemble\"][\"mean\"] - 3 * df_stats[\"ensemble\"][\"std\"],\n",
" df_stats[\"ensemble\"][\"mean\"] + 3 * df_stats[\"ensemble\"][\"std\"],\n",
" df_stats[collection_id][\"mean\"] - 3 * df_stats[collection_id][\"std\"],\n",
" df_stats[collection_id][\"mean\"] + 3 * df_stats[collection_id][\"std\"],\n",
" 1_000,\n",
")\n",
"for model, values in da.groupby(\"model\"):\n",
Expand All @@ -495,12 +517,12 @@
" values,\n",
" weights=np.abs(np.cos(np.deg2rad(values[\"latitude\"]))) if weights else None,\n",
" ).evaluate(x)\n",
" plot_kwargs = {\"color\": \"k\", \"ls\": \"--\"} if model == \"ensemble\" else {}\n",
" plot_kwargs = {\"color\": \"k\", \"ls\": \"--\"} if model == collection_id else {}\n",
" ax.plot(x, y, label=model, **plot_kwargs)\n",
"ax.grid()\n",
"ax.set_xlim(x[[0, -1]])\n",
"ax.set_xlabel(f\"{da.attrs['long_name']} [{da.attrs['units']}]\")\n",
"ax.legend()\n",
"ax.legend(bbox_to_anchor=(1.1, 1.05))\n",
"\n",
"# Add stats\n",
"table = plt.table(\n",
Expand All @@ -525,20 +547,19 @@
"metadata": {},
"outputs": [],
"source": [
"ds_stats = xr.concat(\n",
"df_stats = xr.concat(\n",
" [\n",
" diagnostics.spatial_weighted_statistics(\n",
" ds.drop_sel(model=\"ERA5\"), weights=weights\n",
" ),\n",
" diagnostics.spatial_weighted_statistics(ds.sel(model=\"ERA5\"), weights=True),\n",
" ],\n",
" \"model\",\n",
")\n",
"ds_error = diagnostics.spatial_weighted_errors(\n",
")[variable].to_pandas()\n",
"df_error = diagnostics.spatial_weighted_errors(\n",
" ds.drop_sel(model=\"ERA5\"), ds.sel(model=\"ERA5\"), weights=weights\n",
")\n",
"df_stats_and_error = xr.merge([ds_stats, ds_error])[variable].to_pandas()\n",
"df_stats_and_error"
")[variable].to_pandas()\n",
"df_error"
]
},
{
Expand All @@ -555,6 +576,7 @@
"metadata": {},
"outputs": [],
"source": [
"df_stats_and_error = pd.concat([df_stats, df_error])\n",
"tickRMS = np.linspace(0, df_stats_and_error.loc[\"crmse\"].max(), 5, dtype=int)\n",
"tickSTD = np.linspace(0, df_stats_and_error.loc[\"std\"].max(), 5, dtype=int)\n",
"skill_metrics.taylor_diagram(\n",
Expand Down Expand Up @@ -606,7 +628,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.11"
"version": "3.10.12"
},
"vscode": {
"interpreter": {
Expand Down

0 comments on commit e3f51ab

Please sign in to comment.