🔀 Merge branch 'ice_volume_displacement' (#216)

Closes #216 Calculate rolling time-series of ice volume displacement over subglacial lakes.

Author: weiji14

.github/workflows/python-app.yml

@@ -69,4 +69,4 @@ jobs:
         run: poetry run black . --check
 
       - name: Test with pytest
-        run: poetry run pytest --verbose deepicedrain/
+        run: poetry run pytest --verbose --doctest-modules deepicedrain/

README.md

@@ -18,7 +18,7 @@ in Antarctica using remote sensing and machine learning.
 
 | Along track view of an ATL11 Ground Track | Elevation time-series at Crossover Points |
 |---|---|
-| ![alongtrack_whillans_ix_1080_pt3](https://user-images.githubusercontent.com/23487320/102156291-2210f600-3ee2-11eb-8175-e854b70444df.png) | ![crossover_many_normalized_whillans_ix_2018-10-14_2020-09-30](https://user-images.githubusercontent.com/23487320/102156396-5dabc000-3ee2-11eb-9e42-37463f81058d.png) |
+| ![alongtrack_whillans_ix_1080_pt3](https://user-images.githubusercontent.com/23487320/102156291-2210f600-3ee2-11eb-8175-e854b70444df.png) | ![crossover_anomaly_whillans_ix_2018-10-14_2020-09-30](https://user-images.githubusercontent.com/23487320/102610765-9bcf0b00-4192-11eb-803b-247ed960f9bb.png) |
 
 
 

atlxi_xover.ipynb

@@ -36,9 +36,12 @@
     "import geopandas as gpd\n",
     "import numpy as np\n",
     "import pandas as pd\n",
+    "import pint\n",
+    "import pint_pandas\n",
     "import pygmt\n",
     "import shapely.geometry\n",
-    "import tqdm"
+    "import tqdm\n",
+    "import uncertainties"
    ]
   },
   {
@@ -49,6 +52,8 @@
    },
    "outputs": [],
    "source": [
+    "ureg = pint.UnitRegistry()\n",
+    "pint_pandas.PintType.ureg = ureg\n",
     "tag: str = \"X2SYS\"\n",
     "os.environ[\"X2SYS_HOME\"] = os.path.abspath(tag)\n",
     "client = dask.distributed.Client(\n",
@@ -111,7 +116,7 @@
    "outputs": [],
    "source": [
     "# Choose one Antarctic active subglacial lake polygon with EPSG:3031 coordinates\n",
-    "lake_name: str = \"Subglacial Lake Conway\"\n",
+    "lake_name: str = \"Whillans 12\"\n",
     "lake_catalog = deepicedrain.catalog.subglacial_lakes()\n",
     "lake_ids: list = (\n",
     "    pd.json_normalize(lake_catalog.metadata[\"lakedict\"])\n",
@@ -142,7 +147,11 @@
    "source": [
     "# Subset data to lake of interest\n",
     "placename: str = region.name.lower().replace(\" \", \"_\")\n",
-    "df_lake: pd.DataFrame = region.subset(data=df_dhdt)"
+    "df_lake: cudf.DataFrame = region.subset(data=df_dhdt)  # bbox subset\n",
+    "gdf_lake = gpd.GeoDataFrame(\n",
+    "    df_lake, geometry=gpd.points_from_xy(x=df_lake.x, y=df_lake.y, crs=3031)\n",
+    ")\n",
+    "df_lake: pd.DataFrame = df_lake.loc[gdf_lake.within(lake.geometry)]  # polygon subset"
    ]
   },
   {
@@ -251,9 +260,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
-    "lines_to_next_cell": 2
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "# 2D plot of crossover locations\n",
@@ -318,7 +325,8 @@
     "    stubnames=[\"t\", \"h\"],\n",
     "    j=\"track\",\n",
     ")\n",
-    "df_th = df_th.drop_duplicates(ignore_index=True)"
+    "df_th: pd.DataFrame = df_th.drop_duplicates(ignore_index=True)\n",
+    "df_th: pd.DataFrame = df_th.sort_values(by=\"t\").reset_index(drop=True)"
    ]
   },
   {
@@ -328,9 +336,9 @@
    "outputs": [],
    "source": [
     "# Plot at single location with **maximum** absolute crossover height error (max_h_X)\n",
-    "df_max = df_th.query(expr=\"x == @max_h_X.x & y == @max_h_X.y\").sort_values(by=\"t\")\n",
+    "df_max = df_th.query(expr=\"x == @max_h_X.x & y == @max_h_X.y\")\n",
     "track1, track2 = df_max.track1_track2.iloc[0].split(\"x\")\n",
-    "print(f\"{round(max_h_X.h_X, 2)} metres height change at {max_h_X.x}, {max_h_X.y}\")\n",
+    "print(f\"{max_h_X.h_X:.2f} metres height change at {max_h_X.x}, {max_h_X.y}\")\n",
     "plotregion = np.array(\n",
     "    [df_max.t.min(), df_max.t.max(), *pygmt.info(table=df_max[[\"h\"]], spacing=2.5)[:2]]\n",
     ")\n",
@@ -384,20 +392,95 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Plot all crossover height points over time over the lake area\n",
-    "# with height values normalized to 0 from the first observation date\n",
-    "normfunc = lambda h: h - h.iloc[0]  # lambda h: h - h.mean()\n",
-    "df_th[\"h_norm\"] = df_th.groupby(by=\"track1_track2\").h.transform(func=normfunc)\n",
-    "\n",
+    "# Calculate height anomaly at crossover point as\n",
+    "# height at t=n minus height at t=0 (first observation date at crossover point)\n",
+    "anomfunc = lambda h: h - h.iloc[0]  # lambda h: h - h.mean()\n",
+    "df_th[\"h_anom\"] = df_th.groupby(by=\"track1_track2\").h.transform(func=anomfunc)\n",
+    "# Calculate ice volume displacement (dvol) in unit metres^3\n",
+    "# and rolling mean height anomaly (h_roll) in unit metres\n",
+    "surface_area: pint.Quantity = lake.geometry.area * ureg.metre ** 2\n",
+    "ice_dvol: pd.Series = deepicedrain.ice_volume_over_time(\n",
+    "    df_elev=df_th.astype(dtype={\"h_anom\": \"pint[metre]\"}),\n",
+    "    surface_area=surface_area,\n",
+    "    time_col=\"t\",\n",
+    "    outfile=f\"figures/{placename}/ice_dvol_dt_{placename}.txt\",\n",
+    ")\n",
+    "df_th[\"h_roll\"]: pd.Series = uncertainties.unumpy.nominal_values(\n",
+    "    arr=ice_dvol.pint.magnitude / surface_area.magnitude\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot all crossover height point anomalies over time over the lake area\n",
     "fig = deepicedrain.plot_crossovers(\n",
     "    df=df_th,\n",
     "    regionname=region.name,\n",
-    "    elev_var=\"h_norm\",\n",
+    "    elev_var=\"h_anom\",\n",
     "    outline_points=f\"figures/{placename}/{placename}.gmt\",\n",
     ")\n",
+    "fig.plot(x=df_th.t, y=df_th.h_roll, pen=\"thick,-\")  # plot rolling mean height anomaly\n",
     "fig.savefig(\n",
-    "    f\"figures/{placename}/crossover_many_normalized_{placename}_{min_date}_{max_date}.png\"\n",
+    "    f\"figures/{placename}/crossover_anomaly_{placename}_{min_date}_{max_date}.png\"\n",
+    ")\n",
+    "fig.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Combined ice volume displacement plot\n",
+    "\n",
+    "Showing how subglacial water cascades down a drainage basin!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "lines_to_next_cell": 0
+   },
+   "outputs": [],
+   "source": [
+    "fig = pygmt.Figure()\n",
+    "fig.basemap(\n",
+    "    region=f\"2019-02-28/2020-09-30/-0.3/0.5\",\n",
+    "    frame=[\"wSnE\", \"xaf\", 'yaf+l\"Ice Volume Displacement (km@+3@+)\"'],\n",
+    ")\n",
+    "pygmt.makecpt(cmap=\"davosS\", color_model=\"+c\", series=(-2, 4, 0.5))\n",
+    "for i, (_placename, linestyle) in enumerate(\n",
+    "    iterable=zip(\n",
+    "        [\"whillans_ix\", \"subglacial_lake_whillans\", \"whillans_12\", \"whillans_7\"],\n",
+    "        [\"\", \".-\", \"-\", \"..-\"],\n",
+    "    )\n",
+    "):\n",
+    "    fig.plot(\n",
+    "        data=f\"figures/{_placename}/ice_dvol_dt_{_placename}.txt\",\n",
+    "        cmap=True,\n",
+    "        pen=f\"thick,{linestyle}\",\n",
+    "        zvalue=i,\n",
+    "        label=_placename,\n",
+    "        columns=\"0,3\",  # time column (0), ice_dvol column (3)\n",
+    "    )\n",
+    "fig.text(\n",
+    "    position=\"TL\",\n",
+    "    offset=\"j0.2c\",\n",
+    "    text=\"Whillans Ice Stream Central Catchment active subglacial lakes\",\n",
     ")\n",
+    "fig.legend(position=\"jML+jML+o0.2c\", box=\"+gwhite\")\n",
+    "fig.savefig(\"figures/cascade_whillans_ice_stream.png\")\n",
     "fig.show()"
    ]
   },

atlxi_xover.py

@@ -40,11 +40,16 @@
 import geopandas as gpd
 import numpy as np
 import pandas as pd
+import pint
+import pint_pandas
 import pygmt
 import shapely.geometry
 import tqdm
+import uncertainties
 
 # %%
+ureg = pint.UnitRegistry()
+pint_pandas.PintType.ureg = ureg
 tag: str = "X2SYS"
 os.environ["X2SYS_HOME"] = os.path.abspath(tag)
 client = dask.distributed.Client(
@@ -79,7 +84,7 @@
 
 # %%
 # Choose one Antarctic active subglacial lake polygon with EPSG:3031 coordinates
-lake_name: str = "Subglacial Lake Conway"
+lake_name: str = "Whillans 12"
 lake_catalog = deepicedrain.catalog.subglacial_lakes()
 lake_ids: list = (
     pd.json_normalize(lake_catalog.metadata["lakedict"])
@@ -102,7 +107,11 @@
 # %%
 # Subset data to lake of interest
 placename: str = region.name.lower().replace(" ", "_")
-df_lake: pd.DataFrame = region.subset(data=df_dhdt)
+df_lake: cudf.DataFrame = region.subset(data=df_dhdt)  # bbox subset
+gdf_lake = gpd.GeoDataFrame(
+    df_lake, geometry=gpd.points_from_xy(x=df_lake.x, y=df_lake.y, crs=3031)
+)
+df_lake: pd.DataFrame = df_lake.loc[gdf_lake.within(lake.geometry)]  # polygon subset
 
 
 # %%
@@ -214,7 +223,6 @@
 fig.savefig(f"figures/{placename}/crossover_area_{placename}_{min_date}_{max_date}.png")
 fig.show()
 
-
 # %% [markdown]
 # ### Plot Crossover Elevation time-series
 #
@@ -231,13 +239,14 @@
     stubnames=["t", "h"],
     j="track",
 )
-df_th = df_th.drop_duplicates(ignore_index=True)
+df_th: pd.DataFrame = df_th.drop_duplicates(ignore_index=True)
+df_th: pd.DataFrame = df_th.sort_values(by="t").reset_index(drop=True)
 
 # %%
 # Plot at single location with **maximum** absolute crossover height error (max_h_X)
-df_max = df_th.query(expr="x == @max_h_X.x & y == @max_h_X.y").sort_values(by="t")
+df_max = df_th.query(expr="x == @max_h_X.x & y == @max_h_X.y")
 track1, track2 = df_max.track1_track2.iloc[0].split("x")
-print(f"{round(max_h_X.h_X, 2)} metres height change at {max_h_X.x}, {max_h_X.y}")
+print(f"{max_h_X.h_X:.2f} metres height change at {max_h_X.x}, {max_h_X.y}")
 plotregion = np.array(
     [df_max.t.min(), df_max.t.max(), *pygmt.info(table=df_max[["h"]], spacing=2.5)[:2]]
 )
@@ -279,20 +288,72 @@
 fig.show()
 
 # %%
-# Plot all crossover height points over time over the lake area
-# with height values normalized to 0 from the first observation date
-normfunc = lambda h: h - h.iloc[0]  # lambda h: h - h.mean()
-df_th["h_norm"] = df_th.groupby(by="track1_track2").h.transform(func=normfunc)
+# Calculate height anomaly at crossover point as
+# height at t=n minus height at t=0 (first observation date at crossover point)
+anomfunc = lambda h: h - h.iloc[0]  # lambda h: h - h.mean()
+df_th["h_anom"] = df_th.groupby(by="track1_track2").h.transform(func=anomfunc)
+# Calculate ice volume displacement (dvol) in unit metres^3
+# and rolling mean height anomaly (h_roll) in unit metres
+surface_area: pint.Quantity = lake.geometry.area * ureg.metre ** 2
+ice_dvol: pd.Series = deepicedrain.ice_volume_over_time(
+    df_elev=df_th.astype(dtype={"h_anom": "pint[metre]"}),
+    surface_area=surface_area,
+    time_col="t",
+    outfile=f"figures/{placename}/ice_dvol_dt_{placename}.txt",
+)
+df_th["h_roll"]: pd.Series = uncertainties.unumpy.nominal_values(
+    arr=ice_dvol.pint.magnitude / surface_area.magnitude
+)
 
+# %%
+# Plot all crossover height point anomalies over time over the lake area
 fig = deepicedrain.plot_crossovers(
     df=df_th,
     regionname=region.name,
-    elev_var="h_norm",
+    elev_var="h_anom",
     outline_points=f"figures/{placename}/{placename}.gmt",
 )
+fig.plot(x=df_th.t, y=df_th.h_roll, pen="thick,-")  # plot rolling mean height anomaly
 fig.savefig(
-    f"figures/{placename}/crossover_many_normalized_{placename}_{min_date}_{max_date}.png"
+    f"figures/{placename}/crossover_anomaly_{placename}_{min_date}_{max_date}.png"
+)
+fig.show()
+
+# %%
+
+# %% [markdown]
+# ## Combined ice volume displacement plot
+#
+# Showing how subglacial water cascades down a drainage basin!
+
+# %%
+fig = pygmt.Figure()
+fig.basemap(
+    region=f"2019-02-28/2020-09-30/-0.3/0.5",
+    frame=["wSnE", "xaf", 'yaf+l"Ice Volume Displacement (km@+3@+)"'],
+)
+pygmt.makecpt(cmap="davosS", color_model="+c", series=(-2, 4, 0.5))
+for i, (_placename, linestyle) in enumerate(
+    iterable=zip(
+        ["whillans_ix", "subglacial_lake_whillans", "whillans_12", "whillans_7"],
+        ["", ".-", "-", "..-"],
+    )
+):
+    fig.plot(
+        data=f"figures/{_placename}/ice_dvol_dt_{_placename}.txt",
+        cmap=True,
+        pen=f"thick,{linestyle}",
+        zvalue=i,
+        label=_placename,
+        columns="0,3",  # time column (0), ice_dvol column (3)
+    )
+fig.text(
+    position="TL",
+    offset="j0.2c",
+    text="Whillans Ice Stream Central Catchment active subglacial lakes",
 )
+fig.legend(position="jML+jML+o0.2c", box="+gwhite")
+fig.savefig("figures/cascade_whillans_ice_stream.png")
 fig.show()
 
 # %%

deepicedrain/README.md

@@ -29,6 +29,7 @@ Contents:
 
 - :droplet: lakealgorithms.py - Custom algorithms for detecting and filtering active subglacial lakes
   - find_clusters - Density based clustering algorithm (DBSCAN) to group points into lakes
+  - ice_volume_over_time - Calculate a rolling mean time-series of ice volume displacement over a lake area
 
 - :world_map: vizplots.py - Makes interactive dashboard plots and publication quality figures
   - IceSat2Explorer - Dashboard for interacting with ICESat-2 point clouds on a 2D map

deepicedrain/__init__.py

@@ -11,7 +11,7 @@
     split_tracks,
     wide_to_long,
 )
-from deepicedrain.lake_algorithms import find_clusters
+from deepicedrain.lake_algorithms import find_clusters, ice_volume_over_time
 from deepicedrain.spatiotemporal import (
     Region,
     deltatime_to_utctime,