🔀 Merge branch 'spatiotemporal_cube' (#180)

Closes #180 Spatiotemporal visualization of active subglacial lake surfaces over ICESat-2 cycles

atlxi_lake.py

@@ -44,13 +44,15 @@
 import dask
 import dask.array
 import geopandas as gpd
+import hvplot.xarray
 import numpy as np
 import pandas as pd
 import panel as pn
 import pygmt
 import scipy.spatial
 import shapely.geometry
 import tqdm
+import xarray as xr
 import zarr
 
 import deepicedrain
@@ -339,15 +341,18 @@ 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
-placename: str = "Recovery"  # "Whillans"
-drainage_basins: gpd.GeoDataFrame = drainage_basins.set_index(keys="NAME")
-region: deepicedrain.Region = deepicedrain.Region.from_gdf(
-    gdf=drainage_basins.loc[placename], name="Recovery Basin"
-)
+placename: str = "siple_coast"  # "Recovery"  # "Whillans"
+try:
+    drainage_basins: gpd.GeoDataFrame = drainage_basins.set_index(keys="NAME")
+    region: deepicedrain.Region = deepicedrain.Region.from_gdf(
+        gdf=drainage_basins.loc[placename], name="Recovery Basin"
+    )
+except KeyError:
+    pass
 df_dhdt: cudf.DataFrame = cudf.read_parquet(
     f"ATLXI/df_dhdt_{placename.lower()}.parquet"
 )
@@ -356,29 +361,158 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series:
 # %%
 # Antarctic subglacial lake polygons with EPSG:3031 coordinates
 antarctic_lakes: gpd.GeoDataFrame = gpd.read_file(
-    filename="antarctic_subglacial_lakes.geojson"
+    filename="antarctic_subglacial_lakes_3031.geojson"
 )
 
 # %%
 # Choose one draining/filling lake
-draining: bool = False  # False
-placename: str = "Slessor"  # "Whillans"
+draining: bool = False
+placename: str = "Whillans"  # "Slessor"  # "Kamb"  # "Mercer"  #
 lakes: gpd.GeoDataFrame = antarctic_lakes.query(expr="basin_name == @placename")
-lake = lakes.loc[lakes.maxabsdhdt.idxmin() if draining else lakes.maxabsdhdt.idxmax()]
+lake = lakes.loc[lakes.inner_dhdt.idxmin() if draining else lakes.inner_dhdt.idxmax()]
+# lake = lakes.query(expr="inner_dhdt < 0" if draining else "inner_dhdt > 0").loc[63]
 lakedict = {
-    76: "Subglacial Lake Conway",  # draining lake
-    78: "Whillans IX",  # filling lake
-    103: "Slessor 45",  # draining lake
-    108: "Slessor 23",  # filling lake
+    21: "Mercer 2b",  # filling lake
+    40: "Subglacial Lake Conway",  # draining lake
+    48: "Subglacial Lake Whillans",  # filling lake
+    50: "Whillans IX",  # filling lake
+    63: "Kamb 1",  # filling lake
+    65: "Kamb 12",  # filling lake
+    97: "MacAyeal 1",  # draining lake
+    109: "Slessor 45",  # draining lake
+    116: "Slessor 23",  # filling lake
+    153: "Recovery IX",  # draining lake
+    157: "Recovery 3",  # filling lake
 }
 region = deepicedrain.Region.from_gdf(gdf=lake, name=lakedict[lake.name])
 
+print(lake)
+lake.geometry
+
 # %%
 # Subset data to lake of interest
 placename: str = region.name.lower().replace(" ", "_")
 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("|")]
@@ -445,6 +579,10 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series:
 # Parallelized paired crossover analysis
 futures: list = []
 for rgt1, rgt2 in itertools.combinations(rgts, r=2):
