Recalculate ICESat-2 ATL11 height changes up to 20210715

.gitignore

@@ -36,7 +36,12 @@ ATLXI/df_*.parquet
 ATLXI/ds_*.zarr
 Quantarctica3
 
+# Dask
+**/dask-worker-space/
+
 # Subglacial Lake grid files and figures
 figures/**/*.gif
+figures/**/*.gmt
 figures/**/*.nc
 figures/**/*.png
+figures/**/*.txt

README.md

@@ -12,7 +12,7 @@ in Antarctica using remote sensing and machine learning.
 
 | Ice Surface Elevation trends over Antactica | Active Subglacial Lake fill-drain event |
 |---|---|
-| ![ICESat-2 ATL11 rate of height change over time in Antarctica 2019-03-29 to 2020-12-24](https://user-images.githubusercontent.com/23487320/123902132-65cfd680-d9c0-11eb-88d6-4e0e8c5abc47.png) | ![dsm_whillans_ix_cycles_3-9.gif](https://user-images.githubusercontent.com/23487320/124219379-5ed7ce00-db50-11eb-95d0-f1f660d4d688.gif) |
+| ![ICESat-2 ATL11 rate of height change over time in Antarctica 2019-03-29 to 2021-07-15](https://user-images.githubusercontent.com/23487320/138962441-1288f99c-1bfd-4405-a516-5deee5bb4492.png) | ![dsm_whillans_ix_cycles_3-9.gif](https://user-images.githubusercontent.com/23487320/124219379-5ed7ce00-db50-11eb-95d0-f1f660d4d688.gif) |
 
 ![DeepIceDrain Pipeline Part 1 Exploratory Data Analysis](https://yuml.me/diagram/scruffy;dir:LR/class/[Land-Ice-Elevation|atl06_play.ipynb]->[Convert|atl06_to_atl11.ipynb],[Convert]->[Land-Ice-Height-time-series|atl11_play.ipynb])
 ![DeepIceDrain Pipeline Part 2 Subglacial Lake Analysis](https://yuml.me/diagram/scruffy;dir:LR/class/[Height-Change-over-Time-(dhdt)|atlxi_dhdt.ipynb],[Height-Change-over-Time-(dhdt)]->[Subglacial-Lake-Finder|atlxi_lake.ipynb],[Subglacial-Lake-Finder]->[Crossover-Analysis|atlxi_xover.ipynb])

atlxi_dhdt.py

@@ -7,7 +7,7 @@
 #       extension: .py
 #       format_name: hydrogen
 #       format_version: '1.3'
-#       jupytext_version: 1.11.3
+#       jupytext_version: 1.11.4
 #   kernelspec:
 #     display_name: deepicedrain
 #     language: python
@@ -54,8 +54,8 @@
 import deepicedrain
 
 # %%
-client = dask.distributed.Client(n_workers=16, threads_per_worker=1)
-client
+client = dask.distributed.Client(n_workers=8, threads_per_worker=1)
+print(client)
 
 # %% [markdown]
 # # Select essential points
@@ -138,7 +138,7 @@
 
 # %%
 # Get first and last dates to put into our plots
-min_date, max_date = ("2019-03-29", "2020-12-24")
+min_date, max_date = ("2019-03-29", "2021-07-15")
 if min_date is None:
     min_delta_time = np.nanmin(ds.delta_time.isel(cycle_number=0).data).compute()
     min_utc_time = deepicedrain.deltatime_to_utctime(min_delta_time)
@@ -185,10 +185,10 @@
 
 # %%
 # Save or Load height range data
-# ds_ht.to_zarr(store=f"ATLXI/ds_hrange_time_{placename}.zarr", mode="w", consolidated=True)
+# ds_ht.to_zarr(store=f"ATLXI/ds_hrange_time_{placename}.zarr", mode="w")
 ds_ht: xr.Dataset = xr.open_dataset(
     filename_or_obj=f"ATLXI/ds_hrange_time_{placename}.zarr",
-    chunks={"cycle_number": 7},
+    chunks={"cycle_number": 10},
     engine="zarr",
     backend_kwargs={"consolidated": True},
 )
@@ -241,7 +241,7 @@
 #
 # Performing linear regression in parallel.
 # Uses the [`scipy.stats.linregress`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html) function,
-# parallelized with xarray's [`apply_ufunc`](http://xarray.pydata.org/en/v0.15.1/examples/apply_ufunc_vectorize_1d.html) method
+# parallelized with xarray's [`apply_ufunc`](https://xarray.pydata.org/en/v0.19.0/examples/apply_ufunc_vectorize_1d.html) method
 # on a Dask cluster.
 
 # %%
@@ -262,8 +262,7 @@
     dask="parallelized",
     vectorize=True,
     output_dtypes=[np.float32],
-    output_sizes={"dhdt_parameters": 5},
-    # output_sizes={"slope_ns":1, "intercept":1, "r_value":1, "p_value":1, "std_err":1}
+    dask_gufunc_kwargs={"output_sizes": dict(dhdt_parameters=5)},
 )
 
 # %%
@@ -286,20 +285,21 @@
 
 # %%
 # Do linear regression on single datapoint
-# slope_ns, intercept, r_value, p_value, std_err = nan_linregress(
-#     x=ds.delta_time[:1].data.astype(np.uint64), y=ds.h_corr[:1].data
+# slope_ns, intercept, r_value, p_value, std_err = deepicedrain.nan_linregress(
+#     x=ds.delta_time[:1].astype(np.uint64).data, y=ds.h_corr[:1].data
 # )
 # print(slope_ns, intercept, r_value, p_value, std_err)
 
 # %%
 # Load or Save rate of height change over time (dhdt) data
-# ds_dhdt.to_zarr(store=f"ATLXI/ds_dhdt_{placename}.zarr", mode="w", consolidated=True)
+# ds_dhdt.to_zarr(store=f"ATLXI/ds_dhdt_{placename}.zarr", mode="w")
 ds_dhdt: xr.Dataset = xr.open_dataset(
     filename_or_obj=f"ATLXI/ds_dhdt_{placename}.zarr",
     chunks="auto",
     engine="zarr",
     backend_kwargs={"consolidated": True},
 )
+# ds: xr.Dataset = ds_dhdt  # shortcut for dhdt_maxslp calculation later
 
 # %%
 df_slope: pd.DataFrame = ds_dhdt.dhdt_slope.to_dataframe()
@@ -343,6 +343,69 @@
 
 # %%
 
+# %% [markdown]
+# # Calculate rate of height change over time max slope (dhdt_maxslp)
+#
+# For each ATL11 data point, find the maximum slope (i.e. steepest gradient)
+# for any consecutive paired value within the elevation time-series. Uses a
+# custom `dhdt_maxslp` function, parallelized with xarray's
+# [`apply_ufunc`](https://xarray.pydata.org/en/v0.19.0/examples/apply_ufunc_vectorize_1d.html) method
+# on a Dask cluster.
+#
+# For example, in the plot below, the rate of elevation change over time is
+# greatest from point B to C, so the algorithm will return the dhdt_maxslp
+# value as (elev_C - elev_B) / (time_C - time_B).
+#
+#              ^
+#              |        C
+#              |           D
+#     elev (m) |
+#              |     B        E
+#              |  A              F
+#              -------------------->
+#                      time
+#
+# Note that NaN values are ignored in the calculation. So if point E had a
+# NaN value, the algorithm will calculate dhdt between point F and D.
+
+# %%
+print(f"Running on {len(ds.ref_pt)} ATL11 points")
+
+# %%
+# Do dhdt_maxslp calculation on many datapoints, parallelized using dask
+ds["dhdt_maxslp"]: xr.DataArray = xr.apply_ufunc(
+    deepicedrain.dhdt_maxslp,
+    ds.delta_time.astype(np.uint64),  # x is time in nanoseconds
+    ds.h_corr,  # y is height in metres
+    input_core_dims=[["cycle_number"], ["cycle_number"]],
+    dask="parallelized",
+    vectorize=True,
+    output_dtypes=[np.float32],
+)
+
+# %%
+# Convert dhdt_maxslp units from metres per nanosecond to metres per year
+# 1 year = 365.25 days x 24 hours x 60 min x 60 seconds x 1_000_000_000 nanoseconds
+ds["dhdt_maxslp"] = ds["dhdt_maxslp"] * (365.25 * 24 * 60 * 60 * 1_000_000_000)
+
+# %% time
+# %%time
+# Compute rate of height change over time max slope (dhdt_maxslp).
+# Also include all height and time info
+ds_dhdt: xr.Dataset = ds[["dhdt_slope", "h_corr", "dhdt_maxslp"]].compute()
+
+# %%
+# Load or Save rate of height change over time max slope (dhdt_maxslp) data
+ds_dhdt.to_zarr(store=f"ATLXI/ds_dhdt_maxslp_{placename}.zarr", mode="w")
+
+
+# %%
+# Do dhdt_maxslp calculation on single datapoint
+# dhdt_maxslp = deepicedrain.dhdt_maxslp(
+#     x=ds.delta_time[:1].astype(np.uint64).to_numpy(), y=ds.h_corr[:1].to_numpy()
+# )
+# print(dhdt_maxslp)
+
 # %% [markdown]
 # # Along track plots of subglacial lake drainage/filling events
 #