add maps

bopen · Jan 3, 2024 · 7186896 · 7186896
1 parent 3ec5429
commit 7186896
Showing 1 changed file with 50 additions and 5 deletions.
diff --git a/notebooks/wp4/WIP-maps.ipynb → ...p4/cmip6_sea_ice_concentration_maps.ipynb b/notebooks/wp4/WIP-maps.ipynb → ...p4/cmip6_sea_ice_concentration_maps.ipynb
@@ -9,9 +9,10 @@
     "import math\n",
     "import warnings\n",
     "\n",
+    "import cartopy.crs as ccrs\n",
     "import matplotlib.pyplot as plt\n",
     "import xarray as xr\n",
-    "from c3s_eqc_automatic_quality_control import diagnostics, download\n",
+    "from c3s_eqc_automatic_quality_control import diagnostics, download, plot\n",
     "\n",
     "plt.style.use(\"seaborn-v0_8-notebook\")\n",
     "warnings.filterwarnings(\"ignore\", module=\"cf_xarray\")"
@@ -24,7 +25,16 @@
    "outputs": [],
    "source": [
     "# Time Periods\n",
-    "periods = [slice(\"1985\", \"1994\"), slice(\"1995\", \"2004\"), slice(\"2005\", \"2014\")]\n",
+    "periods = [\n",
+    "    slice(\"1985\", \"1994\"),\n",
+    "    slice(\"1995\", \"2004\"),\n",
+    "    slice(\"2005\", \"2014\"),\n",
+    "]\n",
+    "\n",
+    "\n",
+    "# Regions\n",
+    "regions = (\"northern_hemisphere\", \"southern_hemisphere\")\n",
+    "assert set(regions) <= {\"northern_hemisphere\", \"southern_hemisphere\"}\n",
     "\n",
     "# Choose CMIP6 historical models\n",
     "models = [\n",
@@ -155,7 +165,7 @@
    "outputs": [],
    "source": [
     "cmip6_datasets = {}\n",
-    "for region in (\"northern_hemisphere\", \"southern_hemisphere\"):\n",
+    "for region in regions:\n",
     "    period_datasets = []\n",
     "    for period in periods:\n",
     "        model_datasets = []\n",
@@ -196,7 +206,7 @@
    "outputs": [],
    "source": [
     "satellite_datasets = {}\n",
-    "for region in (\"northern_hemisphere\", \"southern_hemisphere\"):\n",
+    "for region in regions:\n",
     "    period_datasets = []\n",
     "    for period in periods:\n",
     "        print(f\"{region=} {period=}\")\n",
@@ -215,9 +225,44 @@
     "            transform_func_kwargs={\"period\": period},\n",
     "            transform_chunks=False,\n",
     "            chunks={\"year\": 1},\n",
+    "            # Parameters to speed up IO\n",
+    "            concat_dim=\"time\",\n",
+    "            combine=\"nested\",\n",
+    "            data_vars=\"minimal\",\n",
+    "            coords=\"minimal\",\n",
+    "            compat=\"override\",\n",
     "        )\n",
     "        period_datasets.append(ds.expand_dims(period=[f\"{period.start}-{period.stop}\"]))\n",
-    "    cmip6_datasets[region] = xr.concat(period_datasets, \"period\")"
+    "    satellite_datasets[region] = xr.concat(period_datasets, \"period\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for region in regions:\n",
+    "    da_cmip6 = cmip6_datasets[region][\"sic\"].mean(\"model\", keep_attrs=True)\n",
+    "    da_cmip6 = da_cmip6.expand_dims(product=[\"CMIP6\"])\n",
+    "\n",
+    "    da_satellite = satellite_datasets[region][\"sic\"]\n",
+    "    da_satellite = da_satellite.expand_dims(\n",
+    "        product=[request_satellite[1][\"origin\"].replace(\"_\", \"-\").upper()]\n",
+    "    )\n",
+    "    projection = ccrs.Stereographic(\n",
+    "        central_latitude=90 if region.startswith(\"northern\") else -90\n",
+    "    )\n",
+    "    da = xr.concat([da_satellite, da_cmip6], \"product\")\n",
+    "\n",
+    "    plot.projected_map(da, col=\"product\", row=\"period\", projection=projection)\n",
+    "    plt.show()\n",
+    "\n",
+    "    with xr.set_options(keep_attrs=True):\n",
+    "        da_bias = da.diff(\"product\")\n",
+    "    da_bias.attrs[\"long_name\"] = \"Bias of \" + da_bias.attrs[\"long_name\"]\n",
+    "    plot.projected_map(da_bias, row=\"period\", projection=projection, robust=True)\n",
+    "    plt.show()"
    ]
   }
  ],