🔀 Merge branch 'transect_over_3d_surface' (#227)

Closes #227 3D ice surface elevation and trend plot with transect crossing.

atlxi_lake.ipynb

@@ -951,12 +951,12 @@
    ],
    "source": [
     "# Choose one Antarctic active subglacial lake polygon with EPSG:3031 coordinates\n",
-    "lake_name: str = \"Subglacial Lake Conway\"\n",
+    "lake_name: str = \"Whillans 7\"\n",
     "lake_catalog = deepicedrain.catalog.subglacial_lakes()\n",
-    "lake_ids: list = (\n",
+    "lake_ids, transect_id = (\n",
     "    pd.json_normalize(lake_catalog.metadata[\"lakedict\"])\n",
-    "    .query(\"lakename == @lake_name\")\n",
-    "    .ids.iloc[0]\n",
+    "    .query(\"lakename == @lake_name\")[[\"ids\", \"transect\"]]\n",
+    "    .iloc[0]\n",
     ")\n",
     "lake = (\n",
     "    lake_catalog.read()\n",
@@ -983,6 +983,24 @@
     "# Subset data to lake of interest\n",
     "placename: str = region.name.lower().replace(\" \", \"_\")\n",
     "df_lake: cudf.DataFrame = region.subset(data=df_dhdt)\n",
+    "# Get all raw xyz points and one transect line dataframe\n",
+    "track_dict: dict = deepicedrain.split_tracks(df=df_lake.to_pandas())\n",
+    "track_points: pd.DataFrame = (\n",
+    "    pd.concat(track_dict.values())\n",
+    "    .groupby(by=[\"x\", \"y\"])\n",
+    "    .mean()  # z value is mean h_corr over all cycles\n",
+    "    .reset_index()[[\"x\", \"y\", \"h_corr\"]]\n",
+    ")\n",
+    "try:\n",
+    "    _rgt, _pt = transect_id.split(\"_\")\n",
+    "    df_transect: pd.DataFrame = (\n",
+    "        track_dict[transect_id][[\"x\", \"y\", \"h_corr\", \"cycle_number\"]]\n",
+    "        .groupby(by=[\"x\", \"y\"])\n",
+    "        .max()  # z value is maximum h_corr over all cycles\n",
+    "        .reset_index()\n",
+    "    )\n",
+    "except AttributeError:\n",
+    "    pass\n",
     "\n",
     "# Save lake outline to OGR GMT file format\n",
     "outline_points: str = f\"figures/{placename}/{placename}.gmt\"\n",
@@ -1056,7 +1074,38 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# 3D plots of gridded ice surface elevation over time\n",
+    "# 3D plot of mean ice surface elevation (<z>) and rate of ice elevation change (dhdt)\n",
+    "fig = deepicedrain.plot_icesurface(\n",
+    "    grid=f\"figures/{placename}/xyht_{placename}.nc?z_mean\",\n",
+    "    grid_region=grid_region,\n",
+    "    diff_grid=f\"figures/{placename}/xyht_{placename}.nc?dhdt\",\n",
+    "    track_points=track_points.to_numpy(),\n",
+    "    outline_points=outline_points,\n",
+    "    azimuth=157.5,  # 202.5  # 270\n",
+    "    elevation=45,  # 60\n",
+    "    title=f\"{region.name} Ice Surface\",\n",
+    ")\n",
+    "# Plot crossing transect line\n",
+    "fig.plot3d(\n",
+    "    data=df_transect[[\"x\", \"y\", \"h_corr\"]].to_numpy(),\n",
+    "    color=\"yellow2\",\n",
+    "    style=\"c0.1c\",\n",
+    "    zscale=True,\n",
+    "    perspective=True,\n",
+    ")\n",
+    "fig.savefig(f\"figures/{placename}/dsm_{placename}_cycles_{cycles[0]}-{cycles[-1]}.png\")\n",
+    "fig.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "lines_to_next_cell": 2
+   },
+   "outputs": [],
+   "source": [
+    "# 3D plots of gridded ice surface elevation over time (one per cycle)\n",
     "for cycle in tqdm.tqdm(iterable=cycles):\n",
     "    time_nsec: pd.Timestamp = df_lake[f\"utc_time_{cycle}\"].to_pandas().mean()\n",
     "    time_sec: str = np.datetime_as_string(arr=time_nsec.to_datetime64(), unit=\"s\")\n",

atlxi_lake.py

@@ -330,12 +330,12 @@
 
 # %%
 # Choose one Antarctic active subglacial lake polygon with EPSG:3031 coordinates
-lake_name: str = "Subglacial Lake Conway"
+lake_name: str = "Lake 12"
 lake_catalog = deepicedrain.catalog.subglacial_lakes()
-lake_ids: list = (
+lake_ids, transect_id = (
     pd.json_normalize(lake_catalog.metadata["lakedict"])
-    .query("lakename == @lake_name")
-    .ids.iloc[0]
+    .query("lakename == @lake_name")[["ids", "transect"]]
+    .iloc[0]
 )
 lake = (
     lake_catalog.read()
@@ -354,6 +354,24 @@
 # Subset data to lake of interest
 placename: str = region.name.lower().replace(" ", "_")
 df_lake: cudf.DataFrame = region.subset(data=df_dhdt)
+# Get all raw xyz points and one transect line dataframe
+track_dict: dict = deepicedrain.split_tracks(df=df_lake.to_pandas())
+track_points: pd.DataFrame = (
+    pd.concat(track_dict.values())
+    .groupby(by=["x", "y"])
+    .mean()  # z value is mean h_corr over all cycles
+    .reset_index()[["x", "y", "h_corr"]]
+)
+try:
+    _rgt, _pt = transect_id.split("_")
+    df_transect: pd.DataFrame = (
+        track_dict[transect_id][["x", "y", "h_corr", "cycle_number"]]
+        .groupby(by=["x", "y"])
+        .max()  # z value is maximum h_corr over all cycles
+        .reset_index()
+    )
+except AttributeError:
+    pass
 
 # Save lake outline to OGR GMT file format
 outline_points: str = f"figures/{placename}/{placename}.gmt"
@@ -390,7 +408,30 @@
 print(f"Elevation limits are: {z_limits}")
 
 # %%
-# 3D plots of gridded ice surface elevation over time
+# 3D plot of mean ice surface elevation (<z>) and rate of ice elevation change (dhdt)
+fig = deepicedrain.plot_icesurface(
+    grid=f"figures/{placename}/xyht_{placename}.nc?z_mean",
+    grid_region=grid_region,
+    diff_grid=f"figures/{placename}/xyht_{placename}.nc?dhdt",
+    track_points=track_points.to_numpy(),
+    outline_points=outline_points,
+    azimuth=157.5,  # 202.5  # 270
+    elevation=45,  # 60
+    title=f"{region.name} Ice Surface",
+)
+# Plot crossing transect line
+fig.plot3d(
+    data=df_transect[["x", "y", "h_corr"]].to_numpy(),
+    color="yellow2",
+    style="c0.1c",
+    zscale=True,
+    perspective=True,
+)
+fig.savefig(f"figures/{placename}/dsm_{placename}_cycles_{cycles[0]}-{cycles[-1]}.png")
+fig.show()
+
+# %%
+# 3D plots of gridded ice surface elevation over time (one per cycle)
 for cycle in tqdm.tqdm(iterable=cycles):
     time_nsec: pd.Timestamp = df_lake[f"utc_time_{cycle}"].to_pandas().mean()
     time_sec: str = np.datetime_as_string(arr=time_nsec.to_datetime64(), unit="s")
@@ -423,7 +464,7 @@
         "120",
         "-loop",
         "0",
-        f"figures/{placename}/dsm_*.png",
+        f"figures/{placename}/dsm_*cycle_*.png",
         gif_fname,
     ]
 )

atlxi_xover.ipynb

@@ -116,7 +116,7 @@
    "outputs": [],
    "source": [
     "# Choose one Antarctic active subglacial lake polygon with EPSG:3031 coordinates\n",
-    "lake_name: str = \"Whillans 12\"\n",
+    "lake_name: str = \"Lake 12\"\n",
     "lake_catalog = deepicedrain.catalog.subglacial_lakes()\n",
     "lake_ids: list = (\n",
     "    pd.json_normalize(lake_catalog.metadata[\"lakedict\"])\n",
@@ -462,7 +462,7 @@
     "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",
+    "        [\"whillans_ix\", \"subglacial_lake_whillans\", \"lake_12\", \"whillans_7\"],\n",
     "        [\"\", \".-\", \"-\", \"..-\"],\n",
     "    )\n",
     "):\n",

atlxi_xover.py

@@ -84,7 +84,7 @@
 
 # %%
 # Choose one Antarctic active subglacial lake polygon with EPSG:3031 coordinates
-lake_name: str = "Whillans 12"
+lake_name: str = "Lake 12"
 lake_catalog = deepicedrain.catalog.subglacial_lakes()
 lake_ids: list = (
     pd.json_normalize(lake_catalog.metadata["lakedict"])
@@ -335,7 +335,7 @@
 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"],
+        ["whillans_ix", "subglacial_lake_whillans", "lake_12", "whillans_7"],
         ["", ".-", "-", "..-"],
     )
 ):

deepicedrain/atlas_catalog.yaml

@@ -201,10 +201,14 @@ sources:
           ids: [62]
         - lakename: Kamb 34
           ids: [63]
+        - lakename: Kamb 6
+          ids: [60]
         - lakename: Kamb 8
           ids: [61]
         - lakename: MacAyeal 1
           ids: [84]
+        - lakename: MacAyeal 4
+          ids: [83]
         - lakename: Subglacial Lake Mercer
           ids: [15, 19]
         - lakename: Recovery 2
@@ -215,22 +219,36 @@ sources:
           ids: [101]
         - lakename: Slessor 45
           ids: [95]
-        - lakename: Lake 78
-          ids: [16, 46, 48]
         - lakename: Subglacial Lake Conway
           ids: [34, 35]
         - lakename: Subglacial Lake Whillans
           ids: [41, 43, 45]
-        - lakename: Whillans 12
-          ids: [54]
-        - lakename: Whillans 13
-          ids: [39, 40]
+          transect: 0989_pt1
         - lakename: Whillans 6
           ids: [33]
         - lakename: Whillans 7
           ids: [32]
+          transect: 0531_pt1
         - lakename: Whillans IX
           ids: [44]
+          transect: 1080_pt3
+        - lakename: Subglacial Lake Engelhardt 2
+          ids: [42]
+        - lakename: Lake 78
+          ids: [16, 46, 48]
+        - lakename: Lake 12
+          ids: [54]
+          transect: 0593_pt1
+        - lakename: \* 1
+          ids: [36]
+        - lakename: \* 2
+          ids: [39, 40]
+        - lakename: \* 3
+          ids: [56]
+        - lakename: \* 4
+          ids: [69]
+        - lakename: \* 5
+          ids: [130]
   test_data:
     description: 'Sample ICESat-2 datasets for testing purposes'
     args:

deepicedrain/features/subglacial_lakes.feature

@@ -28,7 +28,7 @@ Feature: Mapping Antarctic subglacial lakes
     # | Whillans X   | whillans_upstream   | 3-8    | 157.5   | 45        |
     # | Whillans XI  | whillans_upstream   | 3-8    | 157.5   | 45        |
     # | Whillans IX  | whillans_upstream   | 3-8    | 157.5   | 45        |
-    # | Whillans 12  | whillans_downstream | 3-8    | 157.5   | 45        |
+    # | Lake 12      | whillans_downstream | 3-8    | 157.5   | 45        |
     # | Kamb 8       | whillans_upstream   | 3-8    | 157.5   | 45        |
     # | Kamb 1       | whillans_upstream   | 3-8    | 157.5   | 45        |
     | Kamb 34      | whillans_upstream   | 4-7    | 157.5   | 45        |
@@ -63,4 +63,4 @@ Feature: Mapping Antarctic subglacial lakes
     | Whillans 7               | whillans_upstream   |
     # | Whillans IX              | whillans_upstream   |
     # | Subglacial Lake Whillans | whillans_downstream |
-    | Whillans 12              | whillans_downstream |
+    | Lake 12                  | whillans_downstream |

deepicedrain/spatiotemporal.py

@@ -313,6 +313,7 @@ def spatiotemporal_cube(
     x_var: str = "x",
     y_var: str = "y",
     z_var: str = "h_corr",
+    dhdt_var: str = "dhdt_slope",
     spacing: int = 250,
     clip_limits: bool = True,
     cycles: list = None,
@@ -376,8 +377,8 @@ def spatiotemporal_cube(
         +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84
         +units=m +no_defs', i.e. Antarctic Polar Stereographic EPSG:3031.
     folder : str
-        The folder to keep the intermediate NetCDF file in. Default is to place
-        the files in the current working directory.
+        The folder to keep the intermediate NetCDF files in. Default is to
+        place the files in the current working directory.
 
     Returns
     -------
@@ -413,6 +414,14 @@ def spatiotemporal_cube(
 
     # Create one grid surface for each time cycle
     _placename = f"_{placename}" if placename else ""
+    surface_kwargs = dict(
+        region=grid_region,
+        spacing=spacing,
+        J=f'"{projection}"',  # projection
+        M="3c",  # mask values 3 pixel cells outside/away from valid data
+        T=0.35,  # tension factor
+        V="e",  # error messages only
+    )
     for cycle in tqdm.tqdm(iterable=cycles):
         df_trimmed = pygmt.blockmedian(
             table=table[[x_var, y_var, f"{z_var}_{cycle}"]].dropna(),
@@ -421,17 +430,8 @@ def spatiotemporal_cube(
         )
         outfile = f"{z_var}{_placename}_cycle_{cycle}.nc"
         pygmt.surface(
-            data=df_trimmed.values,
-            region=grid_region,
-            spacing=spacing,
-            J=f'"{projection}"',  # projection
-            L=limits,  # lower and upper limits
-            M="3c",  # mask values 3 pixel cells outside/away from valid data
-            T=0.35,  # tension factor
-            V="e",  # error messages only
-            outfile=outfile,
+            data=df_trimmed.values, outfile=outfile, L=limits, **surface_kwargs
         )
-        # print(pygmt.grdinfo(outfile))
 
     # Move files into new folder if requested
     paths: list = [f"{z_var}{_placename}_cycle_{cycle}.nc" for cycle in cycles]
@@ -448,4 +448,30 @@ def spatiotemporal_cube(
         attrs_file=paths[-1],
     )
 
+    # Extra data variables, calculated using point data and then interpolated to grid
+    if dhdt_var in table.columns:
+        # Mean elevation (<z>)
+        df_zmean: pd.DataFrame = pygmt.blockmedian(
+            table=pd.concat(
+                objs=[table[[x_var, y_var]], z_values.mean(axis="columns")], axis=1
+            ),
+            region=grid_region,
+            spacing=f"{spacing}+e",
+        )
+        dataset["z_mean"]: xr.DataArray = pygmt.surface(
+            data=df_zmean.values, **surface_kwargs
+        )
+        dataset["z_mean"].attrs["long_name"] = f"Mean {z_var}"
+
+        # Rate of elevation change (dhdt)
+        df_dhdt: pd.DataFrame = pygmt.blockmedian(
+            table=table[[x_var, y_var, "dhdt_slope"]],
+            region=grid_region,
+            spacing=f"{spacing}+e",
+        )
+        dataset["dhdt"]: xr.DataArray = pygmt.surface(
+            data=df_dhdt.values, **surface_kwargs
+        )
+        dataset["dhdt"].attrs["long_name"] = f"Rate of elevation change over time"
+
     return dataset