From a4a1e01b7f5573d7756f92eda2156ffae18c6327 Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Mon, 27 Nov 2023 15:06:31 +0100
Subject: [PATCH 01/16] moved combine_grids function from graphics to utils

---
 oggm/graphics.py     | 75 +-------------------------------------------
 oggm/utils/_funcs.py | 72 ++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 73 insertions(+), 74 deletions(-)

diff --git a/oggm/graphics.py b/oggm/graphics.py
index 1241aa3f0..57a4f023a 100644
--- a/oggm/graphics.py
+++ b/oggm/graphics.py
@@ -76,79 +76,6 @@ def surf_to_nan(surf_h, thick):
     return surf_h
 
 
-def combine_grids(gdirs):
-    """ Combines individual grids of different glacier directories to show
-        multiple glaciers in the same plot. The resulting grid extent includes
-        all individual grids completely.
-
-    Parameters
-    ----------
-    gdirs : [], required
-        A list of GlacierDirectories.
-
-    Returns
-    -------
-    salem.gis.Grid
-    """
-
-    new_grid = {
-        'proj': None,
-        'nxny': None,
-        'dxdy': None,
-        'x0y0': None,
-        'pixel_ref': None
-    }
-
-    left_use = None
-    right_use = None
-    bottom_use = None
-    top_use = None
-    dx_use = None
-    dy_use = None
-
-    for gdir in gdirs:
-        # use the first gdir to define some values
-        if new_grid['proj'] is None:
-            new_grid['proj'] = gdir.grid.proj
-        if new_grid['pixel_ref'] is None:
-            new_grid['pixel_ref'] = gdir.grid.pixel_ref
-
-        # find largest extend including all grids completely
-        (left, right, bottom, top) = gdir.grid.extent_in_crs(new_grid['proj'])
-        if (left_use is None) or (left_use > left):
-            left_use = left
-        if right_use is None or right_use < right:
-            right_use = right
-        if bottom_use is None or bottom_use > bottom:
-            bottom_use = bottom
-        if top_use is None or top_use < top:
-            top_use = top
-
-        # find smallest dx and dy for the estimation of nx and ny
-        dx = gdir.grid.dx
-        dy = gdir.grid.dy
-        if dx_use is None or dx_use > dx:
-            dx_use = dx
-        # dy could be negative
-        if dy_use is None or abs(dy_use) > abs(dy):
-            dy_use = dy
-
-    # calculate nx and ny, the final extend could be one grid point larger or
-    # smaller due to round()
-    nx_use = round((right_use - left_use) / dx_use)
-    ny_use = round((top_use - bottom_use) / abs(dy_use))
-
-    # finally define the last values of the new grid
-    if np.sign(dy_use) < 0:
-        new_grid['x0y0'] = (left_use, top_use)
-    else:
-        new_grid['x0y0'] = (left_use, bottom_use)
-    new_grid['nxny'] = (nx_use, ny_use)
-    new_grid['dxdy'] = (dx_use, dy_use)
-
-    return salem.gis.Grid.from_dict(new_grid)
-
-
 def _plot_map(plotfunc):
     """
     Decorator for common salem.Map plotting logic
@@ -213,7 +140,7 @@ def newplotfunc(gdirs, ax=None, smap=None, add_colorbar=True, title=None,
 
         if smap is None:
             if extend_plot_limit:
-                grid_combined = combine_grids(gdirs)
+                grid_combined = utils.combine_grids(gdirs)
                 mp = salem.Map(grid_combined, countries=False,
                                nx=grid_combined.nx)
             else:
diff --git a/oggm/utils/_funcs.py b/oggm/utils/_funcs.py
index 0470a91b2..beecd63d1 100644
--- a/oggm/utils/_funcs.py
+++ b/oggm/utils/_funcs.py
@@ -24,6 +24,7 @@
 import shapely.geometry as shpg
 from shapely.ops import linemerge
 from shapely.validation import make_valid
+from salem.gis import Grid
 
 # Optional libs
 try:
@@ -538,6 +539,77 @@ def recursive_valid_polygons(geoms, crs):
     assert np.all([type(geom) == shpg.Polygon for geom in new_geoms])
     return new_geoms
 
+def combine_grids(gdirs):
+    """ Combines individual grids of different glacier directories to show
+        multiple glaciers in the same plot. The resulting grid extent includes
+        all individual grids completely.
+
+    Parameters
+    ----------
+    gdirs : [], required
+        A list of GlacierDirectories.
+
+    Returns
+    -------
+    salem.gis.Grid
+    """
+
+    new_grid = {
+        'proj': None,
+        'nxny': None,
+        'dxdy': None,
+        'x0y0': None,
+        'pixel_ref': None
+    }
+
+    left_use = None
+    right_use = None
+    bottom_use = None
+    top_use = None
+    dx_use = None
+    dy_use = None
+
+    for gdir in gdirs:
+        # use the first gdir to define some values
+        if new_grid['proj'] is None:
+            new_grid['proj'] = gdir.grid.proj
+        if new_grid['pixel_ref'] is None:
+            new_grid['pixel_ref'] = gdir.grid.pixel_ref
+
+        # find largest extend including all grids completely
+        (left, right, bottom, top) = gdir.grid.extent_in_crs(new_grid['proj'])
+        if (left_use is None) or (left_use > left):
+            left_use = left
+        if right_use is None or right_use < right:
+            right_use = right
+        if bottom_use is None or bottom_use > bottom:
+            bottom_use = bottom
+        if top_use is None or top_use < top:
+            top_use = top
+
+        # find smallest dx and dy for the estimation of nx and ny
+        dx = gdir.grid.dx
+        dy = gdir.grid.dy
+        if dx_use is None or dx_use > dx:
+            dx_use = dx
+        # dy could be negative
+        if dy_use is None or abs(dy_use) > abs(dy):
+            dy_use = dy
+
+    # calculate nx and ny, the final extend could be one grid point larger or
+    # smaller due to round()
+    nx_use = round((right_use - left_use) / dx_use)
+    ny_use = round((top_use - bottom_use) / abs(dy_use))
+
+    # finally define the last values of the new grid
+    if np.sign(dy_use) < 0:
+        new_grid['x0y0'] = (left_use, top_use)
+    else:
+        new_grid['x0y0'] = (left_use, bottom_use)
+    new_grid['nxny'] = (nx_use, ny_use)
+    new_grid['dxdy'] = (dx_use, dy_use)
+
+    return Grid.from_dict(new_grid)
 
 def multipolygon_to_polygon(geometry, gdir=None):
     """Sometimes an RGI geometry is a multipolygon: this should not happen.

From 816cc25d19ca2d7900f4a2dac5cb69706220a6f7 Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Mon, 27 Nov 2023 15:33:55 +0100
Subject: [PATCH 02/16] introduced 'combine_distributed_thickness'
 functionality - not finished yet

---
 oggm/sandbox/distribute_2d.py | 42 ++++++++++++++++++++++++++++++++++-
 1 file changed, 41 insertions(+), 1 deletion(-)

diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py
index d8fdcdb76..44cdfe473 100644
--- a/oggm/sandbox/distribute_2d.py
+++ b/oggm/sandbox/distribute_2d.py
@@ -9,7 +9,7 @@
 import xarray as xr
 from scipy import ndimage
 from scipy.stats import mstats
-from oggm.core.gis import gaussian_blur
+from oggm.core.gis import gaussian_blur, dem_for_combined_grid, CombinedNcdfFile
 from oggm.utils import ncDataset, entity_task
 import matplotlib.pyplot as plt
 
