Merge pull request #33 from adehecq/amaury

Few updates for first RAGMAC submission

ragmac_xdem/mass_balance.py

@@ -97,7 +97,7 @@ def fill_ddem(
     return filled_ddem
 
 
-def fill_ddem_local_hypso(ddem, ref_dem, roi_mask, plot=True, outfig=None):
+def fill_ddem_local_hypso(ddem, ref_dem, roi_mask, roi_outlines, plot=True, outfig=None):
     """
     Function to fill gaps in ddems using a local hypsometric approach.
     """
@@ -159,13 +159,19 @@ def fill_ddem_local_hypso(ddem, ref_dem, roi_mask, plot=True, outfig=None):
         ax3.tick_params(axis="x", colors=p3.patches[0].get_facecolor(), **tkw)
 
         # ddem before interpolation
+        bounds = roi_outlines.bounds
+        pad = 2e3
         ax2 = plt.subplot(132)
-        ddem.show(ax=ax2, cmap="coolwarm_r", vmin=-50, vmax=50)
+        roi_outlines.ds.plot(ax=ax2, facecolor="none", edgecolor="k", zorder=2)
+        ddem.show(ax=ax2, cmap="coolwarm_r", vmin=-50, vmax=50, zorder=1)
+        plt.xlim(bounds.left - pad, bounds.right + pad)
+        plt.ylim(bounds.bottom - pad, bounds.top + pad)
         plt.title("dDEM before interpolation")
 
         # ddem before interpolation
         ax3 = plt.subplot(133, sharex=ax2, sharey=ax2)
-        ddem_filled.show(ax=ax3, cmap="coolwarm_r", vmin=-50, vmax=50)
+        roi_outlines.ds.plot(ax=ax3, facecolor="none", edgecolor="k", zorder=2)
+        ddem_filled.show(ax=ax3, cmap="coolwarm_r", vmin=-50, vmax=50, zorder=1)
         plt.title("dDEM after interpolation")
 
         if outfig is None:

ragmac_xdem/utils.py

@@ -172,7 +172,7 @@ def best_dem_cover(dem_path_list: list, init_stats: pd.Series) -> list[str, floa
     # Select highest ROI coverage
     best = stats_subset.sort_values(by="roi_cover_orig").iloc[-1]
 
-    return best.dem_path, best.roi_cover_orig
+    return best.dem_path, best.roi_cover_orig, best.dem_date
 
 
 def list_pairs(validation_dates):
@@ -263,8 +263,9 @@ def dems_selection(
         else:
             assert init_stats is not None, "`init_stats` must be provided for mode 'best'"
             final_dem_list = []
-            for group in output_list:
-                selected_dem, _ = best_dem_cover(group, init_stats)
+            for date, group in zip(validation_dates, output_list):
+                selected_dem, _, dem_date = best_dem_cover(group, init_stats)
+                print(f"For date {date}, selected DEM from date {dem_date}")
                 final_dem_list.append(
                     [
                         selected_dem,

scripts/main_experiment1.py

@@ -73,11 +73,17 @@
 
     # -- Select different processing modes -- #
     if args.mode == "median":
-        selection_opts = {"mode": "close", "dt": 365, "months": [8, 9, 10]}
+        selection_opts = {"mode": "close", "dt": 400, "months": [8, 9, 10]}
         merge_opts = {"mode": "median"}
         outdir = os.path.join(exp["processed_data"]["directory"], "results_median")
         downsampling = 1
-        method = "DEMdiff"
+        method = "DEMdiff_median"
+    elif args.mode == "best":
+        selection_opts = {"mode": "best", "dt": 400, "months": [8, 9, 10]}
+        merge_opts = {"mode": "median"}
+        outdir = os.path.join(exp["processed_data"]["directory"], "results_best")
+        downsampling = 1
+        method = "DEMdiff_best"
     elif args.mode == "shean":
         selection_opts = {"mode": "subperiod", "dt": 365}
         downsampling = 1
@@ -91,7 +97,7 @@
         outdir = os.path.join(exp["processed_data"]["directory"], "results_knuth")
         method = "TimeSeries2"
     else:
-        raise ValueError("`mode` must be either of 'median','shean' or 'knuth'")
+        raise ValueError("`mode` must be either of 'median', 'best', 'shean' or 'knuth'")
 
     # create output directories
     if not os.path.exists(coreg_dir):
@@ -111,6 +117,8 @@
     # -- Select DEMs to be processed -- #
     print("\n### DEMs selection ###")
     validation_dates = exp["validation_dates"]
+    if selection_opts["mode"] == 'best':
+        selection_opts["init_stats"] = init_stats
     groups = utils.dems_selection(dems_files, validation_dates=validation_dates, **selection_opts)
 
     # -- Postprocess DEMs i.e. coregister, filter etc -- #
@@ -126,7 +134,7 @@
         plot=True,
         method="mp",
     )
-    print(f"--> Coregistered DEMs saved in {outdir}")
+    print(f"--> Coregistered DEMs saved in {coreg_dir}")
 
     # Temporarily downsample DEM for speeding-up testing
     if downsampling > 1:
@@ -170,7 +178,7 @@
         print(pair_id)
         fig_fn = os.path.join(outdir, f"{pair_id}_mb_fig.png")
         ddem_filled, ddem_bins = mb.fill_ddem_local_hypso(
-            ddems[pair_id], ref_dem, roi_mask, plot=True, outfig=fig_fn
+            ddems[pair_id], ref_dem, roi_mask, roi_outlines, plot=True, outfig=fig_fn
         )
         ddems_filled[pair_id] = ddem_filled
 

scripts/main_experiment2.py

@@ -68,11 +68,17 @@
 
     # -- Select different processing modes -- #
     if args.mode == "median":
-        selection_opts = {"mode": "close", "dt": 365, "months": [8, 9, 10]}
+        selection_opts = {"mode": "close", "dt": 400, "months": [8, 9, 10]}
         merge_opts = {"mode": "median"}
         outdir = os.path.join(baltoro_paths["processed_data"]["directory"], "results_median")
         downsampling = 1
-        method = "DEMdiff"
+        method = "DEMdiff_median"
+    elif args.mode == "best":
+        selection_opts = {"mode": "best", "dt": 400, "months": [8, 9, 10]}
+        merge_opts = {"mode": "median"}
+        outdir = os.path.join(baltoro_paths["processed_data"]["directory"], "results_best")
+        downsampling = 1
+        method = "DEMdiff_best"
     elif args.mode == "shean":
         selection_opts = {"mode": "subperiod", "dt": 365}
         downsampling = 1
@@ -106,6 +112,8 @@
     # -- Select DEMs to be processed -- #
     print("\n### DEMs selection ###")
     validation_dates = baltoro_paths["validation_dates"]
+    if selection_opts["mode"] == 'best':
+        selection_opts["init_stats"] = init_stats
     groups = utils.dems_selection(dems_files, validation_dates=validation_dates, **selection_opts)
 
     # -- Postprocess DEMs i.e. coregister, filter etc -- #
@@ -167,7 +175,7 @@
         print(pair_id)
         fig_fn = os.path.join(outdir, f"{pair_id}_mb_fig.png")
         ddem_filled, ddem_bins = mb.fill_ddem_local_hypso(
-            ddems[pair_id], ref_dem, roi_mask, plot=True, outfig=fig_fn
+            ddems[pair_id], ref_dem, roi_mask, roi_outlines, plot=True, outfig=fig_fn
         )
         ddems_filled[pair_id] = ddem_filled