Added a notebook for Mass Uncertainty Project

notebooks/Uncertainties_Mass_Accounting.ipynb

@@ -0,0 +1,172 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0271f4a8-e7ff-4bc3-b60b-9f0f155c2df4",
+   "metadata": {},
+   "source": [
+    "# Uncertainties in Greenland mass loss due to uncertainties in accounting\n",
+    "\n",
+    "### Group Project International Summer School in Glaciology, McCarthy, June 7-17 2024\n",
+    "\n",
+    "## **How much mass is the Greenland Ice Sheet losing?**\n",
+    "\n",
+    "Published numbers vary not only due to the use of different data sets and methods, but also due to use of different ice sheet and basin outlines. Sometimes, glaciers and ice caps disconnected from the main ice sheet are included, sometimes they are not. This makes comparison challenging, and a community consensus is needed. In this project, we will explore different methods for basin delineation.\n",
+    "\n",
+    "<img src=\"gris_basins_tweet.png\" alt=\"Tweet by Hester Jiskoot on the plethora of Greenland basins\" style=\"width: 500px;\"/>"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "31d56b26-e60c-40d6-a86b-617a80252209",
+   "metadata": {},
+   "source": [
+    "## Method 1: Delineation using flowlines"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a65db713-9898-4bce-b451-253aeb4431d3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import xarray as xr\n",
+    "ds = xr.open_dataset(\"/mnt/storstrommen/data/its-live/GRE_G0240_0000.nc\", decode_times=False)\n",
+    "jak_ds = ds.sel(x=slice(-226_000, 120_000), y=slice(-2_120_000, -2_380_000))\n",
+    "jak_ds[[\"v\", \"vx\", \"vy\"]].to_netcdf(\"data/its_live_jak.nc\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2d972f4d-f427-4a13-b528-0e21be9575fb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import xarray as xr\n",
+    "import geopandas as gp\n",
+    "from shapely import Point\n",
+    "import pylab as plt\n",
+    "\n",
+    "from glacier_flow_tools.interpolation import velocity\n",
+    "from glacier_flow_tools.pathlines import (\n",
+    "    compute_pathline,\n",
+    "    series_to_pathline_geopandas_dataframe,\n",
+    ")\n",
+    "from glacier_flow_tools.utils import register_colormaps\n",
+    "register_colormaps()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7112896a-f5cb-48d2-8dae-72a94cb9a8f2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hmin = 1\n",
+    "hmax = 25.0\n",
+    "tol = 1\n",
+    "start_time = 0\n",
+    "end_time = 5_000\n",
+    "\n",
+    "ds = xr.open_dataset(\"data/its_live_jak.nc\")\n",
+    "Vx = np.squeeze(ds[\"vx\"].to_numpy())\n",
+    "Vy = np.squeeze(ds[\"vy\"].to_numpy())\n",
+    "\n",
+    "# Reverse for backward pathlines\n",
+    "Vx = -Vx\n",
+    "Vy = -Vy\n",
+    "\n",
+    "x = ds[\"x\"].to_numpy()\n",
+    "y = ds[\"y\"].to_numpy()\n",
+    "\n",
+    "geom = [Point(-188_000, -2_255_000), Point(-182_000, -2_290_000)]\n",
+    "starting_points_df = gp.GeoDataFrame(geometry=geom, crs=\"EPSG:3413\")\n",
+    "\n",
+    "pathlines = []\n",
+    "for _, df in starting_points_df.iterrows():\n",
+    "    pathline = compute_pathline(\n",
+    "        [*df.geometry.coords[0]],\n",
+    "        velocity,\n",
+    "        f_args=(Vx, Vy, x, y),\n",
+    "        hmin=hmin,\n",
+    "        hmax=hmax,\n",
+    "        tol=tol,\n",
+    "        start_time=start_time,\n",
+    "        end_time=end_time,\n",
+    "        notebook=True,\n",
+    "        progress=True,\n",
+    "    )\n",
+    "    pathlines.append(pathline)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "617aaa27-ad8d-4b37-ba04-ada0446919f9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Combine the pathlines into a GeoPandas DataFrame\n",
+    "\n",
+    "result = pd.concat(        [\n",
+    "    series_to_pathline_geopandas_dataframe(s.drop(\"geometry\", errors=\"ignore\"), pathlines[k])\n",
+    "    for k, s in starting_points_df.iterrows()\n",
+    "]\n",
+    "                  ).reset_index(drop=True)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7ca1f905-1e8e-4cb8-b783-f4b88fd99dc6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ratio = (y.max() - y.min()) / (x.max() - x.min())\n",
+    "fig, ax = plt.subplots(1, figsize=(12, 12 * ratio))\n",
+    "ds[\"v\"].plot(ax=ax, cmap=\"speed_colorblind\", vmin=10, vmax=1500)\n",
+    "result.plot(markersize=25,\n",
+    "            color=\"w\",\n",
+    "            legend=False,\n",
+    "            missing_kwds={},\n",
+    "            ax=ax,\n",
+    "            edgecolors=\"k\",\n",
+    "            linewidths=0.5,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "587ef9f4-5dd5-416b-b938-7dec1fdca2a9",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}