@@ -416,3 +416,43 @@ def distribute_thickness_from_simulation(gdir,
         return ds, out_df
 
     return ds
+
+def combine_distributed_thickness(gdirs, output_folder, source=None):
+    """
+    This function takes a list of glacier directories that have a distributed_thickness
+    Parameters
+    ----------
+    gdirs - list of GlacierDirectories
+    output_folder - folder where the intermediate files and the final combined distributed thickness
+                    will be stored
+    source - DEM source
+
+    Returns
+    -------
+
+    """
+    # create a combined salem.Grid object, which serves as canvas/boundaries of the combined glacier region
+    combined_grid = utils.combine_grids(gdirs)
+
+    #create gridded_data.nc
+    with CombinedNcdfFile(combined_grid, output_folder) as nc:
+        # is there still something else needed in the gridded_data.nc? add it here or use function
+        pass
+
+    # retrieve and project a topography/DEM for the defined region
+    dem_for_combined_grid(combined_grid, output_folder, source=source)
+
+    # add 'write_dem_to_gridded_data' function and use it here
+
+
+    # add individual add individual distributed thicknesses to gridded data
+    # for gdir in gdirs:
+    #     with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
+    #         ds = ds.load()
+    #
+    #     r_data = combined_grid.map_gridded_data(ds[f'simulation_distributed_thickness_{ssp}_projection'],
+    #                                             grid=gdir.grid)
+    #     total_data += r_data.filled(0)
+    # total_data = np.where(total_data == 0, np.NaN, total_data)
+
+

From 2296a14abc9b816986d902627074b1373406d737 Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Mon, 27 Nov 2023 16:00:00 +0100
Subject: [PATCH 03/16] split up and redesigned the 'define_glacier_region' to
 also use parts of it for the new function 'dem_for_combined_grid'

---
 oggm/core/gis.py | 327 +++++++++++++++++++++++++++++++++--------------
 1 file changed, 230 insertions(+), 97 deletions(-)

diff --git a/oggm/core/gis.py b/oggm/core/gis.py
index ea9ad9c35..8b39d7e0e 100644
--- a/oggm/core/gis.py
+++ b/oggm/core/gis.py
@@ -297,76 +297,52 @@ def glacier_grid_params(gdir):
 
     return utm_proj, nx, ny, ulx, uly, dx
 
-
-@entity_task(log, writes=['glacier_grid', 'dem', 'outlines'])
-def define_glacier_region(gdir, entity=None, source=None):
-    """Very first task after initialization: define the glacier's local grid.
-
-    Defines the local projection (Transverse Mercator or UTM depending on
-    user choice), centered on the glacier.
-    There is some options to set the resolution of the local grid.
-    It can be adapted depending on the size of the glacier.
-    See ``params.cfg`` for setting these options.
-
-    Default values of the adapted mode lead to a resolution of 50 m for
-    Hintereisferner, which is approx. 8 km2 large.
-
-    After defining the grid, the topography and the outlines of the glacier
-    are transformed into the local projection.
-
+def check_dem_source(source, extent_ll, rgi_id=None):
+    """
+    This function can check for multiple DEM sources and is in charge of the error handling if a requested source is
+    not available for the given glacier/extent
     Parameters
     ----------
-    gdir : :py:class:`oggm.GlacierDirectory`
-        where to write the data
-    entity : geopandas.GeoSeries
-        the glacier geometry to process - DEPRECATED. It is now ignored
-    source : str or list of str, optional
-        If you want to force the use of a certain DEM source. Available are:
-          - 'USER' : file set in cfg.PATHS['dem_file']
-          - 'SRTM' : http://srtm.csi.cgiar.org/
-          - 'GIMP' : https://bpcrc.osu.edu/gdg/data/gimpdem
-          - 'RAMP' : https://nsidc.org/data/nsidc-0082/versions/2/documentation
-          - 'REMA' : https://www.pgc.umn.edu/data/rema/
-          - 'DEM3' : http://viewfinderpanoramas.org/
-          - 'ASTER' : https://asterweb.jpl.nasa.gov/gdem.asp
-          - 'TANDEM' : https://geoservice.dlr.de/web/dataguide/tdm90/
-          - 'ARCTICDEM' : https://www.pgc.umn.edu/data/arcticdem/
-          - 'AW3D30' : https://www.eorc.jaxa.jp
-          - 'MAPZEN' : https://registry.opendata.aws/terrain-tiles/
-          - 'ALASKA' : https://www.the-cryosphere.net/8/503/2014/
-          - 'COPDEM30' : Copernicus DEM GLO30 https://spacedata.copernicus.eu/web/cscda/cop-dem-faq
-          - 'COPDEM90' : Copernicus DEM GLO90 https://spacedata.copernicus.eu/web/cscda/cop-dem-faq
-          - 'NASADEM':  https://doi.org/10.5069/G93T9FD9
-    """
-
-    utm_proj, nx, ny, ulx, uly, dx = glacier_grid_params(gdir)
+    source
+    extent_ll
+    gdir
 
-    # Back to lon, lat for DEM download/preparation
-    tmp_grid = salem.Grid(proj=utm_proj, nxny=(nx, ny), x0y0=(ulx, uly),
-                          dxdy=(dx, -dx), pixel_ref='corner')
-    minlon, maxlon, minlat, maxlat = tmp_grid.extent_in_crs(crs=salem.wgs84)
+    Returns
+    -------
 
-    # Open DEM
-    # We test DEM availability for glacier only (maps can grow big)
+    """
+# We test DEM availability for glacier only (maps can grow big)
     if isinstance(source, list):  # when multiple sources are provided, try them sequentially
         for src in source:
-            source_exists = is_dem_source_available(src, *gdir.extent_ll)
+            source_exists = is_dem_source_available(src, *extent_ll)
             if source_exists:
                 source = src  # pick the first source which exists
                 break
     else:
-        source_exists = is_dem_source_available(source, *gdir.extent_ll)
-        
+        source_exists = is_dem_source_available(source, *extent_ll)
     if not source_exists:
-        raise InvalidWorkflowError(f'Source: {source} is not available for '
-                                   f'glacier {gdir.rgi_id} with border '
-                                   f"{cfg.PARAMS['border']}")
-    dem_list, dem_source = get_topo_file((minlon, maxlon), (minlat, maxlat),
-                                         rgi_id=gdir.rgi_id,
-                                         dx_meter=dx,
-                                         source=source)
-    log.debug('(%s) DEM source: %s', gdir.rgi_id, dem_source)
-    log.debug('(%s) N DEM Files: %s', gdir.rgi_id, len(dem_list))
+        if rgi_id is None:
+            extent_string = f"the grid extent of longitudes {extent_ll[0]} and latitudes {extent_ll[1]}"
+        else:
+            extent_string = f"the glacier {rgi_id} with border {cfg.PARAMS['border']}"
+        raise InvalidWorkflowError(f'Source: {source} is not available for {extent_string}')
+    return source
+
+
+def reproject_dem(dem_list, dem_path, source, grid_data):
+    """
+    This function
+    Parameters
+    ----------
+    dem_list
+    dem_path
+    source
+    grid_data - holds necessary grid dimension data that is needed for reprojection
+
+    Returns
+    -------
+
+    """
 
     # Decide how to tag nodata
     def _get_nodata(rio_ds):
@@ -392,19 +368,23 @@ def _get_nodata(rio_ds):
 
     # Use Grid properties to create a transform (see rasterio cookbook)
     dst_transform = rasterio.transform.from_origin(
-        ulx, uly, dx, dx  # sign change (2nd dx) is done by rasterio.transform
+        grid_data['ulx'], grid_data['uly'], grid_data['dx'], grid_data['dx']
+        # sign change (2nd dx) is done by rasterio.transform
     )
 
     # Set up profile for writing output
     profile = dem_dss[0].profile
     profile.update({
-        'crs': utm_proj.srs,
+        'crs': grid_data['utm_proj'].srs,
         'transform': dst_transform,
         'nodata': nodata,
-        'width': nx,
-        'height': ny,
+        'width': grid_data['nx'],
+        'height': ['ny'],
         'driver': 'GTiff'
     })
+    profile.pop('blockxsize', None)
+    profile.pop('blockysize', None)
+    profile.pop('compress', None)
 
     # Could be extended so that the cfg file takes all Resampling.* methods
     if cfg.PARAMS['topo_interp'] == 'bilinear':
@@ -415,12 +395,9 @@ def _get_nodata(rio_ds):
         raise InvalidParamsError('{} interpolation not understood'
                                  .format(cfg.PARAMS['topo_interp']))
 
-    dem_reproj = gdir.get_filepath('dem')
-    profile.pop('blockxsize', None)
-    profile.pop('blockysize', None)
-    profile.pop('compress', None)
-    with rasterio.open(dem_reproj, 'w', **profile) as dest:
-        dst_array = np.empty((ny, nx), dtype=dem_dss[0].dtypes[0])
+
+    with rasterio.open(dem_path, 'w', **profile) as dest:
+        dst_array = np.empty((grid_data['ny'], grid_data['nx']), dtype=dem_dss[0].dtypes[0])
         reproject(
             # Source parameters
             source=dem_data,
@@ -430,7 +407,7 @@ def _get_nodata(rio_ds):
             # Destination parameters
             destination=dst_array,
             dst_transform=dst_transform,
-            dst_crs=utm_proj.srs,
+            dst_crs=grid_data['utm_proj'].srs,
             dst_nodata=nodata,
             # Configuration
             resampling=resampling)
@@ -439,6 +416,120 @@ def _get_nodata(rio_ds):
     for dem_ds in dem_dss:
         dem_ds.close()
 
+
+def dem_for_combined_grid(grid, fpath, source=None):
+    minlon, maxlon, minlat, maxlat = grid.extent_in_crs(crs=salem.wgs84)
+    extent_ll = [[minlon, maxlon], [minlat, maxlat]]
+    grid_data = {
+        'utm_proj': grid.proj,
+        'dx': grid.dx,
+        'ulx': grid.x0,
+        'uly': grid.y0,
+        'nx': grid.nx,
+        'ny': grid.ny
+    }
+
+    source = check_dem_source(source, extent_ll)
+
+    dem_list, dem_source = get_topo_file((minlon, maxlon), (minlat, maxlat),
+                                         dx_meter=grid_data['dx'],
+                                         source=source)
+    log.debug('(%s) DEM source: %s', f'extent(lon, lat) {extent_ll}', dem_source)
+    log.debug('(%s) N DEM Files: %s', f'extent(lon, lat) {extent_ll}', len(dem_list))
+
+    # further checks if given fpath exists?
+    dem_path = os.path.join(fpath, 'dem.tif')
+    reproject_dem(dem_list, fpath, source, grid_data)
+
+    """
+    currently not needed for the multiple glacier visualisation task.
+    
+    # Glacier grid
+    x0y0 = (ulx + dx / 2, uly - dx / 2)  # To pixel center coordinates
+    glacier_grid = salem.Grid(proj=utm_proj, nxny=(nx, ny), dxdy=(dx, -dx),
+                              x0y0=x0y0)
+    glacier_grid.to_json(gdir.get_filepath('glacier_grid'))
+
+    # Write DEM source info
+    gdir.add_to_diagnostics('dem_source', dem_source)
+    source_txt = DEM_SOURCE_INFO.get(dem_source, dem_source)
+    with open(gdir.get_filepath('dem_source'), 'w') as fw:
+        fw.write(source_txt)
+        fw.write('\n\n')
+        fw.write('# Data files\n\n')
+        for fname in dem_list:
+            fw.write('{}\n'.format(os.path.basename(fname)))
+    """
+
+@entity_task(log, writes=['glacier_grid', 'dem', 'outlines'])
+def define_glacier_region(gdir, entity=None, source=None):
+    """Very first task after initialization: define the glacier's local grid.
+
+    Defines the local projection (Transverse Mercator or UTM depending on
+    user choice), centered on the glacier.
+    There is some options to set the resolution of the local grid.
+    It can be adapted depending on the size of the glacier.
+    See ``params.cfg`` for setting these options.
+
+    Default values of the adapted mode lead to a resolution of 50 m for
+    Hintereisferner, which is approx. 8 km2 large.
+
+    After defining the grid, the topography and the outlines of the glacier
+    are transformed into the local projection.
+
+    Parameters
+    ----------
+    gdir : :py:class:`oggm.GlacierDirectory`
+        where to write the data
+    entity : geopandas.GeoSeries
+        the glacier geometry to process - DEPRECATED. It is now ignored
+    source : str or list of str, optional
+        If you want to force the use of a certain DEM source. Available are:
+          - 'USER' : file set in cfg.PATHS['dem_file']
+          - 'SRTM' : http://srtm.csi.cgiar.org/
+          - 'GIMP' : https://bpcrc.osu.edu/gdg/data/gimpdem
+          - 'RAMP' : https://nsidc.org/data/nsidc-0082/versions/2/documentation
+          - 'REMA' : https://www.pgc.umn.edu/data/rema/
+          - 'DEM3' : http://viewfinderpanoramas.org/
+          - 'ASTER' : https://asterweb.jpl.nasa.gov/gdem.asp
+          - 'TANDEM' : https://geoservice.dlr.de/web/dataguide/tdm90/
+          - 'ARCTICDEM' : https://www.pgc.umn.edu/data/arcticdem/
+          - 'AW3D30' : https://www.eorc.jaxa.jp
+          - 'MAPZEN' : https://registry.opendata.aws/terrain-tiles/
+          - 'ALASKA' : https://www.the-cryosphere.net/8/503/2014/
+          - 'COPDEM30' : Copernicus DEM GLO30 https://spacedata.copernicus.eu/web/cscda/cop-dem-faq
+          - 'COPDEM90' : Copernicus DEM GLO90 https://spacedata.copernicus.eu/web/cscda/cop-dem-faq
+          - 'NASADEM':  https://doi.org/10.5069/G93T9FD9
+    """
+
+    utm_proj, nx, ny, ulx, uly, dx = glacier_grid_params(gdir)
+
+    # Back to lon, lat for DEM download/preparation
+    tmp_grid = salem.Grid(proj=utm_proj, nxny=(nx, ny), x0y0=(ulx, uly),
+                          dxdy=(dx, -dx), pixel_ref='corner')
+    minlon, maxlon, minlat, maxlat = tmp_grid.extent_in_crs(crs=salem.wgs84)
+
+    # Open DEM
+    source = check_dem_source(source, gdir.extent_ll, rgi_id=gdir.rgi_id)
+
+    dem_list, dem_source = get_topo_file((minlon, maxlon), (minlat, maxlat),
+                                         rgi_id=gdir.rgi_id,
+                                         dx_meter=dx,
+                                         source=source)
+    log.debug('(%s) DEM source: %s', gdir.rgi_id, dem_source)
+    log.debug('(%s) N DEM Files: %s', gdir.rgi_id, len(dem_list))
+
+    grid_data = {
+        'utm_proj': utm_proj,
+        'dx': dx,
+        'ulx': ulx,
+        'uly': uly,
+        'nx': nx,
+        'ny': ny
+    }
+
+    reproject_dem(dem_list, gdir.get_filepath('dem'), source, grid_data)
+
     # Glacier grid
     x0y0 = (ulx+dx/2, uly-dx/2)  # To pixel center coordinates
     glacier_grid = salem.Grid(proj=utm_proj, nxny=(nx, ny), dxdy=(dx, -dx),
@@ -531,6 +622,50 @@ def read_geotiff_dem(gdir):
     return topo
 
 
+def prepareNcdfFile(ncdf_object):
+    """
+    This function takes care of the basic assembly of a gridded netCDF file, before other parameters can be added to it.
+    Parameters
+    ----------
+    ncdf_object
+
+    Returns
+    -------
+
+    """
+    if os.path.exists(ncdf_object.fpath):
+        # Already there - just append
+        ncdf_object.nc = ncDataset(ncdf_object.fpath, 'a', format='NETCDF4')
+        return ncdf_object.nc
+
+    # Create and fill
+    nc = ncDataset(ncdf_object.fpath, 'w', format='NETCDF4')
+
+    nc.createDimension('x', ncdf_object.grid.nx)
+    nc.createDimension('y', ncdf_object.grid.ny)
+
+    nc.author = 'OGGM'
+    nc.author_info = 'Open Global Glacier Model'
+    nc.pyproj_srs = ncdf_object.grid.proj.srs
+
+    x = ncdf_object.grid.x0 + np.arange(ncdf_object.grid.nx) * ncdf_object.grid.dx
+    y = ncdf_object.grid.y0 + np.arange(ncdf_object.grid.ny) * ncdf_object.grid.dy
+
+    v = nc.createVariable('x', 'f4', ('x',), zlib=True)
+    v.units = 'm'
+    v.long_name = 'x coordinate of projection'
+    v.standard_name = 'projection_x_coordinate'
+    v[:] = x
+
+    v = nc.createVariable('y', 'f4', ('y',), zlib=True)
+    v.units = 'm'
+    v.long_name = 'y coordinate of projection'
+    v.standard_name = 'projection_y_coordinate'
+    v[:] = y
+
+    ncdf_object.nc = nc
+    return nc
+
 class GriddedNcdfFile(object):
     """Creates or opens a gridded netcdf file template.
 
@@ -544,39 +679,37 @@ def __init__(self, gdir, basename='gridded_data', reset=False):
             os.remove(self.fpath)
 
     def __enter__(self):
+        return prepareNcdfFile(self)
 
-        if os.path.exists(self.fpath):
-            # Already there - just append
-            self.nc = ncDataset(self.fpath, 'a', format='NETCDF4')
-            return self.nc
-
-        # Create and fill
-        nc = ncDataset(self.fpath, 'w', format='NETCDF4')
-
-        nc.createDimension('x', self.grid.nx)
-        nc.createDimension('y', self.grid.ny)
+    def __exit__(self, exc_type, exc_value, exc_traceback):
+        self.nc.close()
 
-        nc.author = 'OGGM'
-        nc.author_info = 'Open Global Glacier Model'
-        nc.pyproj_srs = self.grid.proj.srs
 
-        x = self.grid.x0 + np.arange(self.grid.nx) * self.grid.dx
-        y = self.grid.y0 + np.arange(self.grid.ny) * self.grid.dy
+class CombinedNcdfFile(object):
+    """
+    Creates or opens a gridded netcdf file template independent of the gdir-logic.
 
-        v = nc.createVariable('x', 'f4', ('x',), zlib=True)
-        v.units = 'm'
-        v.long_name = 'x coordinate of projection'
-        v.standard_name = 'projection_x_coordinate'
-        v[:] = x
+    The other variables have to be created and filled by the calling
+    routine.
+    """
 
-        v = nc.createVariable('y', 'f4', ('y',), zlib=True)
-        v.units = 'm'
-        v.long_name = 'y coordinate of projection'
-        v.standard_name = 'projection_y_coordinate'
-        v[:] = y
+    def __init__(self, combined_grid, fpath, basename='gridded_data', reset=False):
+        """
+        Parameters
+        ----------
+        combined_grid: an already combined salem.Grid of several glaciers.
+        fpath: file path to a directory where the netCDF file of a combined grid is stored.
+        basename: name of the output .nc file
+
+        reset
+        """
+        self.grid = combined_grid
+        self.fpath = os.path.join(fpath, basename)
+        if reset and os.path.exists(self.fpath):
+            os.remove(self.fpath)
 
-        self.nc = nc
-        return nc
+    def __enter__(self):
+        return prepareNcdfFile(self)
 
     def __exit__(self, exc_type, exc_value, exc_traceback):
         self.nc.close()

From 874b2d0e8cde7e6e37126970131471dd912efb85 Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Mon, 27 Nov 2023 16:13:53 +0100
Subject: [PATCH 04/16] bugfix

---
 oggm/core/gis.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/oggm/core/gis.py b/oggm/core/gis.py
index 8b39d7e0e..03afb3041 100644
--- a/oggm/core/gis.py
+++ b/oggm/core/gis.py
@@ -379,7 +379,7 @@ def _get_nodata(rio_ds):
         'transform': dst_transform,
         'nodata': nodata,
         'width': grid_data['nx'],
-        'height': ['ny'],
+        'height': grid_data['ny'],
         'driver': 'GTiff'
     })
     profile.pop('blockxsize', None)
@@ -395,7 +395,6 @@ def _get_nodata(rio_ds):
         raise InvalidParamsError('{} interpolation not understood'
                                  .format(cfg.PARAMS['topo_interp']))
 
-
     with rasterio.open(dem_path, 'w', **profile) as dest:
         dst_array = np.empty((grid_data['ny'], grid_data['nx']), dtype=dem_dss[0].dtypes[0])
         reproject(

From f764bab4088f7141cbf72551fca351882c631927 Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Tue, 28 Nov 2023 10:21:31 +0100
Subject: [PATCH 05/16] bugfix

---
 oggm/core/gis.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/oggm/core/gis.py b/oggm/core/gis.py
index 03afb3041..eca253437 100644
--- a/oggm/core/gis.py
+++ b/oggm/core/gis.py
@@ -438,7 +438,7 @@ def dem_for_combined_grid(grid, fpath, source=None):
 
     # further checks if given fpath exists?
     dem_path = os.path.join(fpath, 'dem.tif')
-    reproject_dem(dem_list, fpath, source, grid_data)
+    reproject_dem(dem_list, dem_path, source, grid_data)
 
     """
     currently not needed for the multiple glacier visualisation task.

From 628ef97bad2518b65168945c7117b069dd363f19 Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Tue, 28 Nov 2023 10:21:54 +0100
Subject: [PATCH 06/16] added comments to the combine_distributed_thickness
 function

---
 oggm/sandbox/distribute_2d.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py
index 44cdfe473..dd1d23e7a 100644
--- a/oggm/sandbox/distribute_2d.py
+++ b/oggm/sandbox/distribute_2d.py
@@ -434,15 +434,16 @@ def combine_distributed_thickness(gdirs, output_folder, source=None):
     # create a combined salem.Grid object, which serves as canvas/boundaries of the combined glacier region
     combined_grid = utils.combine_grids(gdirs)
 
-    #create gridded_data.nc
+    #create gridded_data.nc - maybe use this further down, where we can already add the Topography to the gridded_data
     with CombinedNcdfFile(combined_grid, output_folder) as nc:
-        # is there still something else needed in the gridded_data.nc? add it here or use function
+        # need to check which other data is necessary for the plotting of the distributed_thickness? or which other
+        # data should be included? glacier outlines etc.?
         pass
 
     # retrieve and project a topography/DEM for the defined region
     dem_for_combined_grid(combined_grid, output_folder, source=source)
 
-    # add 'write_dem_to_gridded_data' function and use it here
+    # add 'write_dem_to_gridded_data' function and add dem to the gridded data
 
 
     # add individual add individual distributed thicknesses to gridded data

From 83781d450b98e80bf83ab532094fca3f66f93f6e Mon Sep 17 00:00:00 2001
From: Patrick <patrick.schmitt@uibk.ac.at>
Date: Tue, 28 Nov 2023 13:20:48 +0100
Subject: [PATCH 07/16] cleaning and refractoring

---
 oggm/core/gis.py              | 299 +++++++++++++++++-----------------
 oggm/sandbox/distribute_2d.py |  31 ++--
 oggm/utils/_funcs.py          |   2 +
 3 files changed, 173 insertions(+), 159 deletions(-)

diff --git a/oggm/core/gis.py b/oggm/core/gis.py
index eca253437..3d81a475c 100644
--- a/oggm/core/gis.py
+++ b/oggm/core/gis.py
@@ -297,22 +297,33 @@ def glacier_grid_params(gdir):
 
     return utm_proj, nx, ny, ulx, uly, dx
 
+
 def check_dem_source(source, extent_ll, rgi_id=None):
     """
-    This function can check for multiple DEM sources and is in charge of the error handling if a requested source is
-    not available for the given glacier/extent
+    This function can check for multiple DEM sources and is in charge of the
+    error handling if a requested source is not available for the given
+    glacier/extent
+
     Parameters
     ----------
-    source
-    extent_ll
-    gdir
+    source : str or list of str
+        If you want to force the use of a certain DEM source. If list is
+        provided they are checked in order. For a list of available options
+        check docstring of oggm.core.gis.define_glacier_region.
+    extent_ll : list
+        The longitude and latitude extend which should be checked for. Should
+        be provided as extent_ll = [[minlon, maxlon], [minlat, maxlat]].
+    rgi_id : str
+        The RGI-ID of the glacier, only used for a descriptive error message in
+        case.
 
     Returns
     -------
-
+    String of the working DEM source.
     """
-# We test DEM availability for glacier only (maps can grow big)
-    if isinstance(source, list):  # when multiple sources are provided, try them sequentially
+
+    # when multiple sources are provided, try them sequentially
+    if isinstance(source, list):
         for src in source:
             source_exists = is_dem_source_available(src, *extent_ll)
             if source_exists:
@@ -322,26 +333,36 @@ def check_dem_source(source, extent_ll, rgi_id=None):
         source_exists = is_dem_source_available(source, *extent_ll)
     if not source_exists:
         if rgi_id is None:
-            extent_string = f"the grid extent of longitudes {extent_ll[0]} and latitudes {extent_ll[1]}"
+            extent_string = (f"the grid extent of longitudes {extent_ll[0]} "
+                             f"and latitudes {extent_ll[1]}")
         else:
-            extent_string = f"the glacier {rgi_id} with border {cfg.PARAMS['border']}"
-        raise InvalidWorkflowError(f'Source: {source} is not available for {extent_string}')
+            extent_string = (f"the glacier {rgi_id} with border "
+                             f"{cfg.PARAMS['border']}")
+        raise InvalidWorkflowError(f"Source: {source} is not available for "
+                                   f"{extent_string}")
     return source
 
 
-def reproject_dem(dem_list, dem_path, source, grid_data):
+def reproject_dem(dem_list, dem_source, dst_grid_prop, output_path):
     """
-    This function
+    This function reprojects a provided DEM to the destination grid and saves
+    the result to disk.
+
     Parameters
     ----------
-    dem_list
-    dem_path
-    source
-    grid_data - holds necessary grid dimension data that is needed for reprojection
+    dem_list : list of str
+        list with path(s) to the DEM file(s)
+    dem_source : str
+        DEM source string
+    dst_grid_prop : dict
+        Holds necessary grid properties. Must contain 'utm_proj', 'dx', 'ulx',
+        'uly', 'nx' and 'ny.
+    output_path : str
+        Filepath where to store the reprojected DEM.
 
     Returns
     -------
-
+    None
     """
 
     # Decide how to tag nodata
@@ -349,7 +370,7 @@ def _get_nodata(rio_ds):
         nodata = rio_ds[0].meta.get('nodata', None)
         if nodata is None:
             # badly tagged geotiffs, let's do it ourselves
-            nodata = -32767 if source == 'TANDEM' else -9999
+            nodata = -32767 if dem_source == 'TANDEM' else -9999
         return nodata
 
     # A glacier area can cover more than one tile:
@@ -368,18 +389,19 @@ def _get_nodata(rio_ds):
 
     # Use Grid properties to create a transform (see rasterio cookbook)
     dst_transform = rasterio.transform.from_origin(
-        grid_data['ulx'], grid_data['uly'], grid_data['dx'], grid_data['dx']
+        dst_grid_prop['ulx'], dst_grid_prop['uly'], dst_grid_prop['dx'],
+        dst_grid_prop['dx']
         # sign change (2nd dx) is done by rasterio.transform
     )
 
     # Set up profile for writing output
     profile = dem_dss[0].profile
     profile.update({
-        'crs': grid_data['utm_proj'].srs,
+        'crs': dst_grid_prop['utm_proj'].srs,
         'transform': dst_transform,
         'nodata': nodata,
-        'width': grid_data['nx'],
-        'height': grid_data['ny'],
+        'width': dst_grid_prop['nx'],
+        'height': dst_grid_prop['ny'],
         'driver': 'GTiff'
     })
     profile.pop('blockxsize', None)
@@ -395,8 +417,9 @@ def _get_nodata(rio_ds):
         raise InvalidParamsError('{} interpolation not understood'
                                  .format(cfg.PARAMS['topo_interp']))
 
-    with rasterio.open(dem_path, 'w', **profile) as dest:
-        dst_array = np.empty((grid_data['ny'], grid_data['nx']), dtype=dem_dss[0].dtypes[0])
+    with rasterio.open(output_path, 'w', **profile) as dest:
+        dst_array = np.empty((dst_grid_prop['ny'], dst_grid_prop['nx']),
+                             dtype=dem_dss[0].dtypes[0])
         reproject(
             # Source parameters
             source=dem_data,
@@ -406,7 +429,7 @@ def _get_nodata(rio_ds):
             # Destination parameters
             destination=dst_array,
             dst_transform=dst_transform,
-            dst_crs=grid_data['utm_proj'].srs,
+            dst_crs=dst_grid_prop['utm_proj'].srs,
             dst_nodata=nodata,
             # Configuration
             resampling=resampling)
@@ -416,10 +439,31 @@ def _get_nodata(rio_ds):
         dem_ds.close()
 
 
-def dem_for_combined_grid(grid, fpath, source=None):
+def get_dem_for_grid(grid, fpath, source=None, rgi_id=None):
+    """
+    Fetch a DEM from source, reproject it to the extent defined by grid and
+    saves it to disk.
+
+    Parameters
+    ----------
+    grid : :py:class:`salem.gis.Grid`
+        Grid which defines the extent and projection of the final DEM
+    fpath : str
+        The output filepath for the final DEM.
+    source : str or list of str
+        If you want to force the use of a certain DEM source. If list is
+        provided they are checked in order. For a list of available options
+        check docstring of oggm.core.gis.define_glacier_region.
+    rgi_id : str
+        The RGI-ID of the glacier.
+
+    Returns
+    -------
+    tuple: (list with path(s) to the DEM file(s), data source str)
+    """
     minlon, maxlon, minlat, maxlat = grid.extent_in_crs(crs=salem.wgs84)
     extent_ll = [[minlon, maxlon], [minlat, maxlat]]
-    grid_data = {
+    grid_prop = {
         'utm_proj': grid.proj,
         'dx': grid.dx,
         'ulx': grid.x0,
@@ -428,37 +472,28 @@ def dem_for_combined_grid(grid, fpath, source=None):
         'ny': grid.ny
     }
 
-    source = check_dem_source(source, extent_ll)
+    source = check_dem_source(source, extent_ll, rgi_id=rgi_id)
 
     dem_list, dem_source = get_topo_file((minlon, maxlon), (minlat, maxlat),
-                                         dx_meter=grid_data['dx'],
+                                         rgi_id=rgi_id,
+                                         dx_meter=grid_prop['dx'],
                                          source=source)
-    log.debug('(%s) DEM source: %s', f'extent(lon, lat) {extent_ll}', dem_source)
-    log.debug('(%s) N DEM Files: %s', f'extent(lon, lat) {extent_ll}', len(dem_list))
+
+    if rgi_id is not None:
+        log.debug('(%s) DEM source: %s', rgi_id, dem_source)
+        log.debug('(%s) N DEM Files: %s', rgi_id, len(dem_list))
 
     # further checks if given fpath exists?
-    dem_path = os.path.join(fpath, 'dem.tif')
-    reproject_dem(dem_list, dem_path, source, grid_data)
+    if fpath[-4:] != '.tif':
+        # we add the default filename
+        fpath = os.path.join(fpath, 'dem.tif')
 
-    """
-    currently not needed for the multiple glacier visualisation task.
-    
-    # Glacier grid
-    x0y0 = (ulx + dx / 2, uly - dx / 2)  # To pixel center coordinates
-    glacier_grid = salem.Grid(proj=utm_proj, nxny=(nx, ny), dxdy=(dx, -dx),
-                              x0y0=x0y0)
-    glacier_grid.to_json(gdir.get_filepath('glacier_grid'))
+    reproject_dem(dem_list=dem_list, dem_source=dem_source,
+                  dst_grid_prop=grid_prop,
+                  output_path=fpath)
+
+    return dem_list, dem_source
 
-    # Write DEM source info
-    gdir.add_to_diagnostics('dem_source', dem_source)
-    source_txt = DEM_SOURCE_INFO.get(dem_source, dem_source)
-    with open(gdir.get_filepath('dem_source'), 'w') as fw:
-        fw.write(source_txt)
-        fw.write('\n\n')
-        fw.write('# Data files\n\n')
-        for fname in dem_list:
-            fw.write('{}\n'.format(os.path.basename(fname)))
-    """
 
 @entity_task(log, writes=['glacier_grid', 'dem', 'outlines'])
 def define_glacier_region(gdir, entity=None, source=None):
@@ -506,28 +541,10 @@ def define_glacier_region(gdir, entity=None, source=None):
     # Back to lon, lat for DEM download/preparation
     tmp_grid = salem.Grid(proj=utm_proj, nxny=(nx, ny), x0y0=(ulx, uly),
                           dxdy=(dx, -dx), pixel_ref='corner')
-    minlon, maxlon, minlat, maxlat = tmp_grid.extent_in_crs(crs=salem.wgs84)
 
-    # Open DEM
-    source = check_dem_source(source, gdir.extent_ll, rgi_id=gdir.rgi_id)
-
-    dem_list, dem_source = get_topo_file((minlon, maxlon), (minlat, maxlat),
-                                         rgi_id=gdir.rgi_id,
-                                         dx_meter=dx,
-                                         source=source)
-    log.debug('(%s) DEM source: %s', gdir.rgi_id, dem_source)
-    log.debug('(%s) N DEM Files: %s', gdir.rgi_id, len(dem_list))
-
-    grid_data = {
-        'utm_proj': utm_proj,
-        'dx': dx,
-        'ulx': ulx,
-        'uly': uly,
-        'nx': nx,
-        'ny': ny
-    }
-
-    reproject_dem(dem_list, gdir.get_filepath('dem'), source, grid_data)
+    dem_list, dem_source = get_dem_for_grid(grid=tmp_grid,
+                                            fpath=gdir.get_filepath('dem'),
+                                            source=source, rgi_id=gdir.rgi_id)
 
     # Glacier grid
     x0y0 = (ulx+dx/2, uly-dx/2)  # To pixel center coordinates
@@ -621,94 +638,80 @@ def read_geotiff_dem(gdir):
     return topo
 
 
-def prepareNcdfFile(ncdf_object):
-    """
-    This function takes care of the basic assembly of a gridded netCDF file, before other parameters can be added to it.
-    Parameters
-    ----------
-    ncdf_object
-
-    Returns
-    -------
-
-    """
-    if os.path.exists(ncdf_object.fpath):
-        # Already there - just append
-        ncdf_object.nc = ncDataset(ncdf_object.fpath, 'a', format='NETCDF4')
-        return ncdf_object.nc
-
-    # Create and fill
-    nc = ncDataset(ncdf_object.fpath, 'w', format='NETCDF4')
-
-    nc.createDimension('x', ncdf_object.grid.nx)
-    nc.createDimension('y', ncdf_object.grid.ny)
-
-    nc.author = 'OGGM'
-    nc.author_info = 'Open Global Glacier Model'
-    nc.pyproj_srs = ncdf_object.grid.proj.srs
-
-    x = ncdf_object.grid.x0 + np.arange(ncdf_object.grid.nx) * ncdf_object.grid.dx
-    y = ncdf_object.grid.y0 + np.arange(ncdf_object.grid.ny) * ncdf_object.grid.dy
-
-    v = nc.createVariable('x', 'f4', ('x',), zlib=True)
-    v.units = 'm'
-    v.long_name = 'x coordinate of projection'
-    v.standard_name = 'projection_x_coordinate'
-    v[:] = x
-
-    v = nc.createVariable('y', 'f4', ('y',), zlib=True)
-    v.units = 'm'
-    v.long_name = 'y coordinate of projection'
-    v.standard_name = 'projection_y_coordinate'
-    v[:] = y
-
-    ncdf_object.nc = nc
-    return nc
-
 class GriddedNcdfFile(object):
     """Creates or opens a gridded netcdf file template.
 
     The other variables have to be created and filled by the calling
     routine.
     """
-    def __init__(self, gdir, basename='gridded_data', reset=False):
-        self.fpath = gdir.get_filepath(basename)
-        self.grid = gdir.grid
+    def __init__(self, gdir=None, grid=None, fpath=None,
+                 basename='gridded_data', reset=False):
+        """
+        Parameters
+        ----------
+        gdir : :py:class:`oggm.GlacierDirectory`
+            The glacier directory. If provided it defines the filepath and the
+            grid used for the netcdf file. It overrules grid and fpath if all
+            are provided.
+        grid : :py:class:`salem.gis.Grid`
+            Grid which defines the extend of the netcdf file. Needed if gdir is
+            not provided.
+        fpath : str
+            The output filepath for the netcdf file. Needed if gdir is not
+            provided.
+        basename : str
+            The filename of the resulting netcdf file
+        reset : bool
+            If True, a potentially existing file will be deleted.
+        """
+
+        if gdir is not None:
+            self.fpath = gdir.get_filepath(basename)
+            self.grid = gdir.grid
+        else:
+            if grid is None or fpath is None:
+                raise InvalidParamsError('If you do not provide a gdir you must'
+                                         'define grid and fpath! Given grid='
+                                         f'{grid} and fpath={fpath}.')
+            self.fpath = os.path.join(fpath, basename)
+            self.grid = grid
+
         if reset and os.path.exists(self.fpath):
             os.remove(self.fpath)
 
     def __enter__(self):
-        return prepareNcdfFile(self)
+        if os.path.exists(self.fpath):
+            # Already there - just append
+            self.nc = ncDataset(self.fpath, 'a', format='NETCDF4')
+            return self.nc
 
-    def __exit__(self, exc_type, exc_value, exc_traceback):
-        self.nc.close()
+        # Create and fill
+        nc = ncDataset(self.fpath, 'w', format='NETCDF4')
 
+        nc.createDimension('x', self.grid.nx)
+        nc.createDimension('y', self.grid.ny)
 
-class CombinedNcdfFile(object):
-    """
-    Creates or opens a gridded netcdf file template independent of the gdir-logic.
+        nc.author = 'OGGM'
+        nc.author_info = 'Open Global Glacier Model'
+        nc.pyproj_srs = self.grid.proj.srs
 
-    The other variables have to be created and filled by the calling
-    routine.
-    """
+        x = self.grid.x0 + np.arange(self.grid.nx) * self.grid.dx
+        y = self.grid.y0 + np.arange(self.grid.ny) * self.grid.dy
 
-    def __init__(self, combined_grid, fpath, basename='gridded_data', reset=False):
-        """
-        Parameters
-        ----------
-        combined_grid: an already combined salem.Grid of several glaciers.
-        fpath: file path to a directory where the netCDF file of a combined grid is stored.
-        basename: name of the output .nc file
+        v = nc.createVariable('x', 'f4', ('x',), zlib=True)
+        v.units = 'm'
+        v.long_name = 'x coordinate of projection'
+        v.standard_name = 'projection_x_coordinate'
+        v[:] = x
 
-        reset
-        """
-        self.grid = combined_grid
-        self.fpath = os.path.join(fpath, basename)
-        if reset and os.path.exists(self.fpath):
-            os.remove(self.fpath)
+        v = nc.createVariable('y', 'f4', ('y',), zlib=True)
+        v.units = 'm'
+        v.long_name = 'y coordinate of projection'
+        v.standard_name = 'projection_y_coordinate'
+        v[:] = y
 
-    def __enter__(self):
-        return prepareNcdfFile(self)
+        self.nc = nc
+        return nc
 
     def __exit__(self, exc_type, exc_value, exc_traceback):
         self.nc.close()
@@ -794,7 +797,7 @@ def process_dem(gdir):
         utils.clip_min(smoothed_dem, 0, out=smoothed_dem)
 
     # Write to file
-    with GriddedNcdfFile(gdir, reset=True) as nc:
+    with GriddedNcdfFile(gdir=gdir, reset=True) as nc:
 
         v = nc.createVariable('topo', 'f4', ('y', 'x',), zlib=True)
         v.units = 'm'
@@ -903,7 +906,7 @@ def proj(x, y):
     gdir.write_pickle(geometries, 'geometries')
 
     # write out the grids in the netcdf file
-    with GriddedNcdfFile(gdir) as nc:
+    with GriddedNcdfFile(gdir=gdir) as nc:
 
         if 'glacier_mask' not in nc.variables:
             v = nc.createVariable('glacier_mask', 'i1', ('y', 'x', ),
@@ -1037,7 +1040,7 @@ def simple_glacier_masks(gdir):
                               .format(gdir.rgi_id))
 
     # write out the grids in the netcdf file
-    with GriddedNcdfFile(gdir) as nc:
+    with GriddedNcdfFile(gdir=gdir) as nc:
 
         if 'glacier_mask' not in nc.variables:
             v = nc.createVariable('glacier_mask', 'i1', ('y', 'x', ),
@@ -1751,7 +1754,7 @@ def project(x, y):
                            .format(gdir.rgi_id))
 
     # write out the grids in the netcdf file
-    with GriddedNcdfFile(gdir, reset=True) as nc:
+    with GriddedNcdfFile(gdir=gdir, reset=True) as nc:
 
         v = nc.createVariable('topo', 'f4', ('y', 'x', ), zlib=True)
         v.units = 'm'
diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py
index dd1d23e7a..13b9f0ef6 100644
--- a/oggm/sandbox/distribute_2d.py
+++ b/oggm/sandbox/distribute_2d.py
@@ -9,7 +9,7 @@
 import xarray as xr
 from scipy import ndimage
 from scipy.stats import mstats
-from oggm.core.gis import gaussian_blur, dem_for_combined_grid, CombinedNcdfFile
+from oggm.core.gis import gaussian_blur, get_dem_for_grid, GriddedNcdfFile
 from oggm.utils import ncDataset, entity_task
 import matplotlib.pyplot as plt
 
@@ -417,31 +417,40 @@ def distribute_thickness_from_simulation(gdir,
 
     return ds
 
+
 def combine_distributed_thickness(gdirs, output_folder, source=None):
     """
-    This function takes a list of glacier directories that have a distributed_thickness
+    This function takes a list of glacier directories that have a
+    distributed_thickness.
+
     Parameters
     ----------
-    gdirs - list of GlacierDirectories
-    output_folder - folder where the intermediate files and the final combined distributed thickness
-                    will be stored
-    source - DEM source
+    gdirs : list of :py:class:`oggm.GlacierDirectory` objects
+        The glacier directories which should be combined
+    output_folder : str
+        Folder where the intermediate files and the final combined gridded data
+        should be stored
+    source : str
+        DEM source
 
     Returns
     -------
 
     """
-    # create a combined salem.Grid object, which serves as canvas/boundaries of the combined glacier region
+    # create a combined salem.Grid object, which serves as canvas/boundaries of
+    # the combined glacier region
     combined_grid = utils.combine_grids(gdirs)
 
-    #create gridded_data.nc - maybe use this further down, where we can already add the Topography to the gridded_data
-    with CombinedNcdfFile(combined_grid, output_folder) as nc:
-        # need to check which other data is necessary for the plotting of the distributed_thickness? or which other
+    #create gridded_data.nc - maybe use this further down, where we can already
+    # add the Topography to the gridded_data
+    with GriddedNcdfFile(grid=combined_grid, fpath=output_folder) as nc:
+        # need to check which other data is necessary for the plotting of the
+        # distributed_thickness? or which other
         # data should be included? glacier outlines etc.?
         pass
 
     # retrieve and project a topography/DEM for the defined region
-    dem_for_combined_grid(combined_grid, output_folder, source=source)
+    get_dem_for_grid(combined_grid, output_folder, source=source)
 
     # add 'write_dem_to_gridded_data' function and add dem to the gridded data
 
diff --git a/oggm/utils/_funcs.py b/oggm/utils/_funcs.py
index beecd63d1..e3b13c729 100644
--- a/oggm/utils/_funcs.py
+++ b/oggm/utils/_funcs.py
@@ -539,6 +539,7 @@ def recursive_valid_polygons(geoms, crs):
     assert np.all([type(geom) == shpg.Polygon for geom in new_geoms])
     return new_geoms
 
+
 def combine_grids(gdirs):
     """ Combines individual grids of different glacier directories to show
         multiple glaciers in the same plot. The resulting grid extent includes
@@ -611,6 +612,7 @@ def combine_grids(gdirs):
 
     return Grid.from_dict(new_grid)
 
+
 def multipolygon_to_polygon(geometry, gdir=None):
     """Sometimes an RGI geometry is a multipolygon: this should not happen.
 

From 6eedc1aed2c8c0838badda16b3f10e1b0bf0c00d Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Tue, 5 Dec 2023 11:29:44 +0100
Subject: [PATCH 08/16] added 'unwrapped' option to the entity_task decorator
 function

---
 oggm/utils/_workflow.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/oggm/utils/_workflow.py b/oggm/utils/_workflow.py
index 96a9b56c4..45822402a 100644
--- a/oggm/utils/_workflow.py
+++ b/oggm/utils/_workflow.py
@@ -515,6 +515,8 @@ def _entity_task(gdir, *, reset=None, print_log=True,
                 return out
 
         _entity_task.__dict__['is_entity_task'] = True
+        # adds the possibility to use a function, decorated as entity_task, without its decoration.
+        _entity_task.unwrapped = task_func
         return _entity_task
 
 

From 071150072d49d902f77ac1b50242a457941e2005 Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Tue, 5 Dec 2023 11:30:28 +0100
Subject: [PATCH 09/16] worked on 'combined_distributed_thickness'

---
 oggm/sandbox/distribute_2d.py | 35 ++++++++++++++++++-----------------
 1 file changed, 18 insertions(+), 17 deletions(-)

diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py
index 13b9f0ef6..a46c3816b 100644
--- a/oggm/sandbox/distribute_2d.py
+++ b/oggm/sandbox/distribute_2d.py
@@ -1,6 +1,7 @@
 import logging
 import warnings
 
+import os
 import oggm
 import oggm.cfg as cfg
 from oggm import utils
@@ -9,7 +10,7 @@
 import xarray as xr
 from scipy import ndimage
 from scipy.stats import mstats
-from oggm.core.gis import gaussian_blur, get_dem_for_grid, GriddedNcdfFile
+from oggm.core.gis import gaussian_blur, get_dem_for_grid, GriddedNcdfFile, process_dem
 from oggm.utils import ncDataset, entity_task
 import matplotlib.pyplot as plt
 
@@ -418,7 +419,7 @@ def distribute_thickness_from_simulation(gdir,
     return ds
 
 
-def combine_distributed_thickness(gdirs, output_folder, source=None):
+def combine_distributed_thickness(gdirs, output_folder, suffix='', source=None):
     """
     This function takes a list of glacier directories that have a
     distributed_thickness.
@@ -441,28 +442,28 @@ def combine_distributed_thickness(gdirs, output_folder, source=None):
     # the combined glacier region
     combined_grid = utils.combine_grids(gdirs)
 
-    #create gridded_data.nc - maybe use this further down, where we can already
-    # add the Topography to the gridded_data
-    with GriddedNcdfFile(grid=combined_grid, fpath=output_folder) as nc:
-        # need to check which other data is necessary for the plotting of the
-        # distributed_thickness? or which other
-        # data should be included? glacier outlines etc.?
-        pass
+    # check if output_folder exists
+    if not os.path.exists(output_folder):
+        os.makedirs(output_folder)
 
     # retrieve and project a topography/DEM for the defined region
     get_dem_for_grid(combined_grid, output_folder, source=source)
+    # processes dem(smoothing etc.) and creates gridded_data and adds DEM data
+    process_dem.unwrapped(gdir=None, grid=combined_grid, fpath=output_folder)
 
     # add 'write_dem_to_gridded_data' function and add dem to the gridded data
 
 
     # add individual add individual distributed thicknesses to gridded data
-    # for gdir in gdirs:
-    #     with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
-    #         ds = ds.load()
-    #
-    #     r_data = combined_grid.map_gridded_data(ds[f'simulation_distributed_thickness_{ssp}_projection'],
-    #                                             grid=gdir.grid)
-    #     total_data += r_data.filled(0)
-    # total_data = np.where(total_data == 0, np.NaN, total_data)
+    total_data = np.zeros((combined_grid.ny, combined_grid.nx))
+    for gdir in gdirs:
+        with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
+            ds = ds.load()
+
+        r_data = combined_grid.map_gridded_data(ds['distributed_thickness'],
+                                                grid=gdir.grid)
+        total_data += r_data.filled(0)
+    total_data = np.where(total_data == 0, np.NaN, total_data)
+    return total_data
 
 

From df7012425ea1a23fc54834d6c4d1e8086ae34f67 Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Tue, 5 Dec 2023 11:31:55 +0100
Subject: [PATCH 10/16] adapted read_geotiff_dem

---
 oggm/core/gis.py | 21 ++++++++++++++++++---
 1 file changed, 18 insertions(+), 3 deletions(-)

diff --git a/oggm/core/gis.py b/oggm/core/gis.py
index 3d81a475c..ec048268a 100644
--- a/oggm/core/gis.py
+++ b/oggm/core/gis.py
@@ -619,19 +619,34 @@ def rasterio_to_gdir(gdir, input_file, output_file_name,
                     dst.write(dest, indexes=i)
 
 
-def read_geotiff_dem(gdir):
-    """Reads (and masks out) the DEM out of the gdir's geotiff file.
+def read_geotiff_dem(gdir=None, fpath=None):
+    """Reads (and masks out) the DEM out of the gdir's geotiff file or a geotiff file given by the fpath variable.
 
     Parameters
     ----------
     gdir : :py:class:`oggm.GlacierDirectory`
         the glacier directory
+    fpath : 'str'
+        path to a gdir
 
     Returns
     -------
     2D np.float32 array
     """
-    with rasterio.open(gdir.get_filepath('dem'), 'r', driver='GTiff') as ds:
+    if gdir is not None:
+        dem_path = gdir.get_filepath('dem')
+    else:
+        if fpath is not None:
+            if fpath[-4:] != '.tif':
+                # we add the default filename
+                fpath = os.path.join(fpath, 'dem.tif')
+            dem_path = fpath
+        else:
+            raise InvalidParamsError('If you do not provide a gdir you must'
+                                     f'define a fpath! Given fpath={fpath}.')
+
+
+    with rasterio.open(dem_path, 'r', driver='GTiff') as ds:
         topo = ds.read(1).astype(rasterio.float32)
         topo[topo <= -999.] = np.NaN
         topo[ds.read_masks(1) == 0] = np.NaN

From 44658d6b2eb71ace5c15b9b03dca0d1eef6c2328 Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Tue, 5 Dec 2023 11:32:55 +0100
Subject: [PATCH 11/16] spelling mistake

---
 oggm/core/gis.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/oggm/core/gis.py b/oggm/core/gis.py
index ec048268a..55f6132fc 100644
--- a/oggm/core/gis.py
+++ b/oggm/core/gis.py
@@ -669,7 +669,7 @@ def __init__(self, gdir=None, grid=None, fpath=None,
             grid used for the netcdf file. It overrules grid and fpath if all
             are provided.
         grid : :py:class:`salem.gis.Grid`
-            Grid which defines the extend of the netcdf file. Needed if gdir is
+            Grid which defines the extent of the netcdf file. Needed if gdir is
             not provided.
         fpath : str
             The output filepath for the netcdf file. Needed if gdir is not

From 32b0a2f49964d12129a4e5b0ae1c467bd326cf97 Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Tue, 5 Dec 2023 11:35:55 +0100
Subject: [PATCH 12/16] adapted process_dem

---
 oggm/core/gis.py | 66 +++++++++++++++++++++++++++++++++---------------
 1 file changed, 46 insertions(+), 20 deletions(-)

diff --git a/oggm/core/gis.py b/oggm/core/gis.py
index 55f6132fc..c2d29aced 100644
--- a/oggm/core/gis.py
+++ b/oggm/core/gis.py
@@ -6,6 +6,7 @@
 import os
 import logging
 import warnings
+import json
 from packaging.version import Version
 from functools import partial
 
@@ -733,7 +734,7 @@ def __exit__(self, exc_type, exc_value, exc_traceback):
 
 
 @entity_task(log, writes=['gridded_data'])
-def process_dem(gdir):
+def process_dem(gdir=None, grid=None, fpath=None):
     """Reads the DEM from the tiff, attempts to fill voids and apply smooth.
 
     The data is then written to `gridded_data.nc`.
@@ -743,13 +744,39 @@ def process_dem(gdir):
     gdir : :py:class:`oggm.GlacierDirectory`
         where to write the data
     """
-
-    # open srtm tif-file:
-    dem = read_geotiff_dem(gdir)
-
-    # Grid
-    nx = gdir.grid.nx
-    ny = gdir.grid.ny
+    if gdir is not None:
+        # open srtm tif-file:
+        dem = read_geotiff_dem(gdir)
+        # Grid
+        dem_grid = gdir.grid
+        grid_name = gdir.rgi_id
+    else:
+        if grid is None or fpath is None:
+            raise InvalidParamsError('If you do not provide a gdir you must'
+                                     'define grid and fpath! Given grid='
+                                     f'{grid} and fpath={fpath}.')
+        dem = read_geotiff_dem(fpath=fpath)
+        grid_name = 'custom_grid'
+        dem_grid = grid
+        diagnostics_dict = dict()
+    # Grid parameters
+    nx = dem_grid.nx
+    ny = dem_grid.ny
+    dx = dem_grid.dx
+    xx, yy = dem_grid.ij_coordinates
+
+    def _log_and_diagnostics(pnan, log_string='extra'):
+        if log_string != 'extra':
+            perc_string = 'dem_invalid_perc'
+        else:
+            perc_string = 'dem_extrapol_perc'
+        log.info(grid_name + f': DEM needed {log_string}polation.')
+        if gdir is not None:
+            gdir.add_to_diagnostics(f'dem_needed_{log_string}polation', True)
+            gdir.add_to_diagnostics(perc_string, len(pnan[0]) / (nx * ny))
+        else:
+            diagnostics_dict[f'dem_needed_{log_string}polation'] = True
+            diagnostics_dict[perc_string] = len(pnan[0]) / (nx * ny)
 
     # Correct the DEM
     valid_mask = np.isfinite(dem)
@@ -759,8 +786,7 @@ def process_dem(gdir):
     if np.any(~valid_mask):
         # We interpolate
         if np.sum(~valid_mask) > (0.25 * nx * ny):
-            log.info('({}) more than 25% NaNs in DEM'.format(gdir.rgi_id))
-        xx, yy = gdir.grid.ij_coordinates
+            log.info('({}) more than 25% NaNs in DEM'.format(grid_name))
         pnan = np.nonzero(~valid_mask)
         pok = np.nonzero(valid_mask)
         points = np.array((np.ravel(yy[pok]), np.ravel(xx[pok]))).T
@@ -770,15 +796,13 @@ def process_dem(gdir):
                                  method='linear')
         except ValueError:
             raise InvalidDEMError('DEM interpolation not possible.')
-        log.info(gdir.rgi_id + ': DEM needed interpolation.')
-        gdir.add_to_diagnostics('dem_needed_interpolation', True)
-        gdir.add_to_diagnostics('dem_invalid_perc', len(pnan[0]) / (nx * ny))
+
+        _log_and_diagnostics(pnan, log_string='inter')
 
     isfinite = np.isfinite(dem)
     if np.any(~isfinite):
         # interpolation will still leave NaNs in DEM:
         # extrapolate with NN if needed (e.g. coastal areas)
-        xx, yy = gdir.grid.ij_coordinates
         pnan = np.nonzero(~isfinite)
         pok = np.nonzero(isfinite)
         points = np.array((np.ravel(yy[pok]), np.ravel(xx[pok]))).T
@@ -788,13 +812,11 @@ def process_dem(gdir):
                                  method='nearest')
         except ValueError:
             raise InvalidDEMError('DEM extrapolation not possible.')
-        log.info(gdir.rgi_id + ': DEM needed extrapolation.')
-        gdir.add_to_diagnostics('dem_needed_extrapolation', True)
-        gdir.add_to_diagnostics('dem_extrapol_perc', len(pnan[0]) / (nx * ny))
+        _log_and_diagnostics(pnan, log_string='extra')
 
     if np.min(dem) == np.max(dem):
         raise InvalidDEMError('({}) min equal max in the DEM.'
-                              .format(gdir.rgi_id))
+                              .format(grid_name))
 
     # Clip topography to 0 m a.s.l.
     if cfg.PARAMS['clip_dem_to_zero']:
@@ -802,7 +824,7 @@ def process_dem(gdir):
 
     # Smooth DEM?
     if cfg.PARAMS['smooth_window'] > 0.:
-        gsize = np.rint(cfg.PARAMS['smooth_window'] / gdir.grid.dx)
+        gsize = np.rint(cfg.PARAMS['smooth_window'] / dx)
         smoothed_dem = gaussian_blur(dem, int(gsize))
     else:
         smoothed_dem = dem.copy()
@@ -811,8 +833,12 @@ def process_dem(gdir):
     if cfg.PARAMS['clip_dem_to_zero']:
         utils.clip_min(smoothed_dem, 0, out=smoothed_dem)
 
+    if gdir is None:
+        with open(os.path.join(fpath, 'dem_diagnostics.json'), 'w') as f:
+            json.dump(diagnostics_dict, f)
+
     # Write to file
-    with GriddedNcdfFile(gdir=gdir, reset=True) as nc:
+    with GriddedNcdfFile(gdir=gdir, grid=dem_grid, fpath=fpath, reset=True) as nc:
 
         v = nc.createVariable('topo', 'f4', ('y', 'x',), zlib=True)
         v.units = 'm'

From 1a1be8e98b0d3c6940a398ec8ecc286be14931af Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Tue, 5 Dec 2023 19:39:10 +0100
Subject: [PATCH 13/16] now using the correct data (gridded_simulation)

---
 oggm/sandbox/distribute_2d.py | 58 ++++++++++++++++++++++++++++++-----
 1 file changed, 51 insertions(+), 7 deletions(-)

diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py
index a46c3816b..71811d3c4 100644
--- a/oggm/sandbox/distribute_2d.py
+++ b/oggm/sandbox/distribute_2d.py
@@ -453,17 +453,61 @@ def combine_distributed_thickness(gdirs, output_folder, suffix='', source=None):
 
     # add 'write_dem_to_gridded_data' function and add dem to the gridded data
 
-
-    # add individual add individual distributed thicknesses to gridded data
-    total_data = np.zeros((combined_grid.ny, combined_grid.nx))
+    def _get_gridded_simulation_path(gdir, suffix=''):
+        """
+        can insert a suffix in the gridded_simulation path
+        Parameters
+        ----------
+        gdir - :py:class:`oggm.GlacierDirectory`
+        An OGGM glacier directory
+        suffix - str, optional
+        suffix chosen for the simulation run
+
+        Returns
+        -------
+
+        """
+        gridded_sim_path = gdir.get_filepath('gridded_simulation')
+        return gridded_sim_path[:-3] + suffix + gridded_sim_path[-3:]
+
+    # combine individual distributed thicknesses
+    combined_dt = None
+    time_coord = None
     for gdir in gdirs:
-        with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
+        with xr.open_dataset(_get_gridded_simulation_path(gdir, suffix)) as ds:
             ds = ds.load()
 
+        # assuming that the same simulation ran the same time for each glacier
+        # => taking the simulation duration from the first gdir
+        if combined_dt is None:
+            time_coord = ds.time
+            combined_dt = np.zeros((len(time_coord), combined_grid.ny, combined_grid.nx))
+
         r_data = combined_grid.map_gridded_data(ds['distributed_thickness'],
                                                 grid=gdir.grid)
-        total_data += r_data.filled(0)
-    total_data = np.where(total_data == 0, np.NaN, total_data)
-    return total_data
+        combined_dt += r_data.filled(0)
+    combined_dt = np.where(combined_dt == 0, np.NaN, combined_dt)
+
+    # add combined data to the gridded_data
+    # with xr.open_dataset(os.path.join(output_folder, 'gridded_data.nc')) as ds:
+    #     ds = ds.load()
+    # ds = ds.assign_coords({'time': time_coord.values})
+    # ds = ds.assign({'combined_distributed_thickness': (('time', 'y', 'x'), combined_dt)})
+    # return ds
+
+    with GriddedNcdfFile(gdir=None, fpath=output_folder, grid=combined_grid) as gridded_data:
+        # add time dimension to dataframe
+        gridded_data.createDimension('time', len(time_coord))
+        v = gridded_data.createVariable('time', 'f4', ('time',), zlib=True)
+        v.units = 'years'
+        v.long_name = 'time coordinate of the simulation'
+        v.standard_name = 'time'
+        v[:] = time_coord.values
+        # add combined distributed thickness
+        v = gridded_data.createVariable('combined_distributed_thickness', 'f4', ('time', 'y', 'x',), zlib=True)
+        v.units = 'm'
+        v.long_name = 'combined distributed thickness'
+        v[:] = combined_dt
+
 
 

From 09cb8311fd42bc8ec299a4222bb9c419c3a07951 Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Tue, 5 Dec 2023 19:40:57 +0100
Subject: [PATCH 14/16] creating gridded_data.nc for combined grid with file
 ending

---
 oggm/core/gis.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/oggm/core/gis.py b/oggm/core/gis.py
index c2d29aced..b57fac355 100644
--- a/oggm/core/gis.py
+++ b/oggm/core/gis.py
@@ -689,7 +689,7 @@ def __init__(self, gdir=None, grid=None, fpath=None,
                 raise InvalidParamsError('If you do not provide a gdir you must'
                                          'define grid and fpath! Given grid='
                                          f'{grid} and fpath={fpath}.')
-            self.fpath = os.path.join(fpath, basename)
+            self.fpath = os.path.join(fpath, basename + '.nc')
             self.grid = grid
 
         if reset and os.path.exists(self.fpath):

From b63a72212c57a384ea400fc54722aa2bd7069183 Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Wed, 6 Dec 2023 11:44:04 +0100
Subject: [PATCH 15/16] merged RGI7 related changes with 'multiple glacier
 visualization' code changes

---
 oggm/core/gis.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/oggm/core/gis.py b/oggm/core/gis.py
index c24d09333..cc91d2776 100644
--- a/oggm/core/gis.py
+++ b/oggm/core/gis.py
@@ -443,7 +443,7 @@ def _get_nodata(rio_ds):
         dem_ds.close()
 
 
-def get_dem_for_grid(grid, fpath, source=None, rgi_id=None):
+def get_dem_for_grid(grid, fpath, source=None, gdir=None):
     """
     Fetch a DEM from source, reproject it to the extent defined by grid and
     saves it to disk.
@@ -467,6 +467,11 @@ def get_dem_for_grid(grid, fpath, source=None, rgi_id=None):
     """
     minlon, maxlon, minlat, maxlat = grid.extent_in_crs(crs=salem.wgs84)
     extent_ll = [[minlon, maxlon], [minlat, maxlat]]
+    if gdir is not None:
+        rgi_id = gdir.rgi_id
+    else:
+        rgi_id = None
+
     grid_prop = {
         'utm_proj': grid.proj,
         'dx': grid.dx,
@@ -479,7 +484,7 @@ def get_dem_for_grid(grid, fpath, source=None, rgi_id=None):
     source = check_dem_source(source, extent_ll, rgi_id=rgi_id)
 
     dem_list, dem_source = get_topo_file((minlon, maxlon), (minlat, maxlat),
-                                         rgi_id=rgi_id,
+                                         gdir=gdir,
                                          dx_meter=grid_prop['dx'],
                                          source=source)
 
@@ -548,7 +553,7 @@ def define_glacier_region(gdir, entity=None, source=None):
 
     dem_list, dem_source = get_dem_for_grid(grid=tmp_grid,
                                             fpath=gdir.get_filepath('dem'),
-                                            source=source, rgi_id=gdir.rgi_id)
+                                            source=source, gdir=gdir)
 
     # Glacier grid
     x0y0 = (ulx+dx/2, uly-dx/2)  # To pixel center coordinates

From eda6a10edd8a2bdd1ab611540d184dc1732a5d54 Mon Sep 17 00:00:00 2001
From: afisc <mail@alexfischer.eu>
Date: Thu, 21 Dec 2023 22:42:39 +0100
Subject: [PATCH 16/16] adapted combine_distributed_thickness

---
 oggm/sandbox/distribute_2d.py | 104 +++++++++++++++++++---------------
 1 file changed, 57 insertions(+), 47 deletions(-)

diff --git a/oggm/sandbox/distribute_2d.py b/oggm/sandbox/distribute_2d.py
index 71811d3c4..3e47ef9ea 100644
--- a/oggm/sandbox/distribute_2d.py
+++ b/oggm/sandbox/distribute_2d.py
@@ -451,52 +451,48 @@ def combine_distributed_thickness(gdirs, output_folder, suffix='', source=None):
     # processes dem(smoothing etc.) and creates gridded_data and adds DEM data
     process_dem.unwrapped(gdir=None, grid=combined_grid, fpath=output_folder)
 
-    # add 'write_dem_to_gridded_data' function and add dem to the gridded data
-
-    def _get_gridded_simulation_path(gdir, suffix=''):
-        """
-        can insert a suffix in the gridded_simulation path
-        Parameters
-        ----------
-        gdir - :py:class:`oggm.GlacierDirectory`
-        An OGGM glacier directory
-        suffix - str, optional
-        suffix chosen for the simulation run
-
-        Returns
-        -------
-
-        """
-        gridded_sim_path = gdir.get_filepath('gridded_simulation')
-        return gridded_sim_path[:-3] + suffix + gridded_sim_path[-3:]
-
-    # combine individual distributed thicknesses
-    combined_dt = None
-    time_coord = None
-    for gdir in gdirs:
-        with xr.open_dataset(_get_gridded_simulation_path(gdir, suffix)) as ds:
-            ds = ds.load()
-
-        # assuming that the same simulation ran the same time for each glacier
-        # => taking the simulation duration from the first gdir
-        if combined_dt is None:
-            time_coord = ds.time
-            combined_dt = np.zeros((len(time_coord), combined_grid.ny, combined_grid.nx))
-
-        r_data = combined_grid.map_gridded_data(ds['distributed_thickness'],
-                                                grid=gdir.grid)
-        combined_dt += r_data.filled(0)
-    combined_dt = np.where(combined_dt == 0, np.NaN, combined_dt)
-
-    # add combined data to the gridded_data
-    # with xr.open_dataset(os.path.join(output_folder, 'gridded_data.nc')) as ds:
-    #     ds = ds.load()
-    # ds = ds.assign_coords({'time': time_coord.values})
-    # ds = ds.assign({'combined_distributed_thickness': (('time', 'y', 'x'), combined_dt)})
-    # return ds
-
-    with GriddedNcdfFile(gdir=None, fpath=output_folder, grid=combined_grid) as gridded_data:
-        # add time dimension to dataframe
+
+    def _combine_variable_on_common_grid(var_name, source_file='gridded_data', time_dim=False, suffix=''):
+        combined_data = None
+        time_coord = None
+        for gdir in gdirs:
+            with xr.open_dataset(gdir.get_filepath(source_file, filesuffix=suffix)) as ds:
+                ds = ds.load()
+
+            # assuming that the same simulation ran the same time for each glacier
+            # => taking the simulation duration from the first gdir
+            if combined_data is None:
+                if time_dim:
+                    if time_coord is None:
+                        time_coord = ds.time
+                    else:
+                        if time_coord == ds.time:
+                            continue
+                        else:
+                            raise ValueError('The simulation times of the simulations that you try to combine are '
+                                             'different from each other.\n'
+                                             f'The Glacier with the rgi_id {gdir.rgi_id}\ has {len(ds.time)} timesteps,'
+                                             f'not matching with {gdirs[0].rgi_id}\'s {len(time_coord)} timesteps.')
+                    combined_data = np.zeros((len(time_coord), combined_grid.ny, combined_grid.nx))
+                else:
+                    combined_data = np.zeros((combined_grid.ny, combined_grid.nx))
+            r_data = combined_grid.map_gridded_data(ds[var_name],
+                                                    grid=gdir.grid)
+            combined_data += r_data.filled(0)
+        combined_data = np.where(combined_data == 0, np.NaN, combined_data)
+
+        return combined_data, time_coord
+
+    combined_distributed_thickness, time_coord = _combine_variable_on_common_grid('distributed_thickness',
+                                                                      source_file='gridded_simulation', time_dim=True,
+                                                                      suffix=suffix)
+    combined_glacier_mask, _ = _combine_variable_on_common_grid('glacier_mask')
+    # combined_topo_smoothed = _combine_variable_on_common_grid('topo_smoothed')
+    combined_bedrock, _ = _combine_variable_on_common_grid('bedrock', source_file='gridded_simulation', suffix=suffix)
+    combined_glacier_ext, _ = _combine_variable_on_common_grid('glacier_ext')
+
+
+    with utils.ncDataset(os.path.join(output_folder, 'gridded_data.nc'), 'a') as gridded_data:
         gridded_data.createDimension('time', len(time_coord))
         v = gridded_data.createVariable('time', 'f4', ('time',), zlib=True)
         v.units = 'years'
@@ -504,10 +500,24 @@ def _get_gridded_simulation_path(gdir, suffix=''):
         v.standard_name = 'time'
         v[:] = time_coord.values
         # add combined distributed thickness
+
         v = gridded_data.createVariable('combined_distributed_thickness', 'f4', ('time', 'y', 'x',), zlib=True)
         v.units = 'm'
         v.long_name = 'combined distributed thickness'
-        v[:] = combined_dt
+        v[:] = combined_distributed_thickness
+
+        v = gridded_data.createVariable('combined_glacier_mask', 'f4', ('y', 'x',), zlib=True)
+        v.units = 'm'
+        v.long_name = 'combined glacier mask'
+        v[:] = combined_glacier_mask
 
+        v = gridded_data.createVariable('combined_bedrock', 'f4', ('y', 'x',), zlib=True)
+        v.units = 'm'
+        v.long_name = 'combined bedrock'
+        v[:] = combined_bedrock
 
+        v = gridded_data.createVariable('combined_glacier_ext', 'f4', ('y', 'x',), zlib=True)
+        v.units = 'm'
+        v.long_name = 'combined glacier extent'
+        v[:] = combined_glacier_ext