add p-value

bopen · May 30, 2024 · a08ad3b · a08ad3b
1 parent d9d31c9
commit a08ad3b
Showing 1 changed file with 68 additions and 12 deletions.
diff --git a/notebooks/wp5/glacier_mass_balance_maps.ipynb b/notebooks/wp5/glacier_mass_balance_maps.ipynb
@@ -28,6 +28,7 @@
     "import matplotlib.colors as mcolors\n",
     "import matplotlib.pyplot as plt\n",
     "import xarray as xr\n",
+    "import xskillscore as xs\n",
     "from c3s_eqc_automatic_quality_control import download\n",
     "\n",
     "plt.style.use(\"seaborn-v0_8-notebook\")"
@@ -97,7 +98,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def compute_sum_uncertainty_slope_and_acceleration(ds):\n",
+    "def compute_coeff_and_pvalue(cumsum, degree):\n",
+    "    coeff = cumsum.polyfit(\"time\", degree)\n",
+    "    (fit,) = xr.polyval(cumsum[\"time\"], coeff).data_vars.values()\n",
+    "    p_value = xs.pearson_r_p_value(cumsum.chunk(time=-1), fit.chunk(time=-1), \"time\")\n",
+    "    return coeff, p_value\n",
+    "\n",
+    "\n",
+    "def compute_time_statistics(ds):\n",
     "    with xr.set_options(keep_attrs=True):\n",
     "        cumsum = ds[\"Glacier\"].cumsum(\"time\").drop_vars(\"time\")\n",
     "        ds[\"Uncertainty\"] = (ds[\"Uncertainty\"] ** 2).sum(\"time\") ** (1 / 2)\n",
@@ -108,22 +116,35 @@
     "        da.attrs[\"long_name\"] = f\"Sum of {da.attrs['long_name']}\"\n",
     "\n",
     "    # Linear\n",
-    "    ds[\"Slope\"] = cumsum.polyfit(\"time\", 1)[\"polyfit_coefficients\"].sel(degree=1)\n",
+    "    coeff, p_value = compute_coeff_and_pvalue(cumsum, 1)\n",
+    "    ds[\"Slope\"] = coeff[\"polyfit_coefficients\"].sel(degree=1)\n",
     "    ds[\"Slope\"].attrs = {\n",
     "        \"long_name\": f\"Trend of {da.attrs['long_name']}\",\n",
     "        \"units\": f\"{da.attrs['units']} yr$^{-1}$\",\n",
     "    }\n",
+    "    ds[\"Pvalue1\"] = p_value\n",
+    "    ds[\"Pvalue1\"].attrs = {\n",
+    "        \"long_name\": f\"Linear p-value of {da.attrs['long_name']}\",\n",
+    "    }\n",
     "\n",
     "    # Quadratic\n",
-    "    ds[\"Acceleration\"] = 2 * cumsum.polyfit(\"time\", 2)[\"polyfit_coefficients\"].sel(\n",
-    "        degree=2\n",
-    "    )\n",
+    "    coeff, p_value = compute_coeff_and_pvalue(cumsum, 2)\n",
+    "    ds[\"Acceleration\"] = 2 * coeff[\"polyfit_coefficients\"].sel(degree=2)\n",
     "    ds[\"Acceleration\"].attrs = {\n",
     "        \"long_name\": f\"Quadratic Trend of {da.attrs['long_name']}\",\n",
     "        \"units\": f\"{da.attrs['units']} yr$^{-2}$\",\n",
     "    }\n",
+    "    ds[\"Pvalue2\"] = p_value\n",
+    "    ds[\"Pvalue2\"].attrs = {\n",
+    "        \"long_name\": f\"Quadratic p-value of {da.attrs['long_name']}\",\n",
+    "    }\n",
+    "\n",
+    "    return ds\n",
     "\n",
-    "    return ds"
+    "\n",
+    "def compute_spatial_statistics(ds):\n",
+    "    with xr.set_options(keep_attrs=True):\n",
+    "        return ds.sum((\"latitude\", \"longitude\"))"
    ]
   },
   {
@@ -141,16 +162,25 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "chunks = {\"hydrological_year\": 1}\n",
     "ds = download.download_and_transform(\n",
     "    collection_id,\n",
     "    request,\n",
-    "    chunks={\"hydrological_year\": 1},\n",
-    "    transform_func=compute_sum_uncertainty_slope_and_acceleration,\n",
+    "    chunks=chunks,\n",
+    "    transform_func=compute_time_statistics,\n",
     "    transform_chunks=False,\n",
     ")\n",
+    "ds_timeseries = download.download_and_transform(\n",
+    "    collection_id,\n",
+    "    request,\n",
+    "    chunks=chunks,\n",
+    "    transform_func=compute_spatial_statistics,\n",
+    ")\n",
+    "\n",
     "# Customize some attributes\n",
-    "for da in ds.data_vars.values():\n",
-    "    da.attrs[\"long_name\"] = da.attrs[\"long_name\"].replace(\"_\", \" \").title()"
+    "for obj in (ds, ds_timeseries):\n",
+    "    for da in obj.data_vars.values():\n",
+    "        da.attrs[\"long_name\"] = da.attrs[\"long_name\"].replace(\"_\", \" \").title()"
    ]
   },
   {
@@ -252,19 +282,45 @@
     "        \"label\": \"Mass change trend (Gt yr$^{-1}$)\",\n",
     "        \"title\": \"Linear trends (slopes) of the mass changes of the glaciers.\",\n",
     "    },\n",
+    "    \"Pvalue1\": {\n",
+    "        \"title\": \"Linear p-value.\",\n",
+    "    },\n",
     "    \"Acceleration\": {\n",
     "        \"cmap\": \"coolwarm_r\",\n",
     "        \"vmin\": -0.005,\n",
     "        \"vmax\": +0.005,\n",
     "        \"label\": \"Total mass change error (Gt yr$^{-2}$)\",\n",
     "        \"title\": \"Quadratic trends (accelerations) of the mass changes of the glaciers.\",\n",
     "    },\n",
+    "    \"Pvalue2\": {\n",
+    "        \"title\": \"Quadratics p-value.\",\n",
+    "    },\n",
     "}\n",
-    "\n",
     "for var_name, kwargs in column_kwargs.items():\n",
     "    ax = plot_map(gdf, var_name, **kwargs)\n",
     "    plt.show()"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "17",
+   "metadata": {},
+   "source": [
+    "## Plot timeseries"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "18",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for da in ds_timeseries.data_vars.values():\n",
+    "    da.plot()\n",
+    "    plt.grid()\n",
+    "    plt.show()"
+   ]
   }
  ],
  "metadata": {
@@ -283,7 +339,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.8"
+   "version": "3.11.9"
   }
  },
  "nbformat": 4,