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