add optional custom grid for merging gridded data (OGGM#1779)

* add optional custom grid for merging gridded data

* fix bug when reprojecting relatively small grid to a larger one

* fix bug when merging 3D data

* added whats-new

Author: pat-schmitt

docs/whats-new.rst

@@ -47,6 +47,10 @@ Enhancements
 - New standard glacier directories now use a new reference lookup
   table to decide on the topo data to use (:pull:`1781`).
   By `Fabien Maussion <https://github.com/fmaussion>`_
+- Added the possibility to provide a custom grid  in 
+  ``workflow.merge_gridded_data``. If no grid is provided, the default is to
+  merge all grids of the provided gdirs (:pull:`1779`).
+  By `Patrick Schmitt <https://github.com/pat-schmitt>`_
 
 Bug fixes
 ~~~~~~~~~

oggm/core/gis.py

@@ -1992,6 +1992,8 @@ def reproject_gridded_data_variable_to_grid(gdir,
             total_before = (np.nansum(data.values, axis=sum_axis) *
                             ds.salem.grid.dx ** 2)
 
+            total_before_is_zero = np.isclose(total_before, 0, atol=1e-6)
+
             if smooth_radius != 0:
                 if smooth_radius is None:
                     smooth_radius = np.rint(cfg.PARAMS['smooth_window'] /
@@ -2008,12 +2010,40 @@ def reproject_gridded_data_variable_to_grid(gdir,
 
             total_after = (np.nansum(r_data, axis=sum_axis) *
                            target_grid.dx ** 2)
+            total_after_is_zero = np.isclose(total_after, 0, atol=1e-6)
+
+            # if a relatively small grid is reprojected into a larger grid, it
+            # could happen that no data is assigned at all
+            no_data_after_but_before = np.logical_and(total_after_is_zero,
+                                                      ~total_before_is_zero)
+            if np.any(no_data_after_but_before):
+                # we just assign the maximum value to one grid point and use the
+                # factor for conserving the total value
+                def _assign_max_value(data_provided, data_target):
+                    j_max, i_max = np.unravel_index(
+                        np.nanargmax(data_provided.values), data_provided.shape)
+                    oi_max, oj_max = target_grid.center_grid.transform(
+                        i_max, j_max, crs=gdir.grid.center_grid, nearest=True)
+                    data_target[oj_max, oi_max] = data_provided[j_max, i_max]
+                    return data_target
+
+                if r_data.ndim == 3:
+                    for i in range(r_data.shape[0]):
+                        if no_data_after_but_before[i]:
+                            r_data[i, :, :] = _assign_max_value(data[i, :, :],
+                                                                r_data[i, :, :])
+                else:
+                    r_data = _assign_max_value(data, r_data)
+
+                # and recalculate the total after again
+                total_after = (np.nansum(r_data, axis=sum_axis) *
+                               target_grid.dx ** 2)
 
             # only preserve total if there is some data before
             with warnings.catch_warnings():
                 # Divide by zero is fine
                 warnings.filterwarnings("ignore", category=RuntimeWarning)
-                factor = np.where(np.isclose(total_before, 0, atol=1e-6),
+                factor = np.where(total_before_is_zero,
                                   0., total_before / total_after)
 
             if len(data.dims) == 3:

oggm/sandbox/distribute_2d.py

@@ -490,6 +490,7 @@ def distribute_thickness_from_simulation(gdir,
 def merge_simulated_thickness(gdirs,
                               output_folder=None,
                               output_filename=None,
+                              output_grid=None,
                               simulation_filesuffix='',
                               years_to_merge=None,
                               keep_dem_file=False,
@@ -514,6 +515,9 @@ def merge_simulated_thickness(gdirs,
     output_folder
     output_filename : str
         Default is gridded_simulation_merged{simulation_filesuffix}.
+    output_grid : salem.gis.Grid
+        You can provide a custom grid on which the distributed data should be
+        merged on.
     simulation_filesuffix : str
         the filesuffix of the gridded_simulation file
     years_to_merge : list | xr.DataArray | None
@@ -559,6 +563,7 @@ def _calc_bedrock_topo(fp):
             gdirs,
             output_folder=output_folder,
             output_filename=f"{output_filename}_topo_data",
+            output_grid=output_grid,
             input_file='gridded_data',
             input_filesuffix='',
             included_variables=['glacier_ext',
@@ -607,6 +612,7 @@ def _calc_bedrock_topo(fp):
                 gdirs,
                 output_folder=output_folder,
                 output_filename=f"{output_filename}_{year}_{month:02d}",
+                output_grid=output_grid,
                 input_file='gridded_simulation',
                 input_filesuffix=simulation_filesuffix,
                 included_variables=[('simulated_thickness',
@@ -632,6 +638,7 @@ def _calc_bedrock_topo(fp):
             gdirs,
             output_folder=output_folder,
             output_filename=output_filename,
+            output_grid=output_grid,
             input_file=['gridded_data', 'gridded_simulation'],
             input_filesuffix=['', simulation_filesuffix],
             included_variables=[['glacier_ext',

oggm/tests/test_workflow.py

@@ -414,6 +414,28 @@ def test_merge_gridded_data():
     assert_allclose(df['inv_volume_km3'].sum(), inv_volume_gridded_merged,
                     rtol=1e-6)
 
+    # test providing a custom grid
+    default_grid = utils.combine_grids(gdirs)
+    custom_grid_dict = default_grid.to_dict()
+    custom_grid_dict['nxny'] = (custom_grid_dict['nxny'][0] - 2,
+                                custom_grid_dict['nxny'][1] - 2)
+    custom_grid = salem.gis.Grid.from_dict(custom_grid_dict)
+
+    ds_merged_custom = workflow.merge_gridded_data(
+        gdirs,
+        output_grid=custom_grid,
+        included_variables='distributed_thickness',
+        reset=True)
+
+    assert ds_merged_custom.salem.grid.nx < ds_merged.salem.grid.nx
+    assert ds_merged_custom.salem.grid.ny < ds_merged.salem.grid.ny
+
+    # total volume should be the same with a custom grid
+    inv_volume_gridded_merged_custom = (ds_merged_custom.distributed_thickness.sum() *
+                                        ds_merged_custom.salem.grid.dx ** 2) * 1e-9
+    assert_allclose(inv_volume_gridded_merged_custom, inv_volume_gridded_merged,
+                    rtol=1e-6)
+
 
 @pytest.mark.slow
 @pytest.mark.graphic

oggm/workflow.py

@@ -848,6 +848,7 @@ def _recursive_merging(gdirs, gdir_main, glcdf=None, dem_source=None,
 @global_task(log)
 def merge_gridded_data(gdirs, output_folder=None,
                        output_filename='gridded_data_merged',
+                       output_grid=None,
                        input_file='gridded_data',
                        input_filesuffix='',
                        included_variables='all',
@@ -884,6 +885,10 @@ def merge_gridded_data(gdirs, output_folder=None,
         should be stored. Default is cfg.PATHS['working_dir']
     output_filename : str
         The name for the resulting file. Default is 'gridded_data_merged'.
+    output_grid : salem.gis.Grid
+        You can provide a custom grid on which the gridded data should be
+        merged on. If None, a combined grid of all gdirs will be constructed.
+        Default is None.
     input_file : str or list
         The file(s) which should be merged. If a list is provided the data of
         all files is merged into the same dataset. Default is 'gridded_data'.
@@ -970,9 +975,10 @@ def merge_gridded_data(gdirs, output_folder=None,
         # into a list of lists
         included_variables = [included_variables]
 
-    # create a combined salem.Grid object, which serves as canvas/boundaries of
-    # the combined glacier region
-    combined_grid = utils.combine_grids(gdirs)
+    if output_grid is None:
+        # create a combined salem.Grid object, which serves as canvas/boundaries of
+        # the combined glacier region
+        output_grid = utils.combine_grids(gdirs)
 
     if add_topography:
         # ok, lets get a DEM and add it to the final file
@@ -982,11 +988,11 @@ def merge_gridded_data(gdirs, output_folder=None,
         else:
             dem_source = None
             dem_gdir = gdirs[0]
-        gis.get_dem_for_grid(combined_grid, output_folder,
+        gis.get_dem_for_grid(output_grid, output_folder,
                              source=dem_source, gdir=dem_gdir)
         # unwrapped is needed to execute process_dem without the entity_task
         # overhead (this would need a valid gdir)
-        gis.process_dem.unwrapped(gdir=None, grid=combined_grid,
+        gis.process_dem.unwrapped(gdir=None, grid=output_grid,
                                   fpath=output_folder,
                                   output_filename=output_filename)
         if not keep_dem_file:
@@ -999,7 +1005,7 @@ def merge_gridded_data(gdirs, output_folder=None,
             if os.path.exists(fpath):
                 os.remove(fpath)
 
-    with gis.GriddedNcdfFile(grid=combined_grid, fpath=output_folder,
+    with gis.GriddedNcdfFile(grid=output_grid, fpath=output_folder,
                              basename=output_filename) as nc:
 
         # adding the data of one file after another to the merged dataset
@@ -1036,9 +1042,9 @@ def merge_gridded_data(gdirs, output_folder=None,
                 dim_lengths = []
                 for dim in dims:
                     if dim == 'y':
-                        dim_lengths.append(combined_grid.ny)
+                        dim_lengths.append(output_grid.ny)
                     elif dim == 'x':
-                        dim_lengths.append(combined_grid.nx)
+                        dim_lengths.append(output_grid.nx)
                     else:
                         if slice_of_var is not None:
                             # only keep selected part of the variable
@@ -1088,7 +1094,7 @@ def merge_gridded_data(gdirs, output_folder=None,
 
                 kwargs_reproject = dict(
                     variable=var,
-                    target_grid=combined_grid,
+                    target_grid=output_grid,
                     filename=in_file,
                     filesuffix=in_filesuffix,
                     use_glacier_mask=use_glacier_mask,