docs: def stack(srcfiles, dstfile, resampling, targetres, dstnodata, srcnodata=None, shapefile=None, layernames=None,
            sortfun=None, separate=False, overwrite=False, compress=True, cores=4):
      """
      function for mosaicking, resampling and stacking of multiple raster files into a 3D data cube
  
      Parameters
      ----------
      srcfiles: list
          a list of file names or a list of lists; each sub-list is treated as a task to mosaic its containing files
      dstfile: str
          the destination file or a directory (if `separate` is True)
      resampling: {near, bilinear, cubic, cubicspline, lanczos, average, mode, max, min, med, Q1, Q3}
          the resampling method; see `documentation of gdalwarp <https://www.gdal.org/gdalwarp.html>`_.
      targetres: tuple or list
          two entries for x and y spatial resolution in units of the source CRS
      srcnodata: int, float or None
          the nodata value of the source files; if left at the default (None), the nodata values are read from the files
      dstnodata: int or float
          the nodata value of the destination file(s)
      shapefile: str, Vector or None
          a shapefile for defining the spatial extent of the destination files
      layernames: list
          the names of the output layers; if `None`, the basenames of the input files are used; overrides sortfun
      sortfun: function
          a function for sorting the input files; not used if layernames is not None.
          This is first used for sorting the items in each sub-list of srcfiles;
          the basename of the first item in a sub-list will then be used as the name for the mosaic of this group.
          After mosaicing, the function is again used for sorting the names in the final output
          (only relevant if `separate` is False)
      separate: bool
          should the files be written to a single raster stack (ENVI format) or separate files (GTiff format)?
      overwrite: bool
          overwrite the file if it already exists?
      compress: bool
          compress the geotiff files?
      cores: int
          the number of CPU threads to use; this is only relevant if `separate` is True, in which case each
          mosaicing/resampling job is passed to a different CPU
  
      Returns
      -------
  
      Notes
      -----
      This function does not reproject any raster files. Thus, the CRS must be the same for all input raster files.
      This is checked prior to executing gdalwarp. In case a shapefile is defined, it is internally reprojected to the
      raster CRS prior to retrieving its extent.
      
      Examples
      --------
      
      .. code-block:: python
      
          from pyroSAR.ancillary import groupbyTime, find_datasets, seconds
          from spatialist.raster import stack
          
          # find pyroSAR files by metadata attributes
          archive_s1 = '/.../sentinel1/GRD/processed'
          scenes_s1 = find_datasets(archive_s1, sensor=('S1A', 'S1B'), acquisition_mode='IW')
          
          # group images by acquisition time
          groups = groupbyTime(images=scenes_s1, function=seconds, time=30)
          
          # mosaic individual groups and stack the mosaics to a single ENVI file
          # only files overlapping with the shapefile are selected and resampled to its extent
          stack(srcfiles=groups, dstfile='stack', resampling='bilinear', targetres=(20, 20),
                srcnodata=-99, dstnodata=-99, shapefile='site.shp', separate=False)
      """
      # perform some checks on the input data
      
      if len(dissolve(srcfiles)) == 0:
          raise RuntimeError('no input files provided to function raster.stack')
      
      if layernames is not None:
          if len(layernames) != len(srcfiles):
              raise RuntimeError('mismatch between number of source file groups and layernames')
      
      if not isinstance(targetres, (list, tuple)) or len(targetres) != 2:
          raise RuntimeError('targetres must be a list or tuple with two entries for x and y resolution')
      
      if len(srcfiles) == 1 and not isinstance(srcfiles[0], list):
          raise RuntimeError('only one file specified; nothing to be done')
      
      if resampling not in ['near', 'bilinear', 'cubic', 'cubicspline', 'lanczos',
                            'average', 'mode', 'max', 'min', 'med', 'Q1', 'Q3']:
          raise RuntimeError('resampling method not supported')
      
      projections = list()
      for x in dissolve(srcfiles):
          try:
              projection = Raster(x).projection
          except RuntimeError as e:
              print('cannot read file: {}'.format(x))
              raise e
          projections.append(projection)
      
      projections = list(set(projections))
      if len(projections) > 1:
          raise RuntimeError('raster projection mismatch')
      elif projections[0] == '':
          raise RuntimeError('could not retrieve the projection from any of the {} input images'.format(len(srcfiles)))
      else:
          srs = projections[0]
      ##########################################################################################
      # read shapefile bounding coordinates and reduce list of rasters to those overlapping with the shapefile
      
      if shapefile is not None:
          shp = shapefile.clone() if isinstance(shapefile, Vector) else Vector(shapefile)
          shp.reproject(srs)
          ext = shp.extent
          arg_ext = (ext['xmin'], ext['ymin'], ext['xmax'], ext['ymax'])
          for i, item in enumerate(srcfiles):
              group = item if isinstance(item, list) else [item]
              if layernames is None and sortfun is not None:
                  group = sorted(group, key=sortfun)
              group = [x for x in group if intersect(shp, Raster(x).bbox())]
              if len(group) > 1:
                  srcfiles[i] = group
              elif len(group) == 1:
                  srcfiles[i] = group[0]
              else:
                  srcfiles[i] = None
          shp.close()
          srcfiles = list(filter(None, srcfiles))
      else:
          arg_ext = None
      ##########################################################################################
      # set general options and parametrization
      
      dst_base = os.path.splitext(dstfile)[0]
      
      options_warp = {'options': ['-q'],
                      'format': 'GTiff' if separate else 'ENVI',
                      'outputBounds': arg_ext, 'multithread': True,
                      'dstNodata': dstnodata,
                      'xRes': targetres[0], 'yRes': targetres[1],
                      'resampleAlg': resampling}
      
      if overwrite:
          options_warp['options'] += ['-overwrite']
      
      if separate and compress:
          options_warp['options'] += ['-co', 'COMPRESS=DEFLATE', '-co', 'PREDICTOR=2']
      
      options_buildvrt = {'outputBounds': arg_ext}
      
      if srcnodata is not None:
          options_warp['srcNodata'] = srcnodata
          options_buildvrt['srcNodata'] = srcnodata
      ##########################################################################################
      # create VRT files for mosaicing
      
      for i, group in enumerate(srcfiles):
          if isinstance(group, list):
              if len(group) > 1:
                  base = group[0]
                  # in-memory VRT files cannot be shared between multiple processes on Windows
                  # this has to do with different process forking behaviour
                  # see function spatialist.ancillary.multicore and this link:
                  # https://stackoverflow.com/questions/38236211/why-multiprocessing-process-behave-differently-on-windows-and-linux-for-global-o
                  vrt_base = os.path.splitext(os.path.basename(base))[0] + '.vrt'
                  if platform.system() == 'Windows':
                      vrt = os.path.join(tempfile.gettempdir(), vrt_base)
                  else:
                      vrt = '/vsimem/' + vrt_base
                  gdalbuildvrt(group, vrt, options_buildvrt)
                  srcfiles[i] = vrt
              else:
                  srcfiles[i] = group[0]
          else:
              srcfiles[i] = group
      ##########################################################################################
      # define the output band names
      
      # if no specific layernames are defined, sort files by custom function
      if layernames is None and sortfun is not None:
          srcfiles = sorted(srcfiles, key=sortfun)
      
      # use the file basenames without extension as band names if none are defined
      bandnames = [os.path.splitext(os.path.basename(x))[0] for x in srcfiles] if layernames is None else layernames
      
      if len(list(set(bandnames))) != len(bandnames):
          raise RuntimeError('output bandnames are not unique')
      ##########################################################################################
      # create the actual image files
      
      if separate:
          if not os.path.isdir(dstfile):
              os.makedirs(dstfile)
          dstfiles = [os.path.join(dstfile, x) + '.tif' for x in bandnames]
          jobs = [x for x in zip(srcfiles, dstfiles)]
          if not overwrite:
              jobs = [x for x in jobs if not os.path.isfile(x[1])]
              if len(jobs) == 0:
                  print('all target tiff files already exist, nothing to be done')
                  return
          srcfiles, dstfiles = map(list, zip(*jobs))
          
          multicore(gdalwarp, cores=cores, multiargs={'src': srcfiles, 'dst': dstfiles}, options=options_warp)
      else:
          if len(srcfiles) == 1:
              options_warp['format'] = 'GTiff'
              if not dstfile.endswith('.tif'):
                  dstfile = os.path.splitext(dstfile)[0] + '.tif'
              gdalwarp(srcfiles[0], dstfile, options_warp)
          else:
              # create VRT for stacking
              vrt = '/vsimem/' + os.path.basename(dst_base) + '.vrt'
              options_buildvrt['options'] = ['-separate']
              gdalbuildvrt(srcfiles, vrt, options_buildvrt)
              
              # warp files
              gdalwarp(vrt, dstfile, options_warp)
              
              # edit ENVI HDR files to contain specific layer names
              with envi.HDRobject(dstfile + '.hdr') as hdr:
                  hdr.band_names = bandnames
                  hdr.write()