issue_comments

I opened https://github.com/pydata/xarray/issues/1850 to discuss xarray-contrib.

I've also contributed to developing a python package (xrft) for fft keeping the awareness of the metadata of multidimensional xarray datasets.

I would also be very interested in seeing your codes @lamorton. Overall, I think the xarray community could really benefit from some kind of centralized contrib package which has a low barrier to entry for these kinds of functions. So far, I suspect there has been a large amount of code duplication for routine tasks like the fft, since I have also written a function for that.

@gajomi I can find a place to upload what I have. I foresee some difficulty making a general wrapper due to the issue of naming conventions, but I like the idea too.

Edit: Here's what I have so far ... YMMV, it's still kinda rough. https://github.com/lamorton/SciPyXarray

I've written wrappers for svd, fft, psd, gradient, and specgram, for starts

@lamorton I really like the suggestion from @shoyer about submodules for throwing wrappers from other libraries, but in the meantime I think I might like very much to check out your implementation of fft and gradient in particular if these are somewhere public. I have been hacking at at least the latter and other functions in the numpy/scipy scope.

Hello. I discovered xarray a few days ago, and find it very useful for my job. Integral along a coordinate is one of few things which I found missing so far.

Adding my +1 without offering to do the work. :) This would be very welcome!

+1 for integrate. I found this thread when having the same problem.

Perhaps a new package would be in order?

I would also be very happy to include many of these in a submodule inside xarray, e.g., xarray.scipy for wrappers of the scipy API. This would make it easier to use internal methods like apply_ufunc (though hopefully that will be public API soon).

If you give a mouse a cookie, he'll ask for a glass of milk. There are a whole slew of Numpy/Scipy functions that would really benefit from using xarray to organize input/out. I've written wrappers for svd, fft, psd, gradient, and specgram, for starts. Perhaps a new package would be in order?

By the way, the cumtrapz implementation I pasted above matches the scipy version when initial=0, which I also think would be a more sane default for integration.

