TomNicholas
commented
3 years ago Update on this: in a PR to xhistogram we have a rough proof-of-principle for a dask-parallelized, axis-aware implementation of N-dimensional histogram calculations, suitable for eventually integrating into xarray.
We still need to complete the work over in xhistogram, but for now I want to suggest what I think the eventual API should be for this functionality within xarray:
Top-level function
xhistogram's xarray API is essentially one histogram function, which accepts one or more xarray DataArrays. Therefore we it makes sense to add histogram or hist as a top-level function, similar to the existing cov, corr, dot or polyval.
New methods
We could also add a datarray.hist() method for 1-D histograms and possibly a dataset.hist(vars=['density', 'temperature']) for quickly making N-D histograms.
The existing plot.hist method (the tricky bit)
There is already a da.plot.hist() method, which is a paper-thin wrapper around matplotlib.pyplot.hist, which flattens the dataarray before plotting the result. It would be nice if this internally dispatched to the new da.hist() method before plotting the result, but pyplot.hist does both the bincounting and the plotting, so it might not be simple to do that.
This is also potentially related to @dcherian 's PR for facets and hue with hist, in that a totally consistent treatment would use the axis-aware histogram algorithm to calculate each separate facet histogram by looping over non-binned dimensions. Again the problem is that AFAIK matplotlib doesn't offer a quick way to plot a histogram without recomputing the bins and counts. Any suggestions here?
Adding an optional dim argument to da.plot.hist() doesn't really make sense unless we also add faceting, because otherwise the only shape of result that da.plot.hist() could actually plot is one where we have binned over the entire flattened array.
It would also be nice if
da.hist().plot.hist()
was identical to
da.plot.hist()
which requires the format of the output of da.hist() to be compatible with da.plot.hist().
(We shouldn't need any kind of new da.plot.hist2d() method because as the xhistogram docs show you can already make very nice 2D histogram plots with da.plot().)
Signature
xhistogram adds bin coordinates (the bin centers) to the output dataarray, named after the quantities that were binned.
Following what we currently have, a minimal top-level function signature looks like
def hist(*datarrays, bins=None, dim=None, weights=None, density=False):
"""
Histogram applied along specified dimensions.
If any of the supplied arguments are dask arrays it will use `dask.array.blockwise`
internally to parallelize over all chunks.
datarrays : xarray.DataArray objects
Input data. The number of input arguments determines the dimensionality of
the histogram. For example, two arguments prodoce a 2D histogram.
bins : int or array_like or a list of ints or arrays, or list of DataArrays, optional
If a list, there should be one entry for each item in ``args``.
The bin specification:
* If int, the number of bins for all arguments in ``args``.
* If array_like, the bin edges for all arguments in ``args``.
* If a list of ints, the number of bins for every argument in ``args``.
* If a list arrays, the bin edges for each argument in ``args``
(required format for Dask inputs).
* A combination [int, array] or [array, int], where int
is the number of bins and array is the bin edges.
* If a list of DataArrays, the bins for each argument in ``args``
The DataArrays can be multidimensional, but must not have any
dimensions shared with the `dim` argument.
When bin edges are specified, all but the last (righthand-most) bin include
the left edge and exclude the right edge. The last bin includes both edges.
A ``TypeError`` will be raised if ``args`` contains dask arrays and
``bins`` are not specified explicitly as a list of arrays.
dim : tuple of strings, optional
Dimensions over which which the histogram is computed. The default is to
compute the histogram of the flattened array.
weights : array_like, optional
An array of weights, of the same shape as `a`. Each value in
`a` only contributes its associated weight towards the bin count
(instead of 1). If `density` is True, the weights are
normalized, so that the integral of the density over the range
remains 1. NaNs in the weights input will fill the entire bin with
NaNs. If there are NaNs in the weights input call ``.fillna(0.)``
before running ``histogram()``.
density : bool, optional
If ``False``, the result will contain the number of samples in
each bin. If ``True``, the result is the value of the
probability *density* function at the bin, normalized such that
the *integral* over the range is 1. Note that the sum of the
histogram values will not be equal to 1 unless bins of unity
width are chosen; it is not a probability *mass* function.
"""Weights could also possibly be set via the .weighted() method that we already have for other operations.
Checklist
Desired features in order to fully deprecate xhistogram:
- [ ] axis-aware (ability to loop over dimensions instead of binning over them)
- [ ] optional dask parallelization across all dimensions
- [ ] weights
- [ ]
accept dask-aware bin arrays? - [ ] accept multi-dimensional bins arguments? (see https://github.com/xgcm/xhistogram/issues/28)
- any others?
cc @dougiesquire @gjoseph92
Illviljan
aaronspring
dougiesquire
dcherian
On today's dev call, we discussed the possible role that numpy_groupies could play in xarray (#4540). I noted that many of the use cases for advanced grouping overlap significantly with histogram-type operations. A major need that we have is to take [weighted] histograms over some, but not all, axes of DataArrays. Since groupby doesn't allow this (see #1013), we started the standalone xhistogram package.
Given the broad usefulness of this feature, I suggested that we might want to deprecate xhistogram and move the histogram function to xarray. We may want to also reimplement it using numpy_groupies, which I think is smarter than our implementation in xhistogram.
I've opened this issue to keep track of the idea.