Fix various issues and tune some filtering parameters

ragmac_xdem/dem_postprocessing.py

@@ -766,9 +766,7 @@ def fname_func(pair_id):
             pair_id = pair_ids[count]
             print(f"\nProcessing pair {pair_id}")
 
-            dems_list = groups[count]
             dem_dates = utils.get_dems_date(dems_list)
-
             time_stamps = np.array(matplotlib.dates.date2num(dem_dates))
 #             time_stamps = np.array([utils.date_time_to_decyear(i) for i in dem_dates])
 
@@ -842,7 +840,7 @@ def fname_func(pair_id):
             min_date = np.percentile(ds.time, 2)
             max_date = np.percentile(ds.time, 98)
             time_delta_max = int((max_date - min_date).astype('timedelta64[D]') / np.timedelta64(1, 'D'))
-            time_delta_min = int(time_delta_max * 0.1)
+            time_delta_min = min(4*365, int(time_delta_max * 0.5))
             print("Min 2 percentile date:",np.datetime_as_string(min_date, unit='D'))
             print("Max 98 percentile date:",np.datetime_as_string(max_date, unit='D'))
             print("Time delta between dates:",time_delta_max, 'days')

ragmac_xdem/main.py

@@ -86,12 +86,12 @@ def main(case: dict, mode: str, run_name: str, sat_type: str = "ASTER", nproc: i
         selection_opts = {"mode": "subperiod", "dt": 365}
         downsampling = 1
         merge_opts = {"mode": "TimeSeries3"}
-    elif mode == "TimeSeries4":
+    elif mode == "TimeSeries_full":
         selection_opts = {"mode": "subperiod", "dt": 10958} #~30 years
         downsampling = 1
         merge_opts = {"mode": "TimeSeries3"}
     else:
-        raise ValueError("`mode` must be either of 'DEMdiff_autoselect', 'DEMdiff_median', 'TimeSeries', 'TimeSeries2' or 'TimeSeries3'")
+        raise ValueError("`mode` must be either of 'DEMdiff_autoselect', 'DEMdiff_median', 'TimeSeries', 'TimeSeries2', 'TimeSeries3' or TimeSeries_full")
 
     # Get run parameters
     run = runs[run_name]

ragmac_xdem/mass_balance.py

@@ -106,13 +106,21 @@ def fill_ddem_local_hypso(ddem, ref_dem, roi_mask, roi_outlines, filtering=True)
     `interp_residuals` the difference between input and interpolated ddem,
     `frac_obs` the fraction of observation for each feature in roi_outlines
     """
+    # Filter large outliers - as in Dussaillant et al. 2019 (DOI: 10.1038/s41561-019-0432-5)
+    if filtering:
+        slope = xdem.terrain.slope(ref_dem)
+        ddem_filt = ddem.copy()
+        ddem_filt.data.mask[slope.data > 45] = True
+    else:
+        ddem_filt = ddem
+
     # Calculate mean elevation change within elevation bins
     # TODO: filter pixels within each bins that are outliers
-    ddem_bins = xdem.volume.hypsometric_binning(ddem.data[roi_mask], ref_dem.data[roi_mask])
+    ddem_bins = xdem.volume.hypsometric_binning(ddem_filt.data[roi_mask], ref_dem.data[roi_mask], bins=100)
 
     # Filter outliers in bins
     if filtering:
-        ddem_bins_filtered = ddem_bins_filtering(ddem_bins, verbose=True)
+        ddem_bins_filtered = ddem_bins_filtering(ddem_bins, verbose=True, nmad_fact=-1)
     else:
         ddem_bins_filtered = ddem_bins.copy()
 

scripts/main_experiment1.py

@@ -70,7 +70,7 @@
     if args.mode is not None:
         all_modes = [args.mode]
 
-    if args.qc or args.mode == "TimeSeries3":
+    if args.qc or args.mode == "TimeSeries3" or args.mode == "TimeSeries_full":
         # Launch dask cluster for computation on larger than memory arrays.
         client = io.dask_start_cluster(args.nproc)
 

scripts/main_experiment2.py

@@ -73,7 +73,7 @@
     if args.run is not None:
         all_runs = [args.run]
 
-    if args.qc or args.mode == "TimeSeries3":
+    if args.qc or args.mode == "TimeSeries3" or args.mode == "TimeSeries_full":
         # Launch dask cluster for computation on larger than memory arrays.
         client = io.dask_start_cluster(args.nproc)