Yes, I agree with both of you that we should fix initial=0. (I don't know if I would even bother with adding the option.)

As far as implementation is concerned. Is there any performance downside to using xarrays shift operators versus delving deeper into dask with map_blocks, etc? I looked into using dasks cumreduction function, but am not sure it is possible to implement the trapezoid method in that way without changing dask.

From a performance perspective, it would be totally fine to implement this either in terms of high level xarray operations like shift/sum/cumsum (manipulating full xarray objects) or in terms of high level dask.array operations like dask.array.cumsum (manipulating dask arrays). I would whatever is easiest. I'm pretty sure there is no reason why you need to get into dask's low-level API like map_blocks and cumreduction.

I usually agree that using too many (or any) switches within functions is not ideal. However, I think this is more important for low level or internal routines. For user facing interfaces, I think it is okay. After all, many numpy and scipy functions have convenient switches that control the return values.

By the way, the cumtrapz implementation I pasted above matches the scipy version when initial=0, which I also think would be a more sane default for integration.

As far as implementation is concerned. Is there any performance downside to using xarrays shift operators versus delving deeper into dask with map_blocks, etc? I looked into using dasks cumreduction function, but am not sure it is possible to implement the trapezoid method in that way without changing dask. On Mon, Mar 20, 2017 at 8:48 AM Fabien Maussion notifications@github.com wrote:

An argument against a single function is that the shape of the returned array is different in each case. Also, cumtrapz https://docs.scipy.org/doc/scipy-0.10.1/reference/generated/scipy.integrate.trapz.html has an inital keyword which changes the shape of the returned array. It is currently set to None per default, but should be set to 0 per default IMO.

I this is not a problem, I also like to have one single function for integration (simpler from a user perspective).

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/pydata/xarray/issues/1288#issuecomment-287749414, or mute the thread https://github.com/notifications/unsubscribe-auth/ABUokrAdpysuufZxHdLSdc1nseH9PtOkks5rnnWYgaJpZM4MOCxc .

An argument against a single function is that the shape of the returned array is different in each case. Also, cumtrapz has an inital keyword which changes the shape of the returned array. It is currently set to None per default, but should be set to 0 per default IMO.

I this is not a problem, I also like to have one single function for integration (simpler from a user perspective).

Sorry for letting this lapse.

Yes, we absolutely want this functionality in some form.

My concern is that this doesn't feel like functionality that inherently belongs as a method on a DataArray--if doesn't need to be a method, it shouldn't be. In numpy and scipy, these are separate functions and I think they work fine that way.

This is a fair point, and I agree with you from a purist OO-programming/software-engineering perspective (TensorFlow, for example, takes this approach). But with xarray, we have been taking a different path, putting methods on objects for the convenience of method chaining (like pandas). So from a consistency perspective, I think it's fine to keep these as methods. This is somewhat similar even to NumPy, where a number of the most commonly used functions are also methods.

Perhaps allow generic extension of da.integrate by letting the method keyword of da.integrate accept a function as an argument that performs the actual integration?

I don't see a big advantage to adding such an extension point. Almost assuredly it's less text and more clear to simply write ds.pipe(my_integrate, 'x') or my_integrate(ds, 'x') rather than ds.integrate('x', my_integrate).

Maybe this could be implemented by adding an optional cumulative flag.

I normally don't like adding flags for switching functionality entirely but maybe that would make sense here if there's enough shared code (e.g., simply substituting cumsum for sum). The alternative is something like cum_integrate which sounds kind of awkward and is one more additional method.

One thing that can be useful to do before writing code is to write out a docstring with all the bells and whistles we might eventually add. So let's give that a shot here and see if integrate still makes sense: ``` integrate(dim, method='trapz', cumulative=False)

Arguments

dim : str or DataArray DataArray or reference to an existing coordinate, labeling what to integrate over. cumulative : bool, optional Whether to do a non-cumulative (default) or cumulative integral. method : 'trapz' or 'simps', optional Whether to use the trapezoidal rule or Simpson's rule. ```

I could also imagine possibly adding a bounds or limits argument that specifies multiple limits for controlling multiple integrals at once (e.g., dim='x' and bounds=[0, 10, 20, 30, 40, 50] would result in an x dimension of length 5). This would certainly be useful for some of my current work. But maybe we should save this sort of add for later...

I would also like to see an integrate function. I have had one monkey patched in my own xarray routines for a while now. Also wanted: cumtrapz and friends. Maybe this could be implemented by adding an optional cumulative flag. This shouldn't be too hard to do. For example, in the following cumtrapz implementation all that would need to be changed is the final cumsum call.

```python def cumtrapz(A, dim): """Cumulative Simpson's rule (aka Tai's method)

Notes
-----
Simpson rule is given by
    int f (x) = sum (f_i+f_i+1) dx / 2
"""
x = A[dim]
dx = x - x.shift(**{dim:1})
dx = dx.fillna(0.0)
return ((A.shift(**{dim:1}) + A)*dx/2.0)\
      .fillna(0.0)\
      .cumsum(dim)

```

The issue is that certain types of gridded data (such as output from numerical models) should actually not be integrated with the trapezoidal rule but rather should use the native finite volume discretization for their computational grid.

• We are aiming for the 20% of functionality that covers 80% of use cases, not the long tail.
• We don't want implementations of any complex numerical methods in xarray (like NumPy rather than SciPy).

I can see the problems down the road that @rabernat brings up. Say you have a high-order finite volume discretization and some numerical implementation of high-order integration for that gridding. What would your interface be? You could write it as new_integrate(da, dim, domain) but then it may be confusing to have da.integrate be different (and less accurate).

That might bring us back to the algorithmically descriptive name trapz, but then what about @shoyer's point that given a fixed gridding, da.integrate is the most readable choice of name? Perhaps allow generic extension of da.integrate by letting the method keyword of da.integrate accept a function as an argument that performs the actual integration?

👍 for the functionality (both intergrate and gradient) that work with DataArray. My concern is that this doesn't feel like functionality that inherently belongs as a method on a DataArray--if doesn't need to be a method, it shouldn't be. In numpy and scipy, these are separate functions and I think they work fine that way.

Another way to look at it is that methods are there to encapsulate some kind of manipulation of internal state or to ensure that some kind of invariant is maintained. I don't see how integrate is doing any of this for DataArray--seems like everything integrate would do would be doing can be accomplished using the public API. So really what you're buying is doing this: python da.integrate(dim='x', method='trapezoidal') instead of python integrate(da, dim='x', method='trapezoidal`) If you want to see what the pathological case of putting everything as a method for convenience looks like, go look at all the plot methods on matplotlib's Axes class. Pay special attention to the tangled web of stuff that comes from having ready access to the class's internals.

My real preference would just to have this work: python ds = xr.tutorial.load_dataset('rasm') np.trapz(ds['Tair'], axis='x') but I have no idea what that would take, so I'm perfectly fine with xarray gaining its own implementation.

I like the integrate idea. Nothing further to add not already covered nicely via the above concerns by @rabernat and responses by @shoyer.

As usual @rabernat raises some excellent points!

I weakly prefer not to use the name integrate and instead keep the standard scipy names because they make clear the numerical algorithm that is being applied.

Yes, this is a totally valid concern, if a user might expect integrate to be calculating something different.

One point in favor of calling this integrate is that the name is highly searchable, which provides an excellent place to include documentation about how to integrate in general (including links to other packages, like pangeo's vector calculus package). But we know that nobody reads documentation ;).

But where does it end? Why not implement the rest of the scipy.ode module?

Looking at the rest of scipy.integrate, in some ways the functionality of trapz/cumtrapz/simps is uniquely well suited for xarray: they are focused on data ("given fixed samples") rather than solving a system of equations ("given a function").

numpy.gradient feels very complementary as well, so I could see that as also in scope, but there are similar concerns for the name. There might be some value in complementary names for integrals/gradients.

As a community we need to develop a roadmap that clearly defines the scope of xarray.

I doubt we'll be able to come up with hard and fast rules, but maybe we can enumerate some principles, e.g.,

• Features should be useful to users in multiple fields.
• Features should be primarily about working with labeled data.
• We are aiming for the 20% of functionality that covers 80% of use cases, not the long tail.
• We don't want implementations of any complex numerical methods in xarray (like NumPy rather than SciPy).
• Sometimes it's OK to include a feature in xarray because it makes logical sense with the rest of the package even if it's slightly domain specific (e.g., CF-conventions for netCDF files).

I weakly prefer not to use the name integrate and instead keep the standard scipy names because they make clear the numerical algorithm that is being applied.

Yes, this is a totally valid concern, if a user might expect integrate to be calculating something different. One point in favor of calling this integrate is that the name is highly searchable, which provides an excellent place to include documentation about how to integrate in general (including links to other packages, like pangeo's vector calculus package). But we know that nobody reads documentation ;).

integrate would allow to do things like:

da.integrate(how='rectangle')

da.integrate(how='trapezoidal')

And I'm fine with integrate if that is the consensus here.

Having an xarray wrapper on trapz or cumtrapz would definitely be useful for many users. I weakly prefer not to use the name integrate and instead keep the standard scipy names because they make clear the numerical algorithm that is being applied. The issue is that certain types of gridded data (such as output from numerical models) should actually not be integrated with the trapezoidal rule but rather should use the native finite volume discretization for their computational grid. The goal of our hypothetical pangeo vector calculus package is to implement integrals and derivatives in such a context. A built-in xarray integration function would apply in cases where the data is assumed to be continuous, and where no auxiliary information about the grid (beyond the coordinates) is available.

I will also make the same comment I always make when such feature requests are raised: yes, it always seems desirable to add new features to xarray on a function-by-function basis. But where does it end? Why not implement the rest of the scipy.ode module? And why stop there? As a community we need to develop a roadmap that clearly defines the scope of xarray. Once apply is stable, it might not be that hard to wrap a large fraction of the scipy library. But maybe that should live in a separate package.

An integrate method is probably better for the reason you describe---it's more obvious. I believe the name trapz came from Matlab originally.

With a general integrate, it's probably also useful to allow optional input arguments for lower_bound and upper_bound as a convenience for integrating over a subset of the data instead of the user doing that in a slice. If those arguments aren't given, they would default to all of the data.

I agree that the API should mostly copy the mean/sum reduce methods (and in fact the implementation could probably share much of the logic). But there's still a question of whether the API should expose multiple methods like DataArray.trapz/DataArray.simps or a single method like DataArray.integrate (with method='simps'/method='trapz').

As long as there isn't something else we'd want to reserve the name for, I like the sound of integrate a little better, because it's more self-descriptive. trapz is only obvious if you know the name of the NumPy method. In contrast, integrate is the obvious way to approximate an integral. I would only hold off on using integrate if there is different functionality that comes to mind with the same.

It looks like SciPy implements Simpson's rule with the same API (see scipy.integrate.simps), so that would be easy to support, too. Given how prevalent SciPy is these days, I would have no compunctions about making scipy required for this method and defaulting to method='simps' for DataArray.integrate.

It would be useful to have dask.array versions of these functions, too, but that's not essential for a first pass. The implementation of trapz is very simple, so this would be quite easy to add to dask.

CC @spencerahill @rabernat @lesommer in case any of you have opinions about this

+1 for integrate

The cumulative integral is of very frequent use in atmospheric sciences, too : https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.integrate.cumtrapz.html

I don't at the moment see a reason to use a different API than DataArray.mean or DataArray.sum. DataArrays assume a default spacing of 1 if coordinates are not given, which is exactly what np.trapz does. So the API for trapz might look like:

DataArray.trapz(dim=None, axis=None, skipna=None, keep_attrs=False, **kwargs)