+    # skip if same referencegroundtrack but different laser pair
+    # as they are parallel and won't cross
+    if rgt1[:4] == rgt2[:4]:
+        continue
     track1 = track_dict[rgt1][["x", "y", "h_corr", "utc_time"]]
     track2 = track_dict[rgt2][["x", "y", "h_corr", "utc_time"]]
     future = client.submit(
@@ -511,14 +649,20 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series:
 pygmt.makecpt(cmap="batlow", series=[sumstats[var]["25%"], sumstats[var]["75%"]])
 # Map frame in metre units
 fig.basemap(frame="f", region=plotregion, projection="X8c")
-# Plot actual track points
+# Plot actual track points in green
 for track in tracks:
-    fig.plot(x=track.x, y=track.y, color="green", style="c0.01c")
+    tracklabel = f"{track.iloc[0].referencegroundtrack} {track.iloc[0].pairtrack}"
+    fig.plot(
+        x=track.x,
+        y=track.y,
+        pen="thinnest,green,.",
+        style=f'qN+1:+l"{tracklabel}"+f3p,Helvetica,darkgreen',
+    )
 # Plot crossover point locations
 fig.plot(x=df.x, y=df.y, color=df.h_X, cmap=True, style="c0.1c", pen="thinnest")
-# PLot lake boundary
+# Plot lake boundary in blue
 lakex, lakey = lake.geometry.exterior.coords.xy
-fig.plot(x=lakex, y=lakey, pen="thin")
+fig.plot(x=lakex, y=lakey, pen="thin,blue,-.")
 # Map frame in kilometre units
 fig.basemap(
     frame=[
@@ -530,7 +674,7 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series:
     projection="X8c",
 )
 fig.colorbar(position="JMR", frame=['x+l"Crossover Error"', "y+lm"])
-fig.savefig(f"figures/crossover_area_{placename}.png")
+fig.savefig(f"figures/{placename}/crossover_area_{placename}_{min_date}_{max_date}.png")
 fig.show()
 
 
@@ -587,14 +731,14 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series:
 # Plot dashed line connecting points
 fig.plot(x=df_max.t, y=df_max.h, pen=f"faint,blue,-")
 fig.savefig(
-    f"figures/crossover_point_{placename}_{track1}_{track2}_{min_date}_{max_date}.png"
+    f"figures/{placename}/crossover_point_{placename}_{track1}_{track2}_{min_date}_{max_date}.png"
 )
 fig.show()
 
 # %%
 # Plot all crossover height points over time over the lake area
-fig = deepicedrain.vizplots.plot_crossovers(df=df_th, regionname=region.name)
-fig.savefig(f"figures/crossover_many_{placename}_{min_date}_{max_date}.png")
+fig = deepicedrain.plot_crossovers(df=df_th, regionname=region.name)
+fig.savefig(f"figures/{placename}/crossover_many_{placename}_{min_date}_{max_date}.png")
 fig.show()
 
 # %%
@@ -603,10 +747,15 @@ def find_clusters(X: cudf.core.dataframe.DataFrame) -> cudf.core.series.Series:
 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)
 
-fig = deepicedrain.vizplots.plot_crossovers(
-    df=df_th, regionname=region.name, elev_var="h_norm"
+fig = deepicedrain.plot_crossovers(
+    df=df_th,
+    regionname=region.name,
+    elev_var="h_norm",
+    elev_filter=3 * abs(df.h_X).median(),
+)
+fig.savefig(
+    f"figures/{placename}/crossover_many_normalized_{placename}_{min_date}_{max_date}.png"
 )
-fig.savefig(f"figures/crossover_many_normalized_{placename}_{min_date}_{max_date}.png")
 fig.show()
 
 # %%

deepicedrain/README.md

@@ -17,6 +17,7 @@ Contents:
   - Region - Bounding box data class structure that has xarray subsetting capabilities and more!
   - deltatime_to_utctime - Converts GPS time from an epoch (default is 2018 Jan 1st) to UTC time
   - lonlat_to_xy - Reprojects longitude/latitude EPSG:4326 coordinates to x/y EPSG:3031 coordinates
+  - spatiotemporal_cube - Interpolates a time-series point cloud into an xarray.Dataset data cube
 
 - :card_file_box: extraload.py - Convenience functions for extracting, transforming and loading data
   - array_to_dataframe - Turns a 1D/2D numpy/dask array into a tidy pandas/dask dataframe table

deepicedrain/__init__.py

@@ -11,6 +11,7 @@
     deltatime_to_utctime,
     lonlat_to_xy,
     point_in_polygon_gpu,
+    spatiotemporal_cube,
 )
 from deepicedrain.vizplots import IceSat2Explorer, plot_alongtrack, plot_crossovers