From 1321424b32d2e3ef81ef1bf515f02d9b21b4d5e4 Mon Sep 17 00:00:00 2001 From: Mattia Almansi Date: Wed, 10 Apr 2024 17:42:29 +0200 Subject: [PATCH] implement maps --- notebooks/wp5/satellite_sea_ice_type.ipynb | 277 ++++++++++++++++----- 1 file changed, 218 insertions(+), 59 deletions(-) diff --git a/notebooks/wp5/satellite_sea_ice_type.ipynb b/notebooks/wp5/satellite_sea_ice_type.ipynb index 03d2280..eaef87f 100644 --- a/notebooks/wp5/satellite_sea_ice_type.ipynb +++ b/notebooks/wp5/satellite_sea_ice_type.ipynb @@ -23,11 +23,13 @@ "metadata": {}, "outputs": [], "source": [ + "import cartopy.crs as ccrs\n", + "import cmocean\n", "import matplotlib.dates as mdates\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import xarray as xr\n", - "from c3s_eqc_automatic_quality_control import download\n", + "from c3s_eqc_automatic_quality_control import download, plot\n", "\n", "plt.style.use(\"seaborn-v0_8-notebook\")" ] @@ -47,9 +49,14 @@ "metadata": {}, "outputs": [], "source": [ - "start_year = 1991\n", - "stop_year = 2020\n", + "# Define time periods\n", + "periods = [\n", + " slice(1991, 2000),\n", + " slice(2001, 2010),\n", + " slice(2011, 2020),\n", + "]\n", "\n", + "# Define path to NERSC data\n", "NERSC_PATH = \"/data/wp5/mangini_fabio/nersc_ice_age_v2p1\"" ] }, @@ -70,26 +77,38 @@ "source": [ "collection_id = \"satellite-sea-ice-edge-type\"\n", "\n", - "common_request = {\n", - " \"cdr_type\": \"cdr\",\n", - " \"variable\": \"sea_ice_type\",\n", - " \"region\": \"northern_hemisphere\",\n", - " \"version\": \"3_0\",\n", - " \"day\": [f\"{day:02d}\" for day in range(1, 32)],\n", - "}\n", "\n", - "requests = [\n", - " common_request\n", - " | {\n", - " \"year\": [str(year) for year in range(start_year, stop_year)],\n", - " \"month\": [f\"{month:02d}\" for month in range(10, 13)],\n", - " },\n", - " common_request\n", - " | {\n", - " \"year\": [str(year + 1) for year in range(start_year, stop_year)],\n", - " \"month\": [f\"{month:02d}\" for month in range(1, 5)],\n", - " },\n", - "]" + "def get_requests(year_start, year_stop):\n", + " common_request = {\n", + " \"cdr_type\": \"cdr\",\n", + " \"variable\": \"sea_ice_type\",\n", + " \"region\": \"northern_hemisphere\",\n", + " \"version\": \"3_0\",\n", + " \"day\": [f\"{day:02d}\" for day in range(1, 32)],\n", + " }\n", + " return [\n", + " common_request\n", + " | {\n", + " \"year\": [str(year) for year in range(year_start, year_stop)],\n", + " \"month\": [f\"{month:02d}\" for month in range(10, 13)],\n", + " },\n", + " common_request\n", + " | {\n", + " \"year\": [str(year + 1) for year in range(year_start, year_stop)],\n", + " \"month\": [f\"{month:02d}\" for month in range(1, 5)],\n", + " },\n", + " ]\n", + "\n", + "\n", + "request_timeseries = get_requests(\n", + " year_start=min([period.start for period in periods]),\n", + " year_stop=max([period.stop for period in periods]),\n", + ")\n", + "\n", + "requests_periods = {}\n", + "for period in periods:\n", + " label = f\"{period.start}-{period.stop}\"\n", + " requests_periods[label] = get_requests(period.start, period.stop)" ] }, { @@ -107,37 +126,54 @@ "metadata": {}, "outputs": [], "source": [ + "def get_nersc_data(ds_cds, nersc_path):\n", + " if nersc_path is None:\n", + " nersc_path = NERSC_PATH\n", + "\n", + " paths = ds_cds[\"time\"].dt.strftime(\n", + " f\"{nersc_path}/%Y/arctic25km_sea_ice_age_v2p1_%Y%m%d.nc\"\n", + " )\n", + " ds_nersc = xr.open_mfdataset(\n", + " set(paths.values),\n", + " concat_dim=\"time\",\n", + " combine=\"nested\",\n", + " data_vars=\"minimal\",\n", + " coords=\"minimal\",\n", + " compat=\"override\",\n", + " )\n", + " ds_nersc = ds_nersc.rename(x=\"xc\", y=\"yc\")\n", + " ds_nersc = ds_nersc.assign_coords(\n", + " {coord: ds_cds[coord] for coord in (\"xc\", \"yc\", \"longitude\", \"latitude\")}\n", + " )\n", + " return ds_nersc\n", + "\n", + "\n", "def get_nersc_multiyear_ice(\n", - " time,\n", + " ds_cds,\n", " use_fyi,\n", " age_threshold,\n", " conc_threshold,\n", " nersc_ice_age_path,\n", "):\n", - " if nersc_ice_age_path is None:\n", - " nersc_ice_age_path = NERSC_PATH\n", - " paths = time.dt.strftime(\n", - " f\"{nersc_ice_age_path}/%Y/arctic25km_sea_ice_age_v2p1_%Y%m%d.nc\"\n", - " )\n", - " ds_age = xr.open_mfdataset(set(paths.values)).rename(x=\"xc\", y=\"yc\")\n", + " ds_nersc = get_nersc_data(ds_cds, nersc_ice_age_path)\n", "\n", " if age_threshold is not None:\n", " assert not use_fyi\n", " assert conc_threshold is None\n", - " return ds_age[\"sia\"] > age_threshold\n", + " return ds_nersc[\"sia\"] > age_threshold\n", "\n", " assert conc_threshold is not None\n", "\n", - " conc_myi = ds_age[\"conc_2yi\"]\n", + " conc_myi = ds_nersc[\"conc_2yi\"]\n", " n = 3\n", - " while (varname := f\"conc_{n}yi\") in ds_age.variables:\n", - " conc_myi += ds_age[varname]\n", + " while (varname := f\"conc_{n}yi\") in ds_nersc.variables:\n", + " conc_myi += ds_nersc[varname]\n", " n += 1\n", "\n", " if not use_fyi:\n", " return conc_myi > conc_threshold\n", "\n", - " conc_fyi = ds_age[\"conc_1yi\"]\n", + " conc_fyi = ds_nersc[\"conc_1yi\"]\n", " return ((conc_myi + conc_fyi) > conc_threshold) & (conc_myi > conc_fyi)\n", "\n", "\n", @@ -145,6 +181,11 @@ " return grid_cell_area * da.sum(dim=dim)\n", "\n", "\n", + "def get_classification_mask(ds, use_ambiguous):\n", + " da = ds[\"ice_type\"]\n", + " return (da >= 3) if use_ambiguous else (da == 3)\n", + "\n", + "\n", "def compute_sea_ice_evaluation_diagnostics(\n", " ds,\n", " use_ambiguous,\n", @@ -157,16 +198,14 @@ " grid_cell_area = (dx**2) * 1.0e-6 # 10^6 km2\n", "\n", " # Masks\n", - " da_cds = ds.cf[\"sea_ice_classification\"]\n", - " da_cds = (da_cds >= 3) if use_ambiguous else (da_cds == 3)\n", + " da_cds = get_classification_mask(ds, use_ambiguous)\n", " da_nersc = get_nersc_multiyear_ice(\n", - " ds[\"time\"],\n", + " ds,\n", " use_fyi=use_fyi,\n", " age_threshold=age_threshold,\n", " conc_threshold=conc_threshold,\n", " nersc_ice_age_path=None,\n", " )\n", - " da_nersc = da_nersc.drop_vars([\"xc\", \"yc\"])\n", "\n", " # Fill variables\n", " units = \"$10^6$km$^2$\"\n", @@ -200,7 +239,34 @@ " \"units\": units,\n", " \"long_name\": \"Integrated ice type error\",\n", " }\n", - " return xr.Dataset(dataarrays)" + " return xr.Dataset(dataarrays)\n", + "\n", + "\n", + "def compute_multiyear_ice_percentage(da):\n", + " da = da.groupby(\"time.month\").map(\n", + " lambda da: 100 * da.sum(\"time\") / da.sizes[\"time\"]\n", + " )\n", + " da.attrs = {\n", + " \"units\": \"%\",\n", + " \"long_name\": \"Multi-year sea ice percentage\",\n", + " }\n", + " return da.to_dataset(name=\"percentage\")\n", + "\n", + "\n", + "def compute_cds_multiyear_ice_percentage(ds, use_ambiguous):\n", + " da = get_classification_mask(ds, use_ambiguous)\n", + " return compute_multiyear_ice_percentage(da)\n", + "\n", + "\n", + "def compute_nersc_multiyear_ice_percentage(ds, use_fyi, age_threshold, conc_threshold):\n", + " da = get_nersc_multiyear_ice(\n", + " ds,\n", + " use_fyi=use_fyi,\n", + " age_threshold=age_threshold,\n", + " conc_threshold=conc_threshold,\n", + " nersc_ice_age_path=None,\n", + " )\n", + " return compute_multiyear_ice_percentage(da)" ] }, { @@ -218,20 +284,61 @@ "metadata": {}, "outputs": [], "source": [ + "nersc_kwargs = {\n", + " \"use_fyi\": True,\n", + " \"age_threshold\": None,\n", + " \"conc_threshold\": 0.15,\n", + "}\n", + "kwargs = {\n", + " \"chunks\": {\"year\": 1},\n", + " \"concat_dim\": \"time\",\n", + " \"combine\": \"nested\",\n", + " \"data_vars\": \"minimal\",\n", + " \"coords\": \"minimal\",\n", + " \"compat\": \"override\",\n", + "}\n", + "\n", + "# NERSC Maps\n", "datasets = []\n", - "for use_ambiguous in (True, False):\n", - " print(f\"{use_ambiguous=}\")\n", + "for period, requests in requests_periods.items():\n", + " print(f\"NERSC Maps: {period=}\")\n", " ds = download.download_and_transform(\n", " collection_id,\n", " requests,\n", + " transform_func=compute_nersc_multiyear_ice_percentage,\n", + " transform_func_kwargs=nersc_kwargs,\n", + " transform_chunks=False,\n", + " **kwargs,\n", + " )\n", + " datasets.append(ds.expand_dims(period=[period]))\n", + "ds_nersc = xr.concat(datasets, \"period\")\n", + "\n", + "# CDS Maps\n", + "datasets = []\n", + "for period, requests in requests_periods.items():\n", + " for use_ambiguous in (True, False):\n", + " print(f\"CDS Maps: {period=} {use_ambiguous=}\")\n", + " ds = download.download_and_transform(\n", + " collection_id,\n", + " requests,\n", + " transform_func=compute_cds_multiyear_ice_percentage,\n", + " transform_func_kwargs={\"use_ambiguous\": use_ambiguous},\n", + " transform_chunks=False,\n", + " **kwargs,\n", + " )\n", + " datasets.append(ds.expand_dims(use_ambiguous=[use_ambiguous], period=[period]))\n", + "ds_cds = xr.merge(datasets)\n", + "\n", + "# Timeseries\n", + "datasets = []\n", + "for use_ambiguous in (True, False):\n", + " print(f\"Timeseries: {use_ambiguous=}\")\n", + " ds = download.download_and_transform(\n", + " collection_id,\n", + " request_timeseries,\n", " transform_func=compute_sea_ice_evaluation_diagnostics,\n", - " transform_func_kwargs={\n", - " \"use_ambiguous\": use_ambiguous,\n", - " \"use_fyi\": True,\n", - " \"age_threshold\": None,\n", - " \"conc_threshold\": 0.15,\n", - " },\n", - " chunks={\"year\": 1},\n", + " transform_func_kwargs=nersc_kwargs | {\"use_ambiguous\": use_ambiguous},\n", + " **kwargs,\n", " )\n", " datasets.append(ds.expand_dims(use_ambiguous=[use_ambiguous]).compute())\n", "ds_timeseries = xr.concat(datasets, \"use_ambiguous\")" @@ -264,19 +371,15 @@ " return ds.set_index(time=tuple(coords)).unstack(\"time\")\n", "\n", "\n", - "def plot_against_monthday(\n", - " ds,\n", - " cmap=\"viridis\",\n", - " **kwargs,\n", - "):\n", + "def plot_against_monthday(ds, cmap=\"viridis\", **kwargs):\n", " defaults = {\n", " \"row\": \"variable\",\n", " \"x\": \"time\",\n", " \"hue\": \"year\",\n", " \"add_legend\": False,\n", + " \"figsize\": (10, 10),\n", + " \"sharey\": False,\n", " }\n", - " kwargs = defaults | kwargs\n", - "\n", " ds = rearrange_year_vs_monthday(ds)\n", "\n", " da = ds.to_array()\n", @@ -290,7 +393,7 @@ "\n", " colors = plt.get_cmap(cmap, da.sizes[\"year\"]).colors\n", " with plt.rc_context({\"axes.prop_cycle\": plt.cycler(color=colors)}):\n", - " facet = da.plot(**kwargs)\n", + " facet = da.plot(**(defaults | kwargs))\n", "\n", " for ax, sel_dict in zip(facet.axs.flatten(), facet.name_dicts.flatten()):\n", " ax.grid()\n", @@ -309,6 +412,24 @@ " *_, variable = label.get_text().split()\n", " long_name = ds[variable].attrs[\"long_name\"].replace(\"from\", \"from\\n\")\n", " label.set_text(long_name)\n", + " return facet\n", + "\n", + "\n", + "def plot_maps(da, coastline_color=\"limegreen\", **kwargs):\n", + " defaults = {\n", + " \"row\": \"period\",\n", + " \"col\": \"month\",\n", + " \"projection\": ccrs.Stereographic(central_latitude=90.0),\n", + " \"cmap\": cmocean.cm.ice,\n", + " }\n", + " da = da.sel(\n", + " xc=slice(-2.5e3, 2.5e3),\n", + " yc=slice(2.5e3, -2.5e3),\n", + " month=[10, 11, 12, 1, 2, 3, 4],\n", + " )\n", + " facet = plot.projected_map(da, **(defaults | kwargs))\n", + " for ax in facet.axs.flatten():\n", + " ax.coastlines(color=coastline_color, lw=1)\n", " return facet" ] }, @@ -327,11 +448,49 @@ "metadata": {}, "outputs": [], "source": [ - "facet = plot_against_monthday(\n", - " ds_timeseries, col=\"use_ambiguous\", sharey=False, figsize=(10, 10)\n", - ")\n", + "facet = plot_against_monthday(ds_timeseries, col=\"use_ambiguous\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "15", + "metadata": {}, + "source": [ + "## Plot NERSC" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16", + "metadata": {}, + "outputs": [], + "source": [ + "facet = plot_maps(ds_nersc[\"percentage\"])\n", "plt.show()" ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "## Plot CDS" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "for use_ambiguous, ds in ds_cds.groupby(\"use_ambiguous\", squeeze=False):\n", + " facet = plot_maps(ds[\"percentage\"])\n", + " facet.fig.suptitle(f\"{use_ambiguous=}\")\n", + " plt.show()" + ] } ], "metadata": {