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": {