regionalise ocean color

bopen · Jul 24, 2023 · f0d0ff3 · f0d0ff3
1 parent 240f47c
commit f0d0ff3
Showing 1 changed file with 36 additions and 23 deletions.
diff --git a/notebooks/wp5/ocean_color.ipynb b/notebooks/wp5/ocean_color.ipynb
@@ -26,7 +26,7 @@
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\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, utils\n",
     "\n",
     "plt.style.use(\"seaborn-v0_8-notebook\")"
    ]
@@ -48,7 +48,7 @@
    "source": [
     "# Time period\n",
     "start = \"1997-09\"\n",
-    "stop = \"2023-03\""
+    "stop = \"1998-08\""
    ]
   },
   {
@@ -96,11 +96,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def masked_global_mean(ds):\n",
+    "def regionalised_spatial_weighted_mean(ds, lon_slice, lat_slice):\n",
     "    da = ds[\"chlor_a\"]\n",
+    "    da = utils.regionalise(da, lon_slice=lon_slice, lat_slice=lat_slice)\n",
     "    da = da.where((da > 1.0e-3) & (da < 1.0e2))\n",
     "    with xr.set_options(keep_attrs=True):\n",
-    "        da = 10 ** diagnostics.spatial_weighted_mean(np.log10(da), weights=False)\n",
+    "        da = 10 ** diagnostics.spatial_weighted_mean(np.log10(da))\n",
     "    da.attrs[\"long_name\"] = da.attrs[\"long_name\"].replace(\" (not log-transformed)\", \"\")\n",
     "    return da.to_dataset()"
    ]
@@ -120,12 +121,19 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ds = download.download_and_transform(\n",
-    "    collection_id,\n",
-    "    requests,\n",
-    "    transform_func=masked_global_mean,\n",
-    "    chunks={\"year\": 1, \"month\": 1},\n",
-    ")\n",
+    "lon_slice = slice(-180, 180)\n",
+    "datasets = []\n",
+    "for lats in zip(range(90, -90, -30), range(60, -120, -30)):\n",
+    "    print(f\"{lats=}\")\n",
+    "    ds = download.download_and_transform(\n",
+    "        collection_id,\n",
+    "        requests,\n",
+    "        transform_func=regionalised_spatial_weighted_mean,\n",
+    "        transform_func_kwargs={\"lon_slice\": lon_slice, \"lat_slice\": slice(*lats)},\n",
+    "        chunks={\"year\": 1, \"month\": 1},\n",
+    "    )\n",
+    "    datasets.append(ds.expand_dims(latitudes=[str(lats)]))\n",
+    "ds = xr.concat(datasets, \"latitudes\")\n",
     "da_daily = ds[\"chlor_a\"]\n",
     "da_monthly = da_daily.resample(time=\"M\").mean()\n",
     "da_seasonality = da_daily.groupby(\"time.month\").mean()"
@@ -147,10 +155,13 @@
    "outputs": [],
    "source": [
     "window = 60\n",
-    "da_daily.plot(label=\"daily\")\n",
-    "da_daily.rolling(time=window, min_periods=1).mean().plot(label=f\"{window}-day running\")\n",
+    "for i, (lats, da) in enumerate(da_daily.groupby(\"latitudes\")):\n",
+    "    da.plot(color=f\"C{i}\", ls=\"-\", label=f\"{lats}: daily\")\n",
+    "    da.rolling(time=window, min_periods=1).mean().plot(\n",
+    "        color=f\"C{i}\", ls=\"--\", label=f\"{lats}: {window}-day running\"\n",
+    "    )\n",
     "plt.grid()\n",
-    "_ = plt.legend()"
+    "_ = plt.legend(bbox_to_anchor=(1, 1))"
    ]
   },
   {
@@ -169,12 +180,13 @@
    "outputs": [],
    "source": [
     "window = 6\n",
-    "da_monthly.plot(label=\"monthly\")\n",
-    "da_monthly.rolling(time=window, min_periods=1).mean().plot(\n",
-    "    label=f\"{window}-month running\"\n",
-    ")\n",
+    "for i, (lats, da) in enumerate(da_monthly.groupby(\"latitudes\")):\n",
+    "    da.plot(color=f\"C{i}\", ls=\"-\", label=f\"{lats}: monthly\")\n",
+    "    da.rolling(time=window, min_periods=1).mean().plot(\n",
+    "        color=f\"C{i}\", ls=\"--\", label=f\"{lats}: {window}-month running\"\n",
+    "    )\n",
     "plt.grid()\n",
-    "_ = plt.legend()"
+    "_ = plt.legend(bbox_to_anchor=(1, 1))"
    ]
   },
   {
@@ -193,13 +205,14 @@
    "outputs": [],
    "source": [
     "window = 6\n",
-    "da_seasonality.plot(label=\"seasonality\", marker=\"o\")\n",
-    "da_seasonality.rolling(month=window, min_periods=1).mean().plot(\n",
-    "    label=f\"{window}-month running\", marker=\"o\"\n",
-    ")\n",
+    "for i, (lats, da) in enumerate(da_seasonality.groupby(\"latitudes\")):\n",
+    "    da.plot(color=f\"C{i}\", ls=\"-\", label=f\"{lats}: seasonality\", marker=\"o\")\n",
+    "    da.rolling(month=window, min_periods=1).mean().plot(\n",
+    "        color=f\"C{i}\", ls=\"--\", label=f\"{lats}: {window}-month running\", marker=\"o\"\n",
+    "    )\n",
     "plt.title(f\"From {start} to {stop}\")\n",
     "plt.grid()\n",
-    "_ = plt.legend()"
+    "_ = plt.legend(bbox_to_anchor=(1, 1))"
    ]
   }
  ],