new steric

Author: malmans2

notebooks/wp5/steric_sea_level_ora5.ipynb

@@ -24,8 +24,9 @@
    "outputs": [],
    "source": [
     "import gsw_xarray as gsw\n",
+    "import matplotlib.pyplot as plt\n",
     "import xarray as xr\n",
-    "from c3s_eqc_automatic_quality_control import download"
+    "from c3s_eqc_automatic_quality_control import diagnostics, download"
    ]
   },
   {
@@ -44,7 +45,7 @@
    "outputs": [],
    "source": [
     "year_start = 2004\n",
-    "year_stop = 2004  # 2023"
+    "year_stop = 2023"
    ]
   },
   {
@@ -102,12 +103,31 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def compute_insitu_density(ds):\n",
+    "def compute_densities(ds):\n",
     "    p = gsw.p_from_z(-ds[\"deptht\"], ds[\"latitude\"])\n",
     "    SA = gsw.SA_from_SP(ds[\"vosaline\"], p, ds[\"longitude\"], ds[\"latitude\"])\n",
     "    CT = gsw.CT_from_pt(SA, ds[\"votemper\"])\n",
-    "    rho = gsw.rho(SA, CT, p)\n",
-    "    return xr.merge([p, SA, CT, rho])"
+    "    rho = gsw.rho(SA, CT, p).rename(\"rho\")\n",
+    "    trho = gsw.rho(SA.mean(\"time\"), CT, p).rename(\"trho\")\n",
+    "    srho = gsw.rho(SA, CT.mean(\"time\"), p).rename(\"hrho\")\n",
+    "    return diagnostics.spatial_weighted_mean(xr.merge([rho, trho, srho]))\n",
+    "\n",
+    "\n",
+    "def compute_ssl(rho, rho0=1025):\n",
+    "    grouped_rho = rho.groupby(\"time.month\")\n",
+    "    delta_rho = grouped_rho - grouped_rho.mean()\n",
+    "    ssl = -(delta_rho.fillna(0) / rho0).integrate(\"deptht\")\n",
+    "    match rho.name:\n",
+    "        case \"rho\":\n",
+    "            prefix = \"\"\n",
+    "        case \"trho\":\n",
+    "            prefix = \"thermo\"\n",
+    "        case \"hrho\":\n",
+    "            prefix = \"halo\"\n",
+    "        case _:\n",
+    "            raise NotImplementedError(f\"{rho.name=}\")\n",
+    "    ssl.attrs = {\"long_name\": f\"{prefix}steric sea level\".title(), \"units\": \"m\"}\n",
+    "    return ssl.rename(rho.name.replace(\"rho\", \"ssl\"))"
    ]
   },
   {
@@ -127,7 +147,7 @@
    "source": [
     "ds_reanalysis = download.download_and_transform(\n",
     "    *request_reanalysis,\n",
-    "    transform_func=compute_insitu_density,\n",
+    "    transform_func=compute_densities,\n",
     "    **download_kwargs,\n",
     ")"
    ]
@@ -137,7 +157,7 @@
    "id": "11",
    "metadata": {},
    "source": [
-    "## Download and transform satellite"
+    "## Quick and dirty plot"
    ]
   },
   {
@@ -146,6 +166,31 @@
    "id": "12",
    "metadata": {},
    "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(1, 1, figsize=(12, 5))\n",
+    "for rho in ds_reanalysis.values():\n",
+    "    da = compute_ssl(rho, rho0=1025)\n",
+    "    da.plot(label=da.long_name)\n",
+    "ax.set_ylabel(\"[m]\")\n",
+    "ax.legend()\n",
+    "ax.grid()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "13",
+   "metadata": {},
+   "source": [
+    "## Download and transform satellite"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "14",
+   "metadata": {},
+   "outputs": [],
    "source": [
     "ds_satellite = download.download_and_transform(*request_satellite, **download_kwargs)"
    ]