From 6c37eaf373329f87d9d9e62cce9e7f03e3892532 Mon Sep 17 00:00:00 2001 From: Wei Ji Date: Wed, 13 Jan 2021 09:40:19 +1300 Subject: [PATCH 1/2] :beers: Transect line over 3D ice surface elevation plots Spatial context matters, hence this series of hacks to overlay a transect line over our 3D ice surface elevation plots! The metadata of subglacial_lakes in atlas_catalog.yaml now includes a 'transect' field for the chosen transect to plot. The transect line is now plotted in bright yellow, and lake outline is now a cyan dash-dot line (same as in crossover plot), and the 3D ice_surface elevation plot uses colorbar shading (I) too! Moved some code out of conftest.py to test_subglacial_lake_crossover to ensure that the track points were not subset to within the lake only. Also atlas_catalog following some literature review. Whillans 12 was renamed to Lake 12, and Kamb 6 and MacAyeal 4 were added. Replaced Whillans 13 with *2, and labelled other potentially errorneous lakes from *1 to *5. --- atlxi_lake.ipynb | 26 +++++++++++++--- atlxi_lake.py | 26 +++++++++++++--- atlxi_xover.ipynb | 4 +-- atlxi_xover.py | 4 +-- deepicedrain/atlas_catalog.yaml | 30 +++++++++++++++---- .../features/subglacial_lakes.feature | 4 +-- deepicedrain/tests/conftest.py | 20 +++++++------ .../tests/test_subglacial_lake_animation.py | 14 ++++++++- .../tests/test_subglacial_lake_crossover.py | 9 ++++++ deepicedrain/vizplots.py | 16 ++++------ 10 files changed, 113 insertions(+), 40 deletions(-) diff --git a/atlxi_lake.ipynb b/atlxi_lake.ipynb index 4faa4d2..368b65a 100644 --- a/atlxi_lake.ipynb +++ b/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 = \"Lake 12\"\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,14 @@ "# 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 transect line xyz data\n", + "try:\n", + " _rgt, _pt = transect_id.split(\"_\")\n", + " df_transect: pd.DataFrame = deepicedrain.split_tracks(\n", + " df=df_lake.query(expr=f\"referencegroundtrack == {int(_rgt)}\").to_pandas()\n", + " )[transect_id][[\"x\", \"y\", \"h_corr\", \"cycle_number\"]]\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", @@ -1073,6 +1081,16 @@ " elevation=45, # 60\n", " title=f\"{region.name} at Cycle {cycle} ({time_sec})\",\n", " )\n", + " # Plot crossing transect line\n", + " _xyz = df_transect.query(expr=f\"cycle_number == {cycle}\")[[\"x\", \"y\", \"h_corr\"]]\n", + " if len(_xyz) > 0:\n", + " fig.plot3d(\n", + " data=_xyz.to_numpy(),\n", + " color=\"yellow2\",\n", + " style=\"c0.1c\",\n", + " zscale=True,\n", + " perspective=True,\n", + " )\n", " fig.savefig(f\"figures/{placename}/dsm_{placename}_cycle_{cycle}.png\")\n", "fig.show()" ] diff --git a/atlxi_lake.py b/atlxi_lake.py index 82564b6..fc0b9cf 100644 --- a/atlxi_lake.py +++ b/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,14 @@ # Subset data to lake of interest placename: str = region.name.lower().replace(" ", "_") df_lake: cudf.DataFrame = region.subset(data=df_dhdt) +# Get transect line xyz data +try: + _rgt, _pt = transect_id.split("_") + df_transect: pd.DataFrame = deepicedrain.split_tracks( + df=df_lake.query(expr=f"referencegroundtrack == {int(_rgt)}").to_pandas() + )[transect_id][["x", "y", "h_corr", "cycle_number"]] +except AttributeError: + pass # Save lake outline to OGR GMT file format outline_points: str = f"figures/{placename}/{placename}.gmt" @@ -407,6 +415,16 @@ elevation=45, # 60 title=f"{region.name} at Cycle {cycle} ({time_sec})", ) + # Plot crossing transect line + _xyz = df_transect.query(expr=f"cycle_number == {cycle}")[["x", "y", "h_corr"]] + if len(_xyz) > 0: + fig.plot3d( + data=_xyz.to_numpy(), + color="yellow2", + style="c0.1c", + zscale=True, + perspective=True, + ) fig.savefig(f"figures/{placename}/dsm_{placename}_cycle_{cycle}.png") fig.show() diff --git a/atlxi_xover.ipynb b/atlxi_xover.ipynb index 6703f84..e1b31e1 100644 --- a/atlxi_xover.ipynb +++ b/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", diff --git a/atlxi_xover.py b/atlxi_xover.py index 116ad87..4448902 100644 --- a/atlxi_xover.py +++ b/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"], ["", ".-", "-", "..-"], ) ): diff --git a/deepicedrain/atlas_catalog.yaml b/deepicedrain/atlas_catalog.yaml index 840e315..1cda34b 100644 --- a/deepicedrain/atlas_catalog.yaml +++ b/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: diff --git a/deepicedrain/features/subglacial_lakes.feature b/deepicedrain/features/subglacial_lakes.feature index ed69d26..6cc110b 100644 --- a/deepicedrain/features/subglacial_lakes.feature +++ b/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 | diff --git a/deepicedrain/tests/conftest.py b/deepicedrain/tests/conftest.py index e2b2041..fdbef4c 100644 --- a/deepicedrain/tests/conftest.py +++ b/deepicedrain/tests/conftest.py @@ -58,10 +58,10 @@ def lake_altimetry_data(lake_name: str, location: str, context) -> pd.DataFrame: # Get lake outline from intake catalog 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] ) context.lake: pd.Series = ( lake_catalog.read() @@ -79,12 +79,14 @@ def lake_altimetry_data(lake_name: str, location: str, context) -> pd.DataFrame: # Subset data to lake of interest context.placename: str = context.lake_name.lower().replace(" ", "_") df_lake: cudf.DataFrame = context.region.subset(data=dataframe) # 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(context.lake.geometry) - ] # polygon subset + # Get transect line xyz data + try: + _rgt, _pt = transect_id.split("_") + context.df_transect: pd.DataFrame = deepicedrain.split_tracks( + df=df_lake.query(expr=f"referencegroundtrack == {int(_rgt)}") + )[transect_id][["x", "y", "h_corr", "cycle_number"]] + except AttributeError: + pass # Save lake outline to OGR GMT file format os.makedirs(name=f"figures/{context.placename}", exist_ok=True) diff --git a/deepicedrain/tests/test_subglacial_lake_animation.py b/deepicedrain/tests/test_subglacial_lake_animation.py index 42d8fae..ee9bca8 100644 --- a/deepicedrain/tests/test_subglacial_lake_animation.py +++ b/deepicedrain/tests/test_subglacial_lake_animation.py @@ -5,7 +5,6 @@ import os import subprocess -import geopandas as gpd import numpy as np import pandas as pd import pygmt @@ -108,6 +107,19 @@ def visualize_grid_in_3D( elevation=elevation, title=f"{context.region.name} at Cycle {cycle} ({time_sec})", ) + # Plot crossing transect line + if hasattr(context, "df_transect"): + _xyz = context.df_transect.query(expr=f"cycle_number == {cycle}")[ + ["x", "y", "h_corr"] + ] + if len(_xyz) > 0: + fig.plot3d( + data=_xyz.to_numpy(), + color="yellow2", + style="c0.1c", + zscale=True, + perspective=True, + ) fig.savefig( f"figures/{context.placename}/dsm_{context.placename}_cycle_{cycle}.png" ) diff --git a/deepicedrain/tests/test_subglacial_lake_crossover.py b/deepicedrain/tests/test_subglacial_lake_crossover.py index 04ba488..0eb0039 100644 --- a/deepicedrain/tests/test_subglacial_lake_crossover.py +++ b/deepicedrain/tests/test_subglacial_lake_crossover.py @@ -5,6 +5,7 @@ import itertools import dask +import geopandas as gpd import numpy as np import pandas as pd import pint @@ -39,6 +40,14 @@ def crossover_height_anomaly( df_lake: pd.DataFrame, client: dask.distributed.Client, context ) -> pd.DataFrame: + # Subset data to lake of interest + 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(context.lake.geometry) + ] # polygon subset + # Run crossover analysis on all tracks track_dict: dict = deepicedrain.split_tracks(df=df_lake) rgts, tracks = track_dict.keys(), track_dict.values() diff --git a/deepicedrain/vizplots.py b/deepicedrain/vizplots.py index 7af357d..d352cf1 100644 --- a/deepicedrain/vizplots.py +++ b/deepicedrain/vizplots.py @@ -280,16 +280,11 @@ def _plot_crossover_area( # Map frame in metre units fig.basemap(frame="+n", region=plotregion, projection=f"X{plotsize}c") - # Plot lake boundary in blue - fig.plot(data=outline_points, region=plotregion, pen="thin,blue,-.") + # Plot lake boundary in cyan + fig.plot(data=outline_points, region=plotregion, pen="thin,cyan2,-.") # Plot crossover point locations pygmt.makecpt(cmap=cmap, series=elev_range, reverse=True) - fig.plot( - data=df, - style="d0.1c", - cmap=True, - # pen="thinnest", - ) + fig.plot(data=df, style="d0.1c", cmap=True) # Map frame in kilometre units with pygmt.config(FONT_ANNOT_PRIMARY=f"{plotsize+2}p", FONT_LABEL=f"{plotsize+2}p"): fig.basemap( @@ -603,6 +598,7 @@ def plot_icesurface( ) fig.colorbar( cmap=True, + I=True, # shading position="JMR+o1c/0c+w7c/0.5c+n+mc", frame=['x+l"Elevation (m)"', "y+lm"], perspective=True, @@ -617,7 +613,7 @@ def plot_icesurface( zscale=True, perspective=True, ) - # Plot lake boundary outline as yellow dashed line + # Plot lake boundary outline as cyan dashed line if outline_points is not None: with pygmt.helpers.GMTTempFile() as tmpfile: pygmt.grdtrack( @@ -639,7 +635,7 @@ def plot_icesurface( fig.plot3d( data=tmpfile.name, region=grid_region, - pen="1.5p,yellow2,-", + pen="thicker,cyan2,-.", zscale=True, perspective=True, ) From faad68860ed8c4e6dcfa3d0e9d0f044e9ce9e4d7 Mon Sep 17 00:00:00 2001 From: Wei Ji Date: Wed, 13 Jan 2021 16:29:06 +1300 Subject: [PATCH 2/2] :clown_face: Single 3D plot of mean ice surface elevation and dhdt trend It's tough to animate on paper, so let's take the average and trend instead! The spatiotemporal_cube function will now calculate mean elevation () and rate of elevation change (dhdt) and include it in the output xarray.Dataset. The plot_icesurface code has been modified slightly to handle dhdt as well as dh, and we ensure the roma colormap is centered on zero. Getting all the ICESat-2 raw xyz points and one transect line was complicated (read: needs refactoring). Yellow transect is now plotted only on the new static 3D plot with dhdt trend, leaving the older 3D plots (which becomes a GIF) pretty much intact as before. --- atlxi_lake.ipynb | 63 ++++++++++++++----- atlxi_lake.py | 55 +++++++++++----- deepicedrain/spatiotemporal.py | 50 +++++++++++---- deepicedrain/tests/conftest.py | 18 ++++-- .../tests/test_subglacial_lake_animation.py | 40 +++++++----- deepicedrain/vizplots.py | 35 ++++++++--- 6 files changed, 189 insertions(+), 72 deletions(-) diff --git a/atlxi_lake.ipynb b/atlxi_lake.ipynb index 368b65a..4be6b12 100644 --- a/atlxi_lake.ipynb +++ b/atlxi_lake.ipynb @@ -951,7 +951,7 @@ ], "source": [ "# Choose one Antarctic active subglacial lake polygon with EPSG:3031 coordinates\n", - "lake_name: str = \"Lake 12\"\n", + "lake_name: str = \"Whillans 7\"\n", "lake_catalog = deepicedrain.catalog.subglacial_lakes()\n", "lake_ids, transect_id = (\n", " pd.json_normalize(lake_catalog.metadata[\"lakedict\"])\n", @@ -983,12 +983,22 @@ "# 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 transect line xyz data\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 = deepicedrain.split_tracks(\n", - " df=df_lake.query(expr=f\"referencegroundtrack == {int(_rgt)}\").to_pandas()\n", - " )[transect_id][[\"x\", \"y\", \"h_corr\", \"cycle_number\"]]\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", @@ -1064,7 +1074,38 @@ "metadata": {}, "outputs": [], "source": [ - "# 3D plots of gridded ice surface elevation over time\n", + "# 3D plot of mean ice surface elevation () 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", @@ -1081,16 +1122,6 @@ " elevation=45, # 60\n", " title=f\"{region.name} at Cycle {cycle} ({time_sec})\",\n", " )\n", - " # Plot crossing transect line\n", - " _xyz = df_transect.query(expr=f\"cycle_number == {cycle}\")[[\"x\", \"y\", \"h_corr\"]]\n", - " if len(_xyz) > 0:\n", - " fig.plot3d(\n", - " data=_xyz.to_numpy(),\n", - " color=\"yellow2\",\n", - " style=\"c0.1c\",\n", - " zscale=True,\n", - " perspective=True,\n", - " )\n", " fig.savefig(f\"figures/{placename}/dsm_{placename}_cycle_{cycle}.png\")\n", "fig.show()" ] diff --git a/atlxi_lake.py b/atlxi_lake.py index fc0b9cf..b977dca 100644 --- a/atlxi_lake.py +++ b/atlxi_lake.py @@ -354,12 +354,22 @@ # Subset data to lake of interest placename: str = region.name.lower().replace(" ", "_") df_lake: cudf.DataFrame = region.subset(data=df_dhdt) -# Get transect line xyz data +# 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 = deepicedrain.split_tracks( - df=df_lake.query(expr=f"referencegroundtrack == {int(_rgt)}").to_pandas() - )[transect_id][["x", "y", "h_corr", "cycle_number"]] + 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 @@ -398,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 () 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") @@ -415,16 +448,6 @@ elevation=45, # 60 title=f"{region.name} at Cycle {cycle} ({time_sec})", ) - # Plot crossing transect line - _xyz = df_transect.query(expr=f"cycle_number == {cycle}")[["x", "y", "h_corr"]] - if len(_xyz) > 0: - fig.plot3d( - data=_xyz.to_numpy(), - color="yellow2", - style="c0.1c", - zscale=True, - perspective=True, - ) fig.savefig(f"figures/{placename}/dsm_{placename}_cycle_{cycle}.png") fig.show() @@ -441,7 +464,7 @@ "120", "-loop", "0", - f"figures/{placename}/dsm_*.png", + f"figures/{placename}/dsm_*cycle_*.png", gif_fname, ] ) diff --git a/deepicedrain/spatiotemporal.py b/deepicedrain/spatiotemporal.py index ddf45bb..c34de09 100644 --- a/deepicedrain/spatiotemporal.py +++ b/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 () + 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 diff --git a/deepicedrain/tests/conftest.py b/deepicedrain/tests/conftest.py index fdbef4c..a7bb935 100644 --- a/deepicedrain/tests/conftest.py +++ b/deepicedrain/tests/conftest.py @@ -79,12 +79,22 @@ def lake_altimetry_data(lake_name: str, location: str, context) -> pd.DataFrame: # Subset data to lake of interest context.placename: str = context.lake_name.lower().replace(" ", "_") df_lake: cudf.DataFrame = context.region.subset(data=dataframe) # bbox subset - # Get transect line xyz data + # Get all raw xyz points and one transect line dataframe + track_dict: dict = deepicedrain.split_tracks(df=df_lake) + context.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("_") - context.df_transect: pd.DataFrame = deepicedrain.split_tracks( - df=df_lake.query(expr=f"referencegroundtrack == {int(_rgt)}") - )[transect_id][["x", "y", "h_corr", "cycle_number"]] + context.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 diff --git a/deepicedrain/tests/test_subglacial_lake_animation.py b/deepicedrain/tests/test_subglacial_lake_animation.py index ee9bca8..c77d299 100644 --- a/deepicedrain/tests/test_subglacial_lake_animation.py +++ b/deepicedrain/tests/test_subglacial_lake_animation.py @@ -90,6 +90,31 @@ def visualize_grid_in_3D( print(f"Elevation limits are: {z_limits}") + # %% + # 3D plot of mean ice surface elevation () and rate of ice elevation change (dhdt) + fig = deepicedrain.plot_icesurface( + grid=f"figures/{context.placename}/xyht_{context.placename}.nc?z_mean", + grid_region=grid_region, + diff_grid=f"figures/{context.placename}/xyht_{context.placename}.nc?dhdt", + track_points=context.track_points.to_numpy(), + outline_points=context.outline_points, + azimuth=azimuth, + elevation=elevation, + title=f"{context.region.name} Ice Surface", + ) + # Plot crossing transect line + if hasattr(context, "df_transect"): + fig.plot3d( + data=context.df_transect[["x", "y", "h_corr"]].to_numpy(), + color="yellow2", + style="c0.1c", + zscale=True, + perspective=True, + ) + fig.savefig( + f"figures/{context.placename}/dsm_{context.placename}_cycles_{context.cycles[0]}-{context.cycles[-1]}.png" + ) + # 3D plots of gridded ice surface elevation over time for cycle in tqdm.tqdm(iterable=context.cycles): time_nsec: pd.Timestamp = df_lake[f"utc_time_{cycle}"].mean() @@ -107,19 +132,6 @@ def visualize_grid_in_3D( elevation=elevation, title=f"{context.region.name} at Cycle {cycle} ({time_sec})", ) - # Plot crossing transect line - if hasattr(context, "df_transect"): - _xyz = context.df_transect.query(expr=f"cycle_number == {cycle}")[ - ["x", "y", "h_corr"] - ] - if len(_xyz) > 0: - fig.plot3d( - data=_xyz.to_numpy(), - color="yellow2", - style="c0.1c", - zscale=True, - perspective=True, - ) fig.savefig( f"figures/{context.placename}/dsm_{context.placename}_cycle_{cycle}.png" ) @@ -144,7 +156,7 @@ def animated_dsm_gif(context) -> str: "120", "-loop", "0", - f"figures/{context.placename}/dsm_*.png", + f"figures/{context.placename}/dsm_*cycle_*.png", gif_fname, ] ) diff --git a/deepicedrain/vizplots.py b/deepicedrain/vizplots.py index d352cf1..16ea291 100644 --- a/deepicedrain/vizplots.py +++ b/deepicedrain/vizplots.py @@ -470,7 +470,7 @@ def plot_crossovers( def plot_icesurface( grid: str or xr.DataArray = None, grid_region: tuple or np.ndarray = None, - diff_grid: xr.DataArray = None, + diff_grid: str or xr.DataArray = None, diff_grid_region: tuple or np.ndarray = None, track_points: pd.DataFrame = None, outline_points: str or pd.DataFrame = None, @@ -544,11 +544,21 @@ def plot_icesurface( ## Bottom plot # Normalized ice surface elevation change grid - if diff_grid.min() == diff_grid.max(): - # add some tiny random noise to make plot work - np.random.seed(seed=int(elevation)) - diff_grid = diff_grid + abs(np.random.normal(scale=1e-32, size=diff_grid.shape)) - pygmt.makecpt(cmap="roma", series=diff_grid_region[-2:]) + try: + if diff_grid.min() == diff_grid.max(): + # add some tiny random noise to make plot work + np.random.seed(seed=int(elevation)) + diff_grid = diff_grid + abs( + np.random.normal(scale=1e-32, size=diff_grid.shape) + ) + except AttributeError: + pass + try: + series = diff_grid_region[-2:] + except TypeError: + series = pygmt.grdinfo(grid=diff_grid, T="1+s")[2:-3] + finally: + pygmt.makecpt(cmap="roma", series=series) fig.grdview( grid=diff_grid, projection="X10c", @@ -569,8 +579,13 @@ def plot_icesurface( ) fig.colorbar( cmap=True, - position="JMR+o1c/0c+w7c/0.5c+n+mc", - frame=['x+l"Elevation Change (m)"', "y+lm"], + position="JMR+o1c/0c+w7c/0.5c+n", + frame=[ + 'x+l"Elevation Trend"', + "y+lm/yr", + ] + if "?dhdt" in diff_grid + else ['x+l"Elevation Change"', "y+lm"], perspective=True, ) @@ -599,8 +614,8 @@ def plot_icesurface( fig.colorbar( cmap=True, I=True, # shading - position="JMR+o1c/0c+w7c/0.5c+n+mc", - frame=['x+l"Elevation (m)"', "y+lm"], + position="JMR+o1c/0c+w7c/0.5c+n", + frame=['x+l"Elevation"', "y+lm"], perspective=True, )