code: def meanstack(infn: Path, Navg: int, ut1: Optional[datetime]=None,
                method: str='mean') -> Tuple[np.ndarray, Optional[datetime]]:
  
      infn = Path(infn).expanduser()
  # %% parse indicies to load
      if isinstance(Navg, slice):
          key = Navg
      elif isinstance(Navg, int):
          key = slice(0, Navg)
      elif len(Navg) == 1:
          key = slice(0, Navg[0])
      elif len(Navg) == 2:
          key = slice(Navg[0], Navg[1])
      else:
          raise ValueError(f'not sure what you mean by Navg={Navg}')
  # %% load images
      
      if infn.suffix == '.h5':
          img, ut1 = _h5mean(infn, ut1, key, method)
      elif infn.suffix == '.fits':
          with fits.open(infn, mode='readonly', memmap=False) as f:  # mmap doesn't work with BZERO/BSCALE/BLANK
              img = collapsestack(f[0].data, key, method)
      elif infn.suffix == '.mat':
          img = loadmat(infn)
          img = collapsestack(img['data'].T, key, method)  # matlab is fortran order
      else:  # .tif etc.
          img = imageio.imread(infn, as_gray=True)
          if img.ndim in (3, 4) and img.shape[-1] == 3:  # assume RGB
              img = collapsestack(img, key, method)
  
      return img, ut1
docs: some methods handled individually to improve efficiency with huge files