Feature request: Compute cross-correlation (similar to pd.Series.corr()) of gridded data

[hrishikeshac]

As a earth scientist regularly dealing with 3D data (time, latitude, longitude), I believe it would be great to be able to perform cross-correlation on DataArrays by specifying the axis. It's usage could look like: a.corr(b, axis = 0). It would be even more useful if the two arrays need not have the same dimensions (e.g. 'b' could be a time series).

Currently, the only way to compute this that I am aware of, is by looping through each grid, converting the time series to pd.Series(), and then computing the correlation. This takes a long time.
Would also appreciate suggestions to a faster algorithm.

[shoyer]

The first step here is to find a library that implements the desired functionality on pure NumPy arrays, ideally in a vectorized fashion. Then it should be pretty straightforward to wrap in xarray.

[rabernat]

I agree this would be very useful. But it is also feature creep. There is an extremely wide range of such functions that could hypothetically be put into the xarray package. (all of scipy.signal for example) At some point the community should decide what is the intended scope of xarray itself vs. packages built on top of xarray.

[serazing]

I agree with @rabernat in the sense that it could be part of another package (e.g., signal processing). This would also allow the computation of statistical test to assess the significance of the correlation (which is useful since correlation may often be misinterpreted without statistical tests).

[shoyer]

That said, correlation coefficients are a pretty fundamental operation for working with data. I could see implementing a basic corr in xarray and referring to a separate signal processing package for more options in the docstring.

[rabernat]

To be clear, I am not say that this does not belong in xarray.

I'm saying that we lack clear general guidelines for how to determine whether a particular function belongs in xarray. The criterion of a "pretty fundamental operation for working with data" is a good starting point. I would add:

• used across a wide range of scientific disciplines
• clear, unambiguous / uncontroversial definition
• numpy implementation already exists

corr meets all of these criteria. Many others (e.g. interpolation, convolution, curve fitting) do as well. Expanding xarray beyond the numpy ufuncs opens the door to supporting these things. I'm just saying it should be conscious, deliberate decision, given the limits on developer time.

Many of these things will be pretty trivial once .apply() is here. So perhaps it's not a big deal.

[fmaussion]

I'll chime in here to ask a usage question: what is the recommended way to compute correlation maps with xarray? I.e. I have a dataarray of dims (time, lat, lon) and I'd like to correlate every single grid point with a timeseries of dim (time) to get a correlation map of dim (lat, lon). My current strategy is a wonderfully unpythonic double loop over lons and lats, and I wonder if there's better way?

[hrishikeshac]

FYI @shoyer @fmaussion , I had to revisit the problem and ended up writing a function to compute vectorized cross-correlation, covariance, regression calculations (along with p-value and standard error) for xr.DataArrays. Essentially, I tried to mimic scipy.stats.linregress() but for multi-dimensional data, and included the ability to compute lagged relationships. Here's the function and its demonstration; please feel free to incorporate it in xarray if deemed useful: https://hrishichandanpurkar.blogspot.com/2017/09/vectorized-functions-for-correlation.html

[sebhahn]

@hrishikeshac I was just looking for a function doing a regression between two datasets (x, y, time), so thanks for your function! However, I'm still wondering whether there is a much faster C (or Cython) implementation doing these kind of things?

[max-sixty]

I'm up for adding .corr to xarray

What do want this to look like? It's a bit different from most xarray functions, which either return the same shape or reduce one dimension.

• The basic case here would take a n x m array and return an m x m correlation matrix. We could easily wrap https://docs.scipy.org/doc/numpy/reference/generated/numpy.corrcoef.html
• Another case would be take two similarly sized arrays (with the option of broadcasting) and return an array with one dimension reduced. For example 200 x 10 and 200, return a 10 array.
• I need to think about how those extrapolate to multiple dimensions

Should I start with the first case and then we can expand as needed?

[shoyer]

I tend to view the second case as a generalization of the first case. I would also hesitate to implement the n x m array -> m x m correlation matrix version because xarray doesn't handle repeated dimensions well.

I think the basic implementation of this looks quite similar to what I wrote here for calculating the Pearson correlation as a NumPy gufunc:
http://xarray.pydata.org/en/stable/dask.html#automatic-parallelization

The main difference is that we might naturally want to support summing over multiple dimensions at once via the dim argument, e.g., something like:

# untested!
def covariance(x, y, dim=None):
    return xarray.dot(x - x.mean(dim), y - y.mean(dim), dim=dim)

def corrrelation(x, y, dim=None):
    # dim should default to the intersection of x.dims and y.dims
    return covariance(x, y, dim) / (x.std(dim) * y.std(dim))

If you want to achieve the equivalent of np.corr on an array with dimensions ('n', 'm') with this, you just write something like correlation(x, x.rename({'m': 'm2'}), dim='n').

[hrishikeshac]

Sometime back I wrote a package based on xarray regarding this. I would be happy to be involved in implementing it in xarray as well, but I am new to contributing to such a large-scale project and it looks a bit intimidating!

[max-sixty]

@hrishikeshac if you'd like to contribute, we can help you along - xarray is a v welcoming project!

And from mvstats it looks like you're already up to speed

Let us know