Implementation of FFT geophysical filters in Harmonica

[santisoler]

Hi everyone!

I'm Santiago Soler, one of the core developers of the Fatiando a Terra project: open-source tools for geophysics. We have been quite interested in the development of xrft, specially for managing FFT-based processes and filters that we are planning to include in Harmonica, a library for processing, forward and inverse modelling gravity and magnetic data.

We are very interested in using xrft for managing the whole FFT and inverse FFT processes that these functions need. Our potential fields data (gravity and magnetic) is usually given as a 2D regular grids whose coordinates are projected geographical coordinates: easting and northing (Cartesian coordinates in meters). See this gallery example to get a better picture of them. In the old fatiando package we used to have some implementations of these methods using numpy.fft directly, although we are very interested in a FFT implementation that harness xarray objects, since we are using DataArrays and Datasets as our default structure for our 2D grids.
We therefore find xrft the perfect solution to our problems. We really appreciate you developers for taking the time to build this great project!

I personally started drafting an implementation of a directional derivatives of a regular grid using FFTs (fatiando/harmonica#238) using xrft under the hood for computing the FFTs and their inverse.
On this process, I started getting some user experience of xrft and wanted to share them with you.
We also found that we might have some special needs and we would like to know if you would want to include solutions to them into xrft. We are absolutely willing to help with the implementation and maintenance of those.

Differences between dft and fft

We found that xrft hosts two different functions for handling FFT computations: dft and fft.
We found out that fft is a wrapper of dft that passes hard-coded values to certain parameters.
IMHO, it sounds slightly confusing at the beginning, and in fact I had to surf through the code to finally understand the difference between them.

Is there any strong decision for this design?

Mathematical docstring of dft

The docstring of dft states that it returns the Fourier transform computed as

daft = \mathbb{F}(da - \overline{da})

where I assume that the \overline states the mean value of the array. Am I right?

If so, the function only removes the mean value if we pass detrend="linear", right?

True Phase and True Amplitude

I've been struggling to understand what true_amplitude and true_phase actually mean.
I think it would help to add some mathematical background on the intended behavior of enabling and disabling them.

Coordinates shift

Most of the FFT-related functions we will implement on Harmonica can be summarized in the following workflow:

• compute FFT of a 2D regular grid
• apply a filter in the frequency domain
• compute inverse FFT to obtain the modified grid on the same coordinates.

While implementing the derivatives (fatiando/harmonica#238) I've noticed that xrft.xrft.ifft returns a DataArray whose coordinates are centered around zero. Although this might not be a big deal in the general case, for the Harmonica use cases it's extremely important that the output grid has the same coordinates as the original one, since they mean actual geographical locations. I tried passing a tuple to the lag parameter with the western and southern coordinates of the original region, but ifft returns a warning saying that it ignores this parameter.

So far, we manage to implement this feature on our side: https://github.com/fatiando/harmonica/blob/4e865f00f0464e24d8f1cbb8b3c825ab70aa899d/harmonica/filters/fft.py#L34-L85

Nevertheless we don't intent to maintain wrappers of the xrft functions in the future.

Would you like to see something like this implemented in xrft?

Plans for future padding functions

In the future we will be very interested in having functions to apply padding to our 2D grids before passing them to fft. We've been following the implementation of xarray.pad, although it currently fills with nans the extended coordinates which cause a lot of problems when trying to compute FFTs (see pydata/xarray#3868). Since xarray is not intended to be used exclusively with regular grids, extending the coordinates of a DataArray is not trivial. But since xrft assumes regular grids, will you be interested in implementing specific padding functions in xrft?

---

Thanks again for all your work on this subject! We really appreciate it.
And sorry for my long issue. If you think it might be better to split this Issue into several ones, please let me know and I'll split it so we can chat about every topic in a more organized way.

Cheers!

cc @leouieda @aguspesce @andieie

[roxyboy]

Hi @santisoler, thanks for your interest in xrft and detailed feedback. We are happy to know that there are people using our package.

Differences between dft and fft

We found that xrft hosts two different functions for handling FFT computations: dft and fft.
We found out that fft is a wrapper of dft that passes hard-coded values to certain parameters.
IMHO, it sounds slightly confusing at the beginning, and in fact I had to surf through the code to finally understand the difference between them.

Is there any strong decision for this design?

Part of it was to be compatible with prior releases. We decided to use the naming of fft for similar behavior to numpy.fft while as the functionality of passing the xarray metadata was left for dft.

Mathematical docstring of dft

The docstring of dft states that it returns the Fourier transform computed as

daft = \mathbb{F}(da - \overline{da})

where I assume that the \overline states the mean value of the array. Am I right?

If so, the function only removes the mean value if we pass detrend="linear", right?

detrend="constant" removes the mean while detrend="linear" removes the linear trend either in 1D or 2D.

True Phase and True Amplitude

I've been struggling to understand what true_amplitude and true_phase actually mean.
I think it would help to add some mathematical background on the intended behavior of enabling and disabling them.

The documentation page (https://xrft.readthedocs.io/en/latest/) has some examples with mathematical/theoretical explanation but feel free let us know if you have ideas on how to improve it.

Coordinates shift

Most of the FFT-related functions we will implement on Harmonica can be summarized in the following workflow:

• compute FFT of a 2D regular grid
• apply a filter in the frequency domain
• compute inverse FFT to obtain the modified grid on the same coordinates.

While implementing the derivatives (fatiando/harmonica#238) I've noticed that xrft.xrft.ifft returns a DataArray whose coordinates are centered around zero. Although this might not be a big deal in the general case, for the Harmonica use cases it's extremely important that the output grid has the same coordinates as the original one, since they mean actual geographical locations. I tried passing a tuple to the lag parameter with the western and southern coordinates of the original region, but ifft returns a warning saying that it ignores this parameter.

The lag parameter is only valid when using xrft.idft and not xrft.ifft.

Plans for future padding functions

In the future we will be very interested in having functions to apply padding to our 2D grids before passing them to fft. We've been following the implementation of xarray.pad, although it currently fills with nans the extended coordinates which cause a lot of problems when trying to compute FFTs (see pydata/xarray#3868). Since xarray is not intended to be used exclusively with regular grids, extending the coordinates of a DataArray is not trivial. But since xrft assumes regular grids, will you be interested in implementing specific padding functions in xrft?

This may be related also to issue #146 ?

[santisoler]

Hi @roxyboy ! Thanks for the quick response.

Part of it was to be compatible with prior releases. We decided to use the naming of fft for similar behavior to numpy.fft while as the functionality of passing the xarray metadata was left for dft.

I understand. My question is: does it worth having two different functions that are intended to do the same thing with slight differences?

IMHO, I find it a little bit confusing. Nevertheless, if you find necessary to maintain both functions, I would suggest to improve their docstrings to make it very clear about their differences.

detrend="constant" removes the mean while detrend="linear" removes the linear trend either in 1D or 2D.

Right! My concern was about the difference between the math description on the docstring and the default behavior of the function. By default dft sets detrend=None, which corresponds to not removing any trend from the array. I would suggest that the docstrings should match the default behavior of the function.

The documentation page (https://xrft.readthedocs.io/en/latest/) has some examples with mathematical/theoretical explanation but feel free let us know if you have ideas on how to improve it.

I like the documentation! I would just add mathematical descriptions of what the dft function is doing when passing true_amplitude=True and true_phase=True.
For example, Numpy has an extended explanation about their fft method: https://numpy.org/doc/stable/reference/routines.fft.html#module-numpy.fft

The lag parameter is only valid when using xrft.idft and not xrft.ifft.

Is there any reason why it wouldn't work with ifft?

This may be related also to issue #146 ?

Yes! I'll bring my ideas on #146 then.

---

Thanks again for the response! We are really looking forward to use xrft as the base tool for managing FFTs on Harmonica.

[roxyboy]

I understand. My question is: does it worth having two different functions that are intended to do the same thing with slight differences?

IMHO, I find it a little bit confusing. Nevertheless, if you find necessary to maintain both functions, I would suggest to improve their docstrings to make it very clear about their differences.

Maybe @lanougue can chime to remind me of the discussion we had but I think the idea was that in order to avoid confusion between numpy.fft and xarray.dft, we decided to keep xrft.fft to maintain a functionality similar to numpy.fft.

Right! My concern was about the difference between the math description on the docstring and the default behavior of the function. By default dft sets detrend=None, which corresponds to not removing any trend from the array. I would suggest that the docstrings should match the default behavior of the function.

Noted. Will do this in a new PR.

I like the documentation! I would just add mathematical descriptions of what the dft function is doing when passing true_amplitude=True and true_phase=True.
For example, Numpy has an extended explanation about their fft method: https://numpy.org/doc/stable/reference/routines.fft.html#module-numpy.fft

Ok, I'll add this to the documentation.

Is there any reason why it wouldn't work with ifft?

xrft.fft and xrft.ifft sheds all the coordinate data (as we'd expect from numpy.fft) so we decided that lag should only work with dft which utilizes the coordinate metadata.

[lanougue]

Hi @santisoler,
The initial idea of xrft was to extend fft's to xarrays handling frequency coordinates. However, since xarrays also contain coordinates information, it is possible to correctly perform a Fourier Transform (meaning : with the correct amplitude and correct phase). xrft.fft is the historical function that simply mimic numpy.fft (with frequency vector) while xrft.dft takes care of correct amplitude and phasing of the output. Enabling the correct ampltude and correct phasing compared to xrft.fft is possible with flags "true_amplitude" and "true_phase".
We beleive that maintaining xrft.fft was important for the users for retro-compatibility and to reassure users that would still prefer fft anyway.
Maybe, in the future, only xrft.dft will be kept and documentation will explain that "old fft" behaviour could be recovered with true_amplitude=False and true_phase =False.
I agree that documentation should be more explicit (mea culpa, I am slow in reviewing @roxyboy's work)

"lag flag" could also work with ifft but would not have any signification. fft works together with ifft and dft works together with idft. (However, only dft/idft correctly handle the phase).

[lanougue]

@santisoler,
The "lag" parameter exists in order to be able to get the idft onto the right range of coordinates (same as initial input). However, from a user point of view: do you think that it could be easier to provide the original signal's coordinates instead of the lag?

[leouieda]

@lanougue just jumping in here 🙂 Would it be possible to store this information in the Dataset somehow so that users don't have to pass it separately? This could be through attributes or maybe just keeping the space/time-domain coordinates in the frequency domain Dataset so that they can be restored later.

[lanougue]

Hi @leouieda
I think that keeping the full coordinates would be a waste of memory space since they can easily be regenerated (and will be anyway when calling idft). However, the lag is a single value and would require almost no memory.

It could be possible to store the "lag" (evaluated during dft call) but I would prefer rather not for several reasons:

• It is only useful when calling both dft and then idft. Most of the users only use the direct dft to look at the spectrum for example. This lag information is useless for them.
• Even if the lag evaluated when calling dft is the more probable value used in idft, it is not guaranted.

As I suggested, I think it could be useful to enable passing the initial dataarray coordinates to idft in order to make idft to automattically evaluate the lag.

Maybe @rabernat and @roxyboy have other pros/cons arguments.

[rabernat]

Thanks so much @santisoler and @leouieda for this really useful feedback. We are 💯 committed to making xrft a useful, stable dependency for downstream projects such as yours. I'll chime in with a few of my thoughts.

Part of it was to be compatible with prior releases. We decided to use the naming of fft for similar behavior to numpy.fft while as the functionality of passing the xarray metadata was left for dft.

IMO we should strongly consider deprecating one or the other. dft is not a standard function name, while fft appears in many different packages. I think the first step would be to give these duplicate functions the same API (together with a warning that one will be dropped).

The documentation page (https://xrft.readthedocs.io/en/latest/) has some examples with mathematical/theoretical explanation but feel free let us know if you have ideas on how to improve it.

Our documentation basically only has examples. What is missing is the standard narrative documentation that describes the operation and theory of xrft in detail.

I think that keeping the full coordinates would be a waste of memory space since they can easily be regenerated (and will be anyway when calling idft). However, the lag is a single value and would require almost no memory.

IMO, we should not be concerned about memory footprint for coordinates. If you are using Xarray, you already have decided you want to accept the extra memory overhead in exchange for keeping careful track of coordinates.

• Most of the users only use the direct dft to look at the spectrum for example. This lag information is useless for them.

@lanougue: How do you know this? Have we done a user survey? Not to my knowledge. Here we have some actual users coming and requesting this feature. We can see from the harmonica workaround that this is a feature that users need. We should not be so quick to say no. As I said above, I don't think that conserving memory is a valid reason. People use xrft (and python in general) for its convenience, not its memory efficiency.

In general, I think that xrft should support full "round-tripping" of data, such that ifft(fft(data)) returns a dataarray with the original coordinates of data.

I think it could be useful to enable passing the initial dataarray coordinates to idft in order to make idft to automattically evaluate the lag

Rather than placing this burden of keeping track of coordinate information on the user, I propose we use attributes and / or extra non-dimension coordinates to provide the necessary information about lag within the dataarray returned by fft.

[roxyboy]

IMO we should strongly consider deprecating one or the other. dft is not a standard function name, while fft appears in many different packages. I think the first step would be to give these duplicate functions the same API (together with a warning that one will be dropped).

dft is the legacy naming in xrft so if we are to drop either fft or dft, for compatibility with prior versions, I vote for dropping fft.

Our documentation basically only has examples. What is missing is the standard narrative documentation that describes the operation and theory of xrft in detail.

PR #151 adds some of the the operation and theory, which I think we can merge.

IMO, we should not be concerned about memory footprint for coordinates. If you are using Xarray, you already have decided you want to accept the extra memory overhead in exchange for keeping careful track of coordinates.
In general, I think that xrft should support full "round-tripping" of data, such that ifft(fft(data)) returns a dataarray with the original coordinates of data.
Rather than placing this burden of keeping track of coordinate information on the user, I propose we use attributes and / or extra non-dimension coordinates to provide the necessary information about lag within the dataarray returned by fft.

I'll work on this in another PR.

[roxyboy]

IMO we should strongly consider deprecating one or the other. dft is not a standard function name, while fft appears in many different packages. I think the first step would be to give these duplicate functions the same API (together with a warning that one will be dropped).

dft is the legacy naming in xrft so if we are to drop either fft or dft, for compatibility with prior versions, I vote for dropping fft.

Or perhaps the other way around is better: dft will be a wrapper function of fft with the latter being the default and whenever dft is called, it'll raise a deprecation warning.

[lanougue]

I did not say "no" ! I only started a discussion :)
I like the round-tripping argument :)
@roxyboy, I can start this PR too if you need

[rabernat]

Sure, sorry if I mischaracterized your comments @lanougue. My apologies. I appreciate everyone's viewpoint on this discussion.

[roxyboy]

dft is the legacy naming in xrft so if we are to drop either fft or dft, for compatibility with prior versions, I vote for dropping fft.

Or perhaps the other way around is better: dft will be a wrapper function of fft with the latter being the default and whenever dft is called, it'll raise a deprecation warning.

I made a PR #155 to start the deprecation cycle of dft (coined as byebyedft)

[roxyboy]

I like the round-tripping argument :)
@roxyboy, I can start this PR too if you need

Sure, if you could start a PR about keeping the original coordinates as an attribute and automating the lag computation, that'd be great @lanougue :)

[lanougue]

I like the round-tripping argument :)
@roxyboy, I can start this PR too if you need

Sure, if you could start a PR about keeping the original coordinates as an attribute and automating the lag computation, that'd be great @lanougue :)

Ok, I'll initiate it !

[santisoler]

Thanks everyone for the amazing responses.

I left a few thoughs after reading all the comments in this Issue.

In general, I think that xrft should support full "round-tripping" of data, such that ifft(fft(data)) returns a dataarray with the original coordinates of data.
by @rabernat

I totally agree on this. I think the fft and ifft functions should resemble the mathematical description of the Fourier transform and its inverse as much as possible. And since the F^-1 * F operator is equal to the identity, I cannot agree more to support the round-tripping.

xrft.fft is the historical function that simply mimic numpy.fft (with frequency vector) while xrft.dft takes care of correct amplitude and phasing of the output. Enabling the correct ampltude and correct phasing compared to xrft.fft is possible with flags "true_amplitude" and "true_phase".
We believe that maintaining xrft.fft was important for the users for retro-compatibility and to reassure users that would still prefer fft anyway.
by @lanougue

IMHO, if any user would like to apply the FFT to a xr.DataArray regardless of the frequency arrays or its coordinates, they can always go back to numpy.fft instead of having to go with xrft. Therefore I think that xrft should be focused on the users that are specially worried of keeping track of the frequencies and coordinates arrays.

Or perhaps the other way around is better: dft will be a wrapper function of fft with the latter being the default and whenever dft is called, it'll raise a deprecation warning.
by @roxyboy

I agree with keeping fft and deprecating dft. FFT stands for Fast Fourier Transform, which actually stands for a particular algorithm for computing discrete Fourier transforms that it's much more faster than carrying out the entire summation by iteration. The first time I saw dft I thought that xrft was actually using a different algorithm than FFT. Therefore I agree xrft should keep fft and ifft as the default functions to apply Fourier transformations.

Rather than placing this burden of keeping track of coordinate information on the user, I propose we use attributes and / or extra non-dimension coordinates to provide the necessary information about lag within the dataarray returned by fft.
by @rabernat

If we are all in the same page regarding the round-tripping property of fft, I think this feature is a must. Probably it's easier to add the lag coordinates as attrs, although I'm worried if the attrs could get deleted or modified by any intermediate process.

@leouieda: I think you mention something like this at some point, could you bring some light on this?

It think coordinates and not-index coordinates cannot be easily removed or modified. Nevertheless which option we choose, I think we should stick with the safer one.

Finally but not least, I'm willing to help opening PRs or reviewing code. Feel free to assign me any task you want me to help with!

[lanougue]

Hi @santisoler
We are doing PRs according to your suggestions ! Thanks. Concerning the accidental removal of the lag attribute during intermediate processes (in the frequency domain). I proposed (in the PR) to store the lag as an attribute of the coordinates (and not of the dataarray). This should ensure that it will not be deleted during intermediate operations on this dataarray. (it is also actually an attribute of the coordinates per definition).
The lag will be changed/removed only if operations are realized on the coordinates themselves. In this case, the lag should probably be changed anyway.

[lanougue]

I think the fft and ifft functions should resemble the mathematical description of the Fourier transform and its inverse as much as possible.

Actually this question raised up some time ago. Should xrft.fft mimic the numpy.fft (because users are used to it) or should it be the closest as possible to the mathematical Fourier transform, as you suggest.
At time, the chosen option was to make xrft.fft to mimic numpy.fft and to make xrft.dft the closest to the mathematical Fourier Transform (I agree that the dft name was not a very good choice but it was for historical reasons).

We now need to choose how xrft.fft should behave. I also raised this question in PR #155 too because it is of central interest

[santisoler]

Hi @santisoler
We are doing PRs according to your suggestions ! Thanks. Concerning the accidental removal of the lag attribute during intermediate processes (in the frequency domain). I proposed (in the PR) to store the lag as an attribute of the coordinates (and not of the dataarray). This should ensure that it will not be deleted during intermediate operations on this dataarray. (it is also actually an attribute of the coordinates per definition).
The lag will be changed/removed only if operations are realized on the coordinates themselves. In this case, the lag should probably be changed anyway.

This is much more clever! I totally agree, so feel free to ignore my previous proposals regarding where to place the lags.