code: def plot_histograms_extra_dims(dataobj, keys, **kwargs):
      if isinstance(dataobj, pd.DataFrame):
          df_raw = dataobj
      else:
          df_raw = dataobj.data
      if kwargs.get('subquery', False):
          df = df_raw.query(kwargs.get('subquery'))
      else:
          df = df_raw
      split_timestamps = True
      if split_timestamps:
          group_timestamps = df.groupby('timestep')
          N_ts = len(group_timestamps.groups.keys())
      else:
          group_timestamps = ('all', df)
          N_ts = 1
      columns = keys
      N_c = len(columns)
      plot_log10 = kwargs.get('log10plot', False)
      if plot_log10:
          transformers = ['lin', 'log10']
          N_log10 = 2
      else:
          transformers = ['lin', ]
          N_log10 = 1
      Nx_max = kwargs.get('Nx', 4)
      N = N_ts * N_c * N_log10
      Nx = min(Nx_max, N)
      Ny = int(np.ceil(N / Nx))
      size_x = 5 * Nx / 2.54
      size_y = 5 * Ny / 2.54
      fig, axes = plt.subplots(Ny, Nx, figsize=(size_x, size_y), sharex=True,
                               sharey=True)
      axes = np.atleast_2d(axes)
      index = 0
      for ts_name, tgroup in group_timestamps:
          for column in columns:
              for transformer in transformers:
                  subdata_raw = tgroup[column].values
                  subdata = subdata_raw[~np.isnan(subdata_raw)]
                  subdata = subdata[np.isfinite(subdata)]
                  if transformer == 'log10':
                      subdata_log10_with_nan = np.log10(subdata[subdata > 0])
                      subdata_log10 = subdata_log10_with_nan[~np.isnan(
                          subdata_log10_with_nan)
                      ]
                      subdata_log10 = subdata_log10[np.isfinite(subdata_log10)]
                      subdata = subdata_log10
                  ax = axes.flat[index]
                  ax.hist(
                      subdata,
                      _get_nr_bins(subdata.size),
                  )
                  ax.set_xlabel(
                      units.get_label(column)
                  )
                  ax.set_ylabel('count')
                  ax.xaxis.set_major_locator(mpl.ticker.MaxNLocator(3))
                  ax.tick_params(axis='both', which='major', labelsize=6)
                  ax.tick_params(axis='both', which='minor', labelsize=6)
                  ax.set_title("timestep: %d" % ts_name)
                  index += 1
      for ax in axes[:, 1:].flat:
          ax.set_ylabel('')
      for ax in axes[:-1, :].flat:
          ax.set_xlabel('')
      fig.tight_layout()
      return fig