diff --git a/atlxi_lake.ipynb b/atlxi_lake.ipynb
index 2e8bdb6..f4aa4d7 100644
--- a/atlxi_lake.ipynb
+++ b/atlxi_lake.ipynb
@@ -42,6 +42,7 @@
"import dask\n",
"import dask.array\n",
"import geopandas as gpd\n",
+ "import hvplot.xarray\n",
"import numpy as np\n",
"import pandas as pd\n",
"import panel as pn\n",
@@ -49,6 +50,7 @@
"import scipy.spatial\n",
"import shapely.geometry\n",
"import tqdm\n",
+ "import xarray as xr\n",
"import zarr\n",
"\n",
"import deepicedrain"
@@ -905,12 +907,12 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Select a lake to examine"
+ "# Select a subglacial lake to examine"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 13,
"metadata": {
"lines_to_next_cell": 2
},
@@ -932,7 +934,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
@@ -944,11 +946,41 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 15,
"metadata": {
"lines_to_next_cell": 1
},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "basin_name Whillans\n",
+ "num_points 3276\n",
+ "outer_dhdt 0.377611\n",
+ "outer_std 1.42712\n",
+ "outer_mad 0.083278\n",
+ "inner_dhdt 2.12607\n",
+ "maxabsdhdt 8.65662\n",
+ "refgtracks 74|135|196|266|327|388|516|577|638|769|830|101...\n",
+ "geometry POLYGON ((-444681.0351994165 -545210.454471163...\n",
+ "Name: 50, dtype: object\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 15,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"# Choose one draining/filling lake\n",
"draining: bool = False\n",
@@ -977,7 +1009,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 16,
"metadata": {
"lines_to_next_cell": 2
},
@@ -988,6 +1020,178 @@
"df_lake: cudf.DataFrame = region.subset(data=df_dhdt)"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Create an interpolated ice surface elevation grid for each ICESat-2 cycle"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "100%|██████████| 6/6 [00:03<00:00, 1.71it/s]\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Generate gridded time-series of ice elevation over lake\n",
+ "cycles: tuple = (3, 4, 5, 6, 7, 8)\n",
+ "os.makedirs(name=f\"figures/{placename}\", exist_ok=True)\n",
+ "ds_lake: xr.Dataset = deepicedrain.spatiotemporal_cube(\n",
+ " table=df_lake.to_pandas(),\n",
+ " placename=placename,\n",
+ " cycles=cycles,\n",
+ " folder=f\"figures/{placename}\",\n",
+ ")\n",
+ "ds_lake.to_netcdf(path=f\"figures/{placename}/xyht_{placename}.nc\", mode=\"w\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Elevation limits are: (273.2418518066406, 310.9927062988281)\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Get 3D grid_region (xmin/xmax/ymin/ymax/zmin/zmax),\n",
+ "# and calculate normalized z-values as Elevation delta relative to Cycle 3\n",
+ "grid_region = pygmt.info(table=df_lake[[\"x\", \"y\"]].to_pandas(), spacing=\"s250\")\n",
+ "z_limits: tuple = (float(ds_lake.z.min()), float(ds_lake.z.max())) # original z limits\n",
+ "grid_region: np.ndarray = np.append(arr=grid_region, values=z_limits)\n",
+ "\n",
+ "ds_lake_norm: xr.Dataset = ds_lake - ds_lake.sel(cycle_number=3).z\n",
+ "z_norm_limits: tuple = (float(ds_lake_norm.z.min()), float(ds_lake_norm.z.max()))\n",
+ "grid_region_norm: np.ndarray = np.append(arr=grid_region[:4], values=z_norm_limits)\n",
+ "\n",
+ "print(f\"Elevation limits are: {z_limits}\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# 3D plots of gridded ice surface elevation over time\n",
+ "azimuth: float = 157.5 # 202.5 # 270\n",
+ "elevation: float = 45 # 60\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",
+ "\n",
+ " grdview_kwargs = dict(\n",
+ " cmap=True,\n",
+ " zscale=0.1, # zscaling factor, default to 10x vertical exaggeration\n",
+ " surftype=\"sim\", # surface, image and mesh plot\n",
+ " perspective=[azimuth, elevation], # perspective using azimuth/elevation\n",
+ " # W=\"c0.05p,black,solid\", # draw contours\n",
+ " )\n",
+ "\n",
+ " fig = pygmt.Figure()\n",
+ " # grid = ds_lake.sel(cycle_number=cycle).z\n",
+ " grid = f\"figures/{placename}/h_corr_{placename}_cycle_{cycle}.nc\"\n",
+ " pygmt.makecpt(cmap=\"lapaz\", series=z_limits)\n",
+ " fig.grdview(\n",
+ " grid=grid,\n",
+ " projection=\"X10c\",\n",
+ " region=grid_region,\n",
+ " shading=True,\n",
+ " frame=[\n",
+ " f'SWZ+t\"{region.name}\"', # plot South, West axes, and Z-axis\n",
+ " 'xaf+l\"Polar Stereographic X (m)\"', # add x-axis annotations and minor ticks\n",
+ " 'yaf+l\"Polar Stereographic Y (m)\"', # add y-axis annotations and minor ticks\n",
+ " f'zaf+l\"Elevation (m)\"', # add z-axis annotations, minor ticks and axis label\n",
+ " ],\n",
+ " **grdview_kwargs,\n",
+ " )\n",
+ "\n",
+ " # Plot lake boundary outline\n",
+ " # TODO wait for plot3d to plot lake boundary points at correct height\n",
+ " df = pd.DataFrame([region.bounds()]).values\n",
+ " points = pd.DataFrame(\n",
+ " data=[point for point in lake.geometry.exterior.coords], columns=(\"x\", \"y\")\n",
+ " )\n",
+ " df_xyz = pygmt.grdtrack(points=points, grid=grid, newcolname=\"z\")\n",
+ " fig.plot(\n",
+ " data=df_xyz.values,\n",
+ " region=grid_region,\n",
+ " pen=\"1.5p,yellow2\",\n",
+ " Jz=True, # zscale\n",
+ " p=f\"{azimuth}/{elevation}/{df_xyz.z.median()}\", # perspective\n",
+ " # label='\"Subglacial Lake X\"'\n",
+ " )\n",
+ "\n",
+ " # Plot normalized elevation change\n",
+ " grid = ds_lake_norm.sel(cycle_number=cycle).z\n",
+ " if cycle == 3:\n",
+ " # add some tiny random noise to make plot work\n",
+ " grid = grid + np.random.normal(scale=1e-8, size=grid.shape)\n",
+ " pygmt.makecpt(cmap=\"vik\", series=z_norm_limits)\n",
+ " fig.grdview(\n",
+ " grid=grid,\n",
+ " region=grid_region_norm,\n",
+ " frame=[\n",
+ " f'SEZ2+t\"Cycle {cycle} at {time_sec}\"', # plot South, East axes, and Z-axis\n",
+ " 'xaf+l\"Polar Stereographic X (m)\"', # add x-axis annotations and minor ticks\n",
+ " 'yaf+l\"Polar Stereographic Y (m)\"', # add y-axis annotations and minor ticks\n",
+ " f'zaf+l\"Elev Change (m)\"', # add z-axis annotations, minor ticks and axis label\n",
+ " ],\n",
+ " X=\"10c\", # xshift\n",
+ " **grdview_kwargs,\n",
+ " )\n",
+ "\n",
+ " fig.savefig(f\"figures/{placename}/dsm_{placename}_cycle_{cycle}.png\")\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Make a animated GIF of changing ice surface from the PNG files\n",
+ "gif_fname: str = (\n",
+ " f\"figures/{placename}/dsm_{placename}_cycles_{cycles[0]}-{cycles[-1]}.gif\"\n",
+ ")\n",
+ "# !convert -delay 120 -loop 0 figures/{placename}/dsm_*.png {gif_fname}"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# HvPlot 2D interactive view of ice surface elevation grids over each ICESat-2 cycle\n",
+ "dashboard: pn.layout.Column = pn.Column(\n",
+ " ds_lake.hvplot.image(x=\"x\", y=\"y\", clim=z_limits, cmap=\"gist_earth\", data_aspect=1)\n",
+ " # * ds_lake.hvplot.contour(x=\"x\", y=\"y\", clim=z_limits, data_aspect=1)\n",
+ ")\n",
+ "dashboard.show(port=30227)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Along track plots of ice surface elevation change over time"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
diff --git a/atlxi_lake.py b/atlxi_lake.py
index d008d44..edebad3 100644
--- a/atlxi_lake.py
+++ b/atlxi_lake.py
@@ -44,6 +44,7 @@
import dask
import dask.array
import geopandas as gpd
+import hvplot.xarray
import numpy as np
import pandas as pd
import panel as pn
@@ -51,6 +52,7 @@
import scipy.spatial
import shapely.geometry
import tqdm
+import xarray as xr
import zarr
import deepicedrain
@@ -339,7 +341,7 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series:
# %% [markdown]
-# # Select a lake to examine
+# # Select a subglacial lake to examine
# %%
# Save or load dhdt data from Parquet file
@@ -393,6 +395,124 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series:
df_lake: cudf.DataFrame = region.subset(data=df_dhdt)
+# %% [markdown]
+# ## Create an interpolated ice surface elevation grid for each ICESat-2 cycle
+
+# %%
+# Generate gridded time-series of ice elevation over lake
+cycles: tuple = (3, 4, 5, 6, 7, 8)
+os.makedirs(name=f"figures/{placename}", exist_ok=True)
+ds_lake: xr.Dataset = deepicedrain.spatiotemporal_cube(
+ table=df_lake.to_pandas(),
+ placename=placename,
+ cycles=cycles,
+ folder=f"figures/{placename}",
+)
+ds_lake.to_netcdf(path=f"figures/{placename}/xyht_{placename}.nc", mode="w")
+
+# %%
+# Get 3D grid_region (xmin/xmax/ymin/ymax/zmin/zmax),
+# and calculate normalized z-values as Elevation delta relative to Cycle 3
+grid_region = pygmt.info(table=df_lake[["x", "y"]].to_pandas(), spacing="s250")
+z_limits: tuple = (float(ds_lake.z.min()), float(ds_lake.z.max())) # original z limits
+grid_region: np.ndarray = np.append(arr=grid_region, values=z_limits)
+
+ds_lake_norm: xr.Dataset = ds_lake - ds_lake.sel(cycle_number=3).z
+z_norm_limits: tuple = (float(ds_lake_norm.z.min()), float(ds_lake_norm.z.max()))
+grid_region_norm: np.ndarray = np.append(arr=grid_region[:4], values=z_norm_limits)
+
+print(f"Elevation limits are: {z_limits}")
+
+# %%
+# 3D plots of gridded ice surface elevation over time
+azimuth: float = 157.5 # 202.5 # 270
+elevation: float = 45 # 60
+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")
+
+ grdview_kwargs = dict(
+ cmap=True,
+ zscale=0.1, # zscaling factor, default to 10x vertical exaggeration
+ surftype="sim", # surface, image and mesh plot
+ perspective=[azimuth, elevation], # perspective using azimuth/elevation
+ # W="c0.05p,black,solid", # draw contours
+ )
+
+ fig = pygmt.Figure()
+ # grid = ds_lake.sel(cycle_number=cycle).z
+ grid = f"figures/{placename}/h_corr_{placename}_cycle_{cycle}.nc"
+ pygmt.makecpt(cmap="lapaz", series=z_limits)
+ fig.grdview(
+ grid=grid,
+ projection="X10c",
+ region=grid_region,
+ shading=True,
+ frame=[
+ f'SWZ+t"{region.name}"', # plot South, West axes, and Z-axis
+ 'xaf+l"Polar Stereographic X (m)"', # add x-axis annotations and minor ticks
+ 'yaf+l"Polar Stereographic Y (m)"', # add y-axis annotations and minor ticks
+ f'zaf+l"Elevation (m)"', # add z-axis annotations, minor ticks and axis label
+ ],
+ **grdview_kwargs,
+ )
+
+ # Plot lake boundary outline
+ # TODO wait for plot3d to plot lake boundary points at correct height
+ df = pd.DataFrame([region.bounds()]).values
+ points = pd.DataFrame(
+ data=[point for point in lake.geometry.exterior.coords], columns=("x", "y")
+ )
+ df_xyz = pygmt.grdtrack(points=points, grid=grid, newcolname="z")
+ fig.plot(
+ data=df_xyz.values,
+ region=grid_region,
+ pen="1.5p,yellow2",
+ Jz=True, # zscale
+ p=f"{azimuth}/{elevation}/{df_xyz.z.median()}", # perspective
+ # label='"Subglacial Lake X"'
+ )
+
+ # Plot normalized elevation change
+ grid = ds_lake_norm.sel(cycle_number=cycle).z
+ if cycle == 3:
+ # add some tiny random noise to make plot work
+ grid = grid + np.random.normal(scale=1e-8, size=grid.shape)
+ pygmt.makecpt(cmap="vik", series=z_norm_limits)
+ fig.grdview(
+ grid=grid,
+ region=grid_region_norm,
+ frame=[
+ f'SEZ2+t"Cycle {cycle} at {time_sec}"', # plot South, East axes, and Z-axis
+ 'xaf+l"Polar Stereographic X (m)"', # add x-axis annotations and minor ticks
+ 'yaf+l"Polar Stereographic Y (m)"', # add y-axis annotations and minor ticks
+ f'zaf+l"Elev Change (m)"', # add z-axis annotations, minor ticks and axis label
+ ],
+ X="10c", # xshift
+ **grdview_kwargs,
+ )
+
+ fig.savefig(f"figures/{placename}/dsm_{placename}_cycle_{cycle}.png")
+fig.show()
+
+# %%
+# Make a animated GIF of changing ice surface from the PNG files
+gif_fname: str = (
+ f"figures/{placename}/dsm_{placename}_cycles_{cycles[0]}-{cycles[-1]}.gif"
+)
+# !convert -delay 120 -loop 0 figures/{placename}/dsm_*.png {gif_fname}
+
+# %%
+# HvPlot 2D interactive view of ice surface elevation grids over each ICESat-2 cycle
+dashboard: pn.layout.Column = pn.Column(
+ ds_lake.hvplot.image(x="x", y="y", clim=z_limits, cmap="gist_earth", data_aspect=1)
+ # * ds_lake.hvplot.contour(x="x", y="y", clim=z_limits, data_aspect=1)
+)
+dashboard.show(port=30227)
+
+# %% [markdown]
+# ## Along track plots of ice surface elevation change over time
+
# %%
# Select a few Reference Ground tracks to look at
rgts: list = [int(rgt) for rgt in lake.refgtracks.split("|")]